© ESO, 2014

1. Introduction

Circumstellar dust in debris disks reveals the location and dynamical state of larger bodies and thus sheds light on the architecture of planetary systems in the aftermath of planet formation (see Wyatt 2008 for a review). The dust can be studied by observing its infrared and (sub-)millimeter emission, as well as the stellar radiation it scatters, and is usually found at large distances from the star (tens of AUs, Carpenter et al. 2009). Recently, interferometric observations have found excess near-infrared (NIR) emission emanating from the innermost parts of several debris disks, which has been interpreted as thermal emission from hot (>1000 K) dust (Ciardi et al. 2001; Absil et al. 2006, 2008, 2009; di Folco et al. 2007; Akeson et al. 2009; Defrère et al. 2011, 2012, see Table 1 for an overview). This material is known as hot exozodiacal dust. Its origin, hence what it can tell us about planet formation, is still unclear. In this work, we investigate one possible scenario for explaining hot exozodiacal dust.

Dust grains in debris disks have relatively short lifetimes because of their destruction by collisions and removal by radiation forces. To detect these grains around mature stars therefore implies the existence of a mechanism that continuously replenishes them. Cold dust populations at large distances from the star can be maintained by a collisional cascade grinding down much larger bodies that act as a reservoir of mass (Backman & Paresce 1993). Closer to the star, however, the pace at which material is processed by collisions is much higher, so that the lifetime of a debris belt in collisional equilibrium is much shorter there (Dominik & Decin 2003; Wyatt et al. 2007). For this reason, hot exozodiacal dust cannot be explained by in situ planetesimal belts (Wyatt et al. 2007; Lebreton et al. 2013)¹, and a different mechanism is needed to replenish it, the lifetime of the dust needs to be extended by some process, or both.

Many of the systems that exhibit NIR excess also feature a debris belt at a large distance from the star (see Table 1). Inward transport of material from an outer belt may therefore be a natural explanation for the existence of hot exozodiacal dust. A possible transportation mechanism is Poynting–Robertson (P–R) drag (see, e.g., Burns et al. 1979). Because P–R drag acts on a timescale that is much longer than that of collisions, it is sometimes not favored as a possible mechanism for maintaining exozodiacal dust (e.g., Absil et al. 2006). However, as long as there are no mechanisms that prevent inward migration, a small amount of dust is always transported to the innermost part of the disk (Wyatt 2005), where it produces a NIR signal.

Morphological models of exozodiacal dust disks, constrained by the NIR observations, indicate that the hot dust is concentrated in a sharply peaked ring, whose inner boundary is determined by dust sublimation (Defrère et al. 2011; Mennesson et al. 2013; Lebreton et al. 2013). The process of dust sublimation may therefore play an important role in shaping exozodiacal clouds. Kobayashi et al. (2009) find that the interplay between P–R drag and dust sublimation can lead to a local enhancement of dust in the sublimation zone, leading to radial distributions of dust reminiscent of those found by the morphological models. However, they only investigate this pile-up effect for drag-dominated systems, where collisions are unimportant, and it is unclear what happens to the phenomenon if collisions are taken into account.

In this work, we examine whether it is possible to maintain a pile-up of dust in the sublimation zone of a collisionally active debris disk and whether such a pile-up could explain the exozodiacal NIR emission observed very close to some stars. To do this, we compute the steady-state distribution of dust in the inner parts of debris disks by solving the continuity equation, considering collisions, P–R drag, and sublimation. First, we find an analytical solution, using a number of simplifying assumptions (Sect. 2). Subsequently, we solve the continuity equation numerically using a debris disk model that for the first time treats the complex interaction of collisions, P–R drag, and sublimation in a self-consistent way (Sect. 3). From the obtained steady-state dust distributions, we compute emission spectra to compare with observational data (Sect. 4). We discuss our findings in Sect. 5, and give conclusions in Sect. 6. Details of the numerical techniques employed by the debris disk model are given in Appendix A, model verification tests are described in Appendix B, and the post-processing of model output into useful physical quantities in described in Appendix C.

2. Analytical constraints

In this section we analytically investigate the distribution of dust in the inner regions of debris disks. We focus on the radial distribution of material, by assuming (1) that the disk is axisymmetric and (2) that all particles have the same size. Throughout this work, radial distributions are expressed in terms of vertical geometrical optical depth, which is defined as the surface density of cross-section². Under the two assumptions listed above, it is given by ${\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{\sigma n}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{2}\mathit{\pi r}}\mathrm{·}$(1)Here, r is the radial distance from the central star, σ the cross-section of a particle, and n(r) the one-dimensional number density (i.e., the particle number density integrated over disk height and azimuth).

We consider three processes affecting the evolution of dust particles in debris disks: collisions, P–R drag, and sublimation. The strategy for our analytical estimates is as follows. First, we review the balance between P–R drag and collisions (without considering sublimation) and calculate the inward flux of material due to these two effects (Sect. 2.1). Subsequently, we consider the interplay of P–R drag and sublimation (ignoring collisions), which can lead to the pile-up of dust in the sublimation zone (Sect. 2.2). Finally, we investigate whether collisions can be neglected in the innermost parts of a debris disk and estimate the radial distribution of dust in a disk where all three processes are in operation (Sect. 2.3). At the end of the section, we briefly summarize our findings (Sect. 2.4).

2.1. Poynting–Robertson drag and collisions

The balance between P–R drag and collisions was studied analytically by Wyatt (2005). Because of its importance to the present study, we summarize the main arguments of this work here. The model assumes that (1) there is a source of dust at distance r₀ with a geometrical optical depth of τ_geo(r₀); (2) dust particles follow circular orbits; (3) collisions are always destructive; and (4) all dust grains have the same size. These assumptions lead to simple expressions for the timescales on which P–R drag and collisions typically act.

The P–R drag timescale t_PR is defined as the time it takes for a particle on a circular orbit to spiral from a given distance r to the central star. It is given by (e.g., Burns et al. 1979) ${\mathit{t}}_{\mathrm{PR}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathit{c}{\mathit{r}}^{\mathrm{2}}}{\mathrm{4}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathit{\beta }}\mathit{,}$(2)where c is the speed of light, G the gravitational constant, M_⋆ the stellar mass, and β the ratio of the norms of the direct radiation pressure force and the gravitational force on a particle (β = | F_rp/F_g |)³.

The collisional timescale t_coll indicates the average time between two collisions for a given particle. Wyatt (2005) finds ${\mathit{t}}_{\mathrm{coll}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{orb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{4}\mathit{\pi }{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathit{,}$(3)where t_orb is the orbital period of a circular orbit, given by ${\mathit{t}}_{\mathrm{orb}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }\sqrt{{\mathit{r}}^{\mathrm{3}}\mathit{/}\mathrm{\left(}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathrm{\right)}}$. This equation is valid for particles that are on circular orbits and whose most important collisional partners have similar sizes⁴.

2.1.1. The radial distribution due to P–R drag and collisions

Under the assumptions listed above, there is an analytical solution to the continuity equation, balancing the migration of particles due to P–R drag with the destruction of dust by collisions. Wyatt (2005) finds that the steady-state solution is $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \frac{{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{1}\mathrm{+}\mathrm{4}{\mathit{\eta }}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{-}\sqrt{\mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{0}}}\mathrm{)}}\mathit{,} \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{,}\\ {\mathit{\eta }}_{\mathrm{0}}& \mathrm{=}& \frac{\mathit{c}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}{\mathit{r}}_{\mathrm{0}}\mathrm{\right)}}{\mathrm{2}\mathit{\beta }}\sqrt{\frac{{\mathit{r}}_{\mathrm{0}}}{\mathit{G}{\mathit{M}}_{\mathit{\star }}}}\mathrm{·}\end{array}$The parameter η₀ characterizes the density of the parent belt. It is defined such that for η₀ = 1 the collisional and P–R drag timescales are equal at r₀. Disks with η₀> 1 are collision-dominated at r₀, while disks with η₀< 1 are drag-dominated at r₀. Most debris disks with observed outer belts have η₀> 10 (Wyatt 2005).

The balance between P–R drag and collisions is self-limiting: a denser parent belt produces more dust drifting inwards, but this dust also suffers more mutual collisions that eliminate grains on their way in. This results in a maximum geometrical optical depth profile for η₀ ≫ 1 (i.e., a very dense parent belt) of $\begin{array}{c}\max \end{array}\mathrm{\left[}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\right(}\mathit{r}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\frac{\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}}\hspace{0.17em}\mathit{\beta }}{\mathrm{2}\mathit{c}\mathrm{(}\sqrt{{\mathit{r}}_{\mathrm{0}}}\mathrm{-}{\sqrt{\mathit{r}}}^{\mathrm{)}}}\mathit{,} \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{.}$(5)Figure 1 shows examples of this radial profile for different parent belt locations and host star masses, all for dust grains with β = 0.5. The value for β was set to the blowout limit: particles with β> 0.5 leave the system on hyperbolic paths after they are released from a large parent body on a circular orbit. Therefore, P–R drag is the most efficient for β = 0.5, and for a given system, this value corresponds to the maximum τ_geo profile⁵.

Fig. 1 Maximum geometrical optical depth as a function of distance to the star (Eq. (5)). The profiles are derived from the analytical model of Wyatt (2005), in which dust is produced by a source at radius r0 and subsequently migrates inward due to P–R drag, while suffering destruction from mutual collisions. The solid lines correspond to solar-mass stars, the dashed lines to 2 M⊙ stars. Profiles for different parent belt locations are shown in different colors. All profiles assume dust grains with β = 0.5.

2.1.2. The inward flux of material

Since dust can pile up close to the star due to sublimation (to be discussed in Sect. 2.2.2), τ_geo may exceed the upper limit given by Eq. (5) in the sublimation zone. The material that piles up, however, is supplied from farther out by P–R drag. To investigate the properties of the pile-up, it is therefore useful to assess the inward flux of material.

In the case of a uniform grain size, the inward particle flux due to P–R drag (i.e., the number of particles passing through a ring at radius r per unit of time) can be expressed as ${\mathit{\varphi }}_{\mathrm{PR}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathit{n}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{ṙ}\mathrm{PR}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(6)where ṙ_PR is the P–R drag velocity, counted positively towards r> 0. If the particle orbits are circular, radial migration due to P–R drag is described by (e.g., Burns et al. 1979) $\mathit{ṙ}\mathrm{PR}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathrm{2}\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathit{\beta }}{\mathit{cr}}\mathrm{·}$(7)The grain parameter β can be expressed as (e.g., Burns et al. 1979) $\mathit{\beta }\mathrm{=}\frac{{\mathit{L}}_{\mathit{\star }}}{\mathrm{4}\mathit{\pi cG}{\mathit{M}}_{\mathit{\star }}}\frac{{\mathit{Q}}_{\mathrm{pr}}\mathit{\sigma }}{\mathit{m}}\mathrm{·}$(8)Here, L_⋆ is the stellar luminosity, m the mass of an individual dust grain, and Q_pr the radiation pressure efficiency averaged over the stellar spectrum. Combining Eqs. (1), (6)–(8) and multiplying by particle mass m gives the inward (collision-limited) P–R drag mass flux (cf. Rafikov 2011) $\mathit{Ṁ}\mathrm{PR}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{L}}_{\mathit{\star }}}{{\mathit{c}}^{\mathrm{2}}}{\mathit{Q}}_{\mathrm{pr}}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}$(9)The maximum mass flux can now be found by substituting Eq. (5) for τ_geo(r) in Eq. (9), which yields $\begin{array}{c}\max \end{array}\left[\mathit{Ṁ}\mathrm{PR}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]\mathrm{=}\frac{\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}}\hspace{0.17em}{\mathit{L}}_{\mathit{\star }}\mathit{\beta }{\mathit{Q}}_{\mathrm{pr}}}{\mathrm{2}{\mathit{c}}^{\mathrm{3}}\mathrm{(}\sqrt{{\mathit{r}}_{\mathrm{0}}}\mathrm{-}{\sqrt{\mathit{r}}}^{\mathrm{)}}}\mathit{,} \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{.}$(10)Figure 1 shows that max [ τ_geo(r) ] levels off in the innermost part of the system (r ≪ r₀), at the radii where exozodiacal dust is seen. The mass flux corresponding to this plateau is $\begin{array}{ccc}\begin{array}{c}\max \end{array}\mathrm{\left[}\mathit{Ṁ}\mathrm{PR}\mathrm{\right(}\mathit{r}\mathrm{=}\mathrm{0}\mathrm{\left)}\mathrm{\right]}& & \mathrm{\approx }\mathrm{5.6}\mathrm{\times }{\mathrm{10}}^{-13}\hspace{0.17em}(\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathrm{1}\mathrm{M}}_{\mathrm{\odot }}}{)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}(\frac{{\mathit{L}}_{\mathit{\star }}}{{\mathrm{1}\mathrm{L}}_{\mathrm{\odot }}})\\ & & \mathrm{\times }(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathrm{1}\mathrm{AU}}{)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}(\frac{{\mathit{Q}}_{\mathrm{pr}}}{\mathrm{1}})\hspace{0.17em}(\frac{\mathit{\beta }}{\mathrm{0.5}}){\mathit{M}}_{\mathrm{\oplus }}{\mathrm{yr}}^{-1}\mathit{.}\end{array}$(11)The maximum mass flux only depends on grain properties through Q_pr and β⁶. The radiation pressure efficiency Q_pr must obey 0 ≤ Q_pr ≤ 2, and particles with β> 1 are always unbound. Therefore, Eq. (11) with Q_pr = 2 and β = 1 gives a solid upper limit on the inward mass flux due to P–R drag, unless one of the model assumptions does not hold (e.g., collisions are non-destructive).

