Modeling Collisional Cascades: The Numerical Method

# Modeling Collisional Cascades In Debris Disks: The Numerical Method

András Gáspár Dimitrios Psaltis Feryal Özel George H. Rieke Steward Observatory, University of Arizona, Tucson, AZ 85721
agaspar@as.arizona.edu, dpsaltis@as.arizona.edu, fozel@as.arizona.edu, grieke@as.arizona.edu
Alan Cooney1 Department of Physics, University of Arizona, Tucson, AZ 85721
acooney@physics.arizona.edu
1affiliation: also at Steward Observatory, University of Arizona, Tucson, AZ 85721
###### Abstract

We develop a new numerical algorithm to model collisional cascades in debris disks. Because of the large dynamical range in particle masses, we solve the integro-differential equations describing erosive and catastrophic collisions in a particle-in-a-box approach, while treating the orbital dynamics of the particles in an approximate fashion. We employ a new scheme for describing erosive (cratering) collisions that yields a continuous set of outcomes as a function of colliding masses. We demonstrate the stability and convergence characteristics of our algorithm and compare it with other treatments. We show that incorporating the effects of erosive collisions results in a decay of the particle distribution that is significantly faster than with purely catastrophic collisions.

methods: numerical – circumstellar matter – planetary systems – infrared: stars
journal: The Astrophysical Journal, 749:14 (22pp), 2012 April 10

## 1. Introduction

More than 700 extrasolar planets have been identified to date in over 500 planetary systems. Most of these planets were discovered via radial velocity measurements. As a result, only a handful of them are less than 10 Earth masses; the vast majority are gas giants resembling Jupiter. They are also in extremely close orbits to their host stars, making these systems dramatically different from ours. A large number of additional candidate transiting systems have been found recently with the Kepler mission (Borucki et al., 2010).

In contrast to the great majority of known exoplanet systems, our own solar system has a complex configuration with gas giants at significant distances from their central star and rocky planets/asteroids within the giant planet zone. The direct detection of rocky planets and planetesimals around other stars is only feasible under very rare circumstances. One of the most productive approaches is indirectly via the thermal emission of their planetary debris dust belts. Ever since the discoveries with IRAS (Aumann et al., 1984; Backman & Paresce, 1993), we know that extrasolar systems may harbor disks of dust/debris that are generated by planetesimal collisions. The dust reprocesses the stellar light and emits it as thermal radiation in the 10-1000 wavelength range. A prototypical example of such a system is Fomalhaut, where a planet is shepherding the star’s debris disk resolved in both scattered light (Kalas et al., 2008) and in infrared emission (Holland et al., 2003; Stapelfeldt et al., 2004; Marsh et al., 2005). Debris disks highlight the constituents of planetary systems that are many to hundreds of AU away from their stars.

With the launch of the Spitzer Space Telescope, many observations have been obtained to detect and possibly to resolve debris disks in the infrared regime. Debris disks have been probed around all types of stars, both in stellar clusters and in the field. These observations showed that even though debris disks are common around stars of all spectral types, they are more likely to be detected in the earlier stages of stellar evolution (Wyatt, 2008). We have also learned that debris disks may be located close to or far from their central stars, that there are systems with multiple debris rings (such as our solar system), and that there can be wide varieties of mineralogical compositions within the disks (Carpenter et al., 2009). Debris disk studies are now a major component of the Herschel observing program (Matthews et al., 2010; Eiroa et al., 2010), which will provide substantial advances in our understanding of their outer zones.

Interpreting these results demands theoretical insights in a variety of areas. For example, attempts have been made to understand the evolution of debris disks as a function of stellar type by studying them in stellar clusters of different ages. As concluded in Gáspár et al. (2009), solar-type stars in the field (Beichman et al., 2006; Trilling et al., 2008; Carpenter et al., 2008) and in clusters (Gorlova et al., 2006, 2007; Siegler et al., 2007) may show a faster decay trend compared to that observed for earlier-type stars (Rieke et al., 2005; Su et al., 2006), although the difference is subtle and needs confirmation. The decay trends of the fractional luminosity () show a large range in values. Spangler et al. (2001) find a decay when fitting ISO/IRAS data, while Greaves & Wyatt (2003) get a much shallower decay . The majority of surveys however find a decay (Liu et al., 2004; Moór et al., 2006; Rieke et al., 2005). A better theoretical understanding is needed to sort out these results and to provide testable hypotheses that can be compared with the observations.

Only a handful of debris disks have been resolved; for the majority, we only know the integrated infrared excess emission. Finding the underlying spatial distribution of the debris in these disks is not straightforward, as any spectral energy distribution (SED) can be modeled with a degenerate set of debris rings at different distances. Although much of the uncertainty is associated with the optical constants of the grains, another under-appreciated issue is the grain size distribution. Collisional models can reduce the number of free parameters in the SED models by determining the stable size distribution of particles in the disks.

Observations of resolved debris disks also have raised questions that can best be addressed by theoretical models. For example, Spitzer MIPS images have shown a significant extended halo of dust around Vega (Su et al., 2005), both at 24 and 70 . Initial calculations hypothesized the halo around Vega to be a result of a high outflow of dust due to radiation pressure from a recent high-mass collisional event (Su et al., 2005), while Müller et al. (2010) model it as a result of weakly bound particles on highly eccentric orbits. Further modeling and deep observations of additional systems will help distinguish these two possibilities.

In this paper, we describe a new algorithm for modeling debris disks, in which we refine the physics and numerical methods used in collisional cascade models. In §2, we briefly outline previous models and introduce the basics of our algorithm. In §3 we detail our numerical methods, followed in §4 by our approach for including simplified dynamics. In the last section we compare our numerical algorithm to previous ones and discuss in detail the differences between the codes and the effects those differences have on the outcome of the collisional cascades. We also supplement our paper with an extended appendix that covers the numerical methods and the verification tests of our code.

## 2. The physical and numerical challenges of modeling debris disks

Collisional cascades have been studied both analytically and using collisional integro-differential numerical models. The classic analytic models of Dohnanyi (1969), Hellyer (1970), and Bandermann (1972) took into account both erosive and catastrophic collisional outcomes, assumed a material strength that was independent of the particle mass, a particle mass distribution with no cut-offs, and a constant interaction velocity. These models yielded steady state power-law mass distribution indices of -11/6. This result was in general agreement with the measured size distribution of asteroids in the solar system. More recently, analytic models by Dominik & Decin (2003) and Wyatt et al. (2007) showed that the fractional infrared luminosity in a collision-dominated steady-state system decays following a power-law. As introduced later in this section, Löhne et al. (2008) find a different value for this decay timescale. Wyatt et al. (2007) also derived a maximum mass and fractional luminosity as a function of age and distance from the central star, which they then used to classify systems with possible recent transient events.

However, numerical models are needed to expand on these results. In the particular case of our Solar System, sophisticated numerical models were developed to track the evolution of the largest asteroids (Greenberg et al., 1978). They have been further improved to reproduce the observed wavy structure in the size distribution at the very highest masses (e.g., O’Brien & Greenberg, 2005; Bottke et al., 2005). These models yield power-law distributions that deviate from the classic solution of Dohnanyi (1969), with certain regions steeper than it and others shallower. Using a steeper or shallower distribution and extrapolating it to dust sizes can result in substantial offsets in the number of particles and thus in the infrared emission originating from them for a given planetesimal mass. Conversely, the particle size distribution affects the underlying disk mass calculated from the observed infrared emission.

