# Observing planet-disk interaction in debris disks

S. Ertel, S. Wolf, J. Rodmann

► To cite this version:

S. Ertel, S. Wolf, J. Rodmann. Observing planet-disk interaction in debris disks. Astronomy and Astrophysics - A&A, 2012, 544, 10.1051/0004-6361/201219236. insu-03612414

HAL Id: insu-03612414

https://insu.hal.science/insu-03612414

Submitted on 18 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Observing planet-disk interaction in debris disks

S. Ertel<fnref n="1" />,<fnref n="2" />, S. Wolf<fnref n="2" />, and J. Rodmann<fnref n="3" />,<fnref n="4" />

Received 17 March 2012 / Accepted 30 May 2012

## ABSTRACT

**Context.** Debris disks are commonly considered to be a by-product of planet formation. Structures in debris disks induced by planet-disk interaction are promising to provide valuable constraints on the existence and properties of embedded planets.

**Aims.** We investigate the observability of structures in debris disks induced by planet-disk interaction with future facilities in a systematic way. High-sensitivity, high angular resolution observations with large (sub-)mm interferometers and large space-based telescopes operating in the near- to mid-infrared wavelength range are considered.

**Methods.** The observability of debris disks with the Atacama Large Millimeter/submillimeter Array (ALMA) is studied on the basis of a simple analytical disk model. Furthermore, *N*-body simulations are used to model the spatial dust distribution in debris disks under the influence of planet-disk interaction. From these simulations, images at optical scattered light to millimeter thermal re-emission are computed. Available information about the expected capabilities of ALMA and the *James Webb* Space Telescope (JWST) are used to investigate the observability of characteristic disk structures with these facilities through spatially resolved imaging.

**Results.** Our simulations show that planet-disk interaction can result in prominent structures in the whole considered wavelength range. The exact result depends on the configuration of the planet-disk system and on the observing wavelength which provides the opportunity of detecting and characterizing extrasolar planets in a range of masses and radial distances from the star that is not accessible to other techniques. Facilities that will be available in the near future at both considered wavelength ranges are shown to provide the capabilities to spatially resolve and characterize structures in debris disks that arise because of planet-disk interaction. Limitations are revealed and suggestions for possible instrument setups and observing strategies are given. In particular, ALMA is limited by its sensitivity to surface brightness, which requires a trade-off between sensitivity and spatial resolution. Space-based mid-infrared observations will be able to detect and spatially resolve regions in debris disks even at a distance of several tens of AU from the star, where the emission from debris disks in this wavelength range is expected to be low.

**Conclusions.** Both ALMA and the planned space-based near- to mid-infrared telescopes will provide unprecedented capabilities to study planet-disk interaction in debris disks. In particular, a combination of observations at both wavelengths will provide very strong constraints on the planetary/planetary systems.

**Key words.** techniques: high angular resolution – techniques: interferometric – planets and satellites: detection – planet-disk interactions – infrared: planetary systems – submillimeter: planetary systems

## 1. Introduction

The dust detected in debris disks is thought to be removed from those systems by the stellar radiation on time scales that are short compared to their ages. This means that the dust must be transient, or more likely continuously replenished by ongoing collisions of bigger objects such as planetesimals left over from the planet formation process (for a recent review see, e.g., Krivov 2010). The presence of planets and debris disks is hence thought to be correlated. In a system with a debris disk and one or more planets, one would expect gravitational interaction between the dust grains and the planet, trapping them into resonance (Wyatt 2006; Wolf et al. 2007; Stark & Kuchner 2008, 2009). This results in structures in the disk that may be observable, which in turn provides a method to infer and characterize planets that are in a regime of masses, brightnesses, and radial distances from the star that is not accessible through other techniques such as radial velocity measurements or direct imaging (e.g., Udry & Mayor 2008; Marois et al. 2008; Kalas et al. 2008).

Clumpy structures in debris disks have been observed in several cases (e.g., *ε Eri*, Greaves et al. 1998; AU Mic, Liu 2004; HD 107146, Corder et al. 2009; Hughes et al. 2011). In this

work, we investigate the observability of structures in debris disks in a systematic way. We set up several initial conditions for the planetary mass and orbit as well as for the dust distribution in the disk from which we simulate the spatial dust distribution using *N*-body simulations. Our results from the dynamical simulations are consistent with earlier works (e.g., Holmes et al. 2003; Wyatt 2006; Stark & Kuchner 2008). In contrast to these works which focused mostly on the phenomenological and theoretical description of the structures as well as on the processes of planet-disk interaction, the goals of the present work are the evaluation of the observability and the development of strategies for observations of planet disk interaction. Thus, we use the *N*-body simulations as a tool to produce realistic structures in debris disks, but discuss these structures only briefly emphasizing on the consequences of different planet-disk configurations on the observability. To make predictions on the feasibility to detect and spatially resolve characteristic structures, available information on the capabilities of the Atacama Large Millimeter/submillimeter Array (ALMA) and *James Webb* Space Telescope (JWST) are used representative for near-future facilities.

We describe the approach used in our *N*-body simulations and for the image generation in Sect. 2. The results from the

<footnote n="1">¹ UJF Grenoble 1/CNRS-INSU, Institut de Planétologie et d'Astrophysique de Grenoble (IPAG) UMR 5274, 38041 Grenoble, France
e-mail: steve.ertel@obf.ujf-grenoble.fr</footnote>

<footnote n="2">² Christian-Albrechts-Universität zu Kiel, Institut für Theoretische Physik und Astrophysik, Leibnizstraße 15, 24098 Kiel, Germany</footnote>

<footnote n="3">³ European Space Agency, Space Environment and Effects Section, Keplerlaan 1, PO Box 299, 2200 AG Noordwijk, The Netherlands</footnote>

<footnote n="4">⁴ Georg-August-Universität Göttingen, Institut für Astrophysik, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany</footnote>
dynamical simulations are presented and briefly discussed in Sect. 3. In Sects. 4 to 6, we evaluate how large-scale disk structures may be observed with different facilities. Conclusions are drawn in Sect 7.

## 2. Modeling planet-disk interaction in debris disks

In this section, we describe the initial conditions employed for modeling planet-disk interaction. We simulate the dynamical evolution of a large ensemble of dust particles in the gravitational potential of a central star under the influences of planetary perturbations, radiation pressure, Poynting-Robertson effect, and stellar-wind drag. For these simulations, the code MODUST (Rodmann 2006) is employed. Electromagnetic (Lorentz) forces on charged dust particles are neglected (Gustafson 1994; Holmes et al. 2003). Mutual collisions of dust particles are not considered in the dynamical modeling. A brief description of the methods and numerical approach to simulate the dynamical evolution of the dust grains as well as of the approach to simulate images from the resulting density distributions can be found in the Appendix.

### 2.1. Initial conditions

We use several initial conditions (planetary orbit and mass as well as initial dust distribution) selected by analogy to our solar system and to other known debris disks in order to study the structures that can be expected depending on the configuration of the respective system. Since the simulation of the resulting dust distribution is very time-consuming, we are not able to perform a dense sampling of a large, high dimensional parameter space. Therefore, we focus on the following parameter regime:

- one central star (no multiple systems);

- a solar-type star as central star, photospheric emission realized by a black body with $T_{\text{eff}} = 5778$ K, $L_{\star} = 1.0 L_{\odot}$, $M_{\star} = 1.0 M_{\odot}$, $\xi = 0.35$ (Gustafson 1994);

- one planet (considered to dominate the dynamics of the system);

- one initial dust disk (no initial multi-ring systems);

- planetary orbit and initial dust disk are coplanar.

Astronomical silicate with a bulk density of $2.7 \text{ g/cm}^3$ (Draine & Lee 1984; Weingartner & Draine 2001) is employed for the chemical composition of the dust. In the following, the initial conditions for the dust distribution (spatial distribution and grain size distribution) and for the planetary orbit are described and motivated.

### 2.1.1. Treatment of dust creation

We assume the dust to be produced through collisions in a disk of planetesimals. It is then redistributed because of the effects of stellar radiation and wind, and the gravitational interaction with the planet (Eq. (A.1)). This is realized by placing the initial dust distribution at the same position as the planetesimal disk postulated to produce the dust. Particles that are lost because they are sublimated or ejected from the system are replaced by new ones from this reservoir of initial dust grains. Furthermore, the grain size distribution of all grains in the system applied follows a power-law with exponent -3.5 as expected from an equilibrium collisional cascade (Dohnanyi 1969). This approach is used to mimic dust production through collisions of the parent bodies. Note that the local grain size distribution may be significantly

different due to the redistribution of the grains through dynamical interaction with the planet as well as through the effect of radiation pressure and Poynting-Robertson drag (Moro-Martín & Malhotra 2002, 2003).

### 2.1.2. Grain size distribution

The grain size can be expressed by the corresponding value of $\beta$. The correlation between these two quantities is given by

$$ \beta = \frac{3L_{\star}}{16\pi c G M_{\star} \rho a} \approx 575 \left(\frac{\rho}{\text{kg m}^{-3}}\right)^{-1} \left(\frac{a}{\mu\text{m}}\right)^{-1} \left(\frac{L_{\star}}{L_{\odot}}\right) \left(\frac{M_{\star}}{M_{\odot}}\right)^{-1}, \quad (1) $$

where $a$ is the grain size, $Q_{\text{pr}}$ is the efficiency of radiation pressure on a grain, $L_{\star}$ and $M_{\star}$ are the luminosity and mass of the star, $c$ is the speed of light, $G$ is the gravitational constant, and $\rho$ is the bulk density of the dust grains which are assumed to be spherical and compact. The approximation was given by Wyatt et al. (1999) and is valid in the geometrical optics approximation. Since the critical parameter for the computations is $\beta$, it is used in MODUST directly to describe the grain properties (instead of grain size and optical properties). Around a solar-type star, the approximation gives good results for $\beta \le 0.5$ (Rodmann 2006). The lower boundary of the applied grain size distribution is defined by $\beta = 0.5$ ($a \approx 0.43 \mu\text{m}$), because smaller grains (higher values of $\beta$) are ejected from the system by radiation pressure. The upper boundary of the distribution is set to a grain size of $2 \text{ mm} (\beta \approx 1.06 \times 10^{-4}$ for the applied chemical composition). The contribution of larger grains to the simulated observations is neglected.