2.2. P–R drag and sublimation

In the preceding, we found that P–R drag supplies a small but non-zero amount of dust to the innermost parts of a debris disk. As this material approaches the central star, it is heated by stellar radiation. Eventually, the dust grains become so hot that they start sublimating. We now review the evolution of these particles considering P–R drag and sublimation, but ignoring collisions.

2.2.1. Dust sublimation formalism

For a spherical dust grain in a gas-free environment, the rate at which the grain radius s changes is given by (e.g., Kobayashi et al. 2008) $\frac{\mathrm{d}\mathit{s}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\frac{{\mathit{P}}_{\mathrm{v}}\mathrm{\left(}\mathit{T}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{d}}}\sqrt{\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{u}}}{\mathrm{2}\mathit{\pi }{\mathit{k}}_{\mathrm{B}}\mathit{T}}}\mathrm{·}$(12)Here, P_v is the phase-equilibrium vapor pressure, ρ_d the bulk density of the dust, μ the molecular weight of dust molecules, m_u the atomic mass unit, k_B the Boltzmann constant, and T the temperature of the dust. This theoretical sublimation rate is sometimes lowered to comply with experimental results, parameterized in a sticking efficiency or accommodation coefficient. Here, we ignore this weak effect, by assuming a sticking efficiency of unity.

The temperature dependence of P_v is given by (Kobayashi et al. 2009) ${\mathit{P}}_{\mathrm{v}}\mathrm{\left(}\mathit{T}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{0}}\exp \left(\mathrm{-}\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{u}}\mathit{H}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}\right)\mathit{,}$(13)where P₀ is a normalization constant and H the latent heat of sublimation. By assuming P₀ is constant, we neglect a small temperature dependence beyond the exponential.

Sublimation parameters depend on material and can be determined by laboratory measurements. The material we consider in this study is carbonaceous dust. This choice is motivated by the proximity of hot exozodiacal dust to its host star, which suggests a very refractory material, like carbon (Mennesson et al. 2013; Lebreton et al. 2013). Specifically, we use the sublimation parameters of graphite, for which many laboratory measurements are available. For the molecular weight, we use μ = 36.03, reflecting that graphite sublimation typically releases clusters of three carbon atoms at the temperatures and pressures relevant to this work (Abrahamson 1974). The parameters for C₃ sublimation are P₀ = 2.95 × 10¹⁴ dyn cm^-2 and H = 2.15 × 10¹¹ erg g^-1 (Zavitsanos & Carlson 1973). For the bulk density of the material, we use ρ_d = 1.8 g cm^-3.

For our analytical estimates, we approximate the grain temperature T by its black-body temperature ${\mathit{T}}_{\mathrm{bb}}\mathrm{=}{\left(\frac{{\mathit{L}}_{\mathit{\star }}}{\mathrm{16}\mathit{\pi }{\mathit{\sigma }}_{\mathrm{SB}}{\mathit{r}}^{\mathrm{2}}}\right)}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}$(14)where σ_SB is the Stefan–Boltzmann constant. In reality, the grain temperature is a function of size. Temperatures are generally higher than T_bb for particles that are smaller than the typical wavelength of the stellar radiation, because these grains do not cool efficiently. For simplicity, we ignore this effect here and investigate it further in our numerical calculations, which use realistic grain temperatures (see Sect. 3.2.2).

In the black-body approximation, the sublimation rate ṡ becomes independent of grain size. This leads to a simple expression for the sublimation timescale t_subl, defined as the time it takes for a spherical dust grain to disappear, which is ${\mathit{t}}_{\mathrm{subl}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{-}\frac{\mathit{s}}{\mathit{ṡ}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{·}$(15)This estimate assumes that the grain temperature remains constant (at T_bb) throughout the sublimation process. As a sublimating particle becomes smaller, the black-body approximation is bound to become inaccurate. For sufficiently large particles, however, the true sublimation time is dominated by the black-body regime, and Eq. (15) provides a good approximation.

2.2.2. The pile-up of dust in the sublimation zone

As dust grains become smaller due to sublimation, their β-ratio changes. As a result of this, the interplay between P–R drag and sublimation can lead to a pile-up of dust in the sublimation zone. This phenomenon was studied in detail by Kobayashi et al. (2008, 2009, 2011). Here, we give a brief explanation of the pile-up mechanism.

When dust grains migrate inward owing to P–R drag, their temperature gradually increases. At some point, the grains are heated to the point where sublimation becomes substantial, and their sizes start decreasing significantly. For particles larger than the peak wavelength of the stellar spectrum λ_⋆⁷, the radiation pressure efficiency is roughly constant at Q_pr ≈ 1, and therefore the approximation β ∝ s^-1 holds (see Eq. (8)). As the particles lose mass, radiation pressure therefore becomes more important (relative to stellar gravity), which has the effect of increasing the semi-major axes and eccentricities of their orbits, compensating for the decrease terms due to P–R drag. This effectively slows down the inward migration of the dust grains, leading to an accumulation of dust in the sublimation zone. Eventually, the dust grains either sublimate completely or their β ratios increase to the point where they become unbound and are blown out of the system.

Accumulation of dust occurs when the decrease of semi-major axis due to P–R drag is compensated by the increase of semi-major axis due to sublimation. This happens approximately at the radial distance where the timescale of P–R drag equals that of sublimation (Kobayashi et al. 2008)⁸, We denote this distance with r_pile, and determine its value by solving t_PR(r_pile) = t_subl(r_pile). Assuming that the dust in the pile-up has the black-body temperature T_pile = T_bb(r_pile), which holds for s ≳ λ_⋆, this equation can be written as $\frac{\mathrm{12}{\mathit{\sigma }}_{\mathrm{SB}}}{{\mathit{c}}^{\mathrm{2}}}\frac{{\mathit{Q}}_{\mathrm{pr}}}{{\mathit{P}}_{\mathrm{0}}}{\mathit{T}}_{\mathrm{pile}}^{\mathrm{4}}\mathrm{=}\exp \left(\mathrm{-}\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{u}}\mathit{H}}{{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{pile}}}\right)\sqrt{\frac{\mathit{\mu }{\mathit{m}}_{\mathrm{u}}}{\mathrm{2}\mathit{\pi }{\mathit{k}}_{\mathrm{B}}{\mathit{T}}_{\mathrm{pile}}}}\mathit{,}$(16)which shows that T_pile is independent of stellar parameters, the bulk density of the dust, and grain size (Q_pr is nearly constant for s ≳ λ_⋆), and only depends on material properties (cf. Kobayashi et al. 2008, 2009, 2011). Using Q_pr = 1 and the sublimation parameters of graphite given in Sect. 2.2.1, we numerically solve Eq. (16) to find T_pile ≈ 2020 K, hence ${\mathit{r}}_{\mathrm{pile}}\mathrm{\approx }\mathrm{0.019}\hspace{0.17em}(\frac{{\mathit{L}}_{\mathit{\star }}}{\mathrm{1}{\mathit{L}}_{\mathrm{\odot }}}{)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}(\frac{{\mathit{T}}_{\mathrm{pile}}}{\mathrm{2020}\mathrm{K}}{)}^{-2}\mathrm{AU}\mathit{.}$(17)The pile-up distance is independent of grain size because larger particles take longer to sublimate, but also migrate more slowly because of P–R drag. This holds as long as the black-body approximation is valid, so the grain temperature in the pile-up is independent of grain size.

2.2.3. Conditions for dust pile-up

Kobayashi et al. (2011) list two conditions for the accumulation of dust to be substantial: (1) sufficiently high values of β need to be reached as dust grains sublimate; and (2) the orbital eccentricities of the dust grains need to be low enough when they enter the sublimation zone.

Considering an inward stream of dust grains with a range of sizes, the particles contributing the most to the pile-up are those with the highest β. These particles have the strongest P–R drag drift rates, so per unit of time more of them arrive in the sublimation zone, where their inward drift is canceled as a result of sublimation. In addition, for the typical size distribution resulting from a collisional cascade, the total cross-section is dominated by the smallest particles. Since dust grains with β> 0.5 are typically blown out as soon as they are created, particles that are barely bound (β ≈ 0.5) before the onset of substantial sublimation are the most important for the pile-up. A requirement for dust migration to slow down is that β increases as the grain size decreases. For a given system, however, β reaches a maximum value β_max at s ~ λ_⋆, because smaller particles have lower Q_pr. Considering that relatively high values of β are needed for an efficient pile-up, low luminosity stars do not have significant pile-ups. Kobayashi et al. (2009, 2011) set the limit at β_max> 0.5, and derive a lower limit on the stellar luminosity for significant pile-up, assuming that Q_pr ≈ 1 holds for s ≳ λ_⋆. This limit is ${\mathit{L}}_{\mathit{\star }}\gtrsim \mathrm{0.5}\hspace{0.17em}\left(\frac{{\mathit{M}}_{\mathit{\star }}}{\mathrm{1}{\mathit{M}}_{\mathrm{\odot }}}\right)\hspace{0.17em}{\left(\frac{{\mathit{T}}_{\mathit{\star }}}{\mathrm{5}\mathrm{\times }{\mathrm{10}}^{\mathrm{3}}\mathrm{K}}\right)}^{-1}\hspace{0.17em}\left(\frac{{\mathit{\rho }}_{\mathrm{d}}}{\mathrm{1.0}\mathrm{g}\mathrm{c}{\mathrm{m}}^{-3}}\right){\mathit{L}}_{\mathrm{\odot }}\mathit{,}$(18)where T_⋆ is the effective temperature of the central star.

The pile-up of dust in the sublimation zone is highly dependent on the eccentricity of the dust as it enters the sublimation zone (Kobayashi et al. 2008). For the pile-up mechanism to produce a significant enhancement, the eccentricity of the dust particles in the sublimation zone must be very low (e ≲ 10^-2; Kobayashi et al. 2008, 2011). Particles with higher orbital eccentricities do not spend enough time in the sublimation zone before they are blown out (or sublimate completely) to contribute significantly to the dust enhancement.

When dust particles are created in collisions in the parent belt, they are put on eccentric orbits with their periastron in the parent belt. Particles released from circular orbits will acquire orbital eccentricities of $\mathit{e}\mathrm{=}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\mathrm{·}$(19)The particles that are the most important for the pile-up are the ones with β ≈ 0.5. These barely bound particles initially follow very elliptic orbits, with eccentricities close to unity.

The initial eccentricity of the β ≈ 0.5 particles is far too high for any significant pile-up to occur. However, as the dust grains migrate inward owing to P–R drag, their orbits are circularized. The eccentricity evolution of dust particles experiencing P–R drag is coupled to their orbital size evolution according to (Wyatt & Whipple 1950) $\frac{{\mathit{a}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{)}}{{\mathit{a}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}_{\mathrm{0}}^{\mathrm{2}}\mathrm{)}}\mathrm{=}{\left(\frac{{\mathit{e}}_{\mathrm{1}}}{{\mathit{e}}_{\mathrm{0}}}\right)}^{\mathrm{4}\mathit{/}\mathrm{5}}\mathit{,}$(20)where (a₀,e₀) are the initial semi-major axis and eccentricity, which evolve into (a₁,e₁). This coupling can be used to place a lower limit on the distance of the source region, if significant pile-up is to occur, using the maximum allowed eccentricity in the sublimation zone for efficient pile-up. Written in terms of periastron distance q = a(1 − e), Eq. (20) becomes $\frac{{\mathit{q}}_{\mathrm{1}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{1}}\mathrm{)}}{{\mathit{q}}_{\mathrm{0}}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{e}}_{\mathrm{0}}\mathrm{)}}\mathrm{=}{\left(\frac{{\mathit{e}}_{\mathrm{1}}}{{\mathit{e}}_{\mathrm{0}}}\right)}^{\mathrm{4}\mathit{/}\mathrm{5}}\mathrm{·}$(21)Substituting r₀ and r_pile for q₀ and q₁, respectively, and using e₀ ≈ 1 and e₁ ≲ 10^-2, yields a lower limit on the source radius ${\mathit{r}}_{\mathrm{0}}\gtrsim \mathrm{20}{\mathit{r}}_{\mathrm{pile}}\mathit{,}$(22)for a significant enhancement in the dust density to be possible. Other mechanisms may help in the circularization of the orbits of small particles, decreasing the limit on r₀. An example is the drag force from small amounts of gas that are present in the disk (Takeuchi & Artymowicz 2001).

