Nearinfrared emission from sublimating dust in collisionally active debris disks^{⋆}
^{1}
Astronomical Institute “Anton Pannekoek”, University of
Amsterdam,
PO Box 94249,
1090 GE
Amsterdam,
The Netherlands
email:
r.vanlieshout@uva.nl
^{2}
Department of Astrophysics/IMAPP, Radboud University
Nijmegen, PO Box
9010, 6500 GL
Nijmegen, The
Netherlands
^{3}
Leiden Observatory, Leiden University,
PO Box 9513, 2300 RA
Leiden, The
Netherlands
Received: 14 June 2013
Accepted: 8 April 2014
Context. Hot exozodiacal dust is thought to be responsible for excess nearinfrared (NIR) emission emanating from the innermost parts of some debris disks. The origin of this dust, however, is still a matter of debate.
Aims. We test whether hot exozodiacal dust can be supplied from an exterior parent belt by Poynting–Robertson (P–R) drag, paying special attention to the pileup of dust that occurs owing to the interplay of P–R drag and dust sublimation. Specifically, we investigate whether pileups still occur when collisions are taken into account, and if they can explain the observed NIR excess.
Methods. We computed the steadystate distribution of dust in the inner disk by solving the continuity equation. First, we derived an analytical solution under a number of simplifying assumptions. Second, we developed a numerical debris disk model that for the first time treats the complex interaction of collisions, P–R drag, and sublimation in a selfconsistent way. From the resulting dust distributions, we generated thermal emission spectra and compare these to observed excess NIR fluxes.
Results. We confirm that P–R drag always supplies a small amount of dust to the sublimation zone, but find that a fully consistent treatment yields a maximum amount of dust that is about 7 times lower than that given by analytical estimates. The NIR excess due to this material is much less (≲10^{3} for Atype stars with parent belts at ≳1 AU) than the values derived from interferometric observations (~10^{2}). Pileup of dust still occurs when collisions are considered, but its effect on the NIR flux is insignificant. Finally, the crosssection in the innermost regions is clearly dominated by barely bound grains.
Key words: zodiacal dust / circumstellar matter
Appendices are available in electronic form at http://www.aanda.org
© 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 nearinfrared (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.
NIR interferometric detections of hot exozodiacal dust, together with associated outer debris belt locations.
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)^{1}, 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 pileup effect for dragdominated 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 pileup of dust in the sublimation zone of a collisionally active debris disk and whether such a pileup could explain the exozodiacal NIR emission observed very close to some stars. To do this, we compute the steadystate 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 selfconsistent way (Sect. 3). From the obtained steadystate 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 postprocessing 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 crosssection^{2}. Under the two assumptions listed above, it is given by (1)Here, r is the radial distance from the central star, σ the crosssection of a particle, and n(r) the onedimensional 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 pileup 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_{0} with a geometrical optical depth of τ_{geo}(r_{0}); (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) (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} )^{3}.
The collisional timescale t_{coll} indicates the average time between two collisions for a given particle. Wyatt (2005) finds (3)where t_{orb} is the orbital period of a circular orbit, given by . This equation is valid for particles that are on circular orbits and whose most important collisional partners have similar sizes^{4}.
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 steadystate solution is The parameter η_{0} characterizes the density of the parent belt. It is defined such that for η_{0} = 1 the collisional and P–R drag timescales are equal at r_{0}. Disks with η_{0}> 1 are collisiondominated at r_{0}, while disks with η_{0}< 1 are dragdominated at r_{0}. Most debris disks with observed outer belts have η_{0}> 10 (Wyatt 2005).
The balance between P–R drag and collisions is selflimiting: 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 η_{0} ≫ 1 (i.e., a very dense parent belt) of (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^{5}.
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 r_{0} and subsequently migrates inward due to P–R drag, while suffering destruction from mutual collisions. The solid lines correspond to solarmass 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. 

Open with DEXTER 
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 pileup, 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 (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) (7)The grain parameter β can be expressed as (e.g., Burns et al. 1979) (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 (collisionlimited) P–R drag mass flux (cf. Rafikov 2011) (9)The maximum mass flux can now be found by substituting Eq. (5) for τ_{geo}(r) in Eq. (9), which yields (10)Figure 1 shows that max [ τ_{geo}(r) ] levels off in the innermost part of the system (r ≪ r_{0}), at the radii where exozodiacal dust is seen. The mass flux corresponding to this plateau is (11)The maximum mass flux only depends on grain properties through Q_{pr} and β^{6}. 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 nondestructive).
2.2. P–R drag and sublimation
In the preceding, we found that P–R drag supplies a small but nonzero 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 gasfree environment, the rate at which the grain radius s changes is given by (e.g., Kobayashi et al. 2008) (12)Here, P_{v} is the phaseequilibrium 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) (13)where P_{0} is a normalization constant and H the latent heat of sublimation. By assuming P_{0} 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_{3} sublimation are P_{0} = 2.95 × 10^{14} dyn cm^{2} and H = 2.15 × 10^{11} 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 blackbody temperature (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 blackbody 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 (15)This estimate assumes that the grain temperature remains constant (at T_{bb}) throughout the sublimation process. As a sublimating particle becomes smaller, the blackbody approximation is bound to become inaccurate. For sufficiently large particles, however, the true sublimation time is dominated by the blackbody regime, and Eq. (15) provides a good approximation.
2.2.2. The pileup 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 pileup 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 pileup 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 λ_{⋆}^{7}, 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 semimajor 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 semimajor axis due to P–R drag is compensated by the increase of semimajor 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)^{8}, 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 pileup has the blackbody temperature T_{pile} = T_{bb}(r_{pile}), which holds for s ≳ λ_{⋆}, this equation can be written as (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 (17)The pileup 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 blackbody approximation is valid, so the grain temperature in the pileup is independent of grain size.
2.2.3. Conditions for dust pileup
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 pileup 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 crosssection 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 pileup. 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 pileup, low luminosity stars do not have significant pileups. Kobayashi et al. (2009, 2011) set the limit at β_{max}> 0.5, and derive a lower limit on the stellar luminosity for significant pileup, assuming that Q_{pr} ≈ 1 holds for s ≳ λ_{⋆}. This limit is (18)where T_{⋆} is the effective temperature of the central star.
The pileup 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 pileup 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 (19)The particles that are the most important for the pileup 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 pileup 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) (20)where (a_{0},e_{0}) are the initial semimajor axis and eccentricity, which evolve into (a_{1},e_{1}). This coupling can be used to place a lower limit on the distance of the source region, if significant pileup is to occur, using the maximum allowed eccentricity in the sublimation zone for efficient pileup. Written in terms of periastron distance q = a(1 − e), Eq. (20) becomes (21)Substituting r_{0} and r_{pile} for q_{0} and q_{1}, respectively, and using e_{0} ≈ 1 and e_{1} ≲ 10^{2}, yields a lower limit on the source radius (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_{0}. An example is the drag force from small amounts of gas that are present in the disk (Takeuchi & Artymowicz 2001).
2.3. Pileup 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 solarmass 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: r_{pile} is the radius for which the sublimation timescale equals the P–R drag timescale, r_{crit} is the radius for which the P–R drag timescale equals the collisional timescale, and r_{0} 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). 

