MODELING COLLISIONAL CASCADES IN DEBRIS DISKS: STEEP DUST-SIZE DISTRIBUTIONS

We set up a reference debris disk as a basis for comparison to all other model runs. This model consists of a moderately dense debris disk situated at 25 AU around an A0 spectral-type star, with a width and height of 2.5 AU. This radial distance ensures a moderate evolution speed. The peak emission is in the mid-infrared, with a Rayleigh–Jeans tail in the far-infrared regime, which is the primary imaging window for the Herschel Space Telescope. The total mass in the debris disk is 1 M_⊕, distributed between minimum and maximum particle masses that correspond to radii of 5 nm and 1000 km, when assuming a bulk density of 2.7 g cm⁻³. We summarize the disk parameters of the reference model in Table 1. We evolve the reference model for 10 Gyr.

In Figure 1, we show the evolution of the particle distribution, plotting it at six different points in time up to 10 Gyr. On the vertical axis we plot log₁₀[n(m) × m²], which can be related to the "mass/bin" that is frequently used in other simulations. Even though the number densities decrease with increasing particle masses, the mass distribution increases toward the larger masses in this representation, as long as the mass distribution slope is smaller than 2.

Zoom In Zoom Out Reset image size

The smallest particles reach collisional equilibrium first, roughly at 1 Myr, followed by larger particle sizes as the system evolves. After 50–100 Myr of evolution, the upper, gravity-dominated part of the distribution (m > 10¹³ kg) also reaches equilibrium. The distribution maintains its slope for masses below 10¹⁰ kg, which roughly corresponds to a planetesimal radius of 100 m. The kink in the distribution at the upper end is due to the change in the strength curve slope (O'Brien & Greenberg 2005; Bottke et al. 2005).

Structures in the distribution slope, such as waves, may in principle occur at the low-mass end when assuming softer material properties or higher collision velocities (see Section 3). The distribution may also acquire some curvature (see Paper I). Because of these effects, we evaluate the average slope of the distribution by fitting a power law over a large-mass range, but one that remains below the kink in the distribution. Specifically, we fit the slope between masses 10⁻⁸ and 10⁴ kg, which roughly correspond to sizes of 0.1 mm and 1 m.

We next examine the dependence of the quasi-steady-state distribution slope on the initial conditions, i.e., the initial-mass distribution slope η₀ and the initial total mass in the disk M_tot. Figure 2 shows the evolution of the particle-mass distribution and its slope as a function of these initial conditions. The left panels show the mass distribution after 10 Gyr for different values of the two input parameters, while the right panels show the evolution of the dust-mass distribution slope. Variations in the initial slope, η₀, do not affect the final slope value, although the high-mass end evolves differently or reaches equilibrium at different timescales. A distribution with less dust initially (η₀ < 1.87) also takes more time to reach equilibrium. A shallow distribution with an initial slope of η₀ = 1.5 takes as much as ∼1 Gyr to reach equilibrium, although such initial distribution slopes are unlikely. As shown by Löhne et al. (2008), the evolution of the particle-mass distribution is scalable by the total mass (or number densities of particles), which is what we see in the bottom two panels of Figure 2. All systems with different initial masses reach the same equilibrium mass distribution, but on different timescales. More massive systems evolve on shorter timescales, thus reaching their equilibrium more quickly, while less massive systems evolve more slowly.

Zoom In Zoom Out Reset image size

The parameters that describe the outcomes of collisions, in principle, should be roughly the same for all collisional systems. They are the fragmentation constants and the parameters of the strength curve (Benz & Asphaug 1999). To investigate their effects on the evolution of the particle-mass distribution, we vary their nominal values and evolve the models to the same 10 Gyr, as we did for the reference model.

We give here a detailed analysis of the effects of varying only 5 of the 12 parameters (α, b, Q_sc, s, S), as the remaining seven (γ, β_X, f_Y, f^max_X, f_M, g, G) have no significant effects (see Table 1 for the description of these parameters).

In Figure 3, we plot the resulting mass distributions when varying the cratered mass parameters α and b (Koschny & Grün 2001a, 2001b). The parameter α is the total scaling and b is the exponent of the projectile's kinetic energy in the equation of the cratered mass

$$\begin{equation} M_{\rm cr} = \alpha \left(\frac{\mu V^2}{2}\right)^b\;. \end{equation} \tag{ 1 }$$

As can be seen in Figure 3, the resulting mass distributions depend on the values of α and b only in the gravity-dominated regime. At these larger masses, our model is incomplete, because we do not include aggregation. When increasing α, i.e., basically softening the materials or increasing the effects of erosions, the number of eroded particles in the gravity-dominated regime increases rapidly. A similar effect can be observed when increasing the value of b. However, within reasonable values of α and b, the variation of the equilibrium particle-mass distribution slope in the dust mass regime is negligible.

Zoom In Zoom Out Reset image size

In Figure 4, we plot the resulting distributions and the evolution of the dust distribution slope when we vary the parameters Q_sc, S, and s. These are all variables in the tensile strength curve, which is given as (Benz & Asphaug 1999)

$$\begin{equation} Q^{\ast }(a) = 10^{-4} Q_{\rm sc}\left[S \left(\frac{a}{{\rm cm}}\right)^{s} + G \rho \left(\frac{a}{{\rm cm}}\right)^{g}\right] \frac{{\rm J g}}{{\rm erg}\,{\rm kg}}. \end{equation} \tag{ 2 }$$

The variable Q_sc is a global scaling factor, S is the scaling of the strength-dominated regime, s is the power dependence on particle size of the strength-dominated regime, G is the scaling of the gravity-dominated regime, and g is the power dependence on particle size of the gravity-dominated regime. Variations in the gravity-dominated regime of the curve (G and g) do not have significant effects on the equilibrium dust-mass distribution, so we do not consider these parameters further. The tensile strength curve has been extensively studied for decades. However, as it is dependent on various material properties and the collisional velocity (Stewart & Leinhardt 2009; Leinhardt & Stewart 2012), its parameters do not have universally applicable values. Determining the tensile strength curve at large and small sizes is also extremely difficult experimentally.

Zoom In Zoom Out Reset image size

The slope s of the strength curve in the strength-dominated regime depends on the Weibull flaw-size distribution. Its measured values range anywhere between −0.7 and −0.3 (Holsapple et al. 2002). Steeper values of s make smaller materials harder to disrupt, which results in a steeper dust distribution slope. At s = 1.0, the smallest particles are hard enough to resist catastrophic disruption even when the projectile mass equals the target mass. This results in a mass distribution with a slope equal to the redistribution slope γ = 1.83 at the smallest scales, while in the fitted region it is significantly steeper. At s = 0.6, the smallest particles are still able to destroy each other and generate a dust distribution slope that is close to 1.91.

The scaling constants of the tensile strength curve are the dominant parameters in the evolution toward the quasi-steady-state distribution. When reducing the complete tensile strength curve scaling Q_sc, wave structures form more easily, as a particle becomes capable of affecting the evolution of particles much larger than itself (see Paper I). When upscaling the tensile strength curve, the quasi-steady-state distribution slope starts to resemble the redistribution slope, as it is the particle redistributions that leads the evolution of the particle-mass distribution. When varying the scaling of only the strength side of the curve S, similar effects can be seen.

We have shown that drastic offsets in the collision parameter values result in slight changes to the quasi-steady-state mass distribution slope. We conclude that, for all reasonable values of the collisional parameters, the quasi-steady-state dust-mass distribution slope is larger or equal to 1.88.

There are a number of parameters that can change from one collisional system to another: the material density ρ, the minimum and the maximum particle mass in the system m_min and m_max, the radial distance R, height h, and width ΔR of the disk, and the spectral type of the central star. All these parameters affect three properties of the collisional model: the blowout mass, the collisional velocity, and the number density of particles. Varying these parameters will change the timescale of the evolution and affect the quasi-steady-state distribution slope. In this subsection, we analyze the effects of varying the radial distance on the equilibrium mass distribution through the dependence of the collisional velocity on radius. Modifying either disk parameters ΔR and h or the spectral type of the star would have similar effects. We defer discussion of the variations in the timescales to a future paper.

In the left panel of Figure 5, we show the effects of varying the radial distance, R, on the mass distribution evolved to 10 Gyr. Decreasing the radial distance will increase the collisional velocity, resulting in the appearance of waves at the small-mass end of the distribution. It also generates a much more pronounced kink at the high-mass end. On the other hand, when the velocity is decreased at large radii, the low-mass end of the distribution starts resembling the redistribution function, as smaller particles are not destroyed due to the lower energy collisions. Moreover, no kink is produced at the high-mass end. In the right panel of Figure 5, we show the evolution of the particle-mass distribution slope as a function of the collisional velocity. For high-velocity collisions, the waves render the fitting of a single-mass distribution slope ill constructed but the underlying slope of the wavy-mass distribution is slightly steeper than for the smaller collisional velocity case.

Zoom In Zoom Out Reset image size

We discuss here the effects of three non-physical variables that appear in the numerical algorithm. They are: the mass ratio δ between neighboring grid points and the parameters of the large particle collisional cross-section smoothing formula, Θ and p (see Paper I). In Figure 6, we plot the distribution of the model with varying values of δ at 10 Gyr (left panel) and the evolution in the slope of the dust distribution (right panel). The evolution of the dust distribution is affected by the number of grid points we use, converging at δ = 1.13; this corresponds to 801 grid points between our m_min and m_max mass range. Using a lower number of grid points leads to errors in the numerical integration for the redistribution, leading to an offset larger than 7% for the smallest particles in the system, when only using half as many grid points. We find that the dust distribution slope is practically independent of the smoothing variables Θ and p.

Zoom In Zoom Out Reset image size

In our model calculations, the dust distributions in the vast majority of cases reach quasi-steady state by 10–20 Myr and only in a few cases do they take somewhat longer. The characteristic time is less than 100 Myr for all realistic cases in agreement with previous studies (e.g., Wyatt et al. 2007; Löhne et al. 2008). This shows that, apart from second generation debris disks, the majority of debris disks around stars of ages over 100 Myr are most likely to be in collisional quasi-steady state, at least for the smallest particles (<1 cm). However, young and extended systems, such as β Pic (Smith & Terrile 1984; Vandenbussche et al. 2010), might not be near quasi-steady state at the outer parts of their disks, where interaction timescales are longer.

One of the most surprising results of the wide range of numerical models we computed is the robustness of the quasi-steady-state distribution. Varying the values of the model parameters does not result in significant changes in the slope of the distributions. In Figure 7, we show the effect on the equilibrium slope of the dust-mass distribution function of varying each model parameter. We order the parameters on the horizontal axis as a function of decreasing magnitude of their effects. Parameters that had to be varied by many orders of magnitude to have effects on the equilibrium slope are analyzed in log space. The plot shows that the dominant parameter, by far, is the slope of the strength curve in the strength-dominated regime. This is followed by variables that affect the collisional velocity (ΔR and h) and the power b of the erosive cratered mass formula (Equation (1)). The plot also shows that neither of our arbitrarily chosen collision prescription constants have any significant effects on the outcome of the collisional cascade. The model runs predict an equilibrium dust-mass distribution slope of η = 1.88 ± 0.02 (η_a = 3.65 ± 0.05), taking the maximum offsets originating from the 50% variation in the model parameters test as our error.

Zoom In Zoom Out Reset image size

Waves appear superposed on the mass distribution due to the discontinuity at the low-mass cutoff. This behavior has been explained before as a result of the buildup of particles at the boundary (Campo Bagatin et al. 1994; Wyatt et al. 2011); however, examining the properties of the collisional equation (see Paper I) yields a different picture. In Figure 8, we plot the dμ integrand of the collisional equation

$$\begin{eqnarray} \frac{{d}n(m)}{{d}t} &\propto & \int _{m_{\rm min}}^{m_{\rm max}} {d}\mu \Biggl \lbrace - m^{-\eta } \mu ^{-\eta } \left(m^{\frac{1}{3}} + \mu ^{\frac{1}{3}}\right)^2 \Biggr. \nonumber \\ && + \mu ^{-\eta }M^{-\eta }\left(\mu ^{\frac{1}{3}} + M^{\frac{1}{3}}\right)^2 \delta [X(\mu ,M)-m] \nonumber \\ && + \int _{\mu }^{m_{\rm max}} {d}M \mu ^{-\eta }M^{-\eta }\left(\mu ^{\frac{1}{3}} + M^{\frac{1}{3}}\right)^{2} \nonumber \\ &&\Biggl. \times A(\mu ,M)H[Y(\mu ,M)-m]\Biggr \rbrace \;, \end{eqnarray} \tag{ 3 }$$

which basically gives the amount of variations produced by a μ mass projectile for a certain m mass, either by μ removing m or adding it as an X fragment. In regimes marked with a "+" the effective changes in the number densities are positive, while in regimes marked with a "−" they are negative. With solid black lines we show the integrands for the m masses where the first dip of the waves are situated for the corresponding velocities. It is apparent that for these masses the peaks of the removal integrands are located exactly at the cutoff. The fact that they are peaks is obvious when comparing with the integrands of larger masses. If the distributions were continuous, this effect would smooth out, but with the cutoff located at the peak of the integrand, relatively more will be removed of mass m at the first dip than of other masses. We also note that the additive peaks of the minimum (cutoff) mass, shown with red solid curves, are not correlated with the location of the dip in the distribution (marked with a black solid tick mark), but shift as a function of the collision velocity. At higher collisional velocity, where waves are stronger, the largest number of cutoff particles are actually produced by the particles in the dip, which would provide a negative feedback according to the traditional picture, canceling the waves. Our analysis shows that positive feedback is not the dominant effect in the production of the waves.

Zoom In Zoom Out Reset image size

The location of the peak of the removal term (Term I in Paper I) is the determining factor in the formation of the waves, which can be given by solving

$$\begin{eqnarray} 0 & = &\partial \left\lbrace \mu ^{-\eta }\left[m^{-\eta }\left(\mu ^{\frac{1}{3}}+m^{\frac{1}{3}}\right)^2 \right.\right. \nonumber \\ && \left.\left.\left.- M^{-\eta }\left(\mu ^{\frac{1}{3}}+M^{\frac{1}{3}}\right)^2 \right]\right\rbrace \Biggl /\Biggr.\partial \mu \right|_{\mu =m_{\rm min}}\;, \end{eqnarray} \tag{ 4 }$$

where m gives the location of the first dip and M is given so that X(μ, M) ≡ m. Unfortunately, an analytic solution cannot be given, as the derivative is transcendent, but is easily solvable numerically. The variable that determines the solution will be X, which depends on the collisional velocity, the tensile strength law, and the value of m_min. For a given system, where m_min and the tensile strength law do not change as a function of radial location, the single variable determining the wavy structure is the collisional velocity. The minimum mass will change as a function stellar types and some material/velocity dependence of the tensile strength law can be expected between systems.

In Section 2, we have shown that variations in the proxies for the collisional velocity (such as disk radius) can in fact result in wavy size distributions, even within a narrow debris ring. Extended debris disks will be even more likely to have higher collisional velocities, as the particles with higher β values (but less than 0.5) get placed on high eccentricity orbits. These small particles from the inner rings will collide with the particles in the external rings with increased velocities due to the non-zero collisional angles (Thébault & Augereau 2007). To study the effects of variations only in the collisional velocity, we set the collisional velocity directly within our code to specific values, not changing any other parameters. We present the results from these runs in Figure 9. The figure shows that waves start to appear at collisional velocity values of 3 km s⁻¹ and above.

Zoom In Zoom Out Reset image size

The value of the orbital velocity for a particle on an elliptical orbit originating from R_in with a certain β value at R_out radius (where it will collide with a particle on a circular orbit) is

$$\begin{equation} v_{\rm e} = \sqrt{{\rm G}M_{\ast }(1-\beta)\left[\frac{2}{R_{\rm out}}-\frac{1-2\beta }{(1-\beta)R_{\rm in}}\right]}\;. \end{equation} \tag{ 5 }$$

The angle between the orbital velocities can be expressed as

$$\begin{equation} \zeta = \cos ^{-1}\left[\frac{(2ae)^2-R_{\rm out}^2-(2a-R_{\rm out})^2}{2R_{\rm out}(R_{\rm out} - 2a)}\right], \end{equation} \tag{ 6 }$$

where

$$\begin{eqnarray} a&=&R_{\rm in}\frac{1-\beta }{1-2\beta } \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} e&=&\frac{\beta }{1-\beta }\;. \end{eqnarray} \tag{ 8 }$$

