© ESO, 2011

1. Introduction

The vast majority of stars are born with a circumstellar disk. Observations show that in young clusters the percentage of stars with disks is close to 100%. The typical life times for disks appear to be around 6 Myr (Haisch et al. 2001). Circumstellar disks around young stars are believed to be the sites of planet formation. Solids in the disk grow to form larger objects, planetesimals, and eventually planets. Eventually, after the gas disappears, a planetary system will remain, containing planets and, if there is room in the dynamical structure, a debris belt.

Observations of protoplanetary disks at ever higher spatial resolution have provided promising indirect indications of ongoing planet formation (e.g. Thalmann et al. 2010; Jang-Condell & Kuchner 2010). Another strong indication of the presence of planets comes from structure observed in the dust that is still visible at later times in debris disks (Wilner et al. 2002; Wyatt 2003, 2008).

However, since most of the disks are difficult to resolve, studies also focus on understanding the system’s spectral energy distribution (SED). Comparing the shape of the SED to radiative transfer models gives information on the geometry of the disk and on the temperatures and particles found in the different parts of the circumstellar disk (Dullemond et al. 2007).

Over the past two decades, studying SEDs has improved our knowledge of circumstellar disks significantly and has made us aware of structures as flaring disks, puffed up inner rims, gaps, and holes (e.g. Dullemond et al. 2001, 2007). Holes and gaps are considered to be possible signposts of planet formation (Armitage 2010).

Recent studies of (pre-)transitional disks indicate a possible new structure. Several young stellar objects (YSOs) appear to show signs of an optically thin population of dust with a large scale height close to the star. These structures appear to be optically thin, but they still absorb and reprocess a significant fraction of the stellar light, which requires the dust to be distributed in a structure with high scale height, and not compressed toward the midplane of the disk. Such a large height of the structures is puzzling, since disks that are in hydrostatic equilibrium should be quite flat close to the star (Kenyon & Hartmann 1987; Chiang & Goldreich 1997). In this paper we study the possibility that the large height of the halo-like structure may be reached in a dynamical way by scattering planetesimals. This idea is inspired by the fact that planetesimals in the Solar System (in the asteroid belt and the Kuiper belt) can have substantial inclinations, caused by interactions with planets. If this mechanism is found to be capable of explaining the observed optically thin structures, their presence could be used to infer the existence of planets in these systems. This work will focus on the question of whether a inwardly migrating planet could be responsible for creating such an optically thin halo-like structure.

Section 2 discusses three YSOs that show this new structure. Section 3 discusses high inclinations in our own Solar System, and quantitatively formulates a mechanism that could cause highly-inclined orbits in a system with one planet. In Sect. 4 numerical simulations of the effect of a migrating planet on a planetesimal population are presented. The grinding down of the resulting planetesimal population to form small dust grains is covered in Sect. 5. Finally, a discussion and conclusions are formulated in Sects. 6 and 7.

2. Optically thin dust

Recently, a number of transitional disks have indications of possessing optically thin dust with large scale height: LkCa 15, HD 142527 and HD 163296. These systems are discussed below.

2.1. LkCa 15

LkCa 15 is a low-mass pre-main-sequence star located at about 140 pc in the Taurus star-forming region (Espaillat et al. 2007). Espaillat et al. (2007) compared different models to SEDs obtained with the Infra-red Array Camera (IRAC) aboard Spitzer to study LkCa 15. An optically thick outer disk is found to have an inner radius of 46 AU, containing about 0.1   M_⊙. The inner region of the disk is a bit more complex and Espaillat et al. provide two options which could explain the shape of the SED in the near-IR: an optically thick disk from 0.12 to 0.15 AU combined with 4 × 10^-11   M_⊙ of optically thin dust between 0.15 and 4 AU, or, 5 × 10^-11   M_⊙ of optically thin dust between 0.12 and 4 AU. The first option is said to give the better fit. Both scenarios include an optically thin part, but it is stressed no dust can exist further out than  ~5 AU since the contribution at 20   μm becomes too strong (Espaillat et al. 2007).

In 2009, Mulders et al. investigate the inner region of LkCa 15 in more detail and find an optically thick inner disk will influence the outer disk as well. Because the optically thick inner disk will obscure the outer disk, the outer disk needs to be blown up in the vertical direction in order for it to fit the observed SED. Alternatively, the inner disk could be optically thin, but this is only possible if the material is distributed over a large solid angle, virtually making it a halo. Mulders et al. (2010) used the 2D radiative transfer code MCMax (Min et al. 2009) to model the observed SED. Mulders et al. conclude it is hard to explain the large scale height of the inner dust material. One offered explanation is the violent scattering and disruption of planetesimals, which would end up in them colliding and filling the inner region with dust.

Espaillat et al. (2010) use new observations of the near-IR excess and find it to have the spectral shape of a blackbody, leading them to prefer the option of an optically thick inner disk and a blown-up outer disk. The inner disk is found to lie within 0.19 AU and the upper limit for the mass it contains is found to be 2 × 10^-4   M_⊙.

Whether the inner disk of LkCa 15 contains an optically thin dust component remains a matter of debate, but if optically thin dust is responsible for the near-IR flux, the scale height of the material has to be very large.

2.2. HD 142527

Another interesting system is that of HD 142527. The star was categorized as an F6III star by Houk (1978) and later classified as a Herbig star (Waelkens et al. 1996). HD 142527 is about 20 times as luminous as the sun, has a mass of 2.2 ± 0.3   M_⊙, and shows a huge IR-flux (F_IR = 0.92F_∗) which indicates there must be a considerable amount of circumstellar material (Verhoeff et al. 2011).

Verhoeff et al. used MCMax to model the circumstellar disk and find a model which accurately describes Spitzer and ISO IR spectra. A complicated structure is found for the system; an inner and an outer disk as well as an optically thin halo-like structure. The inner disk stretches from 0.3 to 30 AU and has a puffed up inner rim which shadows the rest of the disk. The small dust grains in the inner disk have a total mass of 2.5 × 10^-9   M_⊙. It is argued most of the other material is already locked up in planetesimals or even planets. The degree of sedimentation also indicates an advanced state of evolution (Verhoeff et al. 2011).

Then there is the optically thin halo, which stretches from 0.3 to 30 AU and has an estimated dust mass of 1.3 × 10^-10   M_⊙. The origin of this halo is unclear and again the dynamical scattering of planetesimals is offered as a speculative explanation.

There is a gap in the disk stretching from 30 to 130 AU. Verhoeff et al. explain why photo-evaporation is unlikely to be the cause of this gap and speculates about the presence of Jupiter or Earth-like planets as the origin of the gap. Lastly there is the outer disk, located between 120 and 200 AU, and with a rim height of 60 AU. Such a large rim height can be reached in hydrostatic equilibrium, but only if the disk is illuminated by the nearly unattenuated stellar radiation. A physical explanation for the high rim therefore requires the warm dust to remain optically thin and, therefore, to be spread out to large scale heights. The total mass in small dust grains is estimated to be 1.0 × 10^-3   M_⊙.

2.3. HD 163296

HD 163296 is an isolated Herbig Ae star located at about 122 pc, which shows a NIR excess indicating a circumstellar disk Hillenbrand et al. (1992). In 2010, Benisty et al. present long-baseline spectro-interferometric observations using the AMBER instrument on VLTI. They find that the inner rim of a disk alone cannot reproduce their observations in the NIR and need an additional emission component between 0.10 and 0.45 AU, with a temperature of 1600 K and an optical depth of  ~0.2.

Benisty et al. discusses whether hot gas could be held responsible for the emission between 0.10 and 0.45 AU but conclude that this is probably not the case. They argue that more detailed, non-LTE models that include the transition from optically thin to optically thick layers in a dust-free environment are needed to completely rule out the possibility of hot gas being the main contributor to the NIR flux.

The existence of refractory dust grains is explored as another possibility for the emission. An optically thin layer of dust grains is assumed to lie between an inner and outer radius, which are both inside the inner rim of the optically thick disk. The vertical optical depth of the dust is assumed to be proportional to 1/r, with r the distance to the star. Benisty et al. find the observed SED is reproduced best with the dust between 0.1 and 0.45 AU. The total mass in small dust between is calculated and found to be 8.7 × 10^-8   M_⊕ for Graphite grains of 0.05−2   μm, and 9.5 × 10^-7   M_⊕ for Iron grains of 0.02−2   μm in size.

3. Direct scattering to reach highly-inclined orbits

3.1. Highly-inclined orbits in our Solar System

