ROCKY PLANET FORMATION: QUICK AND NEAT

To develop predictions for the magnitude of the IR excess produced during rocky planet formation, we identify likely epochs of dust production. In all of the scenarios discussed above, the growth of pebbles and larger planetesimals with r ≈ 1 m to 1000 km takes place in an opaque gaseous environment where debris mixes with preexisting small particles (see Dauphas & Chaussidon 2011, and references therein). Because we cannot distinguish "primordial" dust from debris, we ignore this phase of dust production.

As solids grow from planetesimals into protoplanets and then planets, there are three modes of dust formation. Mergers of planetesimals into protoplanets typically produce modest amounts of debris (e.g., Wetherill & Stewart 1993; Weidenschilling et al. 1997; Kenyon & Luu 1999). As protoplanets grow, they stir the orbits of leftover planetesimals that never become incorporated into a planet. Among the planetesimals, high-velocity collisions then begin to produce numerous smaller particles with sizes ranging from 1 μm to tens of kilometers (Greenberg et al. 1978; Wetherill & Stewart 1993; Kenyon & Bromley 2004; Weidenschilling 2010; Raymond et al. 2011). This debris fuels a collisional cascade, where solids are ground down into submicron-sized particles that are ejected by radiation pressure. Eventually, giant impacts produce larger and larger protoplanets. Debris from giant impacts adds to the debris from collisions of leftover planetesimals and fuels the collisional cascade (Jackson & Wyatt 2012; Genda et al. 2015b).

In all three scenarios, the growth of protoplanets, collisions of leftover planetesimals, and giant impacts generate debris. For simplicity, we ignore the modest amount of debris generated from the growth of protoplanets and derive separate estimates for dust produced from destructive collisions of leftover planetesimals (Section 4.2) and from giant impacts (Section 4.3).

In classical planet formation theory, semianalytical and numerical calculations of collisions among ensembles of rocky planetesimals within 2–3 au report debris production ranging from roughly 5% to almost 50% of the initial mass in solid material (Greenberg et al. 1978; Wetherill & Stewart 1993; Kenyon & Bromley 2004, 2005; Leinhardt & Richardson 2005; Chambers 2008; Weidenschilling 2010; Kenyon & Bromley 2016). The typical mass in the debris is 10%–20% of the initial mass. Significant numbers of destructive collisions begin early, at ∼0.01–0.1 Myr, and last until 10–100 Myr (see also Morbidelli et al. 2012; Raymond et al. 2014; Quintana et al. 2016, and references therein).

In these calculations, the timing of debris production overlaps epochs when we expect a significant decay in the surface density Σ_g of the gaseous disk (Hartmann et al. 1998; Haisch et al. 2001; Mamajek 2009; Williams & Cieza 2011; Alexander et al. 2014, pp. 475–96). When Σ_g is large, gas drag forces debris to spiral into the central star. As Σ_g declines, the system retains a larger and larger fraction of the debris. For young stars with no gaseous disk at ages of 10–20 Myr, we conservatively estimate that 5%–10% of the initial mass in solids is converted into debris. The formation of an Earth-mass planet thus generates 0.05–0.10 ${M}_{\oplus }$ in debris, which is 10–20 times larger than the mass required to produce a detectable 24 μm excess.

To predict debris production from giant impacts, we first consider a simple analytical model. In this approach, we compare the collision energy to the binding energy of a pair of protoplanets (see also Davis et al. 1985; Davis & Ryan 1990; Housen & Holsapple 1990, 1999). We assume conservatively that all collisions are head-on (impact parameter b = 0); more glancing collisions typically yield more debris (e.g., Leinhardt & Stewart 2012). For two protoplanets with mass m and radius r, the center-of-mass collision energy is ${Q}_{c}={v}_{\mathrm{imp}}^{2}/8$ (e.g., Wetherill & Stewart 1993; Kenyon et al. 2014; Kenyon & Bromley 2016), where the impact velocity⁵ is

$$\begin{eqnarray}&&{v}_{\mathrm{imp}}^{2}\approx {v}_{\mathrm{rel}}^{2}+{v}_{\mathrm{rel}}^{2}+{v}_{\mathrm{esc}}^{2}.\end{eqnarray} \tag{ 4 }$$

Here v_rel is the velocity of each protoplanet relative to a circular orbit and v_esc is the mutual escape velocity of the pair at the moment of the collision, ${v}_{\mathrm{esc}}=\sqrt{2G(m+m)/(r+r)}\,=\sqrt{2{Gm}/r}$. The fraction of material ejected during a collision is

$$\begin{eqnarray}&&{f}_{\mathrm{ej}}=0.5\displaystyle \frac{{Q}_{c}}{{Q}_{D}^{\star }},\end{eqnarray} \tag{ 5 }$$

where ${Q}_{D}^{\star }$ is the binding energy (e.g., Agnor et al. 1999; Benz & Asphaug 1999; Canup & Asphaug 2001; Leinhardt & Richardson 2005; Leinhardt & Stewart 2009; Genda et al. 2012, 2015a). When Q_c ≈ ${Q}_{D}^{\star }$, half the mass of the merged pair is ejected to infinity.

Protoplanets typically have relative velocities of 25%–75% of the escape velocity of the largest protoplanets. We set ${v}_{\mathrm{imp}}\approx {{ev}}_{K}$, where v_K is the velocity of a circular orbit and e is the eccentricity. When protoplanets have masses ranging from a lunar mass to an Earth mass, e ≈ 0.1–0.2 (e.g., Chambers 2001; Raymond et al. 2004, 2005; Kenyon & Bromley 2006; Kokubo et al. 2006; Kokubo & Genda 2010; Chambers 2013; Quintana et al. 2016).

For parameters appropriate for rocky objects at 1 au (see Appendix A.1), collisions between equal-mass protoplanets with e ≈ 0.05–0.2 eject roughly 10% of the total mass (Figure 2). Setting the relative velocity equal to 50% of the escape velocity of the protoplanet yields similar results. Although the accretion history of an Earth-mass object is complicated (e.g., Kenyon & Bromley 2006; Kokubo et al. 2006; Chambers 2013), final assembly requires several collisions of sub-Earth-mass objects. If each collision loses roughly 10% of the initial mass, the total amount of lost mass exceeds 0.1 ${M}_{\oplus }$. This mass is roughly an order of magnitude larger than the mass required to produce a detectable debris disk (Section 2).

Zoom In Zoom Out Reset image size

Detailed n-body and SPH calculations support this simple estimate (e.g., Agnor & Asphaug 2004; Asphaug et al. 2006; Raymond et al. 2011; Genda et al. 2012; Jackson & Wyatt 2012; Stewart & Leinhardt 2012; Chambers 2013; Genda et al. 2015b). In SPH simulations, dust production depends on b and v_rel. Head-on collisions with small b and v_rel yield little or no dust. When v_rel is large or the collision is oblique (b ≳ 0.3–0.4), dust production is substantial. Averaged over a complete n-body simulation of an ensemble of growing protoplanets, protoplanet collisions disperse 10%–20% of the initial mass into small fragments. Thus, giant impacts involved in the formation of a single Earth-mass planet yield at least 0.1–0.2 ${M}_{\oplus }$ in debris.

To derive predicted detection rates for warm debris disks, we adopt an initial disk mass M₀ and a model for the time evolution of the total mass ${M}_{d}$ and cross-sectional area ${A}_{d}$ of the debris produced through collisions of leftover planetesimals and giant impacts. In a real debris disk, stochastic collisions add and remove debris; thus, ${M}_{d}$ and ${A}_{d}$ decrease over the long term and increase and decrease on short timescales (e.g., Grogan et al. 2001; Kenyon & Bromley 2004; Weidenschilling 2010; Jackson & Wyatt 2012; Genda et al. 2015b; Kenyon & Bromley 2016). Instead of following this evolution in detail, we consider an analytical model where ${A}_{d}$ and ${M}_{d}$ decline monotonically with time (see Appendix A.1). With this conservative assumption, we derive a lower limit for the expected IR excess from the debris at any time t. If the analytical model predicts much larger IR excesses than observed, then more extensive numerical calculations of planet formation will also yield much larger excesses than observed.

As outlined in Appendix A.1, the long-term evolution of the debris in the analytical model depends on the surface density of solids, the size of the largest object, and the orbital eccentricity of debris particles. The model assumes that collisions among objects with radii smaller than ${r}_{\max }$ are completely destructive, producing debris with particle sizes smaller than ${r}_{\max }$. The cascade of collisions maintains a power-law size distribution with N(r) ∝ r^−3.5 from the smallest size ${r}_{\min }$ up to ${r}_{\max }$. Radiation pressure ejects smaller particles. Protoplanets with radii larger than ${r}_{\max }$ are ignored. With destructive collisions for all particles smaller than ${r}_{\max }$ and ejection of material with r ≲ ${r}_{\min }$, the mass in solids steadily declines with time.