The grain size distribution is sampled with 50 grain sizes between the lower and the upper grain size, distributed logarithmically (to properly sample the thermal re-emission of different grain sizes for the final images). For each size bin, the dust spatial distribution is simulated. A total number of 1000 test particles is used in each size bin. The system evolves over five times the corresponding Poynting-Robertson time scale of grains of the corresponding $\beta$ at the outer edge of the disk, but not longer than the assumed age of the system (see motivation of the single runs, Sect. 2.1.4). At the end of each run, 100 snapshots of the particle distribution (Rodmann 2006) are taken in the reference frame corotating with the planet to virtually increase the number of particles in each size bin. The snapshots are equally distributed over the last 10% of each run time (see description of the different runs, Sect. 2.1.4).

### 2.1.3. Spatial dust distribution

Gaussian distributions are applied to set the initial eccentricities and inclinations of the test particles. At the start of the integration, the majority of dust grains have orbital eccentricities between 0 and ~0.4 (rms of the Gaussian distribution =0.15). The full width at half maximum of the inclination distribution is set to 10°. These values are consistent with measurements on dust parent bodies in the Kuiper belt of the solar system (Jewitt et al. 1996; Vitense et al. 2010). The semi-major axes of the particles are distributed following a power-law (Wolf & Hillenbrand 2003) from a lowest to a highest value considered as the inner and outer radius of the disk.

### 2.1.4. Explored parameter space

The parameter space explored is motivated by analogy to our solar system and to other known debris disks. We use the position
<table><caption>**Table 1.** Initial conditions for the different runs with MODUST.</caption><thead><tr><th>Run</th><th><i>R</i><sub>in</sub> [AU]</th><th><i>R</i><sub>out</sub> [AU]</th><th>α</th><th><i>M</i><sub>pl</sub> [M<sub>J</sub>]</th><th><i>a</i><sub>pl</sub> [AU]</th><th><i>e</i><sub>pl</sub></th></tr></thead><tbody><tr><td>Ia</td><td>5.0</td><td>5.5</td><td>0.0</td><td>1.0</td><td>5.0</td><td>0.0</td></tr><tr><td>Ib</td><td>5.0</td><td>5.5</td><td>0.0</td><td>1.0</td><td>3.15</td><td>0.0</td></tr><tr><td>Ic</td><td>50.0</td><td>55.0</td><td>0.0</td><td>1.0</td><td>50.0</td><td>0.0</td></tr><tr><td>Id</td><td>50.0</td><td>55.0</td><td>0.0</td><td>1.0</td><td>31.5</td><td>0.0</td></tr><tr><td>IIa</td><td>70.0</td><td>250.0</td><td>-0.5</td><td>0.5</td><td>70.0</td><td>0.0</td></tr><tr><td>IIb</td><td>70.0</td><td>250.0</td><td>-0.5</td><td>1.0</td><td>70.0</td><td>0.0</td></tr><tr><td>IIc</td><td>70.0</td><td>250.0</td><td>-0.5</td><td>5.0</td><td>70.0</td><td>0.0</td></tr><tr><td>IId</td><td>70.0</td><td>250.0</td><td>-0.5</td><td>1.0</td><td>70.0</td><td>0.1</td></tr><tr><td>IIIa</td><td>35.0</td><td>210.0</td><td>-0.5</td><td>0.5</td><td>70.0</td><td>0.0</td></tr><tr><td>IIIb</td><td>35.0</td><td>210.0</td><td>-0.5</td><td>1.0</td><td>70.0</td><td>0.0</td></tr><tr><td>IIIc</td><td>35.0</td><td>210.0</td><td>-0.5</td><td>5.0</td><td>70.0</td><td>0.0</td></tr><tr><td>IIId</td><td>35.0</td><td>210.0</td><td>-0.5</td><td>1.0</td><td>70.0</td><td>0.1</td></tr></tbody></table>

and shape of the dust disk as well as the position, mass, and eccentricity of the planet as parameters for our initial conditions. The different combinations of values explored for these parameters are listed in Table 1 and described and motivated below.

A first sequence of initial conditions (Ia to Id) places the planet in 1:1 or in 2:1 resonance with the initial planetesimal belt. The initial dust disk is a narrow ring ($R_{out} = 1.1 R_{in}$) with constant surface density and two realizations of $R_{in}$ (5 AU for runs Ia and Ib and 50 AU for runs Ic and Id). The planet has a mass of 1 $M_J$ and is on a circular orbit. The two assumed positions of the dust ring are motivated by our solar system (Asteroid belt and Kuiper belt) and typical positions of the dust in other known debris disks around Sun-like stars (e.g., $\epsilon$ Eri, Backman et al. 2009; HD 105). An age (run time of the simulation) of 50 Myr is assumed for these runs (e.g., the approximate age of HD 105; Apai et al. 2008).

The second sequence of initial conditions (IIa to IIId) is motivated by the HD 107146 debris disk (Ardila et al. 2004; Ertel et al. 2011). This disk is particularly well-suited for studying structures in debris disks because of its nearly face-on orientation and the large amount of complementary data available. Furthermore, (sub-)mm observations of the disk have been published by Corder et al. (2009) and Hughes et al. (2011), showing clumpy structures in the surface brightness that could be the signposts of planet-disk interaction<fnref n="1" />. Here, we employ a set of initial conditions where the planet is orbiting at the inner edge of an extended disk with a large inner hole. The inner and outer radius of the disk are set to $R_{in} = 70$ AU and $R_{out} = 250$ AU. The radial surface density distribution is described by a power-law with an index of $\gamma = -0.5$ (e.g., Wolf & Hillenbrand 2003). This model is similar to the results of detailed multi-wavelength modeling of the HD 107146 debris disk (Ertel et al. 2011), although a more complex radial density distribution has been found there. However, such a complex distribution might be the result of the interaction of the disk with a possible planet and the temporal evolution of the disk. We explore different parameters for the planet. It is on a circular orbit and has a mass of 0.5 $M_J$, 1.0 $M_J$, and 5.0 $M_J$ in the runs IIa, IIb, and IIc. It is on an eccentric orbit ($e = 0.1$) and has a mass of 1.0 $M_J$ in run IIId. The age of the system is increased compared to that employed in our sequence I to account for the slightly higher age of HD 107146 (e.g., Roccatagliata et al. 2009). An age of 100 Myr is employed.

A third sequence of models (IIIa to IIId) places the planet within a broad disk. Therefore we modify the initial conditions of our sequence II, so that the disk is now closer to the star, ranging from $R_{in} = 35$ AU to $R_{out} = 210$ AU. The inner radius is comparable to that of the Kuiper belt. The power-law index of the radial dust distribution and the parameters of the planet remain the same.

## 2.2. Image creation

We employ Mie theory to compute the optical properties of the dust using the tool miex (Wolf & Voshchinnikov 2004). To create images from the derived dust distributions, we compute the thermal re-emission and the scattered stellar light for each dust particle at a set of observing wavelengths. The fluxes are then projected onto the plane of the sky considering the scattering phases. This can be performed for arbitrary orientations of the disk (assuming an optically thin debris disk). The image extension in one dimension used is 301 pixels for sequence I (resulting in a pixel resolution of 0.1 AU for runs Ia and Ib and 1.0 AU for runs Ic and Id) and 251 pixels for sequence II and sequence III (resulting in a pixel resolution of 2.0 AU). The images for different grain sizes are then weighted following the grain size distribution and added to compose a final image for each model at each wavelength.

## 2.3. Comparison to previous simulations

In the images produced from our simulations (Figs. 1 and 2), it is possible to identify the structures described by Kuchner & Holman (2003) using simple geometrical arguments. The general goal and approach of our simulations is similar to those presented by Stark & Kuchner (2008). However, our models cover a different region of the parameter space (more massive planets, farther away from the star, larger dust grains with smaller $\beta$). While Stark & Kuchner (2008) concentrated on the effects of terrestrial planets on exozodiacal dust clouds, the scope of the present work is to search for structures induced by giant planets in debris disks. Furthermore, we attempt to simulate observations of the structures found by our modeling. This requires a good sampling of the grain size dependent emissivity of the dust and, consequently, of the grain size distribution. The grains in Stark & Kuchner (2008) are assumed to be no larger than ~120 µm and to be produced by an external source (i.e., a planetesimal belt at considerable distance from the modeled region), which is a reasonable assumption for exozodiacal dust. In the present case, the dust-producing planetesimal belt has to be taken into account. Accordingly, larger grains have to be considered. Furthermore, these grains cannot be neglected, since we attempt to simulate observations at (sub-)mm wavelengths, where millimeter-sized dust grains significantly contribute to the dust emission. Thus, a total of 50 grain size bins is used in the present work in contrast to the five bins in Stark & Kuchner (2008), and we consider grains up to a size of 2 mm.

Nonetheless, the results from both studies can be compared qualitatively. In particular, for models that reproduce similar parameter regimes the results are very similar and differences can be attributed to the above described differences in the approaches and to the differences in the nature of the images themselves. The regime of high planetary mass, far from the star in the parameter space considered by Stark & Kuchner (2008) can be compared to the regime of low planetary mass, close to the star in the parameter space considered in the present work, while images at short

<footnote n="1">¹ However, it is doubtful whether these structures are real, since they are only visible at a signal-to-noise ratio (S/N) of ≈5 to 6 (2 to 3 above the average surface brightness of the disk) and their position in the two observations by Corder et al. (2009) and Hughes et al. (2011) seems to be inconsistent.</footnote>
<figure><img src="image_2.png" /></figure>

**Fig. 1.** Simulated images from model sequence I for a face-on orientation of the disk at different wavelengths. The wavelength is indicated at the top of each column. The star is at the center of each image. The position of the planet is indicated by the red cross in each image. The flux is given in arbitrary units and displayed in a logarithmic stretch from zero to peak value. The region marked by a red circle in some images has been attenuated by a factor *f* to properly display the dynamical range of the images. The value of *f* is then indicated in the lower-left corner of each image. If no value of *f* is given, no attenuation has been applied. Note that a disk without a planet would show no azimuthal structures, but a smooth, featureless radial surface brightness distribution (besides inner and outer disk radius).