How do we build a population of planetesimals on highly-inclined orbits from a nearly flat transitional disk? Hints for a mechanism can be found in the Kuiper belt in our own Solar System. The Kuiper belt can be split up into three dynamical classes (Jewitt & Luu 2000): the classical Kuiper belt, the resonant Kuiper belt, and the scattered Kuiper belt. Classical Kuiper belt Objects (KBOs) have low-eccentricity, low-inclination orbits beyond 40 AU. The resonant KBOs are trapped in a mean motion resonance (MMR) with Neptune, and the scattered ones have probably suffered a close encounter with Neptune in the past and show a wide spread in inclination. Scattered KBOs tend to have eccentric orbits with large semi-major axes. Studies on the distribution of inclinations of bodies in these three populations (Brown 2001; Gulbis et al. 2010) show they truly are distinct classes and it is argued that the inclinations in the three classes have different origins.

Observational studies show resonant KBOs to have inclinations up to 30° (Brown 2001; Gulbis et al. 2010). They are commonly believed to have been trapped in an MMR with Neptune when the planet was migrating outward, having their eccentricity and inclination raised during the process. This mechanism explains Pluto’s peculiar orbit (Malhotra 1993), and inclinations of up to  ~15° in general for bodies outside Neptune’s orbit (Malhotra 1995). However, only a small fraction of the resonant KBOs in the simulations have their inclinations raised in excess of 10% (Malhotra et al. 2000). Other mechanisms have been proposed to account for the observed inclinations, including sweeping secular resonances (Li et al. 2008), and stirring by large planetesimals (Morbidelli & Valsecchi 1997), but no definitive answer has been given.

Scattered KBOs show a wider spread in semi-major axes and their inclinations are usually greater than 10° with a peak around 20° (Gulbis et al. 2010). It is this broad distribution in orbital parameters that indicates that their orbits are probably the result of a close encounter with Neptune (Brown 2001). Hahn & Malhotra (2005) conclude from simulations that about 90% of the KBOs in the scattered disk might not have suffered a close encounter, but rather obtained their orbits from Neptune’s resonances during the migration epoch. Their model is unable to reproduce the observed abundance of KBOs with i > 15° however.

From studying the Kuiper belt, it appears there are 2 clear mechanisms that are able to increase an object’s inclination drastically; resonance trapping and close encounters. The degree of inclination-raising by resonance trapping is determined, and therefore limited by, the migration distance of the planet. Malhotra (1995) shows a migration Neptune over a distance of 7 AU to cause an increase in eccentricity of 0.25 and in inclination of about 10°. The inclination that can be reached by direct scattering is limited by the relative velocity between the planet and the object and the escape velocity at the planet’s surface, as we shall see in Sect. 3.2. It appears direct scattering is a more promising way to get bodies onto highly-inclined orbits, and it is this path that we will pursue in the rest of this work.

3.2. Direct planetesimal scattering

Suppose we have a planet of mass M_pl and radius R_pl orbiting a star of mass M_∗, on a circular orbit with Keplerian velocity ${\mathit{v}}_{\mathrm{K}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}{\mathit{a}}_{\mathrm{pl}}}$. We assume a planetesimal of mass m ≪ M_pl on an orbit with negligible inclination is approaching the planet with an impact parameter b and a relative velocity v_rel.

If b is small enough, the close encounter can be treated as a 2 body interaction, without having to consider the central star. During this interaction, the magnitude of the relative velocity is conserved but the direction is changed by a scattering angle of (e.g. Weidenschilling 1975) $\mathit{\theta }\mathrm{=}\mathrm{2}\arctan \left[\frac{\mathit{G}{\mathit{M}}_{\mathrm{pl}}}{{\mathit{v}}_{\mathrm{rel}}^{\mathrm{2}}\mathit{b}}\right]\mathrm{·}$(1)We normalize distances to the radius of the planet and velocities to the Keplerian velocity $\mathit{v̂}\mathit{i}\mathrm{=}\frac{{\mathit{v}}_{\mathit{i}}}{{\mathit{v}}_{\mathrm{K}}}\mathit{,}\mathit{b̂}\mathrm{=}\frac{\mathit{b}}{{\mathit{R}}_{\mathrm{pl}}}\mathrm{·}$(2)This allows us to rewrite Eq. (1) as follows $\mathit{\theta }\mathrm{=}\mathrm{2}\arctan \left[\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{v̂}\mathrm{esc}}{\mathit{v̂}\mathrm{rel}}\right)}^{\mathrm{2}}\mathit{b̂}-1\right]\mathit{,}$(3)where v_esc = (2GM_pl/R_pl)^1/2 is the escape velocity at the surface of the planet.

This equation shows that a scattering angle of 90° is reached for ${\left(\frac{{\mathit{v}}_{\mathrm{rel}}}{{\mathit{v}}_{\mathrm{esc}}}\right)}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{b̂}}\mathrm{·}$(4)Thus, for velocities ${\mathit{v}}_{\mathrm{rel}}\mathit{>}{\mathit{v}}_{\mathrm{esc}}\mathit{/}\sqrt{\mathrm{2}}$, an impact parameter smaller than the physical size of the planet is required. Such a planetesimal will crash into the planet. Due to gravitational focussing any planetesimal with an impact parameter smaller than $\mathit{b̂}\min \mathrm{=}\sqrt{\mathrm{1}\mathrm{+}{\left(\frac{\mathit{v̂}\mathrm{esc}}{\mathit{v̂}\mathrm{rel}}\right)}^{\mathrm{2}}}\mathit{,}$(5)will be lost in a physical collision with the planet (Safronov 1966).

Thus, for a planet to be able to significantly scatter a planetesimal out of the orbital plane, the relative velocity of the encounter needs to be about a factor 2 lower than the escape velocity at the planet’s surface.

Suppose the planet is able to scatter the planetesimal out of the orbital plane by an angle of 90°. In the frame of the planet, the planetesimal’s total velocity is now pointing perpendicular to the orbital plane. In the stellar frame however, the planetesimal is also moving with the planet at the Keplerian velocity, and we can calculate the new inclination using $\mathit{i}\mathrm{=}\arctan \left[\frac{{\mathit{v}}_{\mathrm{rel}}}{{\mathit{v}}_{\mathrm{K}}}\right]\mathrm{·}$(6)This simple calculation shows that an inclination of 45° can be reached for v_rel = v_K.

Summing up, planetesimals can only be scattered into highly-inclined orbits for a particular combination of v_rel, v_esc and v_K. This condition can be written as ${\mathit{v}}_{\mathrm{K}}\lesssim {\mathit{v}}_{\mathrm{rel}}\lesssim \frac{{\mathit{v}}_{\mathrm{esc}}}{\mathrm{2}}\mathrm{·}$(7)We can write the relative velocity of the collision as a function of orbital parameters of both objects. Recall that the planet moves on a circular orbit with semi-major axis a_pl. We define the orbit of the planetesimal by its eccentricity e, semi-major axis a, and inclination i relative to the planet. The relative velocity is then given by (Weidenschilling 1975) $\mathit{v̂}\begin{array}{}\mathrm{2}\\ \mathrm{rel}\end{array}\mathrm{=}\mathrm{3}\mathrm{-}\mathrm{2}\sqrt{\frac{\mathit{a}}{{\mathit{a}}_{\mathrm{pl}}}\left(\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}\right)}\cos \mathit{i}\mathrm{-}\frac{{\mathit{a}}_{\mathrm{pl}}}{\mathit{a}}\mathrm{·}$(8)Recalling i = 0 and taking a = a_pl, we need an eccentricity of e ≥ (3/4)^1/2 ≈ 0.87 to get v_rel ≥ v_K. Particles on orbits larger than that of the planet reach similar relative velocities for lower inclinations. However, if the particle is located inside the planetary orbit, for instance in the 1:2 MMR, it has a/a_pl ≈ 0.63 and only reaches v_rel = v_K for an eccentricity of 0.97.

Since this study focusses on inward planetary migration, most of the material will have a < a_pl. It is obvious from Eq. (8) that this means that the planetesimals need to be on (very) eccentric orbits. Section 4.3 will explore an explanation of these highly eccentric orbits.

The remainder of this section will calculate post-encounter inclinations in more detail and study the conditions for which planetesimals are scattered out of the system on hyperbolic orbits. For simplicity, we assume the orbits of the planetesimals to have e = 1, causing the velocities of the planet and the planetesimal to be perpendicular at the point of close encounter. Figure 1 shows the different quantities involved. Figure 1a depicts the stellar frame. In this case the comet velocity points in the negative x-direction, the planet moves in the direction of the positive y-axis, and the z-axis is perpendicular to the orbital plane.

Figure 1c shows the scattering in the frame of the planet. In this frame the planetesimal’s velocity points in the negative x′-direction, and has a size v_rel determined by $\mathit{v̂}\mathrm{rel}\mathrm{=}\sqrt{\mathrm{1}\mathrm{+}\mathit{v̂}\begin{array}{}\mathrm{2}\\ \mathrm{c}\end{array}}\mathit{.}$(9)The shaded area depicts the plane in which the scattering over an angle of θ occurs, which is the (y′ − z)-plane. Unit vectors x and x′ (and y and y′) lie in the orbital plane, but are rotated with respect to each other by an angle of $\mathit{\alpha }\mathrm{=}\arctan \mathrm{\left(}\mathit{v̂}\mathrm{c}\mathrm{\right)}\mathit{.}$(10)as depicted in Fig. 1b.