2.3. Pile-up in a collisionally active disk

So far, we have looked separately at the balance between P–R drag and collisions (without considering sublimation) and at the balance between P–R drag and sublimation (without considering collisions). We now investigate under what conditions this pairwise approach is justified, and combine the previous findings to estimate the distribution of dust in a debris disk in which all three processes are operational.

Fig. 2 Characteristic timescale as function of radial distance for sublimation (Eq. (15)), P–R drag (Eq. (2)), and mutual collisions (minimum timescale, Eq. (23)), for β = 0.5 particles in a debris disk around a solar-mass star with a very dense (η0 ≫ 1) parent belt located at 30 AU. The gray vertical lines indicate the radial distances used by the analytical model: rpile is the radius for which the sublimation timescale equals the P–R drag timescale, rcrit is the radius for which the P–R drag timescale equals the collisional timescale, and r0 is the location of the parent belt. For the sublimation timescale, we assume that the dust grains are solid spheres of graphite (see Sect. 2.2.1 for the values of the sublimation parameters) with a radius of s = 0.64 μm (the size corresponding to β = 0.5).

2.3.1. Do collisions interfere with dust pile-up?

Since collisions might interfere with the process of dust pile-up, the results of Sect. 2.2 are only valid if collisions do not play an important role at distances where sublimation becomes significant. We now investigate whether collisions can indeed be neglected in the inner regions of debris disks by comparing the characteristic timescales of the three processes as a function of radius.

In the case of a very dense parent belt (η₀ ≫ 1), τ_geo(r) is given by Eq. (5). Putting this in Eq. (3) results in the minimum collisional timescale $\begin{array}{c}\min \end{array}\left[{\mathit{t}}_{\mathrm{coll}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\right]\mathrm{=}\frac{\mathit{c}{\mathit{r}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathit{\beta G}{\mathit{M}}_{\mathit{\star }}}\mathrm{(}\sqrt{{\mathit{r}}_{\mathrm{0}}}\mathrm{-}{\sqrt{\mathit{r}}}^{\mathrm{)}}\mathit{,} \mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{.}$(23)In Fig. 2 this timescale is compared with the sublimation and P–R drag timescales for a debris disk around a solar-mass star with a dense parent belt at 30 AU, consisting of barely bound (β = 0.5) dust grains and using the sublimation parameters of graphite. In this example system, the collisional timescale is longer than the other timescales at r_pile, implying that P–R drag dominates collisions in the sublimation zone, and collisions do not interfere with dust pile-up. We now check whether this a general result or under which conditions it is the case.

We let r_crit be the radial distance at which collisional and P–R drag timescales are equal. Solving t_PR(r_crit) = t_coll(r_crit) for r_crit, with τ_geo(r) given by Eq. (4), yields $\frac{{\mathit{r}}_{\mathrm{crit}}}{{\mathit{r}}_{\mathrm{0}}}\mathrm{=}{\left(\frac{\mathrm{4}{\mathit{\eta }}_{\mathrm{0}}\mathrm{+}\mathrm{1}}{\mathrm{5}{\mathit{\eta }}_{\mathrm{0}}}\right)}^{\mathrm{2}}\mathrm{·}$(24)For η₀ = 1 we recover r_crit = r₀, as expected. From the limiting case η₀ ≫ 1 (i.e., a very dense parent belt), we find ${\mathit{r}}_{\mathrm{crit}}\mathit{>}\mathrm{0.64}{\mathit{r}}_{\mathrm{0}}\mathit{.}$(25)This shows that far enough inward from the parent belt, P–R drag always dominates collisions. Furthermore, if Eq. (22) holds we find r_crit ≳ 12.8r_pile, which means that in systems where the parent belt is distant enough for significant dust pile-up to occur, P–R drag dominates collisions in the sublimation zone, and collisions are so infrequent there that they do not interfere with the pile-up process.

The pile-up of dust means that τ_geo increases around r_pile, locally decreasing the collisional timescale. However, a significant pile-up requires Eq. (22) to hold, in which case the ratio between the minimum collisional timescale and the P–R drag timescale is found to satisfy $\min {}^{\mathrm{[}}\mathit{t}_{\mathrm{coll}}\mathrm{(}{\mathit{r}}_{\mathrm{pile}}{\mathrm{\left)}}^{\mathrm{\right]}}\mathit{/}{\mathit{t}}_{\mathrm{PR}}\mathrm{\left(}{\mathit{r}}_{\mathrm{pile}}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathrm{(}\sqrt{{\mathit{r}}_{\mathrm{0}}\mathit{/}{\mathit{r}}_{\mathrm{pile}}}\mathrm{-}\mathrm{1}\mathrm{)}\gtrsim \mathrm{13.9}$. To overcome this difference, the pile-up would have to raise τ_geo by the same factor. Since the τ_geo enhancement factor is never found to be greater than about 10 (Kobayashi et al. 2009), we expect the disk to remain drag (or sublimation) dominated inside r_crit.

2.3.2. Estimating the pile-up magnitude

The efficiency of dust pile-up was studied in detail by Kobayashi et al. (2009, 2011), who give formulae for the resulting enhancement in particle number density and geometrical optical depth. Here, we present a simple order-of-magnitude estimate of the maximum amount of material in the pile-up, which can be used to assess whether pile-ups can explain the observed NIR excess emission.

Dust grains reside in the pile-up for roughly one sublimation timescale, after which they are either completely sublimated or their size is reduced so much that they are blown out of the system. Since the pile-up occurs roughly at r_pile, defined such that t_subl(r_pile) = t_PR(r_pile), the dust stays in the pile-up for about one P–R drag timescale. Given this, the total number of particles in the pile-up is ${\mathit{N}}_{\mathrm{pile}}\mathrm{=}{\mathit{\varphi }}_{\mathrm{PR}}\mathrm{\left(}{\mathit{r}}_{\mathrm{pile}}\mathrm{\right)}{\mathit{t}}_{\mathrm{PR}}\mathrm{\left(}{\mathit{r}}_{\mathrm{pile}}\mathrm{\right)}\mathit{.}$(26)To describe the radial profile in terms of geometrical optical depth, it is necessary to specify the radial width of the pile-up Δr_pile. In reality, Δr_pile depends on the orbital eccentricities of the particles in the pile-up and on differences in pile-up distance for particles of different sizes that contribute. Since this is beyond the scope of this work, we keep the relative pile-up width Δr_pile/r_pile as a free parameter.

Combining Eqs. (1), (6), and (26), we find that the geometrical optical depth due to the material in the pile-up is given by ${\mathit{\tau }}_{\mathrm{geo}\mathit{,}\hspace{0.17em}\mathrm{pile}}\mathrm{=}\frac{\mathit{\sigma }{\mathit{N}}_{\mathrm{pile}}}{\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{pile}}\mathrm{\Delta }{\mathit{r}}_{\mathrm{pile}}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{pile}}}{\mathrm{2}\mathrm{\Delta }{\mathit{r}}_{\mathrm{pile}}}{\mathit{\tau }}_{\mathrm{geo}\mathit{,}\hspace{0.17em}\mathrm{base}}\mathrm{\left(}{\mathit{r}}_{\mathrm{pile}}\mathrm{\right)}\mathit{,}$(27)where τ_{geo, base}(r) denotes the base level of dust due to P–R drag and collisions given by Eq. (4). Kobayashi et al. (2009) find that sublimating dust particles slowly move outward. Therefore, we assume that the pile-up extends from r_pile outward, overlapping with the inward migrating material. The complete geometrical optical depth profile, which includes the effects of collisions, P–R drag, and sublimation, can then be formally described by (cf. Kobayashi et al. 2011) ${\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{for}\mathit{r}\mathit{<}{\mathit{r}}_{\mathrm{pile}}\\ & \mathrm{for}{\mathit{r}}_{\mathrm{pile}}\mathrm{\le }\mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{pile}}\mathrm{+}\mathrm{\Delta }{\mathit{r}}_{\mathrm{pile}}\\ & \mathrm{for}{\mathit{r}}_{\mathrm{pile}}\mathrm{+}\mathrm{\Delta }{\mathit{r}}_{\mathrm{pile}}\mathit{<}\mathit{r}\mathrm{\le }{\mathit{r}}_{\mathrm{0}}\mathit{.}\end{array}$(28)Using Eq. (5) for τ_{geo, base}(r) gives the maximum profile. Figure 3 shows this maximum profile for different values of Δr_pile/r_pile.

Fig. 3 Maximum geometrical optical depth profiles of debris disks with collisions, P–R drag, and sublimation (Eq. (28) with τgeo, base(r) given by Eq. (5)), for different values of relative pile-up width Δrpile/rpile shown in different colors. The solid lines correspond to disks around solar-mass stars, the dashed lines to disks around 2 M⊙ stars. All profiles assume a parent belt at r0 = 30 AU, dust grains with β = 0.5, and rpile ≈ 0.019 AU, and rpile ≈ 0.064 AU, for the M⋆ = 1 M⊙, and M⋆ = 2 M⊙ cases, respectively, corresponding to spherical graphite grains with β = 0.5.

We can make a rough estimate of the geometrical optical depth enhancement factor f_τgeo of the pile-up (i.e., how much higher τ_geo in the pile-up is compared to the base level of the inner disk), as a function of the radial width of the pile-up Δr_pile: ${\mathit{f}}_{{\mathit{\tau }}_{\mathrm{geo}}}\mathrm{=}\frac{{\mathit{\tau }}_{\mathrm{geo}\mathit{,}\hspace{0.17em}\mathrm{pile}}\mathrm{+}{\mathit{\tau }}_{\mathrm{geo}\mathit{,}\hspace{0.17em}\mathrm{base}}}{{\mathit{\tau }}_{\mathrm{geo}\mathit{,}\hspace{0.17em}\mathrm{base}}}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{pile}}}{\mathrm{2}\mathrm{\Delta }{\mathit{r}}_{\mathrm{pile}}}\mathrm{+}\mathrm{1.}$(29)For a pile-up width of Δr_pile ≈ 0.05r_pile (Kobayashi et al. 2011), this gives f_τgeo ≈ 11, which is comparable to the highest values found by Kobayashi et al. (2009, see their Fig. 6) for carbonaceous dust grains around early F-type stars.

