SMACK: A NEW ALGORITHM FOR MODELING COLLISIONS AND DYNAMICS OF PLANETESIMALS IN DEBRIS DISKS

Our full numerical model uses an N-body integrator, REBOUND, to solve the equations of motion of the planetesimals and detect collisions, combined with a collision resolution algorithm, SMACK, to calculate the effects of collisions on the velocities and size distributions of the planetesimals. REBOUND (Rein & Liu 2012) is a multi-purpose code originally designed for studying collisional dynamics in planetary rings, freely available under an open-source license from https://github.com/hannorein/rebound. While the Grigorieva et al. (2007) model and LIPAD both detect collisions using a 2D grid, REBOUND detects collisions in 3D, without using a grid; it contains a Barnes–Hut tree module to calculate self-gravity and detect collisions, parallelized with MPI and OpenMP.

In the standard version of REBOUND, collisions are resolved using an instantaneous collision model with a normal coefficient of restitution. This approximation handles a variety of problems in planetary ring dynamics. Debris disks are in a very different physical regime than planetary rings; collision speeds in debris disks are much higher than in rings (km s⁻¹ versus mm s⁻¹) and optical depths are lower (∼10⁻⁴ versus ∼1). So rather than allowing each REBOUND code body to represent one planetesimal, we use each body as either a planet with mass or as a massless "superparticle" representing a group of planetesimals on similar trajectories.

Superparticle methods have already been used extensively to model planetesimal formation within gas disks (e.g., Michikoshi et al. 2009; Rein et al. 2010; Zsom & Dullemond 2008; Johansen et al. 2012; Charnoz & Taillifet 2012). They have also been applied to model planetesimals in the solar system and in debris disks. Charnoz & Morbidelli (2003) used a superparticle approach to study how the dynamics of particles ejected from the Jupiter-Saturn region affected their size distributions, although they did not include the feedback from the collisions on the particle dynamics. Grigorieva et al. (2007) used cylindrical superparticles 5 AU in diameter, fixed to the disk midplane, to model collisional avalanches in debris disks over spans of ∼40 orbital periods. The LIPAD model (Levison et al. 2012) uses a superparticle approach; these authors refer to their superparticles as "tracers."

In SMACK, the superparticles represent collections of planetesimals with a range of sizes, as in Charnoz & Morbidelli (2003). However, in Charnoz & Morbidelli (2003), the evolution of the size distribution of each superparticle is calculated only after the entire dynamical integration of the superparticles is complete. In contrast, SMACK calculates the size evolution at each timestep of the N-body integration, allowing us to model the feedback between the collisions and the dynamics of the superparticles. Each superparticle in SMACK is characterized by an incremental size distribution with logarithmic size bins.

In general, two colliding, fragmenting bodies produce a spray of daughter particles with different trajectories. Simulating this distribution of trajectories in great detail would require SMACK to create new superparticles with each collision, quickly increasing the number of bodies tracked by REBOUND with every collision. Instead, SMACK approximates the outcomes of collisions while keeping the number of integrator bodies constant. When REBOUND detects an encounter between two superparticles, it passes the velocities and size distributions of the overlapping superparticles to SMACK. SMACK returns two superparticles with different velocities and size distributions and REBOUND continues its dynamical integrations using these modified superparticles.

The essence of SMACK is simple. In SMACK, the fragments produced by collisions are swapped between superparticles, so that the new size distributions n_a and n_b are given by

$$\begin{equation} n_a(i) = n_A(i) - P_A(i) + F_B(i) \end{equation} \tag{ 1 }$$

$$\begin{equation} n_b(i) = n_B(i) - P_B(i) + F_A(i), \end{equation} \tag{ 2 }$$

where P_A(i) is the number of parent body particles in size bin i in superparticle A that are lost due to collisions, F_B(i) is the number of daughter particles in size bin i produced by colliding particles in superparticle B, and n_A and n_B are the size distributions of the parent superparticles. The detailed forms of P(i) and F(i) used in SMACK are given in Section 2.3.