Open with DEXTER 
2.3.1. Do collisions interfere with dust pileup?
Since collisions might interfere with the process of dust pileup, 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 (η_{0} ≫ 1), τ_{geo}(r) is given by Eq. (5). Putting this in Eq. (3) results in the minimum collisional timescale (23)In Fig. 2 this timescale is compared with the sublimation and P–R drag timescales for a debris disk around a solarmass 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 pileup. 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 (24)For η_{0} = 1 we recover r_{crit} = r_{0}, as expected. From the limiting case η_{0} ≫ 1 (i.e., a very dense parent belt), we find (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 pileup to occur, P–R drag dominates collisions in the sublimation zone, and collisions are so infrequent there that they do not interfere with the pileup process.
The pileup of dust means that τ_{geo} increases around r_{pile}, locally decreasing the collisional timescale. However, a significant pileup 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 . To overcome this difference, the pileup 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 pileup magnitude
The efficiency of dust pileup 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 orderofmagnitude estimate of the maximum amount of material in the pileup, which can be used to assess whether pileups can explain the observed NIR excess emission.
Dust grains reside in the pileup 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 pileup occurs roughly at r_{pile}, defined such that t_{subl}(r_{pile}) = t_{PR}(r_{pile}), the dust stays in the pileup for about one P–R drag timescale. Given this, the total number of particles in the pileup is (26)To describe the radial profile in terms of geometrical optical depth, it is necessary to specify the radial width of the pileup Δr_{pile}. In reality, Δr_{pile} depends on the orbital eccentricities of the particles in the pileup and on differences in pileup distance for particles of different sizes that contribute. Since this is beyond the scope of this work, we keep the relative pileup 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 pileup is given by (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 pileup 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) (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 pileup width Δr_{pile}/r_{pile} shown in different colors. The solid lines correspond to disks around solarmass stars, the dashed lines to disks around 2 M_{⊙} stars. All profiles assume a parent belt at r_{0} = 30 AU, dust grains with β = 0.5, and r_{pile} ≈ 0.019 AU, and r_{pile} ≈ 0.064 AU, for the M_{⋆} = 1 M_{⊙}, and M_{⋆} = 2 M_{⊙} cases, respectively, corresponding to spherical graphite grains with β = 0.5. 

Open with DEXTER 
We can make a rough estimate of the geometrical optical depth enhancement factor f_{τgeo} of the pileup (i.e., how much higher τ_{geo} in the pileup is compared to the base level of the inner disk), as a function of the radial width of the pileup Δr_{pile}: (29)For a pileup 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 Ftype stars.
For comparison with observations, it suffices to assume that all the pileup material is located at r_{pile} (i.e., Δr_{pile}/r_{pile} ≪ 1). This gives the highest possible temperature to all pileup 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 pileup. 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 blackbody grains), it can be approximated by the fraction of the star that is covered by dust (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 pileup is (31)In the limit of r_{0} ≫ r_{pile}, this becomes (32)This is only the fractional luminosity due to the dust in the pileup. 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 nonzero inward mass flux ofdust in the inner disk, which is selflimited by collisions(Eq. (11)).

2.
A pileup 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 pileup.

4.
Given that the pileup of dust occurs around the radial distance where the sublimation timescale equals the P–R drag timescale and that sublimating dust resides in the pileup 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 selfconsistently 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 pileup 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 selfconsistently 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 timedependent 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 phasespace coordinates in a continuous fashion (P–R drag and sublimation) as diffusion terms, and processes that abruptly change phasespace coordinates (collisions) as source and sink terms. Formally, this is decribed by (cf. Krivov et al. 2005; Löhne 2008) (33)where n(m,k,t) is the phasespace 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 phasespace 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 phasespace variables are (1) particle mass m; (2) orbital eccentricity e; and (3) an orbital element characterizing the size of the orbit, such as semimajor axis a or periastron distance q = a(1 − e).
The final phasespace variable can be chosen to fit the numerical needs of the problem under investigation. For our study of the pileup 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 semimajor axis. Using q instead of a as phasespace 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 phasespace grid. Since the P–R drag timescale is usually longer than the orbital period, we use the orbitaveraged change rates of the orbital elements, given by (e.g., Burns et al. 1979) For the rate of change in periastron distance, we find
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 (38)In our numerical model, we use realistic dust grain temperatures (as opposed to the blackbody 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 orbitaveraged change rate in orbital elements and mass to be with where f denotes the true anomaly. For the periastron distance, we find For each phasespace bin, quantities and are numerically evaluated using the standard Euler method^{9}. The change rates of the phasespace variables (Eqs. (40), (41), and (45)) are then used in diffusion terms in the continuity equation.
Using orbitaveraged mass loss rates is only correct if the sublimation timescale is longer than the orbital period. For phasespace 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 phasespace coordinates. In the continuity equation, they are described by sink terms at the phasespace 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 phasespace 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 crosssection, 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 semiopening 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 edgeon 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^{10}.
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 , 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 , 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 , 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 . 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 solarmass star and one using a 2 M_{⊙} star. Following the mass–luminosity relation for mainsequence 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_{0} = 30 AU. Lower values of r_{0} 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_{0}) ≈ 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^{11}. This roughly corresponds to η_{0} ~ 10 for both stellar mass cases, which is apparently enough to approximate an inner disk τ_{geo} profile corresponding to η_{0} ≫ 1. For the semiopening 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 , 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, 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 erg g^{1} found in laboratory experiments with highvelocity 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 phasespace grid
Because this problem is computationally very demanding, great care has to be taken in setting up the phasespace grid. Specifically, resolving the pileup 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 nonuniform 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 highresolution, 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 lowresolution, 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 highresolution part consists of 30 bins between s = 0.5 μm and s = 10 μm. The highresolution 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_{0}. This artificial source term adds particles with sizes between s = 1 mm and s = 1 cm, following the powerlaw 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 steadystate 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 (base10) time steps become less than 1%^{12}. 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 pileup (i.e., the dust in the pileup can be isolated by subtracting the presublimation state from the final one).
3.3. Results
The output of each model run is a steadystate 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 crosssection 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 pileup 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 pileup) 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 crosssection 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), Δr_{pile}/r_{pile} = 0.15, r_{0} = 30 AU, and β = 0.5). 