wavelengths – tracing the small particles in the present work – lead to comparable structures. This can be seen from comparison of Fig. 6 (lower left panel) or Fig. 10 of Stark & Kuchner (2008) to Fig. 1 (model Ib at thermal re-emission wavelengths, model Id at 24 µm and at 70 µm) and Fig. 2 (models IIa and IIb at wavelengths ≤70 µm) of the present work.

Wyatt (2006) modeled planet-disk interaction in debris disks using dynamical simulations including a migrating planet. Although our approach is substantially different, both works attempt to model structures in debris disks. Therefore, we briefly compare the two approaches. In contrast to Wyatt (2006), we use no migrating planet to initially trap the planetesimals into resonance, but instead assume the dust to be produced through collisions in a ringlike, featureless disk. Furthermore, we use a much simpler approach to mimic dust creation in the planetesimal belt. On the other hand, Wyatt (2006) investigated the libration of particles from the initial resonances only, while we include the possibility that particles are trapped into other resonances than those in which they are initially (they do not have to be in resonance initially at all). Thus, the migrating planet is not necessary to produce resonant structures with our approach.

The differences in the two approaches basically result in an inverse situation in the resulting model images. In Wyatt (2006), very prominent structures are seen at long wavelengths (tracing large particles that remain in the resonances which they are in initially). At short wavelengths, the disks appear smoother, because the particles emitting efficiently at these wavelengths cannot be held in the initial resonance. In contrast, our approach results in a much smoother appearance of the disk at long wavelengths that trace particles that cannot be trapped into resonance

efficiently, because they are not moving radially through the system. Hence, strong structures at (sub-)mm wavelengths are only found in our work if the initial dust distribution is placed close to a strong resonance with the planet, or if the dust is close enough to the star that even millimeter-sized particles are significantly affected by Poynting-Robertson drag. On the other hand, small particles are easily trapped into resonances. This is because they move radially through the system due to Poynting-Robertson drag and, thus, can “find” those resonances.

**3. Results from the dynamical modeling**

From the computed images, one can draw several conclusions on the nature and strength of the structures, on the requirements for the generation of strong structures, and on the wavelength dependence of structures in spatially resolved images. The results can be seen in Figs. 1 to 3 and are briefly discussed in the following.

**3.1. The face-on case**

For a face-on or nearly face-on oriented disk, the azimuthal and radial disk structure can be observed directly. This allows one to draw strong conclusions on the disk structure from observations directly, without the need of detailed modeling or disentangling the real structure from projection effects. On the other hand, it is impossible to draw any strong conclusion on the vertical structure of the disk.

*Prominent structures at short wavelengths.* Observations at short wavelengths (optical and near-infrared scattered light
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Same as Fig. 1, but for model sequences II and III. The flux is given in arbitrary units and displayed in a logarithmic stretch from zero to 1/2 of the peak value at wavelengths up to 24 µm and to twice the peak value at longer wavelengths. Note that a disk without a planet would show no azimuthal structures, but a smooth, featureless radial surface brightness distribution (besides inner and outer disk radius).</figcaption></figure>

as well as short wavelength thermal re-emission) predominantly trace small grains. These grains are subject to efficient Poynting-Robertson drag which causes them to move radially through the system. Thus, they can easily reach regions in the disk, where they can be trapped into resonance by the planet. This causes prominent structures in the distribution of small grains. The shape and strength of these structures depend on the mass of the planet. In general, a more massive planet is able to keep a larger region clear of dust than a less massive one. On the other hand, less massive planets produce a very prominent, bar-like structure.

*Prominent structures at small distance from the star.* In a similar way as before, particles at small radii are more efficiently affected by Poynting-Robertson drag than particles at larger distance from the star. Thus, the same resonant configuration of

planet and disk results in more prominent structures in the distribution of larger grains (at long thermal re-emission wavelengths) if placed closer to the star.

*Prominent 1:1 resonance at long wavelength thermal re-emission.* Large particles are mostly traced by observations at long wavelengths, because they emit more efficiently at these wavelengths and, less important, are cooler than smaller grains. These grains are not moving significantly in radial direction through the system because of Poynting-Robertson drag (as long as they are not too close to the star, see previous point). Only their eccentricity causes a small, periodic change of their radial distance from the star in addition to the dynamical perturbation by the planet. Hence, they can only be trapped into resonance at the radial position they are initially placed. This results in prominent structures caused by 1:1 resonance with the planet, while
<figure><img src="image_4.png" /></figure>

**Fig. 3.** Same as Figs. 1 and 2, but for an edge-on orientation of the disk. The position of the planet is indicated by the red arrow in each image. The flux is given in arbitrary units and displayed in a logarithmic stretch from zero to peak value for model sequence I, from zero to 1/2 of the peak value at wavelengths up to 24 µm for model sequence II and III, and at longer wavelengths from zero to 1/2 of the peak value for model sequence II and to the peak value for model sequence III. Note that a disk without a planet would show a smooth, featureless radial surface brightness distribution (besides the possible signpost of inner and outer disk radius).

other resonances can be neglected in most cases. The orbital ve-
locity of these large grains at the same distance from the star
like the planet is nearly the same as that of the planet due to the
negligible effect of radiation pressure. This renders the resulting
structures more stable and, consequently, even more prominent.
In the case of a narrow, ringlike shape of the disk which is ini-
tially placed close to the planet, this results in a very prominent
horseshoe structure (at long wavelength thermal re-emission). If
the planet is close to the inner edge or within a broad disk, this
results in a horseshoe-like structure in addition to an otherwise
nearly featureless disk. If the planet is placed far away from the
disk, so that it cannot influence the disk by 1:1 resonance, the
effect of the planet is very small. In any case, the prominence of
the structures is increasing with increasing mass of the planet.

*Prominent gaps.* If the planet is placed within a broad disk, it opens a ringlike gap in the disk. This gap is very prominent at long wavelengths, while at short wavelengths strong

resonant structures dominate the appearance of the disk. At long
wavelengths, a horseshoe structure of large particles trapped
into 1:1 resonance is visible (see above). The width of the gap
and the prominence of the horseshoe structure are increasing
with increasing mass of the planet.

*Faint inner disks.* In all our simulations, we always find a more or
less significant amount of predominantly small grains which the
planet is unable to prevent from moving toward the star. These
particles move inward till they reach their sublimation radius.
This inner disk, consisting of small, warm grains, is particularly
bright at short wavelengths (thermal re-emission and scattered
light).

*Accumulation of dust close to the planet.* There is a small amount of dust accumulated close to the planet. These particles are directly captured by the planet. This results in an increased brightness of the planet (and its dust envelope), which increases the chances for direct detection. This scenario has also been
suggested to explain the high brightness of Fomalhaut b (Kalas et al. 2008).

## 3.2. The edge-on case

If a disk is seen edge-on, one is faced with a number of challenges when deriving its radial, azimuthal, and vertical structure. Since debris disks are optically thin, one integrates all flux on the line of sight. This results in strong degeneracies between radial and azimuthal structures. Furthermore, vertical structures cannot be assigned to a particular radial and azimuthal position in the disk, but only to a distance from the star projected onto the sky plane. However, most of the known, spatially resolved debris disks are seen close to edge-on. This is most likely because of the higher surface brightness due to more emitting material on the line of sight which, results in an observational bias. For that reason, it is particularly interesting to search for prominent structures in edge-on seen debris disks due to planet-disk interaction, which can be unambiguously identified.

**Particles scattered out of the disk midplane.** Particularly for massive planets close to the radial position of the dust production, there is a halo of particles scattered out of the disk midplane by the interaction with the planet. Therefore, the halo consists mostly of small and intermediate-sized grains moving efficiently radially through the disk. At short to intermediate thermal re-emission wavelengths, resonant structures are visible. These structures are particularly prominent at these wavelengths, because they mostly consist of intermediate-sized grains. These grains are sufficiently affected by Poynting-Robertson drag and radiation pressure to interact with the planet and, thus, to be scattered out of the disk midplane. On the other hand, Poynting-Robertson drag is not strong enough to prevent these grains from being trapped into resonance.

**Asymmetries in the disk radial brightness profile.** There is a brightness asymmetry between the two ansae of the disk. The disk ansa containing the planet is slightly brighter than the opposite ansa. This is due to the accumulation of dust in the resonances which results in clumps of higher density (see face-on case).

**Multiple peaks and dips in the disk radial brightness profile.** A multi-ring structure as seen in particular in the long wavelength results from sequence III results in a wavy radial brightness profile along the disk midplane in the edge-on case, while peaks are seen at the position of rings, and dips are seen at the position of the gaps in the disk. The strength of these wavy structures depends on the exact configuration of the system and on the observing wavelength.

## 3.3. Age dependence of the phenomena

The dominating effect in our model that causes the dust grains to change their orbits beside the gravitational interaction with the planet is Poynting-Robertson drag. Therefore, one has to compare the Poynting-Robertson time scale of different grains with the age of the systems. This time scale (the time it takes for a particle to spiral onto the star starting at a circular orbit at distance $R_0$) can be approximated following Gustafson (1994):

$$t_{PR} \approx \frac{400}{\beta} \left(\frac{M_{\star}}{M_{\odot}}\right)^{-1} \left(\frac{R_0}{AU}\right)^2 \text{yr.} \qquad (2)$$

It is short for small grains (large $\beta$) and increases with grain size (decreasing $\beta$). It also depends on the radial distance from

the star. Accordingly, systems with the dust placed closer to the star evolve faster. Furthermore, structures in small grains evolve faster than large grains. The older a system is, the more time there is for larger grains to move significantly in the radial direction, and thereby to be trapped into different resonances. With increasing age of the system, larger grains will have the time to be trapped in prominent resonant structures instead of in the single 1:1 resonant horseshoe structure (in the face-on case). Because these larger grains are traced at longer wavelengths, the transition between the two types of structures will occur at longer wavelengths for older systems.