This swapping of fragments is a first-order approximation of the velocity distribution of fragments in a planetesimal collision. In a collision between two real planetesimals, the fragments are produced in roughly the center-of-mass frame. But an encounter between two superparticles is more complicated; the center of mass of the superparticles is not the same, in general, as the center of mass of any pair of planetesimals represented by the superparticles. If the mass ratio of the parent bodies is higher (lower) than that of the superparticles, the fragments should be launched in a direction skewed toward that of the more (less) massive parent bodies. The swapping described above crudely approximates this physics.

After swapping the daughter planetesimals, SMACK calculates the kinetic energy lost in the collisions and corrects the superparticle velocities to reflect this loss and to conserve momentum. Let v_A and v_B be the magnitudes of the velocities of the parent superparticles in the center-of-momentum frame and m_A and m_B be the total masses of the parent superparticles. Some fraction of the kinetic energies of the parent superparticles will be lost to planetesimal collisions. This fraction depends on the fraction of planetesimals that experience collisions and the amount of kinetic energy lost by each planetesimal in a collision. Using the energy loss fraction and the collision rates for each pair of size bins, SMACK calculates E_A and E_B, the kinetic energy lost to collisions in superparticles A and B, respectively (as described in Section 2.3). Then, the new superparticle velocities must satisfy the energy conservation law

$$\begin{equation} K = \frac{1}{2} m_a v_a^2 + \frac{1}{2} m_b v_b^2 = \frac{1}{2} m_A v_A^2 + \frac{1}{2} m_B v_B^2 - E_A - E_B, \end{equation} \tag{ 3 }$$

where K is the total kinetic energy of the two superparticles and m_a and m_b are the new total masses of each superparticle, calculated from the new post-encounter size distributions given by Equations (1) and (2). The velocities must also satisfy the momentum conservation law

$$\begin{equation} m_a v_a + m_b v_b = 0. \end{equation} \tag{ 4 }$$

The velocities that solve Equations (3) and (4) are

$$\begin{eqnarray} v_a &= \sqrt{\frac{2 m_b K}{m_a (m_a + m_b)}} \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} v_b &= - \frac{m_a}{m_b} v_a. \end{eqnarray} \tag{ 6 }$$

Equations (1), (2), (5), and (6) give the new size distributions and velocities of the superparticles and define the essential SMACK algorithm.

When a superparticle encounter yields no planetesimal collisions (i.e., all the planetesimals pass through the encounter unaffected), the algorithm gives exact results. When all the planetesimals in each superparticle collide, the algorithm produces an outcome that is a good approximation for the dominant size bins. When there is a mix of collisions and pass-throughs, the algorithm compromises, making errors in the distribution of output velocities, but not in the total energy or angular momentum.

For now, SMACK only models one type of collisional outcome: catastrophic collisions, defined as collisions in which the largest fragment is no larger than half the size of the target. Cratering collisions do not have a significant effect on the steady-state size distribution of a collisional cascade (Dohnanyi 1969), so we neglect them for now. Since we are modeling disks in which the impact velocities are high, we also ignore bouncing, sticking, and gravity between planetesimals. Future versions of SMACK may incorporate these effects.

Consider two superparticles, A and B, that are found to overlap during a given timestep. The number of planetesimals from size bin i in superparticle A that collide with planetesimals in size bin j in superparticle B is

$$\begin{equation} c_A(i,j) = n_A(i) \tau _B(i,j), \end{equation} \tag{ 7 }$$

where n_A(i) is the number of planetesimals in size bin i in superparticle A and τ_B(i, j) is the optical depth along the path of a planetesimal in size bin i passing through superparticle B for collisions with planetesimals of size j. SMACK estimates this optical depth as

$$\begin{equation} \tau _B(i,j) \approx n_B (j) \sigma _{ij} l/V, \end{equation} \tag{ 8 }$$

