11institutetext: NOVA LINCS & DI, FCT, Universidade Nova de Lisboa

# A Deductive Verification Framework For Higher Order Programs

Tiago Soares — 52965
Supervisors: Mário Pereira  António Ravara
<EMAIL_ADDRESS>
###### Abstract

In this report, we present the preliminary work developed for our research
project for the APDC (Área Prática de Desenvolvimento Curricular) course. The
main goal of this project is to develop a framework, on top of the Why3 tool,
for the verification of effectful higher-order programs. We use
defunctionalization as an intermediate transformation from higher-order OCaml
implementations into first order ones. The target for our translation is
WhyML, the Why3’s programming language. We believe defunctionalization can be
an interesting route for the automated verification of higher-order programs,
since one can employ off-the-shelf automated program verifiers to prove the
correctness of the generated first-order program. This report also serves to
introduce the reader to the subject of deductive program verification and some
of the tools and concepts used to prove higher order effectful programs.

###### Keywords:

Deductive Software Verification Imperative Programs Defunctionalization Why3
OCaml

##  1 Introduction

Producing error free code for a complex system is, regardless of the
programmer’s skills, a herculean task. Customarily, the only weapons we have
to combat the plethora of bugs rotting away our code are tests and plenty of
caution. But caution can only get you so far and program testing, in
Dijkstra’s words, "can be a very effective way to show the presence of bugs,
but is hopelessly inadequate for showing their absence"[6]. Therefore, in
order to make sure that our code does what we want, we need formal proofs of
its correctness.

There is no one definition of correctness, however, considering it depends on
how we want to specify our program. We might simply want to make sure that
there are no runtime errors involving type safety, arithmetic overflows or
unhandled exceptions. We might also want to prove that the program eventually
terminates. More generally, we might want to prove some relation to its input
and/or the state with its output. For each definition of correctness there are
better ways of proving it. For the class of programs we will analyse, that is,
functional programs, the best alternative is deductive verification.[7]

Dispatching proofs for imperative implementations can be challenging, as they
heavily rely on modifying data structures after they have been created, which
can be hard to model using first order logic. Therefore, writing in a
programming language that stays as close as possible to logical definitions,
that is, a functional language, means we will have an easier time identifying
the necessary assertions to create an adequate proof. OCaml is a good
candidate for writing effectful, verifiable code, not only because it is
primarily a functional language, but also because it allows us to use
imperative constructs, such as mutable records.

To write specifications for code with side effects, we need additional
infrastructure on top of first order logic. One of the most popular is
Separation Logic [17], which allows us to connect physical memory to logical
statements. Although it is very expressive, its intrinsic technicalities
prevent it from being widely adopted. In this work, we will be using GOSPEL
(Generic OCaml SPEcification Language) [2], which is built on top of
Separation Logic, as a higher level alternative to define changes in the heap
and, to prove the specification, we will use Why3 [9], a program verifier
tailored towards automated verification.

But effectful code is only one half of this project; the other being higher-
order functions. As previously said, code with side effects can be difficult
to map onto logical definitions, but this problem becomes exponentially more
complex when these effects originate from first class functions, due to the
way programs are proved. In short, we annotate them with pre and
postconditions and check that an implementation conforms to the provided
specification. This is normally done by employing some technique of
verification conditions generation, such as Dijkstra Weakest Precondition
Calculus [5]. However, if we don’t know the behaviour of a part of our
routine, and this unknown part might alter our program’s state, then it can be
hard to prove if a given program is correct.

Some work has been put into bringing closer automated verification and higher-
order programs with effects, such as the tool Who [12]. Nonetheless, a
framework of that caliber can be complex to implement and maintain. Another
option would be Iris [11], one of the most powerful tools available in the
realm of program verification, but whose power comes at a cost. The logic
employed by it is, albeit very rich, extremely complex and almost inaccessible
to anyone without years of training in the field. In addition, there is no way
to automate proofs with it, this being the main reason why Iris falls out of
the scope of our project.

Most options available for handling this problem present their own hindrances,
such as, for example, their complexity, cross-compatibility and their
practical usefulness [13, Sec. 3.4]. We have put forth a small contribution to
solve these limitations by creating a Why3 plugin which makes use of
defunctionalization [4, 16] to handle the verification of such programs.
Defunctionalization works by translating a higher-order implementation into an
equivalent first-order counterpart. We believe this technique can open the
door to automatic verification of higher-order code, since we can use an
existing first-order verifier to prove the defunctionalized program.

Previous work studied the combination of defunctionalization and software
verification. [14]. However, no prototype was developed and all examples were
made entirely by hand. In this project, we lay the groundwork for a fully-
fledged verification framework for effectful higher-order programs, by
exploring the technique of defunctionalization. To the best of our knowledge,
this is the first practical implementation of defunctionalization as a proof
technique. Our framework will accept as input a higher-order program,
annotated with assertions written in GOSPEL, translate it into a first-order
equivalent WhyML program so that we can use Why3 to prove it satisfies the
assertions.

##  2 Why3

Why3 is a verification framework that allows us to interface with multiple
theorem provers and automated solvers using the same surface syntax and
underlying logic [8, 1]. This means that, if one such solver is not able to
discharge a verification condition, one can easily switch to another which
might be able to complete our proof.

To understand how it works, we will present the defenition of an array list
type written in Why3’s dedicated programming language WhyML:

⬇

type arrayList ’a = {

mutable arr : array ’a;

mutable size : int

}

This type definition features a field arr of type array ’a, the polymorphic
type of arrays defined in Why3 standard library, and an integer size
representing the number of elements in a value of type arrayList.

Although our type declaration has all we need to build a comprehensive set of
functions for our arrayList, from a logical point of view we are missing
crucial information. In particular, that the size attribute cannot be greater
than the length of the array. To express this, we will add the following type
invariant:

⬇

invariant { 0 <= size <= arr.length }

Next, we will write a function that replaces the element at a given index and
returns the old element:

⬇

let set (al: arrayList ’a) (i : int) (e : ’a) : ’a

For simplicity’s sake, we will add a _precondition_ requiring i to be within
the bounds of our list:

⬇

requires { 0 <= i < al.size }

Had we not included such precondition, we would have to include in the body of
set an if..then..else expression to avoid accessing an index out of bounds.
This a common practice of _defensive programming_.

Next, we will need to specify that the result of our function is the old value
at index i and that the new value is e. To do this, we will will introduce two
_postconditions_ , where we use the old tag to refer to our array list before
the function was called and result for our function’s return value

⬇

ensures { result = old al.arr[i] }

ensures { al.arr[i] = e }

Lastly, we also need to specify that every other element in the array, as well
as its size, is unchanged.

⬇

ensures { al.size = old al.size }

ensures { forall j. 0 <= j < al.size /\ not j = i ->

old al.arr[j] = al.arr[j] }

Now we simply write the body of the function, trivial as it is:

⬇

= let ret = al.arr[i] in

al.arr[i] <- e; ret

The full annotated function can be found in annex 0.A. When fed into Why3, the
above set function generates six verification conditions. These include safety
properties (access index in array bounds, which is proved with our invariant)
and functional ones (type invariant preservation and postcondition
entailment). All of these are automatically verified in less then a second by
the Alt-Ergo [10] SMT solver, one of the provers Why3 allows us to interface
with.

##  3 Defunctionalization

Before jumping into our framework’s implementation lets first take a gander at
a few quick (original) examples to understand how defunctionalization works,
how it can be mechanized and how we can apply it to program verification.

###  3.1 A simple example

As previously mentioned, defunctionalization is the practice of turning
higher-order programs into equivalent first-order implementations. To
illustrate this method, we will convert a higher-order OCaml program that
reverses a list. The code will be written in CPS (Continuation Passing Style)
[3], which, instead of letting the program return a value to the function’s
caller, passes the return values as a parameter to a “continuation” function
that is given as a first order value. The implementation is as follows:

⬇

let rec reverse_aux_cps l k = match l with

| [] -> k []

| x::t -> reverse_aux_cps t (fun r -> x :: (k r))

\parlet reverse l = reverse_aux_cps l (fun r -> r)

There are two anonymous functions: the identity function (fun r -> r) and the
function which builds the reverse of the list (fun r -> x::(k r)). In order to
begin our gentle descent towards a first order version of this program, we
introduce a new algebraic data type to represent continuations:

⬇

type ’a kont = Kid | Krev of ’a * ’a kont

The constructors Kid and Krev match, respectively, the identity function and
the one that builds the reverse. Krev‘s two extra arguments represent its free
variables, namely the appended value and the continuation function. Next, we
need to introduce an apply function which will replace all first order
functions calls:

⬇

let rec apply kont arg = match kont with

| Kid -> let r = arg in r

| Krev(x, k) -> let r = arg in x::(apply r k)

The defunctionalization process terminates by replacing, in the original code,
native function application with calls to apply and anonymous functions with
constructors of the kont data type.

⬇

let rec reverse_aux_defun l k = match l with

| [] -> apply [] k

| x::t -> reverse_aux_defun t (Krev(x, k))

\parlet reverse_defun l = reverse_aux_defun l Kid

The full defunctionalized implementation can be found in annex 0.B and if the
reader wishes to understand in more detail how this function would run, we
have included an example of the program running in annex 0.C.

###  3.2 Proving a defunctionalized program

With defunctionalization’s bases covered, let us introduce how it’s used in a
proof context. For this, we consider a simpler function, namely one that
computes a list’s size, also written in CPS. The WhyML implementation is as
follows:

⬇

let rec length_cps (l : list ’a) (k : int -> ’b) : ’b =

match l with

| [] -> k 0

| _ :: t -> length_cps t (fun l -> k (1 + l))

end

let len (l : list ’a) : int =

length_cps l (fun x -> x)

Now we annotate our functions with their respective postconditions. The len
function’s specification is fairly straightforward:

⬇

ensures { length l = result }

where length is a purely logical function, issued from the Why3 standard
library, which can only be used within specifications and shall never be mixed
with executable code. The CPS function’s specification is a bit more
intricate. We could, naively, write something like

⬇

ensures { k (length l) = result }

However, this leaves us with the problem of how to interpret the use of
functions in our logic, since these may have side effects (_e.g._ ,
divergence) which can lead to logical inconsistencies. Therefore, we introduce
a layer of abstraction to allow us to use potentially effectful programs in
our assertions: instead of thinking of functions as a computation of an output
from a given input, such as f: t1 -> t2, we will think of them as a pair of
predicates, namely their precondition (pre : t1 -> prop) and postcondition
(post : t1 -> t2 -> prop). This is strongly inspired by the work of Régis-
Gianas and Pottier [18]. Under this setting, we would provide with length_cps
the following specification:

