THE COLLISIONAL DIVOT IN THE KUIPER BELT SIZE DISTRIBUTION

The Kuiper Belt is a population of planetesimals outside the orbit of Neptune which exhibits high inclinations and eccentricities by many belt members (Trujillo et al. 2001; Brown 2001), and has a very low mass, M ≲ 0.1 M_⊕ (Gladman et al. 2001; Bernstein et al. 2004; Fuentes & Holman 2008). The high relative encounter velocities, v_rel ∼ 1 km s^-1 (Dell'Oro et al. 2001), and infrequent collisions of the largest members (Davis & Farinella 1997; Durda & Stern 2000) make the growth of Eris-sized bodies impossible over the age of the solar system. Accretion in the early stages of planet building must have been in a more dense and quiescent environment allowing large objects to grow (Stern & Colwell 1997).

The size distribution of the Kuiper Belt is the result of the accretionary and collisional processes that have gone on in that region, and therefore provides one of the main constraints on those processes. For a general review of the Kuiper Belt size distribution, see Petit et al. (2008).³ There are three properties which describe the general shape of the Kuiper Belt size distribution:

• 1.  
  the existence of the largest objects, Eris and Pluto;
• 2.  
  the large object size distribution for D ≳ 100 km which is well characterized by a steep power law $\frac{dN}{dR}\propto r^{-q}$ with q ∼ 4.8 (Gladman et al. 2001; Petit et al. 2006; Fraser et al. 2008);
• 3.  
  the existence of a "roll-over" at r ∼ 25–60 km where the size distribution flattens to a much shallower distribution than for larger objects (Bernstein et al. 2004; Fuentes & Holman 2008; Fraser & Kavelaars 2009).

These three properties need to be reproduced by any successful attempt at recreating the accretionary and collisional history of the Kuiper Belt region. While this has not yet been achieved, these properties provide some insight into the general history of the belt.

Accretion simulations such as those of Stern (1996), Kenyon & Bromley (2001), and Kenyon (2002) have demonstrated that objects as large as Eris could have accreted on timescales shorter than the age of the solar system if the mass of the Kuiper Belt was at least two orders of magnitude larger than the current mass, and if the relative encounter velocities, and hence, the eccentricities and inclinations of the objects were significantly lower in the early solar system (see review by Kenyon et al. 2008). The existence of the steep large object size distribution, however, suggests that accretion was a short-lived process (Kenyon 2002; Fraser et al. 2008). Some event must have disrupted accretion—likely the same event which scattered the primordial Kuiper Belt onto orbits with high inclinations and eccentricities as observed today (see a review of proposed processes by Morbidelli et al. 2008)—otherwise the large object size distribution would be too shallow to be compatible with the observed distribution.

The simulations of Stern (1996), Kenyon & Bromley (2001), and Kenyon (2002) cannot reproduce the large roll-over size. In these simulations, interaction velocities remain low, such that objects larger than r ∼ 1 km do not experience disruptive collisions over the age of the solar system. Thus, the size distribution for objects larger than r ∼ 2 km exhibits a steep slope comparable to that observed today, and the accretion break, or roll-over, occurs at radii too small to be compatible with observations. This suggests that after the belt was dynamically excited, it remained massive enough for a time long enough to allow collisional erosion to produce the large roll-over observed today.

A model presented by Kenyon & Bromley (2004) could reproduce the large roll-over if gravitational stirring by Neptune at its current location was included early in the simulations. Accretion occurred at a sufficient rate in these simulations to produce Eris-sized objects. In these simulations, accretion lasted too long, however, resulting in a large object size distribution too shallow to be compatible with observations. In addition, the mass loss due to collisional grinding was insufficient to produce the tenuous modern-day belt (see discussions by Morbidelli et al. 2008 and Kenyon et al. 2008). The Kuiper Belt must have undergone a more rapid excitation and mass depletion than occurred in these calculations.

Pan & Sari (2005) presented an analytical, order of magnitude, collisional evolution model, reminiscent of the ground-breaking work of Dohnanyi (1969). With this model, they demonstrated that the large roll-over size could be produced by collisional grinding on timescales shorter than the age of the solar system. They, however, assumed that the size distribution rolled over to collisional equilibrium. It has been demonstrated that there is a range of objects with sizes smaller than the roll-over, which are preferentially eroded but which are not being replenished from fragments of larger disrupted objects (Kenyon 2002). These objects do not achieve collisional equilibrium. Equilibrium is only achieved at some smaller radius, the exact size of which depends on the density of planetesimals. Thus, the model of Pan & Sari (2005) likely does not predict the correct shape of the size distribution smaller than the roll-over, or the rate at which the roll-over evolves to larger sizes.

Benavidez & Campo Bagatin (2009) present a collisional evolution model of the Kuiper Belt, in which they include collisions from multiple dynamical classes. Using this model, they calculate the collisional evolution of the Kuiper Belt, starting from an initially steep size distribution for all sizes, with slope similar to that observed for large objects. Their results confirm the findings of Pan & Sari (2005); collisional grinding with relative velocities comparable to that observed can disrupt the majority of objects with r ∼ 50–100 km, and in effect produce a roll-over at that size. These calculations, however, likely do not reproduce the collisional history of the Kuiper Belt, as they start from an initial condition not expected from standard accretion processes (Kenyon 2002). Accretion simulations, which produce Eris-sized objects, produce size distributions with an accretion break at r ∼ 1–2 km, not a steep distribution for all sizes.

Here we present a collisional evolution model, which we utilize to calculate the size distribution of the Kuiper Belt from some starting condition. With this model, we wish to determine whether or not collisional erosion could produce the r ∼ 20–50 km roll-over, starting from a size distribution similar to that produced by models of accretion in the outer solar system, but in the dynamically excited conditions of the current day Kuiper Belt. From these calculations, we wish to constrain the shape of the modern day size distribution smaller than the roll-over, and to constrain the mass of the Kuiper Belt at the end of accretion.

In Section 2, we present our collisional evolution model and the planetesimal strength and shattering models we adopt. In Section 3, we present the main results of our calculation, namely that that existence of a break at r ∼ 1 km left-over from accretion causes a divot in the size distribution at larger sizes. In Section 4, we present the consequences of our model, and finish with concluding remarks in Section 5.

Our collisional evolution model uses a formalism similar to that of Ohtsuki et al. (1990) which uses bins of constant mass, but tracks average bin radius.⁴ The model considers a swarm of planetesimals distributed in a finite number of size bins, where bin i contains a number of objects, N_i, with radii r_i < r < r_i + δr_i, where δ ∼ 1.1.⁵ In our model, we assume that all particles have the same relative collision velocity, v_rel. We do not include velocity evolution in our model, and assume that v_rel is a constant during the simulations. Each collision can either result in catastrophic disruption and dispersal if the collision energy, Q, is sufficiently high, or cratering and mass accretion otherwise. We adopt the standard center-of-mass collision energy $Q=\frac{1}{4} \frac{m_i m_j v_{\rm rel}}{m_i+m_j}$. In the following sections, we describe the details of the model and our prescription for the collisional outcomes.

2.1. Disruption and Cratering Model

We adopt the catastrophic disruption energy threshold functional form presented by Benz & Asphaug (1999). This form includes contributions from the tensile strength of the material body, as well as a gravitational binding term, and is given by

$$\begin{equation} Q_D^*=Q_o\left(\frac{r}{1\,{\rm cm}}\right)^a+B\rho \left(\frac{r}{1\,{\rm cm}}\right)^b, \end{equation} \tag{ 1 }$$

where Q*_D is defined as the disruption energy per unit mass of the target body required to shatter and disperse 50% of the target into a spectrum of fragments, ρ is the target density which we assume is constant for all sizes, and Q_o, a, b, and B are constants appropriate for the material properties of Kuiper Belt objects.

Benz & Asphaug (1999) utilized smooth particle hydrodynamic simulations to determine the outcome of disruptive collisions between icy bodies. They calibrated Equation (1) under a large range of impact parameters (velocity, size, ratio, etc.) for both basalt and water ice targets. They found that the disruption energy and hence the constants Q_o, a, b, and B depended on the relative impact velocity. They also found that the mass of the largest remaining fragment, M_LRF, depended linearly on the ratio of the impact energy to the disruption energy, Q/Q*_D. They demonstrated that the mass of the largest remaining fragment can be well represented by

$$\begin{equation} \frac{M_{\rm LRF}}{M_{t}} = \gamma = Y-X\left(\frac{Q}{M_{t}Q_D^*}\right), \end{equation} \tag{ 2 }$$

where M_t is the initial target mass. Examination of Figure 9 from Benz & Asphaug (1999) reveals that Y = 1.0 and X = 0.55 in Equation (2) is an acceptable representation of γ, and we adopt these parameters for our calculations. Stewart & Leinhardt (2009) confirmed the findings of Benz & Asphaug (1999) and from their simulations determined that X ≈ 0.48. The small difference in X found by Stewart & Leinhardt (2009) and the value we adopt here will make an insignificant difference on our results, and will have no effect on our conclusions.

We assume that the distribution of collisional fragments from a disruption is a power law, and is given by $dN/dr = Ar^{-q_D}$, where q_D is the logarithmic slope. It has been shown that in the asteroid belt, the distribution of fragments is better represented by a two slope distribution (Bottke et al. 2005). In our calculations, we considered a range of two sloped models. We found that, for a reasonable range of two-sloped models, the change in the size distribution results was small compared to the changes caused by adjustment of the strength and cratering laws—also much less than the variation introduced into the results from the unknown collisional history. Thus, we choose the power law as a simple, scale-free, one parameter representation that avoids including much of the uncertain physics of fragmentation in our model.

Assuming that the density of planetesimals is constant for all sizes, from Equation (2) we see that the radius of the largest remaining fragment is given by

$$\begin{equation} r_{\rm LRF} = r_t \gamma ^{1/3}, \end{equation} \tag{ 3 }$$

where r_t is the radius of the target body. Assuming there is only one largest remaining fragment, the normalization of the fragment model is given by $A=(q_D-1) \gamma ^{(q_D-1)/3} r_t^{(q_D-1)}$ where q_D > 1. Conserving mass, we have

$$\begin{equation} \int _0^{r_{\rm LRF}} A r^{-q_D} m(r) dr = m(r_t), \end{equation} \tag{ 4 }$$

where m(r) is the mass of an object with radius r. Equation (4) can be solved for q_D to give

$$\begin{equation} q_D=\frac{\gamma +4}{\gamma +1}. \end{equation} \tag{ 5 }$$

As γ>0, for our fragment model, 1 < q_D < 4.

For impacts with energy insufficient to disrupt the target, we consider a cratering model where the excavated crater mass is given by

$$\begin{equation} M_{\rm crat}=\alpha f_{\rm frac} Q, \end{equation} \tag{ 6 }$$

where α is the crater excavation coefficient, which depends on the material properties of the target and impactor, and f_frac is the fraction of kinetic energy that went into shattering the target material. We consider values of α = 10⁻⁸–10⁻⁹ s² cm^-2 (Petit & Farinella 1993).

We assume that the mass fraction of ejected material with velocity greater than v is given by

$$\begin{equation} f_M(v)=\left(\frac{v}{v_{\rm min}}\right)^{-k} \end{equation} \tag{ 7 }$$

where $k=\frac{9}{4}$ (Petit & Farinella 1993). By assuming that a fraction, f_KE, of the impact energy is imparted into kinetic energy of the excavated material (Equation (6)), and requiring conservation of energy, the minimum ejecta velocity is given by $v_{\rm min}=\sqrt{\frac{2}{9} \frac{f_{\rm KE}}{\alpha f_{\rm frac}}}$, with the fragment mass escaping the body given by $M_{\rm esc}=M_{\rm crat}\big(\frac{v_{{\rm esc}\, i,j}}{v_{\rm min}}\big)^{-k}$, where $v_{{\rm esc}\, i,j}=\sqrt{\frac{2G(m_i+m_j)}{r_i+r_j}}$ is the mutual escape velocity of colliding pair i, j (Wetherill & Stewart 1993). Clearly, accretion will occur if M_esc is smaller than the mass of the impactor.

We assume that the crater fragment distribution is a power law with slope q_c = −3.4, and a largest remaining fragment

$$\begin{equation} r_{{\rm LRF},c}=\left(\frac{4+q_c}{4\pi \rho } M_{\rm crat}\right)^{1/3}. \end{equation} \tag{ 8 }$$

Simulations of Benz & Asphaug (1999), Leinhardt & Stewart (2009), and Stewart & Leinhardt (2009) have demonstrated that catastrophic collisions still occur for energies as little as Q ∼ 0.2 Q*_d. If we consider a maximum crater size as a fraction of the target body radius, f_crat, and require continuity between cratering and disruption regimes, the energy threshold for which collisions is just disruptive, as a fraction of Q*_D, is given by $f_{\rm dis}=\frac{1}{X}(1- (1-f_{\rm crat})^3)$. A minimum disruption energy criterion f_dis = 0.2 corresponds to f_crat = 0.04. In this way, events which strip more than ∼10% of the mass of the target are considered disruptive collisions, while those that strip a smaller amount are considered cratering events. In our calculations, we consider maximum crater sizes f_crat = 0.05–0.15 corresponding to minimum catastrophic disruption energy thresholds of 0.4–0.7Q*_D.

We note here that there are many effects other than size that determine disruption strength and collisional outcome, such as impactor size and target porosity (see discussion by Leinhardt et al. 2008). Given the uncertainty in the knowledge of the processes that shaped the Kuiper Belt size distribution, the specifics of these effects are of little import, and we do not consider them in the collisional model given by Equations (1)–(8).

2.2. Collision Rate

From simple cross-section arguments, it can be shown that, if objects from bin i are passing through the swarm of objects from bin j, the number of collisions between the objects from bins i and j, in time step Δt is

$$\begin{equation} N_{{\rm coll}_{i,j}}=\frac{\pi \big(r_i^2 +r_j^2\big) v_{\rm rel} \Delta t}{V} N_iN_j; \end{equation} \tag{ 9 }$$

where v_rel is the relative encounter velocity, and V is the volume occupied by the swarm of particles.

Because a single disruptive collision, or accretion of that object onto a larger target will remove an object from a bin, Equation (9) does not represent the true number of collisions between bins i and j—there can only be as many disruptive collisions as there are objects in a bin to be disrupted. Therefore, in calculating collision rates, one must consider disruptions of bin i by other impactors k ≠ j, or if i is accreted onto a larger object. We adopt a probabilistic approach. We set the probability that an object from bin i is disrupted by an object from bin j (or accreted onto j if r_j > r_i), P_i,j, as the probability that the disruption of the target i was from an impactor from bin j, $P_{D_{i,j}}$, times the probability that an object from bin j was available to cause the disruption, P_j.

$P_{D_{i,j}}$ is simply given by the number of collisions of objects from bin i by those from bin j, given by Equation (9), divided by the number of objects in bin i. That is $P_{D_{i,j}} = \frac{N_{{\rm coll}_{i,j}}}{N_i}$.

The probability that impactor j is available can be found by setting $P_j = 1-\bar{P_j}$, where $\bar{P_j}$ is the probability that an object j is disrupted by or accreted by an object with radius r_k ≠ r_i. That is, $P_j = 1-\frac{1}{N_j}\sum _{k\ne i} N_{{\rm coll}_{j,k}}$. The true number of disruptive collisions from bin i by bin j is then

$$\begin{equation} N_{{\rm disrupt}_{i,j}} = N_i P_{D_{i,j}} P_j. \end{equation} \tag{ 10 }$$

Equation (10) predicts a smaller number of collisions for large objects compared to Equation (9). The difference imposes significant changes in the long-term collisional evolution when the number of disruptive collisions in a time step is ≳1% the number of objects in that bin. If the time step is small enough, then the effect is small.

For completeness, we apply a corrective term to Equation (9), $\big(1+\big(\frac{V_{{\rm esc}\,i,j}}{V_{\rm rel}}\big)^2\big)$, to account for gravitational focusing. For the impact velocities considered in this work, this term is negligibly small, and could be ignored without changing our conclusions.

2.3. Code Details and Nominal Parameters

2.4. Initial Conditions

The goal of this work is to determine the collisional evolution of the Kuiper Belt size distribution after the epoch of accretion. To that end, we consider a "typical" size distribution found from accretion simulations (Kenyon & Luu 1998; Kenyon & Bromley 2001, 2004; Stern 1996; Stern & Colwell 1997). That is, one which has a steep power law for the largest objects with q₁ ≳ 4. At some smaller size, r_b1, the slope becomes shallower, or breaks to q₂ ∼ −1 to 1 (Stern 1996; Stern & Colwell 1997; Kenyon 2002; Kenyon & Bromley 2004)—we utilize a sharp break for simplicity. The "break-radius" r_b1, typically 2 km for the Kuiper Belt (Kenyon 2002), corresponds to the size that objects which are smaller have undergone at least one disruptive collision. At some even smaller size, r_b2 ∼ 0.5 km, collisional equilibrium is reached, and the size-distribution slope turns up to the canonical q₃ = 3.5. Again we assume a sharp slope change. The collisional equilibrium slope typically only deviates at sizes small enough such that radiation effects remove objects from the belt (Thébault et al. 2003). The effect, however, is inconsequential on our results, so we do not consider radiation effects here.

The most recent measurements of the Kuiper Belt size distribution have determined that q₁ ∼ 4.8 (Fuentes & Holman 2008; Fraser & Kavelaars 2009) and break at a radius r_b1 ∼ 20–60 km. The steep slope implies a short accretion timescale (Kenyon 2002; Fraser et al. 2008) and has been preserved from the epoch of accretion (Durda & Stern 2000; Pan & Sari 2005). Thus, for our nominal simulations, we start with an initial large object slope equal to that observed in the modern day belt, i.e., q₁ = 4.8. We start with break radii and slopes inspired by the simulations of Kenyon (2002) and Kenyon & Bromley (2004) which represent the most complete simulations of accretion in the Kuiper Belt range to date, and for our nominal simulations we take q₂ = 0.0–2.0, q₃ = 3.5, r_b1 = 2 km, and r_b2 = 0.5 km. We discuss variations of these initial conditions in Section 4.

Accretion simulations of the region have demonstrated that the Kuiper Belt must have had a few orders of magnitude more mass (see Kenyon 2002, for a review) than currently observed (Trujillo et al. 2001; Gladman et al. 2001; Bernstein et al. 2004; Fuentes & Holman 2008), otherwise Pluto could not have achieved its size over the age of the solar system. The exact initial mass available at the late stages of accretion, however, is unknown. For our nominal simulations, we adopt a 45 M_⊕ Kuiper Belt and vary this mass to test the evolution for different planetesimal densities. We note here that 45 M_⊕ is at the large end of the mass thought to have existed in the primordial Kuiper Belt (Kenyon 2002; Gomes 2003). Thus, the rate of evolution produced in calculation with such a high mass should be interpreted as upper limits to the real evolution. The parameters for our nominal simulations are listed in Table 1.

Table 1. Nominal Simulation Parameters

Parameter	Adopted Value
Disruption Parameters
Qo	7 × 107 erg g-1
B	2.1 erg cm3 g-2
a	−0.45
b	1.19
X	0.55
Y	1
Cratering Parameters
α	10−9 s2 cm-2
fKE	0.1
ffrac	0.8
fcrat	0.05
qc	3.4
k	$\frac{9}{4}$
Miscellaneous Parameters
ρ	1.2 g cm-3
vrel	1 km s-1
rmin	0.01 cm
δ	1.1
Initial Size Distribution Parameters
q1	4.8
q2	0–2
q3	3.5
rb1	2 km
rb2	0.5 km

Download table as:  ASCIITypeset image

In Figure 1, we present the results of our nominal simulations at specific times. The expected behavior of the collisional erosion occurs at the large and small scales. For particles smaller than the initial r_b2, because these particles were initialized with the collisional equilibrium slope, the overall trend is to remain with that slope, and to decrease in mass as more particles are ground away. For the largest (r ≳ 150 km) objects, very few (if any) disruptive collisions occur. Accretion and cratering result in changes in the large object size distribution by less than ±∼1% over the age of the solar system.

Zoom In Zoom Out Reset image size

For objects with r ∼ 10 r_b1, a distinct divot forms in the size distribution. This behavior is caused by the initial break or turnover at r_b1; the decrease in slope causes an absence in the number of disruptors available to shatter objects with r ∼ r_b1. Thus, due to the change in the slope at r_b1, any objects capable of being shattered by impactors with r ∼ r_b1 experience an enhanced collisional disruption rate. This in turn reduces the disruption rate of objects capable of being shattered by the impactors with r ∼ 10 r_b1, which produces a distinct kink where the size distribution rolls away from the initial accretion slope q₁ into the divot.

The location of the divot can be roughly estimated by equating the collision energy of impactors at the accretion break radius r_b1 to the disruption energy of the bodies with the divot radius r_div. Ignoring the tensile strength term in Equation (1), and solving for r_div, we find

$$\begin{equation} r_{\rm div}=\left(4 f_{\rm dis} B \rho \right)^{\frac{-1}{3+b}} \; v^{\frac{2}{3+b}} \; r_{b1}^{\frac{3}{3+b}}. \end{equation} \tag{ 11 }$$

Equation (11) demonstrates that given the accretion break size, r_b1, the radius of the divot, r_div, depends on the strength parameters of the gravity term in Equation (1), B and b, the density of Kuiper Belt objects, ρ, the average interaction velocity, and the energy at which impacts transition from the cratering regime to the disruption regime, given by f_dis. The factor $\left(4 f_{\rm dis} B \rho \right)^{\frac{-1}{3+b}}$, however, is of order unity. Thus, the central divot radius depends mainly on the slope of the gravity strength scaling-law b, and the velocity at which objects interact. These conclusions are confirmed in Figure 2, where we present the results of our nominal simulations with different strength parameters.

Zoom In Zoom Out Reset image size

The rate at which the divot deepens on the slope change at r_b1; the greater the difference between q₁ and q₂, the greater the rate at which the divot depends (see Figure 1). The divot formation timescale can be roughly estimated from Equation (9). The divot is formed primarily from catastrophic collisions of objects larger than the initial break radius r_b1 where the initial size distribution is given by a power law with a slope q₁, i.e., $\frac{dN}{dr}=Cr^{-q_1}$, where C is a normalization constant. Thus, by inserting a power law for the number of smaller objects N_j in Equation (9), integrating overall sizes of object capable of disrupting objects at the divot, i.e., from r_b1 to r_div, and solving for time, we find that the timescale for divot formation t_div is given by

$$\begin{equation} t_{\rm div} = \left(\frac{N_c}{N_{r_{\rm div}}} \right) \frac{V}{\pi v C r_{\rm div}^{3-q_1}}\left[ \frac{1-f^{1-q_1}}{1-q_1} + \frac{1-f^{3-q_1}}{3-q_1}\right], \end{equation} \tag{ 12 }$$

where $\frac{N_c}{N_{r_{\rm div}}}$ is the fraction of objects with radius r_div that have been disrupted, and $f=\frac{r_{b1}}{r_{\rm div}}\sim 0.2$ (see Equation (11)). By adopting the nominal parameters, and saying a divot has formed when 99% of objects at the divot have been disrupted, i.e., $\frac{N_c}{N_{r_{\rm div}}}=0.99$, then the divot formation timescale is t_div ∼ 25 Myr. For this timescale estimate, we have ignored disruptions of objects with radii r_b1 ⩽ r ⩽ r_div. Thus, Equation (12) represents a lower limit of the true divot formation timescale, which is ∼50% longer if the break slope is q₂ = 0 and ∼3 longer if q₂ = 2 (see Figure 1).

Since the probability of a single object having a collision is approximately linearly dependent on the density of particles, the rate at which the divot forms is also linearly dependent on the density of particles (see Equations (9) and (12)). Minimal variation in the results occurs due to order of magnitude variations in the belt density at equivalent points in the simulations. This is demonstrated in Figure 3, where we compare the nominal simulation with a simulation using the same parameters with the density dropped by a factor of 10.

Zoom In Zoom Out Reset image size

Presented in Figure 4 are the results of our nominal simulations when the parameters relevant to the cratering model are changed. With certain strength models, namely those where cratering is responsible for removing a large fraction of object with r ≳ 10 r_b1, small waves can form at those sizes. But for the majority of the simulations, the roll-over is smooth up to the edge of the divot where the density of objects decreases very rapidly. Other than the maximum cratering radius parameter f_dis which affects the divot location, little variation in the size distribution for r > r_b2 results from variations in the cratering parameters. Changes do occur for r < r_b2. As there currently exists no measurement of the shape of the Kuiper Belt size distribution in this size range, we leave the study of these parameters to later works.

Zoom In Zoom Out Reset image size

The divot produced in our calculations could provide a natural explanation for the observed roll-over at r ∼ 25–60 km detected by Bernstein et al. (2004) and confirmed by Fuentes & Holman (2008) and Fraser & Kavelaars (2009). That is, the size distribution deviates away from the steep large-object accretion size distribution into a divot at r_div ∼ 10–15 km. If the observed roll-over is caused by a break left-over from accretion, then, from Equation (11), we find that the accretion break must have occurred at a radius, r_b1 ∼ 0.5 − 3.5 km, for our adopted strength scaling slope b = 1.19. The implied accretion size break r_b1 is consistent with that found in previous accretion models (Kenyon 2002).

Comparison of the ratio of objects larger and smaller than the kink radius, r_k, from the simulations to the roll-over observed in the Kuiper Belt can provide constraint on the mass of the belt at the end of accretion. Assuming that the accretion break slope is q₂ ∼ 0 and that the Kuiper Belt strength scaling is well represented by our assumed value b = 1.19, a 45 M_⊕ Kuiper Belt would produce a divot consistent with observations in ≳40 Myr if r_b1 = 2 km and ≳100 Myr for r_b1 = 1 km. The divot formation timescales are a factor of ∼2 longer for the initial accretion break slope q₂ = 2 and ∼2 shorter for the scaling law of Leinhardt & Stewart (2009), i.e., if objects are a factor of 3 weaker than suggested by Benz & Asphaug (1999). The time estimates also scale linearly with the density of planetesimals at the epoch of the dynamical excitation.

Our Kuiper Belt mass estimates are consistent with accretion simulations which suggest that the primordial Kuiper Belt had a mass ≳30 M_⊕ (Stern & Colwell 1997; Kenyon & Luu 1998; Kenyon 2002). If we assume that the mass depletion could not have occurred more than 1 Gyr after the Kuiper Belt was excited, then the minimum Kuiper Belt mass required to produce a divot compatible with the observations is ∼1–5 M_⊕.

The simulations of Charnoz & Morbidelli (2007) are the first attempt at coupling the dynamical and collisional history of the Kuiper Belt. Using a similar initial mass distribution (45 M_⊕ between 5 and 50 AU) as we have adopted in our simulations, they found that the mass depletion which has produced the low-density modern day Kuiper Belt must have been a result of combined collisional grinding and dynamical effects. Otherwise, the planetesimal populations external to the Kuiper Belt—the Oort cloud and scattered disk (see Gladman et al. 2008)—would be too anemic by at least an order of magnitude. They found that the mass depletion of the Kuiper Belt occurred over ∼100 Myr. They however adopted an initial size distribution that broke from a steep large-object slope to that of collisional equilibrium rather than the accretion model inspired distribution we adopt here, complicating comparison of the two works. Roughly speaking, however, the mass depletion timescale they found would be sufficient to produce a divot, if starting from an accretion-like size distribution. They also found that populations external to the Kuiper Belt underwent significantly enhanced collisional evolution compared to the Kuiper Belt. This suggests that the scattered disk and Oort cloud populations would have very few objects with 1 ≲ r ≲ 20 km.

Our results are consistent with some of the peculiarities of detected Kuiper Belt objects. Only one object with radius r ⩽ 15 km has been detected (Fuentes et al. 2009). Such an observation is consistent with a flat q ∼ 0 power law for objects smaller than the current roll-over. It is suggestive, however, that there is a substantial depletion of objects with the similar size. This is consistent with the idea that the size distribution is not a power law for r ∼ 10–20 km, but rather has a shape like the divot produced in our calculations. The shape of the luminosity function presented by Fuentes et al. (2009) hints that a divot is detected at m(R) ∼ 26.5. The observational evidence of this, however, is extremely sparse. More observations of small r ∼ 10–40 km Kuiper Belt objects are required before the true size distribution shape is known.

We find that objects with sizes similar to the initial accretion break r_b1 do not suffer significant depletion on timescales necessary to form the divot. Our findings demonstrate that the size distribution should exhibit a depletion of objects for r ∼ 10 km compared to objects with r ∼ r_b1. Comparison of our results to that of Benavidez & Campo Bagatin (2009) demonstrates that the shape of the size distribution caused by collisional evolution is highly dependent on the initial size distribution; the break can be formed with different initial conditions, but the location, and interpretation of the break change. Our results demonstrate that the detection of a divot in the size distribution is distinct evidence for the existence of a break in the size distribution at the onset of collisional disruption.

Our calculations demonstrate that, with a detection of the remaining peak at, r_b1, and the divot center, r_div, the strength scaling slope b as well as the accretion break radius r_b1 can be inferred. This would provide a stringent test of numerical impact simulations of icy bodies, as well as a strong constraint on accretion models of the Kuiper Belt region. The current uncertainty on the roll-over size cannot place constraint on accretion or shattering physics of icy planetesimals, other than to say that the observations are consistent with the current understanding of these two processes.

Our findings suggest that the sky density of objects with r ∼ 1 km is ∼10⁶–10⁷ objects per square degree. Detection by serendipitous occultation of star light is currently the only method of detection sensitive to these small Kuiper Belt objects (KBOs). These exciting experiments have currently placed an upper limit of ∼10⁸–10⁹ objects per square degree at this size (Bickerton et al. 2008; Zhang et al. 2008; Bianco et al. 2009), and future experiments will ultimately determine if the roll-over is caused by the mechanism we have suggested here.

Our calculations, and the existence of the steep slope observed for the large object size distribution, suggest that the Kuiper Belt underwent a short period of accretion, where objects as large as Eris were produced. Some dynamical event truncated accretion by increasing interaction velocities to disruption. On a short ∼10–100 Myr timescale, a divot, or depletion of objects with r ∼ 10 km was produced before the dynamical scattering removed most of the mass from the system, freezing the shape of the size distribution, which has undergone minimal evolution to the current day.

Our calculations have provided the first direct link from the "typical" size distribution produced from accretion to the collisionally modified size distribution of the modern day Kuiper Belt. Our calculations have demonstrated that, given a size distribution similar to that produced by accretion, the likely result of collisional evolution in the Kuiper Belt is to produce a divot in the Kuiper Belt size distribution at r ∼ 10–15 km, and a roll-over from the large object accretion slope at r ∼ 25–60 km, compatible with the observed roll-over in the Kuiper Belt size distribution. We conclude that the size distribution is likely not a power law for objects smaller than the roll-over, and exhibits a dearth of objects at the divot location. Assuming that a divot has been formed, a measurement of the divot center radius would provide constraint on the strength scaling law in the gravity regime for KBOs, and constrain the size distribution of r ∼ 1 km objects left-over from the accretion processes in the early solar system.

From our calculations, we find that the Kuiper Belt must have had at least 1–5 M_⊕ of material after its objects were excited onto orbits with eccentricities comparable to those observed, otherwise the roll-over would not be produced on Gyr timescales. If the belt had a mass of ∼45 M_⊕, then a divot would form in ≳40 Myr.

From the results of our calculations, we predict a deep absence of objects at r ∼ 10 km. We also predict that the sky density of objects with r ∼ 1 km is ∼10⁶–10⁷ per square degree, 2–3 orders of magnitude larger than if the size distribution is a shallow power law for objects smaller than the roll-over, as suggested by Fraser & Kavelaars (2009).

A special thank goes to J. J. Kavelaars and Mike Brown without whose guidance this work would have never come to fruition. I also thank David Trilling for his useful advice on what I should be looking for. Finally, I thank Parker, Raggozzine, Schwamb, and Bannister for humoring me during my rants.

This project was funded in part, by the National Science and Engineering Research Council and the National Research Council of Canada. Support for this work was provided in part by NASA through a grant from the program HST-GO-11644 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc. under NASA contract NAS 5-26555. This research used the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency.

Here we present various tests to confirm the correct behavior of the collisional model. As done by Ohtsuki et al. (1990) and Wetherill (1990), the behavior of a collisional evolution model—at least in terms of accretion—can be assessed by comparing numerical simulations of coagulation to the analytic solutions of the discrete coagulation equation given by

$$\begin{equation} \frac{d}{dt} n_k = \frac{1}{2} \sum _{i+j=k} A_{ij}n_i n_j - n_k \sum _{i=1} A_{ik} n_i. \end{equation} \tag{ A1 }$$

Here, $\frac{d}{dt} n_k$ is the time rate of the change of objects in bin k from collisions adding objects to (left term) and removing objects from (right term) bin k summed over all bins. A_ij is the "kernel" and governs the probability of collisions between i and j. This is usually written in terms of mass bins i and j rather than the size bins we consider above. Any successful accretion or collision model must reproduce analytic solutions of Equation (A1). Discrepancies between the numerical and analytic solutions can highlight possible caveats about the results of the numerical models. Ohtsuki et al. (1990) have demonstrated that, in his constant mass bin formalism, the numerical solutions are typically accelerated compared to the analytic counterparts. It is worth noting that, in the mass-batch formalism of Wetherill (1990), the numerical solution lags the analytic one. As we use the constant mass bin formalism, the conclusions about the behavior of Ohtsuki's model are applicable here. We discuss this below. The acceleration decreases with the bin spacing, δ,⁷ and varies with each kernel used. Thus, the comparison of the numerical and analytic solutions can provide a metric for the largest allowable δ. These comparisons also provide a measure of how large a time step can be used before the incorrect evolution is produced.

The simplest kernel is a constant, i.e., A_ij = λ = constant, i.e., the collision probability does not depend on the target mass. This kernel produces smooth orderly growth, and the solution to the coagulation equation for this kernel, i.e., the number of objects in a bin, n_k, is $\frac{n_k}{n_o}=f^2 \left(1-f\right)^{k-1}$, where f = (1 + η/2)⁻¹, and η = λn_ot is the dimensionless time. n_o is the total number of objects at time zero, all of which are of the smallest mass.

The analytic solution and the numerical solution of the coagulation equation with the constant kernel are presented in Figure 5 utilizing a δ = 1.1, n_o = 10¹⁸ particles,⁸ 128 bins, and by not allowing the number of objects in a bin to change by more than 0.1% per time step. As can be seen, the numerical solution matches the analytic solution well. The discrepancies present are similar to those found by Ohtsuki et al. (1990). The numerical solution produces more large objects than in the analytic one. The difference, however, is almost negligible and increases for larger δ.

Zoom In Zoom Out Reset image size

Another kernel for which an analytic solution exists is for collision probabilities proportional to the sum of colliding masses, i.e., A_ij = β(i + j). The analytic solution for this kernel is $n_k=n_o \frac{k^{k-1}}{k!} f\left(1-f\right)^{k-1}e^{-k(1-f)}$, where f = e^−η, and η = βn_ot is the dimensionless time. The analytic and numerical solutions for this kernel are presented in Figures 6 and 7 for n_o = 10¹⁸ particles, 128 bins, and δ = 1.08.⁹ This kernel creates large objects more rapidly than the constant kernel. As such, the numerical solution is more sensitive to the choice in δ than for the constant kernel. As shown by Ohtsuki et al. (1990), a finite δ in the numerical solution creates an acceleration in the coagulation, as well as an overabundance of large objects compared to the analytic solution. Our model exhibits virtually identical behavior. We find that the acceleration is nearly a constant over the simulations (see Figure 7). That is, for δ = 1.08 the shape of the numerical solution best matches the analytic solution at a time ∼12% earlier than expected. This factor is constant over the simulations, and scales with the choice in δ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The kernel proportional to the product of the colliding masses, i.e., A_ij = γ(i*j) produces runaway growth at η = 1 (again η = γn_ot is the dimensionless time). As such, the kernel is difficult to solve analytically, as the time step must vary with mass of the run-away body to recreate the proper evolution. The analytic solution for the non-runaway bodies is $n_k=\frac{n_o (2k)^{k-1}}{k!k}\big(\frac{\eta }{2}\big)^{k-1}e^{-k\eta }$, and the mass of the runaway body is $m=n_o \exp \big[-\int \sum _{k=1}^{N} k^2 \frac{n_k}{n_o} d\eta ^{\prime }\big]$ (Wetherill 1990).

The code presented here was intended for collisional evolution calculations, where the primary process is disruption of objects to ever smaller sizes. The code inefficiently handles the case of runaway growth; the code takes many more steps than necessary, causing the runtime to increase substantially as runaway continues. Our code was run only long enough to demonstrate that the numerical solution produces the correct mass of the runaway body, and distribution of smaller bodies (see Figure 8). Future versions of the software will include modifications to handle runaway growth more efficiently. But as the simulations presented above have minimal accretion of even the largest objects, inclusion of these modifications was unnecessary, and omitted for this work.

Zoom In Zoom Out Reset image size

Presented in Figure 8 is a comparison of the numerical and analytic solutions to the coagulation equation before and after runaway begins using n_o = 10¹⁸ and δ = 1.06, and allowing the number of objects per bin to vary by 0.1% per time step. Much like the previous examples, the numerical solution produces an accelerated evolution compared to the analytic solution. The acceleration in the numerical solution of runaway growth is not constant with time, as is found in other numerical models (see for instance Kenyon & Luu 1998). The acceleration is most pronounced for short times (∼20%whenη = 0.4), and decreases with time.

In the numerical solution, runaway growth occurs at η = 0.976, slightly early compared to the analytic solution. The post-runaway numerical solutions are presented in Figure 8 at η = 0.98 and 1.01, and very accurately match the analytic solutions at η = 1.0003 and 1.0006, respectively. The mass of the runaway body in the numerical solution at these two times is 3.5 × 10⁻⁴, and 5.6 × 10⁻⁴ of the total system mass. The mass of the runaway body in the analytic solution at the times where the analytic and numerical solutions match is 3.5 × 10⁻⁴, and 6 × 10⁻⁴ of the total system mass, almost identical to the numerical solution. The results presented in Figure 8 demonstrate that our code accurately reproduces the analytic solution to the coagulation equation for the case of runaway growth, before and after runaway commences.

The solution to this kernel is much more sensitive to the choice in δ. For δ ⩾ 1.1, the acceleration verges on unacceptably large, and the mass of the run-away body is not correctly reproduced. Therefore, δ = 1.1 is an approximate upper limit on the bin spacing for simulations in which runaway growth occurs, and which utilizes our model.

Here we present tests of the collisional evolution considering various details of our model. Specifically, we test the effects of time step on our collisional evolution calculations and the numerical solutions to the coagulation equation. Additionally, we test our adopted collisional equilibrium forcing over the smallest bins.

The choice in time step modifies the results of the calculations; too large and the results will not accurately produce the correct evolution. While large time steps may repeat the general behavior of simulations which use smaller time steps, significant stochastic variations in the results occur if the time step used is too large. It was found that not allowing the number of objects in a bin to change by more than 1% per time step was sufficient to ensure accuracy in our collisional evolution simulations as well as the solutions to the coagulation equation for kernels A = λ and A = β(i + j) (see Figures 6, 7, and 9). Further tightening this condition did not change the results by more than a few percent, while relaxing this condition further produced incorrect results for all numerical solutions. Note: none of the numerical solutions of the coagulation equation include a probabilistic treatment of the collision rates as given in Equation (10). Therefore, these solutions suffer the most from larger time steps than do the main simulations presented here.

Zoom In Zoom Out Reset image size

For the simulations in which collisional equilibrium was forced for the smallest bins, we tested over how many bins this should be applied. The results are shown in Figure 10 where we present the evolution of a system already in collisional equilibrium, but forcing equilibrium on the smallest 30, 40, and 50 bins.

Zoom In Zoom Out Reset image size

If equilibrium is forced over 30 bins only, a large mass build-up occurs at the small-size end. This build-up increases in magnitude over time, and produces unrealistic results. If collisional equilibrium was forced over 40 and 50 bins, the mass build-up no longer occurs, and the smallest object becomes collisionally depleted, as expected. The minimum number of bins required to prevent the mass build-up varies with the cratering and fragment distributions chosen for the simulations. The rate of depletion does depend on the number of bins over which collisional equilibrium is forced. The results, however, are similar for objects more than ∼2 orders of magnitude larger than the smallest bin. This effect will not affect the main conclusions drawn above about the large object size distribution.

As demonstrated by Thébault et al. (2003), proper treatment of the smallest objects requires modeling of the radiation forces which are ultimately responsible for removing the dust out of the region. Such a treatment is beyond the scope of this work. We choose to force collisional equilibrium over the smallest 40 bins, as a compromise between having an unrealistic small object mass buildup, and potentially removing mass too quickly from the system. All conclusions drawn above have the caveat, that the small object size distribution might not be correctly accounted for here.