For comparison with observations, it suffices to assume that all the pile-up material is located at r_pile (i.e., Δr_pile/r_pile ≪ 1). This gives the highest possible temperature to all pile-up particles, and therefore results in the maximum NIR flux. In the τ_geo profile, however, it would give a singularity at r = r_pile. To avoid this, we instead compute the fractional luminosity of the pile-up. Fractional luminosity is defined as the ratio of the infrared luminosity of the dust L_D to the stellar luminosity. Assuming the disk is radially optically thin and the dust grains have unity absorption and emission efficiencies at all wavelengths (i.e., assuming black-body grains), it can be approximated by the fraction of the star that is covered by dust $\frac{{\mathit{L}}_{\mathrm{D}}}{{\mathit{L}}_{\mathit{\star }}}\mathrm{=}\mathrm{\int }\frac{\mathit{\sigma n}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{4}\mathit{\pi }{\mathit{r}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}$(30)Evaluating this with n(r)dr = N_pile, and using the maximum geometrical optical depth (Eq. (5)), we find that the maximum fractional luminosity due to material in the pile-up is $\begin{array}{c}\max \end{array}\left[{\left(\frac{{\mathit{L}}_{\mathrm{D}}}{{\mathit{L}}_{\mathit{\star }}}\right)}_{\mathrm{pile}}\right]\mathrm{=}\frac{\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}}\hspace{0.17em}\mathit{\beta }}{\mathrm{8}\mathit{c}\mathrm{(}\sqrt{{\mathit{r}}_{\mathrm{0}}}\mathrm{-}{\sqrt{{\mathit{r}}_{\mathrm{pile}}}}^{\mathrm{)}}}\mathrm{·}$(31)In the limit of r₀ ≫ r_pile, this becomes $\begin{array}{c}\max \end{array}\left[{\left(\frac{{\mathit{L}}_{\mathrm{D}}}{{\mathit{L}}_{\mathit{\star }}}\right)}_{\mathrm{pile}}\right]\mathrm{\approx }\mathrm{6.2}\mathrm{\times }{\mathrm{10}}^{-6}\hspace{0.17em}{\left(\frac{{\mathit{M}}_{\mathit{\star }}}{{\mathrm{1}\mathit{M}}_{\mathrm{\odot }}}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}{\left(\frac{{\mathit{r}}_{\mathrm{0}}}{\mathrm{1}\mathrm{AU}}\right)}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\hspace{0.17em}\left(\frac{\mathit{\beta }}{\mathrm{0.5}}\right)\mathrm{·}$(32)This is only the fractional luminosity due to the dust in the pile-up. Material just beyond r_pile is not accounted for, and will increase the total fractional luminosity.

2.4. Summary of analytical findings

The analytical model presented above yields several tentative conclusions about dust in the inner parts of debris disks:

• 1.

  P–R drag gives rise to a small but non-zero inward mass flux ofdust in the inner disk, which is self-limited by collisions(Eq. (11)).

• 2.

  A pile-up of sublimating dust occurs, as long as the star is luminous enough (Eq. (18)), and the parent belt is distant enough (Eq. (22)).

• 3.

  P–R drag dominates collisions in the inner parts of the disk (Eq. (25)), so collisions do not interfere with dust pile-up.

• 4.

  Given that the pile-up of dust occurs around the radial distance where the sublimation timescale equals the P–R drag timescale and that sublimating dust resides in the pile-up for about one sublimation timescale, there is a maximum fractional luminosity that this dust can provide (Eq. (32)).

3. Numerical modeling

The analytical approach used in Sect. 2 contains several simplifying assumptions. Most importantly, we only self-consistently solve the continuity equation for P–R drag and collisions (under the assumptions listed in Sect. 2.1), and afterwards estimate the effect of dust pile-up due to sublimation, assuming the grains reside in the sublimation zone for one P–R drag timescale (see Sect. 2.3.2). To test the impact of these assumptions, we now proceed to solve the continuity equation numerically using a debris disk model that self-consistently handles the effects of stellar gravity, direct radiation pressure, P–R drag, sublimation, and destructive collisions. Our strategy here is to simulate a few specific cases and compare the results to the analytical maximum distributions found in Sect. 2, to assess the validity and generality of these simple expressions. A description of our numerical debris disk model is given in Sect. 3.1, the runs we performed are detailed in Sect. 3.2, and the resulting dust distributions are presented in Sect. 3.3.

3.1. Model description

Our debris disk model closely follows the method developed by Krivov et al. (2005, 2006) and Löhne (2008). We refer the reader to these publications for a detailed description of the method and only provide a brief outline of its principles here. We focus on the changes we made, which are including a time-dependent treatment of dust sublimation (see Sect. 3.1.3) and implementing additional numerical acceleration techniques (see Appendix A). Our code was tested by comparing its predictions to solutions of the equations of motion and sublimation for individual particles and by benchmarking it against the results of Krivov et al. (2006). This verification is described in Appendix B. Two physical processes that are included by Löhne (2008), but not considered by our present code, are stellar wind drag and erosive (cratering) collisions. We discuss the impact they may have on our results in Sect. 5.

3.1.1. Method basics

The method of Krivov et al. (2005) applies the kinetic method of statistical physics to debris disks, simultaneously following the spatial and size distributions of dust and planetesimals in a phase space of orbital elements and particle masses. Using orbital element instead of radial distance to follow the spatial distribution makes it possible to account for particles on eccentric orbits, whose orbits can span a wide range of radial distances. The continuity equation is solved in this phase space, with processes that affect the evolution of a particle’s phase-space coordinates in a continuous fashion (P–R drag and sublimation) as diffusion terms, and processes that abruptly change phase-space coordinates (collisions) as source and sink terms. Formally, this is decribed by (cf. Krivov et al. 2005; Löhne 2008) $\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{t}}\mathrm{\left(}\mathit{m,}{k}\mathit{,t}\mathrm{\right)}\mathrm{=}{\left(\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{source}}\mathrm{-}{\left(\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{t}}\right)}_{\mathrm{sink}}\mathrm{-}\mathrm{div}\left(\mathit{n}\frac{\mathrm{d}\mathrm{\left\{}\mathit{m,}{k}\mathrm{\right\}}}{\mathrm{d}\mathit{t}}\right)\mathit{,}$(33)where n(m,k,t) is the phase-space number density at time t (i.e., the distribution function that describes the state of the disk), k the vector of orbital elements, and { m,k } denotes the vector consisting of m and k. The divergence term represents the diffusion of material in phase space (i.e., transport due to P–R drag and sublimation). For brevity, we omit the arguments (m,k,t) for all terms on the righthand side of the equation.

To make the numerical evaluation of the continuity equation manageable with limited computational capacity, the number of phase-space variables needs to be reduced. By assuming the disk is axisymmetric, the distribution function can be averaged over three of the orbital elements: longitude of the ascending node, argument of the periastron, and true anomaly. This implicitly assumes that collisional timescales are much longer than orbital timescales, which generally holds for debris disks (for a more detailed discussion, see Sect. 3.1.3 of Löhne 2008). A further assumption is that the distribution of particles over inclinations is constant, which allows the averaging of the distribution function over inclination. The three remaining phase-space variables are (1) particle mass m; (2) orbital eccentricity e; and (3) an orbital element characterizing the size of the orbit, such as semi-major axis a or periastron distance q = a(1 − e).

The final phase-space variable can be chosen to fit the numerical needs of the problem under investigation. For our study of the pile-up of dust due to sublimation, we chose periastron distance. A particle on an eccentric orbit experiences most sublimation around the periastron, owing to the strong temperature dependence of sublimation. Since the periastron distance does not evolve for β changes that happen at the periastron, the orbit’s periastron distance changes much more slowly than its semi-major axis. Using q instead of a as phase-space variable therefore has numerical advantages.

In practice, the phase space is divided into a grid of bins, and the distribution function is replaced by a vector listing the number of particles in each bin. Equation (33) then becomes a system of ordinary differential equations. The source, sink, and diffusion terms are discretized, and they determine the rates at which the particle numbers evolve, dependent on the population levels of other bins. We now proceed to describe these terms for each of the physical processes considered by our model.

3.1.2. Poynting–Robertson drag