To predict the time evolution of dust in the terrestrial zone, we consider a disk with initial solid surface density ${{\rm{\Sigma }}}_{s}(a)={x}_{m}{{\rm{\Sigma }}}_{0}{(a/{a}_{0})}^{-3/2}$, inner radius ${a}_{\mathrm{in}}$ ≈ 0.1 au, and outer radius ${a}_{\mathrm{out}}$ ≈ 1–2 au. This inner radius lies between the value adopted in some numerical calculations (0.05 au; e.g., Hansen & Murray 2012, 2013) and the 0.25 au inner boundary for the Burke et al. (2015) analysis of Kepler data. Our results are insensitive to the exact value of ${a}_{\mathrm{in}}$. Disks with scale factor ${x}_{m}$ = 1 and ${{\rm{\Sigma }}}_{0}$ = 10 ${\rm{g}}\,{\mathrm{cm}}^{-2}$ at a₀ = 1 au have the surface density of the MMSN. For typical conditions in a disk with several protoplanets, ${r}_{\max }$ ≈ 100–1000 km and e ≈ 0.1 (Chambers 2008; Raymond et al. 2011; Kenyon & Bromley 2016). If ${r}_{\max }$ and e remain fixed throughout the evolution, the surface density declines roughly linearly with time, ${{\rm{\Sigma }}}_{s}{(t)\propto (1+t/{t}_{c})}^{-1}$, where t_c ∝ ${r}_{\max }$P/${{\rm{\Sigma }}}_{0}$ is the collision time and P is the orbital period. The relative luminosity of the debris disk is then an analytic function of the extent of the disk, the initial mass, and time (Appendix A.1).

At the start of the evolution, L_d/${L}_{\star }$ depends only on ${x}_{m}$ and the radial extent of the disk (Appendix, Equations (A14)–(A15)). The initial dust luminosity scales linearly with ${x}_{m}$ and with radial distance as a^−3/2. When ${a}_{\mathrm{in}}$ = 0.1 au and ${a}_{\mathrm{out}}$ = 2 au, systems with ${x}_{m}$ ≳ 2 × 10⁻⁴ have dust luminosity larger than the nominal detection limit, L_d/${L}_{\star }$ ≳ 10⁻⁴.

Once the evolution begins, the dust luminosity declines on the local collision time. For a disk with Σ_s ∝ a^−3/2, t_c ∝ a³. Although most of the mass is in the outer disk (${M}_{d}\propto {a}_{\mathrm{out}}^{1/2}$), the material closest to the star has the largest brightness per unit surface area. Solids close to the star also have the shortest collision time. Thus, the dust luminosity begins to decline on the collision time of the inner disk, ${t}_{\mathrm{in}}$. The luminosity declines by a factor of roughly two on this timescale (see also Kennedy & Wyatt 2010).

On timescales larger than ${t}_{\mathrm{in}}$, collisions remove material from larger and larger disk radii. At the outer edge of the disk, the collision time is ${t}_{\mathrm{out}}$. On timescales between ${t}_{\mathrm{in}}$ and ${t}_{\mathrm{out}}$, the dust luminosity declines rather slowly with time, L_d/${L}_{\star }$ ∝ t⁻ⁿ with n ≈ 0.3–0.9. Once the evolution time exceeds ${t}_{\mathrm{out}}$, the decline in the dust luminosity follows the decline of a narrow ring, L_d/${L}_{\star }$ ∝ t⁻ⁿ with n = 1.

At late times, the dust luminosity of disks with different ${x}_{m}$ converges (e.g., Wyatt & Dent 2002; Dominik & Decin 2003). Although the initial disk luminosity scales with ${x}_{m}$, disks with large ${x}_{m}$ evolve more rapidly than disks with small ${x}_{m}$. This convergent luminosity depends on ${r}_{\max }$, ${x}_{m}$, and the extent of the disk.

The lower panel of Figure 3 illustrates the evolution of L_d/${L}_{\star }$ at 0.1–50 Myr for 0.1–1 au debris disks with ${r}_{\max }$ = 300 km, e = 0.1, and a range of ${x}_{m}$. Disks with ${x}_{m}$ = 0.3 (30% of the MMSN) evolve from L_d/${L}_{\star }$ ≈ 10⁻¹ at 10³ yr to L_d/${L}_{\star }$ ≈ 10⁻² at 0.1 Myr to L_d/${L}_{\star }$ ≈ 10⁻³ at 3–4 Myr. These disks remain much brighter than the nominal detection limit until stellar ages of 30–40 Myr. At early times (t ≲ 0.1 Myr), the dust luminosity roughly scales with ${x}_{m}$; a disk with ${x}_{m}$ = 0.003 has a dust luminosity 100 times smaller than a disk with ${x}_{m}$ = 0.3. Because disks with smaller ${x}_{m}$ have longer collision times, they evolve more slowly. It takes a disk with ${x}_{m}$ = 0.003 roughly 15–20 Myr to decline from L_d/${L}_{\star }$ ≈ 10⁻³ to L_d/${L}_{\star }$ ≈ 10⁻⁴.

Zoom In Zoom Out Reset image size

The upper panel of Figure 3 shows that larger disks are always brighter (Equation (A15) of the Appendix). With a factor of 8 larger collision time at the outer edge, a disk extending to 2 au declines more slowly than a disk extending to 1 au. For disks with ${x}_{m}$ = 0.3, the larger disk reaches the detection limit at 100 Myr instead of 30–40 Myr. As ${x}_{m}$ decreases, however, the differences in evolution times become smaller. When ${x}_{m}$ = 0.003, the 0.1–2 au disk reaches the detection threshold at 25 Myr instead of 15–20 Myr for the 0.1–1 au disk.

To derive predicted flux ratios for debris disks, we assume that the solids radiate as blackbodies in equilibrium with radiation from the central star. For simplicity, ${L}_{\star }$ = 1 L_⊙ ${R}_{\star }$ = 1 R_⊙, and ${T}_{\star }$ = 5780 K at all times. The disk temperature is then T_d = 280 (a/a₀)^−1/2 K. With no analytic solution for the flux ratio, we divide the disk into a series of annuli, assign a temperature to each annulus, derive the time evolution of the cross-sectional area and emitted flux in each annulus, and add up the fluxes. In the Appendix, we show that numerical integrations of the dust luminosity agree well with analytic results.

Figure 4 summarizes results for 0.1–2 au disks with ${r}_{\min }$ = 1 μm, ${r}_{\max }$ = 300 km, and various ${x}_{m}$. At current sensitivity limits, robust 8–12 μm detections of warm dust at 10 Myr require massive debris disks with ${x}_{m}$ ≳ 0.3; detections at 8 μm are strongly favored over those at 12 μm. Warm dust is much easier to detect at 16–24 μm. Spitzer 16 μm measurements of ∼10 Myr old stars can detect debris disks with ${x}_{m}$ ≳ 0.03. The sensitivity at 24 μm is fairly remarkable: debris disks with ${x}_{m}$ as low as 1% of the MMSN are detectable.

Zoom In Zoom Out Reset image size

We can use these results to interpret statistics for warm debris from Section 2. Among stars younger than 10 Myr in the FEPS sample, 3%, 7%, and 20% of the nonprimordial disks have a detected 24 μm excess above 50%, 30%, and 10% of the stellar photosphere. From Figure 4, these results allow us to conclude that <3% of FEPS sources younger than 10 Myr have a warm excess consistent with ${x}_{m}$ ≳ 0.1. The 3% upper limit arises because cold dust beyond 2 au could produce some of the observed 24 μm excess. Similarly, fewer than 7% of sources younger than 10 Myr have a warm excess consistent with ${x}_{m}$ ≳ 0.03.

The FEPS 16 μm excess statistics place a more stringent limit on ${x}_{m}$. Because no FEPS source at any age (between 3 Myr and 3 Gyr) has a 16 μm excess above 16% of the stellar photosphere, all sources younger than 10 Myr must have ${x}_{m}$ ≲ 0.05. Similarly, all sources younger than 5 Myr have ${x}_{m}$ ≲ 0.03.

For a given ${x}_{m}$, reducing ${r}_{\max }$ reduces the predicted level of IR excess. Figure 5 shows how the 24 μm excess varies with ${r}_{\max }$ for 0.1–2 au debris disks with ${x}_{m}$ = 0.1 and ${r}_{\min }$ = 1 μm. In systems with identical initial masses, swarms of particles with smaller ${r}_{\max }$ have larger initial ${A}_{d}$ and smaller collision times. Although swarms with ${r}_{\max }$ = 3–10 km initially have larger 24 μm excesses than swarms with ${r}_{\max }$ = 1000 km, they also evolve much more rapidly. By 6 Myr (12 Myr), debris disks with ${r}_{\max }$ = 3 km (10 km) reach the nominal detection limit. Disks with ${r}_{\max }$ = 100–1000 km take 40–100 Myr to reach this limit.

Zoom In Zoom Out Reset image size

Numerical calculations of planet formation suggest that the agglomeration of Earth-mass planets in the terrestrial zone creates copious amounts of debris. Current estimates suggest that from 15% to 30% of the initial mass is potentially available to produce the debris signature of rocky planet formation. Debris production is probably independent of the mode—classical, in situ, or pebble—in which planets form. Unless planetesimal formation mechanisms assemble Mars-mass oligarchs directly, leaving behind 10% of the initial mass in debris is inevitable (e.g., Kenyon & Bromley 2016). By generating smaller planetesimals, current formation paths appear to preclude this possibility (Johansen et al. 2015; Simon et al. 2016). All popular models for the formation of Earth-mass planets include a giant impact phase that also generates significant amounts of debris.