It is important to note that no collisions between the dust particles are considered in our dynamical simulations. Structures that appear only after a long time or that are particularly traced by particles that are very abundant may be destroyed by chaotic events like collisions. Furthermore, particles may be destroyed through collisions before they are able to form these structures. The strength of this effect is decreasing with decreasing dust mass. Thus, such structures are a result of our modeling, but their significance might be over estimated as compared to the strength of these structures in real systems (particularly in systems with very massive disks).

## 3.4. Potential for observations and modeling

From the above discussions, we found that planet-disk interaction in debris disks produces structures in the disk that allow one to constrain the parameters of a planet-disk system. However, data at one of the wavelength regimes considered are expected to provide only weak constraints on the actual configuration of the system (e.g., planetary mass, major axis, eccentricity, radial distribution of the parent bodies), especially when taking into account observational effects such as noise, resolution effects, and stellar point spread function (PSF) subtraction uncertainties (see Sects. 4 to 6 for simulated observations and discussion). In the model images themselves, this is particularly obvious, when comparing the images at mid-infrared wavelengths (and shorter) in face-on orientation produced from our runs II and III (e.g., at 24 µm) and the corresponding (sub-)mm images (Fig. 2). While the mid-infrared data allow one to distinguish particularly well between planetary mass and eccentricity (e.g., between run IIa and IIc), the (sub-)mm data allow one to constrain in particular the major axis of the planet and the position of the parent bodies (e.g., compare images from runs IIa and IIIa). Thus, combined observations at both wavelength regimes will be most useful for a detailed analysis of planet-disk systems.

# 4. Evaluating the observability (1): ALMA

In the following sections, we simulate observations of the modeled disk images using different near-future facilities. In addition, we discuss the detectability of the structures with the Hubble Space Telescope (HST) in order to explain why the prominent scattered light structures seen in our simulations have not been detected yet. In this section, we start with ALMA.

## 4.1. General sensitivity study for debris disks

ALMA will allow interferometric observations of debris disks in the (sub-)mm regime with unprecedented sensitivity and spatial resolution. Consequently, it is the most promising instrument to study structures in debris disks in this wavelength regime in the near future. ALMA is a very complex instrument, providing the
<table><caption>**Table 2.** Longest baseline $B_{\text{max}}$, shortest baseline $B_{\text{min}}$, and resolution characteristics of the array configurations used for the simulations.</caption><thead><tr><th rowspan="2">Array number</th><th rowspan="2">B<sub>min</sub> [m]</th><th rowspan="2">B<sub>max</sub> [m]</th><th colspan="2">Band 10</th><th colspan="2">Band 7</th><th colspan="2">Band 6</th><th colspan="2">Band 5</th></tr><tr><th>δ<sub>min</sub> ["]</th><th>δ<sub>max</sub> ["]</th><th>δ<sub>min</sub> ["]</th><th>δ<sub>max</sub> ["]</th><th>δ<sub>min</sub> ["]</th><th>δ<sub>max</sub> ["]</th><th>δ<sub>min</sub> ["]</th><th>δ<sub>max</sub> ["]</th></tr></thead><tbody><tr><td>01</td><td>15</td><td>161</td><td>0.45</td><td>4.8</td><td>1.22</td><td>13.1</td><td>1.60</td><td>17.2</td><td>2.05</td><td>22.0</td></tr><tr><td>03</td><td>15</td><td>260</td><td>0.28</td><td>4.8</td><td>0.75</td><td>13.1</td><td>0.99</td><td>17.2</td><td>1.27</td><td>22.0</td></tr><tr><td>05</td><td>15</td><td>390</td><td>0.19</td><td>4.8</td><td>0.50</td><td>13.1</td><td>0.66</td><td>17.2</td><td>0.85</td><td>22.0</td></tr><tr><td>07</td><td>15</td><td>538</td><td>0.13</td><td>4.8</td><td>0.36</td><td>13.1</td><td>0.48</td><td>17.2</td><td>0.61</td><td>22.0</td></tr><tr><td>09</td><td>15</td><td>703</td><td>0.10</td><td>4.8</td><td>0.28</td><td>13.1</td><td>0.37</td><td>17.2</td><td>0.47</td><td>22.0</td></tr><tr><td>11</td><td>15</td><td>1038</td><td>0.07</td><td>4.8</td><td>0.19</td><td>13.1</td><td>0.25</td><td>17.2</td><td>0.32</td><td>22.0</td></tr><tr><td>13</td><td>24</td><td>1440</td><td>0.05</td><td>3.0</td><td>0.14</td><td>8.2</td><td>0.18</td><td>10.7</td><td>0.23</td><td>13.8</td></tr><tr><td>15</td><td>24</td><td>1811</td><td>0.040</td><td>3.0</td><td>0.108</td><td>8.2</td><td>0.140</td><td>10.7</td><td>0.180</td><td>13.8</td></tr><tr><td>17</td><td>49</td><td>2297</td><td>0.031</td><td>1.5</td><td>0.085</td><td>4.0</td><td>0.112</td><td>5.3</td><td>0.144</td><td>6.7</td></tr><tr><td>19</td><td>49</td><td>3105</td><td>0.023</td><td>1.5</td><td>0.063</td><td>4.0</td><td>0.083</td><td>5.3</td><td>0.106</td><td>6.7</td></tr><tr><td>21</td><td>79</td><td>6068</td><td>0.012</td><td>0.9</td><td>0.032</td><td>2.5</td><td>0.042</td><td>3.3</td><td>0.054</td><td>4.2</td></tr><tr><td>23</td><td>79</td><td>11 457</td><td>0.006</td><td>0.9</td><td>0.017</td><td>2.5</td><td>0.022</td><td>3.3</td><td>0.029</td><td>4.2</td></tr><tr><td>25, 27</td><td>79</td><td>14 444</td><td>0.005</td><td>0.9</td><td>0.014</td><td>2.5</td><td>0.018</td><td>3.3</td><td>0.023</td><td>4.2</td></tr></tbody></table>

user with multiple options that allow one to optimize the instrument set-up for the planned observations. A detailed study of the performance of each configuration is necessary to optimize the observational outcome and efficiency. This study in the context of spatially resolved observations of debris disk is described in the following.

### 4.1.1. Model description

All simulations are carried out for a debris disk around a solar-type star (realized as a black body radiator with $R_{\star} = 1.0 R_{\odot}$, $L_{\star} = 1.0 L_{\odot}$, $T_{\text{eff}} = 5778$ K). A simple, analytical disk model is used to keep the results as general as possible. It consists of a circular ring ($r_{\text{out}} = 1.1 r_{\text{in}}$) of dust with constant surface density (similar to the width of a ring of dust particles on orbits with identical major axis and eccentricity of 0.1, if the directions of the major axes of the orbits are distributed randomly within the plane of the disk) and an opening angle of $10^{\circ}$. Astronomical silicate with a bulk density of $2.7 \text{ g/cm}^3$ (Draine & Lee 1984; Weingartner & Draine 2001) is employed for the chemical composition of the dust. The grain size distribution follows a power-law distribution (exponent -3.5) with lower cut-off size $a_{\text{min}} = 0.45 \mu\text{m}$ (i.e., the blow-out size of the system) and upper cut-off size $a_{\text{max}} = 2 \text{ mm}$. Only observing wavelengths up to 1.6 mm are considered. Therefore, the emission of larger grains is negligible for a sufficiently steep grain size distribution. The inner disk radius is set to 5 AU, 50 AU, and 100 AU in different runs. This covers a variety of known debris disks, e.g., our solar system, $\epsilon$ Eri (Greaves et al. 1998; Backman et al. 2009), HD 107146 (Ardila et al. 2004; Corder et al. 2009; Ertel et al. 2011), Fomalhaut (Kalas et al. 2005), and AU Mic (Augereau & Beust 2006). The dust mass of the disk is set to $10^{-8} M_{\odot}$ for simulating the model images and scaled to different masses later (assuming an optically thin debris disk) to evaluate the sensitivity. These masses can be converted into fractional luminosities following the equation

$$ \frac{L_d}{L_{\star}} = 4.5 \times 10^7 \left(\frac{r}{\text{AU}}\right)^{-2} \frac{M_{\text{dust}}}{M_{\odot}}, \quad (3) $$

i.e., the dust luminosity for a given dust mass and stellar luminosity scales with $r^{-2}$ and is proportional to the dust mass, where $4.5 \times 10^7$ is the factor for $r = 1$ AU given the above model of the system (stellar properties and disk properties). This is because the total luminosity of a given dust species only depends on the stellar spectrum, the total surface of the dust ($\propto$ dust mass), and

its absorption efficiency, assuming thermal equilibrium (no information about the wavelength range is carried in which the majority of this luminosity is emitted, i.e., about the dust temperature). The resulting factors are $1.8 \times 10^6$, $1.8 \times 10^4$, and $4.5 \times 10^3$ for the disk models with $r_{\text{in}} = 5$ AU, 50 AU, and 100 AU in the present study.

Simulations of images from these disk models are performed for face-on and edge-on orientation of the disk. The pixel resolution of the images in AU can then be computed as $2 r_{\text{out}}/101$ (the total extent of the disk is chosen to be 101 pixels in one direction) and in arcseconds by dividing the result by the distance of the system in pc. Different distances are assumed to simulate observations (Sect. 4.1.2). If the total extent of the images as described above is smaller than five times the resolution element (FWHM) of the observations (in particular for the models with $r_{\text{in}} = 5$ AU at large distance and small array extent), the background fluctuation cannot be quantified in a reliable way. Thus, the resolution of the images is lowered in these cases and the empty region around the disk image is increased. This is done dynamically, so that the pixel resolution is 1/20 of the FWHM and the total extent of the images is 200 pixel ($10 \times \text{FWHM}$ of each observation).

### 4.1.2. Simulation of observations

The CASA ALMA simulator<fnref n="2" /> (procedure `simdata`) is used to simulate observations of the model images. The simulations are conducted using the configurations of the ALMA array in full operations provided by the CASA simulator (using only every second configuration from configurations 01 to 27)<fnref n="3" />. These array configurations contain 50 12-m antennas each. The baseline ranges of the array configurations are listed in Table 2. Bands 10, 7, 6, and 5 (central wavelengths of 350 µm, 950 µm, 1250 µm, and 1600 µm)<fnref n="4" /> are used. The exact central wavelength of each bandpass is chosen to avoid strong atmospheric absorption bands. The band width used is 7.5 GHz for all simulations following the ALMA Cycle 0 Technical handbook<fnref n="4" />. The objects are placed at an optimal position in the sky reaching the zenith during observations (RA = 18h, Dec = -23°). Total