Fig. 1 a) Scattering of a planetesimal in the stellar frame. The planetesimal comes in along the x-axis and the planet moves in the positive y-direction. The z-axis is perpendicular to the orbital plane. b) Unit vectors x and x′ (and y and y′) lie in the orbital plane, but are rotated with respect to each other by an angle α. c) The scattering in the frame of the planet. The planetesimal comes in along the x′-axis, with impact parameter b. It is scattered in the (y′ − z)-plane over an angle θ.

As a planetesimal comes in (along the x′-axis) with impact parameter , the scattering angle θ is determined by Eq. (3). We repeat that a comet approaching with will have a physical collision with the planet.

For the post-encounter comet¹ inclination (i) and velocity (v_post), it is not just the size of the impact parameter that matters, but also the direction of this offset with respect to the motion of the planet. To quantify this direction we use φ. If the planetesimal is aiming at a point slightly behind the planet (φ = π), it will be scattered in the planet’s direction of motion, and therefore accelerated. If the planetesimal is aiming at a point directly above or below the planet (φ =  ± π/2), it is scattered out of the orbital plane.

If we split the post-encounter velocity of the comet into a component in the z-direction and one in the orbital plane, and respectively, we can calculate the total velocity to be $\mathit{v̂}\mathrm{post}\mathrm{=}\sqrt{\mathit{v̂}\begin{array}{}\mathrm{2}\\ \mathrm{post}\mathit{,}\hspace{0.17em}\mathit{xy}\end{array}\mathrm{+}\mathit{v̂}\begin{array}{}\mathrm{2}\\ \mathrm{post}\mathit{,}\hspace{0.17em}\mathit{z}\end{array}}\mathit{,}$(11)and the post-encounter inclination using $\mathit{i}\mathrm{=}\arctan \left[\frac{\mathit{v̂}\mathrm{post}\mathit{,}\hspace{0.17em}\mathit{z}}{\mathit{v̂}\mathrm{post}\mathit{,}\hspace{0.17em}\mathit{xy}}\right]\mathrm{·}$(12)Since the escape velocity from the central star is equal to $\sqrt{\mathrm{2}}{\mathit{v}}_{\mathrm{K}}$, comets with $\mathit{v̂}{}_{\mathrm{post}}\mathit{>}\sqrt{\mathrm{2}}$ will leave the system on a hyperbolic orbit.

For a given combination of , , and φ, the resulting inclination and total velocity are uniquely defined and can be calculated using the following prescription:

• 1.

  calculate and from Eqs. (9) and (5);

• 2.

  determine angles α and θ;

• 3.

  in the planet frame, rotate using θ and φ;

• 4.

  switch to the stellar frame by adding the planet’s velocity under the correct angle α;

• 5.

  calculate using Eq. (11) and check if the particle is ejected;

• 6.

  find the comet’s post-encounter inclination from Eq. (12).

Figure 2 shows the post-encounter inclination of a comet coming in with as a function of and φ for four different values of . The regions where the comet is lost because of ejection from the system or a collision with the planet are coloured white.

Fig. 2 The inclination of a comet after an encounter at with a planet with , as a function of impact parameter and φ. The white regions mark the combinations of and φ where $\mathit{v̂}{}_{\mathrm{post}}\mathit{>}\sqrt{\mathrm{2}}$. Regions where b < bmin are also coloured white.

For a given system, ejection is more likely to happen for comets with relatively small impact parameters, and φ ~ π. The latter is caused by these comets coming in behind the planet, and thus being accelerated in the stellar frame during the close encounter. The value of φ that is ideal for reaching high inclinations appears to lie between 0 and π/2, because the comet not only gains velocity in the z-direction, but it is also slowed down in the xy-plane.

The planet with the lowest mass can be seen not to be able to eject bodies, nor can it scatter them to orbits with inclinations  >5°. As the escape velocity at the planet’s surface is increased, three things happen; 1) increases; 2) the region where particles are ejected grows larger (both in φ and in ); and 3) the region where high inclinations are reached shifts towards higher impact parameters. Figure 2 confirms Eq. (7), as an escape velocity higher than the relative velocity is needed to scatter bodies to highly-inclined orbits.

3.3. Eccentricity pumping and migration

The relative velocity used in Fig. 2 can only be reached for very eccentric planetesimal orbits. A possible explanation for high eccentricities is resonance trapping, see for instance Malhotra (1993, 1995); Malhotra et al. (2000). In our own Solar System, the trapping and consequent dragging along of Plutinos in Neptune’s 3:2 resonance, has raised their inclinations from nearly 0 to up to 30° (see Sect. 3.1).