The debris is long-lived. In standard models of planet formation, the debris has a typical maximum size of 100–1000 km (e.g., Kenyon & Bromley 2016). If the debris contains 15%–30% of the initial mass required to assemble an Earth-mass planet, then the IR excess from the debris is detectable at 24 μm for ≲100 Myr (Figure 4; see also Raymond et al. 2011; Genda et al. 2015b).

For stars with ages ≲10 Myr, it is challenging to detect debris from terrestrial planet formation with current sensitivity limits at 8–12 μm. For robust detections, recent surveys require an initial solid mass ≳30% of the MMSN at 8 μm and ≳100% of the MMSN at 12 μm. If the formation of several Earth-mass planets converts 15%–30% of the initial mass into debris, detection at 12 μm is very unlikely: the debris simply is not luminous enough after 10 Myr. At 8 μm, the predicted level of debris from the formation of a single Earth-mass planet is somewhat smaller than the Spitzer limits. If Earth-mass planets are as common as suggested by Kepler (Burke et al. 2015), the low frequency, ≲3% (e.g., Luhman & Mamajek 2012), of warm dust revealed by 8 μm observations of solar-type stars in the Upper Sco association is roughly consistent with theoretical predictions. However, if we include the excesses produced by the formation of more massive planets (≳1 ${M}_{\oplus }$), the expected 8 μm excess is larger than Spitzer sensitivity limits and potentially in conflict with observations.

If our nominal picture of debris production is correct (Figure 4), current 16–24 μm sensitivity limits are low enough to detect every solar-type star ≲10 Myr engaged in forming Earth-mass planets. For 10 Myr old stars, current 16 μm data should detect debris disks with initial masses of 5%–10% of the MMSN; however, the FEPS survey detected none of these systems. Similarly, <3% of FEPS sources younger than 10 Myr have a warm excess consistent with initial solid masses ≳10% of the MMSN.

Based on this analysis, there is a clear discrepancy between the observed frequency of warm debris from terrestrial planet formation (≲3%; Section 2) and the incidence rate of Earth-mass planets derived from Kepler (≳20%; Section 3). Based on Figure 4, the observed frequency of 16–24 μm excess sources and the inferred incidence rate of Earth-mass planets imply that terrestrial planet formation must leave behind a very small fraction of the initial mass in debris-producing solids (≲1% of MMSN).

When the frequency of false positives f_fp among Kepler Earth-mass planet candidates is large, we overestimate f_p. Current analyses infer f_fp ≈ 10%–75% (e.g., Morton & Johnson 2011; Santerne et al. 2012; Fressin et al. 2013; Sliski & Kipping 2014; Colón et al. 2015; Désert et al. 2015; Coughlin et al. 2016; Morton et al. 2016; Santerne et al. 2016). Many studies focus on gas giant planet candidates; however, Désert et al. (2015) consider whether Spitzer observations confirm transits for Earth-mass candidates with orbital periods ≲100 days. For individual systems, they derive f_fp ≈ 1%–42%; the complete sample suggests a typical f_fp ≈ 10%. Curiously, additional observations confirm planets in several candidate systems with large individual f_fp. Although larger samples are required to understand the true f_fp for Earth-mass planet candidates, these data suggest that f_p is not overestimated by a factor of 10 (see also Coughlin et al. 2016; Morton et al. 2016).

If multiplanet systems occur commonly and with a large range in orbital inclination ı, the true f_p is smaller than our estimates. Although transit surveys of such systems may detect only one of two or more possible planets, detecting a transiting planet in the system is more likely because a planet can be detected from multiple directions. As a result, the fraction of stars with planets in the habitable zone f_p can be reduced while maintaining the same average number of Kepler planets per star overall.

Although observations do not yet provide a robust estimate of the frequency of multiplanet systems (e.g., Tremaine & Dong 2012; Fabrycky et al. 2014; Van Eylen & Albrecht 2015; Winn & Fabrycky 2015, and references therein), comparing results from radial velocity and transit surveys provides useful constraints (S. Tremaine 2016, private communication). Transit surveys can detect multiplanet systems only when the range of mutual inclinations is roughly a few degrees or less; however, radial velocity surveys are sensitive to systems with a much broader range of mutual inclinations. The frequency of multiplanet systems from recent Kepler catalogs (0.19–0.20; Fabrycky et al. 2014; Mullally et al. 2015; Coughlin et al. 2016) is 2/3 of the frequency of ∼0.3 derived from radial velocity measurements (Limbach & Turner 2015). The similarity between the two rates implies that Kepler does not significantly underestimate the fraction of multiplanet systems, i.e., mutual inclinations are small and planetary systems are fairly flat.

While this argument applies to mutual inclinations of planets of all masses, additional dynamical considerations can restrict the mutual inclinations of systems of Earth-mass planets. For known multiplanet systems, the typical range in ı is small, ≈1°–5° (Tremaine & Dong 2012; Fabrycky et al. 2014). With e ≈ (2–4)ı (Fabrycky et al. 2014), e ≈ 0.035–0.20. If this range in e is typical of systems with Earth-mass planets, theory provides a guide to estimate a possible error in f_p at a = 0.25–1 au (see also Tremaine 2015). From numerical simulations, systems of several Earth-mass planets have e ≈ 0.01–0.05 and ı ≈ e/2 (e.g., Chambers & Wetherill 1998; Bromley & Kenyon 2006; Raymond et al. 2007; Morishima et al. 2008; Chambers 2013). These systems are dynamically unstable when the distance between the apocenter of the inner orbit and the pericenter of the outer orbit exceeds 10 mutual Hill radii R_H, where ${R}_{{\rm{H}}}=({a}_{\mathrm{in}}+{a}_{\mathrm{out}}){[({m}_{\mathrm{in}}+{m}_{\mathrm{out}})/3{M}_{\star }]}^{1/3}$ and a_in (a_out) is the semimajor axis of the inner (outer) planet (Chambers et al. 1996; Yoshinaga et al. 1999; Fang & Margot 2013; Petrovich 2015; Pu & Wu 2015). When all planets have e = 0.03 (e = 0.10), the maximum number of stable Earth-mass planets within 0.4–1 au is seven to eight (four to five). The maximum reduction in f_p is then a factor of 5–7.

The number of transiting planets in a multiplanet system depends on the distribution of mutual inclinations (e.g., Tremaine & Dong 2012; Fabrycky et al. 2014). To make estimates for ensembles of closely packed Earth-mass planets, we perform a simple Monte Carlo calculation. For adopted dispersions in e and ı, σ_e and σ_ı = σ_e/2, our algorithm establishes a set of Earth-mass planets with orbital separations of 10 mutual Hill radii and semimajor axes between 0.25 and 1 au. After randomly selecting one of these planets to have impact parameter b ≤ 1, the code chooses random deviates for ı, infers impact parameters for the remaining planets, and counts the number of planets with $| b| \,\leqslant$ 1. For Gaussian (Rayleigh) deviates with σ_e ≈ 0.035 and σ_ı ≈ 1°, closely packed systems have seven (six) planets and four (three) transits per system. When σ_e ≈ 0.175 and σ_ı ≈ 5°, tightly packed systems have three to four (two to three) planets and one (one) transit per system. Kepler should detect maximally packed Earth-mass planets on roughly circular orbits as multiplanet systems, but should fail to identify multiple planets in highly eccentric systems. Either way, these results suggest that f_p might be overestimated by factors of 2–4. Factor-of-10 overestimates are unlikely.

Overall, this analysis suggests that overestimates in f_p are unlikely to reduce the fraction of stars with planets in the 0.25–1 au region (∼20%) to the fraction of stars with warm excesses (<3%). Moreover, our estimates are conservative. If there is a population of Earth-mass planets with a ≈ 1–2 au, then f_p is larger than 20%. If the typical σ_e of multiplanet systems is roughly 4σ_ı instead of our adopted 2σ_ı (Fabrycky et al. 2014), then the maximum number of Earth-mass planets in Kepler systems with a single transit is roughly two. Thus, uncertainties in f_p cannot reduce the discrepancy between f_p and f_d.

Alternatively, we may have overestimated the expected IR excess produced by debris. To derive this excess (Section 4), we assume grains with sizes ≳1 μm. In any debris disk model, the dust luminosity is sensitive to the size of the smallest grains, L_d ∝ ${r}_{\min }^{-1/2}$ . Grains with sizes smaller (larger) than our nominal limit of 1 μm are more (less) luminous than predicted. Although ${r}_{\min }$ is clearly a free parameter, observations of comets in the solar system and warm debris around solar-type stars suggest that 1 μm is a reasonable lower limit to the grain size (e.g., Lisse et al. 2006, 2007, 2008; Currie et al. 2011). Thus, increasing ${r}_{\min }$ seems an unlikely way to remove the discrepancy between f_d and f_p.