⬇

ensures { post k (length l) result }

and the continuation defined in length_cps with

⬇

ensures { post k (l + 1) result }

Now we will proceed with our program’s defunctionalization as we did before,
first by defining a type to represent our first order functions:

⬇

type kont = Kid | Klen kont

The next step, that is, the apply function, is quite similiar to our previous
example, which can be found in annex 0.D. Finally, since we have no higher
order functions, we must define our own post predicate, as follows:

⬇

predicate post (k : kont) (length : int) (result : int) =

match k with

| Kid -> result = length

| Klen k -> post k (length + 1) result

end

For the above first-order implementation and specification, all the
verification conditions generated by Why3 are automatically discharged in no
time.

As we can observe, the defunctionalized definition of predicate post is a
simple matter of collecting the predicate’s uses from the original higher-
order implementation. As for the specification of function apply, since this
function simulates native application, we provide it with a postcondition
stating that applying k to arg must be directly related to result. As we will
show, the translation of the post predicate and the generation of the
specification of the apply are indeed amenable to automatic tool translation.

##  4 The GOSPEL to Why3 Pipeline

In spite of the fact defunctionalization can be a bit of a difficult concept
to grasp, once it is understood it becomes clear that we can divide the
process into neatly defined steps and that these steps can very easily be
mechanized. Although we give a complete overview of how OCaml code is
converted in our translation scheme, found in annex 0.E, we will explain in
detail some of the trickier rules and exemplify with some case studies how our
converter fares when confronted with decently complex programs. The full
implementation for our translation scheme, along with some test cases, can be
found in our GitLab
111https://gitlab.com/releaselab/defunctionalization/-/tree/patch-1

###  4.1 Defunctionalization at Work

Explaining in detail every single one of our translation scheme’s minutia
would be a tad redundant, seeing as WhyML and OCaml are very similar, given
they both derive from the ML dialect. Indeed, most of our translations simply
amount to lightly re-organizing the OCaml code to fit WhyML’s syntax, picking
up the free variables and figuring out the type of the expression.
Nevertheless, we will convert a handful of expressions to explain some of the
peculiarities of how we approach defunctionalization before we get into the
relatively more complex case studies.

First we will be looking at the translation of the following expression : x +
y. When converting an OCaml expression we will return, along with its WhyML
translation, its set of free variables and its type. Seeing as this expression
belongs to the class $e1\ o\ e2$, x being $e1$, + being o and y being $e2$, we
will need $e1$ and $e2$’s translation and free variables. Seeing as $e1$ and
$e2$ are simply two variables, we will translating both directly, meaning e1’s
translation will be x and $e2^{\prime}s$ will be y and their set of free
variables will be a singleton set each only containing a tuple with the
variable’s name and its type, which we obtain from a dictionary we refer to in
our translation scheme as varType which binds the name of user defined
variables to their types.

Finally we need only to combine the two translations with o’s, which will
simply be +, producing "x + y", and the two sets of free variables, giving us
{(x, int), (y, int)}. Note however that we do not define this expression’s
type attribute, seeing as it is only used when translating function
applications and we know that an addition can only produce integers, meaning
we will never apply the result of this expression as a function.

It is worth noting that in order to use the addition operator in WhyML we need
to import the int.Int module. To handle this, we import, along with int.Int,
several other libraries so we can naturally use everything in the OCaml
standard library, such as lists and arithmetic operations.

Next, we will bind the previous expression to a lambda and translate it with
the fun arg* : $\tau$ -> e1 construction

fun (x : int) (z : int) : int -> x + y

Important to note that z is not used and y is not defined in the scope of this
expression. First and foremost, we will be adding the function’s arguments’
type to the aforementioned varType dictionary. After computing $e1$’s free
variables, we must remove $x$ and $z$ from the set. Since $e1$ is our last
translated addition (x + y), our new set of free variables is simply
$\\{(y,int)\\}$, since x is defined in the scope of this expression.
Additionally, we will defunctionalize this lambda using the buildApply
function, which takes the lambda’s free variables, type and arguments’ name
and returns its defunctionalized alias. Given the previous defunctionalization
examples, the reader might be expecting something like this

⬇

type kont1 = K int

let apply (k : kont1) (arg1 : int) (arg2 : int) : int =

match k with K y -> let x = arg1 in let z = arg2 in x + y

Although this would be an accurate translation, it is not fully
defunctionalized, due to the fact OCaml technically does not support functions
with more than one argument, which might seem unintuitive seeing as our
original lambda has two. When OCaml parses it, however, it produces a one
argument function that returns another one argument function, such as:

fun (x : int) -> (fun z : int -> x + y)

This means we can partially apply arguments and, instead of throwing an error,
OCaml will simply return a lambda where the applied arguments are now free
variables. Due to this detail, we defunctionalize our functions in the way
OCaml defines them instead of how the user does. Hence, we will first
translate the innermost lambda into K1 y x with the generated
defunctionalization being similar to:

⬇

type kont1 = K1 int int

predicate post1 (k : kont1) (arg :int) (result : int) =

match k with K1 y x -> let z = arg in true

let apply1 (k : kont1) (arg : int) : int =

ensures { post1 k arg result }

match k with

| K1 y x -> let z = arg in y + x

Seeing as we didn’t annotate our original lambda with a postcondition, we fill
in the post predicate with the default value true. Before moving on to the
outermost lambda, we will bind the type int -> int to the function apply1 the
the predicate post1.

Now all that is left is to translate the outermost lambda, which is largely
identical to the first one, except we now bind int -> (int -> int) to the
function apply2 and the predicate post2 and, seeing as the body of this
function will be our previously generated expression, we must append it a
specification that describes its return value so we may use this function in
future specifications. We then conclude that our post predicate is as follows:

⬇

type kont2 = K2 int

predicate post2 (k : kont2) (arg : int) (result : kont1) =

match k with K2 y -> let x = arg in result = K1 y x

In conclusion, we create two apply functions, one that takes the first
argument and returns a continuation and another which accepts that
continuation along with the final argument returning the result of the
original function.

Another detail is that, unlike in our previous examples, we defined two types
(kont1 and kont2) with one constructor each instead of one type with two
constructors. This is due to the fact that we need different apply functions
for each different function type since WhyML does not support generalized
algebraic data types [19] as arguments or return values. We cannot put all our
constructors in the same definition because that would mean, if we had more
than one apply function, their pattern matching wouldn’t be exhaustive and
Why3 wouldn’t be able to dispatch a proof.

Finally, we will set this expression’s type as int -> (int -> int) and reset
the changes made to varType at the start of the translation, that is, remove
the binding of x to int and the binding of z to int. Important to note that
removing a pairing means that, if there was any other value associated with
the key before adding the removed value, that becomes the new value, if not,
the binding now becomes undefined.

To conclude, we will now convert the application of a function to a set of
arguments: (k 3) 4.

First, we will translate the expression within brackets with the $e1\ e2$
rule. To do so, we will get the translation of $e1$ and $e2$, both of which
are directly obtained; their free variables (which are not relevant for this
example), and their type, e2’s is irrelevant, due to the fact it is a literal
value and those are never used as functions, and e1’s is obtained in the same
manner as before: by checking the current mapping for the variable name k. For
the sake of argument, we will assume k is bound to the previous lambda, whose
type is int -> (int -> int).

With both sub expressions now converted, it is time to turn our attention to
our original application, which will be converted as follows:

⬇

apply2 k 3

and it’s type will be int -> int, seeing as we have already applied k’s first
argument.

Next, we do the same thing for the next argument, the only notable difference
is that our expression’s type is now int -> int, which is bound to the integer
1 and thus we will use the apply1 function, giving us:

⬇

apply1 (apply2 k 3) 4

As the one might have noticed, the previous expressions were defunctionalized
even though they weren’t used in a higher order context. This is another quirk
of our translation: any and all function declarations and applications, when
converted, are also defunctionalized, regardless of whether or not they are
used as parameters. The main reason for this is that, when a function is
declared we generally have no way of knowing if it is going to be used as a
first order value. Despite the fact that we are defunctionalizing more than
what we strictly need, translating every function the same way makes the
conversion process more uniform and easier to define and implement.

###  4.2 Extending GOSPEL

Since we have covered the most relevant rules for translating and
defunctionalizing OCaml code, let us take a look at GOSPEL. We won’t go over
its syntax in much detail, seeing as its specifications are quite similar to
Why3’s, although there are some interesting divergences. Namely, GOSPEL
annotations are added as comments beginning with the symbol "@" at the tail
end of our implementations or, for anonymous functions, between the fun
keyword and the function’s arguments like so

⬇

fun [@gospel {| ensures Q |}] args… -> exp

There are some other which are less relevant, for example, having to redefine
the name of each argument and the name for the result and the absence of curly
brackets, but for the most part we will use similar constructions as we did
for Why3 specifications. Due to this, our WhyML conversions will be almost
identical to their GOSPEL counterparts, with one small difference: the post
"predicate". To understand how we have extended GOSPEL to allow usages of
post, let us consider a small example.

⬇

(*@ r = f g x

ensures post (g : int -> int) x r *)

This will be the specification for a function f which takes a function g as a
parameter and an integer x. we have omitted the implementation since we will
be focusing on how we will deal with the post condition. As previously
explained, post is used to avoid directly using non-pure functions (that is,
with side-effects) in our specifications. Because post has to work for any
function applied to however many arguments of any types, it cannot be defined
through WhyML. Instead, as previously shown, we create several posts, one for
each apply function and whenever we encounter a post in our specifications, we
will replace its application with the respective defunctionalized post. That
is why we put predicate in quotation marks; post is not really a predicate nor
a function, but more of a meta-predicate: a predicate which, depending on the
situation, can be any predicate you need. The only restriction for using it is
that the first argument, that is, the function, must have its type explicitly
declared so we know which post to use.

Now that we understand the keyword, let us look at how we convert it. For one
argument functions it is almost direct. For the aforementioned specification,
we would find the post predicate for functions of type int -> int, and replace
the call to post with post2 from our previous translation:

⬇

ensures { post2 g x r }

It gets interesting, however, when we apply post to a function that takes more
than one argument. Let us imagine that the g function was of type int -> (int
-> int) and was applied to an extra argument y, such as, post g x y r. The
translation of such post would be:

⬇

ensures { forall var0 : kont1. post1 g x var0 -> post2 var0 y r }

In simpler terms, if we assume the return value for g when applied to x is
var0, whose type would be kont1, then when var0 is applied to y, it will
return r.

