Variations on Debris Disks III. Collisional Cascades and Giant Impacts in the Terrestrial Zones of Solar-type Stars
aa r X i v : . [ a s t r o - ph . E P ] D ec Variations on Debris Disks III. Collisional Cascades and GiantImpacts in the Terrestrial Zones of Solar-type Stars
Scott J. Kenyon
Smithsonian Astrophysical Observatory, 60 Garden Street, Cambridge, MA 02138 e-mail: [email protected]
Benjamin C. Bromley
Department of Physics, University of Utah, 201 JFB, Salt Lake City, UT 84112 e-mail: [email protected]
ABSTRACT
We analyze two new sets of coagulation calculations for solid particles orbitingwithin the terrestrial zone of a solar-type star. In models of collisional cascades,numerical simulations demonstrate that the total mass, the mass in 1 mm andsmaller particles, and the dust luminosity decline with time more rapidly thanpredicted by analytic models, ∝ t − n with n ≈ r . n ≈ ∼ a . &
30 Myr old stars agrees with our calculations,the observed frequency of warm debris among 5–20 Myr old stars is smaller thanpredicted.
Subject headings: planetary systems – planets and satellites: formation – proto-planetary disks – stars: formation – zodiacal dust – circumstellar matter 2 –
1. INTRODUCTION
Dense disks of debris surround many main sequence stars (e.g., Backman & Paresce1993; Wyatt 2008; Matthews et al. 2014). In these systems, particles with radii r ≈ µ mto 1 cm and temperatures T ≈ ∼ ∼
50% for ages of ∼
10 Myr to .
5% at 0.5–1 Gyr (e.g., Rieke et al. 2005; Currieet al. 2008; Carpenter et al. 2009a). FGK stars have somewhat lower debris disk frequenciesof 10% to 20%, with excesses that decline more slowly with stellar age (Carpenter et al.2009a,b; Kennedy & Wyatt 2013; Matthews et al. 2014). Several older stars also appear tohave substantial and very luminous disks or rings of debris which may decline rapidly withtime (e.g., Song et al. 2005; Rhee et al. 2008; Melis et al. 2010, 2012; Meng et al. 2015).Debris disks are signposts of planet formation (e.g., Kenyon & Bromley 2002b; Zucker-man & Song 2004; Kenyon & Bromley 2008; Raymond et al. 2011, 2012). In current theory,protoplanets grow through mergers of smaller objects. Once they reach sizes of 500–1000 km,protoplanets stir nearby small objects to large velocities. Continued growth enables stirringover larger and larger volumes, generating a cascade of destructive collisions among smallleftovers. This collisional cascade grinds the leftovers into a ring or disk of very small par-ticles with luminosity comparable to the dust luminosities of many debris disks. As thecascade evolves, occasional ‘giant impacts’ between pairs of protoplanets produce clouds ofsmall particles which interact with the rest of the disk (e.g., Agnor et al. 1999; Asphauget al. 2006; Grigorieva et al. 2007; Genda et al. 2012). In some giant impact models, thedust luminosity roughly matches the very luminous disks associated with old G-type stars.Interpreting observations of debris disks requires a robust theory of planet formationwhich predicts the time evolution of luminosity and the structure of solid material from ∼ r max orbiting the central star at high velocity with eccentricity e , inclination i and a range of semimajor axes a (e.g., Wyatt & Dent 2002; Dominik & Decin 2003; Krivovet al. 2006; Wyatt et al. 2007a,b; L¨ohne et al. 2008; Kobayashi & Tanaka 2010; Kennedy& Wyatt 2010; Wyatt et al. 2011). For an initial surface density Σ , several simplifyingassumptions allow these models to predict the time evolution of the size distribution andtotal mass of the debris. Successful applications include matches to the time evolution ofdebris disks around A-type and G-type stars (Matthews et al. 2014). However, the approachfails to account for the frequency of giant impacts at late times.Several numerical approaches place debris disks in the context of the long-term evolution 3 –of a primordial disk of gas and dust (Kenyon & Bromley 2005, 2008, 2010; Weidenschilling2010a; Raymond et al. 2011, 2012; Kobayashi & L¨ohne 2014). These studies adopt aninitial mass and radial structure for the disk, establish prescriptions for the evolution andinteractions of the gas and the solids, and derive a , e , i , Σ, and other properties of the debrisas a function of the initial conditions and the final architecture of the planetary system. Withan accurate treatment for the long-term evolution of the cascade and the ability to followdust production from giant impacts at late times, these calculations can match observationsof the time evolution of the brightness and frequency of debris disks. However, the cpuintensive nature of this approach often limits their ability to identify how the evolutiondepends on uncertain input parameters.Other numerical approaches extract plausible snapshots from detailed planet formationsimulations and model a well-defined subset of the system (Grigorieva et al. 2007; Boothet al. 2009; Mustill & Wyatt 2009; Stark & Kuchner 2009; Jackson & Wyatt 2012; Nesvoldet al. 2013; Kral et al. 2013). Using n -body dynamical calculations combined with semi-analytic results, these investigations derive the geometry and luminosity of the debris as afunction of time. Aside from learning how these features depend on initial conditions, theseinvestigations excel at matching observations of particular systems. These algorithms arealso very cpu intensive and hinge on specifying a fairly small set of initial conditions from avery broad range of possibilities.Weaving these investigations into a more robust theory for planet formation requirescalculations which test and unify the different approaches. Here we consider results fromtwo sets of coagulation calculations for swarms of particles orbiting at 1 AU from a solar-type star. The first set of simulations enables robust tests of analytic models for collisionalcascades. The second set allows us to place collisional cascades in context with more detailedplanet formation simulations. Comparisons between the two sets yield common features anddifferences which are useful for interpreting existing observations.Our study begins in § § § § § § §
2, glance through the figures in §
4, read the discussionin §
5, and finish with the summaries in § §
6. 4 –
2. ANALYTIC MODEL
Building on the collisional cascade models of Dohnanyi (1969) and Hellyer (1970), Wyatt& Dent (2002) and Dominik & Decin (2003) developed analytic models for the long-termcollisional evolution of particles in a debris disk. Objects with radius r , mass m , and massdensity ρ orbit within a cylindrical annulus with width δa centered at a distance a from acentral star with mass M ⋆ and luminosity L ⋆ . All particles have the same orbital eccentricity e and inclination i . Destructive collisions generate a collisional cascade where the mass insolids gradually diminishes with time. For an upper mass limit m max of the largest solidobjects participating in the cascade, the analytic model yields a simple formula for N max ( t ),the number of these large objects as a function of time. Radiation pressure establishes m min ,a lower mass limit for solids with stable orbits around the central star. Between m min and m max , the cascade produces a power-law size distribution, N ( r ) ∝ r − q (e.g., Dohnanyi 1969;Hellyer 1970; Williams & Wetherill 1994; O’Brien & Greenberg 2003; Kobayashi & Tanaka2010; Wyatt et al. 2011). Setting the slope q of this size distribution yields another simpleformula for the time evolution of the dust luminosity L d ( t ). To derive expressions for N max ( t ) and L d ( t ), we adopt a particle-in-a-box model, wherekinetic theory sets the collision rate (Wyatt & Dent 2002; Dominik & Decin 2003). For anannulus with volume V , we consider an ensemble of mono-disperse particles with numberdensity n max = N max /V , collision cross-section σ = 4 πr max , and relative velocity v . Everycollision destroys two particles; the loss rate is ˙ N max = − N max n max σv .To express this rate in terms of basic disk parameters, we define the initial number ofparticles N ,max , initial surface density Σ = N ,max m max / πaδa , the vertical height of theannulus H , and the volume V = 4 πaδaH . For solids orbiting with angular frequency Ω,orbital period P = 2 π/ Ω, vertical velocity v z , and Keplerian velocity v K within the annulus, H ≈ v z / Ω, i ≈ e/ v z ≈ ev K /
2, and v ≈ ev K (e.g., Wetherill & Stewart 1993; Kenyon &Luu 1998).We set the collision time as t c = r max ρP π Σ . (1)In this approach, t c is the time to destroy one of the largest particles in the ensemble.Formally, collisions destroy two particles on the time scale 2 t c . The collision rate thenbecomes ˙ N max = − N max /N ,max t c , (2) 5 –which has a simple solution N max ( t ) = N ,max t/t c . (3)In this derivation, the collision time depends only on properties of the largest objects inthe swarm (see also Dominik & Decin 2003). In a real cascade, collisions with smaller particlesin the swarm also remove mass from the largest particles (e.g., Wyatt & Dent 2002; Dominik& Decin 2003; Wyatt et al. 2007a,b; Kobayashi & Tanaka 2010; Kennedy & Wyatt 2011;Kennedy et al. 2011; Wyatt et al. 2011). In addition to completely destructive collisions withobjects having m & m max , cratering collisions with much smaller particles graduallychip away material from the largest objects. These ‘extra’ collisional processes modify thecollision time by a factor α . α for an equilibrium cascade, we adopt results from Kobayashi & Tanaka(2010, and references therein). If Q ⋆D is the collision energy required to eject half the massfrom a pair of colliding objects and if Q c is the center-of-mass collision energy, the largestobject in the debris has mass m l ≈ . Q c /Q ⋆D ) − . For simplicity, we assume q ≈ s = 6 . s = 38 .
1, and s = 5 . α ≈ . (cid:18) v Q ⋆D (cid:19) − / . (4)For ( v / Q ⋆D ) ≫
1, eqs. 9–13 of Wyatt et al. (2007a) yield similar results. In most collisionalcascades, ( v / Q ⋆D ) &
4; thus, α . ρ = 3 g cm − in a disk with Σ = 10 g cm − (comparableto the minimum mass solar nebula; e.g., Weidenschilling 1977; Hayashi 1981; Chiang &Youdin 2010), t c ≈ × α (cid:16) r max
100 km (cid:17) (cid:18) Σ
10 g cm − (cid:19) − (cid:18) P (cid:19) yr . (5)Collision times are short, ∼ r max . With N ( t ) known, the dust luminosity follows. Destructive collisions of the largestobjects produce a power-law size distribution extending from the largest objects with mass 6 – m max to the smallest objects with mass m min . The total mass M d and cross-sectional area A d of particles within the annulus are simple functions of the minimum size, the maximumsize, and the slope q (e.g., Wyatt & Dent 2002; Dominik & Decin 2003; Kenyon et al. 2014).The stellar energy intercepted by the solids is L d = A d / πa . If m min , m max , and q neverchange, the initial dust luminosity L d, is a simple function of N max, and these parameters.The time evolution of the dust luminosity is then: L d ( t ) = L d, t/t c . (6)Setting r min = 1 µ m and q = 3.5 for the size distribution for an annulus with δa = 0 . a at1 AU: L d, = 7 . × − (cid:18) ρ − (cid:19) − (cid:16) r max
100 km (cid:17) − / (cid:18) r min µ m (cid:19) − / (cid:18) Σ
10 g cm − (cid:19) L ⋆ . (7)The initial dust luminosity is independent of the collision time.At late times ( t ≫ t c ), L d ( t ) ≈ L d, t c /t . Combining L d, and t c into a new normalizationfactor L ′ d, , L d ( t ) = L ′ d, (cid:18) t (cid:19) − , (8)where L ′ d, = 5 . × − α (cid:16) r max
100 km (cid:17) / (cid:18) r min µ m (cid:19) − / (cid:18) P (cid:19) L ⋆ . (9)Despite starting out with much smaller initial dust luminosity, old cascades with large r max are easier to detect than those with small r max (Wyatt 2008; Kobayashi & Tanaka 2010). Atlate times, the ability to detect a cascade is independent of the initial mass involved in thecascade. Through the parameter α , the dynamics of the swarm also set the detectability. r min and r max Maintaining a power-law size distribution for particles with m min . m . m max requiresdestructive collisions among roughly equal mass objects . For impact velocity v , Q c = µv / m + m ), where µ = m m / ( m + m ) is the reduced mass for a pair of colliding Although cratering can also drive a cascade, the largest objects then probably grow with time. Forsimplicity, we focus here on systems where collisions reduce the masses of all objects. We return to this issuein § m and m . For equal mass objects, Q c = v / Q ⋆D = Q b r β b + Q g ρ p r β g (10)where Q b r β b is the bulk component of the binding energy and Q g ρ g r β g is the gravity com-ponent of the binding energy (e.g., Benz & Asphaug 1999; Leinhardt et al. 2008; Leinhardt& Stewart 2009). For rocky objects, Q b ≈ × erg g − cm − β b , β b ≈ − . Q g ≈ − cm − β g , and β g ≈ Q c ≈ Q ⋆D establishes constraints on the particles destroyed by an adopted col-lision velocity. For large objects with v ≈ ev K and β g = 1.35, collisions destroy equal-massparticles with r . r c,max , where r c,max ≈ (cid:16) e . (cid:17) . (cid:16) v K
30 km s − (cid:17) . (cid:18) ρ − (cid:19) − . (cid:18) Q g . − (cid:19) − . km . (11)For small objects, β b = − .
4; collisions destroy equal mass objects with r & r c,min , where r c,min ≈ . (cid:16) e . (cid:17) − . (cid:16) v K
30 km s − (cid:17) − . (cid:18) Q b erg g − (cid:19) . µ m . (12)The first relation establishes that objects with r .
300 km participate in the cascade.Because radiation pressure typically ejects dust grains with r . µ m, the second relationconfirms that all objects smaller than 300 km participate in the cascade. Once the cascadebegins, collisions maintain it forever.Although collisions between equal-mass objects with r > r c,max result in some net growthof one particle, cratering collisions with the rest of the swarm can gradually reduce the massof these objects (e.g., Kobayashi & Tanaka 2010). As with estimates of the α factor for thecollision time (eq. 4), deriving the maximum size of particles which slowly lose mass requiresintegrating the collision rate over the mass distribution. When we examine results of ournumerical simulations, we will illustrate how this maximum size depends on the propertiesof the swarm.Once the cascade begins, the slope q of the power-law size distribution depends on thedetails of the Q ⋆D relation (O’Brien & Greenberg 2003; Kobayashi & Tanaka 2010; Wyatt 8 –et al. 2011). When Q ⋆D is independent of r , q = 3.5. For the standard Q ⋆D ( r ) in eq. 10, thesize distribution consists of two power laws, with q s = (7 − β b / / (2 − β b / ≈ q l = (7 + β g / / (2 + β g / ≈ Q ⋆D ( r ) reaches a minimum, r Q ≈ Solving eqs. 3–6 for the long-term evolution of the mass in solids and the dust luminosityrequires initial conditions for the central star ( M ⋆ , L ⋆ ), the geometry of the annulus ( a , δa , e ,and i ), and the properties of the solid particles ( r min , r max , N ,max , q , and Q ⋆D ). To illustratethe long-term evolution, we consider a swarm of solid material orbiting a 1 M ⊙ star at 1 AU.Fig. 1 shows the time evolution of the dust luminosity for particles with r min = 1 µ m,various r max Σ = 10 g cm − , e = 0.1, i = e /2, and α = 1 orbiting in an annulus with δa =0.2 at a = 1 AU around a central star with M ⋆ = 1 M ⊙ , L ⋆ = 1 L ⊙ . When r max is 10 km, theinitial luminosity is large, but the collision time is short. Within 1 Myr, the luminosity fallsby a factor of ∼
50. After ∼
30 Myr, the relative luminosity falls below 10 − . When r max ismuch larger, longer collision times allow the luminosity to remain large at later times. For r max = 1000 km, it takes 30–40 Myr (300–400 Myr) for the relative luminosity to fall below10 − (10 − ).
3. NUMERICAL MODEL
To perform numerical calculations of collisional cascades, we use
Orchestra , an ensembleof computer codes for the formation and evolution of planetary systems. In addition toother algorithms,
Orchestra includes a multiannulus coagulation code which derives the timeevolution of a swarm of solid particles orbiting a central object (Kenyon & Bromley 2004a,2008, 2012). Here, we use the coagulation code to follow the evolution of solids at 1 AUaround a solar mass star.
To provide the closest representation of the analytic model, we conduct coagulationcalculations within a single annulus with semimajor axis a and width δa around a starwith mass M ⋆ = 1 M ⊙ . Within this annulus, there are M mass batches with characteristic 9 –mass m k and radius r k (Spaute et al. 1991; Wetherill & Stewart 1993; Weidenschilling et al.1997; Kenyon & Luu 1998). Batches are logarithmically spaced in mass, with mass ratio δ ≡ m k +1 /m k . Each mass batch contains N k particles with total mass M k and average mass¯ m k = M k /N k . Particle numbers N k < are always integers. Throughout the calculation,the average mass is used to calculate the average physical radius ¯ r k , collision cross-section,collision energy, and other necessary physical variables. As mass is added and removed fromeach batch, the average mass changes (Wetherill & Stewart 1993).For any δ , numerical calculations lag the result of an ideal calculation with infinite massresolution (see § A.1). In most cases, simulations with δ = 1.05–1.19 yield better solutionsto the evolution of the largest objects ( r & δ = 1.41–2.00.Although simulations with δ = 1.05–1.10 allow better tracking of the gradual reduction inmass of the largest objects during the cascade, the evolution of the size distribution and thedust luminosity is fairly independent of δ . Thus, we consider a suite of calculations with δ = 1.19 (= 2 / ).In this suite of calculations, we follow particles with sizes ranging from a minimumsize r min = 1 µ m to the maximum size r max . The algorithm for assigning material to themass bins extends the maximum size as needed to accommodate the largest particles. Whencollisions produce objects with radii r < r min , this material is lost to the grid. All calculations begin with a swarm of particles with initial maximum size r and massdensity ρ p = 3 g cm − . These particles have initial surface density Σ , total mass M , andhorizontal and vertical velocities v h, and v z, relative to a circular orbit. The horizontalvelocity is related to the orbital eccentricity, e = 1.6 ( v h /v K ) , where v K is the circularorbital velocity. The orbital inclination is sin i = √ v z /v K .For the simulations in this paper, we consider two different initial size distributions forsolid objects. To follow the analytic model as closely as possible, one set of calculationsbegins with a power law size distribution, N ( r ) ∝ r − q and q = 3.5. To study whether ourcalculations produce this equilibrium size distribution, we begin a second set of calculationswith a mono-disperse set of particles. 10 – The mass and velocity distributions of the solids evolve in time due to inelastic collisions,drag forces, and gravitational encounters. As summarized in Kenyon & Bromley (2004a,2008), we solve a coupled set of coagulation equations which treats the outcomes of mutualcollisions between all particles in all mass bins. We adopt the particle-in-a-box algorithm,where the physical collision rate is nσvf g , n is the number density of objects, σ is thegeometric cross-section, v is the relative velocity, and f g is the gravitational focusing factor(Wetherill & Stewart 1993; Kenyon & Luu 1998). The collision algorithm treats collisions inthe dispersion regime – where relative velocities are large – and in the shear regime – whererelative velocities are small (Kenyon & Luu 1998; Kenyon & Bromley 2014). Within anannulus, massive protoplanets on nearly circular orbits are ‘isolated’; these objects can collidewith smaller objects but cannot collide with other isolated objects (Wetherill & Stewart 1993;Kenyon & Luu 1998; Kenyon & Bromley 2015).In every time step, our algorithm requires an integral number of collisions N kl,int betweenmass bins k and l . For each derived N kl , the algorithm compares the fractional part to arandom number drawn for each ( k, l ) pair. When the random number exceeds the fractionalpart, N kl is rounded down; otherwise N kl is rounded up. Comparisons with a Runge-Kuttacode (which allows fractional particles) and tests on systems with analytic solutions verifythis approach (e.g., Kenyon & Luu 1998; Kenyon & Bromley 2015).For these simulations, we consider two approaches to collision outcomes. Adopting Q ⋆D = constant allows us to make the most direct comparison with analytic models. Calcula-tions with the more standard Q ⋆D ( r ) measure the sensitivity of the evolution to changes in Q ⋆D and enable links to other published numerical simulations of collisional cascades. Forthe rocky solids in this study, we adopt Q b = 3 × erg g − cm − β b , β b = − . Q g =0.3 erg g − cm − β g , and β g = 1.35. These parameters are broadly consistent with publishedanalytic calculations and numerical simulations (e.g., Davis et al. 1985; Holsapple 1994; Love& Ahrens 1996; Housen & Holsapple 1999). At small sizes, they agree with results from lab-oratory experiments (e.g., Ryan et al. 1999; Arakawa et al. 2002; Giblin et al. 2004; Burchellet al. 2005).For each pair of colliding planetesimals, the mass of the merged planetesimal is m = m + m − m esc , (13)where the mass of debris ejected in a collision is m esc = 0 . m + m ) (cid:18) Q c Q ∗ D (cid:19) b d . (14) 11 –The exponent b d is a constant of order unity (e.g., Davis et al. 1985; Wetherill & Stewart1993; Kenyon & Luu 1999; Benz & Asphaug 1999; O’Brien & Greenberg 2003; Kobayashi &Tanaka 2010; Leinhardt & Stewart 2012). Here, we consider b d = 1 and 9/8.To place the debris in the grid of mass bins, we set the mass of the largest collisionfragment as m max,d = m l, (cid:18) Q c Q ∗ D (cid:19) − b l m esc , (15)where m l, ≈ b l ≈ b l is large, catastrophic collisionswith Q c & Q ⋆D crush solids into smaller fragments. Lower mass objects have a differential sizedistribution N ( r ) ∝ r − q d . After placing a single object with mass m max,d in an appropriatebin, we place material in successively smaller mass bins until (i) the mass is exhausted or(ii) mass is placed in the smallest mass bin. Any material left over is removed from the grid.For calculations of collisional cascades, we assume that the orbital e and i are constantwith time. When we consider how the growth of solids leads to a collisional cascade, wederive orbital evolution due to collisional damping from inelastic collisions and gravitationalinteractions. For inelastic and elastic collisions, we follow the statistical, Fokker-Planckapproaches of Ohtsuki (1992) and Ohtsuki et al. (2002), which treat pairwise interactions(e.g., dynamical friction and viscous stirring) between all objects. We also compute long-range stirring from distant oligarchs (Weidenschilling 1989).Our solutions to the evolution equations conserve mass and energy to machine accuracy.Typical collisional cascade calculations require a single run on a system with 16–32 cpus;over the 10 timesteps in a typical 2 Gyr run, calculations conserve mass and energy to betterthan a part in 10 . Planet formation calculations usually require 3–4 runs on 64 cpus. Over10 − timesteps, our calculations conserve mass and energy to better than a part in 10 . Before describing our calculations, it is useful to place our numerical approach in contextwith other simulations. Starting with Safronov (1969), kinetic approximations to the particlecollision rate and Fokker-Planck estimates of velocity evolution have a long history (e.g.,Bandermann 1972; Greenberg et al. 1978; Wetherill 1980; Weidenschilling 1980; Greenberget al. 1984; Ohtsuki & Nakagawa 1988; Wetherill & Stewart 1989; Spaute et al. 1991; Barge& Pellat 1991; Wetherill & Stewart 1993; Kokubo & Ida 1996; Weidenschilling et al. 1997;Grogan et al. 2001; Ohtsuki et al. 2002; Kokubo et al. 2006; Chambers 2008; Morbidelli et al. 12 –2009; Kobayashi & Dauphas 2013; Glaschke et al. 2014; Shannon et al. 2015). In addition toour studies (e.g., Kenyon & Bromley 2002b, 2004a, 2008, 2010), Kobayashi & Tanaka (2010),Weidenschilling (2010b), and Kobayashi & L¨ohne (2014) follow the evolution of debris disksusing a similar set of algorithms. Aside from various details of the implementation – e.g.,single annulus (Kobayashi & Tanaka 2010) or multi-annulus (Weidenschilling 2010b) – thesecoagulation codes are similar to our multi-annulus coagulation code.Several other techniques adopt the kinetic approximation for collisions but do not allowvelocity evolution from gravitational interactions between mass bins. In addition to followingthe evolution of particle number in specific mass bins, Krivov et al. (2006) track bins insemimajor axis and orbital eccentricity (see also L¨ohne et al. 2008, 2012; Krivov et al.2013). In these models, radiation pressure modifies the orbital parameters of grains ejectedfrom collisions of larger particles (e.g., Burns et al. 1979). As an interesting variant of theparticle-in-a-box model with discrete mass bins and integral numbers of particles, G´asp´aret al. (2012a) derive the time evolution of the differential number density in a single annuluswith collision rates inferred from the kinetic approximation. Although this technique doesnot track variations in a and e , G´asp´ar et al. (2012a) derive the impact of radiation pressureas in Krivov et al. (2006).A third group of codes employs n -body techniques to evolve a collection of masslesstracer particles (e.g., Grigorieva et al. 2007; Booth et al. 2009; Stark & Kuchner 2009; Jackson& Wyatt 2012; Nesvold et al. 2013; Kral et al. 2013). These codes solve the equations ofmotion for a swarm of massless tracers responding to the gravity and radiation pressure froma central star and the gravity of a few other massive planets orbiting the central star. Thesetracers serve as proxies for a mass distribution of particles; the surface density of tracers isinput for estimating the collision rates between particles of different masses. Although thealgorithms in each implementation are different, aspects of these codes are similar to n -body-based (e.g., Raymond et al. 2011, 2012; Levison et al. 2012) and hybrid (Weidenschilling et al.1997; Bromley & Kenyon 2006, 2011a,b, 2013; Glaschke et al. 2014) codes for calculating theformation and evolution of planetary systems.In this paper, we compare results from single annulus coagulation calculations to pre-dictions of the analytic model (e.g., Wyatt et al. 2011) and (where possible) other numericalcalculations using the kinetic approximation (e.g., Krivov et al. 2006; Kobayashi & Tanaka2010; G´asp´ar et al. 2012a). Exploring outcomes as a function of Q ⋆D , b d , m l, , b l and q d allowus to establish the sensitivity of observable quantities on unknown parameters. These resultsserve as a foundation for future multi-annulus calculations with tracer particles, which willyield comparisons with n -body approaches using tracer super-particles (e.g., Grigorieva et al.2007; Stark & Kuchner 2009; Nesvold et al. 2013; Kral et al. 2013). 13 –
4. NUMERICAL CALCULATIONS4.1. Collisional Cascades Q ⋆D To explore the long-term evolution of collisional cascades, we consider a ‘standard’calculation of solids orbiting a solar-type star. Material with maximum initial radius r =100 km, r min = 1 µ m, and ρ = 3 g cm − orbits with e = 0.1 and i = 0.05 ( ≈ ◦ ) inside anannulus with a = 1 AU and δa = 0.2 AU. The initial size distribution is either a power lawwith q = 3.5 or a mono-disperse ensemble with r = r . All objects have a constant Q ⋆D =6 × erg g − . To assign collisional debris to mass bins, we set q d = 3.5 and examine resultsfor various combinations of b d , m l, , and b l . All calculations ignore Poynting-Robertson drag.For all calculations with constant e and i , the size of the largest object gradually declineswith time (Fig. 2). The timing of this decline depends on the initial size distribution. Whenthe swarm is initially mono-disperse (upper panel), all collisions are catastrophic. As thesefairly infrequent collisions produce smaller objects, cratering begins to chip away at themass in all of the large objects. When the initial size distribution of the swarm is a power-law (lower panel), cratering collisions begin immediately. Thus, r max initially declines morerapidly for calculations with an initial power law size distribution.Once r max begins to decline, the evolution is more dramatic in originally mono-disperseswarms. During the first 0.01–0.1 Myr of evolution, these swarms produce fewer particleswith sizes smaller than 1 µ m which are lost to the grid. The swarms then have relativelymore mass than swarms with an initial power law size distribution. More massive swarmshave higher collision rates, resulting in a faster decline of r max with time.The time evolution of r max also depends on the way we distribute debris into the massbins. When b d = 1 (9/8), cratering (catastrophic) collisions with Q c . Q c &
1) producerelatively more debris. For swarms with a power law size distribution from 1 µ m to 100 km,cratering collisions produce more debris than catastrophic collisions. Thus, calculations with b d = 1 reduce r max more rapidly than those with b d = 9/8. For any b d , the mass of the largestobject in the debris is smaller in calculations with b l > b l = 0. Thus, r max declines more rapidly when b l is larger than zero.For any combination of b d , m l, , and b l , the dramatic decline in r max roughly follows apower-law. From ∼ r max ∝ t − n with n = 0.1–0.2. After 100 Myr, particleswith r max = 100 km have lost from 97% (power law size distribution, b l = 0) to 99.9%(mono-disperse, b l = 1) of their initial mass. 14 –To describe the evolution of the size distribution, we focus first on results for b d = 1 and b l = 0 (Fig. 3–4). Based on previous analyses, we expect a roughly power law size distributionwith q = 3.5 from 1 µ m to r max . To minimize shot noise, we derive the cumulative sizedistribution, N ( > r ), where the predicted slope is q c = 2.5. To illustrate deviations fromthis prediction, we examine the relative cumulative size distribution, N ( > r ) /r − . . Theexpected power law slope of this relative size distribution is then q r = q c − . r < µ m fromthe grid. Compared to a model which includes these very small particles, grains with r ≈ µ m experience fewer destructive collisions, leading to a clear excess of these particles.With an excess of 1–10 µ m grains, particles with r ≈ µ m experience more destructivecollisions, producing a deficit of these particles (e.g., O’Brien & Greenberg 2003; Kral et al.2013).Collisional evolution also leads to excesses and deficits at r & r ≈ r max ,destructive collisions gradually shift mass into lower mass bins. Without a correspondingaddition in mass from bins with r & r max , there is a deficit in the mass distribution at r ≈ r max . This deficit enables an excess at smaller sizes, r ≈
10 km.Among intermediate sized particles (100 µ m to 1–3 km), there is a small wave aboutthe expected power law. This wave is a function of the deficits and excesses of particlesat smaller and larger sizes (see also § A.2; Campo Bagatin et al. 1994; Durda et al. 1998;O’Brien & Greenberg 2003). As cratering reduces the sizes of the largest particles, the peakof the wave at 5–20 km gradually shifts to smaller sizes. Over 100 Myr, the shape of thewavy power law remains roughly constant at smaller sizes.For swarms starting with a mono-disperse size distribution, the system rapidly evolvesto the ‘standard’ size distribution (Fig. 4). After ∼ µ m to 10 km. At 0.1–100 Myr, there is a clear excess of 10–100 km objects whichis not visible in calculations starting with a power law size distribution of small objects. Inmono-disperse calculations, infrequent collisions among the largest objects prevents reachingthe equilibrium size distribution.Overall, all approaches to fragmentation produce wavy structures in the relative cumu-lative size distribution (Fig. 5). However, the amplitudes and positions of the waves dependon b d and b l . When b d = 9/8 (instead of 1), the deficits at 10 µ m and 10–30 km are pro-nounced but somewhat smaller. Adopting b l = 0.75–1.0 reduces the amplitude of the wave 15 –at small sizes but accentuates it at larger sizes.To conclude this discussion of the standard model, Fig. 6 shows the time evolution ofthe dust luminosity L d . When collisional cascades begin with a mono-disperse swarm oflarge objects, L d ≈
0. As collisions produce small particles, L d rises. In models with b l > L d but a smaller peak L d . Although the dust luminosity at late times depends on b d , m l, , and b l , all models follow L d ∝ t − n with n ≈ L d is muchlarger. As the size distribution establishes the standard wavy power law, L d varies slowlybefore settling in on some constant value which depends on b d , m l, , and b l . At late times, L ∝ t − n with n ≈ L d at t & L d is independent of the initial size distribution. Q ⋆D All analytic models for collisional cascades predict that the lifetime depends on theratio v /Q ⋆D (e.g., eq. 4; see also Wyatt et al. 2007a,b; Kobayashi & Tanaka 2010; Wyattet al. 2011). In our numerical approach, collisions between pairs of equal mass particles with v / Q ⋆D & v , Q ⋆D = 6 × − × erg g − ,and v / Q ⋆D = 0.2–200. For comparison, eq. 10 with standard parameters for rocky objectspredicts Q ⋆D ≈ erg g − for 0.1 km particles.Fig. 7 illustrates the evolution of r max for calculations with b d = 1, m l, = 0.2, b l = 1, andeither a mono-disperse (upper panel) or a power-law (lower panel) initial size distribution.When v / Q ⋆D .
1, the largest objects grow to sizes of 500–1000 km in 10–100 Myr. Forlarger ratios, v / Q ⋆D ≈ Q ⋆D and the initial sizedistribution. In this set of calculations, the loss time scales inversely with Q ⋆D for v / Q ⋆D ≈ v / Q ⋆D , the evolution time is fairly independent of Q ⋆D . For any Q ⋆D ,swarms with an initial power-law size distribution begin to lose mass more rapidly thanthose with an initial mono-disperse set of large objects. At late times, however, originallymono-disperse swarms lose more mass (see § Q ⋆D tend tohave somewhat smaller waves than those with smaller Q ⋆D . Otherwise, the level and shapeof the relative cumulative size distributions simply scale with the collision time.Fig. 8 demonstrates that the evolution of the dust luminosity also scales with Q ⋆D . Whenswarms begin with a power law size distribution, they have substantial L d . Destructivecollisions among mono-disperse swarms gradually build a comparable L d over a few collisiontimes. For this suite of calculations, the peak L d is fairly independent of the ratio v / Q ⋆D .At late times, all systems follow power law declines in L d with time. Setting L d ∝ t − n ,this suite of calculations has n = 1.1–1.2. There is some tendency for larger n at early timesand smaller n at late times, but our sample size is too small to verify this behavior.Throughout this evolution, systems with longer collision times have larger L d . For v / Q ⋆D & L d ∝ ( v / Q ⋆D ) − / . When v / Q ⋆D ≈ L d ∝ ( v / Q ⋆D ) − . Once v / Q ⋆D .
1, most of the solid mass ends up in large objects with negligible surface area.The dust luminosity is then close to zero.When Q ⋆D is independent of radius, all collisions either allow objects to grow (large Q ⋆D ) or produce debris (small Q ⋆D ). In real systems, Q ⋆D is a function of particle size with adistinct minimum at intermediate sizes r ≈ r max ina real system, we consider standard calculations using the Q ⋆D relation in eq. 10.Fig. 9 shows results of calculations for initially mono-disperse and power-law size distri-butions with b d = 1, m l, = 0.2, b l = 1, and r = 100, 300, 500, and 1000 km. In swarms witha starting r = 100 km and 300 km, the largest objects slowly lose mass with time. Systemswith smaller r max evolve more rapidly and lose ∼
97% ( r = 100 km) to 70% ( r = 300 km)of their original mass over 1 Gyr. As in § r is larger than 300 km, collisions among equal mass objects produce largermerged objects. The evolution is then very sensitive to the starting size distribution of theswarm. Among originally mono-disperse swarms, the largest objects reach sizes of 4500–6000 km fairly independently of the initial r max . In swarms with r = 500 km and a powerlaw size distribution, cratering collisions are almost frequent enough to prevent growth.After ∼
10 Myr, growth overcomes cratering, but the largest objects only reach sizes of1000–1500 km instead of 4000–5000 km. For swarms with a larger r , cratering is much lessimportant (Fig. 9, purple curves). 17 –These results follow from eq. 11. Among mono-disperse swarms, we expect growth when r exceeds 300 km; our numerical results confirm this estimate. When swarms start witha power-law size distribution, cratering collisions should prevent growth at some larger r .The numerical calculations suggest this limit is roughly 400 km.Compared to ‘normal’ collisional cascades, systems where particles grow have muchsmaller L d (Fig. 10). In modo-disperse swarms, the time scale for L d to rise from zero scaleswith r . When the largest objects cannot grow, collisions distribute debris among smallerobjects. Peak L d is then fairly independent of r (e.g., Fig. 10, thick orange and magentalines). In ensembles of growing large objects, collisions gradually concentrate more and moremass into larger and larger objects. Growth limits debris production. Peak L d is then muchsmaller (e.g., Fig. 10, thick green and violet lines). At late times, the largest objects sweepup all of the remaining small particles; L d then drops dramatically.Systems with initial power-law size distributions exhibit similar behavior (Fig. 10, thinlines). All swarms begin with large L d . As collisions destroy the largest objects ( r = 100 kmand 300 km), L d slowly declines with time. At late times ( t &
10 Myr), the evolution of L d is independent of the initial size distribution. When the largest particles grow with time ( r = 500 km and 1000 km), all collisions produce some debris; continuous debris productionmaintains a slow decline in L d . Eventually growth dominates debris production; L d thenplummets. For r = 1000 km, the substantial fall in L d begins at ∼
300 Myr (somewhatlater than the drop at ∼
100 Myr for a mono-disperse system). When r = 500 km, a closebalance between growth and debris production maintains large L d past 1–2 Gyr. After a fewmore Gyr, the largest objects sweep up debris and L d rapidly falls to zero.Despite the different evolution in r max and L d , calculations with Q ⋆D ( r ) still producewavy size distributions about a standard power-law (Fig. 11). For the fragmentation pa-rameters adopted here, the expected power law has slope q c ≈ § n ( > r ) /r − . . Independent of the initial size distribu-tion, b d , m l, , and b l , model results for r . r & § A.2; Wyatt et al. 2011). For ensembles with r = 100 km, cumulative size distributionsat 10–100 Myr roughly follow a power-law from 0.1 km to 30–50 km with n ( > r ) ∝ r − q c and q c ≈ q r ≈ − − q d In Figs. 2–6, the evolution of the cascade clearly depends on the algorithms for derivingthe amount of debris and for placing this debris in mass bins. To examine how the sizedistribution of the debris impacts the evolution, we consider calculations with a mono-disperse initial size distribution, b d = 1, m l, = 0.2, b l = 1, and q d = 3.0, 3.5, 4.0, 4.5, and5.0. When q d is small (large), most of the debris lands in bins with smaller (larger) averagemass. Based on results from the standard model, we expect the pace of the cascade todepend on q d .Fig. 12 shows the time evolution of r max as a function of q d . For an initially mono-disperse size distribution, catastrophic collisions between the largest objects place debristhroughout the mass grid. When q d ≈
3, debris is concentrated among the smallest particles;collisions remove more mass from the grid. When q d ≈
5, debris is concentrated among thelargest particles with little mass lost from the grid. Concentrating more debris among largeparticles enables a larger rate of mass loss from the largest objects. Thus, r max declines morerapidly when q d = 5 than when q d = 3.By the end of the calculation at 1 Gyr, catastrophic and cratering collisions substantiallyreduce the mass in the largest objects. When q d = 3 (5), r max = 8 km (6 km). Most ofthis reduction occurs during the first 1 Myr of the evolution. Once most of the mass in thegrid has been converted into 1 µ m or smaller particles which are lost to the grid, the timevariation of r max is slow.For each of these calculations, it takes roughly 0.1 Myr to reach the ‘standard’ wavypower-law size distribution. By 1 Myr, calculations with different q d establish various wavypatterns at small particle sizes (Fig. 13). These patterns remain fixed for the rest of theevolution. When q d = 3.5, there is a characteristic pattern of a large-amplitude wave at 1–100 µ m, several waves with diminishing amplitude at larger sizes, and a final wave of modestamplitude at r . r max (see also Campo Bagatin et al. 1994; Durda et al. 1998; O’Brien &Greenberg 2003; Kral et al. 2013). For other values of q d , the pattern for r & q d , thewavelengths seem independent of q d .Despite these curious differences in the size distribution at small sizes, the mass in thegrid is a simple function of q d . At a fixed time, cascades with larger q d have less massthroughout the grid than cascades with smaller q d . The origin of this result is clear. When q d is small (large), collisions place more mass in smaller (larger) mass bins. Smaller particlesremove less mass from the largest particles, which contain most of the mass in the grid. By 19 –the end of the calculations at 1 Gyr, cascades with q d = 3 have 35% to 40% more mass thancascades with q d = 5.These calculations also reveal interesting differences in the evolution of the relative dustluminosity L d as a function of q d (Fig. 14). At the onset of the cascade, infrequent andrandom collisions (e.g., shot noise) set the population of the small particles which containmost of the surface area. Thus, the timing of the abrupt rises in L d is random and fairlyindependent of q . During this phase of evolution, however, cascades with small q d lose moremass from the grid than cascades with large q d . Once cascades establish the ‘standard’ sizedistribution, those with smaller q d have less mass than those with larger q d . Thus, the peak L d is smaller in systems with smaller q d . This difference is substantial: cascades with q d =5 are 2.5 times brighter at peak L d than cascades with q d = 3.The timing of maximum L d also depends on q d . When q d is large, erosion of the largestparticles is more rapid (Fig. 7). More rapid erosion tends to produce a larger ensemble ofsmall particles. Thus, cascades with large q d reach maximum L d earlier than cascades withsmall q d . This conclusion is independent of the timing of the first rise in L d . In Fig. 8,the calculation with q d = 4 rises before those with q d = 4.5 or 5, but still reaches a smallermaximum L d .At late times, L d also correlates with q d . Systems with larger q d have larger L d at latetimes. Although somewhat counterintuitive, this outcome has a simple physical explanation.When q d is large, more erosion of the largest particles yields a larger maximum L d with asomewhat smaller total mass. Cascades with smaller masses and larger q d evolve more slowly(due to fewer collisions) and retain debris for somewhat longer (due to larger q d ). Thus, thesesystems tend to maintain their larger L d at later times.In all calculations, L d declines somewhat more rapidly than predicted by the analyticmodel. Our results follow L d ∝ t − n with n = 1.1–1.2 instead of the analytic n = 1. Theslow decline of r max with time is responsible for this difference. When r max is fixed as inthe analytic model, n must be unity as outlined in §
2. When r max declines with time, thecollision time becomes progressively shorter with time (eq. 1). As the collision time becomesshorter, the luminosity declines more rapidly. Our analysis suggest several clear differences between predictions of the standard ana-lytic model and the results of numerical simulations with Q ⋆D = constant (Figs. 2–8, 12–14).For material orbiting with e = 0.1 inside an annulus with Σ = 10 g cm − at 1 AU, (i) the 20 –radius (mass) of the largest object declines by a factor of 12–16 (2000–4000) in 1 Gyr, (ii)the size distribution is a power law which generally follows the N ( r ) ∝ r − . of the analyticmodel but has large amplitude waves at large and small sizes, and (iii) the evolution of theluminosity at late times is somewhat steeper than the predicted L d ∝ t − .When Q ⋆D is a function of r , the numerical calculations provide other tests of the analyticmodel. For e = 0.1 and i/e = 0.5 at a = 1 AU, swarms with r = 100–400 km follow theevolution of the constant Q ⋆D models fairly closely. Over 2 Gyr, (i) the largest objects losefrom 75% to 99% of their initial mass and (ii) the dust luminosity declines as L d ∝ t − n with n ≈ r . r Q = 0.1 km roughly tracks a powerlaw with the q ≈ r . µ m leads to an excess (deficit) of particles with r ≈ µ m (10–100 µ m). Amonglarger particles with r & r Q , the power law slope is smaller and has a broader range, q ≈ q ≈ r & r , there are large-amplitude waves about the average power law at largesizes.Calculations with r = 500–1000 km evolve differently. When e ≈ M d ∝ t − n with n . L d ∝ t − n with n ≈ § M ∝ t − (e.g., Kobayashi & Tanaka2010; Wyatt et al. 2011). This behavior allows the cascade to maintain an invariant sizedistribution. In our calculations, the shape of the size distribution is nearly fixed in time(Fig. 3–4), matching a basic prediction of the analytic model.When Q ⋆D = constant over a finite range of particle sizes, analytic models predict a setof waves superposed on a power-law size distribution (Wyatt et al. 2011). Our calculationsmatch the predicted slope and the amplitude of the primary wave close to the small size cutoffat 1 µ m (see § A.2). The numerical simulations produce (i) waves with longer wavelengthsand smaller amplitudes at 0.1 mm to 1 km and (ii) a larger wave at 1–100 km which are notpresent in the analytic model. 21 –When Q ⋆D is a function of radius, the numerical calculations match (i) the predictedslope of the power law from 1 µ m to 0.1 km, (ii) the amplitude of the wave at the smallestsizes, (iii) the distinctive rise in the relative cumulative size distribution at 0.1 km, and(iv) the general shape of the relative size distribution from 0.1 km to 10–30 km (see § A.2).However, the level of waviness at 0.1 mm to 0.1 km is larger than predicted; the slope of thesize distribution for r & ξ as the ratio of the mass in smallparticles to the total mass. In models with fixed Q ⋆D , r = 100 km, r min = 1 µ m, and q =3.5, ξ = 10 − . When Q ⋆D is a function of particle size, q = 3.68 (3.04) for r . r Q ( r & r Q );with r Q = 0.1 km, ξ ≈ − .Fig. 15 illustrates the behavior of ξ for a set of cascade models with r = 100 km. Whenthe swarm has an initial power law size distribution, ξ ≈ − . As systems with constant Q ⋆D evolve (Fig. 15, thin black and green curves), the relative dust mass gradually growswith time and reaches ξ ≈ − × − at t = 1–2 Gyr. When Q ⋆D is a function of particlesize (thin orange and magenta lines), ξ declines to 2 − × − at 0.1–1 Myr and then slowlyrecovers to 4 − × − at t = 1–2 Gyr. Swarms with mono-disperse initial size distributionsfollow similar paths. Once collisions begin to produce dust, the relative dust mass growsto ξ ≈ − (for constant Q ⋆D swarms) or ξ ≈ − × − (for swarms with Q ⋆D ( r )). Therelative dust mass then grows slowly with time.For swarms with constant Q ⋆D , gradual reductions in r max explain the evolution of therelative dust mass. After 1–2 Gyr of evolution, r max ≈ q = 3.5,the expected relative dust mass is ξ ≈ − × − . Allowing for small differences in dustmass and total mass due to the wavy size distribution (Fig. 5), the ξ ≈ − × − fromour calculations is close to the expected value.When Q ⋆D is a function of r , the change in the size distribution at large sizes is asimportant as the reduction of r max with time. For swarms with a double power-law sizedistribution and a break at r Q ≈ ξ ≈ (cid:18) − q l − q s (cid:19) (cid:18) r − q s d r − q l max (cid:19) r q s − q l Q , (16)where r d = 1 mm is the radius of the largest dust particle and q s , q l are the slopes of the size 22 –distribution for large and small particles. In our calculations, the slope of the size distributionfor small particles – q s ≈ .
68 – is close to expectation. At large sizes, however, the slope isshallower than predicted, q l ≈ q l ≈ ξ ≈ − × − for r max = 100 km and ξ ≈ − × − for r max = 50 km. Our calculations match these expectations.For the time evolution of ˙ M and L d in the analytic model, we derive an estimate forthe collision time from eq. 5 with α = 1. Adjusting α allows us to match the analytic˙ M = M /t c at t = 0 with the initial production rate for 1 µ m and smaller objects from ournumerical simulations. This mass loss rate then declines with time as˙ M = t − c M (1 + t/t c ) . (17)Adopting an initial power law size distribution with q = 3.5 for Q ⋆D = constant and r ≤ r max yields the initial L d which then evolves as in eq. 6.Fig. 16 compares ˙ M derived from our baseline calculations with predictions of the an-alytic model for α = 0.01, 0.1, and 1.0. At t = 0, the numerical calculations match theanalytic model for α ≈ M declines slightlymore rapidly than the analytic prediction. By the end of the calculation, the mass loss rateis similar to analytic model predictions for α ≈ Q ⋆D (e.g., Fig. 7). At late times, ˙ M ≈ M t c t − . From eq. 5, t c ≈ t α with t ≈ × yr.Substituting this result into the expression for ˙ M yields α ≈ ˙ M t /M t . Numerical resultsfor ˙ M ( t ) thus yield α ( t ). For models with Q ⋆D = 6 − × erg g − , we infer: α ≈ (cid:18) v Q ⋆D (cid:19) − (cid:18) t (cid:19) − . . (18)At 0.1–1 Gyr, the numerical calculations evolve roughly 3 times more slowly than the analyticmodel (e.g., eq. 4) when v / Q ⋆D ≈ v / Q ⋆D ≈ M .Fig. 17 compares the baseline L d ’s for the simulations with the analytic model. At earlytimes, the number of small particles in the calculations oscillates about a constant value as 23 –the system establishes an equilibrium size distribution. This constant is a factor of 4–5 largerthan the prediction of the analytic model. At ∼ L d in the numerical model declines fasterthan predicted by the analytic model, L d from the numerical models eventually falls belowthe analytic prediction.In the numerical models, the excess of small particles at r ≈ r min produces the larger L d . Although these particles remove larger particles from the size distribution, the excesssurface area from the small particles more than compensates for the deficit in surface area atsomewhat larger sizes. Thus, the numerical models have larger L d than the analytic model.Despite the spikes in the mass loss rate, the numerical models show no spikes in L d .When r max ≈ M sprinkles debris throughout the grid, with noimpact on the shape of the size distribution and little impact on the total mass in the grid.Thus, there is negligible change in surface area and L d . In § r max produce clear spikes in L d .Results for models where Q ⋆D depends on particle size yield similar results. The dustluminosity and the mass loss rate from the grid generally follow the trends established incalculation with constant Q ⋆D : (i) a plateau in ˙ M or L d followed by a decline and (ii)occasional spikes in ˙ M (but not L d ) as the system evolves. Our analysis pinpoints several clear differences between the analytic model and detailedcoagulation calculations of collisional cascades. However, real cascades do not occur inisolation. Within planet formation theory, cascades begin as the gravity of growing or fully-formed planets stirs up surrounding smaller particles (e.g., Kenyon & Bromley 2002b; Mustill& Wyatt 2009). As these systems evolve, the collisional cascade converts a fraction f c ofthese small particles into dust grains which are ejected by radiation pressure from the centralstar. Massive planets accrete or dynamically eject the rest of the small particles. Deriving f c and making a more robust link between the theories of collisional cascades and planetformation requires more detailed sets of calculations.To begin to make this link, we consider a suite of planet formation calculations. Withina single annulus, the initial swarm consists of a set of mono-disperse particles with totalmass M d = 0.5 M ⊕ , surface density Σ = 10 g cm − , and initial radius r = 10 n ( n = 0, 1, 2, . . . , 8). Particles have e = 10 − ( r = 1–10 cm), 10 − ( r = 1–100 m), 10 − ( r = 1–10 km),10 − ( r = 1–100 km), and i = e /2. As particles grow, collisional damping, dynamical 24 –friction, and viscous stirring modify e and i for each mass bin. All calculations ignore gasdrag and Poynting-Robertson drag.Although the initial conditions strongly favor growth through mergers, every collisionproduces a modest amount of debris. To specify the amount and size distribution of thedebris, we adopt the standard Q ⋆D ( r ) relation with parameters specified after eq. 10 and set b d = 1, m l, = 0.2, and b l = 1 in eqs. 14–15. Modest changes in these assumptions have littleimpact on the results.To quantify the diversity of outcomes, we perform ∼ Orchestra uses a random number generator to decide whetherto round up (or down) a fractional number of collisions to an integer. Randomness in thecollision rates generates a dispersion in outcomes.Despite some intrinsic diversity, all calculations follow a standard pattern. Initially,particles grow through mergers. As particles grow, collisional damping and gravitationalinteractions modify e and i for each mass bin. Damping reduces e and i for small particles, r . r & e and i forall particles.As e and i evolve, the largest particles step through several stages of growth. At thestart of each calculation, gravitational focusing factors f g are small. Growth is slow andsteady. Damping and dynamical friction enhance f g , enabling a phase of runaway growthwhere the largest particles gain mass much more rapidly than smaller particles. Once viscousstirring dominates collisional damping and dynamical friction, gravitational focusing factorsfor the largest particles drop. Oligarchic growth – where large particles add mass more slowly– begins.Rising viscous stirring rates also initiate the collisional cascade. Collisions among parti-cles close to the minimum in Q ⋆D – at roughly r Q ≈ e and i for small particles, collisionsdestroy particles which are smaller and larger than r Q . With Q ⋆D ∝ r − . for r < r Q and Q ⋆D ∝ r . for r > r Q , it is easier for collisions with fixed kinetic energy to destroy particleswith r ≪ r Q than those with r ≫ r Q . These collisions produce debris with sizes r .
100 cm.As the debris accumulates, there are two possible outcomes (Kenyon & Bromley 2015,and references therein). If collisional damping dominates the loss of particles from destructivecollisions, small particles have very small e and i . Damping then limits the collisional cascade.Large gravitational focusing factors initiate a second phase of runaway growth, where the 25 –largest objects accrete most of the mass in small particles. If damping is ineffective, smallparticles remain at large e and i ; the collisional cascade proceeds unabated.Fig. 18 illustrates the growth of the largest particles as a function of r . From eq. 1,small particles evolve more rapidly than large particles. The ensemble of 1 cm particlesrequires less than a year to reach sizes of 30 cm; after another 20 yr, the largest particleshave radii of 1 km. By ∼ yr, the radii of the largest protoplanets surpass 10 km. As theevolution slows, occasional collisions among these objects eventually yields a single planetwith a radius of roughly 6000 km. In contrast, sets of 100 km or 1000 km particles take morethan 0.1 Myr for the largest objects to double in mass. Despite the delay, protoplanets stillachieve radii of roughly 6000 km on a time scale of 1 Myr.Throughout this evolution, changes in the relative cumulative size distribution follow afairly standard pattern (Fig. 19). For a swarm of particles with r = 10 cm, mergers rapidlyproduce a swarm of particles with r ≈ ∼
100 yr, the size distribution has three mainpieces: (i) an exponential at large sizes, (ii) a plateau at intermediate sizes, and (iii) a debristail at small sizes. Over the next 10 yr, the exponential and the small peak at the smallsize end of the plateau shift to larger and larger sizes. Large particles sweep up the smallestparticles, gradually diminishing the population of the debris tail.From ∼ yr to ∼ yr, the system makes a transition from runaway growth tooligarchic growth. During oligarchic growth, the size distribution develops a small rise at r ≈ r ≈ r Q ≈ r & ∼
100 Myr, the size distribution consists of a steeplyrising piece for r &
10 km and a fairly flat piece for r .
10 km. By 1 Gyr, accretion by 1–2large objects and destruction by the collisional cascade remove all small particles from theannulus.The sharp minimum at r ≈ r damp . Defining e H = ( m max / M ⊙ ) / as the Hill eccentricity ofthe largest object with mass m max , we express particle eccentricities as e rel = e/e H . (19)Particles with r & r damp have e rel ≈ r . r damp have e rel ≈ e rel between small and large particles tempers the collisionalcascade. Although collisions between particles with r & r damp produces substantial debris,collisions among smaller particles often yield larger merged objects. As material cyclesbetween small and large objects, the largest objects accrete solids from both reservoirs.Eventually, ∼
90% of the initial mass ends up in the largest objects. The collisional cascadeconverts only ∼
10% of the initial mass into 1 µ m and smaller particles.Despite significantly different evolution times for systems with different r (Fig. 18), allsystems have similar size distributions at 0.1–1 Gyr (Fig. 20). Every swarm has one or twoobjects which contain nearly all of the mass. For r &
10 km, the cumulative size distributionapproximately follows a power law with N ( > r ) ∝ r − q c and q c ≈ r ≈ r ≈ µ m to 1 cm, the size distributions are steeper with morepronounced waviness. This feature depends on the timing of the last large collision betweenobjects with r &
100 km. Systems with more recent large collisions have larger rises to smallsizes and more waviness along this rise.To examine the observable consequences of this evolution, we consider the dust luminos-ity L d (Fig. 21). At early times, growth concentrates mass into objects with smaller surfacearea; L d declines with time. Once the collisional cascade produces sufficient debris, L d rises.For systems with r . ∼ yr. Among swarms withlarger r , the rise in L d occurs at progressively later times, reaching ∼ r =100 km and 5–10 Myr when r = 1000 km.For our calculations, the peak relative dust luminosity, L d /L ⋆ ≈ − × − is relativelyindependent of r . Close to the peak, stochastic variations in the collision rate produce 50%variations in L d . The magnitude of these variations is independent of r . Repeat calculationswith identical initial conditions (but different random number seeds) yield similar variationsand similar peak L d .When r & r & L d for all swarms gradually declines. As the systems decline,fluctuations – including large spikes in L d – become more pronounced. The magnitude andfrequency of the spikes is independent of r ; however, the spikes are clearly more obvious atlater times. Although the overall evolution of L d results from countless collisions of 1–100 kmobjects, binary collisions between 300–1000 km objects produce the large spikes. These giantimpact events release dust masses comparable to (and sometimes exceeding) the total dustmass in the rest of the swarm.With few large objects in the swarm, impact events are rare and random (see also Gendaet al. 2015). The typical spacing between large collisions is ∼ ∼ & . . Kepler satellite might identify 1–2 transitevents from impact debris during the lifetime of the satellite.To develop a better understanding of the spikes in L d , we define the break radius r brk asthe maximum size for a particle destroyed in a collision between equal mass objects. In ourapproach, Q c = v / Q ⋆D is the collision energy required to eject half the mass of the colliding pair to infinity.As the largest objects grow, their gravity stirs up smaller objects to larger and larger e and i . Particles with larger e collide with larger velocities. For fixed r , the ratio Q c /Q ⋆D growswith time; collisions among pairs of objects with identical r eject more mass at later times.Fixing Q c /Q ⋆D = 1, larger collision velocities result in collisions that destroy pairs of objectswith larger and larger r brk as the system evolves.Fig. 22 shows the time evolution of r brk for swarms with r = 1 cm to 1000 km. After 28 –10 yr, all systems follow the same pattern: r brk gradually grows from roughly 1 km at0.01 Myr to roughly 1000 km at 10 Myr. After 10 Myr, one or two objects contain nearly allof the mass. Stirring of e and i effectively ceases; r brk remains roughly constant with time.This behavior has a profound impact on L d (Fig. 21). When r brk is small, every collisionamong roughly equal mass objects yields a merger (for large objects) or a small amountof debris (for small objects). This debris has little impact on L d . As r brk grows, binarycollisions among objects with r ≈ r brk produce more and more debris. At peak L d , thesecollisions cause modest fluctuations in L d . As L d declines, collisions among larger and largerobjects disperse more and more debris, resulting in more prominent spikes in L d . By theend of the calculation, giant impacts among 1000 km objects lead to factor of 5–10 changesin L d .Repeat calculations confirm this behavior. Modest changes in e and i have littleimpact on the results. Over 2 Gyr of evolution, the collisional cascade converts 10% ±
3% ofthe initial mass into 1 µ m and smaller particles. The rest of the initial mass lies in 1–2 largeobjects with radii of 5000–6000 km. During the early stages of the cascade (0.1–10 Myr),the peak dust luminosity is L d /L ⋆ ≈ − × − . Near peak brightness, factor of 1.5–2 fluctuations in L d are typical. For older systems (10–100 Myr), a slow decline in L d isoccasionally interrupted by factor of 2–10 spikes in the dust luminosity. As collisions destroylarger and larger objects, the amplitudes of these spikes gradually grow with time. After ∼
100 Myr, the dust luminosity drops by several orders of magnitude. The timing of thisdecline has a broad range, 100–300 Myr, and is fairly independent of r and e . Our suite of numerical calculations yields several interesting insights into the evolutionof collisional cascades.For ‘pure’ cascade models with constant e and i , we confirm several aspects of theanalytic model. Throughout the evolution of systems where Q ⋆D in independent of r , thesize distribution generally follows the expected N ( r ) ∝ r − . . When Q ⋆D is a function of r ,the size distribution follows the predicted N ( r ) ∝ r − . for r . r Q ≈ r & r Q ,the typical slope of q ≈ q ≈ § A.2, see also Campo Bagatin et al. 1994; O’Brien & Greenberg 2003; Wyatt et al. 2011;Kral et al. 2013).The long-term evolution of cascades is close to the predictions of the analytic model. In 29 –our calculations, the decline is somewhat faster – L d , M d ∝ t − n and n = 1.1–1.2 – than theanalytic solution, where L d , M d ∝ t − . The ratio ξ of the dust mass to the total mass slowlygrows with time, ξ ∝ t n and n = 0.05–0.15. For fixed e and i , the collision time scale for thenumerical simulations – t c ∝ ( v / Q ⋆D ) − – is reasonably close to the analytic result where t c ∝ ( v/ Q ⋆D ) − / .A major difference between analytic and numerical models for collisional cascades is theevolution of r max . Analytic models assume constant r max . In the numerical calculations, r max either grows ( v / Q ⋆D .
1) or declines ( v / Q ⋆D &
1) with time. After 1 Gyr of evolution,factor of 2–10 changes in r max are typical. In systems where r max declines with time, ourresults yield r max ∝ t − n with n = 0.1–0.2. With t c ∝ r max , we then expect L d ∝ t − n with n = 1.1–1.2. Thus, the gradual decline of r max with time accounts for the more rapid declineof L d with time. A changing r max is likely responsible for other, more subtle differencesbetween the analytic and numerical model.Planet formation simulations yield important insights into the onset and evolution ofreal collisional cascades. At 1 AU, systems of small particles with Σ ≈
10 g cm − evolveinto a few large oligarchs in 0.1–1 Myr. The time scale to produce 2000–4000 km oligarchsis fairly independent of r . For t & r .For collisional cascades, the unique feature of planet formation simulations is the steadyevolution of r brk , the maximum particle size destroyed in collisions of equal mass objects.In analytic models, r brk = r max = constant. In numerical simulations of planet-formingdisks, growing oligarchs gradually stir up the orbits of smaller solid particles; r brk growspersistently with time. As r brk evolves, more and more material becomes involved in thecollisional cascade. Thus, L d and M d decline more slowly with time.The relative dust mass ξ also evolves differently in planet formation calculations. Incascades with r ≈ ξ gradually grows over time (Fig. 15). In cascades with r = 500–1000 km and in planet formation calculations, the relative dust mass is usually muchsmaller ( ξ ≈ − − − ) and fluctuates by a factor of 2–20 when giant impacts add debristo the swarm. As the largest objects sweep up debris, planet formation calculations exhibitlarge drops in ξ which do not occur in cascade calculations.To focus on observable differences between numerical calculations of planet formationand pure cascades, Fig. 23 compares the evolution of L d /L ⋆ for two collisional cascadecalculations and two planet formation calculations. For the cascade calculations, we set r max = 100 km and e = 0.1 (black curve) and r max = 1000 km and e = 0.25 (violet curve).For these cascades, we set the initial mass M = 0.05 M ⊕ , the mass typically lost from the 30 –planet formation calculations. Tracks for the two planet formation calculations repeat thosefrom Fig. 21 for r = 10 cm (magenta curve) and r = 1 km (cyan curve).The four sequences show several similarities in the time evolution of L d . All have roughlythe same peak luminosity, L d /L ⋆ ≈ − × − . The rise time for the 100 km collisionalcascade model is similar to rise times for the two planet formation calculations. The largespikes in L d for the 1000 km collisional cascade model are similar to those in the two planetformation calculations.Several differences in the evolution are also apparent. When r brk ≈
100 km in the twoplanet formation calculations (0.1–1 Myr), L d shows significant fluctuations not visible inthe collisional cascade calculation with r = 100 km. When r brk evolves well past 100 kmin the planet formation calculations, the decline in L d is much slower than in the collisionalcascade calculation with r = 100 km. Finally, when the largest object in the planet formationcalculations has accreted most of the mass in the annulus, L d plummets. In contrast, bothcollisional cascade calculations maintain a t − . decline for much longer time scales.
5. DISCUSSION5.1. Theoretical Issues
To follow the evolution of a swarm of planetesimals into a planetary system, coagulationcalculations include a broad set of physical processes (e.g., Greenberg et al. 1978; Wetherill& Stewart 1993; Weidenschilling et al. 1997; Kenyon & Luu 1999; Kenyon & Bromley 2008;Kobayashi et al. 2010; Ormel et al. 2010; Raymond et al. 2011; Glaschke et al. 2014; Shannonet al. 2015; Johansen et al. 2015, and references therein). Modeling each physical processinvolves several assumptions and choices for various input parameters. Here, we considerhow these decisions impact our results and conclusions.For this study, we ignore interactions between solid particles and a gaseous disk. Al-though gaseous circumstellar disks are a major component in planet formation models (e.g.,Youdin & Kenyon 2013), all analytic models of collisional cascades assume the gas has al-ready dissipated. Neglecting the gas allows us to pinpoint common features in numericalcalculations of collisional cascades and planet-forming disks.Although we consider radiation pressure from the central star to set r min , we do notinclude Poynting-Robertson drag. For 1 µ m particles at 1 AU, the collision time is shorterthan the time scale for Poynting-Robertson drag when L d /L ⋆ & − (e.g., Backman &Paresce 1993; Wyatt 2008). In most of the calculations reported here, the dust luminosity 31 –meets this limit.Unlike several numerical algorithms (e.g., Krivov et al. 2006; G´asp´ar et al. 2012a; Kralet al. 2013), our calculations also neglect the radiation force on fragments ejected from acollision. Defining β as the ratio of the radiation force to the gravitational force, debrisparticles have orbits with semimajor axis a ′ = ( − β − β ) a and eccentricity e ′ = β/ (1 − β ) (e.g.,Burns et al. 1979). Unbound particles with β & r . µ m leave the swarm on a timescale comparable to the orbital period. For swarms with Σ = 10 g cm − at 1 AU, collisionsbetween bound and unbound particles are rare. Particles with β ≈ r ≈ µ m)remain bound to the central star but lie outside the single annulus in our calculations.Despite its absence from current analytic models, it is worth estimating the impact ofhigh- β particles on our planet formation simulations. With larger a and e , these particleshave longer lifetimes than particles inside our single annulus (Wyatt et al. 2010). Fewercollisions tend to increase the population of particles with (i) r & µ m inside the annulusand (ii) r . µ m outside the annulus. Despite producing a steeper size distribution at smallsizes, the impact on L d is probably small (Wyatt et al. 2010). With negligible total mass,the abundance and orbits of this population have little impact on the growth of the largestobjects. We plan multiannulus calculations to investigate these possibilities in more detail.To calculate collisional processes, we specify a variety of parameters for collisional damp-ing, the cross-section, gravitational stirring, and fragmentation. Comprehensive analytic and n -body calculations verify the algorithms for collisional damping, the collision cross-section,and stirring (e.g., Greenzweig & Lissauer 1990; Spaute et al. 1991; Ohtsuki 1999; Lee 2000;Ohtsuki et al. 2002; Goldreich et al. 2004; Levison & Morbidelli 2007; Kobayashi & Tanaka2010; Kenyon & Bromley 2015). Within the fragmentation algorithm, we derive materiallost from ‘cratering’ and ‘catastrophic disruption’ using the expressions in eq. 10, eqs. 14–15,and a power law size distribution for the debris. For cratering collisions where the mass ofone component (the ‘projectile’) is much smaller than the other (the ‘target’), results for theejected mass agree with detailed analyses of laboratory measurements (see the discussion in § n -body simulations (e.g., Benz & Asphaug 1999; Leinhardt& Stewart 2009, 2012). At present, these algorithms are state-of-the-art.For cascade models, Fig. 8 illustrates the impact of changing Q ⋆D on the collision timeand the dust luminosity. For fixed e and i , systems where Q ⋆D is small evolve more rapidlythan systems where Q ⋆D is large. Within a suite of planet formation calculations, the impactof Q ⋆D is less dramatic (e.g., Kenyon & Bromley 2008, 2010, 2012). Although small particleswith r . Q b and β b probably have little impact 32 –on the evolution (e.g., Kenyon & Bromley 2008). Among larger objects with r & r Q , self-gravity limits the ability to change Q ⋆D (e.g., Benz & Asphaug 1999; Leinhardt & Stewart2012). However, even large changes in Q g and β g may not change L d significantly (e.g.,Kenyon & Bromley 2012).Once a few large objects contain most of the mass in the annulus, the approximationsused in our statistical estimates for collision rates and stirring begin to break down. As thelargest protoplanets continue to accrete material from the rest of the swarm, their dynamicalinteractions become stronger and more chaotic. Among these large objects, our calculationsthen tend to overestimate the growth rate and underestimate e and i .To avoid these issues in several previous investigations, we promote the largest objectsinto the n -body component of Orchestra , where we follow the orbits of each large protoplanet(e.g., Kenyon & Bromley 2006; Bromley & Kenyon 2011a; Kenyon & Bromley 2014). In thisstudy, our main focus is on the evolution of the collisional cascade instead of the final orbitsof the largest objects. To save cpu time for a broad suite of simulations, we therefore chosenot to promote the largest protoplanets. When these protoplanets contain nearly all of themass, calculations for the rest of the swarm – including the flow of mass from r brk to r min and the accretion of small objects by the largest objects – remain accurate. The evolutionof the largest objects then have little impact on the evolution of the mass and luminosity ofthe small particles.In multiannulus calculations of terrestrial planet formation, promotion into the n -bodycode is more important (e.g., Kenyon & Bromley 2006; Raymond et al. 2011, 2012). Whenprotoplanets contain more than half the mass, they undergo a period of chaotic growthwhere the largest protoplanets scatter smaller objects throughout the terrestrial zone. Gi-ant impacts are then more frequent and occur at higher velocity than in pure coagulationcalculations (see also Asphaug et al. 2006; Genda et al. 2015). Thus, our single annuluscalculations probably somewhat underestimate the number and magnitude of large spikes indust luminosity on time scales &
10 Myr. Our main conclusion – that the dust luminosityin coagulation calculations declines more slowly than in collisional cascade calculations – isunchanged.
In the past 10–15 yr, there have been many numerical studies of collisional cascadesand planet-forming disks using coagulation and related techniques. Most investigationsconcentrate on the evolution at large distances, 10–200 AU, from the host star (e.g., Kenyon 33 –& Bromley 2002a; Krivov et al. 2006; Th´ebault & Augereau 2007; L¨ohne et al. 2008; G´asp´aret al. 2013; Kral et al. 2013; Kobayashi & L¨ohne 2014). Where possible, we compare ourresults with these studies and the few analyses which focus on the long-term evolution offairly massive disks in the terrestrial zone.In systems with large v / Q ⋆D and short t c , the linear decline of mass with time is abasic feature of the analytic model for collisional cascades (e.g., Wyatt & Dent 2002; Dominik& Decin 2003; Wyatt 2008; Kobayashi & Tanaka 2010). Aside from Kobayashi & Tanaka(2010), results from other studies do not test this prediction: the largest objects (i) areallowed to grow or (ii) have very long collision times (e.g., Krivov et al. 2006; L¨ohne et al.2008; G´asp´ar et al. 2012b, 2013; Kral et al. 2013).For all of our calculations, the disk mass declines with time somewhat more rapidlythan the analytic prediction: M d ∝ t − n with n = 1.1–1.2 instead of 1. In every calculation,catastrophic and cratering collisions slowly reduce the size of the largest objects. As theradius of the largest objects grows smaller, the collision time also decreases (eq. 1). Shortercollision times speed up the evolution, resulting in a larger mass loss rate per unit time.Thus, the slow reduction in r max produces a faster evolution of M d . Calculations with othercodes (e.g., L¨ohne et al. 2008; Kobayashi & Tanaka 2010; Kral et al. 2013) could test thisresult.Tests of the variation of t c with Q ⋆D are currently ambiguous. The numerical simulationsof Kobayashi & Tanaka (2010) confirm the analytic result, t c ∝ ( v /Q ⋆D ) − / . However,L¨ohne et al. (2008) quote t c ∝ ( v /Q ⋆D ) − / . Our result, t c ∝ ( v /Q ⋆D ) − , lies roughlymidway between these two extremes.All calculations produce wavy size distributions (e.g., Krivov et al. 2006; L¨ohne et al.2008; G´asp´ar et al. 2012b; Kral et al. 2013). For small particles and intermediate massparticles with r . r max and theslope of the size distribution as a function of time.When numerical simulations consider the long-term evolution of objects with r & v / Q ⋆D &
1, giant impacts which release copious amounts of dust are inevitable(e.g., Agnor et al. 1999; Agnor & Asphaug 2004; Genda et al. 2015). The giant impact eventsin our calculations are similar to those in Weidenschilling (2010a) and Raymond et al. (2011).To conclude this section, our results for planet formation time scales in the inner solarsystem are similar to those reported previously (e.g., Kominami & Ida 2002; Kenyon & 34 –Bromley 2004b; Nagasawa et al. 2005; Raymond et al. 2005; Kenyon & Bromley 2006; O’Brienet al. 2006; Kokubo et al. 2006; Chambers 2008; Raymond et al. 2014, and references therein).Typically, it takes 10 − yr to produce lunar mass protoplanets. Collisions combine theseobjects into a few 3000–6000 km objects over 10–100 Myr. Because we neglect gas drag, ourcalculations take a factor of ∼ Based on a set of multiannulus coagulation calculations at 0.68–1.32 AU, Kenyon &Bromley (2004b, hereafter KB04) first predicted the detectability of warm dust (T ≈ . In these studies,the predicted mid-IR (10–25 µ m) emission from dust reaches a peak relative flux F d /F ⋆ ≈ F d /F ⋆ ≈ ∼
10 Myr. Superimposed on thisdecline, spikes from giant impacts raise the relative flux by factors of 3–10. Rapid variabilityof dust production on 1–100 yr time scales is another characteristic of these calculations.Recent investigations confirm these results (e.g., Weidenschilling 2010a; Raymond et al.2011, 2012; Jackson & Wyatt 2012; Genda et al. 2015). Once 500–1000 km objects form, de-structive collisions produce copious amounts of debris throughout the terrestrial zone. Giantimpacts, including those capable of producing the Earth-Moon system, generate additionaldebris. As the overall level of debris declines, sporadic spikes in emission from the giantimpacts become more prominent. At late times, these spikes dominate dust production.Although the new calculations in § r .
100 km is somewhatdifferent. By explicitly following the dynamical evolution of 1 µ m to 1 m particles, we nowidentify an evolutionary phase at 1–10 Myr where collisional damping among mm- and cm-sized particles overcomes viscous stirring by large oligarchs. During this epoch, the lackof destructive collisions among small particles effectively halts the cascade. With a slowercascade, these systems lose a smaller fraction of their initial mass and are 2–3 times fainterat peak brightness than the debris disks in KB04 and KB05. However, the disks also remainbrighter for a factor of 2–3 longer, allowing the debris from terrestrial planet formation to KB05 also predicted levels for warm dust emission in the terrestrial zones of A-type stars.
35 –remain detectable among older stars.Scaling the results of the single 0.9–1.1 AU annulus in § F d /F ⋆ ≈ µ m and F d /F ⋆ ≈ µ m. Occasional bright spikes from giant impactscan raise the peak flux by a factor of 3–10. Although the bulk emission usually declines belowdetectable levels at 10–20 Myr, giant impacts occasionally raise the system above currentdetection limits until ∼
100 Myr (see also Raymond et al. 2011, 2012; Jackson & Wyatt2012; Genda et al. 2015).Calculations with r & Spitzer or WISE .Despite the identification of many cold debris disks around solar-type stars, initial at-tempts to detect warm emission with
IRAS and ground-based data were discouraging (e.g.,Weinberger et al. 2004; Mamajek et al. 2004). Song et al. (2005) then discovered a veryluminous debris disk associated with the old field star BD+20 ◦ IRAS also measured asignificant 12 µ m excess associated with HD 23514, a member of the Pleiades cluster; confir-mation as a warm debris disk required data from ISO , Spitzer , and ground-based telescopes(Spangler et al. 2001; Rhee et al. 2008).Over the past decade, data from
AKARI , Spitzer , and
WISE have enabled many surveysfor 12–25 µ m emission from debris orbiting solar-type main sequence stars (e.g., Stauffer et al.2005; Silverstone et al. 2006; Beichman et al. 2006; Siegler et al. 2007; Gorlova et al. 2007;Currie et al. 2007; Meyer et al. 2008; Trilling et al. 2008; Carpenter et al. 2009a,b; Mo´oret al. 2009; Stauffer et al. 2010; Beichman et al. 2011; Chen et al. 2011; Smith et al. 2011;Zuckerman et al. 2011; Luhman & Mamajek 2012; Ribas et al. 2012; Urban et al. 2012;Kennedy & Wyatt 2012; Zuckerman et al. 2012; Fujiwara et al. 2013; Kennedy & Wyatt2013; Cloutier et al. 2014). Current compilations include stars with ages ranging from ∼
10 Myr to & F d /F ⋆ ≈ µ m. In some systems,high quality spectra and imaging data confirm the presence of silicate dust orbiting within1–3 AU of the central star (e.g., Smith et al. 2008; Currie et al. 2011; Olofsson et al. 2012;Smith et al. 2012; Schneider et al. 2013; Ballering et al. 2014). For several others, short timescale variations in F d /F ⋆ suggest rapid changes in the typical sizes and total mass of smallparticles (e.g., Melis et al. 2012; Meng et al. 2012, 2014, 2015). Interferometric observationsfurther reveal 300 K or hotter dust in some systems (e.g., Absil et al. 2013; Ertel et al. 2014;Mennesson et al. 2014; Defr`ere et al. 2015). 36 –Taken together, the complete set of data suggest a fairly low warm dust frequency forsolar-type stars (see also Stauffer et al. 2005; Gorlova et al. 2007; Carpenter et al. 2009a,b;Chen et al. 2011; Luhman & Mamajek 2012; Ribas et al. 2012; Fujiwara et al. 2013; Kennedy& Wyatt 2013; Patel et al. 2014; Cloutier et al. 2014). Among older stars with ages &
100 Myr, only ∼ F d /F ⋆ & µ m. For youngerstars, the frequency is ∼
1% for 100 Myr-old stars and ∼
2% to 3% for 10–20 Myr stars.These statistics are generally consistent with the standard theoretical picture where acollisional cascade and occasional giant impacts generate a large population of small particlesin the terrestrial zone (see also Raymond et al. 2011; Fujiwara et al. 2013; Kennedy & Wyatt2013; Genda et al. 2015). As in Fig. 23, emission from the cascade gradually declines withtime from 1–10 Myr to ∼ Spitzer and
WISE detection limits ( F d /F ⋆ ≈ µ m). Based on current planet detection statistics, we assume 50% of solar-type starshave at least one Earth-mass planet at 0.4–2 AU (e.g., Youdin 2011; Petigura et al. 2013;Foreman-Mackey et al. 2014; Burke et al. 2015). For these systems, the initial surfacedensity distribution adopted in our calculations is a reasonable choice (e.g., Chambers 2001;Kominami & Ida 2004; Nagasawa et al. 2005; Kenyon & Bromley 2006; Kokubo et al. 2006;O’Brien et al. 2006; Chambers 2008; Fogg & Nelson 2009; Lunine et al. 2011; Raymond et al.2011, 2012; Hansen & Murray 2012, 2013).Among this set of calculations, the expected detection frequency f d depends on stellarage and r , the maximum size of solid objects at the start of the calculation. For any r , f d . f d . Spitzer and
WISE .Observations of dust emission for younger stars appear to rule out the possibility thata large fraction of planetary systems produce Earth-mass planets starting from 1000 km orlarger objects. When r & f d . r .
500 km 37 –tend to produce too much dust. Among the youngest stars (1–5 Myr), the predicted f d isa strong function of r , ranging from f d ≈ r .
100 km to f d .
5% for r ≈
500 km. For older stars with ages of 5–10 Myr (10–30 Myr), f d ≈ r . Given current constraints on warm dust for stars with ages of 1–30 Myr,our predictions agree with observations for 20–30 Myr old stars but are too large by factorsof 2–5 for younger stars. The differences between theory and observations for these youngstars are fairly independent of r .Among younger stars, the difference between theory and observations is probably lessextreme than indicated by these statistics. In our simulations, the typical level of dustemission at 1–10 Myr is only a factor of ∼ r & µ m (Plavchan et al. 2005; Wurz 2012),(iv) small particles are much weaker than assumed, (v) collisional damping is more effectiveamong small particles, or (vi) protoplanets accrete small particles more efficiently. Althoughinvestigating these possibilities is beyond the scope of this paper, we plan to analyze theirlikelihood in future studies. Better limits on the frequency of debris disks with relative fluxes3–5 times smaller than current detection limits would help to guide this analysis.
6. SUMMARY
We analyze a suite of coagulation calculations to examine the long-term evolution ofcollisional cascades in a ring of solid material with surface density Σ ≈
10 g cm − at 1 AUfrom a solar-mass star. In pure cascade models where the orbital e and i of the solids remainconstant with time, destructive collisions slowly reduce the size of the largest objects byfactors of 3–20. Over 1–2 Gyr, the total mass, the mass in small particles ( r . t − n , with n = 1.1–1.2. The ratio ξ of the mass in small particlesto the total mass grows by a factor of 2–3 from 1 Myr to 2 Gyr. After 1 Myr, the shapeof the size distribution for r & µ m is fairly constant in time and approximately follows apower law, N ( r ) ∝ r − q . The slope q depends on the form of the relation for the bindingenergy Q ⋆D : for constant Q ⋆D , q = 3.5; for Q ⋆D ( r ), q s = 3.68 at small sizes and q l ≈ q s , the overall level of 38 –waviness, and the approximately invariant shape from 1 µ m to 6–30 km agree with resultsfrom the coagulation calculations. Analytic models underestimate the mass lost from de-structive collisions. Coagulation calculations have smaller total mass, dust mass, and dustluminosity than analytic models. These differences between analytic models and numericalcalculations probably stem from the variation of the size of the largest object with time. Inthe coagulation calculations, large reductions in r max shorten the collision time and speedup the evolution.To place the collisional cascade calculations in context with planet formation models,we examine a suite of coagulation calculations where gravitational interactions among thesolids yield an internally consistent e and i as a function of particle size. For swarms with Σ = 10 g cm − , r = 10 n , and n = 0, 1, 2, . . . , 8, it takes 0.1–1 Myr for the largest objects toreach sizes of 2000–4000 km. As these protoplanets grow, their gravity stirs up the e and i forsmaller particles, initiating a collisional cascade. From 0.1–10 Myr, the sizes of the largestobjects that suffer catastrophic collisions grows from r brk ≈ r brk ≈ r .
100 km, the dust luminosity reaches a peak L d /L ⋆ ≈ − × − and thenslowly declines. During this decline, giant impacts between 300–1000 km and larger objectsrelease large amounts of debris, producing pronounced spikes in L d . Eventually, accretionby large objects and destructive collisions among smaller objects eliminate all of the smallerparticles. The dust luminosity then drops to zero.Compared to analytic and numerical models of pure collisional cascades, the planetformation calculations reveal several stark differences. In planet formation calculations, thepersistent growth in r brk involves larger and larger objects in the cascade; r brk is constant inpure cascades. Until the precipitous drop in L d when protoplanets achieve their final masses,the decline of L d in planet formation models is slower in time. Throughout the evolution,these models also produce a higher frequency of spikes in L d . Despite the similar maximum L d , planet formation models typically have a factor of 3–10 less mass in small particles ( r . r & L d /L ⋆ ≪ − would allow more stringent tests of the calculations.We acknowledge a generous allotment of computer time on the NASA ‘discover’ cluster.Comments from M. Geller, G. Kennedy, J. Najita, and an anonymous referee improved ourpresentation. Portions of this project were supported by the NASA Astrophysics Theory and
Origins of Solar Systems programs through grant NNX10AF35G and the
NASA OuterPlanets Program through grant NNX11AM37G.
A. Appendix
To test the algorithms used in
Orchestra , we compare numerical results with analyticsolutions to the coagulation equation (Kenyon & Luu 1998; Kenyon & Bromley 2015) andpublished results from other investigators (Kenyon & Bromley 2001; Bromley & Kenyon2006; Kenyon & Bromley 2008; Bromley & Kenyon 2011a; Kenyon & Bromley 2015). Here,we examine how
Orchestra performs for collisional cascades.
A.1. The Mass Spacing Parameter
The accuracy of coagulation calculations depends on the mass spacing parameter be-tween adjacent mass bins, δ k = m k +1 /m k (e.g., Wetherill 1990; Kenyon & Luu 1998; Kenyon& Bromley 2015). At the start of each calculation, we set the typical mass m k and theboundaries m k − / and m k +1 / of each mass bin. These parameters remain fixed throughouta calculation. The initial average mass within each bin is ¯ m k = M k /N k ; typically ¯ m k = m k .As the calculation proceeds, collisions add and remove material from all bins; the averagemass in a bin ¯ m k and the average physical radius of particles in the bin ¯ r k = (3 ¯ m k / πρ p ) / then change with time.To illustrate how the evolution depends on δ , we consider an annulus centered at a =1 AU. The annulus has a width ∆ a = 0 . a ; thus the inner and outer radii are a in = 0.9 AUand a out = 1.1 AU. We seed the annulus with particles having mass density ρ p = 3 g cm − , e = 0.1, and i = e /2. The annulus has surface density Σ = 10.6 g cm − and total mass M = 0.5 M ⊕ .Fig. 24 shows the time-variation of r max (the radius of the largest object) for calculationswith a mono-disperse set of 100 km objects and five different values of δ . All curves have thesame general shape: after a brief, ∼ yr period where r max is roughly constant, the size 40 –of the largest objects declines monotonically with time. Tracks with larger δ decline faster.Along each track, the evolution consists of a gradual reduction in r max interspersedwith large jumps to smaller r max . Cratering collisions – where somewhat smaller objectsgradually chip away at the mass of larger objects – produce continuous mass loss from thelargest objects. Thus, the average mass in the largest mass bin falls with time. Eventually,this mass falls below the mass boundary between adjacent bins (e.g., ¯ m k < m k − / ). Objectsin bin k are then placed into bin k −
1. Averaging the mass of the ‘old’ objects in bin k − k yields a new average mass ¯ m k − which is smallerthan the average mass of bin k . Thus, the radius (average mass) of the largest object jumpsdownward. Because the spacing of mass bins scales with δ , calculations with larger δ havelarger jumps than those with smaller δ .Although the mass loss rate from the grid is fairly insensitive to δ (see below), the massof the largest object declines faster in calculations with larger δ . Cratering collisions areresponsible for this difference. For all δ , these collisions are rare. Thus, only a few of thelargest objects suffer substantial mass loss from cratering collisions every time step. When δ is small (1.05–1.10), these objects are placed into the next smallest mass bin; the averagemass of the remaining objects in the mass bin is unchanged. When δ is large (1.41–2.00),the amount of mass loss is not sufficient to place objects into the next smallest mass bin;the average mass of all objects in the bin then decreases. As a result, the average mass ofthe largest objects declines faster when δ = 2 than when δ = 1.05.Despite this difference, other aspects of the evolution are fairly insensitive to δ . Fig. 25shows the relative size distributions at 1 Myr. Each curve follows a standard pattern, withan excess (deficit) of particles at 1–10 µ m (10–100 µ m), a wavy pattern at 0.1 mm to 0.1 km,and then a rise at larger sizes. For δ = 1.05–1.41, there is little difference between the curves.When δ = 2, the amplitude of the waves is larger.Fig. 26 shows a similar pattern at 100 Myr. The three curves for δ = 1.05–1.19 arealmost identical. Although the pattern for δ = 1.41 is very similar, the waviness is morepronounced for particles with sizes of 100 µ m. This waviness is even larger for δ = 2.To explore the evolution of the particle flux through the grid, Fig. 27 shows the produc-tion rate for particles with r < µ m. At early times, collisions among the mono-disperse setof 100 km objects produce little debris with sizes smaller than 1 mm. The mass loss rate fromthe grid is then zero. As each calculation proceeds, destructive collisions disperse more andmore small particles. Collisions among these particles yield particles smaller than 1 µ m. Themass loss rate then rises abruptly. The timing of this rise depends on stochastic variationsin the collision rates among large particles. Thus, the timing of the rise is not sensitive to 41 – δ . The size distribution then takes awhile to reach an equilibrium which depends on δ andon the timing of the rise in the mass loss rate. As this equilibrium is established, the massloss rate rises slowly. Afterwards, the mass loss rate declines, ˙ M ∝ t − α with α ≈ δ = 1.05–1.41, the long term time evolution of ˙ M is independent of δ .When δ = 2.0, the decline of ˙ M with time is more rapid at late times.The evolution of the dust luminosity follows the evolution of ˙ M (Fig. 28). Early on,it takes awhile for destructive collisions to produce a large ensemble of small objects. Oncethese objects exist, the dust luminosity rises very rapidly. As the size distribution of smallparticles adjusts to an equilibrium, the luminosity slowly rises to a clear maximum. Oncethe size distribution reaches an equilibrium, the luminosity declines roughly linearly withtime. Calculations with δ = 1.05–1.41 yield nearly identical solutions for the time variationof the dust luminosity. At late times, the luminosity declines more rapidly when δ = 2.00than when δ = 1.05–1.41.To conclude this discussion, we consider how the time evolution of ˙ M changes in a‘hybrid’ model, where we adopt different values of δ for large and small particles. We adopt δ = 2.0 for particles with r . δ = 1.05, 1.10, 1.19, or 1.41 for particles with r & δ varies linearly with particle size. This approach seeksto capture the better accuracy of small δ models with a smaller expense of cpu time.Fig. 29 shows the results. Aside from the magenta curve (the pure δ = 2.00 from Fig. 27),all of the hybrid models follow the time variation of ˙ M established for models with δ = 1.05for all particles sizes. These results demonstrate that we can follow collisional cascades with δ = 2.0 for small particles and δ . A.2. Analytic Model for Wavy Size Distributions
Wavy size distributions are characteristic of coagulation calculations with a low masscutoff (Campo Bagatin et al. 1994; Durda et al. 1998; O’Brien & Greenberg 2003; Kral et al.2013). To extend previous analytic and numerical studies, Wyatt et al. (2011) explore a newsteady state numerical model which shows how the amplitudes and positions of the wavesdepend on the collision velocity. Here, we derive an analytic solution to their model andbriefly illustrate how the properties of the waves depend on Q ⋆D .Following method I of the Wyatt et al. (2011) framework (see their § r min to r max with indices k = 1 to k = N.Reversing the sense of the labeling from Wyatt et al. (2011), the smallest (largest) bin has k = 1 ( k = N). If the mass loss rate through the grid is independent of particle size, the mass 42 – M k contained in a mass bin is then a simple function of the rate of destructive collisions R k summed over all particles with index i and i < k : M k = C R − k , (A1)where C is an arbitrary constant (see also eq. 15 of Wyatt et al. 2011). The sum R k is R k = C i = k X i =1 ǫ ik M i ( r i + r k ) /r i , (A2)where C is another constant. Recalling that Q c is the collision energy per unit mass and Q ⋆D is the binding energy per unit mass, ǫ ik = Q c < Q ⋆D Q c ≥ Q ⋆D (A3)Thus, R k is the rate of collisions which disperse at least half the mass of the colliding pairof particles.The mass conservation equation in eq. A1 has a recursive solution. For k = 1, R k hasonly one term: R k = 4 C M k /r k . Thus, M k = C / C r k . For k ≥ R k is a sum over termswith known M i and one term with M k : R k = 4 C M k /r k + C i = k − X i =1 ǫ ik M i ( r i + r k ) /r i . (A4)Setting R ki equal to the second term in eq. A4 leads to a quadratic equation,4 C M k + R ki M k − C = 0 , (A5)which has only one possible solution with M k > Q ⋆D , we consider modelswith Q ⋆D = constant and Q ⋆D = Q b r β b + Q g ρ p r β g with the standard fragmentation parameters( § v = 3 km s − , the mass density ρ = 3 g cm − ,and δ r = r i +1 /r i = 10 . = 1.023. The center-of-mass collision energy Q c = v c m i m k / m i + m k ) then establishes the ratio Q c /Q ⋆D and ǫ ik for each pair of particles. Recursive solutionof eq. A5 for k = 1 (1 µ m) to k = 1101 (100 km) yields M k for each particle size. Dividing M k by the expected M k ∝ r αk with α = 1/2 for Q ⋆D = constant and α = 0.32 for Q ⋆D ( r ) yieldsthe relative mass in each bin.Fig. 30 compares results for the two models. When Q ⋆D = constant (violet line), therelative mass generally follows the power law. Waves with modest amplitude start at the 43 –small-size cutoff and gradually diminish in amplitude to large sizes. Relative to the powerlaw, there is a clear excess of particles at the smallest sizes: these bins require extra mass togenerate enough collisions to maintain a constant rate of mass loss per bin.When Q ⋆D is a function of particle size, the system follows a power law from 1 µ mto 0.1 km. For larger sizes, the gravity component of the binding energy dominates; theslope of the size distribution changes from 3.68 to 3.0. Although systems where Q ⋆D is afunction of r still display wavy size distributions, the waves have smaller amplitude andshorter wavelength.For both Q ⋆D relations, the shape of the size distribution derived from numerical simula-tions is remarkably close to the analytic prediction. In particular, (i) the amplitudes of thewaves at 1–10 µ m and (ii) the shape of the rise at 0.1–10 km in systems with Q ⋆D ( r ) are verysimilar in the analytic and numerical calculations. In the numerical simulations, crateringproduces a longer wavelength compared to the analytic model. The gradual reduction in thesize of the largest object modifies the shape of the size distribution at the largest sizes. REFERENCES
Absil, O., Defr`ere, D., Coud´e du Foresto, V., et al. 2013, A&A, 555, A104Agnor, C., & Asphaug, E. 2004, ApJ, 613, L157Agnor, C. B., Canup, R. M., & Levison, H. F. 1999, Icarus, 142, 219Arakawa, M., Leliwa-Kopystynski, J., & Maeno, N. 2002, Icarus, 158, 516Asphaug, E., Agnor, C. B., & Williams, Q. 2006, Nature, 439, 155Asphaug, E., & Benz, W. 1996, Icarus, 121, 225Backman, D. E., & Paresce, F. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I.Lunine, 1253–1304Ballering, N. P., Rieke, G. H., & G´asp´ar, A. 2014, ApJ, 793, 57Bandermann, L. W. 1972, MNRAS, 160, 321Barge, P., & Pellat, R. 1991, Icarus, 93, 270Beichman, C. A., Tanner, A., Bryden, G., et al. 2006, ApJ, 639, 1166Beichman, C. A., Lisse, C. M., Tanner, A. M., et al. 2011, ApJ, 743, 85 44 –Benz, W., & Asphaug, E. 1999, Icarus, 142, 5Booth, M., Wyatt, M. C., Morbidelli, A., Moro-Mart´ın, A., & Levison, H. F. 2009, MNRAS,399, 385Bromley, B. C., & Kenyon, S. J. 2006, AJ, 131, 2737—. 2011a, ApJ, 731, 101—. 2011b, ApJ, 735, 29—. 2013, ApJ, 764, 192Burchell, M. J., Leliwa-Kopysty´nski, J., & Arakawa, M. 2005, Icarus, 179, 274Burke, C. J., Christiansen, J. L., Mullally, F., et al. 2015, ApJ, 809, 8Burns, J. A., Lamy, P. L., & Soter, S. 1979, Icarus, 40, 1Campo Bagatin, A., Cellino, A., Davis, D. R., Farinella, P., & Paolicchi, P. 1994,Planet. Space Sci., 42, 1079Carpenter, J. M., Mamajek, E. E., Hillenbrand, L. A., & Meyer, M. R. 2009a, ApJ, 705,1646Carpenter, J. M., Bouwman, J., Mamajek, E. E., et al. 2009b, ApJS, 181, 197Chambers, J. 2008, Icarus, 198, 256Chambers, J. E. 2001, Icarus, 152, 205Chen, C. H., Mamajek, E. E., Bitner, M. A., et al. 2011, ApJ, 738, 122Chiang, E., & Youdin, A. N. 2010, Annual Review of Earth and Planetary Sciences, 38, 493Cloutier, R., Currie, T., Rieke, G. H., et al. 2014, ApJ, 796, 127Currie, T., Kenyon, S. J., Balog, Z., et al. 2008, ApJ, 672, 558Currie, T., Kenyon, S. J., Rieke, G., Balog, Z., & Bromley, B. C. 2007, ApJ, 663, L105Currie, T., Lisse, C. M., Sicilia-Aguilar, A., Rieke, G. H., & Su, K. Y. L. 2011, ApJ, 734,115Davis, D. R., Chapman, C. R., Weidenschilling, S. J., & Greenberg, R. 1985, Icarus, 63, 30 45 –Defr`ere, D., Hinz, P. M., Skemer, A. J., et al. 2015, ApJ, 799, 42Dohnanyi, J. S. 1969, J. Geophys. Res., 74, 2531Dominik, C., & Decin, G. 2003, ApJ, 598, 626Durda, D. D., Greenberg, R., & Jedicke, R. 1998, Icarus, 135, 431Ertel, S., Absil, O., Defr`ere, D., et al. 2014, A&A, 570, A128Fogg, M. J., & Nelson, R. P. 2009, A&A, 498, 575Foreman-Mackey, D., Hogg, D. W., & Morton, T. D. 2014, ApJ, 795, 64Fujiwara, H., Ishihara, D., Onaka, T., et al. 2013, A&A, 550, A45G´asp´ar, A., Psaltis, D., ¨Ozel, F., Rieke, G. H., & Cooney, A. 2012a, ApJ, 749, 14G´asp´ar, A., Psaltis, D., Rieke, G. H., & ¨Ozel, F. 2012b, ApJ, 754, 74G´asp´ar, A., Rieke, G. H., & Balog, Z. 2013, ApJ, 768, 25Genda, H., Kobayashi, H., & Kokubo, E. 2015, ArXiv e-prints, arXiv:1508.00977Genda, H., Kokubo, E., & Ida, S. 2012, ApJ, 744, 137Giblin, I., Davis, D. R., & Ryan, E. V. 2004, Icarus, 171, 487Glaschke, P., Amaro-Seoane, P., & Spurzem, R. 2014, MNRAS, 445, 3620Goldreich, P., Lithwick, Y., & Sari, R. 2004, ARA&A, 42, 549Gorlova, N., Balog, Z., Rieke, G. H., et al. 2007, ApJ, 670, 516Greenberg, R., Hartmann, W. K., Chapman, C. R., & Wacker, J. F. 1978, Icarus, 35, 1Greenberg, R., Weidenschilling, S. J., Chapman, C. R., & Davis, D. R. 1984, Icarus, 59, 87Greenzweig, Y., & Lissauer, J. J. 1990, Icarus, 87, 40Grigorieva, A., Artymowicz, P., & Th´ebault, P. 2007, A&A, 461, 537Grogan, K., Dermott, S. F., & Durda, D. D. 2001, Icarus, 152, 251Hansen, B. M. S., & Murray, N. 2012, ApJ, 751, 158—. 2013, ApJ, 775, 53 46 –Hayashi, C. 1981, Progress of Theoretical Physics Supplement, 70, 35Hellyer, B. 1970, MNRAS, 148, 383Holsapple, K. A. 1994, Planet. Space Sci., 42, 1067Holsapple, K. A., & Housen, K. R. 2007, Icarus, 191, 586Housen, K. R., & Holsapple, K. A. 1999, Icarus, 142, 21Jackson, A. P., & Wyatt, M. C. 2012, MNRAS, 425, 657Johansen, A., Mac Low, M.-M., Lacerda, P., & Bizzarro, M. 2015, Science Advances, 1,15109Kennedy, G. M., & Wyatt, M. C. 2010, MNRAS, 405, 1253—. 2011, MNRAS, 412, 2137—. 2012, MNRAS, 426, 91—. 2013, MNRAS, 433, 2334Kennedy, G. M., Wyatt, M. C., Sibthorpe, B., et al. 2012, MNRAS, 426, 2115Kennedy, G. M., Wyatt, M. C., Su, K. Y. L., & Stansberry, J. A. 2011, MNRAS, 417, 2281Kenyon, S. J., & Bromley, B. C. 2001, AJ, 121, 538—. 2002a, AJ, 123, 1757—. 2002b, ApJ, 577, L35—. 2004a, AJ, 127, 513—. 2004b, ApJ, 602, L133—. 2005, AJ, 130, 269—. 2006, AJ, 131, 1837—. 2008, ApJS, 179, 451—. 2010, ApJS, 188, 242—. 2012, AJ, 143, 63 47 –—. 2014, AJ, 147, 8—. 2015, ApJ, 806, 42Kenyon, S. J., Currie, T., & Bromley, B. C. 2014, ApJ, 786, 70Kenyon, S. J., & Luu, J. X. 1998, AJ, 115, 2136—. 1999, AJ, 118, 1101Kobayashi, H., & Dauphas, N. 2013, Icarus, 225, 122Kobayashi, H., & L¨ohne, T. 2014, MNRAS, 442, 3266Kobayashi, H., & Tanaka, H. 2010, Icarus, 206, 735Kobayashi, H., Tanaka, H., Krivov, A. V., & Inaba, S. 2010, Icarus, 209, 836Kokubo, E., & Ida, S. 1996, Icarus, 123, 180Kokubo, E., Kominami, J., & Ida, S. 2006, ApJ, 642, 1131Kominami, J., & Ida, S. 2002, Icarus, 157, 43—. 2004, Icarus, 167, 231Kral, Q., Th´ebault, P., & Charnoz, S. 2013, A&A, 558, A121Krivov, A. V., L¨ohne, T., & Sremˇcevi´c, M. 2006, A&A, 455, 509Krivov, A. V., Eiroa, C., L¨ohne, T., et al. 2013, ApJ, 772, 32Lee, M. H. 2000, Icarus, 143, 74Leinhardt, Z. M., & Stewart, S. T. 2009, Icarus, 199, 542—. 2012, ApJ, 745, 79Leinhardt, Z. M., Stewart, S. T., & Schultz, P. H. 2008, in The Solar System Beyond Neptune,ed. Barucci, M. A., Boehnhardt, H., Cruikshank, D. P., & Morbidelli, A. (Universityof Arizona Press, Tucson, AZ), 195–211Levison, H. F., Duncan, M. J., & Thommes, E. 2012, AJ, 144, 119Levison, H. F., & Morbidelli, A. 2007, ArXiv Astrophysics e-prints, astro-ph/0701544L¨ohne, T., Krivov, A. V., & Rodmann, J. 2008, ApJ, 673, 1123 48 –L¨ohne, T., Augereau, J.-C., Ertel, S., et al. 2012, A&A, 537, A110Love, S. G., & Ahrens, T. J. 1996, Icarus, 124, 141Luhman, K. L., & Mamajek, E. E. 2012, ApJ, 758, 31Lunine, J. I., O’brien, D. P., Raymond, S. N., et al. 2011, Advanced Science Letters, 4, 325Mamajek, E. E., Meyer, M. R., Hinz, P. M., et al. 2004, ApJ, 612, 496Matthews, B. C., Krivov, A. V., Wyatt, M. C., Bryden, G., & Eiroa, C. 2014, ArXiv e-prints,arXiv:1401.0743Melis, C., Zuckerman, B., Rhee, J. H., & Song, I. 2010, ApJ, 717, L57Melis, C., Zuckerman, B., Rhee, J. H., et al. 2012, Nature, 487, 74Meng, H. Y. A., Rieke, G. H., Su, K. Y. L., et al. 2012, ApJ, 751, L17Meng, H. Y. A., Su, K. Y. L., Rieke, G. H., et al. 2014, Science, 345, 1032—. 2015, ApJ, 805, 77Mennesson, B., Millan-Gabet, R., Serabyn, E., et al. 2014, ApJ, 797, 119Meyer, M. R., Carpenter, J. M., Mamajek, E. E., et al. 2008, ApJ, 673, L181Mo´or, A., Apai, D., Pascucci, I., et al. 2009, ApJ, 700, L25Morbidelli, A., Bottke, W. F., Nesvorn´y, D., & Levison, H. F. 2009, Icarus, 204, 558Mustill, A. J., & Wyatt, M. C. 2009, MNRAS, 399, 1403Nagasawa, M., Lin, D. N. C., & Thommes, E. 2005, ApJ, 635, 578Nesvold, E. R., Kuchner, M. J., Rein, H., & Pan, M. 2013, ApJ, 777, 144O’Brien, D. P., & Greenberg, R. 2003, Icarus, 164, 334O’Brien, D. P., Morbidelli, A., & Levison, H. F. 2006, Icarus, 184, 39Ohtsuki, K. 1992, Icarus, 98, 20—. 1999, Icarus, 137, 152Ohtsuki, K., & Nakagawa, Y. 1988, Progress of Theoretical Physics Supplement, 96, 239 49 –Ohtsuki, K., Stewart, G. R., & Ida, S. 2002, Icarus, 155, 436Olofsson, J., Juh´asz, A., Henning, T., et al. 2012, A&A, 542, A90Ormel, C. W., Dullemond, C. P., & Spaans, M. 2010, Icarus, 210, 507Patel, R. I., Metchev, S. A., & Heinze, A. 2014, ApJS, 212, 10Petigura, E. A., Howard, A. W., & Marcy, G. W. 2013, Proceedings of the National Academyof Science, 110, 19273Plavchan, P., Jura, M., & Lipscy, S. J. 2005, ApJ, 631, 1161Raymond, S. N., Kokubo, E., Morbidelli, A., Morishima, R., & Walsh, K. J. 2014, Protostarsand Planets VI, 595Raymond, S. N., Quinn, T., & Lunine, J. I. 2005, ApJ, 632, 670Raymond, S. N., Armitage, P. J., Moro-Mart´ın, A., et al. 2011, A&A, 530, A62—. 2012, A&A, 541, A11Rhee, J. H., Song, I., & Zuckerman, B. 2008, ApJ, 675, 777Ribas, ´A., Mer´ın, B., Ardila, D. R., & Bouy, H. 2012, A&A, 541, A38Rieke, G. H., Su, K. Y. L., Stansberry, J. A., et al. 2005, ApJ, 620, 1010Rodriguez, D. R., Duchˆene, G., Tom, H., et al. 2015, MNRAS, 449, 3160Rodriguez, D. R., & Zuckerman, B. 2012, ApJ, 745, 147Ryan, E. V., Davis, D. R., & Giblin, I. 1999, Icarus, 142, 56Safronov, V. S. 1969, Evoliutsiia doplanetnogo oblaka. (Evolution of the ProtoplanetaryCloud and Formation of the Earth and Planets, Nauka, Moscow [Translation 1972,NASA TT F-677] (1969.)Schneider, A., Song, I., Melis, C., et al. 2013, ApJ, 777, 78Shannon, A., Wu, Y., & Lithwick, Y. 2015, ApJ, 801, 15Siegler, N., Muzerolle, J., Young, E. T., et al. 2007, ApJ, 654, 580Silverstone, M. D., Meyer, M. R., Mamajek, E. E., et al. 2006, ApJ, 639, 1138 50 –Smith, R., Jeffries, R. D., & Oliveira, J. M. 2011, MNRAS, 411, 2186Smith, R., Wyatt, M. C., & Dent, W. R. F. 2008, A&A, 485, 897Smith, R., Wyatt, M. C., & Haniff, C. A. 2012, MNRAS, 422, 2560Song, I., Zuckerman, B., Weinberger, A. J., & Becklin, E. E. 2005, Nature, 436, 363Spangler, C., Sargent, A. I., Silverstone, M. D., Becklin, E. E., & Zuckerman, B. 2001, ApJ,555, 932Spaute, D., Weidenschilling, S. J., Davis, D. R., & Marzari, F. 1991, Icarus, 92, 147Stark, C. C., & Kuchner, M. J. 2009, ApJ, 707, 543Stauffer, J. R., Rebull, L. M., Carpenter, J., et al. 2005, AJ, 130, 1834Stauffer, J. R., Rebull, L. M., James, D., et al. 2010, ApJ, 719, 1859Stern, S. A., & Colwell, J. E. 1997, AJ, 114, 841Tanaka, H., Inaba, S., & Nakazawa, K. 1996, Icarus, 123, 450Th´ebault, P., & Augereau, J.-C. 2007, A&A, 472, 169Trilling, D. E., Stansberry, J. A., Stapelfeldt, K. R., et al. 2007, ApJ, 658, 1289Trilling, D. E., Bryden, G., Beichman, C. A., et al. 2008, ApJ, 674, 1086Urban, L. E., Rieke, G., Su, K., & Trilling, D. E. 2012, ApJ, 750, 98Weidenschilling, S. J. 1977, Ap&SS, 51, 153—. 1980, Icarus, 44, 172—. 1989, Icarus, 80, 179—. 1994, Nature, 368, 721—. 2010a, ApJ, 722, 1716—. 2010b, ApJ, 722, 1716—. 2011, Icarus, 214, 671Weidenschilling, S. J., Spaute, D., Davis, D. R., Marzari, F., & Ohtsuki, K. 1997, Icarus,128, 429 51 –Weinberger, A. J., Becklin, E. E., Zuckerman, B., & Song, I. 2004, AJ, 127, 2246Wetherill, G. W. 1980, ARA&A, 18, 77—. 1990, Icarus, 88, 336Wetherill, G. W., & Stewart, G. R. 1989, Icarus, 77, 330—. 1993, Icarus, 106, 190Williams, D. R., & Wetherill, G. W. 1994, Icarus, 107, 117Wurz, P. 2012, in Astrophysics and Space Science Library, Vol. 385, Astrophysics and SpaceScience Library, ed. I. Mann, N. Meyer-Vernet, & A. Czechowski, 161Wyatt, M. C. 2008, ARA&A, 46, 339Wyatt, M. C., Booth, M., Payne, M. J., & Churcher, L. J. 2010, MNRAS, 402, 657Wyatt, M. C., Clarke, C. J., & Booth, M. 2011, Celestial Mechanics and Dynamical Astron-omy, 111, 1Wyatt, M. C., & Dent, W. R. F. 2002, MNRAS, 334, 589Wyatt, M. C., Smith, R., Greaves, J. S., et al. 2007a, ApJ, 658, 569Wyatt, M. C., Smith, R., Su, K. Y. L., et al. 2007b, ApJ, 663, 365Youdin, A. N. 2011, ApJ, 742, 38Youdin, A. N., & Kenyon, S. J. 2013, From Disks to Planets, ed. T. D. Oswalt, L. M. French,& P. Kalas (Dordrecht: Springer Science & Business Media), 1Zuckerman, B., Melis, C., Rhee, J. H., Schneider, A., & Song, I. 2012, ApJ, 752, 58Zuckerman, B., Rhee, J. H., Song, I., & Bessell, M. S. 2011, ApJ, 732, 61Zuckerman, B., & Song, I. 2004, ApJ, 603, 738
This preprint was prepared with the AAS L A TEX macros v5.2.
52 –Table 1. List of Variables
Variable Definition a semimajor axis, radial coordinate A d cross-sectional area of particles b d exponent in relation for debris from collisions b l exponent in relation for mass of largest particle in debris e eccentricity e H Hill eccentricity f c fraction of solids converted into small particles ejected by radiation pressure f d detection frequency of dust emission f g gravitational focusing factor F radiation flux emitted by a particle or a star at a specific wavelength h horizontal velocity H vertical scale height i inclination L d luminosity of particles L ⋆ stellar luminosity m particle mass¯ n average mass of particle in a mass bin m l mass of largest particle in debris m l, coefficient in relation for mass of largest particle in debris m max mass of largest particle m min mass of smallest particle M d total mass in particles˙ M mass loss rate M ⋆ stellar mass n exponent of power law or particle number density N particle number N ( > r ) cumulative number of particles larger than rN max number of largest particles q exponent of power law size distribution q c , q r exponents for cumulative and relative cumulative size distribution q d exponent of power law size distribution for debris Q b , Q g coefficients in Q ⋆D relation Q c center of mass collision energy
53 –Table 1—Continued
Variable Definition Q ⋆D collision energy required to eject 50% of the mass r particle radius¯ r average radius of particle in a mass bin r brk radius of particles destroyed in equal mass collisions r max radius of largest particle r min radius of smallest particle r Q particle size with smallest Q ⋆D t time t c collision time T temperature v relative collision velocity v h horizontal velocity v K orbital velocity v z vertical velocity V volume α correction factor for collision time β ratio of the radiation force to the gravitational force β b , β g exponents in Q ⋆D relation δ mass spacing factor δa width of annulus ξ ratio of mass in small particles ( r < ρ particle mass density σ geometric cross sectionΣ surface densityΩ angular velocityNote. — Variables with a subscript ‘0’ refer to initial conditions; e.g., e is the initial eccentricity
54 – Time (yr)10 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y
10 km100 km1000 km 1 AU, MMSN
Fig. 1.— Time evolution of the relative luminosity L d /L ⋆ derived from the analytic modelfor a collisional cascade. Material orbits a 1 M ⊙ star inside a cylindrical annulus with a = 1 AU, δa = 0.2 AU, and Σ = 10 g cm − . Solid curves plot the evolution of L d /L ⋆ forcascades with r max = 10 km (black curve), 100 km (violet curve), and 1000 km (orangecurve). Cascades with smaller r max have larger initial luminosity which declines on muchshorter time scales. 55 – Time (yr)020406080100
1, 0.2, 0.01, 0.2, 1.01, 0.5, 0.759/8, 0.2, 0.0 pd, Q constant D* M a x i m u m R a d i u s ( k m )
1, 0.2, 0.01, 0.2, 1.01, 0.5, 0.759/8, 0.2, 0.0 md, Q constant D* Fig. 2.— Time evolution of r max in collisional cascades with e = 0.1, Q ⋆D = 6 × erg g − and various initial conditions and input parameters. The legend indicates b d , m l, , and b l for each curve. Upper panel: results for calculations with a single initial particle size (‘md’; r = 100 km); after 10 − yr, the maximum particle size rapidly evolves from 100 kmto 10–20 km. Lower panel: results for calculations with a power law initial size distributionhaving q = 3.5 from 1 µ m to 100 km (‘pd’); the largest particles gradually diminish in sizefrom 100 km to 20–30 km. 56 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r Fig. 3.— Time evolution of the relative cumulative size distribution N ( > r ) /r − . for a col-lisional cascade starting with a mono-disperse (‘md’) set of planetesimals with r = 100 km, Q ⋆D = 6 × erg g − , b d = 1, m l, = 0.2, and b l = 0 (as indicated in the upper left corner forthis and subsequent plots). The legend in the lower right corner indicates the evolution timein Myr. The grey line shows the predicted relative size distribution when N ( > r ) ∝ r − . .Aside from the gradual loss of large objects ( r & -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r Fig. 4.— As in Fig. 3 for a cascade with an initial power-law size distribution (‘pd’, q =3.5) from 1 µ m to 100 km. Compared to calculations starting with a mono-disperse set oflarge objects, the amplitude of the wave at r ≈
10 km is smaller. 58 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r md, 1, 0.2, 0.0md, 1, 0.2, 1.0md, 1, 0.5, 0.75md, 9/8, 0.2 0.0 Fig. 5.— Comparison of n rel ( > r ) at 1 Myr for calculations with a mono-disperse set oflarge planetesimals ( r = 100 km) and Q ⋆D = 6 × erg g − . The legend indicates the typeof initial size distribution and values for the input parameters b d , m l, , and b l . Relative toa power-law size distribution with N ( > r ) ∝ r − . , all curves have a similar morphology at1 Myr: an excess of 1–10 µ m particles, a deficit of 10–100 µ m particles, a wavy behaviorfrom 100 µ m to 10 km, an excess of 10–30 km particles, and a deficit of 30–100 km particles.Although the fine details of the shape depend on the input parameters, the shape for eachset of input parameters is independent of time (see also Figs. 3–4). 59 – -1 Time (yr)10 -8 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y
1, 0.2, 0.01, 0.2, 1.01, 0.5, 0.759/8, 0.2, 0.0 Q constantd
Fig. 6.— Evolution of the relative dust luminosity L d /L ⋆ for the calculations in Fig. 2.Starting from a mono-disperse ensemble of planetesimals (thick lines), L d /L ⋆ rises rapidlyfrom roughly zero to ∼ − and then declines roughly linearly with time. When calculationsbegin with a power-law size distribution of planetesimals (thin lines), L d /L ⋆ fluctuates abouta constant as the size distribution reaches an equilibrium and then declines roughly linearlywith time. Although the evolution of L d /L ⋆ at late times depends on b d , m l, , and b l (asindicated in the legend), it is insensitive to the shape of the initial size distribution. 60 – Time (yr)020406080100 M a x i m u m R a d i u s ( k m ) Fig. 7.— Evolution of r max for cascade models with constant Q ⋆D , b d = 1, m l, = 0.2, b l = 1.0, initially mono-disperse (‘md’, upper panel) or power-law (‘pd’, lower panel) sizedistributions, and various ratios v / Q ⋆D as indicated in the legend on the left side of eachpanel. Systems with larger v / Q ⋆D have shorter collision time scales. 61 – Time (yr)10 -8 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y Fig. 8.— Evolution of L d /L ⋆ for cascade models with various ratios v / Q ⋆D as in Fig. 7.Although the peak luminosity is fairly independent of v / Q ⋆D , systems with larger ratiosreach peak luminosity earlier in time. Independent of the initial size distribution, systemswith smaller v / Q ⋆D ratios have larger relative luminosity at late times. 62 – Time (yr)10 M a x i m u m R a d i u s ( k m ) Q (r), 1, 0.2, 1 D* Fig. 9.— Evolution of r max for cascade models with Q ⋆D ( r ), b d = 1, m l, = 0.2, b l = 1,initially mono-disperse (thick lines) or power-law (thin lines) size distributions, and various r . The boundary between growth and destructive evolution is at r ≈
400 km. 63 – Time (yr)10 -8 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y
100 km300 km500 km1000 km
Fig. 10.— Evolution of L d /L ⋆ for cascade models with various r (in km) as indicated inthe legend and mono-disperse (thick lines) or power-law (thin lines) initial size distributions.When the largest objects grow through mergers, originally mono-disperse swarms have min-imal L d for short periods of time. In swarms with an power law initial size distribution, L d is fairly well-correlated with r . 64 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r
1, 0.2, 0.001, 0.2, 1.001, 0.5, 0.759/8, 0.2, 0.00 100 Myr
Fig. 11.— Relative cumulative size distributions at 100 Myr for cascade models with amono-disperse initial size distribution, Q ⋆D ( r ), r = 100 km, and various b d , m l, , and b l asindicated in the legend. Relative to a power-law size distribution with N ( > r ) ∝ r − . , allcurves have a wavy morphology with distinct excesses of particles at r ≈ µ m and r & Q ⋆D at 0.1 km. Although the fine details of the shape depend on the input parameters, the shapefor each set of input parameters is independent of time for t & Time (yr)020406080100 M a x i m u m R a d i u s ( k m ) Fig. 12.— Evolution of r max for cascade models with a mono-disperse initial size distribution, b d = 1, m l, = 0.2, and b l = 1.0 (as indicated in the legend in the upper right corner) andvarious q d as indicated in the legend in the lower left corner. In simulations with larger q d ,the largest objects reach smaller sizes at late times. 66 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r Fig. 13.— Relative cumulative size distributions at 1 Myr for mono-disperse cascade modelswith b d = 1, m l, = 0.2, b l = 0.1 and various q as indicated in the legend. Aside from themorphology of the waviness for r ≈ µ m to 10–100 cm, the size distribution is independentof q . 67 – Time (yr)10 -8 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y Fig. 14.— Evolution of L d /L ⋆ for the mono-disperse cascade models in Figs. 12– 13. Systemswith larger q d reach larger L d and decline more slowly than systems with smaller q d . 68 – Time (yr)10 -7 -6 -5 -4 R e l a t i v e D u s t M a ss c, 1, 0.2, 0c, 1, 0.2, 1v, 1, 0.2, 0v, 1, 0.2, 1Q constant D* Fig. 15.— Time evolution of the relative dust mass ξ for collisional cascades with mono-disperse (thick lines) and power law (thin lines) initial size distributions. The legend indicatesthe prescription for Q ⋆D (‘c’ for constant and ‘v’ for Q ⋆D ( r )) and values for b d , m l, , and b l .On time scales longer than the collision time, the relative dust mass slowly grows with time.Although ξ depends on Q ⋆D , it is fairly independent of the initial size distribution and otherfragmentation parameters. 69 – Time (yr)10 M a ss L o ss R a t e ( g / s ) Analytic1, 0.2, 0.01, 0.2, 1.01, 0.5, 0.75 Q constant D Fig. 16.— Comparison of the mass loss rate from numerical simulations with the analyticmodel. The legend indicates values of b d , m l, , and b l for simulations with a constant Q ⋆D =6 × erg g − . The black lines plot predictions from the analytic model for α = 1 (upperdashed), 0.1 (solid), and 0.01 (lower dashed). Colored lines plot results for calculationswith an initial power law size distribution. Although the numerical estimates match theanalytic prediction for the mass loss rate with α = 0.1 at early times, they decline with timemore rapidly than the analytic result. Occasional spikes in the mass loss rate result fromoccasional collisions of large objects. 70 – Time (yr)10 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y Analytic1, 0.2, 0.01, 0.2, 1.01, 0.5, 0.75 Q constant D Fig. 17.— Comparison of the relative luminosity L d /L ⋆ from numerical simulations with theanalytic model. Despite the random spikes in the mass loss rate (see Fig. 16), the luminositydeclines smoothly with time. 71 – -3 -1 Time (yr)10 -6 -5 -4 -3 -2 -1 M a x i m u m R a d i u s ( k m ) Fig. 18.— Growth of the largest object in a suite of planet formation calculations with Σ =10 g cm − in an annulus with δa = 0.2 AU at a = 1 AU. Labels indicate the initial radius r of a mono-disperse set of objects. Swarms with small r evolve more rapidly than those withlarge r . After 0.1–1 Myr, all swarms produce a few objects with radii exceeding 3000 km.This final radius is independent of r . 72 – -6 -4 -2 Radius (cm)10 -6 -4 -2
10 yr10 yr10 yr10 yr -6 -4 -2
10 yr10 yr10 yr10 yr
Fig. 19.— Evolution of the relative cumulative size distribution ( N ( > r ) /n r − . , with n = 10 ) for a planet formation calculation with Σ = 10 g cm − and r = 10 cm. In eachpanel, the legend indicates the evolution time in yr. To aid in the comparison at late times,both panels include the size distribution at 10 yr (orange points). From 10 yr to 10 yr(upper panel; black, violet, and green points), the largest object grows from 1 km to 100 kmand sweeps up smaller particles in the debris tail. As objects grow from 100 km to 1000 km(0.01–1 Myr, green, orange, and magenta points), collisions produce a prominent debris tailat sizes of 100 cm and smaller. Continued evolution produces a significant deficit of particlesat sizes ranging from 1 cm to 1 km (1–10 Myr; magenta and maroon points). Eventuallycollisions remove most particles smaller than ∼
100 km (100 Myr; black points). 73 – -4 -2 Radius (cm)10 -4 -3 -2 -1 R e l a t i v e C u m u l a t i v e N u m b e r
100 Myr
Fig. 20.— Comparison of relative cumulative size distributions at 100 Myr for calculationswith Σ = 10 g cm − at 1 AU for various r (in cm) as indicated in the legend. At r &
10 km, all starting points lead to the same size distribution. For 0.1 cm to 10 km, the sizedistributions are nearly identical except for a modest offset which depends on the completecollision history but not on r . The relative abundance of particles with sizes smaller than0.1 cm depends on the recent history of collisions between large objects with radii of 10–1000 km. 74 – Time (yr)10 -6 -5 -4 -3 R e l a t i v e L u m i n o s i t y Fig. 21.— Evolution of the dust luminosity for the calculations in Fig. 18. The legendindicates r (in cm) for each calculation. At early times, particle growth dramatically reducesthe ratio of the surface area to the mass for the entire swarm. The dust luminosity declines.Once the collisional cascade begins, the dust luminosity rises rapidly. Despite different risetimes for calculations with different r , the maximum dust luminosity is fairly independentof r . At late times, the evolution consists of spikes from giant impacts superposed on a slowdecline in L d . 75 – Time (yr)10 B r e a k R a d i u s ( k m ) Fig. 22.— Time evolution of r brk for planet formation simulations with various r as indicatedin the legend. Independent of r , r brk gradually increases with time and then reaches a plateauat r brk ≈ Time (yr)10 -6 -5 -4 R e l a t i v e L u m i n o s i t y cc, 100 kmcc, 1000 kmpf, 1 kmpf, 10 cm Fig. 23.— Comparison of L d /L ⋆ for two collisional cascade calculations (black and violetcurves, with r listed in the legend) and two planet formation calculations from Fig. 21(green and orange curves, with r listed in the legend). At early times ( .
10 Myr), threecalculations produce similar L d . When t = 10–100 Myr, L d from the cascade with r = 100 kmis considerably smaller than from the two planet formation models. Despite negligible L d at early times, the cascade with r = 1000 km achieves a dust luminosity similar to theplanet formation calculations after 30–40 Myr. Once t &
100 Myr, planets have removed allof the smaller particles from the annulus; L d for these simulations then lies well below thepredictions of cascade models. 77 – Time (yr)020406080100 M a x i m u m R a d i u s ( k m ) Fig. 24.— Time evolution of r max in collisional cascades with a mono-disperse ensemble of100 km objects and no velocity evolution. The legend in the lower left corner indicates themass spacing factor δ for each calculation. In the upper right corner, the legend indicatesadopted values for b d , m l, , and b l . Continuous cratering collisions produce a gradual declinein r max . Downward jumps in r max occur when mass loss moves the largest objects into thenext smallest mass bin. Calculations with larger δ produce larger and less frequent jumps. 78 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r Fig. 25.— Relative cumulative size distributions N ( > r ) /r − . for collisional cascades at1 Myr. All of the curves have a similar morphology: (i) an excess for particle sizes r ≈ µ m, (ii) a deficit for r ≈ µ m, (iii) a constant part for r ≈ r & r ≈ δ (as indicated in the legend) have larger fluctuations (waves) about a smooth curve. 79 – -4 -2 Radius (cm)10 -1 R e l a t i v e C u m u l a t i v e N u m b e r Fig. 26.— As in Fig. 25 for 100 Myr. At later times, calculations with smaller δ still havesmaller fluctuations relative to calculations with larger δ . 80 – Time (yr)10 M a ss L o ss R a t e ( g / s ) Fig. 27.— Time evolution of the mass loss rate for collisional cascade calculations with amono-disperse swarm of 100 km particles, no velocity evolution, and different δ (as indicatedin the legend in the lower left corner). The legend in the upper right corner indicates adoptedvalues for b d , m l, , and b l . The mass loss rate measures the production rate of objects with r < µ m from collisions between objects with typical radii of 1–100 µ m. Starting from anensemble of 100 km planetesimals, the production rate is initially zero. Destructive collisionsamong large objects eventually produce sufficient numbers of small particles; the mass lossrate then rises abruptly. Within ∼ Time (yr)10 -8 -7 -6 -5 -4 -3 -2 -1 R e l a t i v e L u m i n o s i t y Fig. 28.— As in Fig. 27 for the dust luminosity L d /L ⋆ . Rare destructive collisions among100 km objects gradually produce debris. When the debris contains sufficient numbers of1–100 µ m objects, the dust luminosity rises. For all calculations, the dust luminosity reachesa peak at 0.1 Myr and then declines roughly linearly with time. 82 – Time (yr)10 M a ss L o ss R a t e ( g / s ) Fig. 29.— As in Fig. 27 for hybrid models where δ varies with particle size. As indicated inthe legend, we adopt one δ for r . δ for r & -4 -2 Radius (cm)10 -1 R e l a t i v e M a ss i n B i n Q constant D* Q (r) D* Fig. 30.— Relative mass distributions derived from the analytic model described by eqs. A1–A5 with v = 3 km s − . The relative mass in each bin is M k /r / ( Q ⋆D = constant, violetcurve) or M k /r . ( Q ⋆D ( r ), cyan curve). Systems with Q ⋆D = constant have size distributionswith larger amplitude and longer wavelength waves than systems with Q ⋆D ( r ). When Q ⋆D⋆D