We also assume disks with power-law size distributions of small particles. In a real collisional cascade, the equilibrium size distribution has pronounced waves relative to our adopted $N(r)\propto {r}^{-3.5}$ power law (Campo Bagatin et al. 1994; O'Brien & Greenberg 2003; Wyatt et al. 2011; Kenyon & Bromley 2016). Adopting an analytic model for the waves that matches numerical simulations (Kenyon & Bromley 2016) yields IR excesses a factor of 2–3 larger than the predictions of the basic analytical model outlined in the Appendix. With this assumption, it is possible to detect debris disks with 10%–30% (8–12 μm), 1% (16 μm), or 0.3% (24 μm) of the MMSN.

We also assume that small grains radiate as perfect blackbodies. In a real debris disk, small (1–10 μm) grains are probably hotter and radiate less efficiently than perfect blackbodies (e.g., Backman & Paresce 1993; Dent et al. 2000; Lisse et al. 2006). If we adopt an extreme model where all grains radiate inefficiently (e.g., with emissivity ∝ (λ/λ₀)^−b, with λ₀ = 1 μm and b ≈ 0.8–1), the IR excess at 8–24 μm is roughly 50% smaller than the blackbody prediction. At 8–12 μm, this prediction has little impact on our ability to detect excesses around 10 Myr stars: blackbody grains are already at or below the nominal detection limits for Spitzer and WISE. At 16–24 μm, systems with initial masses of 10% (16 μm) or 1% (24 μm) of the MMSN still remain above the detection threshold even with an extreme emission model.

Considering the uncertainties in ${r}_{\min }$ and the radiative properties and size distribution of small particles, our estimates for the IR excesses of warm debris disks seem reasonable. Allowing for changes in ${L}_{\star }$ from pre-main-sequence stellar evolution (e.g., Baraffe et al. 1998, 2015; Siess et al. 2000; Bressan et al. 2012; Chen et al. 2014) and corresponding variations in dust temperature also has little impact on predicted IR excess emission at 8–24 μm. Thus, the discrepancy between f_d and f_p remains.

If we accept the currently measured frequencies of Earth-mass planets and warm debris disks, we must then devise a theory that assembles Earth-mass planets with little detectable debris.

One way to achieve this goal is if planet formation is quick. If Earth-mass planets orbiting solar-mass stars reach their final masses during the T Tauri phase, debris ejected from giant impacts or collisions of leftover planetesimals would be rapidly mixed with the optically thick primordial dust left within the gaseous disk. Instead of producing a distinctive "debris disk signature," the debris simply contributes to the already large IR excess of the primordial disk. As viscous evolution and other processes remove the gas, the debris from Earth-mass planet formation leaves with the gas.

In the classical theory, planet formation is not quick. Although protoplanet formation might occur during the T Tauri phase (e.g., Wetherill & Stewart 1993; Weidenschilling 1997; Ohtsuki et al. 2002; Kenyon & Bromley 2006, 2016; Kokubo et al. 2006; Kokubo & Genda 2010), growth then stalls (e.g., Weidenschilling et al. 1997; Chambers 2001, 2008, 2013; Kenyon & Bromley 2006; Raymond et al. 2007, 2014; Genda et al. 2015b; Quintana et al. 2016). Once the gaseous disk dissipates, giant impacts among Mars-mass protoplanets then build Earth-mass planets (Kominami & Ida 2004), generating debris. The gaseous disk plays no further role in the evolution of the debris.

Pebble accretion is not intrinsically quick (e.g., Chambers 2014; Johansen et al. 2015; Levison et al. 2015). Recent studies suggest the rapid formation of 100–1000 km planetesimals at ≲1 Myr. These planetesimals rapidly accrete leftover pebbles. Because these models begin with initial masses comparable to the MMSN, the giant impact phase begins well after the gaseous disk has dissipated. Giant impacts then occur on timescales when debris is easily detected.

The massive disks of solids proposed in the in situ theory (∼10 times the MMSN) enable rapid planet formation (Hansen & Murray 2012, 2013; Hansen 2015). Formation times scale inversely with the mass of solids; thus, protoplanets grow much faster than in the classical or pebble accretion theories. Even in a primordial gaseous disk, the large number of closely packed protoplanets generates a series of giant impacts, leading to multiple Earth-mass planets in a few megayears. The gas can then remove small particles produced in giant impacts. However, the gas probably cannot remove any leftover planetesimals or other 1–100 km particles generated during the giant impact phase. To produce a 24 μm excess that is <10% of the stellar photosphere in a system where the initial mass in solids is ∼10 times the MMSN, the giant impact phase must leave behind less than 0.1%–0.2% of the initial mass in 100–1000 km objects (including any leftover planetesimals). Otherwise, the in situ theory makes too much debris at 10 Myr.

Thus, another way to satisfy the observational constraints is if planet formation is intrinsically neat: Earth-mass planets assemble with negligible dust emission. As described above, the classical, in situ, and pebble accretion scenarios are not neat enough to match observational constraints. To isolate the appropriate physical conditions for a "neat" scenario, we rely on the analytical model outlined in Appendix A.1 and published numerical simulations (e.g., Kenyon & Bromley 2016). To match current sensitivity limits at 16 μm and 24 μm, we must reduce dust production by at least a factor of 10 (i.e., reducing ${x}_{m}$ from 15%–30% to <3%; Figure 4). In all numerical simulations, dust production depends on v_imp/${Q}_{D}^{\star }$, where v_imp is the impact velocity and ${Q}_{D}^{\star }$ is the binding energy. Simple dynamics sets v_imp; changing it significantly is unlikely.

The binding energy is based on analytical and numerical simulations using state-of-the-art equations of state. While increasing ${Q}_{D}^{\star }$ seems unlikely, it is worth exploring the required physics. As a simple comparison, the gravitational binding energy per unit mass of a uniform sphere exceeds the ${Q}_{D}^{\star }$ for basalt only when r ≳ 10 ${R}_{\oplus }$. Thus, the internal degrees of freedom of protoplanets are important in setting ${Q}_{D}^{\star }$. While it may be possible to modify these by factors of 2–3, order-of-magnitude changes seem unlikely (see also Housen & Holsapple 1999; Holsapple et al. 2002, and references therein).

If dust production cannot be limited, the main alternative is to reduce the lifetime of the debris. In this approach, a shorter debris lifetime lessens the likelihood of detection. From Equations (A4) and (A10), lowering ${r}_{\max }$ or increasing v_imp/${Q}_{D}^{\star }$ by a factor of 10 shortens the lifetime of the debris by a similar factor. Having already ruled out factor-of-10 changes in v_imp/${Q}_{D}^{\star }$, we consider the impact of changing ${r}_{\max }$.

In the analytical model, changing ${r}_{\max }$ modifies the detectability of the debris in two ways. For a fixed mass in debris, the cross-sectional area A_d ∝ ${r}_{\max }$^−1/2. Reducing ${r}_{\max }$ therefore increases the initial dust luminosity and shortens the collision time. As shown in Figure 5, ${r}_{\max }$ must be quite small to reduce the predicted warm excess to the level of the 24μm excess displayed by the brightest 20% of ∼10 Myr old stars. If ${x}_{m}$ = 0.1 and ${r}_{\max }$ = 10 km, then at 10 Myr F_d/F_star ∼ 0.12 at 24 μm. Approximately 20% of FEPS sources younger than 10 Myr have an excess brighter than F_d/F_star ∼ 0.1 at 24 μm, a fraction comparable to the fraction of solar-type stars with Earth-like planets.

Although this solution seems attractive, it has several major drawbacks. Numerical simulations of giant impacts already yield 24 μm fluxes much larger than observed (Genda et al. 2015b); generating a smaller likelihood of detection at 10 Myr age in exchange for a much larger initial luminosity is probably a poor trade. For dust generated from planetesimals, numerical simulations suggest that collisional damping of small particles maintains large dust luminosities far longer than suggested by the analytical model (Kenyon & Bromley 2016). Increasing the population of small particles by lowering ${r}_{\max }$ probably exacerbates this problem.

If we cannot substantially alter the observed frequency of terrestrial planets, models for the assembly of terrestrial planets, or the amount of debris generated by terrestrial planets, some process must remove small particles with r ≲ 0.1–1 mm from the terrestrial zone. Possibilities include interactions with the stellar radiation field, the stellar wind, or a remnant gaseous disk. If one of these processes removes small particles faster than the collisional cascade produces them, the predicted IR emission from warm dust is reduced by the requisite factor of ∼10, eliminating the discrepancy between f_d and f_p.

Compared to aerodynamic drag from a remnant gaseous disk, other mechanisms probably have a limited role. Observations suggest that radiation pressure may remove particles smaller than 1 μm, but larger particles are robustly detected (e.g., Lisse et al. 2008; Currie et al. 2011; Matthews et al. 2014, and references therein). For L_d/${L}_{\star }$ ≳ 10⁻⁷, Poynting–Robertson drag is ineffective (Wyatt 2008). In the inner solar system, stellar wind drag is roughly 30% as effective as Poynting–Robertson drag in removing small particles (Gustafson 1994). Observations of the youngest solar-type stars suggest mass-loss rates ≲10 times the mass-loss rate ${\dot{M}}_{\odot }$ of the current Sun (Wood et al. 2014, and references therein). Theoretical studies predict mass-loss rates up to 100 ${\dot{M}}_{\odot }$ (e.g., Cohen & Drake 2014; Airapetian & Usmanov 2016). Both of these estimates fall well below the 1000 ${\dot{M}}_{\odot }$ required for stellar wind drag to remove small particles rapidly when L_d/${L}_{\star }$ ≳ 10⁻⁴. Finally, the sputtering rate for small particles is too small compared to the collision rate (Wurz 2012).

In some circumstances, photophoresis⁶ drives small particles radially outward through an optically thin gaseous disk (e.g., Krauss & Wurm 2005; Herrmann & Krivov 2007; Cuello et al. 2016). When photophoresis is effective, the timescale for radial drift is comparable to (and sometimes shorter than) the timescale for radial drift due to gas drag. However, the drift time is sensitive to the ratio of the heat conductivity to the particle asymmetry factor, which is uncertain (e.g., von Borstel & Blum 2012).

Here we focus on the impact of gas drag, where the physical properties of small particles are less important. Our goal is to identify a range of Σ_g for the gaseous disk that allows giant impacts (0.1%–1% or less of the MMSN) and removes 0.1–1 mm and smaller particles on timescales shorter than the collision time.

Within any gaseous disk, particles weakly bound to the gas rapidly drift inward when their "stopping time" (${t}_{s}={{mv}}_{d}/{F}_{d}$, where v_d is the drift velocity and F_d is the drag force) is comparable to their orbital period P (Adachi et al. 1976; Weidenschilling 1977b). At 1 au, the shortest drift time is 50–100 yr (Figure 6), shorter than the typical collision time of 1000 yr for 1 μm to 1 cm particles. In an optically thin disk, radiation pressure drives small, weakly coupled particles to large a, where the particles are colder and have much smaller dust luminosities (Takeuchi & Artymowicz 2001). Effective radial drift thus eliminates particles that produce an IR excess.⁷

Zoom In Zoom Out Reset image size

To quantify radial drift at 1 au in a low-density gaseous disk, we combine the approaches of Weidenschilling (1977b) and Takeuchi & Artymowicz (2001). As outlined in Appendix A.2, we solve for the radial and azimuthal velocity of particles relative to the gas in a protostellar disk with standard relations for the surface density, midplane temperature, vertical scale height, gas pressure, thermal velocity, and other physical variables (e.g., Kenyon & Hartmann 1987; Chiang & Goldreich 1997; Rafikov 2004; Chiang & Youdin 2010; Armitage 2013; Youdin & Kenyon 2013).

Figure 6 illustrates the impact of a residual gas disk on the radial drift velocity of particles as a function of size at 1 au under MMSN-like conditions. For this example, a disk with the surface density of the MMSN has Σ_g = 2000 ${\rm{g}}\,{\mathrm{cm}}^{-2}$, T = 278 K, and vertical scale height H = 0.03 au. Filled circles indicate inward drift; open symbols indicate outward drift. Other choices for these parameters—e.g., T = 150–300 K and H = 0.02–0.10—change the maximum drift velocity and the particle size for this maximum drift by 25%–50% but do not change the overall trends.

When the disk has Σ_g comparable to the MMSN (Figure 6, black symbols), particles with r ≈ 50 cm have τ_s ≈ 1 and the largest drift velocity. Inflection points in ${v}_{\mathrm{rad}}(r)$ occur when particles enter different drag regimes (e.g., Epstein, Stokes, quadratic). For particle sizes r ≈ 30–40 μm, outward radiation pressure balances inward gas drag; these particles do not drift. Smaller particles drift outward⁸ with maximum velocities slightly smaller than 0.1 cm s⁻¹.

As the surface density of the gaseous disk declines, smaller particles drift more rapidly through the gas. Although the peak in ${v}_{\mathrm{rad}}$ shifts to smaller sizes, the maximum drift velocity is always roughly 8000 cm s⁻¹. The balance between radiation pressure and gas drag remains fixed for 30–40 μm particles. However, in disks with lower Σ_g, smaller particles have larger stopping times. With a weaker drag force, radiation pressure drives these particles outward at larger velocities. Once Σ_g is roughly 0.001% of the MMSN, the outward drift velocities of r ≲ 30 μm particles surpass the maximum inward drift velocities of 100 μm particles.

To judge the impact of radial drift on the collisional cascade, we compare the drift time to the collisional time. For disks with surface densities of 0.001% of the MMSN, 1–10 μm and 100 μm particles at 1 au drift inward or outward on 50–75 yr timescales. In less (more) tenuous gas disks, particle lifetimes are longer (shorter; Figure 6). For comparison, in a collisional cascade at 0.1–2 au with L_d/${L}_{\star }$ ≈ 10⁻⁴, the collision time is roughly 1000 yr. In disks with Σ_g ≈ 0.001% of the MMSN, radial drift from radiation pressure and gas drag removes particles smaller than 100 μm in <100 yr, rapidly depleting the disk of the 1–100 μm particles that compose 90% of the surface area. Gas drag also rapidly transports larger (0.1–1 mm) particles toward the central star.

These two processes have a dramatic effect on debris production. In a standard collisional cascade, mass flows at a roughly constant rate from the largest particles to the smallest particles (e.g., Wyatt 2008; Kobayashi & Tanaka 2010; Wyatt et al. 2011). When the gas drags 0.1–1 mm particles inward, the equilibrium population of these particles drops. Collisions among these objects are then less frequent, depressing the flow of mass to smaller particles. In turn, a smaller mass flow rate lowers the cross-sectional area of the swarm and weakens the infrared excess.

When gas drag and radiation pressure effectively remove small particles inside 2 au, the IR excess drops dramatically. For disks with 0.001%–0.01% of the MMSN, the dust luminosity drops by 2–4 orders of magnitude. If this remnant disk extends close to the stellar photosphere, particles with sizes ≳1 mm either evaporate or fall onto the photosphere. Thus, residual gas disks as dilute as 0.001% of the MMSN can remove the small solids (<100 μm) that make up 90% of the cross-sectional area of the debris, reducing the IR excess at 5–20 Myr to a level consistent with the observational constraints. Less dilute gas disks (0.01% of the MMSN) can reduce the population of small grains <100 μm indirectly, by removing the 1–10 mm solids that generate them through collisions. More dilute gas disks are less efficient in removing solids >10 μm and therefore less able to limit the IR excess.

We conclude that a residual gas disk in the range 0.001%–1% of the MMSN that persists for ∼10 Myr is dilute enough to allow giant impacts to assemble terrestrial planets in classical planet formation, but dense enough to minimize the observable IR excess produced by the resulting debris.

Current observational constraints on the surface density of residual circumstellar gas disks cannot exclude our picture for reducing IR excesses around 10 Myr old solar-type stars. Constraints on the residual gas content of disks from in situ diagnostics are model dependent (Gorti & Hollenbach 2004) and not highly restrictive in this context. For example, analyses of Spitzer IRS emission-line diagnostics for a sample of FEPS targets with little evidence for warm debris disks place a formal limit on the gas surface density Σ_g < 0.01% of the MMSN beyond 1 au (Pascucci et al. 2006). If this upper limit were to extend inside 1 au, the surface density would be low enough to allow giant impacts in the classical planet formation picture (e.g., Kominami & Ida 2002, 2004) and large enough to allow the removal of small particles by gas drag.

Measuring stellar accretion is an alternate way to search for evidence of a residual gas disk. While accretion indicates that residual gas is present, it does not directly measure Σ_g. To obtain a rough scaling between stellar accretion rate ${\dot{M}}_{\star }$ and Σ_g, we can consider the average accretion rates of classical T Tauri stars at 1–10 Myr age. Stellar accretion rates decrease from ∼10⁻⁸ M_⊙ yr⁻¹ at the age of Taurus-Auriga (few megayears) to ∼4 × 10⁻¹⁰M_⊙ yr⁻¹ at 10 Myr (Sicilia-Aguilar et al. 2006), i.e., only a factor of 25. If the fractional decrement in ${\dot{M}}_{\star }$ indicates a similar reduction in Σ_g, an accretion rate of ∼4 × 10⁻¹⁰M_⊙ yr⁻¹ corresponds to a surface density at 1 au of 4 g cm⁻² if the initial state is an α-disk with a surface density of 100 g cm⁻² at 1 au (Hartmann et al. 1998) or 80 g cm⁻² if the initial state is the MMSN with a surface density of 2000 g cm⁻² at 1 au (Weidenschilling 1977a). Nonaccreting sources, with presumably lower Σ_g, make up more than 90% of young stars at 5–10 Myr age (Fedele et al. 2010).

Based on the ${\dot{M}}_{\star }$–Σ scaling relation, it appears that we need to probe very low accretion rates ∼10⁻¹³−10⁻¹²M_⊙ yr⁻¹ to reach Σ_g ≲ 10⁻⁵–10⁻⁴ of the MMSN and to exclude the possibility that 5–10 Myr old disks contain enough residual gas to erase a debris disk signature. For approximately solar-mass stars with ages of 5–10 Myr, chromospheric emission limits our ability to measure stellar accretion rates below ∼10⁻¹⁰M_⊙ yr⁻¹ using U-band excesses (Ingleby et al. 2011a) or hydrogen emission line fluxes (Manara et al. 2013). Hα emission line widths (Natta et al. 2004) and line profile analyses (Muzerolle et al. 1998a, 1998b) can probe lower accretion rates (e.g., Muzerolle et al. 2000; Manara et al. 2013).

From the shape and velocity extent of the Hα line profile (White & Basri 2003; Natta et al. 2004), Riviere-Marichalar et al. (2015) derive a broad range of accretion rates for 11 sources in the 8 Myr old η Cha group. Six sources have ${\dot{M}}_{\star }$ ≳ 10⁻¹⁰M_⊙ yr⁻¹ from at least one spectrum; the other five sources have ${\dot{M}}_{\star }$ ≈ 10⁻¹²–3 × 10⁻¹¹M_⊙ yr⁻¹. While it may not be appropriate to apply the Natta et al. (2004) relation derived for low-mass stars to a set of solar-type stars, these accretion rates are close to the rates required for gas drag to eliminate the IR excesses of debris disks.

There are multiple prospects for detecting a dilute reservoir of gas in solar-type stars. CO fundamental emission (e.g., Najita et al. 2003; Salyk et al. 2011; Doppmann et al. 2016) and UV transitions of molecular hydrogen (e.g., Ingleby et al. 2011b) probe molecular gas within an au of young stars. Among A-type stars with debris disks, several have neutral or ionized gas close to the star (e.g., Hobbs et al. 1985; Ferlet et al. 1987; Welsh et al. 1998; Redfield 2007; Montgomery & Welsh 2012; Cataldi et al. 2014; Kiefer et al. 2014, and references therein). Evaporation of comets (Lagrange et al. 1987; Beust et al. 1990) and vaporization of colliding particles (e.g., Czechowski & Mann 2007) might supply some of this material (see also Kral et al. 2016).

In our analysis, we have considered four main options to resolve the discrepancy between the frequencies of warm debris disks f_d and Earth-mass planets f_p among solar-type stars. Despite clear gaps in our understanding of the physics of terrestrial planet formation, there is no obvious way to modify the theories to limit the detectability of warm debris at 5–20 Myr. Uncertainties in our ability to relate IR excesses to debris from rocky planet formation are unlikely to change f_d. Although overestimating the observed frequency of Earth-mass planets is a plausible cause for the difference between f_d and f_p, the "best" explanation is rapid drift of small particles in a residual gaseous circumstellar disk.

Future observations and theoretical investigations can test all four explanations. We outline several possibilities.

• 1.  
  Search for fainter warm excess signatures in large samples of young solar-type stars. Current sensitivity limits are not much lower than theoretical predictions. If warm excesses are simply a factor of 2–3 fainter than predicted, detecting these excesses with more sensitive surveys (e.g., using MIRI on James Webb Space Telescope [JWST]; Wright et al. 2003; Wells et al. 2006; Bouchet et al. 2015) may be able to distinguish neat scenarios from quick scenarios. For example, placing stricter limits on the level of 16 μm excess (≲10% of the stellar photosphere) would probe initial solid masses in debris of a few percent of the MMSN (Figure 4).
• 2.  
  In Kepler systems with apparent Earth-mass planets, better limits on (i) the false-positive rate, (ii) the frequency of multiple planets, and (iii) the distributions of e and ı in systems with multiple Earth-mass planets would reduce the uncertainties in the fraction of stars that host Earth-mass planets inside 1–2 au. Theoretical studies into the stability of systems with multiple Earth-mass planets would place constraints on the ability of these systems to "hide" from radial velocity and transit observations.
• 3.  
  Hunt for planets around all young (5–20 Myr) solar-type stars. Several recent studies identify massive planets orbiting several pre-main-sequence stars (e.g., David et al. 2016; Donati et al. 2016; Gaidos et al. 2016; Johns-Krull et al. 2016a, 2016b; Mann et al. 2016). Identifying Earth-mass planets orbiting young stars would constrain the timescale of terrestrial planet formation relative to the production of debris. The K2 mission (Howell et al. 2014) might identify more short-period planets within several nearby young stellar associations.⁹ TESS (Ricker et al. 2015; Sullivan et al. 2015) can search for planets with a much larger range of orbital periods.
• 4.  
  Search for evidence of tenuous residual gas disks, at the 10⁻⁵–10⁻² of MMSN level, around young solar-type stars. Direct detection of gas or robust measurement of very low mass accretion rates tests the idea that radiation pressure and aerodynamic drag remove the debris of terrestrial planet formation. In addition to surveys with the NIRSPEC or MIRI spectroscopic instruments on JWST (e.g., Wright et al. 2003; Wells et al. 2006), sensitive optical (e.g., Xu et al. 2016) or radio (e.g., ALMA, VLA) observations might reveal low-mass gaseous disks in the terrestrial zones of 5–20 Myr old solar-type stars.
• 5.  
  Examine quick and neat modes of planet formation. Although our current understanding appears to preclude these ideas, it is important to quantify the ability of planet formation scenarios to assemble planets quickly and neatly in the absence of a long-lived gaseous disk. Future theoretical calculations of terrestrial planet formation should include clear predictions of dust production for comparison with observations (e.g., Raymond et al. 2011; Genda et al. 2015b; Kenyon & Bromley 2016).
• 6.  
  Consider the late stages of protostellar disk evolution in more detail. Current mechanisms for disk dispersal (e.g., photoevaporation and viscous evolution) do not make firm predictions for the structure of gaseous material on timescales of 10–100 Myr (e.g., Gorti et al. 2015, and references therein). For example, including the gas from the evaporation of comets or vaporization of colliding planetesimals in photoevaporation models would help to constrain the long-term evolution of debris and gas in the terrestrial zones of young stars.

Wyatt & Dent (2002) and Dominik & Decin (2003) first developed an analytic theory to track the long-term evolution of collisional cascades (see also Wyatt et al. 2007a, 2007b, 2011; Heng & Tremaine 2010; Kobayashi & Tanaka 2010; Kenyon & Bromley 2016). In this model, particles with sizes ranging from ${r}_{\min }$ to ${r}_{\max }$, orbital period P, and orbital eccentricity e occupy a single annulus with width δa at a distance a from a central star with mass ${M}_{\star }$, radius ${R}_{\star }$, and luminosity ${L}_{\star }$. Particles collide at a rate nσv, where n is the number density, σ is the physical cross-section, and v is the collision velocity. Setting the center-of-mass collision energy¹⁰ Q_c larger than the binding energy ${Q}_{D}^{\star }$ of the largest particles ensures destructive collisions that produce debris. For an initial surface density of solids Σ₀ in the annulus, the collision time is

$$\begin{eqnarray}&&{t}_{c}=\displaystyle \frac{\alpha {r}_{\max }\rho P}{12\pi {{\rm{\Sigma }}}_{0}},\end{eqnarray} \tag{ A1 }$$

where the correction factor is $\alpha \approx {\alpha }_{0}{({v}^{2}/{Q}_{D}^{\star })}^{-n}$ (Wyatt et al. 2007a, 2011; Heng & Tremaine 2010; Kobayashi & Tanaka 2010). The surface density of material in the annulus then evolves as

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{s}(t)=\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{1+t/{t}_{c}}.\end{eqnarray} \tag{ A2 }$$

Other collective properties of the particles—including the total mass, cross-sectional area, and relative luminosity—follow the decline of Σ_s with time (Wyatt & Dent 2002; Dominik & Decin 2003; Wyatt et al. 2007a, 2011; Kobayashi & Tanaka 2010; KB2016).

Here we expand the analytic theory to a disk¹¹ extending from an inner radius ${a}_{\mathrm{in}}$ to an outer radius ${a}_{\mathrm{out}}$. Particles with sizes ranging from ${r}_{\min }$ to ${r}_{\max }$ have an initial surface density

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{s}={x}_{m}{{\rm{\Sigma }}}_{0}{(a/{a}_{0})}^{-p},\end{eqnarray} \tag{ A3 }$$

where Σ₀ = 10 ${\rm{g}}\,{\mathrm{cm}}^{-2}$ is the surface density of solids at a₀ = 1 au in the MMSN and ${x}_{m}$ is a scale factor (Weidenschilling 1977a; Hayashi 1981; Kenyon & Bromley 2008). We assume that the correction factor $\alpha ={\alpha }_{0}{(a/{a}_{0})}^{-n}$. Setting $P={P}_{0}{(a/{a}_{0})}^{3/2}$, the collision time is then ${t}_{c}={t}_{0}{(a/{a}_{0})}^{3/2+p-n}$, where

$$\begin{eqnarray}&&{t}_{0}=\displaystyle \frac{{\alpha }_{0}{r}_{\max }\rho {P}_{0}}{12\pi {x}_{m}{{\rm{\Sigma }}}_{0}}.\end{eqnarray} \tag{ A4 }$$

To derive the evolution of the disk, we assume that the particles have a differential size distribution $N(r)\propto {r}^{-3.5}$ in each annulus (e.g., Dohnanyi 1969; Williams & Wetherill 1994; Tanaka et al. 1996; Kobayashi & Tanaka 2010). Other choices lead to similar results (Wyatt et al. 2011). We can then relate the cross-sectional area of the swarm to the mass in each annulus:

$$\begin{eqnarray}&&{dA}=\displaystyle \frac{3}{4\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}{dM}.\end{eqnarray} \tag{ A5 }$$

Setting ${dM}={{\rm{\Sigma }}}_{s}2\pi {ada}$, $\tilde{a}=a/{a}_{0}$, and $d\tilde{a}={da}/{a}_{0}$, the general result of the analytic model in Equation (A2) and our relation for t_c yields

$$\begin{eqnarray}{dA}(t) & = & \displaystyle \frac{3\pi {x}_{m}{{\rm{\Sigma }}}_{0}{a}_{0}^{2}{\tilde{a}}^{5/2-n}d\tilde{a}}{2\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}\\ & & \times {({\tilde{a}}^{-n+p+3/2}+t/{t}_{0})}^{-1}.\end{eqnarray} \tag{ A6 }$$

To relate this evolution to an observable quantity, we set the relative luminosity ${dL}(t)/{L}_{\star }={dA}(t)/4\pi {a}^{2}$ = ${dA}(t)/4\pi {\tilde{a}}^{2}{a}_{0}^{2}$. We then have

$$\begin{eqnarray}\displaystyle \frac{{dL}(t)}{{L}_{\star }} & = & \displaystyle \frac{3{x}_{m}{{\rm{\Sigma }}}_{0}{\tilde{a}}^{1/2-n}d\tilde{a}}{8\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}\\ & & \times {({\tilde{a}}^{-n+p+3/2}+t/{t}_{0})}^{-1}.\end{eqnarray} \tag{ A7 }$$

Integrating Equations (A6)–(A7) over $\tilde{a}$ yields the time evolution of the cross-sectional area and relative luminosity for material between ${\tilde{a}}_{\mathrm{in}}={a}_{\mathrm{in}}/{a}_{0}$ and ${\tilde{a}}_{\mathrm{out}}={a}_{\mathrm{out}}/{a}_{0}$.

Applying this theory to a real disk requires specifying parameters for the particles (${r}_{\min },$ ${r}_{\max }$, and ρ) and the disk (α₀, n, ${x}_{m}$, and p). Radiation pressure from the central star ejects particles with radii smaller than ${r}_{\min }$ (Burns et al. 1979); here we set ${r}_{\min }$ = 1 μm. For equal-mass objects with relative velocity v ≈ e v_K, the collision energy is Q_c = v²/8. The binding energy is

$$\begin{eqnarray}&&{Q}_{D}^{\star }={Q}_{b}{r}^{{\beta }_{b}}+{Q}_{g}{\rho }_{p}{r}^{{\beta }_{g}},\end{eqnarray} \tag{ A8 }$$

where Q_br^β_b is the bulk component of the binding energy and ${Q}_{g}{\rho }_{g}{r}^{{\beta }_{g}}$ is the gravity component of the binding energy (e.g., Benz & Asphaug 1999; Leinhardt et al. 2008; Leinhardt & Stewart 2009). Adopting parameters appropriate for rocky objects and setting Q_c ≈ ${Q}_{D}^{\star }$ yields (KB2016)

$$\begin{eqnarray}{r}_{c,\max } & \approx & 300{\left(\displaystyle \frac{e}{0.1}\right)}^{1.48}{\left(\displaystyle \frac{{v}_{K}}{30{{\rm{km}}{\rm{s}}}^{-1}}\right)}^{1.48}\\ & & \times {\left(\displaystyle \frac{\rho }{3{{\rm{g}}{\rm{cm}}}^{-3}}\right)}^{-0.74}{\left(\displaystyle \frac{{Q}_{g}}{{{\rm{0.3\; erg\; g}}}^{-1}}\right)}^{-0.74}\,\mathrm{km}.\end{eqnarray} \tag{ A9 }$$

Based on several simulations of planet formation in the terrestrial zone, e ≈ 0.1 is a reasonable choice (e.g., Weidenschilling et al. 1997; Kenyon & Bromley 2006; Raymond et al. 2007; Chambers 2013). Setting ${r}_{\max }$ = 300 km, the collision time at 1 au is

$$\begin{eqnarray}{t}_{0} & \approx & 2.4\times {10}^{5}\,{\alpha }_{0}\,{x}_{m}^{-1}\,\left(\displaystyle \frac{{r}_{\max }}{300\,\mathrm{km}}\right)\\ & & \times {\left(\displaystyle \frac{{{\rm{\Sigma }}}_{0}}{10{\rm{g}}{\mathrm{cm}}^{-2}}\right)}^{-1}\left(\displaystyle \frac{{P}_{0}}{1\,\mathrm{yr}}\right)\,\mathrm{yr}.\end{eqnarray} \tag{ A10 }$$

For a standard protostellar disk with ${{\rm{\Sigma }}}_{s}\propto {a}^{-p}$, we adopt p = 3/2 (e.g., Birnstiel et al. 2010; Windmark et al. 2012; Testi et al. 2014). In the standard analytic model, $\alpha \propto {({v}^{2}/{Q}_{D}^{\star })}^{-n}$ with n = 5/6 (Wyatt et al. 2007a, 2007b; Kobayashi & Tanaka 2010). However, comprehensive numerical simulations suggest n = 1 (KB2016). Here we consider α₀ = 1 and n = 0 (v = constant throughout the disk) or n = 1 (e = constant throughout the disk). Adopting $\tilde{t}$ = t/t₀, the luminosity evolution is then

$$\begin{eqnarray}\displaystyle \frac{L(t)}{{L}_{\star }} & = & \displaystyle \frac{{x}_{m}{{\rm{\Sigma }}}_{0}}{4\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}\\ & & \times \left[{\tan }^{-1}\left(\displaystyle \frac{{\tilde{a}}_{\mathrm{out}}^{3/2}}{\sqrt{\tilde{t}}}\right)-{\tan }^{-1}\left(\displaystyle \frac{{\tilde{a}}_{\mathrm{in}}^{3/2}}{\sqrt{\tilde{t}}}\right)\right]{\tilde{t}}^{-1/2}\end{eqnarray} \tag{ A11 }$$

for n = 0 and

$$\begin{eqnarray}&&\displaystyle \frac{L(t)}{{L}_{\star }}=\displaystyle \frac{3\sqrt{2}{x}_{m}{{\rm{\Sigma }}}_{0}}{32\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}{\tilde{t}}^{-3/4}{\rm{\Lambda }}(\tilde{a},\tilde{t})\end{eqnarray} \tag{ A12 }$$

for n = 1, where

$$\begin{eqnarray}{\rm{\Lambda }}(\tilde{a},\tilde{t}) & = & \left[\mathrm{ln}\left(\displaystyle \frac{\tilde{a}+{\tilde{t}}^{1/2}+\sqrt{2\tilde{a}{\tilde{t}}^{1/2}}}{\tilde{a}+{\tilde{t}}^{1/2}-\sqrt{2\tilde{a}{\tilde{t}}^{1/2}}}\right)\right.\\ & & +2{\tan }^{-1}(\sqrt{2\tilde{a}/{\tilde{t}}^{1/2}}+1)\\ & & +{2{\tan }^{-1}(\sqrt{2\tilde{a}/{\tilde{t}}^{1/2}}-1)]}_{{\tilde{a}}_{\mathrm{in}}}^{{\tilde{a}}_{\mathrm{out}}}.\end{eqnarray} \tag{ A13 }$$

At t = 0, solutions for any n yield the same initial luminosity:

$$\begin{eqnarray}&&\displaystyle \frac{{L}_{\mathrm{init}}}{{L}_{\star }}=\displaystyle \frac{{x}_{m}{{\rm{\Sigma }}}_{0}}{4\rho }{\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}({\tilde{a}}_{\mathrm{in}}^{-3/2}-{\tilde{a}}_{\mathrm{out}}^{-3/2}).\end{eqnarray} \tag{ A14 }$$

In terms of our adopted parameters:

$$\begin{eqnarray}\displaystyle \frac{{L}_{\mathrm{init}}}{{L}_{\star }} & \approx & 0.5{x}_{m}{\left(\displaystyle \frac{{r}_{\min }}{1\mu {\rm{m}}}\right)}^{-1/2}{\left(\displaystyle \frac{{r}_{\max }}{300\mathrm{km}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{a}_{\mathrm{in}}}{0.1\mathrm{au}}\right)}^{-3/2}\left[1-{\left(\displaystyle \frac{{a}_{\mathrm{in}}}{{a}_{\mathrm{out}}}\right)}^{3/2}\right].\end{eqnarray} \tag{ A15 }$$

When ${x}_{m}$ ≳ 2 × 10⁻⁴, the initial luminosity of the debris disk exceeds our nominal detection limit, ${L}_{d}/{L}_{\star }$ ≳ 10⁻⁴. At late times, we derive a single expression when n = 0 or n = 1:

$$\begin{eqnarray}\displaystyle \frac{{L}_{d}(t\gg {t}_{0})}{{L}_{\star }} & \simeq & \displaystyle \frac{(2n+1){x}_{m}{{\rm{\Sigma }}}_{0}}{4\rho }\\ & & \times ({\tilde{a}}_{\mathrm{out}}^{-n+3/2}-{\tilde{a}}_{\mathrm{in}}^{-n+3/2}){\left(\displaystyle \frac{1}{{r}_{\min }{r}_{\max }}\right)}^{1/2}{\tilde{t}}^{-1}.\end{eqnarray} \tag{ A16 }$$

As in the standard model for a single annulus, the luminosity declines linearly with time when t ≫ t₀. For our adopted parameters and n = 0,

$$\begin{eqnarray}\displaystyle \frac{{L}_{d}(t\gg {t}_{0})}{{L}_{\star }} & \approx & 0.04{x}_{m}{\left(\displaystyle \frac{{r}_{\min }}{1\mu {\rm{m}}}\right)}^{-1/2}{\left(\displaystyle \frac{{r}_{\max }}{300\mathrm{km}}\right)}^{-1/2}\\ & & \times {\left(\displaystyle \frac{{a}_{\mathrm{out}}}{2\mathrm{au}}\right)}^{-3/2}\left[1-{\left(\displaystyle \frac{{a}_{\mathrm{out}}}{{a}_{\mathrm{in}}}\right)}^{3/2}\right]{\tilde{t}}^{-1}.\end{eqnarray} \tag{ A17 }$$

With ${x}_{m}$ = 1 and t = 10 Myr, L_d/${L}_{\star }$ ≈ 10⁻³. The relative luminosity reaches the nominal detection limit of L_d/${L}_{\star }$ ≈ 10⁻⁴ at roughly 100 Myr.

Figure 7 illustrates results for a disk with a_in = 0.1 au and a_out = 2 au. The two sets of curves compare the analytic model with numerical results derived from dividing the disk into 101 annuli and summing dL(t)/${L}_{\star }$ (Equation (A7)) from each annulus. There is superb agreement between the two approaches. For any time, the difference between the analytic and numerical results for n = 0 or n = 1 is less than 0.1%. At 10 Myr, all of the model curves lie above the approximate detection limit of L_d/${L}_{\star }$ ≈ 10⁻⁴.

Zoom In Zoom Out Reset image size

To derive the predicted evolution of disk fluxes at specific wavelengths, we assume that all particles radiate as blackbodies at the equilibrium temperature, ${dF}(t)={dA}(t){B}_{\nu }(T)$, where $T=280\,{(a/{a}_{0})}^{-1/2}$ K. Although ${dF}(t)$ might be integrable, the result is very messy. Given the good agreement between the analytical and numerical results for L_d(t)/${L}_{\star }$, we integrate dF(t) numerically.

As outlined in the main text, drag within a remnant gaseous disk can rapidly remove small particles from a debris disk around 10 Myr old stars. To quantify this drift, we adopt a model gaseous disk (e.g., Adachi et al. 1976; Weidenschilling 1977b; Kenyon & Hartmann 1987; Chiang & Goldreich 1997; Takeuchi & Artymowicz 2001; Rafikov 2004; Chiang & Youdin 2010; Armitage 2013; Youdin & Kenyon 2013) with surface density, ${{\rm{\Sigma }}}_{g}={{\rm{\Sigma }}}_{g,0}{(a/1\mathrm{au})}^{-p}$, midplane temperature ${T}_{g}={T}_{g,0}{(a/1\mathrm{au})}^{-q}$, and vertical scale height $H/a={H}_{0}{(a/1\mathrm{au})}^{s}$. The midplane density is ${\rho }_{g}={{\rm{\Sigma }}}_{g}/2H$. Setting the sound speed ${c}_{s}^{2}=(\gamma {k}_{{\rm{B}}}{T}_{g}/\mu {{\rm{m}}}_{{\rm{H}}})$, where γ is the ratio of specific heats, k_B is Boltzmann's constant, μ is the mean molecular weight, and m_H is the mass of a hydrogen atom, the gas pressure in the midplane is ${P}_{g}=\rho {c}_{s}^{2}/\gamma$. The gas has mean free path $\lambda =\mu {{\rm{m}}}_{{\rm{H}}}{c}_{s}/{\rm{\Omega }}{{\rm{\Sigma }}}_{g}{\sigma }_{{H}_{2}}$ and viscosity $\nu =\lambda {v}_{t}{\rho }_{g}/3$, where ${\sigma }_{{H}_{2}}\simeq {10}^{-15}\,{{\rm{cm}}}_{2}$ is the collision cross section for H₂ and ${v}_{t}^{2}=(8{k}_{{\rm{B}}}{T}_{g}/\pi \mu {{\rm{m}}}_{{\rm{H}}})$ is the thermal velocity.

For the calculations in this paper we adopt p = 1, q = 0.5, H₀ = 0.03, s = 1/8, γ = 1.4, and μ = 2.0. Other choices for the exponents—p = 0.5–1.5, q = 0.3–0.7, and s = 0.0–0.3—have modest impact on the results. In all configurations, the ratio of the mass density of solids to the mass density of gas is small.

Particles with radii r ≤ 9λ/4 are in the Epstein regime, where the drag force is ${F}_{d}=\pi {r}^{2}\rho {v}_{d}{({v}_{t}^{2}+{v}_{d}^{2})}^{1/2}$. As in Takeuchi & Artymowicz (2001), the term ${({v}_{t}^{2}+{v}_{d}^{2})}^{1/2}$ provides an estimate of the drag force in the subsonic and supersonic regimes. For larger particles, all drift is in the subsonic regime, where the drag force is ${F}_{d}={C}_{d}\pi {r}^{2}\rho {v}_{d}^{2}/2$ and C_d is a drag coefficient that depends on the Reynolds number $\mathrm{Re}\,=\,2\rho {v}_{d}r/\nu$ (e.g., Equations (8a)–(8c) in Weidenschilling 1977b).

The main parameters in this model are the gravity of the central star $g={{GM}}_{\star }/{a}^{2}$, the residual gravity ${\rm{\Delta }}g={\rho }_{g}^{-1}{dP}/{da}$, and the stopping time ${t}_{s}={{mv}}_{d}/{F}_{d}$ (Adachi et al. 1976; Weidenschilling 1977b; Takeuchi & Artymowicz 2001). The gas orbits at a slightly smaller velocity than the local circular velocity. Defining $\eta =-{\rm{\Delta }}g/g$, the gas velocity is

$$\begin{eqnarray}&&{v}_{g}={v}_{K}{(1-\eta )}^{1/2}.\end{eqnarray} \tag{ A18 }$$

The maximum velocity of drifting particles is roughly

$$\begin{eqnarray}&&{v}_{d,\max }\approx \eta {v}_{K}/2.\end{eqnarray} \tag{ A19 }$$

In an optically thin disk, radiation pressure drives small particles radially outward (e.g., Burns et al. 1979; Takeuchi & Artymowicz 2001). Defining β as the ratio of the radiation pressure force to the gravitational force,

$$\begin{eqnarray}&&\beta =\displaystyle \frac{3{L}_{\star }{Q}_{{pr}}}{16\pi {{cGM}}_{\star }r{\rho }_{s}},\end{eqnarray} \tag{ A20 }$$

where c is the speed of light. Although Takeuchi & Artymowicz (2001) also consider the impact of Poynting–Robertson drag, this drag has little impact when the collision time is short. Thus, we ignore Poynting–Robertson drag in this discussion.

Weidenschilling (1977b) and Takeuchi & Artymowicz (2001) divide the drift velocity of solid particles relative to the gas into a radial component ${v}_{\mathrm{rad}}$ and an azimuthal component v_θ, with ${v}_{d}^{2}={v}_{\mathrm{rad}}^{2}+{v}_{\theta }^{2}$. Defining the angular velocity ${\rm{\Omega }}={({{GM}}_{\star }/{a}^{3})}^{1/2}$ and setting τ_s = t_sΩ, the two components of the drift are

$$\begin{eqnarray}&&{v}_{\mathrm{rad}}=(\eta -\beta ){v}_{K}/({\tau }_{s}+{\tau }_{s}^{-1})\end{eqnarray} \tag{ A21 }$$

and

$$\begin{eqnarray}&&{v}_{\theta }={\tau }_{s}{v}_{\mathrm{rad}}/2.\end{eqnarray} \tag{ A22 }$$

Three parameters—β, η, and τ_s—divide drift into four major regimes. Large particles poorly coupled to the gas (τ_s ≫ 1) follow circular orbits with ${v}_{\mathrm{rad}}$ ≈ 0 and v_θ ≈ η v_K. Intermediate-size particles weakly coupled to the gas (τ_s ≈ 1) have maximal drift velocities and somewhat smaller angular velocities relative to the gas. Small particles (τ_s ≪ 1) are entrained in the gas. When η ≳ β (η ≲ β) for small objects, ${v}_{\mathrm{rad}}$ is positive (negative); particles drift inward (outward).

To identify these regimes in gaseous disks at 1 au, we develop an iterative technique to solve Equations (A21)–(A22). For each particle size, the algorithm adopts a drag regime, derives t_s and τ_s, infers ${v}_{\mathrm{rad}}$ and v_θ, and then verifies the drag regime. This process repeats until the final drag regime is identical to the initial drag regime. When particles are small enough to experience Epstein drag, the algorithm loops through an initial iteration to derive a consistent v_d for the subsonic and supersonic regimes. Comparisons with results in Weidenschilling (1977b) and Takeuchi & Artymowicz (2001) confirm the accuracy of our approach.

Figure 8 illustrates the variation in drift velocity with semimajor axis for particles with sizes of 1–100 μm in a gaseous disk with a surface density of 0.001% of the MMSN. Close to the star, radiation pressure drives particles to larger a (open circles). Far from the star, the gas drags particles inward (filled circles). At some intermediate a, small particles find an equilibrium distance a_eq where radiation pressure balances gas drag. For the disk model adopted here, this equilibrium is at

$$\begin{eqnarray}&&{a}_{\mathrm{eq}}\approx \displaystyle \frac{40\,\mu {\rm{m}}}{r}\,\mathrm{au}.\end{eqnarray} \tag{ A23 }$$

Zoom In Zoom Out Reset image size

In solar-type stars with ${r}_{\min }$ ≈ 1 μm, small particles with r ≈ 1–10 μm produce most of the IR excess. When the disk's surface density is 0.001% of the MMSN, these small particles have a_eq ≳ 4 au and T_d ≲ 140 K. These particles are then too cold to produce a warm debris disk.