<footnote n="2">² http://casa.nrao.edu/</footnote>

<footnote n="3">³ For a detailed explanation of all parameters used in this procedure, see the corresponding manuals at</footnote>

<footnote>http://casaguides.nrao.edu/index.php?</footnote>

<footnote>title=Simulating_Observations_in_CASA</footnote>

<footnote n="4">⁴ http://almascience.eso.org/document-and-tools</footnote>
<figure><img src="image_5.png" /></figure>

**Fig. 4.** Contribution of phase noise to the total noise in our ALMA observations (Ruge, priv. comm.). The noise levels are computed as the standard deviation of all pixels in simulated observations of empty images (all pixels set to 0.0). The quantity $\sigma_{th}$ represents the standard deviation in simulations with only thermal noise considered, while $\sigma_{th+ph}$ is derived from simulations considering thermal and phase noise. Ratios higher than 1.0 may occur because of uncertainties in noise estimates and because of different noise contributions that in part cancel out each other. The curves for Band 5 and Band 6 are very close to 1.0 and cannot be distinguished in the plot. In any case, the contribution of phase noise is lower than 5%.

observing time is 8 h, while the single integration time is 60 s. Thermal noise is added to the simulated visibilities. Good, but realistic weather conditions are assumed (precipitable water vapor = 0.65 mm, ground temperature = 269 K<fnref n="5" />). From a comparison with thermal noise, we find phase noise to be negligible (Fig. 4) and, therefore, neglect this contribution to the total noise. For image reconstruction, natural weighting of the visibilities is applied. Deconvolution is performed with 500 iterations and a threshold of 0.01 mJy/beam. Single pointing observations are simulated rather than mosaicking, because most of the simulated images fit into one field of view of a single pointing observation. Mosaicking would require significantly more observing time, increasing the total time on target to an unrealistic amount considering the high pressure on ALMA expected.

For the simulated observations of the model images, the disks are placed at 11 different distances distributed logarithmically from 10 pc to 100 pc. The dust mass (total disk flux assuming an optically thin disk) is scaled to different values. The S/N is estimated as the ratio between peak flux in the simulated observations of a target and the rms of the background fluctuation. A linear fit is performed on the distribution of S/N over disk mass considering only data with S/N > 20. From this fit, the disk mass needed to reach an S/N of 10 is computed. This is considered to be a robust detection of the disk.

### 4.1.3. Results

To properly interpret the results, one has to consider the following specifics of an interferometric observation:

- The field of view of an interferometric observation using one single pointing is limited by the FWHM of an observation with one antenna (the primary beam). For the 12-m antennas of ALMA, this is $\lambda/12$ m.

- The spatial resolution $\delta_{min}$ of an observation (the FWHM of the synthesized beam) is determined by the maximum baseline $B_{max}$ following $\delta_{min} = \lambda/B_{max}$.

- The largest angular scale $\delta_{max}$ of a target that will not be filtered out (resolved out) by the interferometric observations is determined by the minimum baseline $B_{min}$ following $\delta_{max} \approx 0.6\lambda/B_{min}^6$.

The derived values of $B_{min}$ and $B_{max}$ as well as the resolution characteristics $\delta_{min}$ and $\delta_{max}$ for the different array configurations and bands considered are listed in Table 2. With the derived values one can now explain the features seen in the sensitivity maps that are shown in Fig. 5. First, there are a number of configurations in the explored parameter space that cannot be observed with the strategy assumed or where the approach used to estimate the sensitivity produces erroneous results:

- In particular very extended, nearby disks, cannot be observed at short wavelengths without mosaicking, since the spatial extent of the disks is larger than the field of view of the observations. The configurations where this is the case are marked by a black cross in an empty (white) pixel in Fig. 5. Note that parts of the disk can still be observed with one single pointing (Boley et al. 2012).

- For the highest resolution (most extended arrays), in particular nearby, extended disks are resolved out. These regions are marked in Fig. 5 by an "x".

- For high spatial resolution, the FWHM of the synthesized beam is similar to, or smaller than, the pixel size of the synthetic images. This results in heavy residuals from image reconstruction and does not allow one to properly quantify the S/N in the simulated maps. Note that this is a limitation of the approach, not a real limitation of ALMA. However, at such high resolution, the disks would have to be extremely massive (bright) to be detected at all (see the discussion of the sensitivity below).

In the regions of the parameter space at which the simulations are not affected by the above limitations, one can evaluate the sensitivity. The following behavior in the parameter space can be found:

- In general, the highest sensitivity to disk mass is reached in Band 6 (central wavelength 1250 µm).

- For a given spatial resolution, the sensitivity is only slightly decreasing with increasing distance as long as the disk is spatially resolved ($FWHM \le 2 r_{out}$). This can be explained by the fact that the surface brightness of the disk (e.g., in mJy/AU²) is decreasing with the distance squared, but the beam covers an increasing area of the disk, depending on the scales of the dominating disk structures.

- For a given distance, the S/N is increasing with decreasing spatial resolution as long as the disk is spatially resolved. This is because the beam covers a larger area of the disk at very similar sensitivity (e.g., in mJy/beam). As long as the FWHM is smaller than the dominating structure of the disk (e.g., the ring width), the sensitivity is increasing with the square of the FWHM.

- As soon as the disk is spatially unresolved, the sensitivity drops with the distance squared for a given spatial resolution.

- In general, disks seen edge-on have a higher surface brightness. Thus, spatially resolved observations of an edge-on

<footnote n="5">⁵ http://www.eso.org/gen-fac/pubs/astclim/</footnote>

<footnote n="6">⁶ ALMA Cycle 0 Proposers Guide, http://almascience.eso.org/document-and-tools</footnote>
<figure><img src="image_6.png" /></figure>

**Fig. 5.** Maps of the 10σ sensitivity of ALMA to face-on seen (**left**) and edge-on seen (**right**) debris disks as a function of the distance and array configuration (spatial resolution) for different wavelengths of observations. The three rows refer to different disk extents given by $r_{\text{in}}$ of 5 AU (**top**), 50 AU (**middle**), and 100 AU (**bottom**). Disk mass can be converted into fractional luminosity following Eq. (3). For details about the model images and about the simulation of observations see Sects. 4.1.1 and 4.1.2. The detectable disk mass is color-coded and can be converted into fractional luminosity following Eq. (3). White pixels with symbols mark positions where the simulations failed or observations are impossible for various reasons denoted by the symbols as follows (see Sect. 4.1.3 for more details): cross – the spatial extent of the disk is larger than the field of view of a single pointing, × – the disk is “resolved out” (the radial width of the ring is larger than the largest scale detectable), open square – artifacts of the image reconstruction caused by deficiencies of the model image dominate the simulated observations.

seen disk will result in a higher S/N than the same obser-
vations of the same disk seen face-on. In addition to that,
there is no significant difference between imaging of edge-
on and face-on disks (for the employed disk models). This
can be explained by the fact that for spatially resolved imag-
ing only the peak S/N is considered, which depends on the
scale and brightness of the brightest structures in the disk.
This is in both cases the radial width of the ring and, for disk
seen edge-on, the vertical height (which is of the same or-
der in our model). For spatially unresolved observations, the
orientation of the disk has no effect at all.

The above qualitative discussion leads to an important recom-
mendation for optimal spatially resolved imaging of debris disks
with ALMA. As long as the structures of interest are bright
enough, one should always use a spatial resolution similar to the
scales of interest. Most importantly, one should avoid to over-
resolve the structures, because sensitivity then decreases with
the square of the spatial resolution.

ALMA will not only increase our knowledge about the
known spatially resolved debris disks owing to the higher sensi-
tivity and spatial resolution. Because of the nearly constant sen-
sitivity to dust mass in spatially resolved imaging at distances
up to ≈100 pc (for disks with radial extents of few tens of AU),
ALMA will significantly increase the sample of debris disks spa-
tially resolved at millimeter wavelengths toward more distant
systems compared to earlier instruments.

## 4.2. Observability of planet-disk interaction

From our dynamical modeling, we select several runs that re-
sult in prominent structures at wavelengths ≥350 µm in face-on
and/or edge-on orientation. The total flux in the selected im-
ages is scaled to the total flux of selected, well-studied debris
disks (Table 3). The CASA simulator is used to predict results

of observations. The optimum array configuration is selected by
the resulting S/N and spatial resolution in the simulated images.

4.2.1. Selection of runs and preparation of the images

The runs for which ALMA simulations are carried out are se-
lected according to the prominence of the structures seen in the
(sub-)mm:

- runs Ib, Ic, IIc, and IIId are found to represent the prominent structures in the synthetic images well for face-on orientation of the disk (Figs. 1 and 2);

- run IIIa is selected as an example for structures seen in edge-on orientation of the disk (Fig. 3).

From these runs, observations are simulated at Band 6 (central
wavelength 1250 µm, found in Chapt. 4.1 to be the most sensitive
one for the intended observations) as follows:

- known debris disks with radial dust distributions similar to the simulated ones are selected as reference objects. Systems around solar-type stars are considered, because solar-type stars are considered for the dynamical simulations as well;

- a representative photometric measurement (in the (sub-)mm, where possible) is used to scale the total flux of the simulated images. Therefore, the total flux of the disk at the wavelength of these observations is computed;

- the disk is placed at the distance of the reference debris disk. A flux-scaling factor is derived by comparing the flux of the modeled disk with the measured flux at the reference wavelength. This scaling factor is applied to the image at Band 6 to properly scale the flux in this image to a realistic value;

- the stellar contribution of the host star of the reference debris disk is added to the center of the synthetic image.