where σ_ij is the combined collisional cross section of a planetesimal in size bin i and a planetesimal in size bin j, l is the path length traveled by a planetesimal in superparticle A, and V is the volume of superparticle B. The path length l is the distance traveled by superparticle A through the disk relative to the local Keplerian flow since its last encounter. SMACK estimates this distance as l = v_ABt_enc, where v_AB is the magnitude of the relative velocity of superparticles A and B and t_enc is the time since superparticle A's last encounter. The superparticle volume V is the parameter used by REBOUND to determine when a collision occurs. The superparticle volume must be chosen carefully to minimize both computation time (see Section 2.4) and numerical heating (see Section 3.2). We run numerical heating tests before every new simulation to select the optimal superparticle size. The collisional cross section is purely geometric because gravitational effects between planetesimals are ignored. The cross section is given by

$$\begin{equation} \sigma _{ij} = \frac{\pi }{4} (D(i) + D(j))^2, \end{equation} \tag{ 9 }$$

where D(i) and D(j) are the diameters of planetesimals in size bins i and j, respectively.

We aspire to model a system in which each planetesimal involved in a fragmenting collision loses some fraction f_KE of its kinetic energy. This fraction f_KE is a parameter of the code. The total desired energy loss for size bin i in superparticle A during its encounter with superparticle B is calculated by SMACK using

$$\begin{equation} E_A(i) = \displaystyle \sum \limits _{j} f_{{\rm KE}} \frac{1}{2} m(i) v_A(i,j)^2 c_A(i,j), \end{equation} \tag{ 10 }$$

where v_A(i, j) is the magnitude of the velocity of superparticle A in the center-of-momentum frame of a planetesimal in size bin i in A and a planetesimal in size bin j in B.

Some of the colliding planetesimals counted in Equation (7) may shatter, creating a distribution of smaller fragments. SMACK determines which colliding planetesimals will fragment by comparing the collision energy of each pair of colliding planetesimals with the minimum kinetic energy needed for a catastrophic collision. In the center-of-mass frame, the kinetic energy of two colliding planetesimals with masses m(i) and m(j) is

$$\begin{equation} E_{{\rm col}} = \frac{1}{2} \frac{m(i) m(j)}{m(i) + m(j)} v_{{\rm rel}}^2, \end{equation} \tag{ 11 }$$

where v_rel is the magnitude of the relative velocity of the two planetesimals. In laboratory experiments, Hartmann (1980) found that approximately half of this collisional kinetic energy is partitioned into each colliding body, regardless of the mass ratio of the two, so a planetesimal in size bin i will shatter in a collision with a planetesimal in size bin j if

$$\begin{equation} \frac{1}{2} E_{{\rm col}}(i,j,v_{{\rm rel}}) \ge E_{{\rm min}}(i), \end{equation} \tag{ 12 }$$

where E_min(i) is the minimum energy needed to shatter a planetesimal in size bin i.

The minimum shattering energy, E_min(i), is a combination of the gravitational binding energy of the planetesimal and its internal impact strength. We use the minimum energy criteria derived by Durda (1993),

$$\begin{equation} E_{{\rm min}}(i) = \frac{1}{f_{{\rm KE}}} \left(0.822 \frac{G m(i)^2}{D(i)} + \frac{1}{6} \pi S D(i)^3\right), \end{equation} \tag{ 13 }$$

where G is the gravitational constant, m(i) is the mass of planetesimals in size bin i, D(i) is the diameter of planetesimals in size bin i, and S is the impact strength of the planetesimals. We use f_KE = 0.1 (Fujiwara 1982) and a size-independent impact strength of S = 3 × 10⁶ J m⁻³ (Greenberg et al. 1977).

If Equation (12) holds true, the loss of parent bodies in size bin i due to fragmentation in collisions with size bin j is equal to the number of collisions between planetesimals in i and j:

$$\begin{equation} P_A(i,j) =c_A(i,j). \end{equation} \tag{ 14 }$$