##  5 Case studies

With the most important rules detailed and understood, we will examine how our
tool handles relatively complex programs. The following case studies, unlike
the ones presented in Section 3, are not original and were presented by one of
the project’s supervisors in their precursor paper on defunctionalized
proofs[14].

###  5.1 Height of a tree

Firstly, we will look at at a function that computes the height of a binary
tree with the following signature

⬇

let height_tree_cps (t: int tree): int = …

(*@ r = height_tree t

ensures (height t) = r *)

An interesting peculiarity of this program is that it only accepts a tree of
integers, even though, intuitively, a height function should work for any type
of tree. However, we won’t be dealing with any kind of polymorphism, for that
would add certain complications when building our constructors [15]. Another
thing worth mentioning is that OCaml does not have a binary tree type,
although WhyML does. This means that, even though the OCaml program we will
present technically does not compile on its own unless the programmers defines
a tree exactly as it’s done in WhyML, it will result in a valid program after
the conversion if we import the tree.Tree module. Even so, one could wonder
why we simply don’t define the tree type in our code instead of going in this
roundabout way. The reason is that we are going to need a logical function
from Why3 standard library, namely the height function, in our specification,
which we can only apply in conjunction with WhyML’s tree type. This solution,
flawed as it may be, will be sufficient to prove our program.

First let us define the recursive auxiliary function that performs the
calculation

⬇

let rec height_tree (t: int tree) (k: int -> int) : int

= match t with

| Empty -> k 0

| Node((lt : int tree), _, (rt : int tree)) ->

height_tree lt (fun [@gospel

{| ensures post (k : (integer -> integer))

(1 + max hl (height rt)) result |}]

(hl : int) : int ->

height_tree rt (fun [@gospel

{| ensures post (k : (integer -> integer))

(1 + max hl hr) result |}] (hr : int) : int -> k (1 + max hl hr)))

(*@ r = height_tree t k

ensures post (k : (integer -> integer)) (height t) r*)

This one is a bit trickier than our previous examples, due to the nested
continuations. In short, we create two continuations, one for each tree
branch, which will calculate their heights. Finally, they will compare the two
and return the highest value plus 1. Then, all our main function has to do is
run the recursive auxiliary function with the identity function.

⬇

let height_tree_cps (t: int tree): int =

height_tree t

(fun [@gospel {|ensures result = x|}] (x : int) : int -> x)

(*@ r = height_tree t

ensures (height t) = r*)

The fully annotated program and its conversion can be found in our GitLab and
can be proved using the CVC4 prover, which dispatches the proof in less than a
second.

###  5.2 Small-step interpreter

The final example we will present is an interpreter for a small language that
consists entirely of integer literals and subtractions thereof:

⬇

type exp = Const of integer | Sub of exp * exp

To define how expressions are reduced, we will present two rules, first we
assert that a $Sub$ node with two $Const$ can be reduced to a single $Const$:

$Sub(Const\ v1)(Const\ v2)\xrightarrow{e}Const(v1-v2)$

Additionally, we will also need an inference rule that allows us to properly
parse larger expressions:

$\displaystyle C[e]\xrightarrow{e}C[e^{\prime}]e\xrightarrow{e}e^{\prime}$

C representing the following reduction context:

$\begin{array}[]{crl}C&::=&\square\\\
&|&C[\mathtt{Sub}\;\square\;\mathtt{e}]\\\
&|&C[\mathtt{Sub}\;(\mathtt{Const}\>\mathtt{v})\;\square]\end{array}$

Simply put, $C[e]$ represents a context where the expression $e$ has
completely replaced the symbol $\square$, meaning our inference rule reads: if
an expression $e$ can be reduced to a simpler $e^{\prime}$, then any context
$C[e]$ can be reduced to $C[e^{\prime}]$

With our interpreter clearly described let us move on to how we represent our
contexts. Although the most natural translation would be of an algebraic data
type with three constructors, one for each of our rules, we will represent
them as first order functions that take an expression and return an
expression, leading to much cleaner and compact code.

Next, we will implement a function that quite directly translates our
inference rule by taking a context $C$ and an expression $e$ which returns a
context and a new expression making sure that $C[e]=C^{\prime}[e^{\prime}]$,
$C^{\prime}$ and $e^{\prime}$ being the returned context and expression
respectively.

⬇

let rec decompose_term (e: exp) (c : exp -> exp) :

((exp -> exp) * exp)

= match e with

| Sub (Const (v1 : int), Const (v2 : int)) -> (c, e)

| Sub (Const (v : int), (e : exp)) ->

decompose_term e (fun [@gospel

{| ensures post (c : (exp -> exp)) (Sub (Const v) x) result |}]

(x : exp) : exp -> c (Sub(Const v, x)))

\par| Sub ((e1 : exp), (e2 : exp)) -> decompose_term e1

(fun [@gospel {| ensures post (c : (exp -> exp)) (Sub x e2) result |}]

(x : exp) : exp -> c (Sub(x, e2)))

\par(*@ c_res, e_res = decompose_term e c

requires not (is_value e)

ensures is_redex e_res &&

forall res. post (c : (exp -> exp)) e res ->

post (c_res : (exp -> exp)) e_res res*)

The first thing one might notice is the fact we only match e with Subs and no
Const. When pattern matching is non-exhaustive, WhyML will add an extra branch
of the following type |_ -> absurd and attempt to prove that it is never
reached. In this case, it will meet this obligation knowing that e is not a
value, in other words, not a constant. This knowledge comes from the is_value
predicate, which we do not define above. We haven’t added its implementation,
nor the others we use, seeing as their quite simple. Moreover, functions with
preconditions are handled a bit differently: any function presenting these
will be converted directly instead of defunctionalized, seeing as there are no
interesting case studies of stateless continuations with preconditions.
Therefore we will assume that any function with these will never be used as a
first order value and therefore has no need in being defucntionalized, thereby
simplifying our proof obligations.

Now that we have a function that satisfies $C^{\prime}[e^{\prime}]=C[e]$ we
need one for the condition $C^{\prime}[e^{\prime}]=e$. For this purpose, all
we need is to apply our original function with a context so that $C[e]=e$, in
other words, the identity function.

⬇

let decompose (e: exp) : (exp -> exp) * exp =

decompose_term e

(fun [@gospel {| ensures result = x |}] (x : exp) : exp -> x)

(*@ c_res, e_res = decompose e

requires not (is_value e)

ensures is_redex e_res &&

post (c_res : (exp -> exp)) e_res e *)

With decomposition out of the way, all that’s left is to properly evaluate our
expressions. To do this, we will define a function which takes an expression
and keeps decomposing and performing head reductions, that is, applying the
first rule we presented that reduces a $Sub$ into a $Const$, until it reaches
a constant:

⬇

let red (e : exp) : int =

match e with

| Const (v : int) -> v

| _ -> match decompose e with

((c : exp -> exp), (r : exp)) ->

let r : exp = head_reduction r in red (c r)

(*@ r = red e

ensures r = eval e *)

When converted, this program, whose translation also available in our GitLab,
is once again successfully proved, this time by the Alt-ergo prover, in less
then a second. However to reach it, we need to supply it with the following
lemma, which says that for any two expressions, if they evaluate to the same
value and are applied to the same context, their evaluations will stay equal
relative to one another.

⬇

lemma post_eval:

forall c: exp -> exp, arg1: exp, arg2: exp, r1: exp, r2 : exp.

eval arg1 = eval arg2 -> post0 c arg1 r1 -> post0 c arg2 r2 ->

eval r1 = eval r2

Since this lemma is much broader than our target function’s post condition,
there could be some value in making it part of GOSPEL’s (and, naturally,
Why3’s) standard library, along with other general axioms regarding higher
order functions, so that the programmer has a fairly decent safety net when
their provers are unable to reach the required proof obligations. Of course we
couldn’t put it as is, given that it depends on the eval function, which is
specific to our case study. Nonetheless, we could turn it into a kind of meta-
lemma where we take a logical function of type ’a -> ’a as an argument. To use
such a "meta"-lemma to prove our transitive property would be just a matter of
instantiating it with the proper ’eval’-like function.

Due to the fact that GOSPEL does not yet support lemmas, it will have to be
placed directly in our translation. An interesting limitation for this lemma
is that it could cease to work if we introduced more first order functions of
type exp -> exp that are divorced from the contexts we have described. To
handle this, our tool would have to know which specific continuations the user
is referring to and create a lemma tailored to them. A hypothetical
translation that solves this problem can be found in annex 0.F.

##  6 Limitations

Due to time constraints, we could only hope to produce a rough prototype that
could be used less as a fully fledged tool and more of a proof of concept.
Naturally, what we created presents itself with the following non-exhaustive
list of restrictions, along with those previously mentioned in Section 5:

  * •

To prove an OCaml program, we first needed a way to convert the source code
into a more digestible format, such as an abstract syntax tree. Two choices
emerged; the Typedtree module and the Parsetree module, both from OCaml’s
compiler-libs library; the main difference between the two being that the
Typedtree has already inferred the types for all expressions and the Parsetree
only has the type information that the user puts in the source code. Seeing
as, to build our defunctionalized constructors, we need to know the types of
the free variables, our options were to simply use the Typedtree or,
alternatively, the Parsetree, with the caveat that the user would have to
explicitly add the types for all variables defined in ’let’ bindings, match
cases and function definitions. We decided on using the later seeing as it
would be the easiest option to convert programs, although it does pose a
restriction on the type of programs we can convert.

Due to this choice, our program does not actually type-check the original code
and leaves this to WhyML. This means that, theoretically speaking, we could
have a situation where a badly typed program is converted and the resulting
WhyML program is correctly typed, leading to a nonsensical specification.
Although we believe this situation is impossible, we have no formal proof that
this is the case.

  * •

In order to make the best of the time given, we decided, instead of
translating it in its entirety, to focus our efforts onto a subset of the
OCaml language. Although we are missing some very useful syntactic sugar such
as module declarations and optional arguments, the constructions allowed by
our converter (_e.g._ , match expressions, let..in, if..then..else, lambda
expressions) give us plenty of breathing room to write rich, functional code.

  * •

Due to the fact our apply functions will be mutually recursive, proving
termination, while not impossible [14] ,would be quite difficult, therefore,
we only prove partial correctness; in other words, we do not prove that these
functions are terminating.

  * •

Defunctionalization can only be used knowing a priori which functions will be
used as first order values. Therefore, if a function takes, for example, a
lambda of type int -> int, and no functions of this type are defined, then our
converter throws an error.

  * •