The collisional velocity is then

$$\begin{equation} v_{\rm coll} = \sqrt{\left(v_{\rm o}-v_{\rm e}\right)^2 + \left(v_{\rm e}\sin \zeta \right)^2}\;, \end{equation} \tag{ 9 }$$

where v_o is the orbital velocity for the particle on a circular orbit at R_out.

We plot the values of this collisional velocity for particles originating from R_in = 10 AU and 50 AU in systems around a solar- and early-type star in Figure 10. We shade the collisional velocity region above our 3 km s⁻¹ limit, where waves start to appear. As we can see, nearby rings at ΔR = 1 AU have low collisional velocities, reassuring that our prescription for a constant collisional velocity value in narrow debris rings is reasonable. We can also see that once we depart from the narrow ring assumption (ΔR > 10 AU), the collisional velocities do increase.

Zoom In Zoom Out Reset image size

Around a solar-type star, particles originating from an inner radius of 10 AU are able to achieve collisional velocities of ∼3 km s⁻¹; however, the particle-size range that is able to achieve this limit is very narrow, between 0.4 and 0.6 μm. Taking into account the dilution of these small particles within a limited narrow size range, we conclude that extended debris disks with inner boundaries outside of 10–15 AU around solar-type stars should be well approximated by a simple power-law mass (size) distribution function. In Paper I, we show a comparison run to Löhne et al. (2008) and Wyatt et al. (2011), assuming an extended disk between radii of 7.5 and 15 AU around a solar-type star. The runs by Wyatt et al. (2011) and our code show negligible/small amplitude waves, while Löhne et al. (2008) show moderate amplitude waves, likely due to their full three-dimensional modeling of the system, taking into account smaller particles originating from regions inward of 10 AU. These results confirm our analytic analysis on the conditions required to initiate waves in the dust particle-mass distribution.