A planetesimal orbiting a star is in an MMR with a planet if their mean motions are locked, and the resonant argument is librating (see Murray & Dermott 1999, for a detailed treatment of the resonant argument). The locations of the mean motion resonances of a planet orbiting a star at a_pl are ${\mathit{a}}_{\mathrm{(}\mathit{p}\mathrm{+}\mathit{q}\mathrm{)}\mathrm{:}\mathit{p}}\mathrm{=}{\mathit{a}}_{\mathrm{pl}}{\left(\frac{\mathit{p}\mathrm{+}\mathit{q}}{\mathit{p}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}\mathit{,}$(13)with p and q integers, and the ratio of an MMR depicting the ratio of the orbital periods, so that a 1:2 resonance lies inside the planet’s orbit while a 2:1 resonance lies outside. Resonances with |q| = 1 are called first-order resonances and are often described with j instead of p.

When the planet is migrating, the locations of the MMRs will change and bodies trapped within a resonance can be dragged along. This will not only change the dragged-along body’s semi-major axis, but also various other orbital parameters. Malhotra (1993) showed Pluto’s peculiar orbit could be explained by it being stuck in Neptune’s 3:2 resonance, as Neptune migrated out by  ~7 AU.

Wyatt (2003) provides us with an expression for the eccentricity of a planetesimal orbit, that has been dragged in a (p + q):p resonance by a planet ${\mathit{e}}^{\mathrm{2}}\mathrm{=}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{+}\left(\frac{\mathit{q}}{\mathit{p}\mathrm{+}\mathit{q}}\right)\ln \frac{\mathit{a}}{{\mathit{a}}_{\mathrm{0}}}\mathit{,}$(14)where a₀ and a are the initial and current semi-major axes of planet. Supposing a planetesimal gets captured in the 2:1 resonance with e₀ = 0, it would have to be dragged to a/a₀ ⋍ 5 to reach an eccentricity of 0.9.

Resonance capture is not guaranteed however. As the resonance of a planet moves in (or out), because of the planet migrating, and it encounters a planetesimal, the probability of capture depends on both eccentricities, the migration rate, and planet mass (Mustill & Wyatt 2011). Mustill & Wyatt conclude resonance capture is guaranteed for low eccentricity and low migration rates, impossible at fast migration rates and low eccentricities, and possible for fast migration rates and high eccentricities. The transitions between these regimes depend on the planet mass. Furthermore, the critical migration rate for capture increases with j, so if the particle survives the passing of the 2:1 resonance, it may still be captured in another resonance closer to the planet with higher j. It is also pointed out that planetesimals that are not caught, experience a small jump in e up or down depending on the migration rate. Mustill & Wyatt (2011) calculate capture probabilities as a function of migration rate for different MMRs, migration rates, and planetary masses. For instance, a Neptune-mass planet is found to capture all bodies in its 3:2 resonance for migration rates  ≲ 4 × 10^-6   AU   yr^-1.

Resonance trapping followed by eccentricity-pumping appears to be needed to ensure high relative velocities. Another argument for an inwardly migrating planet is the fact it can collect and carry material inwards. A stationary planet can only scatter planetesimals in its direct vicinity. Since the optically thin structures are located close to their central stars and have considerable mass in small dust grains (see Sect. 2) this build-up of material by an inwardly migrating planet is crucial.

Assume a planetesimal disk with a mass distribution equal to the mass in solids in a minimum mass Solar Nebula (MMSN). The surface density of solids in a MMSN disk is given by (Hayashi 1981) ${\mathrm{\Sigma }}_{\mathrm{s}}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{0}}{\left(\frac{\mathit{r}}{\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\hspace{0.17em}\mathrm{g}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-2}\mathit{.}$(15)Σ₀ = 7.1 for r < 2.7 AU because ices cannot survive in the hot inner regions. Outside the snowline Σ₀ = 30.

A stationary planet of Neptune-mass on a circular orbit at 5 AU can scatter planetesimals within a region the size of its Hill sphere (e.g. Armitage 2010) $\mathrm{\Delta }\mathit{r}\mathrm{=}{\mathit{R}}_{\mathrm{H}}\mathrm{=}\mathit{a}{\left(\frac{{\mathit{M}}_{\mathrm{pl}}}{\mathrm{3}{\mathit{M}}_{\mathrm{\ast }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{3}}\mathrm{·}$(16)The mass of planetesimals within a ring of this size is given by $\mathrm{\Delta }\mathit{m}\mathrm{=}\mathrm{2}\mathit{\pi a}{\mathrm{\Sigma }}_{\mathrm{p}}\mathrm{\Delta }\mathit{r}\mathit{.}$(17)Using Σ_p = 30   g   cm^-2, we find Δm ≃ 4.6   M_⊕.

Assume now that the planet migrated from 30 to 5 AU at such a rate, it collected all material in between. We can integrate Eq. (15) to find Δm ≃ 46   M_⊕. This is an increase by a factor of 10 with respect to the mass scattered by a stationary planet at 5 AU. The difference becomes even bigger when looking at planets closer to the star, since Δm in Eq. (17) is proportional to a².

4. Simulations

4.1. Introduction

The scenario for creating the highly-inclined orbits is now as follows. A planet with v_esc ≳ v_K and a planetesimal population start out on (nearly) circular orbits and in the orbital plane (Fig. 3a). Due to for example the presence of gas the planet starts migrating inwards, and collecting material in mean motion resonances. Planetesimals that are stuck in an MMR, have their eccentricities raised significantly. Figure 3b shows the system after some migration has taken place. All planetesimal orbits in this figure have the same semi-major axis (they are presumed stuck in the same MMR) but differ in eccentricity. As the planet migrates even further, the increasing eccentricities cause the orbits of the planetesimals to be planet-crossing (Fig. 3c). At this point, the inclinations are only a couple of degrees. The relative velocity between the planet and the planetesimals grows together with the eccentricity of the planetesimal orbit and at some point this will lead to violent scattering events. As the planet has nearly reached the star it will leave behind a scattered population of planetesimals, which show (depending on the planet’s escape velocity) a large spread in inclinations, Fig. 3d. This population of scattered planetesimals, over time, will collide and grind down to form the small dust that causes the optical depth. In this section numerical simulations are presented which show the resonance capture, eccentricity-pumping, and scattering of planetesimals happening simultaneously.

Fig. 3 Illustration of the mechanism to reach a planetesimal population with high orbits, as explored in this paper. a) both the planet orbit (dotted) and planetesimal orbits (solid) start out in the orbital plane. b) As the planet has migrated in a some distance, the planetesimals are stuck in an MMR. As they are dragged along their orbits become increasingly eccentric. c) At some point the planetesimal orbits become so eccentric, the bodies become planet-crossing. As the eccentricity rises even further, the relative velocity between the planet and the planetesimals grows. d) The high relative velocity compared to the local Keplerian velocity enables the planet to violently scatter the planetesimals throughout the system, to orbits with high inclinations.

4.2. Setup

For the integrations of the migrating planet the software package Swifter is used². Swifter is written by Kaufmann and based upon the Swift package developed by Levison & Duncan. Within Swifter, the user is given the choice between 7 integrating schemes, including Wisdom-Holman Mapping (Wisdom & Holman 1991), the Regularized Mixed Variable Symplectic (RMVS) integrator (Levison & Duncan 1994), and the Symplectic Massive Body Algorithm (SyMBA) (Duncan et al. 1998; Levison & Duncan 2000). In addition, Swifter also includes a Bulirsch-Stoer integrator (Press et al. 1992). In this work SyMBA was used most of the time, since it can handle close encounters correctly and is significantly faster than the Bulirsch-Stoer method. The results of the SyMBA integrator are checked by comparing a few cases to a Bulirsch-Stoer integration with a very small timestep.

During the integrations, a planet is pushed to migrate inwards from 30 AU, while a population of 100 massless planetesimals orbits the star at circular orbits between 10 and 15 AU. The planetesimals all start out with e = i = 0, whereas the planet is given a small inclination of 0.72°, to ensure φ can take on all values between 0 and 2π, and an eccentricity of 0.012. The planet mass is varied between 0.1 and 20 Neptune masses.

For simplicity, the migration is not done in a self-consistent fashion. We want to focus on the effect of the migrating planet on the planetesimal population, to see if the resonance capture, eccentricity pumping, and scattering can take place. A self-consistent migration would vastly increase the complexity of the problem. Instead, the planet is pushed within SyMBA to have it migrate at a constant rate ξ, which is varied between 7.3 × 10^-6   AU   yr^-1 and 7.3 × 10^-4   AU   yr^-1. These migration rates are based upon estimates for migration rates and time scales for type I migration (Tanaka et al. 2002), type II migration (Armitage 2010), and migration through a planetesimal disk (Armitage 2010). Since the migration in this work is artificial, the system does not contain any gas or massive planetesimals and, therefore, does not distinguish between the various types of migration. Section 6 will discuss what effect the presence of gas or other drivers for migration might have on the results.

For SyMBA to be able to handle close encounters between massless bodies and planets, a Hill-sphere radius, R_H, needs to be specified for every planet. This radius determines when SyMBA will shorten the integration timestep. Since the Hill-sphere radius is proportional to a planet’s semi-major axis (Eq. (16)), it will decrease in time for a planet that migrates inwards. For simplicity, we will treat R_H as a constant. This results in the size of the Hill sphere being overestimated as the planet migrates in, slowing down the integrations somewhat.

The system is evolved long enough for the planet to pass through the planetesimal population, and end up close to the star,  <1 AU, where it has a negligible effect on the planetesimals. The timestep for the integrations ranges between 0.5 days and 10^-1 yr, depending on the smallest perihelion distance in the simulation, making sure all orbits are resolved. Within SyMBA, the timestep is decreased when a planetesimal enters the Hill sphere of the planet. During the simulations, planetesimals are removed if they:

• come within a distance r_min from the origin. We use r_min = 0.05 AU;

• have a > r_max. In our simulations r_max = 1000 AU is used;

• come within R_pl. In this work R_pl = R_Nep is used.

The choice for r_min comes from the fact that very eccentric orbits with perihelion distances smaller than 0.05 AU are hard to resolve without using a very small timestep, and will lead to non-physical behaviour. Furthermore, bodies that come this close to the central star are likely to be evaporated.

4.3. Resonance capture

An example of resonance capture in the simulations is illustrated in Fig. 4. The figure shows a Neptune-mass planet being pushed artificially to migrate in from 30 AU at a rate of 7.3 × 10^-5   AU   yr^-1. A massless asteroid is orbiting the solar-mass star at 13 AU on an initially circular orbit. The asteroid survives the passing of the 1:2 (j = 1) MMR but is eventually captured by the 3:4 (j = 3) MMR. For slower migration rates (or higher planet masses) the asteroid gets trapped in the 1:2 resonance already, while for faster migration rates (or lower planet masses) the asteroid only gets captured very close to the planet.

Figure 5 shows the evolution of a massless body as it is dragged along in a migrating planet’s 1:2 resonance. The integration has been performed using two different integrators, to check the consistency of SyMBA for highly eccentric orbits. The test particle is captured at t ≈ 0.13 Myr and eventually lost to the central star after about 0.32 Myr. During this time, the 10   M_Nep mass planet has travelled  ~14 AU. The results of the SyMBA and Bulirsch-Stoer integrators are very similar, and in both cases the massless body is eventually lost to the central star when it reaches e ≈ 0.99. The expected increase in eccentricity, as predicted by Eq. (14) is also plotted and seen to agree well with the simulations. The discrepancy between Eq. (14) and the simulations towards higher e comes from the fact Eq. (14) is a first-order approximation.

Fig. 4 A Neptune-mass planet is migrating in with ξ = 7.3 × 10-5   AU   yr-1. The asteroid survives the passing of the 1:2 resonance but is captured in the 3:4 resonance and dragged inward.

Fig. 5 Evolution of the eccentricity of a massless body stuck in the 1:2 internal mean motion resonance of a 10   MNep planet migrating inward at a rate of 7.3 × 10-5   AU   yr-1. The timestep for both the SyMBA and Bulirsch-Stoer integration was set at 0.5 days. The dotted line shows the expected increase in eccentricity as predicted by Eq. (14). The test particle is captured at t ≈ 0.13 Myr and eventually lost to the central star after about 0.32 Myr.

4.4. Inclinations

Figure 6a shows the evolution of the planetesimal swarm as a function of time. The left panels show eccentricity, the right panels inclination, and time passes from top to bottom. The large circle shows the location of the planet and the small circles represent massless planetesimals. Planetesimals above the solid lines are on a planet-crossing orbit. The dotted lines show the 1:2 and 2:3 MMR, where the 1:2 MMR is furthest from the planet. The planet is 0.1   M_Nep and migrating at a rate of 7.3 × 10^-5   AU   yr^-1. The planetesimals start out between 10 and 15 AU with e = i = 0. After 0.15 Myr (top panels), both resonances have reached the planetesimal population, and both appear unable to capture planetesimals. The passing of the resonances through the swarm of planetesimals is seen to have a very small effect on their orbital parameters. After 0.3 Myr the planet itself has migrated through the population of massless bodies, and has scattered most of them via direct encounters. Since the planetesimal orbits were not excited before the encounter, the relative velocity remained small, and the scattering events proved unable to get the planetesimal inclinations higher than a couple of degrees.The resulting population, shown after 0.4 Myr in the bottom panels, is very similar to the one after 0.3 Myr, i.e., the planet has little effect on the bodies once it has migrated inside the population. During this simulation, 2% of the planetesimals is accreted by the migrating planet.

Fig. 6 Evolution of the eccentricities (left) and inclinations (right) in a planetesimal swarm as a planet (large circle) migrates inwards from 30 AU at a rate of 7.3 × 10-5   AU   yr-1. The planet mass used are 0.1   MNep for a), and 2   MNep for b). Planetesimals (small circles) above the solid lines are on a planet crossing orbit, and the dotted lines show the location of the 1:2 and 2:3 MMR.