P–R drag affects the orbits of particles, circularizing them, and making them smaller. These effects are accounted for in the model by diffusion terms in the continuity equation that move particles to adjacent bins in the phase-space grid. Since the P–R drag timescale is usually longer than the orbital period, we use the orbit-averaged change rates of the orbital elements, given by (e.g., Burns et al. 1979) $\begin{array}{ccc}{⟨\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{PR}}& \mathrm{=}& \mathrm{-}\frac{\mathit{\beta G}{\mathit{M}}_{\mathit{\star }}}{\mathit{ca}}\frac{\mathrm{2}\mathrm{+}\mathrm{3}{\mathit{e}}^{\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathit{,}\\ {⟨\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{PR}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{5}\mathit{\beta G}{\mathit{M}}_{\mathit{\star }}}{\mathrm{2}\mathit{c}{\mathit{a}}^{\mathrm{2}}}\frac{\mathit{e}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$For the rate of change in periastron distance, we find $\begin{array}{ccc}{⟨\frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{PR}}& \mathrm{=}& \frac{\mathit{\partial q}}{\mathit{\partial a}}{⟨\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{PR}}\mathrm{+}\frac{\mathit{\partial q}}{\mathit{\partial e}}{⟨\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{PR}}\\ & \mathrm{=}& \mathrm{-}\frac{\mathit{\beta G}{\mathit{M}}_{\mathit{\star }}}{\mathrm{2}\mathit{cq}}\frac{\mathrm{(}\mathrm{4}\mathrm{-}\mathit{e}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}{\mathrm{)}}^{\mathrm{3}}}{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{·}\end{array}$

3.1.3. Sublimation

The formalism that we use for dust sublimation is described in Sect. 2.2.1. For a spherical dust particle, it gives a mass loss rate of $\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}{\mathit{P}}_{\mathrm{v}}\mathrm{\left(}\mathit{T}\mathrm{\right)}{\mathit{s}}^{\mathrm{2}}\sqrt{\frac{\mathrm{8}\mathit{\pi \mu }{\mathit{m}}_{\mathrm{u}}}{{\mathit{k}}_{\mathrm{B}}\mathit{T}}}\mathrm{·}$(38)In our numerical model, we use realistic dust grain temperatures (as opposed to the black-body temperatures used in Sect. 2). The method for computing these temperatures is explained in Sect. 3.2.2.

Since the sublimation rate strongly depends on grain temperature, which varies along the path of an eccentric orbit, the mass loss rate needs to be averaged over the orbit. As described qualitatively in Sect. 2.2.2, the change in β ratio associated with mass loss induces changes in orbital elements. Kobayashi et al. (2009) derive the orbit-averaged change rate in orbital elements and mass to be $\begin{array}{ccc}{⟨\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{d}\ln \mathit{\beta }}{\mathrm{d}\ln \mathit{m}}\left(\frac{\mathrm{1}\mathrm{+}{\mathit{e}}^{\mathrm{2}}}{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\mathit{\psi ̅}\mathrm{m}\mathrm{+}\frac{\mathrm{2}\mathit{e}}{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}\mathit{\phi ̅}\mathrm{m}\right)\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\frac{\mathit{a}}{\mathit{m}}\mathit{,}\\ {⟨\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{d}\ln \mathit{\beta }}{\mathrm{d}\ln \mathit{m}}\mathrm{(}\mathit{e}\mathit{\psi ̅}\mathrm{m}\mathrm{+}\mathit{\phi ̅}\mathrm{m}\mathrm{)}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\frac{\mathrm{1}}{\mathit{m}}\mathit{,}\\ {⟨\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}& \mathrm{=}& \mathrm{-}\mathit{\psi ̅}\mathrm{m}\mathit{,}\end{array}$with $\begin{array}{ccc}\mathit{\psi ̅}\mathrm{m}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}}\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{t}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{f,}\\ \mathit{\phi ̅}\mathrm{m}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}}\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{t}}\cos \mathit{f}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}{\mathrm{)}}^{\mathrm{3}\mathit{/}\mathrm{2}}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\cos \mathit{f}{\mathrm{)}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{f,}\end{array}$where f denotes the true anomaly. For the periastron distance, we find $\begin{array}{ccc}{⟨\frac{\mathrm{d}\mathit{q}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}& \mathrm{=}& \frac{\mathit{\partial q}}{\mathit{\partial a}}{⟨\frac{\mathrm{d}\mathit{a}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}\mathrm{+}\frac{\mathit{\partial q}}{\mathit{\partial e}}{⟨\frac{\mathrm{d}\mathit{e}}{\mathrm{d}\mathit{t}}⟩}_{\mathrm{subl}}\\ & \mathrm{=}& \mathrm{-}\frac{\mathrm{d}\ln \mathit{\beta }}{\mathrm{d}\ln \mathit{m}}\frac{\mathit{\psi ̅}\mathrm{m}\mathrm{-}\mathit{\phi ̅}\mathrm{m}}{\mathrm{1}\mathrm{+}\mathit{e}}\frac{\mathit{\beta }}{\mathrm{1}\mathrm{-}\mathit{\beta }}\frac{\mathit{q}}{\mathit{m}}\mathrm{·}\end{array}$For each phase-space bin, quantities and are numerically evaluated using the standard Euler method⁹. The change rates of the phase-space variables (Eqs. (40), (41), and (45)) are then used in diffusion terms in the continuity equation.

Using orbit-averaged mass loss rates is only correct if the sublimation timescale is longer than the orbital period. For phase-space bins for which this does not hold (small particles close to the star), we compute an equilibrium population of particles from the product of their sublimation timescale and the sum of their gain terms. This implicitly assumes that particles are created on their orbits with a uniform distribution over true anomaly.

3.1.4. Collisions

Collisions are different from P–R drag and sublimation in that they cause abrupt rather than smooth changes in phase-space coordinates. In the continuity equation, they are described by sink terms at the phase-space coordinates of targets and projectiles and by source terms at the coordinates of the resulting fragments. Here, we give a summary of the way our model handles collisions. More thorough descriptions of the treatment of collisions, including all relevant equations, are given by Krivov et al. (2006) and Löhne (2008).

For each pair of phase-space bins, we determine collision rates for a range of relative orbit orientations (differences in the longitudes of the periastra of the two orbits). The collision rate is the product of the target and projectile number densities, their relative velocity, the collisional cross-section, and the effective volume of interaction. Krivov et al. (2006) found analytical expressions for these factors in two dimensions as a function of the orbital elements and masses corresponding to both bins and the relative orientation of the orbits. To account for the third dimension, a correction is applied based on the semi-opening angle of the disk, equivalent to the maximum inclinations of the particles. This correction assumes that the disk is relatively flat, consistent with observations of resolved edge-on systems.

To save computational power, we ignore collisions involving unbound particles. This is a valid approximation if the radial geometrical optical depth of the system is much smaller than unity (true for most debris disks), because then the blowout timescale of such particles is much shorter than their collisional timescale¹⁰.

The nature of a collision is determined by the impact energy available per unit of mass. We only consider catastrophic collisions, defined as destructive events in which the largest fragment contains at most half of the mass of the more massive of the two impactors. The threshold for these catastrophic collisions is the critical specific energy for dispersal ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$, which incorporates the fact that fragments may reassemble after destruction, and generally depends on particle size. If the specific energy of a collision is higher than this threshold, the impact destroys both bodies, and their mass is distributed over a swarm of fragments. At specific energies just below ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$, collisions are erosive. Such cratering collisions, however, are not considered in our present model, so if the specific energy of an impact is lower than ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$, no collision is considered to occur.

A catastrophic collision results in a range of fragments with different masses and orbits. The fragments are distributed over particle masses according to a single power law, up to a maximum fragment mass, which is determined by the kinetic energy of the impact and the material strength of the target. The maximum fragment mass is at most half of the mass of the target and projectile combined, but it can also be less, if the specific energy involved in the collision is more than ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$. The amount of particles that end up in each mass bin (up to the maximum fragment mass) is computed by integrating the fragment mass distribution. Particles with masses below the lowest mass bin (i.e., that fall off the grid) are considered lost due to immediate blowout or vaporization. For each fragment mass bin, new orbital elements are calculated using the conservation of momentum, and taking into account direct radiation pressure (i.e., the values of β of the fragments), using Eqs. (19) and (20) of Krivov et al. (2006). This assumes that the fragments are not launched away from the collision with any velocity. These orbital elements are rounded to the nearest periastron distance and eccentricity bins for each fragment mass bin.

3.2. Setup of the model runs

3.2.1. Stellar and disk parameters

Hot exozodiacal dust has been detected around stars with spectral types ranging from A to K (see Table 1). To focus on this range of stellar types, we did one model run with a solar-mass star and one using a 2 M_⊙ star. Following the mass–luminosity relation for main-sequence stars (L_⋆ ∝ M_⋆^3.5; Allen 1976), we set the stellar luminosities corresponding to these stellar masses to L_⋆ = 1 L_⊙ and L_⋆ = 11.31 L_⊙, respectively.

Both runs use a parent belt radius of r₀ = 30 AU. Lower values of r₀ can in principle yield higher dust levels in the innermost regions (see Eq. (5)), but parent belts closer to the star are generally not dense enough to provide these large amounts of dust, because they do not survive the intense collisional grinding (Dominik & Decin 2003; Wyatt et al. 2007). In addition, many of the observed outer belts are located at tens of AUs (see Table 1). The level of dust in the source region is set to τ_geo(r₀) ≈ 5 × 10^-5, chosen such that (i.e., iterated until) the geometrical optical depth in the inner regions does not become any higher by increasing the level of dust at the source¹¹. This roughly corresponds to η₀ ~ 10 for both stellar mass cases, which is apparently enough to approximate an inner disk τ_geo profile corresponding to η₀ ≫ 1. For the semi-opening angle of the disk we use ε = 8.5°.

3.2.2. Material properties

We consider carbonaceous dust particles with a density of ρ_d = 1.8 g cm^-3. To compute the optical properties of these grains with different radii, we use the DHS method of Min et al. (2005) with an irregularity parameter of f_max = 0.8. This method simulates the properties of irregularly shaped particles. For the material we use amorphous carbon with the refractive index data taken from Preibisch et al. (1993). These optical properties are used to compute β(s) and dust temperatures. The dust temperatures are computed by solving the balance between absorption and thermal emission as a function of grain size and distance to the central star. For the sublimation properties, we use those of graphite, given in Sect. 2.2.1.

Modeling collisions requires a prescription for the specific energy threshold for dispersal ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$, which is generally found to depend on size. In our numerical simulation, we consider particles with radii up to 1 cm. (see Sect. 3.2.4). For bodies smaller than s ~ 100 m, ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}$ is often described by a power law with a negative exponent (e.g., Benz & Asphaug 1999). However, such a prescription predicts unrealistically high values for the small particles that we consider. Therefore, following Heng & Tremaine (2010), we use the constant value of ${\mathit{Q}}_{\mathrm{D}}^{\mathit{\star }}\mathrm{=}{\mathrm{10}}^{-7}$ erg g^-1 found in laboratory experiments with high-velocity collisions of small particles (Flynn & Durda 2004). We follow Krivov et al. (2006) in setting the collisional fragment mass distribution to n(m) ∝ m^{−11 /6}, and in assuming the maximum fragment mass scales with specific impact energy to the power −1.24 (Fujiwara et al. 1977).

3.2.3. The phase-space grid

Because this problem is computationally very demanding, great care has to be taken in setting up the phase-space grid. Specifically, resolving the pile-up requires high resolution in the sublimation zone and at small particles sizes, where radiation pressure becomes important. To achieve this with limited computational resources, we designed a non-uniform grid that has a higher resolution where it is required.

The eccentricity grid contains ten logarithmic bins between e = 0 and e = 1, with the lowest bin at e = 10^-4. In addition, there are two linearly spaced eccentricity bins between e = 1 and e = 2 for hyperbolic orbits, as well as two bins between e = −2 and e = −1 to account for “anomalous” hyperbolic orbits followed by β> 1 particles (see Krivov et al. 2006).

The periastron distance grid consists of two parts: (1) a high-resolution, linear grid of 21 bins is used to cover the sublimation zone (0.01 AU <q< 0.03 AU for the M_⋆ = 1 M_⊙ run, and 0.05 AU <q< 0.1 AU for the M_⋆ = 2 M_⊙ run). (2) A low-resolution, logarithmic grid covers the rest of the disk out to about q = 100 AU, with 60 bins in the M_⋆ = 1 M_⊙ run and 50 bins in the M_⋆ = 2 M_⊙ run. Care was taken to place one bin exactly at q = 30 AU, which was used as the source region.

The mass grid has 48 logarithmically spaced bins in both runs with higher resolution at the smaller sizes (β ≳ 0.05) and a maximum mass corresponding to s = 1 cm. For the M_⋆ = 1 M_⊙ run, the high-resolution part consists of 30 bins between s = 0.5 μm and s = 10 μm. The high-resolution part of the M_⋆ = 2 M_⊙ run has 36 bins between s = 2 μm and s = 100 μm.

3.2.4. Simulation strategy

Owing to computational limitations, the largest particles we consider have a radius of 1 cm. In reality, the size distribution in the parent belt extends up to planetesimals of tens to hundreds of kilometers. To account for the fragmentation of these larger bodies, we include a source of dust at r₀. This artificial source term adds particles with sizes between s = 1 mm and s = 1 cm, following the power-law size distribution n(s) ∝ s^-3.5. The size of these grains is chosen such that the effect of radiation pressure on them is very small (β< 10^-4). Therefore, their eccentricity distribution follows that of the parent bodies. We add the source particles at eccentricities ranging from e = 0 to e = 0.1.

The artificial supply of large dust particles is balanced by the loss of material due to blowout and sublimation. Therefore, solving the continuity equation results in a steady-state distribution function. We start the integration without any material in the disk (only the artificial supply is acting) and let the model run until steady state is reached. This is considered to be the case when relative changes in the radial and size distributions between logarithmic (base-10) time steps become less than 1%¹². Initially, we only consider collisions and P–R drag, and the sublimation module of the code is switched off. At this stage, particles migrate inward due to P–R drag until they reach the inner edge of the grid. Once steady state is reached, sublimation is switched on, and we let the distribution function settle into a new steady state. This procedure is necessary because sublimation forces the time step to become very short. With sublimation switched on from the start, the computation would take unnecessarily long. Additionally, it allows us to isolate the effect of dust pile-up (i.e., the dust in the pile-up can be isolated by subtracting the pre-sublimation state from the final one).

3.3. Results

The output of each model run is a steady-state distribution of particles in the phase space of orbital elements and masses. To analyze this output, we convert it into a radial geometrical optical depth profile (Sect. 3.3.1) and size distributions (in terms of cross-section density per unit size decade A) at different radial locations (Sect. 3.3.2). The conversion from raw model output to these quantities is detailed in Appendix C.

3.3.1. Radial distribution

Figure 4 shows the geometrical optical depth profiles derived from the numerical model runs, together with the analytical maxima given by Eq. (28), using a pile-up width of Δr_pile/r_pile = 0.15, chosen to match the numerical profile. Generally, there is good correspondence between the numerical results and the analytical maxima, but there are some important differences. The profiles have roughly the same shape, with a slightly steeper slope close to the source region in the numerical results. For both cases of stellar mass, the base level of τ_geo in the inner disk (i.e., away from the pile-up) is a factor of about 7 lower in the numerical profiles. This discrepancy is a result of the the assumption in the analytical model that all orbits are circular. In the numerical model, the small particles that contribute most to the cross-section are released in the parent belt on eccentric orbits. They therefore suffer a higher rate of destruction by collisions.

Fig. 4 Geometrical optical depth profiles of debris disks with a dense parent belt at 30 AU around stars of M⋆ = 1 M⊙ (solid lines) and M⋆ = 2 M⊙ (dashed lines). The black lines show the end results of the numerical simulations. In green are the maximum τgeo profiles as given by the analytical model (Eq. (28), with τgeo, base(r) given by Eq. (5), Δrpile/rpile = 0.15, r0 = 30 AU, and β = 0.5).

As predicted, switching on sublimation leads to an accumulation of dust. The pile-ups are located very close to r_pile, as determined for graphite grains, using black-body temperatures and Q_pr = 1 (Eq. (17), which gives r_pile ≈ 0.019 AU for the M_⋆ = 1 M_⊙ case, and r_pile ≈ 0.064 AU for the M_⋆ = 2 M_⊙ case). The temperature of the dust (as computed using the full optical properties) at the inner edge of the disk and in the pile-up is between 2000 and 2100 K. In both runs, the pile-up has a geometrical optical depth enhancement factor of about f_τgeo ≈ 3. The fractional luminosities due to the pile-ups are (L_D/L_⋆)_pile ≈ 7.5 × 10^-8 for the M_⋆ = 1 M_⊙ case and (L_D/L_⋆)_pile ≈ 1.1 × 10^-7 for the M_⋆ = 2 M_⊙ case. Both are about a factor 15 lower than the maxima given by Eq. (32). Given that the base levels of τ_geo in the numerical profiles are a factor of about 7 lower than the analytical maxima, however, the discrepancy is only about a factor of 2. In short, the pile-up mechanism is found to be somewhat less efficient than predicted by the analytical estimates.

The outer disk (r> 30 AU) is not the focus of this work. Nevertheless, its radial profile is relevant, since it reflects the status of the balance between collisions and P–R drag. To first order, the geometrical optical depth profile of the outer disk can be characterized by a power law τ_geo ∝ r^−α. Strubbe & Chiang (2006) derive the theoretical values of α = 1.5 and α = 2.5 for collision and P–R drag-dominated disks, respectively. We find slopes of α ≈ 2.0 for both runs, consistent with the outer slope found by Vitense et al. (2010) for the Edgeworth–Kuiper Belt, and interpreted as the sign of a disk that is in between drag and collision-dominated. This indicates that the density of the parent belt (characterized by η₀ ~ 10) is insufficient to make the outer disk completely collision-dominated.

Fig. 5 Size distributions at different radial distances for the M⋆ = 1 M⊙ run before switching on sublimation. The quantity on the vertical axis, A, is the cross-section density per unit size decade (see Appendix C.2), which is horizontal if all sizes contribute equally to the cross-section. The gray band indicates the slope of a size distribution that follows the classical Dohnanyi (1969) power law (n(s) ∝ s-3.5, hence A ∝ s-0.5). The vertical lines mark the particle sizes corresponding to relevant values of β.

3.3.2. Size distribution

The size distribution results of the two runs are similar in many ways, so we only discuss the M_⋆ = 1 M_⊙ run here. Figure 5 shows how the size distribution changes with radial distance when only considering P–R drag and collisions (i.e., before sublimation is switched on in the model). In the parent belt (r = 30 AU), it follows the classical Dohnanyi (1969) power law (n(s) ∝ s^-3.5, which is valid for an infinite collisional cascade with self-similar collisions), from the blowout radius upwards. Particles with β> 0.5 are depleted by about three orders of magnitude in terms of collective cross-section. Superimposed on the power law is a well-known wave pattern related to the discontinuity in the size distribution at the blowout size (see, e.g., Campo Bagatin et al. 1994; Durda & Dermott 1997; Gáspár et al. 2012). The first bump (i.e., the one at β ≈ 0.5) in the size distribution at r = 30 AU does not extend far above the power-law prediction. The reason for this may be that particles with β ≳ 0.1 experience more destructive collisions, because their eccentricities are significantly higher than those of the parent bodies, which are distributed over the range 0 <e< 0.1 (Sect. 3.2.4, see also Eq. (19), and cf. Fig. 5 of Krivov et al. 2006).

Inward from the parent belt, the size distribution seems to become steeper, which is expected from the dependence of the radial profile on β (Eq. (5)). However, this effect is difficult to isolate, because the profile is distorted by the wave pattern, which increases in amplitude and “wavelength” with decreasing radial distance (cf. Fig. 7 of Krivov et al. 2006). The prominent wave pattern indicates that collisions are still important for larger particles in the inner disk (while P–R drag dominates for β ≈ 0.5 particles there, see Fig. 2). In the innermost parts of the disk (r ≲ 1 AU), particles with β ≈ 0.5 clearly dominate the cross-section, contributing at least three orders of magnitude more than any other size. At r = 0.1 AU, the local slope of the size distribution between the blowout size and s ≈ 10 μm is approximately A ∝ s^-7.5, equivalent to n(s) ∝ s^-10.5. This means that not only the cross-section, but also the mass is dominated by barely bound grains in the innermost parts of the disk. Interestingly, such steep size distributions are also invoked to explain NIR interferometric observations of hot exozodiacal dust (Defrère et al. 2011; Mennesson et al. 2013; Lebreton et al. 2013). The drop in the size distribution from β ≈ 0.5 to β ≈ 1 is also much more pronounced in the inner disk than in the parent belt.

Sublimation only has a significant effect on the size distribution around r ≈ r_pile (i.e., in the pile-up). Figure 6 shows the size distribution in the pile-up, before and after sublimation is switched on. It clearly demonstrates that sublimation enhances the density of particles with 0.5 ≲ β< 1 around r = r_pile. This indicates that the pile-up consists mostly of particles that started with β ≈ 0.5 before active sublimation and lost mass due to sublimation, increasing their β. The rest of the size distribution does not change significantly. Size distributions at larger radial distances are largely unaffected by sublimation.

Fig. 6 Size distribution at r = 0.02 AU (the location of pile-up) for the M⋆ = 1 M⊙ run, before and after sublimation is switched on.

Fig. 7 SEDs of stars and their debris disks at a distance of 10 pc. The disk SED is shown for different dust distributions with the numerical results in black and the analytical distributions in green. Dust distributions including the pile-up of dust in the sublimation zone are shown with dashed lines, and distributions in which the pile-up is excluded are shown with solid lines, but the spectra largely overlap. Disk SEDs are calculated according to Eqs. (46) and (47). The stellar photosphere is indicated by a dotted Planck curve. The vertical gray areas mark the NIR H and K bands in which hot exozodiacal dust is observed.

4. Comparison with observations

To assess whether the pile-up effect can explain the observed NIR excess, we computed the spectral energy distributions (SEDs) of the dust distributions found in the previous sections. We only calculated the emission spectrum of the dust and ignored the (viewing-angle-dependent) contribution of scattered light, since thermal emission was found to dominate scattered light at the wavelengths in which hot exozodiacal dust is detected (Absil et al. 2006). For the analytical dust distributions, we assumed the dust has a black-body temperature (Eq. (14)). The flux density of the dust, expressed as a function of wavelength λ, can then be computed as ${\mathit{F}}_{\mathit{\nu ,}\mathrm{D}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathit{\lambda }}^{\mathrm{2}}}{\mathit{c}{\mathit{d}}^{\mathrm{2}}}\mathrm{\int }\mathit{r}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{B}}_{\mathit{\lambda }}\mathrm{\left(}{\mathit{T}}_{\mathrm{bb}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{r,}$(46)where d is the distance to the source, and B_λ(T) the Planck function. For the numerically determined dust distributions, we use realistic dust temperatures and optical properties (see Sect. 3.2.2), which leads to ${\mathit{F}}_{\mathit{\nu ,}\mathrm{D}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }{\mathit{\lambda }}^{\mathrm{2}}}{\mathit{c}{\mathit{d}}^{\mathrm{2}}}\underset{}{{\mathrm{\int }}^{\mathit{s}}}\underset{}{{\mathrm{\int }}^{\mathit{r}}}\mathit{r}{\mathit{Q}}_{\mathrm{abs}}\mathrm{\left(}\mathit{s}\mathrm{\right)}\hspace{0.17em}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{s,r}\mathrm{\right)}{\mathit{B}}_{\mathit{\lambda }}\mathrm{\left[}\mathit{T}\mathrm{\right(}\mathit{s,r}\mathrm{\left)}\mathrm{\right]}\hspace{0.17em}\mathrm{d}\mathit{s}\hspace{0.17em}\mathrm{d}\mathit{r}\mathit{.}$(47)Here, Q_abs is the absorption efficiency of the grains (equal to their emission efficiency) and τ_geo(s,r) is the geometrical optical depth profile of dust with grain radius s.