A complete numerical model of collisional cascades would follow outcomes from all types of collisions, include a kinematic description of the system, incorporate coagulation below certain thresholds, and do all this with high numerical fidelity. Although such a model has not yet been built because of its complexity, there are a number of approaches in the literature that model collisional cascades down to particles of micron size, each with distinctive strengths and weaknesses.

The collisional code ACE has been used in many studies (Krivov et al., 2000, 2005, 2006, 2008; Löhne et al., 2008; Müller et al., 2010). It follows the evolution of the particle size distributions as well as the spatial distribution of the dust in debris disks. The code initially only accounted for collisions resulting in catastrophic outcomes, while the latest version (Krivov et al., 2008; Müller et al., 2010) includes erosive (cratering) events as well. The collisional outcome prescriptions are based on the Dohnanyi (1969) particle-in-a-box model, but with a more elaborate description of material strengths in collision outcomes as well as the radiation force blowout. A strength of the code is that it calculates the dynamical evolution of the systems, as well. Since following the dynamical evolution of a system makes large demands on computer memory space and CPU speed, the code can only model the size distribution with a low number of mass grid points; it originally used a first order Euler Ordinary Differential Equation (ODE) solving algorithm, but has been modified to include a more precise one (A. Krivov, priv. comm.). Krivov et al. (2006) and Müller et al. (2010) applied this algorithm to debris disks in general and to the specific example of the Vega system. They followed the orbital paths of fragments and placed special emphasis on radiation effects. Löhne et al. (2008) modeled debris disk evolution around solar-type stars and found, both with analytic and numerical analysis, that the majority of physical quantities, such as the mass and the infrared luminosity, decrease with time as to . This is in contrast with the observed decay found by some obervations (Liu et al., 2004; Moór et al., 2006; Rieke et al., 2005). However, the population synthesis verification tests in Löhne et al. (2008) yield good agreement with the latest Spitzer observations.

Thébault et al. (2003, 2007) study the evolution of extended debris disks with a particle-in-a-box algorithm. They include both catastrophic and erosive collisions and employ resolution and numerical methods similar to the ones implemented in the ACE code. They model the extended disk structure by dividing the disk into separate, but interacting rings. Their model does not include dynamics.

Campo Bagatin et al. (1994) showed that a series of wave patterns is produced in the mass distribution of particles when a low-mass cutoff is enforced, such as in the case of a radiation force blowout limit. This signature is produced by the other debris disk numerical models as well (Thébault et al., 2003; Thébault & Augereau, 2007; Krivov et al., 2006; Löhne et al., 2008; Wyatt et al., 2011). However, the conditions under which these waves are produced have not been completely analyzed. Wyatt et al. (2011) do show that the amplitude and wavelength of the waves is collisional velocity dependent. Such strong features in the particle size distribution are not observed in the dust collected within the solar system. The interplanetary dust flux model of Grun et al. (1985), which used in situ satellite measurements of the micro-meteroid flux in the solar system, and the terrestrial particle flux measurements of the LDEF satellite (Love & Brownlee, 1993) only show a single peak at in the dust distribution. However, these measurements detected particles that were brought inward from the outer parts of the solar system via Poynting-Robertson drag and particles removed from the inner parts of the solar system via radiation force blowout. Their results are reflections of more than a single parent distribution and of multiple physical effects.

The dynamical code of Kuchner & Stark (2010) models the evolution and 3D structure of the Kuiper belt, with a Monte Carlo algorithm and a simple treatment of particle collisions. Their models predict that grain-grain collisions are important even in a low density debris ring such as our Kuiper belt.

Kuchner & Stark (2010) and Dullemond & Dominik (2005) both emphasize the strong effects of fragmentation in their models. Müller et al. (2010) point out that including erosive (cratering) events is necessary for their models to reproduce the observed surface brightness profiles of Vega. Thébault & Augereau (2007) also show that a complete collisional treatment will result in significant deviations from the classic power-law solution.

Our goal is to set up a numerical model that places special emphasis on investigating these issues. Our new empirical description of collisional outcomes avoids discontinuities between erosive and catastrophic collisions and thus enables a more stable and accurate calculation. We also solve the full scattering integral, thus ensuring mass conservation and the propagation of the largest remnants of collision outcomes. Finally, we use second order integration and fourth order ODE solving methods to improve the numerical accuracy. Below, we outline the physical and numerical techniques we will employ. In the Appendix, we present verification tests of these treatments.

### 2.1. Collisional outcomes

In collision theory, two types of outcomes are generally distinguished: catastrophic and erosive (the latter also known as cratering). For an erosive collision (EC), the collisional energy is relatively small, resulting in one big fragment whose mass is close to the original target mass. In the case of a catastrophic collision (CC), both colliding bodies are completely destroyed.

We illustrate these outcomes in Figure 1. In the first panel, we plot the outcome distribution of an erosive collision, where there is a single large fragment and a distribution of dust at much lower masses. The fragment is over half the mass of the original target mass . The redistributed mass is equal to the cratered mass plus the projectile mass. The largest fragment in the distribution, , is arbitrarily set to be 20% of the cratered mass. In §3.3 we elaborate on the validity of this arbitrary value.

In the second panel of Figure 1, we plot the outcome at the boundary case between catastrophic and erosive collisions, where the single largest fragment is exactly half of the original target mass . The redistributed mass is equal to the other half of the target mass plus the projectile mass. The largest fragment in the distribution, , is arbitrarily set to be 20% of the cratered mass here as well (10% of the target mass).

Finally, in the third panel of Figure 1, we plot the outcome of a super-catastrophic collision, where the target and projectile masses are equal. The mass of the single largest fragment is given by the relation of Fujiwara et al. (1977). The redistributed mass is equal to plus the projectile mass. The largest fragment in the distribution, , is arbitrarily set to be at .

In reality, there is no strict boundary between catastrophic and erosive collisions (Holsapple et al., 2002). The outcomes between these two extreme scenarios should be continuous. In laboratory experiments, however, it is easier to test the extreme outcomes. In our model, we use the laboratory experiments to describe the extreme solutions and connect them with simple interpolations throughout the parameter space. We revise the currently used models to include an fragment for both erosive and catastrophic collisions as a separate new gain term. In this treatment, the placement of the fragments is grid size independent, further improving precision and guaranteeing the accurate downward propagation of these fragments. We are able to express the loss term in a much simpler form, including collisions from both regimes. Previous models only included a full loss term for catastrophic collisions and removed fractions of particles for erosive collisions.

The slope of the power-law particle redistribution has been studied extensively. Dohnanyi (1969) used a single power-law value from the largest mass to the smallest. Later experiments have shown that a double (or even a triple) power-law distribution is a more likely outcome (see, e.g., Davis & Ryan, 1990). This has led to the widespread use of a double power-law for the redistribution in numerous collision models. We conducted numerical tests that demonstrated that there is negligible difference in using a wide range of slopes with a single power-law. The fact that the outcome does not depend on the bimodality of the redistribution has also been proven by Thébault et al. (2003). Therefore, we have used the simplest method of redistributing the fragmented particles with a single power-law slope from the second largest fragment downwards, scaled to conserve mass. As a nominal value, we will use -11/6 as the redistribution slope. This is close to the initial value of -1.8 used by Dohnanyi (1969).

### 2.2. Incorporating the complete redistribution integral

The classic solution to the collisional evolution of an asteroid system involves solving the Smoluchowski (1916) integro-differential equation. This was first done by Dohnanyi (1969). Because erosive collisions remove only a small part of the target mass in a collision, Dohnanyi (1969) expressed the erosive removal term in a differential form. This is not appropriate for our case. Our system has well defined boundaries; thus a continuity equation cannot be used. The locality of the collisional outcomes in phase space is not certain either.