If the inequality in Equation (12) does not hold, the planetesimal in size bin i will not shatter and will not be counted as loss, so P_A(i, j) = 0. The planetesimal in size bin j may fragment, however, depending on whether E_min(j) satisfies Equation (12); SMACK performs these calculations by looping through one index at a time. The total loss in size bin i in superparticle A is the sum of the losses due to each size bin in superparticle B:

$$\begin{equation} P_A(i) = \displaystyle \sum \limits _{j} P_A(i,j). \end{equation} \tag{ 15 }$$

Collisions produce fragments in power-law size distributions, which we represent with incremental logarithmic bins. The individual daughter particle distribution resulting from fragmentation in size bin i is

$$\begin{equation} F_A(i,j) = P_A(i,j) \kappa _AD(i)^{\alpha }. \end{equation} \tag{ 16 }$$

The fragment distribution index, α, is a parameter of the algorithm. The constant κ_A is calculated for every collision for each superparticle such that the largest fragment in each distribution is one half the mass of the parent planetesimal that produces the fragments. Again, we find the total fragment gain from size bin i in superparticle A by summing over the size bins in superparticle B that collide with i:

$$\begin{equation} F_A(i) = \displaystyle \sum \limits _{j} F_A(i,j). \end{equation} \tag{ 17 }$$

If the total optical depth for a given size bin is ever greater than one, i.e.,

$$\begin{equation} \tau _B(i) = \displaystyle \sum \limits _{j} \tau _B(i,j) > 1, \end{equation} \tag{ 18 }$$

the loss in that size bin given by Equation (15) could be greater than the number of planetesimals in that bin. We generally try to avoid this situation by choosing the number and volume of the superparticles so that the superparticle encounter time is less than the collision time for the small grains. However, it sometimes occurs anyway as the disk evolves and more small planetesimals are produced. This situation tends to arise for the smallest planetesimals first, since they collide most frequently, and becomes more common with increasing superparticle encounter times and with increasing optical depths.

To address this issue, we adopted a variable-timestep method. If the maximum optical depth for all the size bins in a superparticle is ever found to be greater than 1, SMACK divides the path length traveled by the superparticle since its last encounter into segments such that the maximum optical depth in any size bin is ⩽1. After each segment, the optical depths are recalculated and the energy and planetesimal losses are calculated as described above. The relative velocity of the superparticles is updated after each segment to accurately reflect the energy lost to collisions, but the superparticle trajectories are not changed until the final segment.

This "small particle loop" allows us to run SMACK uninterrupted in regions of unexpectedly high density with fewer superparticles. However, it also introduces noise into the velocity evolution of the superparticles, as the trajectories are not updated during the loop. We therefore check each astrophysical simulation by running it with a few different superparticle sizes to ensure that we get the same result and that the small particle loop is not introducing excess noise.

While REBOUND can be adapted to model small dust grains by adding addition forces to the integrator such as radiation pressure and Poynting–Robertson drag, SMACK cannot simultaneously model dust grains and planetesimals within the same superparticles if they are subject to different forces. We therefore model only grains larger than 1 mm, which do not experience significant radiative forces during their collisional lifetimes, even in a very sparse disk like a zodiacal cloud. The current version of SMACK thus cannot directly simulate disk images at optical wavelengths, but is meant for simulating high-resolution millimeter and sub-millimeter images from instruments like ALMA.

To compare our models with images and photometry of known debris disks, we need to know the face-on optical depth, τ_disk, of the models and how it relates to τ_SP, the total cross section of the planetesimals in the superparticle divided by the cross section of the superparticle, $\pi r_s^2$. Equation (8) implies that the density of an individual superparticle models the local disk density. So, τ_disk is related to the individual optical depth of each superparticle by a linear filling factor, f_SP, where

$$\begin{equation} \tau _{{\rm disk}} = f_{{\rm SP}} \tau _{{\rm SP}}. \end{equation} \tag{ 19 }$$