Figure 7 shows the debris disk spectra corresponding to dust distributions found from the numerical modeling, as well as the analytical maximal dust profiles and the stellar spectra. The analytical maxima are computed from Eq. (28), with τ_{geo, base}(r) given by Eq. (5), r₀ = 30 AU, β = 0.5, and all the pile-up material at r_pile (i.e., Δr_pile/r_pile ≪ 1)¹³. Singularities at r = r₀ are removed by imposing the condition τ_geo ≤ 0.01, which corresponds to η₀ ≈ 2200. The numerical dust distributions without pile-up were created by subtracting the isolated pile-up dust from the final profile (see Sect. 3.2.4). The stellar spectra are given by black-body curves, with L_⋆ = 1 L_⊙, T_⋆ = 5780 K for the M_⋆ = 1 M_⊙ star, and L_⋆ = 11.31 L_⊙, T_⋆ = 7730 K for the M_⋆ = 2 M_⊙ star, following the main-sequence relations L_⋆ ∝ M_⋆^3.5 and T_⋆ ∝ L_⋆^0.12 (Allen 1976). The distance is arbitrarily set to 10 pc.

The synthetic spectra are similar to those described by Wyatt (2005), but are truncated at short wavelengths because there is no material with a higher temperature than the sublimation temperature. Apart from an overall shift in flux, the differences between the analytical and numerical SEDs are minor. In the solar-mass star run, the peak of emission of the numerical result is shifted toward shorter wavelengths with respect to the analytical SED. This is because the grains that dominate the emission in the parent belt are significantly hotter than the black-body temperature.

Only considering thermal emission, the NIR flux originates almost exclusively in the inner 1 AU of the debris disk. The pile-up does not add a significant amount of flux. The theoretical SEDs display NIR flux ratios between the disk and star of about F_ν,D/F_ν,⋆ ~ 10^-4. This is much less than the observed flux ratios that indicate hot exozodiacal dust, which are typically on the order of 1% (see Table 1). Furthermore, the NIR radiation is accompanied by mid-IR flux at a comparable level, which is incompatible with observed excess spectra (e.g., Akeson et al. 2009). The pile-up does not contain enough material to create a bump in thermal flux in the NIR.

To generalize the above results, we investigate how the NIR flux ratio depends on parent belt distance and stellar type using the analytical model. In Fig. 8, we show the analytical maximum NIR flux ratios for four different stellar types. As in the above analysis, the disk flux is calculated from Eq. (46), using the analytical maximum dust distribution, and the stars are approximated by black bodies. The M0V, G2V, A5V, and B5V stars correspond to M_⋆ = 0.5, 1, 2, and 4 M_⊙, respectively, with L_⋆ ∝ M_⋆^3.5 and T_⋆ ∝ L_⋆^0.12 (Allen 1976). Comparing the analytical maximum NIR flux ratios to the observed ones demonstrates that P–R drag does not provide enough material to the inner disk to explain the observations.

Fig. 8 Flux ratio between the disk and the star at 2.1 μm (K band), as function of parent belt radius r0. Solid lines indicate the flux ratio due to the maximum amount of material moved in by P–R drag, truncated at r = rpile. Ratios including the small effect due to pile-up are shown with dashed lines. Different colors are used for different stellar types. The light shaded area marks the region in parameter space where no significant pile-up is expected to occur because the orbits of small dust grains cannot circularize sufficiently (see Eq. (22)). The dark shaded area marks where the parent belt r0 is closer than the pile-up radius rpile. The observed K-band excess fluxes of main-sequence stars from Table 1 are shown at their respective estimated parent belt distances with a different symbol for each star and coloring according to the star’s spectral type: Vega (the CHARA/FLUOR measurement): circles; β Leo: hexagon; Fomalhaut: squares; β Pic: diamond; η Lep: upward pointing triangles; 110 Her: left-pointing triangle; 10 Tau: right-pointing triangle; τ Cet: pentagon.

5. Discussion

5.1. The size distribution in the inner disk

The size distribution in the inner parts of dense debris disks (i.e., inward from the parent belt) has not been studied many times before. Acke et al. (2012) analytically derive a power-law distribution of n(s) ∝ s^-3, but ignore collisions in the inner disk. In a detailed modeling study of the debris disk around ε Eri, Reidemeister et al. (2011) find a flat profile (n(s) ∝ s^-3, flat in terms of A) for small sizes, and a steeper one (n(s) ∝ s^-3.7) for larger sizes. However, this system is a special case, because ε Eri is not luminous enough to blow small dust grains out of the system.

With our numerical model, we find a size distribution that deviates significantly from a power law. Specifically, in the innermost regions of the disk (r ≲ 1 AU), the cross-section is completely dominated by barely bound (β ≈ 0.5) particles. This is the result of two effects: (1) Since the efficiency of P–R drag depends on β, larger particles tend to stay closer to the parent belt (apparent from Eq. (5)). (2) Owing to the high relative velocities in the inner disk¹⁴, the wave in the size distribution is extremely strong. As a consequence, the single size assumption used in our analytical model is a good approximation.

Erosive (cratering) collisions affect the size distribution (Thébault & Augereau 2007). We expect that including this type of collisions in our model would result in a less wavy size distribution, possibly eliminating the second bump in the size distribution (see Sect. 4.3.4 of Löhne 2008). This would mean that β ≈ 0.5 particles dominate the cross-section even more than suggested by Fig. 5. Since erosive collisions present an additional mechanism for destroying dust, we expect that including them in our model would only lower the level of dust in the inner disk, so it does not change the conclusions of this work. Erosive collisions do affect the timescale on which debris disks evolve, but since we compute steady-state dust distributions, this does not have any impact on our results.

It may be surprising that a large number of particles are present in the pile-up with 0.5 ≲ β< 1. Usually (as in Sect. 2), particles with β> 0.5 are assumed to be absent, because they are blown out of the system as soon as they are released from parent bodies on circular orbits. Parent bodies on eccentric orbits can give rise to a population of bound grains with β> 0.5, but the parent body orbits in our model runs are not eccentric enough to produce bound dust grains with β values close to unity. The origin of this high β population is the sublimation of barely bound β ≈ 0.5 grains. As these particles migrate inward from the parent belt, their orbits are circularized by P–R drag. When they arrive in the sublimation zone, their orbital eccentricities are as low as e ≈ 10^-4. Kobayashi et al. (2009) find that the subsequent evolution of the eccentricities during the active sublimation phase can be described by $\mathit{e}\mathrm{=}{\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\beta }}_{\mathrm{1}}}{\mathrm{1}\mathrm{-}\mathit{\beta }}\right)}^{\mathit{\kappa }}{\mathit{e}}_{\mathrm{1}}\mathit{,}$(48)where subscript 1 denotes quantities at the start of substantial sublimation, and the exponent κ can be treated as a constant, which mostly depends on the optical and sublimation properties of the material under consideration. For spherical, black-body graphite particles, we find κ ≈ 10. This shows that particles starting with β₁ ≈ 0.5, and e₁ ≈ 10^-4 will only become unbound (e ≥ 1) when they reach β ≳ 0.8 (corresponding to s ≈ 0.8 μm).

5.2. Stellar wind drag

Neither our analytical nor our numerical model includes the effects of stellar wind. For stars with a strong stellar wind and/or low luminosity, the stellar wind equivalent of P–R drag shortens the inward migration timescale and thus increases the maximum geometrical optical depth profile. Since stellar wind drag works in the same way as P–R drag, its effect can be accounted for by replacing β with (Burns et al. 1979; Minato et al. 2004, 2006) $\begin{array}{ccc}{\mathit{\beta }}_{\mathrm{PR}}& \mathrm{=}& \mathit{\beta }\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\mathit{,}\\ \mathit{\gamma }& \mathrm{=}& \frac{\mathit{Ṁ}\mathit{\star }{\mathit{c}}^{\mathrm{2}}}{{\mathit{L}}_{\mathit{\star }}}\frac{{\mathit{Q}}_{\mathrm{sw}}}{{\mathit{Q}}_{\mathrm{pr}}}\mathit{,}\end{array}$where Ṁ_⋆ is the stellar mass loss rate by stellar wind, and Q_sw is the stellar wind momentum transfer efficiency. For carbonaceous particles orbiting the Sun, Minato et al. (2004, 2006) find γ ~ 1.

Kobayashi et al. (2011) argue that the pile-up scenario could work for Vega if the disk is drag-dominated, which requires γ ≈ 300. While this value is consistent with the upper limit on the mass-loss rate of Vega from radio-continuum observations (Ṁ_⋆ ≲ 10^-10 M_⊙ yr^-1; Hollis et al. 1985), stellar wind models predict much lower mass-loss rates for main-sequence A-type stars (Ṁ_⋆ ≲ 10^-16 M_⊙ yr^-1; Babel 1995). This theoretical mass-loss rate, together with L_⋆ ≈ 10 L_⊙ and Q_sw ≈ Q_pr, gives a ratio between stellar wind drag and P–R drag of γ ≲ 10^-4. We therefore conclude that it is unlikely that stellar wind drag has a significant effect on debris disks around main-sequence A-type stars.