Although we hope to eventually use defunctionalization to prove functions with
side effects, our prototype can only be used to prove pure programs.

##  7 Conclusions and future work

In this report we have proposed a translator that can convert GOSPEL annotated
OCaml programs into WhyML in order to prove their specification, possibly
involving higher order functions. This represents not only the first work
involving the verification of GOSPEL implementations, but also, as far as we
know, the first usage of defunctionalization in the context of function
specification.

While our tool covers an interesting subset of OCaml programs, we have only
scratched the surface of what defunctionalization is capable of. We will save
for future work patching up some of the limitations presented in Section 6,
namely dealing with effectful programs and polymorphism. To tackle the former,
we anticipate the need to by pass the state of the program as an extra
parameter to our posts [14]; for the latter, a much more profound extension to
Why3 is needed, namely extending Why3’s type system with generalized algebraic
data types [15].

## References

  * [1] Bobot, F., Filliâtre, J.C., Marché, C., Paskevich, A.: Why3: Shepherd your herd of provers. In: . (2011)
  * [2] Charguéraud, A., Filliâtre, J., Lourenço, C., Pereira, M.: GOSPEL — Providing OCaml with a Formal Specification Language. In: Formal Methods - The Next 30 Years - Third World Congress. Lecture Notes in Computer Science, vol. 11800, pp. 484–501. Springer (2019), 10.1007/978-3-030-30942-8_29
  * [3] Danvy, O., Millikin, K., Nielsen, L.R.: On One-pass CPS Transformations. Journal of Functional Programming 17(6), 793–812 (2007)
  * [4] Danvy, O., Nielsen, L.R.: Defunctionalization At Work. In: International Conference on Principles and Practice of Declarative Programming (PPDP) (2001). https://doi.org/10.1145/773184.773202
  * [5] Dijkstra, E.W.: Guarded Commands, Nondeterminacy and Formal Derivation of Programs. Commun. ACM 18(8), 453–457 (Aug 1975), 10.1145/360933.360975
  * [6] Dijkstra, E.W.: The Humble Programmer, p. 1972. Association for Computing Machinery, New York, NY, USA (2007), https://doi.org/10.1145/1283920.1283927
  * [7] Filliâtre, J.C.: Deductive Software Verification. International Journal on Software Tools for Technology Transfer (STTT) 13(5), 397–403 (Aug 2011), 10.1007/s10009-011-0211-0
  * [8] Filliâtre, J.C.: One Logic To Use Them All. In: International Conference on Automated Deduction. pp. 1–20. Springer (2013)
  * [9] Filliâtre, J.C., Paskevich, A.: Why3 — Where Programs Meet Provers. In: Programming Languages and Systems. pp. 125–128. Springer Berlin Heidelberg (2013)
  * [10] Iguernelala, M.: Strengthening the Heart of an SMT-Solver: Design and Implementation of Efficient Decision Procedures. Thèse de doctorat, Université Paris-Sud (Jun 2013)
  * [11] Jung, R., Krebbers, R., Jourdan: Iris From The Ground Up: A Modular Foundation For Higher-order Concurrent Separation Logic. Journal of Functional Programming 28 (2018)
  * [12] Kanig, J., Filliâtre, J.C.: Who: A Verifier for Effectful Higher-order Programs. In: Workshop on ML. Edinburgh, Scotland (Aug 2009), 10.1145/1596627.1596634
  * [13] O’Hearn, P.W.: Continuous Reasoning: Scaling the Impact of Formal Methods. In: Symposium on Logic in Computer Science. p. 13–25. IEEE (2018), 10.1145/3209108.3209109
  * [14] Pereira, M.: Desfuncionalizar para Provar. CoRR abs/1905.08368 (2019), http://arxiv.org/abs/1905.08368
  * [15] Pottier, F., Gauthier, N.: Polymorphic typed defunctionalization. SIGPLAN Not. 39(1), 89–98 (Jan 2004). https://doi.org/10.1145/982962.964009, https://doi.org/10.1145/982962.964009
  * [16] Reynolds, J.C.: Definitional Interpreters for Higher-Order Programming Languages. In: Proceedings of the ACM Annual Conference - Volume 2. p. 717–740. ACM ’72 (1972), 10.1145/800194.805852
  * [17] Reynolds, J.C.: Separation logic: A Logic For Shared Mutable Data Structures. In: Symposium on Logic in Computer Science. pp. 55–74. IEEE (2002)
  * [18] Régis-Gianas, Y., Pottier, F.: A Hoare Logic for Call-by-Value Functional Programs. In: International Conference on Mathematics of Program Construction (MPC). vol. 5133, pp. 305–335. Springer (Jul 2008), http://dx.doi.org/10.1007/978-3-540-70594-9_17
  * [19] Xi, H., Chen, C., Chen, G.: Guarded recursive datatype constructors. In: Proceedings of the 30th ACM SIGPLAN-SIGACT symposium on Principles of programming languages. pp. 224–235 (2003)

##  Appendix 0.A Annotated set function

⬇

use array.Array

use int.Int

type arrayList ’a = {

mutable arr : array ’a;

mutable size : int

}

invariant { 0 <= size <= arr.length }

let set (al : arrayList ’a) (i : int) (e : ’a) : ’a

requires { 0 <= i < al.size }

ensures { al.arr[i] = e }

ensures { result = old al.arr[i]}

ensures { al.size = old al.size}

ensures { forall j. 0 <= j < al.size /\ not j = i ->

old al.arr[i] = result.arr[i] }

= let ret = al.arr[i] in

al.arr[i] <- e;

ret

##  Appendix 0.B Defunctionalized reverse function

⬇

type ’a kont = Kid | Krev of ’a * ’a kont

\parlet rec apply kont arg = match kont with

| Kid -> let r = arg in r

| Krev(x, k) -> let r = arg in x::(apply r k)

\parlet rec reverse_aux_defun l k = match l with

| [] -> apply [] k

| x::t -> reverse_aux_defun t (Krev(x, k))

\parlet reverse_defun l = reverse_aux_defun l Kid

##  Appendix 0.C Reverse function stack trace

⬇

reverse_defun [1;2;3]

which then has the following stack trace

⬇

reverse_aux_defun [1;2;3] Kid

\parreverse_aux_defun [2;3] Krev(1, Kid)

\parreverse_aux_defun [3] Krev(2, Krev(1, Kid))

\parreverse_aux_defun [] Krev(3, Krev(2, Krev(1, Kid)))

\parapply [] Krev(3, Krev(2, Krev(1, Kid)))