The filling factor, f_SP, is simply the average number of superparticles that a perpendicular path through the disk would intersect. In other words, f_SP is the inverse of the fraction of the perpendicular path through the disk that would be contained within one superparticle. This factor is approximately

$$\begin{equation} f_{{\rm SP}} \approx \frac{3h}{4r_s}, \end{equation} \tag{ 20 }$$

where h is the full height of the planetesimal distribution and 4r_s/3 is the average length of a chord through a spherical superparticle with radius r_s.

We calculate f_SP before the simulation begins using Equation (20). We generate 10⁶ superparticles with the same orbital element distributions that we will use as the initial conditions. Then, we choose a fiducial perpendicular path through the disk and create a histogram of the positions of the superparticles along that path. We estimate h along that path by normalizing the histogram such that the maximum value is 1 and summing together all the bins. Knowing f_SP allows us to set the total cross section of the planetesimals in the superparticles to yield an initial face-on optical depth of our choosing.

The planetesimal size distributions in the superparticles begin at 1 mm, but dust grains in disks are ground down to ∼1 μm in size before they are removed from the system by radiation pressure. To ensure that we are including the contribution of these smaller dust grains to the cross-sectional area of the planetesimals in the disk, we fit a power law to the size distributions of each superparticle and then extrapolate the size distributions down to 1 μm. We use the extrapolated size distributions when calculating the cross-sectional area of the superparticles during normalization.

Dohnanyi (1969) studied collisional cascades analytically assuming an infinite range of particle masses with each body experiencing a constant impact velocity, assuming a mass-independent material strength. He found that the incremental mass distribution of such a collisional cascade at steady state can be described by a power law with index q = −1.833 (with linear mass bins). The corresponding incremental size distribution with logarithmic bins is a power law with index p = −2.5 (Durda 1993).

Dohnanyi (1969) also found that the index of the equilibrium size distribution is independent of the fragment size distribution index α. Durda (1993) confirmed this in his collisional model and noted that steeper values of α corresponded to a faster convergence of the size distribution to an equilibrium power law. We chose a relatively steep value of α = −2.8 to decrease the computation time required for the size distribution evolution tests.

As an initial test of the collisional evolution simulated by SMACK, we placed 1000 superparticles around a solar mass star with semi-major axes uniformly distributed between 90 AU and 110 AU. The superparticles were given eccentricities uniformly distributed between 0 and e_max = 0.2 and inclinations between 0 and i_max = e_max/2. The other orbital elements were uniformly distributed between 0 and 2π. We distributed the planetesimals in each superparticle among 31 logarithmic size bins ranging from 1 mm in diameter to 1 m using a logarithmic increment of 0.1. In order to keep the average impact velocity constant to match the Dohnanyi (1969) scenario, we turned off the velocity corrections in SMACK by allowing SMACK to update the size distribution of each superparticle after an encounter without updating the superparticle trajectories.

If the assumption of an infinite collisional cascade is relaxed and a cutoff is present in the small-sized end of the size distribution, a wave-like pattern emerges in the size distribution (Campo Bagatin et al. 1994). A real cutoff of this sort occurs in disks at the particle blowout size for the system, producing real wave patterns, but artificial wave patterns can also appear as numerical noise in collisional models with artificial particle size cutoffs. We encountered this wave in our initial tests of the collisional algorithm in SMACK. Our simulated disk had an artificial cutoff in the size distribution at 1 mm.

To remove this artificial wave, we use the method of Durda (1993). For each superparticle encounter, SMACK extrapolates the size distribution for each superparticle to 30 smaller size bins down to 1 μm. The number of collisions due to planetesimals in each of these smaller bins is calculated using Equation (7) and is included in the particle loss and energy loss calculations (Equations (15) and (10)) for each superparticle. This small-particle extrapolation reduces the wave effect of the cutoff to negligible levels.