The reference debris disk for each run, the reference wavelength,
the applied flux at this wavelength, some basic information about
<table><caption>**Table 3.** Reference systems for our simulated ALMA observations.</caption><thead><tr><th>Run</th><th>Reference</th><th>λ<sub>ref</sub><br>[μm]</th><th>d<br>[pc]</th><th>Host star</th><th>Age<br>[Myr]</th><th>R<sub>dust</sub><br>[AU]</th><th>F<sub>disk</sub> at λ<sub>ref</sub><br>[mJy]</th><th>M<sub>disk</sub><br>[M<sub>☉</sub>]</th><th>Ref.</th></tr></thead><tbody><tr><td>Ib</td><td>ε Eri</td><td>24</td><td>3.2</td><td>K2 V</td><td>850</td><td>3.0</td><td>330.0</td><td>5.7 × 10<sup>-11</sup></td><td>1</td></tr><tr><td>Ic</td><td>HD 105</td><td>850</td><td>40.0</td><td>G0 V</td><td>30</td><td>45–120</td><td>10.7</td><td>1.2 × 10<sup>-7</sup></td><td>2</td></tr><tr><td>IIc</td><td>HD 107146</td><td>350</td><td>29.5</td><td>G2 V</td><td>130</td><td>50–250</td><td>319.0</td><td>7.3 × 10<sup>-7</sup></td><td>3, 4</td></tr><tr><td>IIIa</td><td>HD 107146</td><td>350</td><td>29.5</td><td>G2 V</td><td>130</td><td>50–250</td><td>319.0</td><td>3.7 × 10<sup>-7</sup></td><td>3, 4</td></tr><tr><td>IIId</td><td>HD 107146</td><td>350</td><td>29.5</td><td>G2 V</td><td>130</td><td>50–250</td><td>319.0</td><td>5.5 × 10<sup>-7</sup></td><td>3, 4</td></tr></tbody></table>

**References.** (1) Backman et al. (2009); (2) Nilsson et al. (2010); (3) Corder et al. (2009); (4) Ertel et al. (2011).

the reference system, and the resulting dust mass when scaling the model disk to the observed flux are shown in Table 3. It is important to note that the model scaled to the observed flux of the reference debris disk at one wavelength is not expected to reproduce the whole spectral energy distribution of the disk properly. However, the reference debris disks are selected by their similarity to the model systems in the radial dust distribution and Band 6 is clearly in the Rayleigh-Jeans regime of the dust re-emission. Accordingly, the flux at this wavelength is considered to be at least realistic for a debris disk similar to the reference disk.

The reference debris disks used for the scaling of the model are motivated as follows:

- ε Eri for run Ib: the inner debris ring (~3 AU from the star; Backman et al. 2009) is used. For this ring, the 24 µm flux is used as reference value to scale the flux in all model images, since the other disk components are expected to have significant contribution to the flux at longer wavelengths.

- HD 105 for run Ic: HD 105 is a massive debris disk around a solar-type star. Most of the dust is expected to be concentrated at a distance of 40 AU to 50 AU from the star (Hillenbrand et al. 2008).

- HD 107146 for runs IIc, IIIa, and IIId: the work on the HD 107146 disk (Ertel et al. 2011) was the motivation for these runs (Sect. 2.1).

### 4.2.2. Parameters for the CASA simulations

The simulations are carried out using all configurations of the ALMA array in full operations provided by the CASA simulator. Observing conditions and parameters for the simulations are identical to the simulations in Sect. 4.1. Band 6 (central wavelength 1250 µm) is used. Single pointing observations are simulated rather than mosaicking, since all simulated images fit into one field of view of a single pointing at Band 6.

### 4.2.3. Results

We find from the simulations of observations at different spatial resolution (different maximum baseline) that a trade-off between sensitivity and spatial resolution is necessary (see also Sect. 4.1.3). This is particularly well illustrated in the example of run IIId (Fig. 6, top). The bar-like structure within the (outer) gap in the disk is only visible at the lowest resolutions, while one needs a higher resolution to clearly resolve the inner gap and to clearly separate the innermost peak from the inner ring of dust at 1" to 2" from the star. For the five model disks considered for the ALMA simulations, the best array configurations based on high S/N and high resolution are selected. The model image resulting from run Ib is found to be too faint to be detected by any of the

simulated observations. For the other four model images, the results are shown in Fig. 6 (bottom). The selection of the optimum array configuration for each case is described in the following. The results are consistent with the suggestions for optimal spatially resolved imaging of debris disks derived in Sect. 4.1.3.

For run Ic, the structures are very smooth and no substructure is seen. Thus, one can use the peak S/N as a good tracer for the significance of the structures detected. Based on an additional, visible inspection of the images, a peak S/N of 22 (array 19, FWHM = 0.11 × 0.15) is found to give the best results, while the next larger array would result in a peak S/N of 15, but large parts of the disk would be detected at an S/N < 10.

For run IIc, the structure of interest would be the faint horseshoe structure at the inner edge of the disk. This structure can in any case only be observed at an S/N of up to 4. This is reached when the size of the beam is comparable to the size of this structure. The optimum is reached using array 08 (FWHM = 0.62 × 0.78). It is important to note that there is a negative background that is nearly homogeneous over the whole image. This can be attributed to the limitations of ALMA to observe very extended structures and of the simulations as described in Sect. 4.1.3. The signal has to be evaluated as the flux above this homogeneous background.

For run IIId, value is placed on the multi-ring structure, while the bar-like structure is ignored. An example for an array configuration that results in a significant detection of the bar-like structure can be seen in Fig. 6 (top). A peak S/N of 10 (array 11, FWHM = 0.38 × 0.47) is found to result in significant detection of the three disk components (outer and inner ring as well as innermost accumulation of dust), since they have a very similar surface brightness. The bar-like structure is visible at low S/N in these simulated observations, as expected from the previous discussion about the connection between spatial resolution and sensitivity to surface brightness.

For the disk from run IIIa seen edge-on, the structure of interest is the gap in the disk seen as a dip in the radial brightness distribution and the outer disk seen as a secondary peak of the surface brightness beyond this gap. An S/N of 15 (array 16, FWHM = 0.16 × 0.22 – similar to the vertical extent of the disk) in the secondary peak of the surface brightness results in a clear detection of both features and allows one to marginally resolve the vertical extent of the disk.

These examples illustrate the findings in Sect. 4.1.3 for structured disks. ALMA is not able to reach both high sensitivity to surface brightness and high spatial resolution simultaneously. The highest sensitivity at reasonable spatial resolution is reached if the resolution of the observations is similar to the scale of the structures expected, as long as these structures are bright enough.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** *Top:* simulated ALMA observations of the model image resulting from run IIId (face-on orientation) for different array configurations. Only the results from a representative number of array configurations are shown. The number in the upper-right corner of each image denotes the array number. The corresponding beam is displayed in the lower-right corner of each image. *Bottom:* simulated ALMA observations of all selected model images. The optimum array configurations have been selected. The run and the array number are displayed in the upper-right corner, the corresponding beam is displayed in the lower-right corner of each image. All observations were simulated for Band 6 (central wavelength 1250 µm).</figcaption></figure>

**5. Evaluating the observability (2): the Hubble Space Telescope**

From the dynamical simulations, prominent structures in debris disks are found, in particular in scattered light and short-wavelength thermal re-emission. In contrast, no debris disk seen face-on is known to exhibit such structures in scattered light observations. Therefore, we briefly examine the observability of such structures in face-on seen disks with present instruments in scattered light, in particular with the HST with which most of the scattered-light detected debris disks have been discovered. Our example here is the HD 107146 debris disk. It is one of the few debris disks seen (close to) face-on detected in scattered light (Ardila et al. 2004; Ertel et al. 2011).

The HST/ASC image in the F606W ($\lambda_c$ = 0.6 µm) filter presented by Ardila et al. (2004) is used to evaluate the observability of simulated structures in these data. Therefore, the model-subtracted image of Ertel et al. (2011) is employed, which represents the pure PSF subtraction residuals (within the capabilities of the modeling). To these data, simulated images from the dynamical modeling (inclination of 25° from face-on applied from Ardila et al. 2004) are added after scaling them to the scattered-light flux of the fitted disk at the same wavelength. This results in simulated observations of the model images using the same strategy as for the original observations. The results for the model images from runs IIb and IIc are shown in Fig. 7.

While the disk is clearly detected, the structures of the modeled images are corrupted by the strong PSF subtraction residuals. Thus, one would not be able to draw strong conclusions on the azimuthal structure of the disk from these data. Furthermore, there is only a small number of debris disks seen face-on and detected in scattered light. In contrast, the two most prominent debris disks seen edge-on, β Pic and AU Mic, exhibit clumpy and warped structures that might be interpreted as structures

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Simulated HST/ACS coronagraphic observations of the simulated images from runs IIb and IIc at 0.6 µm. The disk is inclined by 25° along the *x*-axis, so that the lower half of the disk points toward the observer.</figcaption></figure>

in the disk induced by planet-disk interaction (e.g., Liu 2004; Golimowski et al. 2006). All this suggests that the structures found by the modeling in this work may be present, but not detected in face-on seen debris disks, yet.

**6. Evaluating the observability (3): space-based near- to mid-infrared telescopes**

Near-future space-based near- to mid-infrared telescopes (e.g., the James Webb Space Telescope, JWST, Gardner et al. 2006, or the Space Infrared Telescope for Cosmology and Astrophysics, SPICA, Swinyard et al. 2009) are expected to provide high-sensitivity, high angular resolution imaging capabilities in the
<table><caption>**Table 4.** Reference systems for our simulated JWST observations.</caption><thead><tr><th>Run</th><th>Reference</th><th><i>F</i><sub>disk</sub><br>[mJy]</th><th><i>F</i><sub>*</sub><br>[mJy]</th><th><i>M</i><sub>disk</sub><br>[M<sub>☉</sub>]</th><th>Ref.</th></tr></thead><tbody><tr><td>Ia</td><td>ε Eri</td><td>330.0</td><td>1726</td><td>8.5 × 10<sup>-12</sup></td><td>1</td></tr><tr><td>Ib</td><td>ε Eri</td><td>330.0</td><td>1726</td><td>7.3 × 10<sup>-12</sup></td><td>1</td></tr><tr><td>IIa</td><td>HD 107146</td><td>19.0</td><td>40.8</td><td>4.4 × 10<sup>-8</sup></td><td>2</td></tr><tr><td>IIc</td><td>HD 107146</td><td>19.0</td><td>40.8</td><td>1.4 × 10<sup>-7</sup></td><td>2</td></tr><tr><td>IId</td><td>HD 107146</td><td>19.0</td><td>40.8</td><td>5.7 × 10<sup>-8</sup></td><td>2</td></tr><tr><td>IIIc</td><td>HD 107146</td><td>19.0</td><td>40.8</td><td>7.4 × 10<sup>-9</sup></td><td>2</td></tr></tbody></table>