If we increase the planet mass, this will lead to a high capture probability. Figure 6b shows a similar calculation, but for a planet mass of 2   M_Nep. After 0.15 Myr the effect of the enhanced planet mass are already visible as the passing of the 1:2 resonance raises the eccentricities of the planetesimals slightly, while the 2:3 resonance appears to be able to capture bodies, and subsequently drag them along. This becomes even more visible after 0.3 Myr, where a significant fraction of the planetesimals are stuck in the 2:3 MMR. As these bodies are dragged along, their eccentricities are raised as described by Eq. (14). As the eccentricity of a planetesimal goes up, it will eventually reach a planet-crossing orbit, allowing it to have a close encounter with the planet. This has already happened for the bodies no longer in the 2:3 resonance. Since the close encounters happen when the planetesimals are on eccentric orbits, the relative velocities involved are higher, and the inclinations resulting from the scattering events are much greater than for Fig. 6a. During this simulation, 1% of the bodies were accreted by the central star, and 1% by the migrating planet.

If we increase the planet mass even further, the capture probability reaches almost unity, and the vast majority of planetesimals get stuck in a resonance far away from the planet. As they are dragged along and their eccentricities go up, they are protected from a close encounter with the planet for some time, since they are so far away. It turns out, almost all these bodies will reach e ≫ 0.9, and are accreted by the central star before they have time to be scattered by the migrating planet. An example of this scenario is shown in Fig. 5. This means, that in the systems where you would expect to find the highest relative velocities (and therefore the highest inclinations) the majority of the planetesimals are actually lost to the central object. A way to get around this, might be to put an additional, stationary planet close to the central star, which will scatter the highly excited planetesimals before they come too close to the central object. This planet could not be too massive, as this would lead to it ejecting the planetesimals on hyperbolic orbits. In this study we have restricted ourselves to single-planet systems, but the study of stellar systems with multiple planets, one of which is migrating, might be very interesting indeed. Another way to protect the planetesimals from accretion by the central star would be to introduce a gas drag, which would dampen the eccentricities. The effect of gas of the results is further discussed in Sect. 6.

From Sect. 4.3 we know that the capture probability also depends on the migration rate of the planet. We have therefore conducted a small parameter study, where the migration rate and planet mass have been varied. In total we carried out 12 simulations with the planet-mass ranging from 0.1−20   M_Nep and ξ between 7.3 × 10^-4−7.3 × 10^-6   AU   yr^-1. Table 1 shows the characteristics of planetesimal populations after migration has taken place, as well as the fraction of planetesimals which have been lost during the integration. Both higher planet mass and lower migration rate appear to ensure high average inclinations. This is the case because; 1) massive planets are more capable of scattering planetesimals into highly-inclined orbits (see Sect. 3.2); and 2) relatively massive planets and low migration rates increase the probability of resonance capture, allowing for more eccentricity-pumping and a higher relative velocity between the planetesimal and the planet during the close encounter. However, if the combination of migration rate and planet mass causes the planetesimals to get stuck in a resonance too far, they are likely to be lost to the central star before having a close encounter with the migrating planet.

5. Dust production

Section 4 has left us with populations of planetesimals and constraints on their orbital parameters. Section 5.1 will calculate the time scales on which these planetesimals should collide and form dust, and Sect. 5.2 will estimate the total mass and surface area of the dust size distribution resulting from the collisional cascade.

5.1. Collision times

Suppose we have N planetesimals of radius R, mass m and a geometrical cross-section for collisions of σ_coll = 4πR², flying around in some volume V. We can define the sweeping time as ${\mathit{t}}_{\mathrm{s}}\mathrm{=}\frac{\mathit{V}}{{\mathit{v}}_{\mathrm{coll}}{\mathit{\sigma }}_{\mathrm{coll}}}\mathit{,}$(18)where v_coll is the collision velocity between the planetesimals. Since the populations from Table 1 show pretty eccentric orbits, the collision velocity will vary a lot between collisions, but for the purpose of this calculation we estimate it as ${\mathit{v}}_{\mathrm{coll}}\mathrm{=}\mathit{\eta }{\mathit{v}}_{\mathrm{K}}\left(\mathrm{⟨}\mathit{a}\mathrm{⟩}\right)\mathit{,}$(19)with η between 0.5 and 2. For two nearly circular orbits η is determined by the vertical component of the individual velocities (Dominik & Decin 2003). But since the orbits in our populations are eccentric and randomly oriented, η will lie closer to $\sqrt{\mathrm{2}}$.

If we assume the volume is wedge-shaped and lies between an inner and an outer radius, r_in and r_out, and has a normalized height of h = H/r, we can write it as $\mathit{V}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}\mathit{h}\mathrm{(}{\mathit{r}}_{\mathrm{out}}^{\mathrm{3}}\mathrm{-}{{\mathit{r}}_{\mathrm{in}}^{\mathrm{3}}}^{\mathrm{)}}\mathit{,}$(20)where h = tani_max.

We have to think carefully about r_in and r_out. One might insert the minimum and maximum values of the semi-major axis found in the resulting populations, but this is not correct since most comets are on pretty eccentric orbits. On the other hand, choosing r_out to be equal to a_max(1 + e_max) will increase the volume drastically even though only one planetesimal is able to get to r_out. It makes more sense to define something like $\begin{array}{ccc}{\mathit{r}}_{\mathrm{out}}& \mathrm{=}\mathrm{(}\mathrm{⟨}\mathit{a}\mathrm{⟩}\mathrm{+}{\mathit{\sigma }}_{\mathit{a}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}\mathrm{⟨}\mathit{e}\mathrm{⟩}\mathrm{)}\mathit{,}& \\ {\mathit{r}}_{\mathrm{in}}& \mathrm{=}\mathrm{(}\mathrm{⟨}\mathit{a}\mathrm{⟩}\mathrm{-}{\mathit{\sigma }}_{\mathit{a}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{⟨}\mathit{e}\mathrm{⟩}\mathrm{)}\mathit{,}& \end{array}$where  ⟨ a ⟩  and  ⟨ e ⟩  can be taken from Table 1. We approximate ${\mathit{\sigma }}_{\mathit{a}}\mathrm{~}\sqrt{\mathrm{⟨}\mathit{a}\mathrm{⟩}}$, which is found to work well for most populations in Table 1.

The collision time for the comets will become ${\mathit{t}}_{\mathrm{coll}}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{s}}}{\mathit{N}}\mathit{,}$(23)with the number of bodies being constrained by the total mass of the population M$\mathit{N}\mathrm{=}\frac{\mathit{M}}{\mathit{m}}\mathrm{=}{\left(\frac{\mathit{R}}{{\mathit{R}}_{\mathrm{\oplus }}}\right)}^{-3}{\left(\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{\oplus }}}\right)}^{-1}\left(\frac{\mathit{M}}{{\mathit{M}}_{\mathrm{\oplus }}}\right)\mathrm{·}$(24)where ρ denotes the density of the planetesimal material. We can now calculate the collision time scale for different populations. A promising case for creating a halo-like structure appears to be M2, which has high inclinations but is still located relatively close to the star. From Table 1 we find  ⟨ a ⟩  = 8.61 AU,  ⟨ e ⟩  = 0.538 and i_max = 35.25°. For this specific configuration of comets, the collision time becomes ${\mathit{t}}_{\mathrm{coll}}^{\mathrm{\left(}\mathit{M}\mathrm{2}\mathrm{\right)}}\mathrm{=}\mathrm{3.7}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}{\left(\frac{\mathit{\eta }}{\mathrm{1.41}}\right)}^{-1}\left(\frac{\mathit{R}}{\mathrm{km}}\right)\left(\frac{\mathit{\rho }}{\mathrm{0.5}{\mathit{\rho }}_{\mathrm{\oplus }}}\right){\left(\frac{\mathit{M}}{\mathrm{50} {\mathit{M}}_{\mathrm{\oplus }}}\right)}^{-1}\hspace{0.17em}\mathrm{yr}\mathit{.}$(25)The collision time scale for the population in M2 is comparable to the migration time itself. This means a number of asteroids might suffer collisions during the migration itself, which we have not considered. The collision time is also an indication of how long the dust can survive in the system, as new dust is continuously formed during this period.