Open with DEXTER 
As predicted, switching on sublimation leads to an accumulation of dust. The pileups are located very close to r_{pile}, as determined for graphite grains, using blackbody 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 pileup is between 2000 and 2100 K. In both runs, the pileup has a geometrical optical depth enhancement factor of about f_{τgeo} ≈ 3. The fractional luminosities due to the pileups 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 pileup 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 dragdominated 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 collisiondominated. This indicates that the density of the parent belt (characterized by η_{0} ~ 10) is insufficient to make the outer disk completely collisiondominated.
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 crosssection density per unit size decade (see Appendix C.2), which is horizontal if all sizes contribute equally to the crosssection. 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 β. 

Open with DEXTER 
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 selfsimilar collisions), from the blowout radius upwards. Particles with β> 0.5 are depleted by about three orders of magnitude in terms of collective crosssection. Superimposed on the power law is a wellknown 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 powerlaw 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 crosssection, 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 crosssection, 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 pileup). Figure 6 shows the size distribution in the pileup, 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 pileup 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 pileup) for the M_{⋆} = 1 M_{⊙} run, before and after sublimation is switched on. 

Open with DEXTER 
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 pileup of dust in the sublimation zone are shown with dashed lines, and distributions in which the pileup 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. 

Open with DEXTER 
4. Comparison with observations
To assess whether the pileup 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 (viewingangledependent) 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 blackbody temperature (Eq. (14)). The flux density of the dust, expressed as a function of wavelength λ, can then be computed as (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 (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_{0} = 30 AU, β = 0.5, and all the pileup material at r_{pile} (i.e., Δr_{pile}/r_{pile} ≪ 1)^{13}. Singularities at r = r_{0} are removed by imposing the condition τ_{geo} ≤ 0.01, which corresponds to η_{0} ≈ 2200. The numerical dust distributions without pileup were created by subtracting the isolated pileup dust from the final profile (see Sect. 3.2.4). The stellar spectra are given by blackbody 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 mainsequence 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 solarmass 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 blackbody temperature.
Only considering thermal emission, the NIR flux originates almost exclusively in the inner 1 AU of the debris disk. The pileup 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 midIR flux at a comparable level, which is incompatible with observed excess spectra (e.g., Akeson et al. 2009). The pileup 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 r_{0}. Solid lines indicate the flux ratio due to the maximum amount of material moved in by P–R drag, truncated at r = r_{pile}. Ratios including the small effect due to pileup 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 pileup 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 r_{0} is closer than the pileup radius r_{pile}. The observed Kband excess fluxes of mainsequence 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: leftpointing triangle; 10 Tau: rightpointing triangle; τ Cet: pentagon. 

Open with DEXTER 
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 powerlaw 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 crosssection 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^{14}, 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 crosssection 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 steadystate 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 pileup 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 (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, blackbody graphite particles, we find κ ≈ 10. This shows that particles starting with β_{1} ≈ 0.5, and e_{1} ≈ 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) 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 pileup scenario could work for Vega if the disk is dragdominated, which requires γ ≈ 300. While this value is consistent with the upper limit on the massloss rate of Vega from radiocontinuum observations (Ṁ_{⋆} ≲ 10^{10} M_{⊙} yr^{1}; Hollis et al. 1985), stellar wind models predict much lower massloss rates for mainsequence Atype stars (Ṁ_{⋆} ≲ 10^{16} M_{⊙} yr^{1}; Babel 1995). This theoretical massloss 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 mainsequence Atype 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 pileup of dust is too inefficient to have an effect on the slope of the SED^{15}. 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 crosssection of the material. To investigate this scenario, the steadystate 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 lowmass 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 pileup 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 always 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 depends on the mass of the central star and distance to the parent belt(Fig. 1). When the production of dust is treatedselfconsistently, 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 Atype stars with parent belts at ≳1 AU (Fig. 8).

3.
The pileup 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 pileup 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 pileup. The fractional luminosity provided by dust in the pileup is relatively low, so the pileup 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 pileup, 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.
Kennedy & Wyatt (2013) find that “warm” exozodiacal dust around solartype stars can be explained by insitu planetesimal belts. This type of exozodiacal dust is detected at midinfrared wavelengths and has a typical temperature of a few hundred K, placing it around 1 AU from the star.
Assuming circular orbits is valid for particles with low β ratios, while the small particles most relevant to our study have high β ratios. We relax this simplifying assumption in our numerical model (Sect. 3).
Equation (3) ignores a factor in the orbital period of radiation pressure affected particles, but this only changes the collisional timescale by a factor of about 0.7 for β = 0.5.
Debris disks around stars with a strong stellar wind have higher values of τ_{geo}. We ignore stellar wind in our present analysis, and discuss its effects in Sect. 5.2.
More precise estimates of the pileup distance are given by Kobayashi et al. (2009, 2011). We use this simple approximation to facilitate the estimate of the pileup magnitude in Sect. 2.3.2.
If the radial geometrical optical depth is higher than ~10^{2}, the disk is subject to dust avalanches (Grigorieva et al. 2007).
The actual input parameter used (which indirectly determines the geometrical optical depth at the source region) is the mass supply term of large dust particles in the source region (see Sect. 3.2.4). It is set to 10^{10}M_{⊕}yr^{1} for the M_{⋆} = 1 M_{⊙} run and 8 × 10^{10}M_{⊕}yr^{1} for the M_{⋆} = 2 M_{⊙} run. Comparing these mass fluxes to those given by Eq. (11) indicates that the vast majority of the material is destroyed in collisions before it reaches the sublimation zone.
Formally, steady state is only reached after ~ 10 Gyr, which is the time it takes for the largest particles we consider to move from the parent belt to the sublimation zone by P–R drag. The barely bound grains that dominate the crosssection, however, already settle into a steady state after ~ 10 Myr, which is short compared to the typical lifetime of a debris disk.
Our approach to computing the SED of a dust profile with pileup is different from that of Kobayashi et al. (2011), who use τ_{geo, pile} and Δr_{pile} as independent input variables. In our estimate, the total amount of dust is fixed, and placing it all at r_{pile} is the most optimistic configuration for detecting the pileup.
In debris disks relative velocities are mostly due to the eccentricities and mutual inclinations of particle orbits. While orbital eccentricities are diminished in the inner disk by P–R drag induced circularization, inclinations are simply inherited from the parent belt. This is handled correctly in our code, since orbital inclinations are parameterized by the (constant) opening angle of the disk.
This corresponds to Eq. (3.107) of Löhne (2008), corrected for a typographical error.
Acknowledgments
We thank the anonymous referee for a very thorough and detailed report, that helped us to improve the quality of this paper.
References
 Abrahamson, J. 1974, Carbon, 12, 111 [CrossRef] [Google Scholar]
 Absil, O., di Folco, E., Mérand, A., et al. 2006, A&A, 452, 237 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Absil, O., di Folco, E., Mérand, A., et al. 2008, A&A, 487, 1041 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Absil, O., Mennesson, B., Le Bouquin, J.B., et al. 2009, ApJ, 704, 150 [NASA ADS] [CrossRef] [Google Scholar]
 Absil, O., Defrère, D., Coudé du Foresto, V., et al. 2013, A&A, 555, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Acke, B., Min, M., Dominik, C., et al. 2012, A&A, 540, A125 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Akeson, R. L., Ciardi, D. R., MillanGabet, R., et al. 2009, ApJ, 691, 1896 [NASA ADS] [CrossRef] [Google Scholar]
 Allen, C. W. 1976, Astrophysical Quantities (London: Athlone) [Google Scholar]
 Babel, J. 1995, A&A, 301, 823 [Google Scholar]
 Backman, D. E., & Paresce, F. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine, 1253 [Google Scholar]
 Benz, W., & Asphaug, E. 1999, Icarus, 142, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Bonsor, A., Augereau, J.C., & Thébault, P. 2012, A&A, 548, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bonsor, A., Kennedy, G. M., Crepp, J. R., et al. 2013, MNRAS, 431, 3025 [NASA ADS] [CrossRef] [Google Scholar]
 Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Campo Bagatin, A., Cellino, A., Davis, D. R., Farinella, P., & Paolicchi, P. 1994, Planet. Space Sci., 42, 1079 [NASA ADS] [CrossRef] [Google Scholar]
 Carpenter, J. M., Bouwman, J., Mamajek, E. E., et al. 2009, ApJS, 181, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C. H., Patten, B. M., Werner, M. W., et al. 2005, ApJ, 634, 1372 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, C. H., Sargent, B. A., Bohac, C., et al. 2006, ApJS, 166, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Ciardi, D. R., van Belle, G. T., Akeson, R. L., et al. 2001, ApJ, 559, 1147 [NASA ADS] [CrossRef] [Google Scholar]
 Defrère, D., Absil, O., Augereau, J.C., et al. 2011, A&A, 534, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Defrère, D., Lebreton, J., Le Bouquin, J.B., et al. 2012, A&A, 546, L9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dent, W. R. F., Walker, H. J., Holland, W. S., & Greaves, J. S. 2000, MNRAS, 314, 702 [NASA ADS] [CrossRef] [Google Scholar]
 di Folco, E., Thévenin, F., Kervella, P., et al. 2004, A&A, 426, 601 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 di Folco, E., Absil, O., Augereau, J.C., et al. 2007, A&A, 475, 243 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531 [NASA ADS] [CrossRef] [Google Scholar]
 Dominik, C., & Decin, G. 2003, ApJ, 598, 626 [NASA ADS] [CrossRef] [Google Scholar]
 Durda, D. D., & Dermott, S. F. 1997, Icarus, 130, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Eiroa, C., Marshall, J. P., Mora, A., et al. 2013, A&A, 555, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flynn, G. J., & Durda, D. D. 2004, Planet. Space Sci., 52, 1129 [NASA ADS] [CrossRef] [Google Scholar]
 Fujiwara, A., Kamimoto, G., & Tsukamoto, A. 1977, Icarus, 31, 277 [NASA ADS] [CrossRef] [Google Scholar]
 Gáspár, A., Psaltis, D., Rieke, G. H., & Özel, F. 2012, ApJ, 754, 74 [NASA ADS] [CrossRef] [Google Scholar]
 Gáspár, A., Rieke, G. H., & Balog, Z. 2013, ApJ, 768, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Greaves, J. S., Wyatt, M. C., Holland, W. S., & Dent, W. R. F. 2004, MNRAS, 351, L54 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Grigorieva, A., Artymowicz, P., & Thébault, P. 2007, A&A, 461, 537 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Haug, U. 1958, Z. Astrophys., 44, 71 [NASA ADS] [Google Scholar]
 Heng, K., & Tremaine, S. 2010, MNRAS, 401, 867 [NASA ADS] [CrossRef] [Google Scholar]
 Hochbruck, M., & Ostermann, A. 2010, Acta Numerica, 19, 209 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hollis, J. M., Chin, G., & Brown, R. L. 1985, ApJ, 294, 646 [NASA ADS] [CrossRef] [Google Scholar]
 Kalas, P., Graham, J. R., & Clampin, M. 2005, Nature, 435, 1067 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kennedy, G. M., & Wyatt, M. C. 2013, MNRAS, 433, 2334 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, H., Watanabe, S.I., Kimura, H., & Yamamoto, T. 2008, Icarus, 195, 871 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, H., Watanabe, S.I., Kimura, H., & Yamamoto, T. 2009, Icarus, 201, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Kobayashi, H., Kimura, H., Watanabe, S.i., Yamamoto, T., & Müller, S. 2011, Earth, Planets, and Space, 63, 1067 [Google Scholar]
 Krivov, A. V., Sremčević, M., & Spahn, F. 2005, Icarus, 174, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Krivov, A. V., Löhne, T., & Sremčević, M. 2006, A&A, 455, 509 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lagage, P. O., & Pantin, E. 1994, Nature, 369, 628 [NASA ADS] [CrossRef] [Google Scholar]
 Lawler, S. M., Beichman, C. A., Bryden, G., et al. 2009, ApJ, 705, 89 [NASA ADS] [CrossRef] [Google Scholar]
 Lebreton, J., van Lieshout, R., Augereau, J.C., et al. 2013, A&A, 555, A146 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Löhne. 2008, Ph.D. Thesis, FriedrichSchillerUniversität Jena, Germany [Google Scholar]
 Marshall, J. P., Krivov, A. V., del Burgo, C., et al. 2013, A&A, 557, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mennesson, B., Absil, O., Lebreton, J., et al. 2013, ApJ, 763, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Min, M., Hovenier, J. W., & de Koter, A. 2005, A&A, 432, 909 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minato, T., Köhler, M., Kimura, H., Mann, I., & Yamamoto, T. 2004, A&A, 424, L13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Minato, T., Köhler, M., Kimura, H., Mann, I., & Yamamoto, T. 2006, A&A, 452, 701 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Morales, F. Y., Werner, M. W., Bryden, G., et al. 2009, ApJ, 699, 1067 [NASA ADS] [CrossRef] [Google Scholar]
 Pantin, E., Lagage, P. O., & Artymowicz, P. 1997, A&A, 327, 1123 [NASA ADS] [Google Scholar]
 Plavchan, P., Werner, M. W., Chen, C. H., et al. 2009, ApJ, 698, 1068 [NASA ADS] [CrossRef] [Google Scholar]
 Preibisch, T., Ossenkopf, V., Yorke, H. W., & Henning, T. 1993, A&A, 279, 577 [NASA ADS] [Google Scholar]
 Rafikov, R. R. 2011, ApJ, 732, L3 [NASA ADS] [CrossRef] [Google Scholar]
 Reidemeister, M., Krivov, A. V., Stark, C. C., et al. 2011, A&A, 527, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Strubbe, L. E., & Chiang, E. I. 2006, ApJ, 648, 652 [NASA ADS] [CrossRef] [Google Scholar]
 Su, K. Y. L., Rieke, G. H., Malhotra, R., et al. 2013, ApJ, 763, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Takeuchi, T., & Artymowicz, P. 2001, ApJ, 557, 990 [NASA ADS] [CrossRef] [Google Scholar]
 Thébault, P., & Augereau, J.C. 2007, A&A, 472, 169 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Trilling, D. E., Bryden, G., Beichman, C. A., et al. 2008, ApJ, 674, 1086 [NASA ADS] [CrossRef] [Google Scholar]
 Vitense, C., Krivov, A. V., & Löhne, T. 2010, A&A, 520, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyatt, M. C. 2005, A&A, 433, 1007 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyatt, M. C. 2008, ARA&A, 46, 339 [NASA ADS] [CrossRef] [Google Scholar]
 Wyatt, M. C., Smith, R., Greaves, J. S., et al. 2007, ApJ, 658, 569 [NASA ADS] [CrossRef] [Google Scholar]
 Wyatt, S. P., & Whipple, F. L. 1950, ApJ, 111, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Zavitsanos, P. D., & Carlson, G. A. 1973, J. Chem. Phys., 59, 2966 [NASA ADS] [CrossRef] [Google Scholar]
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 phasespace 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 timestep sizes, which reduces the number of steps needed per run (Appendix A.1). Furthermore, we calculated as many timeindependent 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 (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^{16}. 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 timestep 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 (A.2)where A^{′} is a diagonal matrix that contains only the diagonal elements of A, while the terms marked “const” contain the offdiagonal 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(A.3)where Δt denotes the timestep size^{17}. This scheme is known as the exponential Euler method (see Hochbruck & Ostermann 2010 for an introductory review). It is suitable for time integration of semilinear problems, which consist of a stiff linear part and a nonstiff nonlinear part. In short, the strategy is to solve the linear part exactly and to approximate the nonlinear part using an explicit integration scheme.
The timestep 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) (A.4)Here, t and p are the target and projectile bin indices, respectively, and B_{jtp} is an entry in the timeindependent 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 phasespace 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 phasespace 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 nonzero 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 nonzero 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 nonzero 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 phasespace 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 semimajor axes as the distance dimension of the phasespace grid, the benchmark runs were performed with a version of our code that also uses the semimajor 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 phasespace 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). 

Open with DEXTER 
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). 