Therefore, to solve the Smoluchowski equation for the problem at hand, we need to solve the full scattering integral. This is complicated numerically, as the integrations must extend over the entire dynamical range of orders of magnitude in mass. To be able to perform accurate integrations over such a large interval, we need to use a large number of grid points and sophisticated numerical methods. To achieve this in a reasonable time, we chose to drop the radial dependence of the various quantities and to solve the equations under a “particle-in-a-box” approximation. With this approach, we lose radial and velocity information but gain accuracy.

### 2.3. The effect of radiation forces

Poynting-Robertson drag can be an effective form of removing particles from the disk, so we include it in our model. However, the strongest and most dominating radiation effect is the removal of particles via radial radiation forces. These act on orbital timescales and can remove or place particles on extremely eccentric orbits. This gives rise to the challenge of incorporating a radial dependent removal term into a particle-in-a-box model that does not carry radial information. In §3.1.1. and §3.1.2, we discuss our approach for incorporating these radiation effects.

Stellar wind drag is an important dust removal effect for late type stars such as in the case of AU Mic (Augereau & Beust, 2006; Strubbe & Chiang, 2006), an M1 spectral-type star. We concentrate on modeling debris disks in early- and solar-type systems, so we chose to neglect the effects of stellar wind drag. We do not take into account the Yarkovsky effect either, as it is small in high-density debris disks compared to the other radiation effects.

## 3. The collisional model

We now discuss our collisional code (CODE-M - COllisional Disk Evolution Model), which solves the system of integro-differential equations that describe the evolution of the number densities of particles of different masses. The code includes outcomes from erosive (cratering) collisions and catastrophic collisions and qualitatively follows the effects originating from radiation forces and Poynting-Robertson drag.

Our system-dependent parameters are: the spectral-type of the central star (which defines the stellar mass and the magnitude of the radiation effects it will have on the particles), the minimum and maximum particle masses ( and , respectively), the radius, width, and height of the debris ring (, , and , respectively), the total mass within the ring (), and the slope of the initial size distribution of the particles (). We estimate the total volume of the narrow ring, , as

 V=2πhRΔR, (1)

which together with defines a mass density.

### 3.1. The evolution equation

In general, the change in the differential number density at any given time for a particle of mass is given by (Smoluchowski, 1916)

 ddtn(m,t)=TPRD+Tcoll, (2)

where is the Poynting-Robertson drag (PRD) term and is the sum of the collisional terms. We define the differential number density of particles such that

 N(t)=∫n(m,t)dm (3)

is the time-dependent total number density of particles within the ring.

Effects such as radiation force blowout and Poynting-Robertson drag are able to deplete the low-mass end of the distribution, which in turn alters the evolution of the disk and more importantly, its infrared signature. Because we do not follow the radial profile of the various debris disk quantities in our algorithm, we can only capture the effects of radiation forces in a simplified way.

#### 3.1.1 Poynting-Robertson drag term

A complete analysis of the effects of the Poynting-Robertson drag is given by Burns et al. (1979), who correct many errors made in previous work. This effect causes the particles to slow in their orbit and follow an inward spiral. Burns et al. (1979) show that the change in the orbital distance can be written as

 dR(m)dt=−2GM∗β(m)cR, (4)

where G is the gravitational constant, c is the speed of light, and is the ratio of radiation to gravitational force experienced by a particle of mass .

We calculate the values as a function of the particle masses, optical constants, and the spectral type of the central star following Gáspár et al. (2008). For the calculations we assume a silicate composition for the particles and a bulk density of g cm. We show the calculated values for a few different spectral type stars in Figure 2.

We use equation (4) to derive an approximate term that captures the effect of Poynting-Robertson drag as (see eq. [2])

 TPRD=−n(m,t)τPRD(m), (5)

where

 τPRD(m)=c2GM∗β(m)RΔR. (6)

The mass dependence of the timescale comes from the mass dependence of the parameter . In principle, once a particle is removed from our collisional system it still radiates in the IR; it just does not take part in the collisional cascade. We keep track of the removal rate of these particles, but do not follow the total amount removed or their infrared emission.

The effects of the radiation force blowout are incorporated in our code with the simplified dynamics treatment introduced in §4, and not by the inclusion of a separate term in the differential equation as are the effects of Poynting-Robertson drag. Removing a particle from the collisional system via radiation force blowout requires roughly an orbital timescale

 τRFB=2π√R3GM∗. (7)

As we will show in §4, under our assumptions a newly created particle of mass gets removed via radiation force blowout if and is unaffected by radiation forces when .

Although the radiation force blowout timescale is not used in our code, in Figure 3 we compare it to the Poynting-Robertson drag timescale around an A0 spectral-type star, assuming a disk width-to-radius ratio of 0.1. The plot shows that within reasonable disk radii estimates, radiation force blowout will always dominate in the domain, while outside of it Poynting-Robertson drag will be the stronger effect. Whether Poynting-Robertson drag is an effective form of removal in the domain depends on the number density of particles in the ring, i.e., the collisional timescale of the system (Wyatt, 2005). The outcomes are similar for all spectral type stars and for realistic values.

### 3.2. The collisional term

The probability of a collision between particles is a function of their number densities and their collisional cross section. We express the collisional cross section for particles of mass and as

 σ(m,m′)=π[r(m)+r(m′)]2, (8)

where is the radius of each particle. We express the differential rate of collisions between the two masses as

 P(m,m′;t) = n(m,t)n(m′,t)Vσ(m,m′) (9) = n(m,t)n(m′,t)Vπ[r(m)+r(m′)]2 = κn(m,t)n(m′,t)Vπ(m13+m′13)2

where is their characteristic collisional velocity, and are the differential number densities for particles of mass and ,

 κ≡(34πρ)23, (10)

and is the bulk mass density of the particles. The number densities in the problem are naturally all time dependent. However, for brevity, hereafter we drop the time dependence from our notation.

The decrease or increase in number density at a certain mass will be determined by three separate events: the removal of particles caused by their interaction with all other particles, the addition of the particles from the interactions of other particles (see §2.1 and Figure 1), and the addition of particles from the redistribution of smaller fragments originating from collisions of other particles.

We express the first event, which describes the removal of particles, as

 ddtn(m)∣∣∣rem=−mmax∫mmindm′P(m,m′). (11)

We completely remove all particles from all grid points if they take part in a collision, even if they are the target objects in erosive collisions.

The second event to be described is the addition of the large fragments. To this end, we need to calculate the mass that will produce a particle of mass when interacting with a particle of mass . We achieve this with a root finding algorithm from the collisional equations presented in the following sections and calculate it only once in the beginning of each run. In equation form,

 ddtn(m)∣∣∣m≡X(m′,M)=μX(m)∫mmindm′P(M,m′). (12)

The lower limit of the integration is the minimum mass in the distribution. We denote the largest mass that can create a particle of mass as the fragment as . Its value can also be calculated via root finding algorithms and has to be calculated for each value of once in the beginning of each run (see Figure 15 in the Appendix).

These first two integrals may catastrophically cancel, meaning that the difference between the two terms may be significantly smaller than the absolute value of each, causing the former to be artificially set to zero when evaluated numerically. It is therefore useful to combine these terms into a single integral in a way that will lessen the probability of catastrophic cancellation:

 TI(m)=−Vπκ{μX(m)∫mmindm′n(m′)[n(m) (m13+m′13)2−n(M)(M13+m′13)2] +mmax∫μX(m)dm′n(m′)n(m)(m13+m′13)2}, (13)