5.2. Dust size distribution

When planetesimals collide at sufficiently high velocities, they will fragment and create smaller bodies. The distributions of these bodies has the shape of a power law $\mathit{f}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{\propto }{\mathit{m}}^{\mathrm{-}\mathit{q}}\mathit{,}$(26)with q = 11/6 (Dohnanyi 1968). Assuming spherical particles we can rewrite this as a size distribution $\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{R}}{\mathit{R}}^{\mathit{\gamma }}\mathit{,}$(27)with γ = 2 − 3q =  −3.5. The proportionality factor f_R depends on the material properties of the initial population and is calculated by Dominik & Decin (2003) as ${\mathit{f}}_{\mathit{R}}\mathrm{=}\mathit{N}\sqrt{\frac{\mathrm{2}\mathit{m}{\mathit{\sigma }}_{\mathrm{coll}}}{\mathrm{(}\mathrm{4}\mathit{\pi }\mathit{/}\mathrm{3}\mathrm{)}\mathit{\rho }{\mathit{ϵ}}_{\mathrm{0}}}}\mathrm{=}\mathit{N}{\mathit{R}}^{\mathrm{5}\mathit{/}\mathrm{2}}\sqrt{\frac{\mathrm{8}\mathit{\pi }}{{\mathit{ϵ}}_{\mathrm{0}}}}\mathit{,}$(28)with ${\mathit{ϵ}}_{\mathrm{0}}\mathrm{=}\mathrm{-}\frac{\mathit{\pi }{\mathit{ϵ}}^{\mathit{\gamma }\mathrm{+}\mathrm{1}}}{\mathit{\gamma }\mathrm{+}\mathrm{1}}\mathit{,}$(29)and $\mathit{ϵ}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{coll}}^{\mathrm{2}}}{\mathrm{4}\mathit{S}}\mathrm{-}\sqrt{\frac{{\mathit{v}}_{\mathrm{coll}}^{\mathrm{4}}}{\mathrm{16}{\mathit{S}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{v}}_{\mathrm{coll}}^{\mathrm{2}}}{\mathrm{2}\mathit{S}}}\mathrm{-}\mathrm{1}\mathit{,}$(30)where S is the binding energy of the material. For asteroid type bodies S = 10⁷   erg   g^-1 (Dominik & Decin 2003).

5.3. Total mass and optical depth

The size distribution can be integrated to yield the total mass ${\mathit{M}}_{\mathrm{dust}}\mathrm{=}{\mathrm{\int }}_{{\mathit{R}}_{\min }}^{{\mathit{R}}_{\max }}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{R}}^{\mathrm{3}}\mathit{\rho f}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{d}\mathit{R,}$(31)in dust grains with radii between R_min and R_max.

The smallest size, R_min, is determined by comparing the gravitational force acting on the dust grain with the radiation pressure $\mathit{\beta }\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{=}\frac{{\mathit{F}}_{\mathrm{rad}}}{{\mathit{F}}_{\mathrm{grav}}}\mathrm{=}\frac{\mathrm{3}{\mathit{L}}_{\mathrm{\ast }}{\mathit{Q}}_{\mathrm{pr}}}{\mathrm{16}\mathit{\pi cG}{\mathit{M}}_{\mathrm{\ast }}\mathit{R\rho }}\mathit{,}$(32)with c the speed of light and Q_pr the radiation pressure coefficient ${\mathit{Q}}_{\mathrm{pr}}\mathrm{\equiv }{\mathit{Q}}_{\mathrm{abs}}\mathrm{+}{\mathit{Q}}_{\mathrm{scat}}\mathrm{(}\mathrm{1}\mathrm{-}\mathrm{⟨}\cos \mathit{\alpha }\mathrm{⟩}\mathrm{)}\mathit{,}$(33)where α denotes the scattering angle. For isotropic forward scattering,  ⟨ cosα ⟩  = 1 (Burns et al. 1979). Particles with β > 1/2 are ejected from the system.

Assuming Q_pr = 1, a dust grain with ρ = 2.76   g   cm^-3 orbiting a sun-like star reaches the critical value of β = 1/2 for a radius of about 0.4   μm. From here on on we will assume 10   μm to be the maximum size of a dust particle still contributing significantly to the optical depth.

Again considering M2, we have v_coll = 14.4   km   s^-1. Together with S = 10⁷   erg   g^-1 this gives ϵ₀ = 1.90 × 10⁴. If we take the total planetesimal mass to be 50   M_⊕, and the material density of the bodies to be 2.76   g   cm^-3, we find the total mass in dust grains between 0.4 and 10   μm to be M_dust = 8.7 × 10^-10   M_⊙. This mass compares well to the mass in small dust grains quoted in Sect. 2, but depends heavily on parameters such as the total planetesimal mass, and is only intended as an order of magnitude estimate.