We tested SMACK's ability to reproduce the convergence to a Dohnanyi (1969) size distribution by running five simulations with different initial indices for the total size distributions for the planetesimals in the disk. For each simulation, we set the vertical optical depth of the disk to 10⁻². We summed the size distribution over all superparticles and fit a power law to the total size distribution at each output timestep. Then, we compared the evolution of the power law indices over time for each of the five simulations. The results are shown in Figure 1. As each collisional cascade evolved, the size distribution index grew or shrank until it reached the Dohnanyi (1969) equilibrium index of −2.5. The size distribution index for each simulation converged to the Dohnanyi (1969) index within 2 × 10⁶ yr. For comparison, the collision time for the smallest planetesimals is ∼10⁵ yr in this simulation.

Zoom In Zoom Out Reset image size

Using a size-dependent breaking strength (e.g., Gáspár et al. 2012b) or a size-dependent velocity distribution (see Pan & Schlichting 2012) can cause breaks in the steady-state distribution of a collisional cascade. We did not attempt to reproduce these phenomena in these initial numerical tests.

Simulations of bouncing collisions suffer from numerical heating associated with the finite size of the superparticles. The superparticles in SMACK do not simply bounce, but the simulations of bouncing collisions can nonetheless inform us about the numerical noise in SMACK. For example, Lithwick & Chiang (2007) found that their grid-based simulations of bouncing collisions contained numerical heating on the scale of their grid size. Likewise, we expect that SMACK cannot model disks where the mean eccentricity of the planetesimals is <r_s/r or the mean inclination of the planetesimals is <r_s/r, where r is the radius of the disk.

We studied the eccentricity evolution of SMACK simulations as a function of superparticle size to search for additional numerical heating effects, like any artificial viscous stirring associated with the finite superparticle size (see Goldreich & Tremaine 1978). We simulationed a planetless ring with radius 90–110 AU and optical depth 5 × 10⁻³ around a solar mass star. We varied the superparticle radius between runs while varying the number of superparticles to keep the superparticle encounter time constant and measured how the mean eccentricity of the planetesimals evolved. We utilized the results of this test to choose the maximum superparticle size to use for the simulation described in Section 4.

At first, we found that simulations using a large range of planetesimal sizes (1 mm–1 m) cooled more slowly and the eccentricity evolution curves never converged for any value of r_s. This situation probably resulted from the coupling between planetesimals of various sizes. To reduce this kind of numerical noise, we decided to restrict the size range of the planetesimals in each superparticle to 1 mm–10 cm for all simulations that included velocity evolution.

Figure 2 shows the simulations with the reduced range of planetesimal sizes. The mean eccentricity of the superparticles in each simulation decreased rapidly at first due to inelastic collisions, then flattened and leveled off at a nonzero eccentricity. As we decreased the superparticle radius, the eccentricity evolution curved converged to a single damping curve, in which the mean eccentricity dropped by a factor of 1/2 in 0.53 Myr or approximately 2τ_col, where τ_col is the collision timescale for the smallest planetesimals. The damping curves converged at a superparticle radius of r_s = 10^−1.3 AU, showing that this superparticle size probably suffices to mitigate numerical heating in our astrophysical simulations.

Zoom In Zoom Out Reset image size

The finite size of the superparticles also leads to diffusion of angular momentum, which can cause a narrow ring to spread on an artificially short timescale (e.g., Lithwick & Chiang 2007). We tested our models for numerical viscous spreading by simulating the evolution of a very narrow (Δr = 0.01 AU) planetless ring at r = 10 AU with superparticle radii of 10^−2.5 AU for 1 Myr. Figure 3 shows the histogram of the number density of superparticles at several times during the simulation. The spreading time of this ring, the time it takes for the width of the ring to double, is τ_spread = 6.67 × 10⁵ yr.

Zoom In Zoom Out Reset image size

Petit & Hénon (1987) solved the diffusion equation for a collisional ring and found that the width of the ring Δr at time t is related to the particle size (in our case, the superparticle size) r_s by