\par3::(apply [] Krev(2, Krev(1, Kid))

\par3::2::(apply [] Krev(1, Kid))

\par3::2::1::(apply [] Kid)

\par3::2::1::[]

Given this, we can see that this function returns the list 3::2::1::[], that
is, [3;2;1], the reverse of [1;2;3]. Interesting to note that for each call to
apply the first argument, namely the list, doesn’t change. With this in mind,
we could have had the apply function take a void value (called unit in ocaml)
as its first argument and had the Kid function simply return the empty list.
We chose to do it this way, despite this, to keep the reverse_cps function as
close to CPS conventions as possible.

##  Appendix 0.D Defunctionalized length

⬇

use int.Int

use list.List

use list.Length

let rec function apply (k : kont) (arg : int) : int

ensures { post k arg result } =

match k with

| Kid -> let x = arg in x

| Klen kont -> let l = arg in apply kont (1 + l)

end

let rec function length_defun (l : list ’a) (k : kont) : int

ensures { post k (length l) result } =

match l with

| Nil -> apply k 0

| Cons _ t -> length_defun t (Klen k)

end

let length_fin (l : list ’a)

ensures { result = length l } =

length_defun l Kid

##  Appendix 0.E Translation Scheme

Symbol |  Construction |  Translation
---|---|---
e | () |  e.trans := "()";
| s |  e.trans := s; e.type := varType[s].top(); e.freeVar := [(s, e.type)];
| [] |  e.trans := "Nil"
| n |  e.trans := n
| (s)?(e1 (,e2)*) |  if is_constructor() then c.trans := format("%s %s…, s, e1.trans, …) else c.trans := format("(%s, …), e1.trans…) e.freeVar := union(c1.vars,…)
| e1 o e2 |  e.trans := format("%s %s %s", e1.trans, o.trans, e2.trans); e.freeVar := union(e1.freeVar, e2.freeVar)
| e1 :: e2 |  e.trans := format("Cons (%s) (%s)", e1.trans, e2.trans); e.freeVar = union(e1.freeVar, e2.freeVar);
| e1;e2 |  e.trans := format("( %s) ; (%s)", e1.trans, e2.trans); e.freeVar := union(e1.freeVar, e2.freeVar); e.type := e2.type
| var_def in e2 |  e.trans := format("%s in %s", var_def.trans, e2.trans); e.freeVar := union(var_def.freeVar, e2.freeVar.except(var_def.vname)); e.type := e2.type; varType[var_def.vname].pop();"
| match e1 with (|c -> e2)+ |  e.trans := format("match %s with |%s -> %s… end", e1.trans, c.trans, e2.trans, …); e.freeVar := union(e1, e2.except(c.vars),..); e.type := e2.type;
| if e1 then e2 else e3 |  e.trans := format("if %s then %s else %s", e1.trans, e2.trans, e3.trans); e.freeVar := union(e1.freeVar, e2.freeVar, e3.freeVar); e.type := e2.type
| fun arg+ : $\tau$ -> e1 |  for all arg do varType[arg.name].push(arg.type); e.freeVar := e1.freeVar.except(arg1,…); let name = buildApply(e.freeVar, e1.trans, arg1, …, $\tau$.trans) in e.trans := format(%s(%s…,), name, e.freeVar[0]…); e.type := ($\tau$, arg1, …); for all arg do varType[arg.name].pop();
| e1 e2 |  e.trans := format("%s %s %s", getApply(e1.type), e1.trans, e2.trans e.freeVar := union(e1.freeVar, e2.freeVar) e.type := consumeArg(e1.type)
t | type (’s1 (,’s)*)? s2 = (| s3 (of s4(* s)*)?)+ |  t.trans := fomat("type %s %s… = |"%s %s… = %s…", s2, s1,…, arg.trans)
type (’s1 (,’s)*)? s2 = arg(;arg)* |  t.trans := fomat("type %s %s… = |"%s %s… = %s…", s2, s1,…, arg.trans)
type (’s1 (,’s)*)? s2 = $\tau$ |  t.trans := format("type %s %s… = %s, s2, s1,…, $\tau$.trans)
c | n |  c.trans := n
| (s)?(c1 (,c2)*) |  if is_constructor() then c.trans := format("%s %s…, s, c1.trans, …) else c.trans := format("(%s, …), c1.trans…) c.vars := union(c1.vars, c2.vars, …)
| arg |  l.trans := arg.trans; l.vars = {(arg.name)}
| [ ] |  c.trans := "Nil"
| c1::c2 |  l.trans := format("Cons %s %s", c1.trans, c2.trans); c1.vars := union(c1.vars, c2.vars)
o | +(.)? |  o.trans := "+"
| -(.)? |  o.trans := "-"
| *(.)? |  o.trans := "*"
| /(.)? |  o.trans := "/"
| && |  o.trans := "&&"
| || |  o.trans := "||"
| = |  o.trans := "="
$\tau$ | s |  $\tau$.trans := s
| s (-> $\tau$)+ |  $\tau$.trans :=  String.format("%s -> %s, s, $\tau$.trans)
| () |  $\tau$.trans :="()"
arg | () |  arg.trans := "()"
| (s : $\tau$) |  arg.trans := format("s : %s", $\tau$.trans); arg.type := $\tau$.trans; arg.var := s
var_def | let (rec)? s arg* : $\tau$ = e1 |  varDef.vname := s; if isRec() then varType[s].push(($\tau$, arg1,…)); if isFun() then buildApply(s, e1.freeVar, arg1, …, $\tau$.trans); var_def.trans := format("let %s = %s in ", s, e1.trans); if not isRec() then varType[s].push((arg1.type,…, $\tau$));
top_level | t;; |  top_level.trans := t.trans
| var_def;; |  top_level.trans := var_def.trans
| e;; | top_level.trans :=  format("let %s = %s", VAR_NAME, e.trans)

##  Appendix 0.F Defunctionalized lemma

⬇

let rec lemma post_eval (k: kont0) (arg1 arg2 r1 r2: exp)

requires { eval arg1 = eval arg2 }

requires { post0 k arg1 r1 && post0 k arg2 r2 }

ensures { eval r1 = eval r2 }

variant { k }

= match k with

| K2 -> ()

| K1 e c ->

post_eval c (Sub arg1 e) (Sub arg2 e) r1 r2

| K0 n c ->

post_eval c (Sub (Const n) arg1) (Sub (Const n) arg2) r1 r2

|_ -> true

end


11institutetext:  European Southern Observatory, Alonso de Córdova 3107,
Casilla 19001, Vitacura, Santiago 19, Chile
11email<EMAIL_ADDRESS>22institutetext: Department of Astronomy, University
of Florida, Gainesville, FL 32611, USA
33institutetext: Institute of Astronomy, University of Cambridge, Madingley
Road, Cambridge CB3 0HA, UK
44institutetext: National Science Foundation, Division of Astronomical
Sciences, 4201 Wilson Blvd., Suite 1045, Arlington, VA 22230, USA
55institutetext: Gemini Observatory, Northern Operations Center, 670 N. A
ohoku Place, Hilo, HI 96720, USA

# Asymmetric Heating of the HR 4796A Dust Ring
Due to Pericenter Glow

Margaret M. Moerchen 1122 Laura J. Churcher 33 Charles M. Telesco 22 Mark
Wyatt 33 R. Scott Fisher & Christopher Packham 445522

###### Abstract

Context.  We have obtained new resolved images of the well-studied HR 4796A
dust ring at 18 and 25 $\mu$m with the 8-meter Gemini telescopes. These images
confirm the previously observed spatial extent seen in mid-IR, near-IR, and
optical images of the source. We detect brightness and temperature asymmetries
such that dust on the NE side is both brighter and warmer than dust in the SW.
We show that models of so-called pericenter glow account for these
asymmetries, thus both confirming and extending our previous analyses. In this
scenario, the center of the dust ring is offset from the star due to
gravitational perturbations of a body with an eccentric orbit that has induced
a forced eccentricity on the dust particle orbits. Models with 2-$\mu$m
silicate dust particles and a forced eccentricity of 0.06 simultaneously fit
the observations at both wavelengths. We also show that parameters used to
characterize the thermal-emission properties of the disk can also account for
the disk asymmetry observed in shorter-wavelength scattered-light images.

Aims.

Methods.

Results.

###### Key Words.:

circumstellar matter, planetary systems, Stars: individual: HR 4796A,
Infrared: planetary systems

##  1 Introduction

HR 4796A is an A0V star with the highest fractional IR luminosity
($L_{IR}$/$L_{\star}$= 5 x $10^{-3}$) yet discovered among debris disks. Its
recently revised distance estimate is 73 pc (van Leeuwen, 2007) (updated from
67 pc [Perryman et al. 1997]), and it has an M-star T Tauri companion with a
separation of 7.7” (Jura et al., 1993). Its age of 8$\pm$2 Myr (Stauffer et
al., 1995) (based on the lithium abundance of HR 4796B) places the disk of HR
4796A in the very interesting early phase of chaotic disk evolution dominated
by collisions (e.g., Kenyon & Bromley 2005).

Based on the dust temperature, Jura et al. (1993) first tentatively suggested
40 AU (now 44 AU) as the main location of dust responsible for the HR 4796A
excess thermal emission. Indeed, the mid-IR discovery images of the source at
CTIO and Keck revealed a highly inclined ringlike disk with a dust
distribution peaking near 70 AU (now 76 AU) (Jayawardhana et al., 1998;
Koerner et al., 1998). Follow-up coronagraphic near-IR images with
$HST$/NICMOS tightly constrained the width of the ring to $\sim$17 AU (now
$\sim$18.5 AU), with the dust population severely depleted both interior and
exterior to this region (Schneider et al., 1999).

Further imaging observations of the thermal emission confirmed a strong peak
in the dust density at 70 AU (Telesco et al., 2000; Wahhaj et al., 2005).
Telesco et al. (2000) also noted a 1.8-$\sigma$ brightness asymmetry at 18
$\mu$m in which the northeast side of the disk is brighter than the southwest
side. In a companion paper, Wyatt et al. (1999) demonstrated that such an
asymmetry could arise from the phenomenon of pericenter glow. A second body in
the system on an eccentric orbit about the star may cause pericenter glow
when, through secular perturbations, it effectively shifts the center of the
dust ring away from the star and closer to the apastron of the perturbing
companion’s orbit. The side of the ring shifted closer to the star becomes
relatively warmer and more luminous (Dermott et al., 1998), giving rise to the
“glow”. At present, while several works have acknowledged the existence of a
brightness asymmetry (e.g., Wahhaj et al. 2005, Debes et al. 2008), no
alternative explanations for the brightness asymmetry have been proposed.

We have obtained mid-IR images of HR 4796A with the Michelle and T-ReCS
cameras at the Gemini North and Gemini South telescopes, respectively. These
permit us to confirm unambiguously the brightness asymmetry, which was
previously characterized at only the 1.8-$\sigma$ level of significance. By
imaging thermally emitting dust at two wavelengths, we also examine the
spatial distribution of dust temperatures within the resolved disk, and we
conclude that the dust population in the NE ansa of the disk is not only
brighter, but also warmer, than the dust in the SW ansa. With these new data
and with the consideration of shorter-wavelength images of scattered light, we
re-examine the possibility that these observed asymmetric characteristics are
caused by pericenter glow.

##  2 Mid-IR images of HR 4796A

We observed HR 4796A with T-ReCS at Gemini South (using the Qb filter) in May
2004 and with Michelle at Gemini North (using the Qa filter) in April 2005.
The central wavelengths for these two filters are 18.1 $\mu$m for Qa and 24.5
$\mu$m for Qb. A PSF and flux standard star were observed close in time to the
target object for both sets of images. We reduced the images with an IDL
reduction package developed at the University of Florida, and we performed
photometric measurements with the Starlink GAIA program. The images in both
bandpasses show that the disk of HR 4796A is clearly resolved, in the form of
an elongated disk $\sim$3” in extent with two prominent peaks in brightness
(Figure 1).

Figure 1: Top: 18.1-$\mu$m image contours drawn at 2-$\sigma$ intervals
starting at 5$\sigma$ (background $\sigma$ = 0.295 mJy/pixel). Bottom:
24.5-$\mu$m image contours drawn at 1-$\sigma$ intervals starting at
3-$\sigma$ (background $\sigma$ = 1.463 mJy/pixel). Note that the 18.1-$\mu$m
image was taken with Michelle (plate scale: 0.1005”/pixel) and the 24.5-$\mu$m
image was taken with T-ReCS (plate scale: 0.089”/pixel). The images have been
scaled to show the same area in square arcseconds, and both images are
photosphere-subtracted.

Table 1 lists results of the sky-subtracted aperture photometry of HR 4796A at
18.1 $\mu$m and 24.5 $\mu$m. The uncertainties reflect only measurement
errors. Photometric uncertainties are driven primarily by variable sky
transmission throughout the night, and lacking multiple standard star
measurements, we adopt typical photometric uncertainties of 15% for both
filters.

Table 1: Photometry of the HR 4796A Disk

$\lambda_{c}$ | Total Flux Density | Photosphere | Disk Flux Density
---|---|---|---
$[\mu$m] | [mJy] | [mJy] | [mJy]
18.1 | 1106 $\pm$ 7 | 48 | 1058 $\pm$ 7
24.5 | 3307 $\pm$ 47 | 33 | 3274 $\pm$ 47

Notes– Uncertainties listed in this table are measurement uncertainties.
Photometric uncertainties are assumed to be 15% for both bandpasses.

The photospheric contribution in the bandpasses of the imaging observations
are based on a 10.8-$\mu$m value derived by Jura et al. (1998) by
extrapolating the $K$-band (2.2 $\mu$m) flux density using a $\nu^{1.88}$
power law as given by the Kurucz (1979) model atmosphere with T = 9500 K and
log $g$ = 4.0. Photospheric values at longer wavelengths were estimated with
the assumption that HR 4796A and Vega have the same 10-to-18 $\mu$m flux
ratio, since they are both A0V stars. The photospheric contribution to the
total flux density of the source is 4% at 18.3 $\mu$m and 1% at 24.5 $\mu$m.

The PSF was scaled such that the integrated flux matched the photosphere
level, and its center was shifted to the position assumed for the star in the
target source images. The maximum intensities of the scaled PSF and of the
disk (prior to the photosphere subtraction) at the stellar position are,
respectively, 1.4 mJy/pixel and 5.0 mJy/pixel at 18.1 $\mu$m and 0.3 mJy/pixel
and 9.8 mJy/pixel at 24.5 $\mu$m. For the 18.3-$\mu$m image, the adopted
stellar position was the centroid peak of the 2-pixel-smoothed central lobe of
emission (which we assume arises from the presence of the star) between the NE
and SW ansae, where the peaks were also determined from a centroid measurement
of the 2-pixel-smoothed image. The stellar centroid peak position in the
18.1-$\mu$m image differs from the actual midpoint between the NE and SW peaks
of emission by $<$1 pixel, which is insignificant as this is approximately the
level of certainty in the determination of the peak location. One pixel in the
18.1-$\mu$m image corresponds to 0.1”, or $\sim$7 AU. For the 24.5-$\mu$m
image, the stellar position is not obvious, and so the midpoint between the
two centroid peaks of the NE and SW ansae was adopted as the stellar position.
The distance from the NE peak to the SW peak in the 2-pixel-smoothed image is
17.3 pixels (1.74”) at 18.1 $\mu$m and 18.2 pixels (1.62”) at 24.5 $\mu$m.

The scaled PSF image in each bandpass was then subtracted from the
corresponding image of HR 4796A. The total photosphere-subtracted flux density
of the source was measured at each wavelength, and the results are given in
Table 1.

##  3 Disk asymmetry

The global dust color temperature is 102 $\pm$ 11 K, based on the
18.1-$\mu$m/24.5$\mu$m flux density ratio determined from aperture photometry
of the entire photosphere-subtracted disk. This temperature is consistent with
the initial color temperature estimated by Jura et al. (1993) based on $IRAS$
photometry (T$\sim$110 K). The color temperature of the dust as a function of
disk radius was also calculated using the approach described below.

Figure 2: Color temperature profile of the HR 4796A disk based on photometric
measurements at 18.1 $\mu$m and 24.5 $\mu$m (images “cross-convolved” to
achieve the same spatial resolution). The NE profile is represented by
circles, and the SW profile is represented by squares. Combined measurement
uncertainties from both bandpasses are smaller than the size of the data
points. The temperatures at the location of the peak dust density in the ring
(labeled) are the key values to note. Figure 3: Brightness profiles of the HR
4796A observations and model at 18.1 $\mu$m (lower curves) and 24.5 $\mu$m
(upper curves). The brightness profile for each observation or disk model was
generated by summing the central 3” along the long axis of the disk. Model
curves are shown for four different angles of pericenter, and we have
concluded that a pericenter angle of zero provides the best fit.

The position angle (PA), east of north, is defined as the angle between the
vertical image axis and the line connecting the two central peaks of emission
of the two ansae: 28.1∘ in the 18.1-$\mu$m image and 29.9∘ in the 24.5-$\mu$m
image. We estimate an uncertainty for the PA of 0.9 degrees, half the
difference between the values measured in the stacked image at each
wavelength. We adopt a PA value of 29∘, the average of the two values. For
convenient display, the images of HR 4796A were rotated counter-clockwise by
an angle of 61∘, to orient the disk plane parallel to the x-axis.

Quantifying the asymmetric structure– The asymmetric structure is quantified
in both the observations and in the model distributions (§4) by finding the
brightness peak for each ansal “lobe” and its offset from the star. We
calculated the following values for each image:

  * •

$D_{mean}$, the mean offset of the lobes’ peak brightness from the center, in
pixels

  * •

$\frac{\Delta D}{D_{mean}}$, the difference between the distance of the lobe
brightness peaks ($D_{NE}-D_{SW}$), divided by the mean offset of the
brightness peaks from the center, given as a percentage value

  * •

$F_{mean}$, the mean peak brightness of the lobes

  * •

$\frac{\Delta F}{F_{mean}}$, the difference in lobe peak brightness,
($F_{NE}-F_{SW}$), divided by the mean peak brightness of the lobe
($F_{mean}$), given as a percentage value

The level of the brightness asymmetry, $\Delta F/F_{mean}$, is of primary
interest, and it is 15.3$\leavevmode\nobreak\ \pm$ 2.6% at 18.1 $\mu$m and
13.0$\leavevmode\nobreak\ \pm$ 3.8% at 24.5 $\mu$m. The uncertainty in $\Delta
F/F_{mean}$ was calculated by taking a circular symmetrical fit to the
observations (i.e., repeating the modeling but imposing zero forced
eccentricity) and then using Monte Carlo methods to add noise to the model at
the observed level (0.007 mJy/pixel at 18.3 $\mu$m, 0.135 mJy/pixel at 24.5
$\mu$m) to work out the asymmetry that would have been detected purely due to
noise.

Assessing the temperature asymmetry– After removal of the photospheric
contribution, the angular resolution of each image was degraded to achieve the
same resolution in both images for accurate spatial comparison in the color
temperature calculation. Using the $gauss$ routine in IRAF, each of the target
images was convolved with a Gaussian profile having the same FWHM as the PSF
from the other bandpass. The PSF FWHM was estimated by a Gaussian fit to the
azimuthally averaged profile. The FWHM values were 5.19 pixels (0.52”) in the
18.1-$\mu$m image and 8.08 pixels (0.72”) in the 24.5-$\mu$m image. The PSF
images themselves were convolved in the same manner, which confirmed that the
resulting resolution (assessed with FWHM measurements) was identical for both
images. The final resolution of the cross-convolved images is $\sim$0.9”,
which corresponds to $\sim$65 AU for the source distance of 73 pc.

The rotated images were cropped to a swath of vertical width
$\Delta$y$\sim$3”. The pixel values in these images were summed along the
y-axis, resulting in a one-dimensional brightness profile for the disk at each
wavelength. From these brightness profiles, a color profile and a color
temperature profile (Figure 2) were constructed to illustrate variations in
temperature along the extent of the disk. As a test of the robustness of the
temperature estimates, the 24.5-$\mu$m brightness profile was shifted by
$\pm$1 pixel ($\sim$6 AU) along the x-axis relative to the 18.1-$\mu$m
brightness profile before the two brightness profiles were combined to
determine the color profile. The two resulting color profiles that
incorporated these offsets yielded color temperature profiles that did not
differ from the results presented in Figure 2 within 60 AU of the star by more
than 2 K.

The NE and SW profiles are overplotted on the same radius scale (along the
x-axis) to better show differences in temperature at the same distances from
the central star. We see no significant variation in color temperature within
a radius of 15 AU from the star. While that result is expected for such a
narrow annulus, this result is in fact attributable principally to the
degradation of resolution that we carried out to make the temperature
estimates. This effect is also seen in the separation of the two ansae peaks,
which decreased for both images following the cross-convolution, from 1.74” to
1.1” at 18.1 $\mu$m, and from 1.62” to 1.58” at 24.5 $\mu$m (where the lesser
effect for the 24.5 $\mu$m image is likely due to being convolved by a smaller
kernel).

There is a clear temperature asymmetry, with the dust in the NE side of the
disk having significantly higher color temperatures at greater than 30 AU from
the star and in particular at the location of the emission peaks corresponding
to the two ansae of the dust ring. The temperature difference peaks at $\sim$9
K near 60 AU. The uncertainties due to measurement error in the temperature
estimates are less than 2 K. The color temperature uncertainties that
incorporate photometric uncertainties are $\pm$8 K, but in this analysis we
are most interested in the relative temperatures of the two ansae and are
therefore most concerned with uncertainties within the images and not absolute
photometric uncertainties.

##  4 Pericenter glow modeling

###  4.1 Pericenter glow scenario and prior work

One motivation for observing debris disks is the potential to reveal as-yet
unseen planetary companions through the disk morphology. A massive orbiting
body can gravitationally perturb the orbits of smaller bodies like dust
particles over long-period (“secular”) timescales, and these perturbations may
be observed in the global characteristics of the disk. For example, a planet
on an eccentric orbit eventually imposes a forced eccentricity on the dust
particles, which results in an offset of the disk’s center of symmetry from
the host star. While this offset may be observable directly, other
manifestations of the offset may be more prominent. In particular, if the dust
population is azimuthally homogenous in size distribution and composition, the
side of the disk (the pericenter side) that is offset toward the star will
experience enhanced heating, increasing both its temperature and brightness.
The resulting asymmetry in brightness between the apocenter and pericenter
sides of the disk is referred to as pericenter glow (Dermott et al., 1998).

The asymmetric morphology of the HR 4796A disk observed in previous thermal
emission images at 18 $\mu$m was well approximated by models that invoked
pericenter glow (Wyatt et al., 1999; Telesco et al., 2000), but the
statistical significance of the asymmetry in those earlier images was
relatively low. We note that the brightness asymmetry was also previously
detected at 10.8 $\mu$m (Telesco et al., 2000), but the much greater
prominence of the starlight at those wavelengths makes such shorter wavelength
observations less useful for our analysis. With new 18.1-$\mu$m images and the
addition of images at 24.5 $\mu$m in this dataset, we have re-examined this
hypothesis by using the same models to reconstruct the observed thermal
emission and color temperature profiles.

###  4.2 Thermal emission reconstruction

Model inputs– The physics of pericenter glow and the associated modeling of it
are outlined in detail by Wyatt et al. (1999), and here we review only the key
parameters involved in that analysis. The model for the disk’s density
distribution of dust was generated with three distinct inputs: the physical
structure of the disk, the combination of the optical properties and size
distribution of the dust particle population, and the orientation of the disk
along the line of sight. We recall that the mass of the perturbing body itself
does not come into play and is not constrained by this approach. The physical
structure of the disk is modelled by a radial density distribution of dust
with defining parameters as follows: $a_{min}$ and $a_{max}$, the inner and
outer semimajor axes of the dust ring (which are roughly equivalent to the
inner and outer radii), $\gamma$, the power-law index of the semimajor axis
distribution, $e_{f}$, the forced eccentricity, and $\sigma_{tot}$, the total
surface area of the dust particles in the disk. The forced eccentricity is
that imposed on all particles in the disk by the planet’s secular
perturbations (Murray & Dermott, 1999); as such, the planet’s orbit is not
affected by its own perturbations. The eccentricities of particle orbits are
the vector sum of this forced eccentricity and their own proper eccentricity.
To first order, the disk remains circular but is just offset. The total
surface area of the dust is effectively a scaling factor that can be adjusted
to approximate the overall brightness level of the dust emission. While the
intrinsic disk offset and asymmetry (i.e., if viewed face-on) are determined
by the forced eccentricity $e_{f}$, the measured offset and asymmetry are
determined by the additional variables of disk inclination, $i$, and the angle
of pericenter, ${\omega}_{f}$. For example, the forced eccentricity required
to replicate a given measured asymmetry increases as the angle of pericenter
deviates from being perpendicular to the line of sight.

The optical properties of dust were calculated using Mie theory, Rayleigh-Gans
theory, and geometric optics (Li & Greenberg, 1997; Augereau et al., 1999;
Bohren & Huffman, 1983) for astronomical silicate spheres the diameter of
which ($D_{typ}$=2 $\mu$m) was constrained by a simultaneous fit to the images
at both wavelengths. The resulting dielectric constant was 6.7x$10^{-5}$, and
the absorption coefficients were 0.54 and 0.30 at 18.1 and 24.5 $\mu$m,
respectively. While a single particle size is used to represent a particle
population that has a range of sizes and composition, it reasonably
encapsulates key mid-IR emission characteristics of the disk. We adopted
values for the stellar luminosity and effective temperature of 21 L⊙ and 9500
K, respectively (e.g., Debes et al., 2008).

The output from the disk radial distribution model and the calculation of the
dust optical properties were combined with the specified orientation to
calculate the line-of-sight brightness and generate images of the model disk.
The model images at each wavelength were rotated to the same position angle as
the observed disk, and each image was convolved by a PSF kernel with the
profile width of the observed PSF reference star in the same filter. In the
generation of color temperature profiles for the model data, the images were
also convolved by the profile width corresponding to the other filter to match
the treatment of the observational data (§3).

Figure 4: Color temperature profiles of the HR 4796A observations and models.
Colors were calculated by summing the flux within a 3”-wide swath along the
long axis of the disk for each image. Model curves are shown for two different
angles of pericenter, and we have concluded that a pericenter angle of zero
provides the best fit.

Optimizing the model– With the images of the model disk, we generated
brightness profiles and color temperature profiles in the same way that we did
with the observational data (§3). We compared the observed and model profiles,
and we adjusted the model parameters until the model disk yielded a good
approximation to the observations. First, we varied the inner and outer disk
radii and the total disk surface area until the brightness profile width and
magnitude were well approximated. We then adjusted the forced eccentricity
until the brightness asymmetry level was matched for a fixed value of the
angle of pericenter. The fit was repeated for several values of the angle of
pericenter to assess the corresponding forced eccentricity necessary to
reproduce the same asymmetry (Table 2). The forced eccentricity and the angle
of pericenter define the magnitude of the radial offset asymmetry in the disk
model, which in turn defines the brightness asymmetry and the dust color
temperature asymmetry. The initial assessment of how changing these parameters
affected the resulting brightness profiles was performed by chi-squared
minimization on the brightness profiles and images. The brightness profiles of
the observations and the model disk with varying values for the angle of
pericenter are shown in Figure 3, and the color temperature profiles are shown
in Figure 4. The quantitative measurements used to assess the asymmetry in the
model and in the observed disk, as described above, are summarized in Table 2.

Table 2: Asymmetry parameters for observations and models

|  | 18.1 |  | 24.5 |
---|---|---|---|---|---
| $e_{f}$ |  $\frac{\Delta D}{D_{mean}}$ [%] |  $\frac{\Delta F}{F_{mean}}$ [%] |  |  $\frac{\Delta D}{D_{mean}}$ [%] |  $\frac{\Delta F}{F_{mean}}$ [%] | reduced $\chi^{2}$
observed | – | 0 | 15.28 |  | 0 | 13.00 | –
model $\omega_{f}=0^{\circ}$ | 0.06 | -11.71 | 15.26 |  | -6.11 | 13.10 | 2.17
model $\omega_{f}=40^{\circ}$ | 0.10 | -10.81 | 15.31 |  | -7.01 | 13.30 | 2.31
model $\omega_{f}=55^{\circ}$ | 0.13 | -11.18 | 15.35 |  | -7.15 | 13.05 | 2.21
model $\omega_{f}=75^{\circ}$ | 0.30 | -14.50 | 10.55 |  | -7.39 | 12.70 | 2.26
111 All models were simultaneously fitted to 18.1- and 24.5-$\mu$m brightness
profile linecuts of 0.5” and share the following parameter values:
$a_{min}$=70 AU, $a_{max}$=84 AU, $\gamma$=-1.5, $D_{typ}$=2 $\mu$m, and
$i$=14.1∘. Also recall that the value $\frac{\Delta D}{D_{mean}}$ for the
observations is 0 because we defined the stellar position at 24.5 $\mu$m as
the midpoint between the brightness peaks of the two ansae (with the measured
difference between the midpoint and the stellar peak at 18.1 $\mu$m to be $<$1
pixel.)

Results of the modelling– We have determined that a forced eccentricity of
0.06 imposed on a ring of dust with $D_{typ}$= 2 $\mu$m spanning 74–84 AU can
reproduce the thermal IR images, fitting the observed brightness asymmetry
well at both 18.1 and 24.5 $\mu$m. The disk density power-law exponent was set
to -1.5. The same result for the characteristic particle size was inferred
previously in order to simultaneously fit 10-$\mu$m and 18-$\mu$m images of HR
4796A (Wyatt et al., 1999; Telesco et al., 2000). We have estimated the
uncertainty in the forced eccentricity to be 0.01, which is the change
required to effect a 1-$\sigma$ change in the brightness asymmetry of the
model at 18.3 $\mu$m (the more limiting case). The angle of pericenter was
0$\pm$30∘ (perpendicular to the line of sight), and the disk inclination was
14.1∘ (Schneider et al., 2009). We note that these results differ from those
found by Wyatt et al. (1999) and Telesco et al. (2000) ($e_{f}$=0.02,
$\omega_{f}$=$75^{\circ}$) likely because of the different levels of observed
asymmetry, which are nonetheless within 3$\sigma$ of one another (5.1 $\pm$
3.2% in the previous work, 15.3 $\pm$ 2.6% in this work).

We have also used the model to quantify the limits set by these data on the
presence of hot dust inside the ring. We compared the photometry of the
photosphere-subtracted images with that of the model within an aperture radius
equal to the size of the PSF FWHM (0.52” at 18.3 $\mu$m, 0.72” at 24.5
$\mu$m), and then we repeated the measurement with the same aperture on the
models, where an unresolved flux component could be added at the position of
the star. Excluding calibration uncertainties, the measured flux density
within a 0.52”-radius aperture on the 18.3-$\mu$m photosphere-subtracted image
was 289 $\pm$ 3 mJy. The same measurement on the model with no added
unresolved flux component was 290 mJy. Therefore, the maximum unresolved
contribution we can add to the model that is 3-$\sigma$ consistent with the
observations is 8 mJy. Repeating the same procedure at 24.5 um with an
aperture of radius 0.72”, we measured a flux density of 1290 $\pm$ 12 mJy for
the photosphere-subtracted images and 1289 mJy for the model. Therefore, the
3-$\sigma$ upper limit of 1326 mJy for this aperture allows the presence of an
unresolved component at a level of 37 mJy. However, these measurements do not
account for the assumed 15% calibration uncertainty. Since a change in
calibration would also result in a corresponding change to the model (to fit
the peaks at the correct level), it is instructive to consider what would have
been derived with a revised calibration factor. In the extreme situation that
the calibration factor was 45% higher, the flux measured in the previously
described apertures at 18.1 and 24.5 $\mu$m respectively would have been 438
$\pm$ 4 mJy and 1881 $\pm$ 17 mJy, with corresponding model fluxes of 290*1.45
= 421 mJy and 1289*1.45 = 1869 mJy, and so a flux excess of 17 $\pm$ 4 and 12
$\pm$ 17 mJy. Therefore, a 3-$\sigma$ deviation in both calibration and
statistical uncertainty puts the limits on the unresolved flux at $<$29 mJy
and $<$63 mJy at 18.1 and 24.5 $\mu$m. Wahhaj et al. (2005) estimated an
exozodiacal contribution of 87 $\pm$ 22 mJy at 24.5 $\mu$m, which is
consistent to within 3$\sigma$ with the upper limit we derive.

###  4.3 Scattered light comparison

A brightness asymmetry has also been detected at wavelengths that are
dominated by light scattered by the dust disk. Images obtained by Schneider et
al. (1999) at 1.1 $\mu$m were suggestive of a brightness imbalance favoring
the NE that was later noted to be at the 10–15% level (Schneider et al.,
2001). The disk ansae were again compared in the analysis of data sets taken
with $HST$ NICMOS at 1.71–2.22 $\mu$m, but the NE ansa was not consistently
brighter in all bandpasses (Debes et al., 2008). Most recently, $HST$ STIS
images at 0.57 $\mu$m show that the SW ansa is 0.74 $\pm$ 0.07 times as bright
as the NE ansa, as measured in “lobes” within $\pm$20∘ of the major axis and
between 0.8” and 1.5” from the star. Schneider et al. (2009) measured an
offset asymmetry, $(D_{NE}-D_{SW})/D_{mean}\leavevmode\nobreak\ =-3.7\%$, such
that the center of the disk is offset from the star by 1.4 $\pm$ 0.4 AU (19
$\pm$ 6 mas) in the plane perpendicular to the line of sight.

Thus, another test for our pericenter glow model is whether it can replicate
the observations in scattered light. To examine this, we modified the model’s
particle properties to correctly calculate the scattered light fluxes, using a
combination of Mie theory, the Rayleigh-Gans approximation, and geometric
optics, which yielded a mean albedo of 0.05. We assumed that the same particle
population (in this case, modeled as a first approximation by 2-$\mu$m
silicate spheres) is responsible for both the thermal emission and the
scattered light reflection. We adopted a Henyey-Greenstein phase function to
model the non-isotropic scattering properties of the small grains (Augereau et
al., 1999), with the form

$f(\alpha)=\frac{1-g^{2}}{1-g^{2}-2gcos(\alpha)}$ (1)

where $\alpha$ is the scattering angle to the line of sight and $g$ is the
asymmetry parameter. The parameter $g$ indicates the relative level of front
or back scattering, such that $g=1$ for 100% forward scattering, $g=-1$ for
100% backward scattering, and $g=0$ for isotropic scattering. We adopted the
value $g=0.16$ as determined from the optical $STIS$ images (Debes et al.,
2008; Schneider et al., 2009), which is consistent within 2$\sigma$ of the
values for $g$ determined by Debes et al. (2008) for the near-IR $NICMOS$
images.

To check that our model produced realistic flux levels, we compared them to
those measured by Schneider et al. (2009) at a wavelength of 0.57 $\mu$m. Due
to the presence of coronagraphic spikes in the image, the flux density
measurement had to be corrected for the missing area, yielding a total flux
9.4$\pm$0.8 mJy. The total flux in our model scattered light image is 10.8
mJy. As a secondary check, we used the formula given by Weinberger et al.
(1999) to predict the surface brightness of a disk in scattered light (in mJy
arcsec-2)

$S=\frac{F\tau\omega}{4\pi\phi^{2}}$ (2)

where $F$ is the flux received from the star in mJy, $\tau$ is the optical
depth of the scattering material, $\omega$ is the albedo (from Mie theory,
where $Q_{sca}$=0.08), and the factor $\tau\omega$ can be approximated by
$L_{disk}/L_{\star}$ if $\omega=1$. The expected flux density calculated with
this relation agrees with the flux density values both from observations and
from our model to within 20% at 0.57, 1.1, and 1.6 $\mu$m.

We also replicated the lobe asymmetry measurement by Schneider et al. (2009)
with our disk model. Within the same apertures (within $\pm$20∘ of the major
axis and between 0.8” and 1.5” from the star), we determine a SW/NE brightness
ratio of 0.92 at 0.57 $\mu$m and 0.93 at 1.1 $\mu$m. Our model values are
therefore consistent within 3$\sigma$ with the Schneider et al. (2009) value
of 0.74$\leavevmode\nobreak\ \pm\leavevmode\nobreak\ $0.07, and are also
consistent with a simple 1/$r^{2}$ attenuation of starlight due to the disk
ansae being at different distances from the star. Our model also replicates
the radial offset asymmetry at a level of 3.3%, which is consistent with the
Schneider et al. (2009) value of 3.7 $\pm$ 1.1%.

Schneider et al. (2009) acknowledge that dynamical perturbations such as those
we propose may be responsible for part of the brightness asymmetry, also
suggesting that a segregation of the particle population or an azimuthal
density variation in a homogenous dust ring could produce the asymmetry.
However, whether such a model could also reproduce the observed temperature
asymmetry has yet to be tested. Additionally, we note that particles in the
disk should be moving slightly faster at pericenter, relative to their
velocity at apocenter, and the collisional timescale at pericenter should be
relatively shorter in turn. Therefore, the population of small particles that
scatter near-IR light may be enhanced due to the increased collision rate near
pericenter, thus yielding a density variation in the dust ring and a greater
brightness imbalance in scattered light. There is the potential to further
develop the pericenter glow model with such an enhanced particle population,
which may help to determine whether factors beyond dynamical perturbations
must be invoked to explain the observed morphology. Such refinements to our
model would bear greater consideration especially if the large near-IR
brightness asymmetry is confirmed.

##  5 Discussion and conclusions

Using new mid-IR images at 18 and 25 $\mu$m, we have confirmed the previously
observed brightness asymmetry discovered by Telesco et al. (2000) in the HR
4796A disk to a high level of statistical significance. We have determined
that a model disk with forced eccentricity of 0.06 that is dynamically induced
by a companion can replicate that brightness asymmetry. After incorporating
the dust scattering properties, we can also reproduce, to within 3$\sigma$,
the level of brightness asymmetry seen in near-IR scattered-light images . We
believe that near-IR imaging in additional bandpasses will further constrain
the level of the asymmetry observed in scattered light, and that such
constraints will assist in confirming whether the pericenter glow model is
sufficient to explain the disk asymmetry observed throughout the spectrum.

The model disk is comprised of particles 2 $\mu$m in diameter ($D_{typ}$), in
agreement with previous analysis of the disk (Wyatt et al., 1999; Telesco et
al., 2000), and the particles are assumed to be spherical silicates. Dust
particles smaller than $\sim$8 $\mu$m in diameter should be blown out roughly
within an orbital timescale ($\sim$400 yr). In reality the dust population is
not this simple, but we have adopted this value to characterize the overall
population. We must consider that the particles may have a more complex set of
optical properties that allows them to be heated to the temperatures we infer
and to remain in the system without being rapidly blown out. More complex dust
properties (e.g., Li & Lunine 2003) could be incorporated into the pericenter
glow model in future iterations.

Finally, we recall that a key result for a system like HR 4796A that exhibits
pericenter glow is the implied presence of a perturbing companion. In fact,
the principal perturbing body may be the nearby M-star companion HR 4796B.
However, the orbital parameters of HR 4796B are unknown and may not be
commensurate with its being the primary perturber. The sharp inner edge of the
ring points to the possibility that the perturbing object responsible for the
asymmetry is also truncating the inner edge, as has been suggested for
Fomalhaut (Kalas et al., 2005; Quillen, 2006). In this case, the putative
planet responsible would be expected to lie just interior to the ring. For
example, a Jupiter-mass planet would be expected to orbit at $\sim$60 AU with
an eccentricity of $\sim$0.06 (maximum projected separation 0.87”), but it
could be closer to the dust ring if it were more massive. Given the ambiguity
in potential perturbers and the fact that the level of forced eccentricity
induced by the companion does not constrain its mass, the continued search for
the putative companion is warranted, if not with the current generation of
telescopes then with the next, such as the space-based $JWST$ or the ground-
based TMT or E-ELT facilities. Kalas et al. (2008) have demonstrated a
successful direct planet detection that was pursued following similar evidence
of a forced stellar offset, as in the case of Fomalhaut, which is nearer and
brighter than HR 4796A.

###### Acknowledgements.

MMM gratefully acknowledges fellowship support from the Michelson Science
Center. This work was performed [in part] under contract with JPL funded by
NASA through the Michelson Fellowship Program. JPL is managed for NASA by the
California Institute of Technology. LJC is grateful for the support of an STFC
studentship. Observations were obtained at the Gemini Observatory, operated by
AURA, Inc., under agreement with the NSF on behalf of the Gemini partnership:
NSF (US), PPARC (UK), NRC (Canada), CONICYT (Chile), ARC (Australia), CNPq
(Brazil), and CONICET (Argentina).

## References

  * Augereau et al. (1999) Augereau, J. C., Lagrange, A. M., Mouillet, D., Papaloizou, J. C. B., & Grorod, P. A. 1999, A&A, 348, 557
  * Bohren & Huffman (1983) Bohren, C. F. & Huffman, D. R. 1983, Absorption and scattering of light by small particles, ed. Bohren, C. F. & Huffman, D. R.
  * Debes et al. (2008) Debes, J. H., Weinberger, A. J., & Schneider, G. 2008, ApJ, 673, L191
  * Dermott et al. (1998) Dermott, S. F., Grogan, K., Holmes, E. K., & Wyatt, M. C. 1998, in Exozodiacal Dust Workshop:, ed. D. E. Backman, L. J. Caroff, S. A. Sandford, & D. H. Wooden, 59–+
  * Jayawardhana et al. (1998) Jayawardhana, R., Fisher, S., Hartmann, L., Telesco, C., Pina, R., & Fazio, G. 1998, ApJ, 503, L79+
  * Jura et al. (1998) Jura, M., Malkan, M., White, R., Telesco, C., Pina, R., & Fisher, R. S. 1998, ApJ, 505, 897
  * Jura et al. (1993) Jura, M., Zuckerman, B., Becklin, E. E., & Smith, R. C. 1993, ApJ, 418, L37+
  * Kalas et al. (2008) Kalas, P., Graham, J. R., Chiang, E., Fitzgerald, M. P., Clampin, M., Kite, E. S., Stapelfeldt, K., Marois, C., & Krist, J. 2008, Science, 322, 1345
  * Kalas et al. (2005) Kalas, P., Graham, J. R., & Clampin, M. 2005, Nature, 435, 1067
  * Kenyon & Bromley (2005) Kenyon, S. J. & Bromley, B. C. 2005, AJ, 130, 269
  * Koerner et al. (1998) Koerner, D. W., Ressler, M. E., Werner, M. W., & Backman, D. E. 1998, ApJ, 503, L83+
  * Kurucz (1979) Kurucz, R. L. 1979, ApJS, 40, 1
  * Li & Greenberg (1997) Li, A. & Greenberg, J. M. 1997, A&A, 323, 566
  * Li & Lunine (2003) Li, A. & Lunine, J. I. 2003, ApJ, 590, 368
  * Murray & Dermott (1999) Murray, C. D. & Dermott, S. F. 1999, Solar system dynamics
  * Perryman et al. (1997) Perryman, M. A. C., Lindegren, L., Kovalevsky, J., Hoeg, E., Bastian, U., Bernacca, P. L., Crézé, M., Donati, F., Grenon, M., van Leeuwen, F., van der Marel, H., Mignard, F., Murray, C. A., Le Poole, R. S., Schrijver, H., Turon, C., Arenou, F., Froeschlé, M., & Petersen, C. S. 1997, A&A, 323, L49
  * Quillen (2006) Quillen, A. C. 2006, MNRAS, 372, L14
  * Schneider et al. (2001) Schneider, G., Nicmos/Eons, T., & Go/8624 Teams. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 244, Young Stars Near Earth: Progress and Prospects, ed. R. Jayawardhana & T. Greene, 203–+
  * Schneider et al. (1999) Schneider, G., Smith, B. A., Becklin, E. E., Koerner, D. W., Meier, R., Hines, D. C., Lowrance, P. J., Terrile, R. J., Thompson, R. I., & Rieke, M. 1999, ApJ, 513, L127
  * Schneider et al. (2009) Schneider, G., Weinberger, A. J., Becklin, E. E., Debes, J. H., & Smith, B. A. 2009, AJ, 137, 53
  * Stauffer et al. (1995) Stauffer, J. R., Hartmann, L. W., & Barrado y Navascues, D. 1995, ApJ, 454, 910
  * Telesco et al. (2000) Telesco, C. M., Fisher, R. S., Piña, R. K., Knacke, R. F., Dermott, S. F., Wyatt, M. C., Grogan, K., Holmes, E. K., Ghez, A. M., Prato, L., Hartmann, L. W., & Jayawardhana, R. 2000, ApJ, 530, 329
  * van Leeuwen (2007) van Leeuwen, F., ed. 2007, Astrophysics and Space Science Library, Vol. 350, Hipparcos, the New Reduction of the Raw Data
  * Wahhaj et al. (2005) Wahhaj, Z., Koerner, D. W., Backman, D. E., Werner, M. W., Serabyn, E., Ressler, M. E., & Lis, D. C. 2005, ApJ, 618, 385
  * Weinberger et al. (1999) Weinberger, A. J., Becklin, E. E., Schneider, G., Smith, B. A., Lowrance, P. J., Silverstone, M. D., Zuckerman, B., & Terrile, R. J. 1999, ApJ, 525, L53
  * Wyatt et al. (1999) Wyatt, M. C., Dermott, S. F., Telesco, C. M., Fisher, R. S., Grogan, K., Holmes, E. K., & Piña, R. K. 1999, ApJ, 527, 918