Open with DEXTER 
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: Postprocessing of model output
The output of a model run are the population levels n of all bins in the threedimensional 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 phasespace 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 (C.1)where P is the particle’s orbital period, and h = 2r sin ε is the disk height, where ε is the semiopening 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 (C.2)The radial velocity ṙ of the particle is given by (C.3)Combining Eqs. (C.1)–(C.3), and considering all particles on bound orbits, gives (cf. Krivov et al. 2005) (C.4)with the integration domain (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 verticallyaveraged particle number density is given by with the integration domain (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) crosssection. It is computed from as To characterize size distributions, we use the quantity A(s,r), defined as the crosssection density per base10 logarithmic unit of size. It is given by This quantity allows for an easy comparison between the relative contributions of particles with different sizes to the total crosssection of the disk. In the size distributions plots (Figs. 5 and 6), a horizontal line means particles of all sizes contribute equally to the crosssection, and equal areas under the curve correspond to equal contributions to the total crosssection.
All Tables
NIR interferometric detections of hot exozodiacal dust, together with associated outer debris belt locations.
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 r_{0} and subsequently migrates inward due to P–R drag, while suffering destruction from mutual collisions. The solid lines correspond to solarmass 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. 

Open with DEXTER  
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 solarmass 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: r_{pile} is the radius for which the sublimation timescale equals the P–R drag timescale, r_{crit} is the radius for which the P–R drag timescale equals the collisional timescale, and r_{0} 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). 