Unfortunately can still suffer from catastrophic cancellation (when is much smaller than , and by definition only for the first integral). We overcome this issue by employing a Taylor-series expanded form of , as given in the Appendix.

The third event in the collisional term is the addition of the power-law fragments back to the distribution. The description of this process is quite simple; however, its precise calculation is not. We write in general

 TII(m)= ∫mmaxmmindμ∫mmaxμdMP(μ,M)A(μ,M)× (14) H[Y(μ,M)−m]m−γ,

where is the projectile mass, is the target mass, is the scaling of the power-law distribution, and is the Heaviside function. The total redistributed mass is

 Mredist.(μ,M)=∫Y(μ,M)0A(μ,M)m−γ+1dm, (15)

where is the largest fragment within the redistribution (i.e., the second largest fragment in the collision, after X; see §2.1). This gives the scaling factor

 A(μ,M)=(2−γ)Mredist.(μ,M)Y2−γ(μ,M). (16)

The precision of this integration depends strongly on the resolution of our grid points, due to the integration limits set by the Heaviside function. We discuss in detail the integration methods we used in the Appendix.

### 3.3. Collision outcomes

The collisional equations can be integrated if the values of , , and are known as a function of the colliding masses. Their values are strongly dependent on the outcome of the collision they originate from, which is determined by the energies of the colliding parent bodies. We show the domains of erosive, interpolated erosive (explained later in the section), and catastrophic collisions as a function of the colliding body masses in Figure 4 for collisional velocities of , 1.0 and 1.5 km s. We introduce the method for calculating collisional velocities from orbital velocities in §4, but it is a good general approximation that the collisional velocity is roughly an order of magnitude smaller than the orbital velocity. These collisional velocities will then correspond to debris ring radii of 100, 25 and 10 AU around an A spectral type star, respectively.

A collision is considered to be catastrophic if

 Q(μ,M)impact≡μV22M≥Q∗(M), (17)

where is the dispersion strength parameter of the target mass , is the projectile mass, and V is the relative velocity of the projectile compared to the parent ring (§3.1). We use the dispersion strength description of Benz & Asphaug (1999) and discuss our choice in the Appendix. Note that, in a more accurate treatment, we would redistribute the relative kinetic energy to both masses and not just to the target mass (i.e., divide by instead of ). We are, however, using the original definition of (as opposed to using the relative kinetic energy) because the values that we will be comparing it to were defined the same way (Benz & Asphaug, 1999) and this definition makes the problem more tractable numerically. We note that some work has indicated that the tensile strength curve itself may be collision velocity dependent (Benz & Asphaug, 1999; Holsapple et al., 2002; Stewart & Leinhardt, 2009), which we currently do not take into account.

In catastrophic collisions both particles are completely destroyed. Based on experimental evidence (Fujiwara et al., 1977; Matsui et al., 1984; Takagi et al., 1984; Holsapple et al., 2002), we will assume that apart from the largest fragment , the total mass is redistributed as a power-law distribution of particles up to a mass that we denote as . We calculate the largest single mass created using the relation (Fujiwara et al., 1977)

 X(μ,M)=M12[μV22MQ∗(M)]−βX. (18)

At the separatrix between catastrophic and erosive collisions, , and , which is exactly what we expect. The factor is measured to be 1.24 by Fujiwara et al. (1977) and this is the fiducial value that we use. Some experiments have shown that the shape and material of the target have an effect on the exact value of (Matsui et al., 1984; Takagi et al., 1984). We will elaborate on the effects of varying in an upcoming paper.

The second largest fragment, , is always a fraction () of the cratered mass, , in the erosive collision domain up to the erosive/catastrophic collision boundary. In the catastrophic collision domain, is a fraction () of the fragment. We interpolate from its value defined by at the separatrix (where , as ) to a specified value at the super catastrophic collision case of as

 fX=exp⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩ln(fY)+ln(fmaxXfY)ln[μ2Q∗(M)MV−2]ln[M2Q∗(M)MV−2]⎫⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪⎭. (19)

Our fiducial values for these fractions are and . We express the remaining mass in catastrophic collisions as

 Mredist(μ,M)=μ+M−X(μ,M), (20)

which is redistributed as a large number of smaller particles.

Erosive collisions are more complicated and less well understood. A collision will be erosive if

 Q(μ,M)impact≡μV22M

As described in §3, an erosive collision will result in a single large fragment, which will be a remnant of the target body, and a distribution of smaller particles. We use the formula of Koschny & Grün (2001a, b), i.e.,

 Mcr=α[μV22]b, (22)

to calculate the total mass cratered from the target, where and are constants, with fiducial values of and . This formula is only valid for small cratered masses; it can lead to artificially high values for the cratered masses (much larger than the target mass) even in the erosive collision domain. When the cratered mass given by this formula is larger than an arbitrarily set fraction of the target mass, we use the following interpolation formula

 Mcr=M×exp⎧⎪ ⎪⎨⎪ ⎪⎩ln(fM)+ln(0.5fM)ln(μV22M/Ql)ln[Q∗D(M)/Ql]⎫⎪ ⎪⎬⎪ ⎪⎭, (23)

where

 Ql=(fMαM1−b)1/b. (24)

We choose an arbitrary fiducial value for of .

In erosive collisions, the single large fragment is expressed as

 X=M−Mcr. (25)

As defined before, the largest fragment of the redistributed mass is a fraction of the cratered mass

 Y(μ,M)=fYMcr(μ,M), (26)

while the redistributed mass is

 Mredist(μ,M)=μ+Mcr. (27)