5.3. Other explanations for hot exozodiacal dust

We find that P–R drag does not provide enough material to the innermost parts of the disk to explain the interferometric detections of NIR excess. There are two additional problems with this scenario, which also serve as clues for solving the hot exozodiacal dust mystery. (1) Hot exozodiacal dust is thought to consist mostly of blowout grains with sizes around 0.01–0.1 μm (Akeson et al. 2009; Defrère et al. 2011; Mennesson et al. 2013; Lebreton et al. 2013), while P–R drag only transports bound grains to the sublimation zone, the smallest of which have sizes of about 1 μm. (2) The dust distribution resulting from the balance between P–R drag and collisions yields an SED with a positive slope in the infrared domain, while observations find negative slopes (e.g., Akeson et al. 2009; Acke et al. 2012), and the pile-up of dust is too inefficient to have an effect on the slope of the SED¹⁵. For these reasons, the origin of hot exozodiacal dust must involve mechanisms that we did not consider in this work.

Our treatment of sublimation assumes spherical dust grains that sublimate uniformly (i.e., layers of material are removed one by one). In reality, dust grains may be aggregates that fall apart during sublimation, abruptly increasing the collective cross-section of the material. To investigate this scenario, the steady-state amount of material can be estimated by multiplying the maximum P–R drag inward mass flux (Eq. (11)) with an estimate of the lifetime of the fragments. We performed this analysis for Fomalhaut and find that P–R drag from a parent belt at 2 AU still does not supply enough material, unless the lifetime of the fragments is significantly longer than what can be expected from sublimation and blowout (Lebreton et al. 2013).

Two mechanisms have recently been proposed that can lead to an extended lifetime for small exozodiacal dust particles. Lebreton et al. (2013) investigate the impediment of blowout due to the presence of gas for the hot dust around Fomalhaut, but find that this requires unrealistically high gas densities. Su et al. (2013) propose that charged nanograins can remain trapped in the magnetic field of the star and qualitatively show that this may help explain the NIR excess of Vega.

Another mechanism for the inward transport of material from a cold outer belt is the inward scattering of material by planets. Bonsor et al. (2012) investigated this scenario and find that it is marginally capable of providing mass influxes compatible with observations. However, this requires relatively contrived planetary system architectures, consisting of closely packed chains of low-mass planets. Furthermore, the scenario was investigated for the inward scattering down to 1 AU, and reaching the sublimation zone is likely to be less efficient.

6. Conclusions

In this work, we investigated hot dust in the inner regions of debris disks, whose presence is suggested by interferometrically resolved excess NIR emission observed in some debris disk systems (Table 1). We tested whether the hot dust can be supplied by P–R drag from a distant parent belt and whether the pile-up of dust in the sublimation zone still occurs if collisions are considered. Our main conclusions follow.

• 1.

  As predicted by Wyatt (2005), P–R drag al-ways brings a small amount of dust from an outer debris beltinto the sublimation zone. The maximum geometrical opticaldepth that can be reached in the innermost parts of the disk de-pends on the mass of the central star and distance to the parent belt(Fig. 1). When the production of dust is treatedself-consistently, this maximum is found to be a factor of about 7lower than the analytical estimate (Fig. 4). Thisis because small dust particles, which are efficiently draggedinward by radiation forces, are also put on highly eccentric orbitsby those radiation forces and therefore suffer more collisionaldestruction.

• 2.

  Dust that reaches the sublimation zone produces some NIR emission, but this excess flux is insufficient to explain the interferometric observation. While the observed excess ratios are ~10^-2, the maximum flux ratio due to material supplied by P–R drag is ≲10^-3 for A-type stars with parent belts at ≳1 AU (Fig. 8).

• 3.

  The pile-up of dust from the interplay of P–R drag and sublimation still occurs when collisions are considered (Fig. 4), as long as the parent belt in which the dust originates is distant enough to allow for sufficient circularization of the orbits, and the central star is luminous enough to blow small dust grains out of the system. Collisions do not interfere with the pile-up process, since in the inner disk, the collisional timescale is longer than the P–R drag timescale for the barely bound grains that are the most important for the pile-up. The fractional luminosity provided by dust in the pile-up is relatively low, so the pile-up does not influence the disk SED significantly (Fig. 7).

• 4.

  In the inner parts of dense debris disks, the cross section is clearly dominated by barely bound (β ≈ 0.5) grains, and the size distribution features a prominent wave pattern, related to the discontinuity in the size distribution at the blowout size (Fig. 5). In the pile-up, there is an enhancement of particles with 0.5 ≲ β< 1 (Fig. 6). These particles are still bound, because of their almost circular orbits at the start of substantial sublimation.

Online material

Appendix A: Numerical techniques

Our model is computationally very demanding, primarily because of the large number of calculations needed for collisions, which scales with the number of phase-space bins cubed. To ensure that the model can produce results in a reasonable amount of time, an effort was made to reduce the amount of calculations needed per run. To achieve this, we used an integration technique that allows large time-step sizes, which reduces the number of steps needed per run (Appendix A.1). Furthermore, we calculated as many time-independent numerical factors as possible before the actual integration, so these factors do not have to be determined at every time step (Appendix A.2).

Appendix A.1: Integration method

After discretization, the distribution function n(m,k,t) is replaced by the vector of population levels n, and Eq. (33) can be treated as the system of linear equations $\frac{\mathrm{d}{n}}{\mathrm{d}\mathit{t}}\mathrm{=}{A}{n}\mathrm{+}{b}\mathit{,}$(A.1)where A is a matrix of coefficients, and b is a constant vector containing the artificial source terms that replenish dust in the parent belt¹⁶. This system of equations suffers from stiffness: the population levels of some bins change very rapidly compared to others mainly due to large differences in collisional timescales, and the time-step size is therefore determined by the stability of the integration method rather than by its accuracy. When using standard explicit integrators, this leads to an impractically small step size that prevents long integrations.

Following Löhne (2008), we resolve the stiffness by writing the differentials as $\frac{\mathrm{d}{n}}{\mathrm{d}\mathit{t}}\mathrm{=}{{A}}^{{\prime }}\hspace{0.17em}{n}\mathrm{+}{\frac{\mathrm{d}{n}}{\mathrm{d}\mathit{t}}|}_{\mathrm{const}}\mathit{,}$(A.2)where A^′ is a diagonal matrix that contains only the diagonal elements of A, while the terms marked “const” contain the off-diagonal parts of A, as well as the constant terms b. Our integration scheme for the jth component of n, for timestep m + 1, reads as${\mathit{n}}_{\mathit{j,m}\mathrm{+}\mathrm{1}}\mathrm{=}{\mathit{n}}_{\mathit{j,m}}\exp \mathrm{(}{\mathit{A}}_{\mathit{jk,m}}^{\mathrm{\prime }}\mathrm{\Delta }{\mathit{t}}^{\mathrm{)}}\mathrm{+}\frac{\mathit{ṅ}\mathit{j,m}{\mathrm{|}}_{\mathrm{const}}}{{\mathit{A}}_{\mathit{jk,m}}^{\mathrm{\prime }}}\mathrm{[}\exp \mathrm{(}{\mathit{A}}_{\mathit{jk,m}}^{\mathrm{\prime }}\mathrm{\Delta }{\mathit{t}}^{\mathrm{)}}\mathrm{-}{\mathrm{1}}^{\mathrm{]}}\mathit{,}$(A.3)where Δt denotes the time-step size¹⁷. This scheme is known as the exponential Euler method (see Hochbruck & Ostermann 2010 for an introductory review). It is suitable for time integration of semi-linear problems, which consist of a stiff linear part and a non-stiff non-linear part. In short, the strategy is to solve the linear part exactly and to approximate the non-linear part using an explicit integration scheme.

The time-step size Δt is determined dynamically from the condition that population levels can never become negative. We use a scheme similar to that of Krivov et al. (2005), but adapted for the exponential Euler method.

Appendix A.2: Precalculation

Only considering collisions (i.e., ignoring diffusion terms), the evolution of the jth component of n can be written as (Löhne 2008) $\mathit{ṅ}\mathit{j}\mathrm{=}\sum _{\mathit{tp}}{\mathit{B}}_{\mathit{jtp}}{\mathit{n}}_{\mathit{t}}{\mathit{n}}_{\mathit{p}}\mathit{.}$(A.4)Here, t and p are the target and projectile bin indices, respectively, and B_jtp is an entry in the time-independent tensor of collisional coefficients B. Specifically, coefficient B_jtp is the rate at which the population level of bin j changes, per particle in bin t, per particle in bin p, combining all considered relative orbit orientations. If j = t or j = p, this is the loss rate of the target or projectile bin, respectively, caused by collisional destruction (assuming that the mass grid resolution is high enough that fragments do not end up in their parent bin). Otherwise, it is a rate at which fragments are created.

For the phase-space grids we use, the entire tensor B is too large to be stored in memory. However, it is very sparse, because (1) not all orbits that are part of the phase-space grid have mutual overlap; (2) overlap may occur for a limited range in relative orbit orientation; (3) impact velocities are not always high enough to cause catastrophic collisions; and (4) only a fraction of all possible masses are created as fragments.

By only storing the non-zero entries of B, it becomes possible to keep it in memory. We store the source and sink terms separately. For the source terms, which are by far the most memory intensive, each non-zero entry requires (1) the index of the target bin; (2) the index of the projectile bin; (3) the index of the bin fragment bin; and (4) the rate at which particles are created in the fragment bin, per particle in the target bin, per particle in the projectile bin. By further manual compression, only (3) and (4) need to be stored for each entry separately. The target bin index (1) only needs to be stored once for each possible target bin, along with the number of possible projectile bins for that target bin. Then, the projectile bin index (2) only needs to be stored once for each possible collisional pair of bins, along with the number of possible fragment bins for that pair. A similar scheme is used for the sink terms.

The compressed version of B is still very large (several gigabytes for the runs presented in this work). A disadvantage of the precalculation technique is therefore that the size of the grid that can be used is restricted by the available amount of memory, whereas a code that recalculates (parts of) B at each time step is only limited by CPU power. A cubic dependence of the non-zero entries of B on the number of mass bins (as opposed to quadratic dependencies on the orbital element variables) motivates the choice of a relatively small mass grid, representing only a part of the collisional cascade.

Appendix B: Model verification

We performed several tests to verify our numerical debris disk model. The P–R drag and sublimation modules of the code were tested by comparing the behavior of the model with independent numerical solutions of the equations of motion and sublimation for individual particles. For this purpose, the collisional module of the code was switched off, and a single bin was filled as the initial setup of the model, corresponding to the orbital elements of the particle. For all these tests, the resulting evolution of the dust distribution (not shown here) matches that of the independent solution. The accuracy of the predictions is limited by numerical diffusion and becomes better with higher resolution (i.e., larger phase-space grids).

To test the collisional module of our code, we benchmarked the model against the code of Krivov et al. (2006), by replicating one of their model runs for the debris disk of Vega as accurately as possible. Since the code of Krivov et al. (2006) uses semi-major axes as the distance dimension of the phase-space grid, the benchmark runs were performed with a version of our code that also uses the semi-major axis (as opposed to periastron distance, used in the rest of this paper). These runs do not include the effects of P–R drag or sublimation, so they can be used to separately test the collisional module of our code by switching off P–R drag and sublimation. Of the various runs presented by Krivov et al. (2006), the specific one that was reproduced is characterized by an initial optical depth profile in the outer disk (beyond 120 AU) that scales as τ_geo(r) ∝ r^-4, an initial eccentricity distribution between 0 and 0.375, and material properties for “rocky” grains. We refer the reader to Krivov et al. (2006) for the specific values used for parameters describing the phase-space grid, the initial setup of the disk, material properties, etc.

Fig. B.1 Evolution of the radial geometrical optical depth profile of the benchmark run, to be compared with Fig. 10 of Krivov et al. (2006).

Fig. B.2 Evolution of the size distribution at r = 100 AU of the benchmark run, to be compared with Fig. 6 of Krivov et al. (2006).

The evolution of the radial and size distribution predicted by the benchmark run are shown in Figs. B.1 and B.2, respectively. The corresponding results of Krivov et al. (2006) are their Figs. 10 and 6, respectively. Comparison of the results reveals good agreement between the two codes, and the remaining discrepancies can be accounted for. Relative to the benchmark, our model predicts (1) lower unbound particle populations and (2) shorter evolutionary timescales. Point (1) is to be expected since the unbound particle densities of Krivov et al. (2006) are computed using the product of their production rate and their disk crossing time, rather than from Eqs. (C.6) (Löhne, priv. comm.). We attribute point (2) to a mislabeling of the time steps in Figs. 6 and 10 of Krivov et al. (2006). Further discrepancies can be due to small differences in the input parameters and numerical techniques used.

Appendix C: Post-processing of model output