We show the same analysis for early-type stars in the right panels of Figure 10. The collisional velocities of particles originating from 10 AU with particles on external circular orbits around an A0 spectral-type star will be significantly larger than our limit for a larger range of particle sizes (5–10 μm). Particles originating from orbits inside of 10 AU will definitely initiate waves in extended disks or external rings. However, particles originating from 50 AU will barely reach collisional velocities at or above our limit and for a limited range of particle sizes. We conclude that extended debris disks with inner boundaries outside of 50–60 AU around early-type stars should be well approximated by a simple power-law mass (size) distribution function.

Our one-dimensional "particle-in-a-box" numerical code adjusts the collisional velocities through the ring thickness and height (ΔR and h), and the resulting orbital inclinations for the particles. Our results are only moderately affected over the relatively narrow range of velocities obtained by varying them (see Figure 7). We modified the code (see above) to simulate the higher-velocity collisions with particles on elliptical orbits to determine a threshold for generation of waves in the particle distribution. However, our model cannot explore the nuances resulting from varying collisional velocities in extended systems (such as varying collisional rates and collisional outcomes), which could also play a role in the possible formation of substructures superposed on the generally steeper particle-size distributions (see, e.g., Krivov et al. 2005; Thébault & Augereau 2007; Löhne et al. 2008; Müller et al. 2010).

As stated at the end of Section 3.1, our generalization assumed a specific system (with same material properties and minimum particle mass). Although we explored varying the minimum particle mass when introducing different spectral-type central stars, both the tensile strength and the minimum mass may vary with material properties as well. In addition to minerals such as basalt, it is likely that grains in the outer zones of debris disks contain significant amounts of ice. The mechanical properties of icy grains are explored by Benz & Asphaug (1999) and Leinhardt & Stewart (2009), yielding roughly a factor of three in difference between the tensile strength of basalt and ice. Given the relative insensitivity of our results to the scaling of the grain strength curve (Q_sc, see Figure 7), these properties are similar enough that our basic conclusions about the grain size distribution and region of applicability of the power-law approximation to it are essentially unchanged even outside of the iceline.