**Notes.** The reference wavelength at which the fluxes are given is 24 µm for both systems. Also note the information given in Table 3 about the two reference systems.

**References.** (1) Backman et al. (2009); (2) Hillenbrand et al. (2008).

optical to near-infrared. In this section, observations with such facilities are predicted based on available information about the expected capabilities of the JWST. Since the focus of this work is on spatially resolved imaging capabilities, the two relevant instruments are the Near InfraRed Camera (NIRCam) and the Mid InfraRed Instrument (MIRI). NIRCam provides imaging capabilities in the 0.6–5.0 µm wavelength range. In this range, coronagraphy is necessary to block the direct stellar radiation. The HST is very successful in coronagraphic imaging of bright debris disks and the JWST will exceed these capabilities due to its higher sensitivity and spatial resolution. However, the performance of coronagraphic instruments depends on several influences such as PSF stability. Consequently, the exact results of such observations are hard to predict without detailed knowledge about the performance of the instrument during science operations. Hence, the focus of the present section is placed on predictions on the capability of debris disk observations through none-coronagraphic imaging as are possible with MIRI.

MIRI will provide imaging capabilities in the wavelength range of 5.0–29 µm. In particular at the long wavelength edge of this range, thermal radiation of debris disks is expected to have a detectable level. With a spatial resolution of ≈1" and high sensitivity, MIRI is the first instrument that will allow for resolved imaging of a large number of debris disks in this wavelength range.

## 6.1. Selection of runs and approach to simulate observations

Several runs are selected from the dynamical simulations to explore the capabilities of MIRI to detect the structures. A summary of the considered runs and the references used for each run is given in Table 4. Since $\epsilon$ Eri is the only known debris disk this close and it is seen close to face-on, only the face-on orientation is used for the simulated observations on the example of the inner disk in this system. For the simulated observations on the example of the HD 107146 debris disk, face-on and edge-on orientations are used, since a number of disks with different orientations are known at distances of few tens of AU (e.g., AU Mic, Kalas et al. 2004; $q^1$ Eri, Liseau et al. 2010; HD 207129, Krist et al. 2010).

Because the star contributes significantly to the flux at mid-infrared wavelengths, a stellar PSF subtraction is expected to be necessary to clearly reveal any resolved disk structure. To evaluate the observability of structures in debris disks, a two-step approach is used:

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Simulated PSF of JWST/MIRI at 25.5 µm. The image is shown in a logarithmic stretch from 0 to 5% of the peak flux to highlight the high order structure. The white contours represent 10% and 50% of the peak flux.</figcaption></figure>

(1) observations of the disk only are simulated by convolving the synthetic disk images with the telescope PSF. From the resulting images, the sensitivity needed is estimated and the capability to spatially resolve the structures is evaluated;

(2) the stellar PSF contribution is estimated and the PSF subtraction accuracy needed is evaluated.

For simulated observations of the disk only, the images at 24 µm resulting from the dynamical simulations are scaled to the 24 µm fluxes of two reference debris disks, $\epsilon$ Eri and HD 107146. This approach is expected to result in a realistic mid-infrared brightness of the disk images, which is important for the simulated observations, but not in a good reproduction of the SED of the system over the whole wavelength range. The disks are placed at the distance of the reference debris disks. The images are convolved with a simulated PSF of MIRI at 25.5 µm (Fig. 8) produced by the software WebbPSF⁶. The resulting sensitivity of the observations is estimated using a sensitivity of 30 µJy (10$\sigma$, $\lambda$ = 25.5 µm, integration time = $10^4$ s)<fnref n="7" /> and scaling it with the square root of the actual integration time (Gardner et al. 2006).

In Fig. 8, one can clearly see the hexagonal structure and the peaks of the PSF at high order. Since the star contributes significantly to the flux in the mid-infrared, one has to perform accurate PSF subtraction to distinguish between real disk structures and these high order PSF structures. To evaluate the accuracy of stellar PSF subtraction necessary, the contribution of the PSF depending on the radial distance from the star has to be computed. Therefore, the PSF used to convolve the disk images is scaled to the total flux of the star in each of the two reference systems. The maximum of the PSF in one-pixel wide radial bins (pixel scale of the simulated disk images) is computed. The result is shown in Fig. 10. To obtain a reliable image of the disk, the contribution of the stellar PSF should be reduced by PSF subtraction (e.g., through observations of a reference star) to a level of ≈1/10 or less of the disk flux at a comparable distance from the star.

<footnote n="7">⁷ http://www.stsci.edu/jwst/science/data-simulation-resources; effective Aug. 2011.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** JWST/MIRI images at 25.5 µm. Leftmost column: images produced from the runs Ia and Ib. The $\epsilon$ Eri inner debris disk is used as reference. Columns two to five: images produced from the runs IIa, IIc, IId, and IIIc. The HD 107146 debris disk is used as reference. The contours in each image represent $3\sigma$, $5\sigma$, $10\sigma$, $20\sigma$, etc. (double every step) S/N levels for a total integration time of 10 s ($\epsilon$ Eri) and 900 s (HD 107146). The images are shown in a logarithmic stretch.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Radial distribution of the PSF structures at 25.5 µm as the maximum flux at a given radial distance. The total flux is scaled to the flux of the $\epsilon$ Eri stellar photosphere. To scale it to the HD 107146 stellar photosphere, one has to divide it by 42.3.</figcaption></figure>

## 6.2. Results

The results from the simulated observations of the disk only are shown in Fig. 9. In both the $\epsilon$ Eri and the HD 107146 case, one can see that the disks are spatially resolved. One can clearly distinguish between the different planet-disk configurations in both the face-on and the edge-on case. Only the results of run IIIc do not show a clearly resolved image of the disk. In contrast, the bright, unresolved inner part of the disk dominates the emission, which results in a bright image of the PSF overlayed on the outer disk structure.

The accuracy of stellar PSF subtraction necessary for the different runs is given in Table 5. This is the result of a comparison of the simulated disk observations (Fig. 9) and the results plotted in Fig. 10. We find that a PSF subtraction accuracy of 1% is sufficient to detect and spatially resolve all simulated debris disks considered in this study (significantly lower accuracy is needed in some cases). Since this uncertainty dominates the total uncertainty in the observations, very accurate photometry is required to scale the reference PSF to the correct level. This represents the minimum requirements to unequivocally de-

<table><caption>**Table 5.** PSF subtraction accuracy necessary for observations of the modeled debris disks.</caption><thead><tr><th>Run</th><th>ρ<sub>disk</sub><br>[arcsec]</th><th>F<sub>disk</sub><br>[mJy<br>beam]</th><th>ρ<sub>peak</sub><br>[arcsec]</th><th>F<sub>PSF</sub><br>[mJy<br>beam]</th><th>σ<br>[%]</th></tr></thead><tbody><tr><td>Ia</td><td>1.0</td><td>21.8</td><td>1.4</td><td>90.0</td><td>2.4</td></tr><tr><td>Ib</td><td>1.0</td><td>55.6</td><td>1.4</td><td>90.0</td><td>6.2</td></tr><tr><td>IIa (fo)</td><td>2.0</td><td>1.3</td><td>1.4</td><td>2.1</td><td>6.2</td></tr><tr><td>IIa (eo)</td><td>2.0</td><td>1.1</td><td>1.4</td><td>2.1</td><td>5.2</td></tr><tr><td>IIc (fo)</td><td>2.6</td><td>0.6</td><td>2.3</td><td>0.2</td><td>30.0</td></tr><tr><td>IIc (eo)</td><td>3.0</td><td>1.3</td><td>3.2</td><td>0.2</td><td>65.0</td></tr><tr><td>IId (fo)</td><td>3.0</td><td>0.2</td><td>3.2</td><td>0.2</td><td>10.0</td></tr><tr><td>IId (eo)</td><td>3.0</td><td>0.9</td><td>3.2</td><td>0.2</td><td>45.0</td></tr></tbody></table>

**Notes.** The quantity $\rho_{\text{disk}}$ denotes the (projected) distance from the star at which the disk flux is measured (the brightest structures of intent), $\rho_{\text{peak}}$ denotes the distance from the star of the closest peak of the stellar PSF$^8$, $F_{\text{disk}}$ and $F_{\text{PSF}}$ are the surface brightness of the disk and the PSF structures at the corresponding distances, and $\sigma$ denotes the accuracy at which the stellar PSF has to be subtracted to reduce the PSF structures to 1/10 of the disk structure (i.e., $1/10 \times F_{\text{disk}}/F_{\text{PSF}}$). Values are given for a face-on (fo) and edge-on (eo) orientation of the disk.

tect structures in the model debris disks without further considerations. However, a sophisticated approach might include proper rotation of the optics, so that edge-on seen disks are imaged between the bright hexagonal wings of the PSF – in regions where the PSF structures are less bright – increasing the depth of the observations. Furthermore, it might be possible in some cases to scale the reference PSF to the high order structures of the PSF in regions of the science observations where no significant signal from the disk is expected.

## 7. Conclusions

We demonstrated that planet-disk interaction may produce detectable structures in the dust distribution of debris disks. This depends on the configuration of the planet-disk system and