The output of a model run are the population levels n of all bins in the three-dimensional phase space of particle mass and orbital elements, as well as their change rates ṅ (which should equal zero once steady state is reached, except for bins corresponding to unbound orbits). While this data is useful for analyzing the orbital characteristics of different classes of particles in the debris disk, it is often more convenient to know about the state of the disk as a function of distance from the star. This is essential, for example, if the results of the model are to be compared with observations. Here, we describe the processing steps that are applied to the model output to derive the quantities used in Sect. 3.3.

Appendix C.1: Conversion from orbital elements to radial distance

To find the radial distribution of matter in the debris disk, the orbital element phase-space distribution function n(m,q,e) (dimension: [g^-1 cm^-1]) needs to be converted to the configuration space distribution function , which denotes the vertically averaged number density of particles with masses [m,m + dm] at distance r (dimension: [g^-1 cm^-3]). This problem was first solved by Haug (1958) for a rotationally symmetric ensemble of particles on Keplerian orbits. Here, we give a brief derivation under the additional assumption that the distribution of inclinations is independent of the distribution of the other orbital elements.

Consider an individual particle on a bound Keplerian orbit that spends dt time to cross a radial annulus with width dr at distance r from the star. The contribution of this particle to is ${\mathrm{𝒩}}_{\mathrm{part}}\mathrm{=}\frac{\mathrm{2}\mathrm{d}\mathit{t}}{\mathit{P}}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi r}\hspace{0.17em}\mathrm{d}\mathit{r}}\frac{\mathrm{1}}{\mathit{h}}\mathit{,}$(C.1)where P is the particle’s orbital period, and h = 2r sin ε is the disk height, where ε is the semi-opening angle of the disk. The explicit factor 2 in the numerator accounts for the fact that the particle passes through this radial annulus twice during each orbit. In terms of orbital elements q and e, the orbital period P, accounting for direct radiation pressure, is given by $\mathit{P}\mathrm{=}\mathrm{2}\mathit{\pi }\sqrt{\frac{{\mathit{q}}^{\mathrm{3}}}{\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}{\mathrm{)}}^{\mathrm{3}}}}\mathrm{·}$(C.2)The radial velocity ṙ of the particle is given by $\frac{\mathrm{d}\mathit{r}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{\pm }\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}}{\mathit{r}}\left[\mathrm{2}\mathrm{-}\frac{\mathit{r}}{\mathit{q}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}\mathrm{-}\frac{\mathit{q}}{\mathit{r}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}\right]}\mathit{.}$(C.3)Combining Eqs. (C.1)–(C.3), and considering all particles on bound orbits, gives (cf. Krivov et al. 2005) $\begin{array}{ccc}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\sin \mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}}\frac{\mathrm{1}}{{\mathit{r}}^{\mathrm{3}}}\underset{}{{\mathrm{\int }}^{\mathit{q}}}\underset{}{{\mathrm{\int }}^{\mathit{e}}}\mathit{n}\mathrm{\left(}\mathit{m,q,e}\mathrm{\right)}{\left[\frac{\mathit{r}}{\mathit{q}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}\right]}^{\mathrm{3}\mathit{/}\mathrm{2}}\\ & & \mathrm{\times }{\left[\mathrm{2}\mathrm{-}\frac{\mathit{r}}{\mathit{q}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}\mathrm{-}\frac{\mathit{q}}{\mathit{r}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{d}\mathit{q}\hspace{0.17em}\mathrm{d}\mathit{e,}\end{array}$(C.4)with the integration domain $\frac{\mathrm{1}\mathrm{-}\mathit{e}}{\mathrm{1}\mathrm{+}\mathit{e}}\mathit{r}\mathrm{\le }\mathit{q}\mathrm{\le }\mathit{r,} \mathrm{0}\mathrm{\le }\mathit{e}\mathit{<}\mathrm{1.}$(C.5)The contribution of unbound particles to is calculated using their production rate ṅ(m,q,e) and the radial velocity with which they leave the system (Löhne, priv. comm.). Assuming all unbound particles are created at the periastron of their orbits, their vertically-averaged particle number density is given by $\begin{array}{ccc}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi rh}}\underset{}{{\mathrm{\int }}^{\mathit{q}}}\underset{}{{\mathrm{\int }}^{\mathit{e}}}\frac{\mathit{ṅ}\mathrm{\left(}\mathit{m,q,e}\mathrm{\right)}}{\mathrm{\left|}\mathit{ṙ}\mathrm{\right(}\mathit{m,q,e,r}\mathrm{\left)}\mathrm{\right|}}\hspace{0.17em}\mathrm{d}\mathit{q}\hspace{0.17em}\mathrm{d}\mathit{e}\\ & \mathrm{=}& \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }\sin \mathit{\epsilon }\sqrt{\mathit{G}{\mathit{M}}_{\mathit{\star }}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\beta }\mathrm{)}{\mathit{r}}^{\mathrm{3}}}}\underset{}{{\mathrm{\int }}^{\mathit{q}}}\underset{}{{\mathrm{\int }}^{\mathit{e}}}\mathit{ṅ}\mathrm{\left(}\mathit{m,q,e}\mathrm{\right)}\\ & & \mathrm{\times }{\left[\mathrm{2}\mathrm{-}\frac{\mathit{r}}{\mathit{q}}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{e}\mathrm{)}\mathrm{-}\frac{\mathit{q}}{\mathit{r}}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{e}\mathrm{)}\right]}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{d}\mathit{q}\hspace{0.17em}\mathrm{d}\mathit{e,}\end{array}$with the integration domain $\mathit{q}\mathrm{\le }\mathit{r,} \mathit{e}\mathrm{\le }\mathrm{-}\mathrm{1}\mathrm{\vee }\mathit{e}\mathrm{\ge }\mathrm{1.}$(C.8)Negative eccentricities correspond to “anomalous” hyperbolic orbits (see Krivov et al. 2006).

In applying this theory to the raw data, we replace the integrals in Eqs. (C.4) and (C.7) with sums; replace n(m,q,e) and ṅ(m,q,e) with n and ṅ, respectively; and evaluate the resulting equations at discrete points of r. For each bin, we sample the phase space it represents using a Monte Carlo method.

Appendix C.2: Derived quantities

For radial dust distribution profiles, we use the vertical geometrical optical depth τ_geo, defined as the surface density of collective (i.e., combining particles of all sizes) cross-section. It is computed from as $\begin{array}{ccc}{\mathit{\tau }}_{\mathrm{geo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{=}& \mathit{h}\underset{}{{\mathrm{\int }}^{\mathit{m}}}\mathit{\sigma }\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{m}\\ & \mathrm{=}& \mathrm{2}\mathit{\pi }\sin \mathrm{\left(}\mathit{\epsilon }\mathrm{\right)}\mathit{r}\underset{}{{\mathrm{\int }}^{\mathit{m}}}{\mathit{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{m}\mathrm{\right)}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{m}\mathit{.}\end{array}$To characterize size distributions, we use the quantity A(s,r), defined as the cross-section density per base-10 logarithmic unit of size. It is given by $\begin{array}{ccc}\mathit{A}\mathrm{\left(}\mathit{s,r}\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{d}\mathit{s}}{\mathrm{d}{\log }_{\mathrm{10}}\mathrm{\left(}\mathit{s}\mathrm{\right)}}\frac{\mathrm{d}\mathit{m}}{\mathrm{d}\mathit{s}}\mathit{\sigma }\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{4}{\mathit{\pi }}^{\mathrm{2}}\ln \mathrm{\left(}\mathrm{10}\mathrm{\right)}{\mathit{\rho }}_{\mathrm{d}}{\mathit{s}}^{\mathrm{5}}\mathrm{𝒩}\mathrm{\left(}\mathit{m,r}\mathrm{\right)}\mathit{.}\end{array}$This quantity allows for an easy comparison between the relative contributions of particles with different sizes to the total cross-section of the disk. In the size distributions plots (Figs. 5 and 6), a horizontal line means particles of all sizes contribute equally to the cross-section, and equal areas under the curve correspond to equal contributions to the total cross-section.

Acknowledgments

We thank the anonymous referee for a very thorough and detailed report, that helped us to improve the quality of this paper.

References

All Tables

All Figures

Fig. 1 Maximum geometrical optical depth as a function of distance to the star (Eq. (5)). The profiles are derived from the analytical model of Wyatt (2005), in which dust is produced by a source at radius r0 and subsequently migrates inward due to P–R drag, while suffering destruction from mutual collisions. The solid lines correspond to solar-mass stars, the dashed lines to 2 M⊙ stars. Profiles for different parent belt locations are shown in different colors. All profiles assume dust grains with β = 0.5.

In the text

Fig. 2 Characteristic timescale as function of radial distance for sublimation (Eq. (15)), P–R drag (Eq. (2)), and mutual collisions (minimum timescale, Eq. (23)), for β = 0.5 particles in a debris disk around a solar-mass star with a very dense (η0 ≫ 1) parent belt located at 30 AU. The gray vertical lines indicate the radial distances used by the analytical model: rpile is the radius for which the sublimation timescale equals the P–R drag timescale, rcrit is the radius for which the P–R drag timescale equals the collisional timescale, and r0 is the location of the parent belt. For the sublimation timescale, we assume that the dust grains are solid spheres of graphite (see Sect. 2.2.1 for the values of the sublimation parameters) with a radius of s = 0.64 μm (the size corresponding to β = 0.5).

In the text

Fig. 3 Maximum geometrical optical depth profiles of debris disks with collisions, P–R drag, and sublimation (Eq. (28) with τgeo, base(r) given by Eq. (5)), for different values of relative pile-up width Δrpile/rpile shown in different colors. The solid lines correspond to disks around solar-mass stars, the dashed lines to disks around 2 M⊙ stars. All profiles assume a parent belt at r0 = 30 AU, dust grains with β = 0.5, and rpile ≈ 0.019 AU, and rpile ≈ 0.064 AU, for the M⋆ = 1 M⊙, and M⋆ = 2 M⊙ cases, respectively, corresponding to spherical graphite grains with β = 0.5.

In the text

	Fig. 4 Geometrical optical depth profiles of debris disks with a dense parent belt at 30 AU around stars of M⋆ = 1 M⊙ (solid lines) and M⋆ = 2 M⊙ (dashed lines). The black lines show the end results of the numerical simulations. In green are the maximum τgeo profiles as given by the analytical model (Eq. (28), with τgeo, base(r) given by Eq. (5), Δrpile/rpile = 0.15, r0 = 30 AU, and β = 0.5).
In the text

Fig. 5 Size distributions at different radial distances for the M⋆ = 1 M⊙ run before switching on sublimation. The quantity on the vertical axis, A, is the cross-section density per unit size decade (see Appendix C.2), which is horizontal if all sizes contribute equally to the cross-section. The gray band indicates the slope of a size distribution that follows the classical Dohnanyi (1969) power law (n(s) ∝ s-3.5, hence A ∝ s-0.5). The vertical lines mark the particle sizes corresponding to relevant values of β.

In the text

	Fig. 6 Size distribution at r = 0.02 AU (the location of pile-up) for the M⋆ = 1 M⊙ run, before and after sublimation is switched on.
In the text

Fig. 7 SEDs of stars and their debris disks at a distance of 10 pc. The disk SED is shown for different dust distributions with the numerical results in black and the analytical distributions in green. Dust distributions including the pile-up of dust in the sublimation zone are shown with dashed lines, and distributions in which the pile-up is excluded are shown with solid lines, but the spectra largely overlap. Disk SEDs are calculated according to Eqs. (46) and (47). The stellar photosphere is indicated by a dotted Planck curve. The vertical gray areas mark the NIR H and K bands in which hot exozodiacal dust is observed.

In the text

Fig. 8 Flux ratio between the disk and the star at 2.1 μm (K band), as function of parent belt radius r0. Solid lines indicate the flux ratio due to the maximum amount of material moved in by P–R drag, truncated at r = rpile. Ratios including the small effect due to pile-up are shown with dashed lines. Different colors are used for different stellar types. The light shaded area marks the region in parameter space where no significant pile-up is expected to occur because the orbits of small dust grains cannot circularize sufficiently (see Eq. (22)). The dark shaded area marks where the parent belt r0 is closer than the pile-up radius rpile. The observed K-band excess fluxes of main-sequence stars from Table 1 are shown at their respective estimated parent belt distances with a different symbol for each star and coloring according to the star’s spectral type: Vega (the CHARA/FLUOR measurement): circles; β Leo: hexagon; Fomalhaut: squares; β Pic: diamond; η Lep: upward pointing triangles; 110 Her: left-pointing triangle; 10 Tau: right-pointing triangle; τ Cet: pentagon.

In the text

	Fig. B.1 Evolution of the radial geometrical optical depth profile of the benchmark run, to be compared with Fig. 10 of Krivov et al. (2006).
In the text

	Fig. B.2 Evolution of the size distribution at r = 100 AU of the benchmark run, to be compared with Fig. 6 of Krivov et al. (2006).
In the text