$$\begin{equation} \Delta r(t) =\left(\frac{36 k N r_s^4 t}{t_{{\rm per}}r}\right)^{1/3}, \end{equation} \tag{ 21 }$$

where N is the number of particles, r is the radius of the ring, t_per is the orbital period at that radius, and k is a dimensionless constant. This equation provides a good fit to the simulation shown in Figure 3 when we chose k = 7.6 × 10⁻⁴. With this constant, we can use Equation (21) to choose a superparticle size that yields acceptable levels of numerical viscosity. For example, if we want to limit the spreading of the ring in this simulation to 10% of its initial width in 10⁷ yr, we must choose a superparticle size less than r_s = 10^−4.12 AU.

Collisions between planetesimals produce fragments smaller than the smallest size bin tracked by the superparticles. The mass of these dust particles is effectively lost to the superparticles and carries away some angular momentum. We measured the angular momentum of the superparticles and the mass lost to small fragments in the planetless ring simulation described in Section 3.2 to verify that these torques balanced and that total angular momentum was conserved by SMACK.

We calculated the time derivative of the angular momentum of the superparticles at each output timestep t as

$$\begin{equation} \left. \frac{dL}{dt} \right|_{{\rm superparticle}} \approx \frac{\displaystyle \sum \limits _{i} L_i (t+\Delta t) - \displaystyle \sum \limits _{i} L_i(t)}{\Delta t}, \end{equation} \tag{ 22 }$$

where Δt is the size of an output timestep and L_i(t) is the magnitude of the angular momentum of superparticle i at time t, given by

$$\begin{equation} L_i(t) = m_i(t) |\textbf {v}_i(t)| |\textbf {r}_i(t)| \sin \theta _i(t), \end{equation} \tag{ 23 }$$

where r_i(t) and v_i(t) are the position and velocity vectors (in the heliocentric frame) of superparticle i, m_i(t) is the total mass of superparticle i, and θ_i(t) is the angle between r_i(t) and v_i(t).

We then set SMACK to output the mass lost in every superparticle encounter and the velocity of that mass. We calculated the time derivative of the angular momentum of the small fragments lost as dust in the same way as we calculated the time derivative of the superparticle angular momentum:

$$\begin{equation} \left. \frac{d\mathcal {L}}{dt} \right|_{{\rm dust}} \approx \frac{\displaystyle \sum \limits _{j} \mathcal {L}_j(t+\Delta t) - \displaystyle \sum \limits _{k} \mathcal {L}_k(t)}{\Delta t}, \end{equation} \tag{ 24 }$$

where $\mathcal {L}_j(t)$ is the angular momentum of the small fragments produced in each of the j collisions during the output timestep t and $\mathcal {L}_k(t)$ are the fragments produced in the k collisions during the previous output timestep.

In Figure 4, we plot dL/dt of the superparticles, $d\mathcal {L}/dt$ of the lost mass, and the total for the two populations. The figure shows that most of the angular momentum lost by the superparticles is gained by the lost mass. The total change in angular momentum for the system varies stochastically over the simulation with a maximum variation of 0.48% of the initial total angular momentum, L₀, per Myr. The systematic change in angular momentum is only 1.79 × 10⁻³ of L₀ over 10 Myr. We also calculated the change in the x-, y-, and z-components of the system's angular momentum and found that the systemic change of each component was −3.13 × 10⁻³L₀, −4.46 × 10⁻²L₀, and 1.79 × 10⁻³L₀ over 10 Myr, respectively.

Zoom In Zoom Out Reset image size

In the absence of velocity evolution, the rate of mass loss in a debris disk from planetesimal collisions is given by

$$\begin{equation} \dot{M}_{{\rm disk}}(t) = -CM_{{\rm disk}}^{2} \end{equation} \tag{ 25 }$$

and the total mass in a debris disk at time t is given by