Open with DEXTER  
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 pileup width Δr_{pile}/r_{pile} shown in different colors. The solid lines correspond to disks around solarmass stars, the dashed lines to disks around 2 M_{⊙} stars. All profiles assume a parent belt at r_{0} = 30 AU, dust grains with β = 0.5, and r_{pile} ≈ 0.019 AU, and r_{pile} ≈ 0.064 AU, for the M_{⋆} = 1 M_{⊙}, and M_{⋆} = 2 M_{⊙} cases, respectively, corresponding to spherical graphite grains with β = 0.5. 

Open with DEXTER  
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), Δr_{pile}/r_{pile} = 0.15, r_{0} = 30 AU, and β = 0.5). 

Open with DEXTER  
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 crosssection density per unit size decade (see Appendix C.2), which is horizontal if all sizes contribute equally to the crosssection. 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 β. 

Open with DEXTER  
In the text 
Fig. 6 Size distribution at r = 0.02 AU (the location of pileup) for the M_{⋆} = 1 M_{⊙} run, before and after sublimation is switched on. 

Open with DEXTER  
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 pileup of dust in the sublimation zone are shown with dashed lines, and distributions in which the pileup 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. 

Open with DEXTER  
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 r_{0}. Solid lines indicate the flux ratio due to the maximum amount of material moved in by P–R drag, truncated at r = r_{pile}. Ratios including the small effect due to pileup 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 pileup 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 r_{0} is closer than the pileup radius r_{pile}. The observed Kband excess fluxes of mainsequence 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: leftpointing triangle; 10 Tau: rightpointing triangle; τ Cet: pentagon. 

Open with DEXTER  
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). 

Open with DEXTER  
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). 

Open with DEXTER  
In the text 