The optical depth along a line of sight dr can be written as $\mathrm{d}\mathit{\tau }\mathrm{=}\mathrm{-}\mathit{\kappa \rho }\mathrm{d}\mathit{r}\mathrm{=}\mathit{n\sigma }\mathrm{d}\mathit{z,}$(34)with κ the opacity and n the number density of the particles responsible for the absorption. Assuming the dust is evenly distributed we can rewrite this as $\mathrm{d}\mathit{\tau }\mathrm{=}\frac{{\mathit{\sigma }}_{\mathrm{dust}}}{\mathit{V}}\mathrm{d}\mathit{r,}$(35)where V is the volume in which the dust is located, and σ_dust is the total combined cross-section for radiation of the dust grains ${\mathit{\sigma }}_{\mathrm{dust}}\mathrm{=}{\mathrm{\int }}_{{\mathit{R}}_{\min }}^{{\mathit{R}}_{\max }}\mathit{\pi }{\mathit{R}}^{\mathrm{2}}\mathit{f}\mathrm{\left(}\mathit{R}\mathrm{\right)}\mathrm{d}\mathit{R}\mathit{.}$(36)Assuming the dust traces the planetesimal population, i.e., destructive collisions happen throughout V, we can integrate Eq. (35) find ${\mathit{\tau }}^{\mathrm{\left(}\mathrm{M}\mathrm{2}\mathrm{\right)}}\mathrm{=}\mathrm{9.6}\mathrm{\times }{\mathrm{10}}^{-3}{\left(\frac{\mathit{R}}{\mathrm{km}}\right)}^{-0.5}{\left(\frac{\mathit{\rho }}{\mathrm{2.76}\hspace{0.17em}\mathrm{g}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}\left(\frac{\mathit{M}}{\mathrm{50} {\mathit{M}}_{\mathrm{\oplus }}}\right)\mathrm{·}$(37)The radial optical depth is indeed  < 1. Again, the derived optical depth depends a lot on the initial parameters. Since we have assumed the dust density is constant in the volume V, the radial optical depth is almost independent of the angle away from the midplane.

The optical depth could be increased by making the planetesimal population more massive, or changing their size. This will however also increase M_dust.

A better way to increase τ might be to confine the volume in which the dust is located. This can be done in two ways; either moving the entire planetesimal population closer to the star, or by realizing the planetesimals have pretty eccentric orbits and most destructive collisions might therefore happen near periapsis. Wyatt et al. (2010) showed, using a Monte-Carlo simulation, that a population of bodies with identical semi-major axes and eccentricities, but random mean longitudes, arguments of pericentre and longitudes of ascending nodes, 90% of destructive collisions happen within a radius r/a = (1 − e²)/(1 − 0.72e), which might be a hint in this direction, even though the bodies in the populations in Table 1 are distributed both in a and in e.

6. Discussion

Comparing the SEDs of several YSOs to radiative transfer models and near-IR interferometric visibilities appears to indicate the presence of an optically thin dust structure with a large scale height (Sect. 2). We have argued that an inwardly migrating planet could explain the observed structure. The scenario is that an inwardly migrating planet collects planetesimals in mean motion resonances, raises their eccentricities, and subsequently scatters them to orbits with greater inclinations (Fig. 3). Over time, the planetesimals collide and form small dust in a collisional cascade.

6.1. The role of eccentricity pumping

An important result of this study is that scattering planetesimals efficiently into highly-inclined orbits requires the relative velocity between the planet and the planetesimal to be high, of order of the local Keplerian speed. This is the reason why in our setup the sequence of eccentricity pumping in resonances followed by dynamical scattering is the most successful strategy.

The resonance capture probability of a migrating planet is a function of migration rate, planet mass, and distance to the star, and increases for more massive planets (Mustill & Wyatt 2011). This puts constraints on the planet mass and migration rate. The process of eccentricity-pumping is well-defined by Malhotra (1995); Wyatt (2003) and Eq. (14) can be used to estimate the distance over which a particle has to be dragged in a particular MMR for its orbit to reach a certain eccentricity. However, if a planetesimal gets stuck in a resonance that is too far from the planet, it is likely to be accreted by the central star after its eccentricity approaches unity. This fate could possibly be averted by introducing gas drag, or a second, stationary planet close to the central star. In this case, the migrating planet would excite and push the planetesimals inwards, while the stationary planet would cause the scattering events. The mass of this second planet should not be too high, as this would lead to the ejection of the excited planetesimals. Systems with Neptune-mass planets within an AU of their parent stars are not uncommon (e.g. Bouchy et al. 2009), and adding such a planet to the simulations presented in this work could well lead to even higher planetesimal inclinations.

Yu & Tremaine (2001) study resonant capture by a Jupiter-mass planet migrating inward using analytic arguments and three-body integrations. They find that in typical systems, the resonance capture and subsequent eccentricity-pumping will most likely result in accretion of the test particle by the central star. This fate can be avoided by having a close encounter with the migrating planet, which eventually leads to an escape from the system. In the study presented here, systems with very massive planets show very similar behaviour.

6.2. Migration

Our numerical simulations show that the basic scenario does work well for a certain range in planet mass and migration rate. The migration rates used are consistent with migration time scales for type I, type II, and migration through a planetesimal disk. However, the migration itself is not done self-consistently in our simulations. Because of this we have not taken into account the effect of the driver for the migration on the planetesimal population.

One possible driver is gas. A planet in a gaseous disk can exchange energy and angular momentum with the gas. If the planet is unable to significantly alter the disk structure it is referred to as type I. In this scenario the exchange of angular momentum can be written as the sum of the torques exerted by the Lindblad resonances of the planet on the disk (Goldreich & Tremaine 1979, 1980). The direction of the migration depends on the size of the interior and exterior torques and is found to be inward for an isothermal disk (Ward 1997). Paardekooper & Mellema (2006) simulated non-isothermal disks and showed the direction of migration depends on the ability of the disk to radiate away generated heat. Tanaka et al. (2002) calculated the sizes the torques for isothermal disks and derived migration time scales of the order of 1 Myr. The migration time scale is inversely proportional to planet mass, massive planets migrate faster (Armitage 2010).

A slightly more massive planet will be able to open up a gap in the disk. In this case the co-rotating torque disappears. The migration rate for type II migration is equal to the gas-inflow velocity and found to be independent of planet mass (Armitage 2010). Type II migration is generally faster than type I and at 5 AU a typical disk results in a migration time scale of 10⁵ yr.

Gas is not necessarily needed for migration though. A planet embedded in a planetesimal belt can exchange angular momentum with these bodies during scattering events. If the planet scatters bodies inward and outward at an equal rate, the net change in the angular momentum of the planet will be zero. If however, the planet usually scatters bodies to larger orbits, it will slowly loose angular momentum and migrate towards the star (Ida et al. 2000). For the change in angular momentum to have a significant effect on the planet’s orbital radius, the mass of the planet has to be lower or equal to the total mass of scattered bodies. Also, for the migration not to halt the planet has to encounter fresh planetesimals continuously. Ida et al. (2000) calculate the migration rate to be $\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{4}\frac{{\mathit{M}}_{\mathrm{disk}}}{{\mathit{M}}_{\mathrm{\odot }}}\frac{{\mathit{a}}_{\mathrm{pl}}}{{\mathit{P}}_{\mathrm{pl}}}\mathit{,}$(38)with the disk mass equal to ${\mathit{M}}_{\mathrm{disk}}\mathrm{=}\mathit{\pi }\mathrm{\Sigma }{\mathit{a}}_{\mathrm{pl}}^{\mathrm{2}}$. Kirsh et al. (2009) study the orbital evolution of an isolated planet in a planetesimal belt and confirm this result, adding that for typical planetesimal surface density profiles the direction of migration will generally be inward.

It is stressed that the planetesimals that are having their eccentricities raised in an interior resonance are not the ones responsible for the migration. Only when these bodies reach very eccentric orbits and suffer close encounters with the planet are they able to contribute to the planetary migration.

Aside from planet-disk interactions, the semi-major axis of a planet can also be changed significantly as a result of interactions with other planets. Raymond et al. (2011) numerically investigated systems with an initially unstable configuration of giant planets within a debris disk. It was found that the dynamical instability could trigger a collisional cascade.

6.3. The influence of gas in the disk

While planetesimal scattering alone can make a planet migrate, the preferred scenario may be to have the planet migrate in the phase when the disk is still gas-rich, while our simulations with forced migration have ignored the effects of gas on the planetesimals. If we attribute the planet’s radial motion to type I, type II, or planetesimal-driven migration in a gas-rich disk, we need to consider the effects of the gas during the migration and scattering phases. The presence of gas will cause aerodynamic gas drag. This drag will have two main implications; it will dampen the eccentricity and inclination of a planetesimal, and, it will have an effect of the probability of capturing a planetesimal in an MMR. Since the capturing of planetesimals followed by the eccentricity-pumping plays a vital role in this work, the effects of gas on these two mechanisms have to be addressed.

The aerodynamic gas drag time scale is given by Adachi et al. (1976)${\mathit{\tau }}_{\mathrm{aero}}\mathrm{=}\frac{\mathrm{8}\mathit{\rho R}}{\mathrm{3}{\mathit{C}}_{\mathrm{D}}{\mathit{\rho }}_{\mathrm{gas}}{\mathit{v}}_{\mathrm{K}}}\mathit{,}$(39)with ρ_gas and ρ the density of the gas and planetesimal material respectively. The drag coefficient C_D is of order unity, and a non-linear function of the planetesimal’s radius R and relative velocity to the gas. From the equation it is clear that gas drag will predominantly effect smaller bodies at relatively small distances to the star.

Capobianco et al. (2011) study planetesimal driven migration in a gas disk and find the eccentricity of a planetesimal is dampened by a factor of order itself on a time scale of ${\mathit{\tau }}_{\mathit{e}}\mathrm{=}\frac{\mathrm{2.3}\mathrm{\times }{\mathrm{10}}^{\mathrm{2}} \mathrm{yr}}{\mathit{f}\mathrm{\left(}\mathit{a,e}\mathrm{\right)}}\left(\frac{\mathrm{1.0}}{{\mathit{f}}_{\mathrm{g}}}\right)\left(\frac{\mathit{R}}{\mathrm{km}}\right){\left(\frac{\mathit{a}}{\mathrm{AU}}\right)}^{\mathrm{13}\mathit{/}\mathrm{4}}\mathit{,}$(40)where f_g = 1 corresponds to a MMSN-disk and f(a,e) is given by $\mathit{f}\mathrm{\left(}\mathit{a,e}\mathrm{\right)}\mathrm{=}\sqrt{{\mathit{\eta }}_{\mathrm{0}}^{\mathrm{2}}\left(\frac{\mathit{a}}{\mathrm{AU}}\right)\mathrm{+}\mathrm{7.5}\mathrm{\times }{\mathrm{10}}^{-5}{\left(\frac{{\mathit{M}}_{\mathrm{pl}}}{{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{\mathrm{2}\mathit{/}\mathrm{3}}{\left(\frac{\mathit{e}}{\mathit{\chi }}\right)}^{\mathrm{2}}}\mathit{,}$(41)with χ = R_H/a_pl and η₀ = 1.95 × 10^-3   AU ^− 1/2. The damping time scale is smallest for small bodies since they couple to the gas better. For a larger planetesimal of 10 km with a semi-major axis of 10 AU, the damping time scale is  ~10^6 − 7 yr and the system approaches the gas-free case (Capobianco et al. 2011). It appears that the eccentricity-damping plays a small role for planetesimals of sizes  >1−10 km. For smaller planetesimals or higher gas densities, gas drag will play an important role. From Eqs. (40) and (41) it is obvious gas drag will affect eccentric orbits more than circular ones. In such a regime, eccentricity-damping could therefore protect planetesimals against accretion by the central star. This would allow for them to be scattered by the migrating planet, even if they are initially captured in a resonance far from the planet.

The coupling of planetesimals to the gas has an effect on the probability of resonance capture as well. Capobianco et al. (2011) define a radius R_trap for which half of the planetesimals are captured. For a MMSN disk this radius corresponds to ${\mathit{R}}_{\mathrm{trap}}\mathrm{\simeq }\mathrm{9.2}\left(\frac{{\mathit{C}}_{\mathrm{D}}}{\mathrm{0.5}}\right){\left(\frac{\mathit{\rho }}{\mathrm{0.5}\hspace{0.17em}\mathrm{g}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{-3}}\right)}^{-1}{\left(\frac{{\mathit{M}}_{\mathrm{p}}}{{\mathit{M}}_{\mathrm{\oplus }}}\right)}^{-0.73}{\left(\frac{\mathit{a}}{\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{11}\mathit{/}\mathrm{4}}\hspace{0.17em}\mathrm{km}\mathit{,}$(42)with C_D the drag coefficient of order unity. For a Neptune-mass planet at 15 AU, R_drag ~ 0.01 km. Assuming the planetesimals are of km-size, the effect of gas on the capture probability can thus be neglected.

Fogg & Nelson (2005) have studied the effect of a giant planet migrating in on the process of terrestrial planet formation, using N-body simulations. The inward migration is implemented artificially with a prescription that is chosen to match type II migration. The disk is assumed to be gas-rich. Consequently, the gas drag causes radial drift of planetesimals, and a damping of eccentricities and inclinations. Fogg & Nelson find that the inward migration causes the planet to capture small bodies in its inner resonances, after which the orbits of these bodies are excited. The excitation of the orbits, followed by direct scattering by the planet, causes the formation of an exterior scattered disk of planetesimals, which looks very similar to the resulting populations in Fig. 6. Depending on the initial disk properties, the exterior disk is found to contain 30–70% of the inner disk material. The gas damping is found to play an important role in the inner regions of the disk. For increasingly mature disks, Fogg & Nelson find a wider spread in the orbital parameters of the scattered disk, caused by the increasingly smaller effect of gas damping.

In a somewhat similar study, with the focus of forming Earth-like planets in the habitable zone, is conducted by Mandell et al. (2007). They find similar results for an inwardly migrating Jupiter-mass planet as Fogg & Nelson (2005), including a scattered exterior disk. They also show that a second, stationary giant planet outside of the migrating one, will remove this scattered disk by ejecting the bodies from the system. A smaller or vanishing gas drag leads to more material ending up in the outer disk. The spread in orbital parameters of the bodies in the scattered disk also increases when less gas is present.

The extended scattered disks in these bodies show similar characteristics as the populations found in Table 1, when gas does not play a big role. However, these studies focus on giant planet migration, in relatively young disks, thus focussing on low-velocity collisions between the planetesimals that allow them to grow and form planet embryos and Earth-like planets. In the study presented in this paper, we are instead focussing on less massive planets, and interested in the orbital parameters of planetesimals in the scattered disk.

After the migration has taken place, the planetesimals collide on a certain time scale to form dust. The dust distribution can be calculated under some approximations and the total mass and radial optical depth are found to be consistent with the observed haloes for certain values of variable parameters such as the total mass in the disk, the initial size of the planetesimals, and the planetesimal material density. The collisional model used in Sect. 5, and the estimates of the optical depth and total mass are very simple and have their limitations. For instance, we assume the dust to trace the comet population, use simple geometrical estimates to define the volume the planetesimals occupy, and neglect the effect of any gas that might be present. Because of the nature of the discussed systems, including the spread in semi-major axes, eccentricities and inclinations, it is difficult to avoid these simplifications without having to make a detailed collisional and radiative transfer model, which is beyond the scope of this paper. If gas is present in the inner regions of the system, it will have a significant effect on the small dust grains since they couple to the gas on short time scales. Gas which is confined to the disk might rapidly remove small dust grains from highly-inclined orbits and reduce the scale height of the optically thin structure. Since we are dealing with pre-transitional disks, it is very well possible the inner regions of the system do not contain any gas at all.

7. Conclusions

The goal of this work was to see if an inwardly migrating planet could be the cause of the observed optically thin halo-like structures. The main findings are:

• An inwardly migrating planet is capable of capturing andtransporting material inward in mean motion resonances. Theprobability of capture depends on planet mass and migration rate.Larger planet masses and lower migration rates increase thecapture probability. The presence of gas on the captureprobability can be ignored for km-sized planetesimals.

• When a trapped planetesimal is pushed inwards, its eccentricity will rise rapidly and in agreement with Eq. (14). During this phase the inclination does not reach values greater than a couple of degrees. The growing eccentricity will eventually put the planetesimal on a planet-crossing orbit while increasing the relative velocity between the two.

• Simulations show that the inward migration of a planet can lead to a population of planetesimals on highly-inclined orbits. The average inclination of the remaining population increases for higher planet masses and lower migration rates. Resonance capture and subsequent eccentricity pumping plays an important role in this outcome.

• The high planetesimal inclinations in these systems result from direct scattering events, and reach significantly higher values if the planetesimal orbits are excited prior to this event.

• In single-planet systems where planetesimals get stuck in an MMR far from the migrating planet, the vast majority of planetesimals is lost to the central object. This fate could perhaps be averted by introducing a second, stationary planet close to the central star. This planet should not be very massive, as close encounters would then lead to ejection of the excited planetesimals.

• The resulting population of planetesimals on highly-inclined orbits have collisional time scales of the order of 10^5 − 6 yr depending on the properties of the individual planetesimals and their combined mass. Since the relative velocities between the planetesimals are high (~10   km   s^-1) the collisions will be destructive and will start a collisional cascade that will grind the bodies down to dust grains. The total mass and optical depth of this dust can be estimated and found to be of the same order as the total mass and optical depth observed in various systems.

Future studies should focus on self-consistent migration and a more careful description of the collisional cascade and optical depth determination to see if this mechanism can work and constrain the parameter-space further. With these future studies, it might be possible to use the optically thin structures as signposts for planet formation and migration. The presence of an optically thin halo-like structure could be induced from an SED, without having to resolve the system. The properties of the halo could then put constraints on the planet mass, migration rate, and distance it has travelled.

Acknowledgments

The authors would like to thank M. Wyatt for useful discussions, and the referee, P. J. Armitage, whose comments have led to a significant improvement of the manuscript.

References

All Tables

All Figures

Fig. 1 a) Scattering of a planetesimal in the stellar frame. The planetesimal comes in along the x-axis and the planet moves in the positive y-direction. The z-axis is perpendicular to the orbital plane. b) Unit vectors x and x′ (and y and y′) lie in the orbital plane, but are rotated with respect to each other by an angle α. c) The scattering in the frame of the planet. The planetesimal comes in along the x′-axis, with impact parameter b. It is scattered in the (y′ − z)-plane over an angle θ.

In the text

	Fig. 2 The inclination of a comet after an encounter at with a planet with , as a function of impact parameter and φ. The white regions mark the combinations of and φ where $\mathit{v̂}{}_{\mathrm{post}}\mathit{>}\sqrt{\mathrm{2}}$. Regions where b < bmin are also coloured white.
In the text

Fig. 3 Illustration of the mechanism to reach a planetesimal population with high orbits, as explored in this paper. a) both the planet orbit (dotted) and planetesimal orbits (solid) start out in the orbital plane. b) As the planet has migrated in a some distance, the planetesimals are stuck in an MMR. As they are dragged along their orbits become increasingly eccentric. c) At some point the planetesimal orbits become so eccentric, the bodies become planet-crossing. As the eccentricity rises even further, the relative velocity between the planet and the planetesimals grows. d) The high relative velocity compared to the local Keplerian velocity enables the planet to violently scatter the planetesimals throughout the system, to orbits with high inclinations.

In the text

	Fig. 4 A Neptune-mass planet is migrating in with ξ = 7.3 × 10-5   AU   yr-1. The asteroid survives the passing of the 1:2 resonance but is captured in the 3:4 resonance and dragged inward.
In the text

Fig. 5 Evolution of the eccentricity of a massless body stuck in the 1:2 internal mean motion resonance of a 10   MNep planet migrating inward at a rate of 7.3 × 10-5   AU   yr-1. The timestep for both the SyMBA and Bulirsch-Stoer integration was set at 0.5 days. The dotted line shows the expected increase in eccentricity as predicted by Eq. (14). The test particle is captured at t ≈ 0.13 Myr and eventually lost to the central star after about 0.32 Myr.

In the text

	Fig. 6 Evolution of the eccentricities (left) and inclinations (right) in a planetesimal swarm as a planet (large circle) migrates inwards from 30 AU at a rate of 7.3 × 10-5   AU   yr-1. The planet mass used are 0.1   MNep for a), and 2   MNep for b). Planetesimals (small circles) above the solid lines are on a planet crossing orbit, and the dotted lines show the location of the 1:2 and 2:3 MMR.
In the text