<footnote>$^7$ Of the two peaks of the PSF next to $\rho_{\text{disk}}$ in radial distance, the brighter one is considered.</footnote>
on the observing wavelength. The detected structures enable one to infer and characterize extrasolar planets in a range of mass and radial distance from the star unattainable by other techniques. In particular, detailed modeling of a combination of high-sensitivity, high spatial resolution observations at mid-infrared wavelengths and (sub-)mm wavelengths is able to put strong constraints on the configuration of the planet-disk system. We demonstrated that HST scattered-light observations are in most cases unable to unambiguously detect such structures, in particular in debris disks seen face-on. In contrast, both ALMA and the JWST will provide the sensitivity and resolution to detect and spatially resolve the spatial dust distribution in debris disks at a level of sensitivity and resolution that allow one to distinguish between different planet-disk configurations. However, we also demonstrated limitations of the instruments. ALMA is unable to reach both high sensitivity to surface brightness and high spatial resolution simultaneously, requiring a sophisticated observing strategy to optimize the outcome of planned observations. For a debris disk with typical shape of the SED, intermediate observing wavelengths around 1 mm and small to intermediate array extents are ideally suited to reach a high S/N and reasonable spatial resolution with ALMA. The highest sensitivity at reasonable spatial resolution is reached if the resolution of the observations is similar to the scale of the structures expected, as long as these structures are bright enough. Mid-infrared observations of debris disks with the JWST will be able to detect and spatially resolve dust in debris disks even at a distance of several tens of AU, where the emission from debris disks in this wavelength range is expected to be low. For such observations, stellar PSF subtraction with an accuracy of a few percent is necessary to unequivocally detect structures in the spatial distribution of the dust.

*Acknowledgements.* We thank J. P. Ruge for valuable discussions and important contributions to the consideration of phase noise. S. Ertel acknowledges financial support from DFG under contract WO 857/7-1 and from the French National Research Agency (ANR) for financial support through contract ANR-2010 BLAN-0505-01 (EXOZODI), and for general support by K. Ertel.

# Appendix A: Orbit integration of mass-less test particles

## A.1. Equation of motion

Because no interaction among the dust grains is considered, one can treat them as mass-less test particles independently of each other. Summarizing the different force terms, the equation of motion governing the dynamics of a dust particle at stlocentric position **r** and velocity **v** can be written as (Burns et al. 1979; Liou et al. 1995; Moro-Martín & Malhotra 2002)

$$ \mathbf{a} = \mathbf{a}_G + \mathbf{a}_R + \mathbf{a}_{PR} + \mathbf{a}_{SW} + \mathbf{a}_{Pl} \quad (A.1) $$

$$ = -\frac{GM_{\star}(1-\beta)}{r^3}\mathbf{r} - \frac{\beta(1+\xi)}{c}\frac{GM_{\star}}{r^2}\left(\frac{\dot{\mathbf{r}}}{r} + \mathbf{v}\right) + \sum_{j=1}^{p} GM_{\text{pl},j} \frac{\mathbf{r}_{\text{pl},j} - \mathbf{r}}{|\mathbf{r}_{\text{pl},j} - \mathbf{r}|^3} \quad (A.2) $$

$$ = -\frac{GM_{\star}}{r^3} \left\{ (1-\beta)\mathbf{r} + \frac{\beta(1+\xi)}{c} (\dot{\mathbf{r}} + \mathbf{r}\mathbf{v}) \right\} + \sum_{j=1}^{p} GM_{\text{pl},j} \frac{\mathbf{r}_{\text{pl},j} - \mathbf{r}}{|\mathbf{r}_{\text{pl},j} - \mathbf{r}|^3}, \quad (A.3) $$

where G is the gravitational constant, $M_{\star}$ is the stellar mass, P is the number of planets (1 in our case), and $M_{\text{pl},j}$ and $\mathbf{r}_{\text{pl},j}$ are the mass and position vector of planet j. The quantity $\beta$ is the ratio of radiation pressure and stellar gravity force (Burns et al. 1979). The strength of the corpuscular stellar wind (consisting mostly of protons) is prescribed by the parameter $\xi$, defined as the ratio of stellar wind drag to Poynting-Robertson drag. Following Gustafson (1994), we adopt the solar value of 0.35.

## A.2. Orbit integration

In this section, we only briefly describe the computational procedure; for a complete description we refer the reader to Kokubo et al. (1998). The orbits of the dust particles and that of the planet are integrated using a fourth-order Hermite scheme, a direct integration method using a predict-evaluate-correct (PEC) algorithm to solve the N-body problem (Makino 1991b; Makino & Aarseth 1992). For a particle i, the acceleration $\mathbf{a}_{0,i}$ and its first derivative with respect to time $\dot{\mathbf{a}}_{0,i}$ are used to predict the future position $\mathbf{r}_i(t)$ and velocity $\mathbf{v}_i(t)$ from the instantaneous values $\mathbf{r}_{0,i}$ and $\mathbf{v}_{0,i}$ at time $t_0 < t$.

The higher-order derivatives $\mathbf{a}_{0,i}^{(2)}$ and $\mathbf{a}_{0,i}^{(3)}$ are obtained directly by Hermite interpolation from $\mathbf{a}$ and $\dot{\mathbf{a}}$ at times $t_0$ and $t$ (Makino & Aarseth 1992). They are used to improve the accuracy of the force polynomial to refine the predicted position and velocity of the particle. For all simulations we iterate the sequence of (1) evaluation of the acceleration and its first derivative and (2) correction of the position and velocity with higher derivatives three times to increase the accuracy of the orbit integration.

For the choice of the individual time step, we use the prescription proposed by Aarseth (1985). It uses the acceleration and its first, second, and third time derivatives to compute the time step at which the integration for particle i should be advanced in time. The (continuous) time step values are then “quantized” to the nearest negative powers of 2 to further economize the computation (Makino 1991a). Then the computation proceeds in blocks of particles that are integrated simultaneously.

We have tested the accuracy of our numerical integration against the semi-analytical solutions of Wyatt & Whipple (1950) for initially circular and elliptical orbit of a dust particle. The temporal evolution of the orbital elements and the Poynting-Robertson decay time agree very well with Wyatt & Whipple (1950). The relative energy error of the orbit integration is typically below $10^{-6}$.

## References

Aarseth, S. J. 1985, in Multiple time scales, eds. J. U. Brackbill, & B. I. Cohen, 377

Apai, D., Janson, M., Moro-Martín, A., et al. 2008, ApJ, 672, 1196

Ardila, D. R., Golimowski, D. A., Krist, J. E., et al. 2004, ApJ, 617, L147

Augereau, J.-C., & Beust, H. 2006, A&A, 455, 987

Backman, D., Marengo, M., Stapelfeldt, K., et al. 2009, ApJ, 690, 1522

Boley, A. C., Payne, M. J., Corder, S., et al. 2012, ApJ, 750, L21

Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1

Corder, S., Carpenter, J. M., Sargent, A. I., et al. 2009, ApJ, 690, L65

Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531

Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89

Ertel, S., Wolf, S., Metchev, S., et al. 2011, A&A, 533, A132

Gardner, J. P., Mather, J. C., Clampin, M., et al. 2006, Space Sci. Rev., 123, 485

Golimowski, D. A., Ardila, D. R., Krist, J. E., et al. 2006, AJ, 131, 3109

Greaves, J. S., Holland, W. S., Moriarty-Schieven, G., et al. 1998, ApJ, 506, L133

Gustafson, B. A. S. 1994, Ann. Rev. Earth Planet. Sci., 22, 553
Hillenbrand, L. A., Carpenter, J. M., Kim, J. S., et al. 2008, ApJ, 677, 630

Holmes, E. K., Dermott, S. F., Gustafson, B. Å. S., & Grogan, K. 2003, ApJ, 597, 1211

Hughes, A. M., Wilner, D. J., Andrews, S. M., et al. 2011, ApJ, 740, 38

Jewitt, D., Luu, J., & Chen, J. 1996, AJ, 112, 1225

Kalas, P., Liu, M. C., & Matthews, B. C. 2004, Science, 303, 1990

Kalas, P., Graham, J. R., & Clampin, M. 2005, Nature, 435, 1067

Kalas, P., Graham, J. R., Chiang, E., et al. 2008, Science, 322, 1345

Kokubo, E., Yoshinaga, K., & Makino, J. 1998, MNRAS, 297, 1067

Krist, J. E., Stapelfeldt, K. R., Bryden, G., et al. 2010, AJ, 140, 1051

Krivov, A. V. 2010, Res. Astron. Astrophys., 10, 383

Kuchner, M. J., & Holman, M. J. 2003, ApJ, 588, 1110

Liou, J.-C., Zook, H. A., & Jackson, A. A. 1995, Icarus, 116, 186

Liseau, R., Eiroa, C., Fedele, D., et al. 2010, A&A, 518, L132

Liu, M. C. 2004, Science, 305, 1442

Makino, J. 1991a, PASJ, 43, 859

Makino, J. 1991b, ApJ, 369, 200

Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141

Marois, C., Macintosh, B., Barman, T., et al. 2008, Science, 322, 1348

Moro-Martín, A., & Malhotra, R. 2002, AJ, 124, 2305

Moro-Martín, A., & Malhotra, R. 2003, AJ, 125, 2255

Nilsson, R., Liseau, R., Brandeker, A., et al. 2010, A&A, 518, A40

Roccatagliata, V., Henning, T., Wolf, S., et al. 2009, A&A, 497, 409

Rodmann, J. 2006, Ph.D. Thesis, University of Heidelberg, http://www.ub.uni-heidelberg.de/archiv/6123

Stark, C. C., & Kuchner, M. J. 2008, ApJ, 686, 637

Stark, C. C., & Kuchner, M. J. 2009, ApJ, 707, 543

Swinyard, B., Nakagawa, T., Merken, P., et al. 2009, Exper. Astron., 23, 193

Udry, S., & Mayor, M. 2008, in Extreme Solar Systems, eds. D. Fischer, F. A. Rasio, S. E. Thorsett, & A. Wolszczan, ASP Conf. Ser., 398, 13

Vitense, C., Krivov, A. V., & Löhne, T. 2010, A&A, 520, A32

Weingartner, J. C., & Draine, B. T. 2001, ApJ, 548, 296

Wolf, S., & Hillenbrand, L. A. 2003, ApJ, 596, 603

Wolf, S., & Voshchinnikov, N. V. 2004, Comp. Phys. Commun., 162, 113

Wolf, S., Moro-Martín, A., & D'Angelo, G. 2007, Planet. Space Sci., 55, 569

Wyatt, M. C. 2006, ApJ, 639, 1153

Wyatt, S. P., & Whipple, F. L. 1950, ApJ, 111, 134

Wyatt, M. C., Dermott, S. F., Telesco, C. M., et al. 1999, ApJ, 527, 918