$$\begin{equation} M_{{\rm disk}}(t) = \frac{M_0}{1+t/\tau }, \end{equation} \tag{ 26 }$$

where M₀ is the initial mass of the disk, τ is a characteristic time at which the disk mass is half its initial mass, and C = 1/M₀τ (Dominik & Decin 2003).

We compared the evolution of the ring described in Section 3.2 with this analytic expression. Figure 5 shows the total mass of the disk during each simulation, with the full SMACK model and also for an identical model with velocity evolution turned off (i.e., the velocities of the superparticles were not updated). The mass-loss curves for each simulation are in good agreement with the Dominik & Decin (2003) prediction for t < τ = 3.7 × 10⁵ yr, but the simulations lose mass at a faster rate than predicted for t > τ. The mass-loss curve for the full SMACK simulation begins to flatten at t ≈ 10⁶ yr as the mean eccentricity of the superparticles is damped (see Figure 2), decreasing the rate of fragmentation and slowing the mass loss. Löhne et al. (2008) also studied this problem with a model that included velocity evolution, but they did not report on this effect. More recently, Gáspár et al. (2012a) used a collisional model to demonstrate a similar mass-loss damping to the damping seen in our simulations in Figure 5.

Zoom In Zoom Out Reset image size

Rings of identical initial mass will have different characteristic times depending on their radii and the mean eccentricities of the planetesimals. Wyatt et al. (2007) analytically derived the following power law dependence of C on the radius of the ring:

$$\begin{equation} C \propto r^{-13/3}. \end{equation} \tag{ 27 }$$

This power law is the product of three contributions: the r^−1/2 dependence of the relative velocities in the ring, the r^−5/6 dependence of the total number of projectiles above the minimum required disruptive size on the impact velocities, and an r⁻³ dependence of the density in the ring (Löhne et al. 2008). Löhne et al. (2008) studied mass loss rates in debris disks and replicated this r^−13/3 dependence in their simulations. In our simulations, the initial conditions of each system are determined by a user-input optical depth, so the density of the ring is independent of its size. This choice eliminates the r⁻³ factor for the ring density, so we predict a power-law dependence on ring radius of

$$\begin{equation} C \propto r^{-4/3}. \end{equation} \tag{ 28 }$$

We measured the mass loss in three simulated rings of radial width Δr = 1 AU at five different radii from r = 1 to 3 AU, evolved for 10⁴ yr. We turned off SMACK's velocity evolution for these simulations by allowing SMACK to update the size distributions of the superparticles in each encounter without changing their velocities. We calculated C for each ring and fit a power law to the relationship between C and r for each ring (Figure 6). The resulting dependence of C on radius in our simulated rings was C∝r^−1.24542, in good agreement with the prediction of C∝r^−4/3.

Zoom In Zoom Out Reset image size

Wyatt et al. (2007) also derived the dependence of C on the mean eccentricity, e, of the ring planetesimals:

$$\begin{equation} C \propto e^{5/3}. \end{equation} \tag{ 29 }$$

Löhne et al. (2008) also tested this mass-loss rate dependence on eccentricity, but found a power law dependence of C∝e^9/4. They argued that the analytical model must be incomplete.

We measured the mass loss-eccentricity relationship in our models by simulating three rings at the same radius (r = 1 AU) with different mean planetesimal eccentricities. Again, we turned off the velocity evolution in SMACK to compare the simulations with the predictions of Wyatt et al. (2007). After measuring the mass loss in each ring, we calculated C for each ring and fit a power law to the results (Figure 7). The fit to the simulated data yields an index of 1.56767, in good agreement with the derived index of 5/3. However, at higher mean eccentricities, the radial excursions of the particles on their eccentric orbits can begin to dominate the width of the ring. At this point, increasing eccentricity will decrease the optical depth of the ring, which slows the mass loss rate dependence on eccentricity and flattens out the C versus e relationship. Figure 7 begins to show this effect in our simulations at higher eccentricities.

Zoom In Zoom Out Reset image size