Thus, the , , and parameters can be summarized as

 X(μ,M) = ⎧⎪⎨⎪⎩M12[μV22MQ∗(M)]−βXin CCM−Mcr(μ,M)in EC (28) Y(μ,M) = {fX(μ,M)X(μ,M)~{}~{}in CCfYMcr(μ,M)~{}~{}in EC (29) Mredist(μ,M) = {μ+M−X(μ,M)~{}~{}in CCμ+Mcr(μ,M)~{}~{}in EC (30)

The is redistributed as a large number of smaller particles in both collision types, with a slope of and a scaling given by equation (16). We choose a redistribution slope of , which is a value close to that given by experimental results (Davis & Ryan, 1990) and is the same as used by Dohnanyi (1969).

We give a list of the variable collisional parameters of our model and their fiducial values in Table 1.

### 3.4. The initial distribution and fiducial parameters

We use the Dohnanyi (1969) steady-state solution of as our initial distribution, where is the slope of the initial distribution and yields an initial number density of , where C is an appropriate scaling constant for the distribution. The exact value of this slope is unknown for all real systems. Fortunately, the convergent solutions and the timescales of reaching a convergent solution are fairly insensitive to this value.

## 4. Simplified Dynamics

For the smallest particles, which we are particularly interested in modeling, radiation forces lead to effects such as reduced collisional probabilities in thin ring disks and increased collisional velocities in extended disks. In this section, we describe our approximate treatment of these effects.

The radiation originating from the central star effectively modifies the mass of the star seen by the particles; the orbits themselves remain conic sections. The effective mass of the star is decreased by a factor of , where is defined in §3.1.1. If

 β(m)≥12 (31)

a newly created particle is generally put on a hyperbolic orbit, which we take as the requirement for radiation force blowout to occur (Kresak, 1976; Burns et al., 1979, and references therein).

The effects of dynamical evolution on the collisional cascade can be traced to the eccentricity of the orbits. To follow the orbital path of a dust grain that has been created in a collision, we assume that the parent bodies were on circular orbits at radius and had . We also assume that the produced grain will be created with very small relative velocity with respect to the parent bodies. Under these conditions, it can be shown that the semi-major axis of the acquired orbit will be

 a(m)=1−β(m)1−2β(m)R. (32)

At the semi-major axis becomes infinite, while at it is equal to the semi-major axis of the colliding particles’ original orbit. The eccentricity () of the orbit can be determined from the fact that the periapsis will equal the original orbital distance

 a(m)[1−eβ(m)]=R, (33)

yielding

 eβ(m)={β(m)/[1−β(m)]if β≤0.5>1if β>0.5. (34)

At , the eccentricity equals 1, and at , it equals zero, consistent with our expectations. Similar derivations can be found elsewhere (e.g., Harwit, 1963; Kresak, 1976). A particle on an eccentric orbit will have a modified probability of interaction with other particles in the parent ring, which we address in §4.2.

### 4.1. Collisional velocities

Lissauer & Stewart (1993) give the velocity of a planetesimal relative to the other planetesimals in the swarm (i.e., the collisional velocity), averaged over an epicycle and over a vertical oscillation as

 V=vorb√54(e2)2+(i2)2, (35)

where is the maximum eccentricity and is the maximum inclination in the system. This equation is valid for a swarm of particles in Rayleigh distributed equilibrium. This condition is true for a system in quasi-collisional equilibrium. We use this equation to estimate the collisional velocity of all particles, setting

 e=ΔR2R (36)

and

 i=h2R. (37)

The smallest particles that are in highly eccentric orbits will have varying velocities along their trajectories. However, when at their periapsis, they will have their original orbital velocities, as by definition they are on eccentric orbits due to their original periapsis velocity. Because of this, in our simplified dynamical treatment we only use a single collisional velocity for all particles, which is described by equation (35).

### 4.2. Reduced collisional probabilities of β critical particles

Particles with less than 0.5, but which are still non-zero, called critical particles, are thought to produce halos around debris disks via the highly eccentric orbits radiation forces place them on (Thébault & Wu, 2008; Müller et al., 2010).

For a particle to go into an eccentric orbit, it must acquire a radial velocity component that is different than zero. In collisions, fragments will be ejected in all directions with a certain velocity distribution. Since the smallest fragments will tend to escape with the highest velocities (e.g., Jutzi et al., 2010), it is a fair question to ask whether thermalization of velocity vectors and their high values is a stronger effect in placing dust particles on higher eccentricity orbits compared to radiation effects that reduce the effective stellar mass.

To answer this question, we need to examine the origin of the particles that contribute to the increase of the differential density in each mass grid point. We calculate only integrating in space, thus calculating the rate of increase of the differential number density of particles with mass that originate from collisions with targets of mass . Our calculations show that there is a pronounced peak in roughly on the same scale (at most one order of magnitude higher) as itself. That is, most particles originate from targets larger in size than the particle itself. The results of Jutzi et al. (2010) clearly show that the velocities acquired by collision fragments at 1/3 sizes are more than an order of magnitude lower than the collisional velocities, meaning that, in the most extreme case, a fragment will receive up to a few tenths of a km s radial velocity compared to its 10-30 km s orbital velocity. We can thus safely say that particles that are created with values similar to tend to be placed on eccentric orbits by the radiation forces rather than being dispersed. These orbits will extend out from the initial debris disk ring, preventing the particles from being destroyed and from them creating other particles.

Our approach to calculating a weighting factor for each particle mass, determined by the fraction of its orbital period it spends in the parent ring is similar to that of Thébault & Wu (2008). The orbital time of a particle in an elliptical orbit as a function of its distance from the center of mass is (Taff, 1985)

 cos−1(a−laeβ)−eβ√1−(a−laeβ)2=(t−t0)√GMa3, (38)

where is the initial time at periapsis, is its distance at time from the center of mass, and we omit the dependences of and for clarity. We estimate the semi-major axis as

 a=R−ΔR/21−eβ, (39)

and we calculate the time needed for a particle to reach the outer edge of the disk at . Dividing by the half of the orbital period gives the weighting factor for each mass as

 w = 1π{cos−1[ΔR(2−eβ)−2Reβeβ(ΔR−2R)] (40) −2 ⎷ΔR(eβ−1)(ΔR−2Reβ)(ΔR−2R)2}

We plot these weighting factors as a function of the particle mass and orbital eccentricity in Figure 5. When analyzing the particle distributions, we only plot the number of particles within the parent ring, which we calculate as

 nring(m)=n(m)w(m). (41)

### 4.3. Reduced collisional probabilities of the largest particles

The very last grid point in the domain of solution will only reduce its number density, as it cannot gain from larger masses either as an fragment or from being part of a redistribution. The full phase space removal, introduced in §3.2, causes its evolution time to become very small compared to all others and leads to a numerical instability. In order to avoid this, we multiply the collisional rates with a weight that smooths to zero for the largest particles

 (42)

for both the projectile and target particle. We chose to be a number a few orders of magnitude larger than and use an arbitrary . The modified collisional rates, therefore, read

 P(m,m′) = Vπκn(m)n(m′)w(m)w(m′)× (43) σw(m)σw(m′)(m13+m′13)2.

We discuss the implications of our choice of the weight function and of its parameters in §5.2.

## 5. Results

As we discussed in §2, collisional cascades in debris disks have been studied extensively in the past decades, with many different analytic and numerical solutions to the problem. To demonstrate the similarities and differences between our model and some earlier ones, we show in the following subsection the results of a few comparison tests. The system variables used by our code for these runs are summarized in Table 2.

We compare our numerical model to three previous well known algorithms, the particle-in-a-box code of Thébault et al. (2003), the dynamical code ACE (Krivov et al., 2000, 2005; Löhne et al., 2008; Müller et al., 2010), and the 1D steady-state solver code of Wyatt et al. (2011). Although we do make an effort to model their systems as accurately as possible, a true benchmark between the codes is impossible. This is due in part to the fact that all models have somewhat different collisional and dynamical prescriptions.

### 5.1. Comparison to Thébault et al. (2003)

A relatively straightforward comparison can be made between CODE-M and the Thébault et al. (2003) model. Although the initial Thébault et al. (2003) model has been subsequently improved in Thébault & Augereau (2007), we chose to compare our results with the former, as they are both particle-in-a-box approaches to the collisional cascade, with some dynamical effects included in a simplified manner.

Thébault et al. (2003) aimed to model the inner 10 AU region of the Pictoris disk, with their reference model being a dense debris ring at 5 AU, with a width of 1 AU and height of 0.5 AU. We adopted their largest particle size of 54 km for the comparison run. However, we adopted a smaller minimum particle size than they did (in our case well below the blowout size), to be able to follow the removed particles more completely. This choice does not affect the actual distribution within the ring. We also had Poynting-Robertson drag turned off. Although both of our models include erosive (cratering) collisions, the Thébault et al. (2003) prescription uses hardness constants () of much softer material than that of our nominal case and a linear relationship between cratered mass and impact energy (the prescription for erosive collisions has been changed in their later paper Thébault & Augereau 2007). For a better agreement, we also model a modified cratered prescription case, where we set , and . With these adjustments our cratering prescriptions agree better; however our interpolation formula is offset compared to Thébault et al. (2003). While ours has a continuous prescription at the CC/EC boundary (i.e., the cratered mass is ), the Thébault et al. (2003) interpolation does not (i.e., the cratered mass is ). Thébault & Augereau (2007) improve on this, by employing an interpolation formula very similar to our equation (23).

Figure 6 compares the evolution of the distribution of particles between the Thébault et al. (2003) nominal case and our runs. In the vertical axes we plot , which is similar to the “mass/bin” value used by Thébault et al. (2003). To make them exact, we divide the Thébault et al. (2003) values by (-1), which places them on the same scale. A few similarities and a few significant differences can be noted. Generally both models show wavy structure - which is a well studied phenomenon (see e.g., Campo Bagatin et al., 1994; Wyatt et al., 2011) - but the exact structure of the waves differs.

Our modified erosive (cratering) prescription model gives a much better agreement with the Thébault et al. (2003) results than our fiducial prescription, in the sense that it yields a deeper first wave in the distribution with a larger wavelength. The offset between the locations of the first dip and the subsequent peak in the two models could likely result from the higher collisional velocities that Thébault et al. (2003) calculate for the smallest particles. Our modified erosion prescription gives good agreement with the Thébault et al. (2003) results for particles larger than a km in size, which is a surprise as the Thébault et al. (2003) erosive constants are for much softer materials than our nominal values. Just above the blow-out regime our model becomes abundant in dust particles, as more and more dust is placed on highly eccentric orbits. Although some smoothing is expected in reality, we do expect the number of dust particles near the blowout limit to increase.

While both our distributions show the typical double power-law feature of quasi steady-state collisional cascades (see e.g., Wyatt et al., 2011) above and below the change in the strength curve, it is masked in the Thébault et al. (2003) model, due to the high amplitude wavy structures. These structures are substantially reduced for the innermost ring of the disk modeled in Thébault & Augereau (2007), but persist in the outer zones.

In Figure 7, we show the differences in the evolution of the disk mass between the nominal case of Thébault et al. (2003) and our models. The left panel shows the evolution of the total disk mass within the debris ring, while the right panel shows the evolution of the dust-to-planetesimal mass ratio. These figures are equivalent to Thébault et al. (2003) figures 2 and 3 (except that these are in plotted in logarithmic scales). Our nominal model predicts a faster decay of the total disk mass, reaching 25% mass loss, while the modified erosive prescription agrees with the Thébault et al. (2003) model and loses 12% of its initial mass. The evolution of the dust-to-planetesimal disk mass differs while the quasi steady-state is being reached, after which all models decay with the same slope. Our nominal case model has an order of magnitude larger dust-to-planetesimal mass ratio at all times compared to the Thébault et al. (2003) model, while our modified erosive collision prescription case is close to it.

There are some easily identified differences between our models. Thébault et al. (2003) use the same Benz & Asphaug (1999) dispersive strength curve as we do, although they do average it to account for impact angle variations. This is an unnecessary step, as the Benz & Asphaug (1999) strength curve is already impact angle averaged, and is corrected in Thébault & Augereau (2007). However, we find that this scaling offset does not have a significant effect on the outcome of the distribution evolution. Thébault et al. (2003) use a double power-law for fragment redistribution, while we use only a single power-law. We find that varying the slope of the single power-law does not have a significant effect on the evolution of the distribution either, so it appears that this difference is also likely not a significant contributing factor to the discrepancies. A noteworthy difference between our models is that Thébault et al. (2003) calculate fragment re-accumulation, while we do not. This is a possible explanation for our discrepancies at high masses, and the offsets we have in the total mass decay.

The most significant difference between the models is that ours uses a single interaction velocity, while Thébault et al. (2003) model the interaction velocity between the critical elliptical orbit smallest particles and the parent ring. This is likely to account for some of the additional offsets for the smallest particles, as higher interaction velocities have been shown to initiate higher amplitude waves (Campo Bagatin et al., 1994; Wyatt et al., 2011). Thébault et al. (2003) also take into account the constant presence of particles smaller than the blowout limit within the parent ring, which we do not. Within denser disks this step may smooth out any features near the blowout limit.

### 5.2. Comparison to Löhne et al. (2008) and Wyatt et al. (2011)

As introduced in §2, the numerical code ACE (Krivov et al., 2000, 2005, 2006, 2008; Löhne et al., 2008; Müller et al., 2010) solves the dynamical evolution of the collisional system as well as the collisional fragmentation, thus a straight comparison to CODE-M cannot be performed. We use their ii-0.3 model (Löhne et al., 2008) for comparison, which is for a relatively wide (7.5-15 AU), extremely dense (1 M total mass with a largest planetesimal size of 74 km) debris ring. We turned off the effects of Poynting-Robertson drag for this comparison model. The initial parameters we assumed are summarized in Table 2. This same system was modeled by Wyatt et al. (2011), whose results we also use for comparison.

In Figure 8, we show the evolution of the dust distribution of the system given by CODE-M, ACE, and Wyatt et al. (2011). As the version of ACE used in Löhne et al. (2008) only modeled catastrophic collisions and the Wyatt et al. (2011) model uses catastrophic collision rates, we also include a CODE-M run in the plot that only models the outcomes of catastrophic collisions. Since the Löhne et al. (2008) values are already downscaled by (), no additional scaling was required of their data. The Wyatt et al. (2011) data points are divided by to convert them to differential number densities. Qualitatively, all distributions agree much better than in the Thébault et al. (2003) comparison (Figure 6); however, there is some scaling offset between the full CODE-M and the other models, especially at large masses.

The wavelengths of the waves roughly agree between CODE-M and ACE, with the single difference being the absence of the strong offset of the first crest in our model; the agreement is also good between CODE-M and Wyatt et al. (2011). The double power-law distribution due to the change from strength to gravity dominated thresholds in the strength curve (Benz & Asphaug, 1999) can be distinguished in all three models, with roughly the same slopes. The ACE and the Wyatt et al. (2011) models maintain their initial number density distribution slope, while CODE-M becomes somewhat steeper for the smallest particles.222In Figure 8 we plot in the y-axis the product , so that a steeper number density slope will show up as a flatter distribution. Our catastrophic-collision-only model has a smaller amplitude and wavelength wave structure than our full model or the other models. The most significant difference between the models is the scaling offset of the full CODE-M model, which we analyze below.

In Figure 9, we show a comparison to figures 1 and 2 of Löhne et al. (2008). In the left panel we show the evolution of (), which is introduced in Equation (11) of Löhne et al. (2008) as

 CL=−˙MdiskM2disk, (44)

This quantity is inversely proportional to the characteristic timescale of the system. As expected, since our system evolves faster, its characteristic timescale is shorter, so the factor for our models is larger. This can be seen in the right panel as well, where we plot the decay of the total mass in our system and that given in figure 2 of Löhne et al. (2008). We adopted the exact strength curve of Löhne et al. (2008) in this run of our model, with the corrections given by Wyatt et al. (2011).

In Figure 10, we show a comparison to the top panel of figure 4 in Löhne et al. (2008), which shows the evolution of the total mass within each of their mass bins. In this plot we show the evolution of the full collisional system, which includes erosive and catastrophic collisions. Since we do not use mass bins, but rather a differential number density griding, we integrate our distribution between 14 grid points for each mass value, which roughly corresponds to a single mass bin of Löhne et al. (2008). Up to roughly a few hundred meters in size (where the strength curve has its minimum) all mass “bins” decay in close parallel slopes to each other after reaching their quasi steady-state around 10,000 yr. This is in contrast to the Löhne et al. (2008) results and agrees more with figure 2 of Wyatt et al. (2011), who model the same system. Our intermediate size planetesimals ( km) show a steeper decay than that modeled by either Löhne et al. (2008) or Wyatt et al. (2011).

The obvious difference between ACE and CODE-M, is that ACE also evolves the dynamics in the system and takes into account the varying collisional velocities in the system from particles that are in elliptical orbits within the parent ring. This could easily explain the offset of the first wave given by ACE in Figure 8.

The increasing offset between the full CODE-M run and the other two models is likely due to the omission of erosive collisions by Löhne et al. (2008) and to using catastrophic collision rates in Wyatt et al. (2011). Kobayashi & Tanaka (2010) have shown earlier erosive (cratering) collisions to be the dominant effect for mass loss in collisionally evolving systems. This effect is demonstrated by the CODE-M model we run with only catastrophic collisions included, which scales exactly with the ACE and Wyatt et al. (2011) models. Since CODE-M does not include aggregation, the collisions of the smallest particles with the largest bodies is not modeled perfectly. We assume the realistic distribution decay to lie between the two models given by CODE-M.

As introduced in §4.3, we artificially reduce the collision cross section of the largest particles in the system to zero in order to avoid numerical instabilities. However, as this is a completely arbitrarily defined numerical necessity, we investigate its effects on the total mass decay, where we expect it to be the strongest. We reproduce Figures 8 & 12, but with varying the values of and in the smoothing function, in Figures 11 & 12. As can be seen in these plots, the variable does not affect either the evolution of the distribution or the total mass decay, as long as it is larger than one. Varying the values of does have an effect on both the evolution of the distribution and the total mass decay. The effect is only on the largest bodies in the system; below a size of one hundred meters, the shapes of the distributions remain unchanged, with the only differences being scaling offsets.

Visual examination of the distributions in Figure 9 hint at a slightly steeper distribution slope for the CODE-M models than the Löhne et al. (2008) and Wyatt et al. (2011) ones. Since the distributions have wavy structures in them, this is difficult to show with a slope fit. Therefore, we calculate the dust-to-planetesimal mass ratios of the distributions. We define the dust sizes to be from 0.1 mm to 10 cm and the planetesimal sizes to be from 10 cm to 100 m. In Figure 13 we plot the evolution of the mass ratios of our models and the mass ratios for the Löhne et al. (2008) and Wyatt et al. (2011) models at 0.5 and 50 Myr. Our models have significantly higher mass ratios than theirs. This is likely a result of the differences between our collisional equations.

## 6. Conclusions

In this paper we present a numerical model of the evolution of the distribution of dust in dense debris disks. We calculate our model with a new numerical code, CODE-M, which we extensively verify and test in the Appendix.

Our collisional model and numerical solution both present improvements to certain aspects of previous numerical collisional cascade models, but also have some limitations. Like most debris disk models (Thébault & Augereau, 2007; Krivov et al., 2006; Löhne et al., 2008; Müller et al., 2010; Wyatt et al., 2011), we solve the Smoluchowski collisional equation (Smoluchowski, 1916), which does not enable the study of stochastic impact events. Our model also does not evolve the dynamical state of the systems; the only code to currently do so is ACE (Krivov et al., 2005; Löhne et al., 2008; Müller et al., 2010). In its current state, our model is also limited to solving the collisional cascade in debris rings, similarly to the initial versions of most collisional models (Krivov et al., 2000; Thébault et al., 2003; Wyatt et al., 2011), some of which have since been expanded to model extended disk structures (Krivov et al., 2005; Thébault & Augereau, 2007). Although our model only uses a single interaction velocity, as we show in §4.1, this is a valid assumption as long as the interactions occur within the debris ring.

We introduce a new prescription for describing erosive collisions, which always takes into account the reduction of the mass in the largest fragment during the cratering event. We introduce a number of approximate interpolations to ensure that our description of erosive and catastrophic collisions yields a continuous set of outcomes as a function of the colliding masses, while being consistent with experimental results at various limits. Moreover, we employ an efficient numerical algorithm that allows us to evaluate the scattering integrals with high precision, considering the enormous dynamic range of masses involved in the calculation.

We compare our code to the previously published numerical models of collisional cascades in debris disks, reassuringly showing general agreement, particularly with Löhne et al. (2008) and Wyatt et al. (2011). Nonetheless, our model shows faster decays than previously published ones (Thébault et al., 2003; Löhne et al., 2008; Wyatt et al., 2011) and also yields slightly higher dust-to-planetesimal mass ratios. We attribute these characteristics to be a result of our accurate treatment of collisional cascades. Our model will be used in a series of upcoming papers to study those aspects of debris disk behavior for which it is uniquely well suited.

We thank Richard Greenberg and David O’Brien for helpful discussions on collisional outcomes and William Hartmann for advice on the smallest particles produced in collisions. We thank Philippe Thébault, Alexander Krivov and Mark Wyatt for their comments on an earlier version of the manuscript. Support for this work was provided by NASA through Contract Number 1255094 issued by JPL/Caltech.

## Appendix A Strength curves

The redistribution outcome of collisions depends almost solely on the energy of the impact and the colliding masses. In experiments it is common to specify the ratio of the kinetic energy of the projectile to the mass of the target. This ratio is known as the specific energy of the impact. Gault & Wedekind (1969) already noticed that the fragment distribution of particles depends on (which they called “rupture energy”) when firing aluminum projectiles into glass targets. Their experiments showed that the fragments will have a power-law distribution, with the largest fragment being a function of the specific energy of the impact. This relationship was first given in equation format in Fujiwara et al. (1977) for basalt targets. They note an offset from the Gault & Wedekind (1969) results, likely due to material strength differences.

Two specific values of are used: (the shattering specific energy) and a somewhat larger (the dispersion specific energy). The value of gives the energy required to shatter the target so that the mass of the largest fragment is no more than half of the original target mass. However, if the target is large enough, then self gravity pulls the fragments back together, leaving a remnant larger than half of the original. The larger gives the value of needed to disperse the fragments, so that the largest remaining piece is half of the original target mass. At lower target masses, where self-gravity can be neglected, . We use in our code and refer to it as .

Determining the value of is difficult, especially for such a large range of particle sizes, from m to km. The values for smaller bodies on the order of a few kilograms are mostly determined from laboratory experiments, while the values for larger bodies are determined from under-surface explosions, observations of large asteroids and with experiments done under very high pressure (Holsapple et al., 2002). However, material strength varies greatly as a function of material type, object size, surface type and the number of shattering events an object has gone through over its lifetime. An object that has gone through many collisions in its lifetime, but still remains in one piece (descriptively called a “rubble pile”) can endure harder collisions, which can actually be absorbed and help to compact the object, rather than dispersing it into smaller particles. This may seem like an important parameter only for larger objects; however, the evolution of larger objects significantly influences the evolution of smaller particles, and thus is important in our study. We also lack experiments done with targets and impactors cooled down to space temperatures of 100-150 K, where one would assume that objects get more brittle and easier to shatter.

Experiments clearly show that is a function of the target mass , meaning that different mass targets will get shattered (with a largest fragment) by different specific energies. Holsapple et al. (2002) reviews experimental and theoretical results on collisions and strength curves. A common result for all of them is a minimum in the strength curve for bodies around 0.3 km in radius, where planetesimals are easiest to disperse (the number of cavities and cracks weakening the bodies increases, while self-gravitation is not dominant yet). As a result, there is a bump in the size distribution of minor planets in the solar system around this size. Smooth particle hydrodynamic (SPH) models give the strength curve for larger bodies, while experiments help to anchor the curve for smaller rocks on the scale of a few cm in radius. It is still not clear whether the power-law shape of the curve can be extrapolated down to micron size particles, where experiments cannot be carried out. To study the collisional evolution of the smallest particles, the exact value of the strength curve must be known. In the absence of any models/experiments currently at those sizes, the best that can be done is a simple extrapolation of the strength curve to those regimes. Stewart & Leinhardt (2009) introduce a velocity-dependent tensile strength curve, that is defined by variables such that it removes ambiguities over material density and projectile-to-target mass ratio. Their tensile strength curve is ideal for low-velocity (1-300 m s) collisions, such as those found during planet formation or at large radii debris disks. However, their universal relationship does not hold for conditions that strongly depart from the catastrophic disruption regime ().

In our models we use the Benz & Asphaug (1999) dispersion strength curve. It is derived from SPH models, represents a reasonable average of all previous strength curves, and is impact angle averaged. This curve can be written as (all units are in SI)

 Q∗(a)=10−4Qsc[S(a1 cm)s+Gρ(a1 cm)g]J gerg−1 kg−1, (A1)

where the fiducial values in the equation are given in Table 1.

## Appendix B Mass conservation of the model

A crucial test of any collisional code is for it to conserve the initial total mass of the system. Since particles are removed at the low mass end, this behavior can be complicated to verify. However, a system can only maintain its total mass numerically, if its collisional equations are mass conserving analytically. Here, we prove that our collisional equation is mass conserving.

The collisional equation can be written as

 dn(m)dt = −∫∞0dm′n(m)n(m′)σ(m,m′) (B1) +∫∞0dμ∫∞μdMn(μ)n(M)σ(μ,M)δ[X(μ,M)−m] +∫∞0dμ∫∞μdMn(μ)n(M)σ(μ,M)R(m;μ,M),

where is the redistribution function to mass from collisions, such that

 ∫∞0dmR(m;μ,M)m=μ+M−X(μ,M), (B2)

and is the Kronecker function. Multiplying Equation (B1) by and integrating over d gives

 dMdt=dn(m)mdt = −∫∞0dm∫∞0dm′n(m)n(m′)σ(m,m′)m (B3) +∫∞0dμ∫∞μdMn(μ)n(M)σ(μ,M)∫∞0dmδ[X(μ,M)−m]m +∫∞0dμ∫∞μdMn(μ)n(M)σ(μ,M)∫∞0dmR(m;μ,M)m,

where

 ∫∞0dmδ[X(μ,M)−m]m=X(μ,M), (B4)

resulting in

 dMdt=−∫∞0dm∫∞0dm′n(m)n(m′)σ(m,m′)m+∫∞0dμ∫∞μdMn(μ)n(M)σ(μ,M)(μ+M) (B5)

The first integral can be separated into two sections as

 ∫∞0dm∫∞0dm′n(m)n(m′)σ(m,m′)m = ∫∞0dm∫m0dm′n(m)n(m′)σ(m,m′)m (B6) +∫∞0dm∫∞mdm′n(m)n(m′)σ(m,m′)m.

Since is a symmetric function, we can swap the limits of integration for and in the second integral of Equation (B6) and, after making a change of variables of and in the first and and in the second integral, the full equation becomes

 dMdt=∫∞0dμ∫∞μdMn(μ)n(M)[σ(μ,M)μ+σ(M,μ)M−σ(μ,M)(μ+M)]. (B7)

Since the collisional cross section is completely symmetric, the integral itself becomes zero, thus proving that our equation is mass conserving.

## Appendix C Numerical evaluation of the model and verification tests

The integro-differential equation presented in §3 must be integrated over 40 orders of magnitude in mass space, contains a double integral whose errors can easily increase if not evaluated carefully, and is bundled in a differential equation that evolves the number densities of dust grains and boulders within the same step. These characteristics demand attention in its numerical evaluation. In the following subsections we explain the numerical methods used to evaluate each integral and the ordinary differential equation (ODE). We also present verification and convergence tests for our code, which explain why such precisions are really necessary.

### c.1. Taylor series expansion of TI

First, we expand equation (13) to use a Taylor series when and . For this, we rewrite in terms of and as

 M=m+m′G(m,m′), (C1)

where equals the cratered mass, is the largest particle created and can be found by root finding algorithms. As written, can be related to the parameter used by Dohnanyi (1969), for which he used a constant value of 130 for collisions. We plot the value of as a function of and in Figure 14, with the thicker solid line giving the contour of . This contour lies at sizes reasonable for experiments in laboratory conditions, which is why Dohnanyi (1969) used a value close to it. The positions of the contours are a strong function of the interaction velocities. The integrand can be written as

 I(m,m′) = f(m′)w(m′)σw(m′)m′−ηa(m)2(f(m)w(m)σw(m)(1+Z)2 (C2)

where and and are dimensionless number densities that can be expressed as

 f(m)=n(m)Cm−η. (C3)

We rewrite this integrand as

 I(m,m′) = f(m′)w(m′)σw(m′)m′−ηa(m)2f(m)w(m)σw(m)((1+Z)2 (C4) −f(M)w(M)σw(M)f(m)w(m)σw(m)[1+G(m,m′)Z3]−η{Z+[1+G(m,m′)Z3]13}2).

The Taylor series for the components are

 (1+Z)2 = 1+2Z+Z2 (C5) ≡ T1

and

 [1+G(m,m′)Z3]−η[Z+3√1+G(m,m′)Z3]2 = 1+2Z+Z2+[2G(m,m′)3−ηG(m,m′)]Z3 (C7) +[2G(m,m′)3−2ηG(m,m′)]Z4−ηG(m,m′)Z5 ≡ T1+T2

Both and are close to 1, while deviates from 1 as approaches . In those cases, the ratio can be expressed as

 σw(M)σw(m)=1+∂σw(m)∂M∣∣∣M=m(M−m)=1−PExp(−mmax−mΘ)Θ[1−Exp(−mmax−mΘ)]m′G(m,m′), (C8)

since we know that . We write this ratio as

 σw(M)σw(m)=1−F. (C9)

The integrand then takes the form

 I(m,m′)=f(m′)w(m′)σw(m′)m′−ηa(m)2f(m)w(m)σw(m)×[T1−(T1+T2)w(M)f(M)w(m)f(m)(1−F)]. (C10)

Rearranging it gives us

 I(m,m′)=f(m′)w(m′)σw(m′)m′−ηa(m)2f(m)w(m)σw(m)×{T1[1−(1−F)w(M)f(M)w(m)f(m)]−T2w(M)f(M)w(m)f(m)(1−F)} (C11)

When

 1−w(M)f(M)w(m)f(m)<10−9, (C12)

we use the approximate formula

 I(m,m′)=f(m′)w(m′)σw(m′)m′−ηa(m)2f(M)w(M)σw(m)×[T1F−T2(1−F)]. (C13)

We use the Taylor series of the components to write the integrand below the limit of (i.e., ). This means that our full integral for the first term () takes the final form

 dfI(m,t)dt=−VπC⎧⎪ ⎪⎨⎪ ⎪⎩μX(m)∫mmindm′ I+mmax∫μX(m)dm′f(m′,t)(m′)−ηf(m,t)(a(m)+a(m′))2⎫⎪ ⎪⎬⎪ ⎪⎭ (C14)

where

In Figure 15, we show the mass of the fragments created when particles of mass and collide. The regions are well defined in our single collisional velocity case. When using collisional velocity that depends on particle size, more than one boundary may exist.

### c.2. Verification of the numerical precision of Ti

To verify the precision of our integrator we set up an equation that is similar to in behavior and that has an analytic solution and compare values given by our code to it. The integral we evaluate both analytically and numerically with our code is

 TVI(m) = ∫μX(m)mmindm′(m′)−η{(m13+m′13)2−[(m+m′Γ)13+m′13]2(m+m′Γm)−η}+ (C15) ∫mmaxμX(m)dm′(m′)−η(m1