{'fraction_non_alphanumeric': 0.030757112582284654, 'fraction_numerical': 0.020423472928114625, 'mean_word_length': 4.384428960920933, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Collisional Divot in the Kuiper belt Size Distribution

This paper presents the results of collisional evolution calculations for the Kuiper belt starting from an initial size distribution similar to that produced by accretion simulations of that region - a steep power-law large object size distribution that breaks to a shallower slope at r ~1-2 km, with collisional equilibrium achieved for objects r ~0.5 km. We find that the break from the steep large object power-law causes a divot, or depletion of objects at r ~10-20 km, which in-turn greatly reduces the disruption rate of objects with r>25-50 km, preserving the steep power-law behavior for objects at this size. Our calculations demonstrate that the roll-over observed in the Kuiper belt size distribution is naturally explained as an edge of a divot in the size distribution; the radius at which the size distribution transitions away from the power-law, and the shape of the divot from our simulations are consistent with the size of the observed roll-over, and size distribution for smaller bodies. Both the kink radius and the radius of the divot center depend on the strength scaling law in the gravity regime for Kuiper belt objects. These simulations suggest that the sky density of r ~1 km objects is ~10^6-10^7 objects per square degree. A detection of the divot in the size distribution would provide a measure of the strength of large Kuiper belt objects, and constrain the shape of the size distribution at the end of accretion in the Kuiper belt.


Introduction
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) 1 There are three properties which describe the general shape of the Kuiper belt size distribution: • the existence of the largest objects, Eris and Pluto • the large object size distribution for D 100 km which is well characterized by a steep power-law dN dR ∝ r −q with q ∼ 4.8 (Gladman et al. 2001;Petit et al. 2006;Fraser et al. 2008) • 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 a 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.

The Model
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 2 . 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. 3 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 = 1 4 m i m j v rel m i +m j . In the following sections we describe the details of the model and our prescription for the collisional outcomes.

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 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 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.  have confirmed the findings of Benz & Asphaug (1999) and from their simulations have determined that X ≈ 0.48. The small difference in X found by  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 were 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 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 where m(r) is the mass of an object with radius r. Equation 4 can be solved for q D to give As γ > 0, for our fragment model, 1 < q D < 4.
For impacts with energy insufficient to disrupt the target, we considering a cratering model where the excavated crater mass is given by where α is the crater excavation coefficient, which depends on the material properties of the target and impactor, and f f rac is the fraction of kinetic energy that went into shattering the target material. We consider values of α = 10 −8 − 10 −9 s 2 cm −2 (Petit & Farinella 1993).
We assume that the mass fraction of ejected material with velocity greater than v is given by where k = 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 min = 2 9 f KE αf f rac , with the fragment mass escaping the body given by 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 Simulations of Benz & Asphaug (1999);  and  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 are just disruptive, as a fraction of Q * D , is given by f dis = 1 X 1 − (1 − f 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 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.

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 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 The probability that impactor j is available can be found by setting P j = 1 −P j whereP 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 − 1 N j k =i N coll j,k . The true number of disruptive collisions from bin i by bin j is then Equation 10 predicts a smaller number of collisions for large objects compared to Equation 9. The difference imposses 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, 1 + to account for gravitational focusing. For the impact velocities considered in this work, this term is negibibly small, and could be ignored without changing our conclusions.

Code Details and Nominal Parameters
In any given time-step, the number of disruptive collisions per bin is first calculated with Equation 10, using an adaptive time-step, such that the number of objects being removed from a bin -either from disruption, or accretion and cratering -is no more than 1% the number of objects in the bin before the time-step. This minimizes code run-times while preserving the correct collisional evolution (see Appendix A). If the probability of collision from Equation 10 is less than unity, then a single collision occurrence is decided randomly. ie. the collision occurs if R < N i P Di,j P j , where R is a randomly generated number with 0 ≤ R ≤ 1.
Formally, if the number of collisions of a bin is less than 10 8 per time-step, the number of collisions is kept as an integer value, and a random number generator is used to determine if the number is rounded up or down. Wetherill & Stewart (1989) has demonstrated that requiring an integer number of collisions per time-step when the number is less than 10 8 is a sufficient condition to produce the correct time averaged collision rates.
For each non-disruptive collision between a large object i and a small object j, the resultant body has a mass m(r i ) + m(r j ) − M esc . An object is removed from bin j, and if the resulting mass is sufficiently large, a body is removed from bin i and added to the bin with the appropriate bin k. For these collisions, their is a difference in mass between the resultant body, and the average mass of the bin it was promoted to. This mass, accreted by bin k, M acc,k (negative if cratering occurs) is accumulated over multiple time-steps. When the accreted (or cratered) mass becomes large enough, objects can "grow" (or shrink) to subsequent bins, ie. when M acc,k > M k±1 − M k (sign set according to accretion or cratering) objects are moved randomly from bin i to k ± 1. This is done in such a way as to reflect the fact that the accreted mass is occurring over all objects in bin k and not just one. Thus, M acc,k |M k±1 −M k | objects are added to bin k ± 1 from bin k when M acc,k N k > R where R is a randomly generated value with 0 ≤ R ≤ 1. This process preserves total system mass.
We consider bins increasing logarithmically in radius. That is r i = r i−1 δ where δ > 1. For our standard simulations, we consider δ = 1.1. It has been shown that a finite δ produces an acceleration in the results (Wetherill (see 1990), and Appendix A). THus, our conclusions drawn here have the caveat that the evolution might occur in nature slightly slower than in our simulations.
Campo Bagatin et al. (1994) has demonstrated that a finite minimum size bin will cause substantial wave-like deviations of the small size distribution away from the expected collisional equilibrium distribution; the smallest particles have no smaller impactors, and are in over abundance compared to larger particles. This over abundance creates an increased collision rate for larger particles, and so on. To prevent such a behavior, we enforce a collisional equilibrium slope for the smallest bins by removing excess objects from those bins. While this method is not the most accurate way of correcting for the finite bin minimum (see Kenyon et al. 2008, for a discussion), it is trivial to implement, and ensures accurate enough small object treatment to not affect our conclusions about the 1 km size distribution. It was found that enforcing collisional equilibrium over the smallest 30 − 40 bins is necessary to prevent the over abundance of the smallest objects (see Appendix A).
Our calculations were performed on an NVIDIA GTX 280 graphics card using the CUDA programming environment 4 . These cards provide 30 processors, each containing 8 processor cores, all running in parallel, and excel at highly parallel calculations, with speeds up to 50 times higher than typical multi-core CPUs. Due to the nature of the calculations however, we are limited to 2 n bins. For our nominal calculations, we use 256 bins with δ = 1.1 for our nominal simulations, allowing us to probe more than 10 orders of magnitude in radius. We considered larger δ with fewer bins to determine the effects of bin width and minimum bin size on our results. We found that changing the bin size did not produce a noticeable affect on our conclusions considering the range of results caused by the uncertain shattering and cratering models. Thus, we adopt δ = 1.1 for all the calculations we present here. We wish to model collisional evolution that occurs in dynamical conditions typical of the modern day Kuiper belt. To that extent, we roughly model the Kuiper belt as a single annulus, with inner and outer radii of 30 and 60 AU. We set a scale height of 20 AU similar to the scale height of the current Kuiper belt. We note that this scale height is significantly larger than that achieved in standard accretion calcuations (see Kenyon 2002, for discussion). Such a scale height can only be achieved from scattering dynamics not included in those accretion models, and which is necessary to shape the primordial Kuiper belt to its current state (see Morbidelli et al. 2008).
It has been shown that accretion calculations which utilize a single annulus, as done here, produce different results than those which utilize many annuli (Kenyon 2002). The calculations we present here are intended to show the collisional evolution of an excited Kuiper belt. Standard multiannulus codes cannot accurately represent the dynamical structure of the Kuiper belt and therefore do not accurately model the collisional evolution of this region (Benavidez & Campo Bagatin 2009). Rather, a proper treatment would involve a combined N-body and collisional evolution calculation, as done in Bottke et al. (2005) and Charnoz & Morbidelli (2007). Given that the dynamical history of the Kuiper belt is not yet known, such a calculation is beyond the scope of this work. As such, the calculations presented should only be interpreted as rough approximations to the true collisional evolution that occurred in the Kuiper belt.
We set v rel = 1 km s −1 which Dell' Oro et al. (2001) has demonstrated is a typical collision velocity in the modern Kuiper belt. At this velocity, targets of r 5 km will be disrupted in collisions with size ratio larger than ∼ 1 : 10 for the impactor and target. Collisions with smaller impactors will typically cause cratering, and accretion will occur -slowly -for only the largest (r 250 km) objects.
For our simulations we define two strength models. For the "strong" model, we adopt the strength parameters of water-ice for a 0.5 km s −1 collision velocity presented by Benz & Asphaug (1999), with Q o = 7×10 7 erg g −1 , B = 2.1 erg cm 3 g −2 , a = −0.45, and b = 1.19.  however, has demonstrated that the above destruction criterion is likely a factor of ∼ 3 too high for objects in the Kuiper belt. Thus, we adopt the strength parameters suggested by , which have the same slopes a and b, as the strong model, but with a factor of 3 less disruption energy at all sizes for our "weak" model. We utilize the strong model for our nominal simulations, and discuss the effects of using the weak model in Sections 3 and 4.
We set the maximum crater size, f Crat = 0.05, and partition the collision energy into fracturing and kinetic energy as f f rac = 0.8 and f KE = 0.1. Following Petit & Farinella (1993), for our nominal model, we set the cratering velocity fragment k = 9 4 , the cratering efficiency α = 10 −9 s 2 cm −2 corresponding to "strong" surfaces or low cratering efficiency, and the fragment size distribution slope q c = 3.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 2001Stern 1996;Stern & Colwell 1997). That is, one which has a steep power-law for the largest objects with q 1 4. At some smaller size, r b1 , the slope becomes shallower, or breaks to q 2 ∼ −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 = 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 1 ∼ 4.8 (Fuentes & Holman 2008;Fraser & Kavelaars 2009) and breaks 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, ie. q 1 = 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 2 = 0.0 − 2.0, q 3 = 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). This, 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.

Results
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 results in changes in the large object size distribution by less than ± ∼ 1% over the age of the Solar system.
For objects with r ∼ 10 r b1 , a distinct divot forms in the size distribution. This behavior is caused by the initial break or turn-over 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 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 1 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 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 (4f dis Bρ) −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.
The rate at which the divot deepens depends on the slope change at r b1 ; the greater the difference between q 1 and q 2 , the greater the rate at which the divot deepens (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 slope q 1 , ie. 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 over all sizes of object capable of disrupting objects at the divot, ie. from r b1 to r div , and solving for time, we find that the timescale for divot formation t div is given by where Nc Nr div is the fraction of objects with radius r div that have been disrupted, and f = r b1 r div ∼ 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, ie. Nc Nr 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. The true divot formation timescale is ∼ 50% longer if the break slope is q 2 = 0 and ∼ 3 longer if q 2 = 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 occur 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.
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 10r 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.

Discussion
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 2 ∼ 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. r b1 = 2. The divot formation timescales are a factor of ∼ 2 longer for the initial accretion break slope q 2 = 2 and ∼ 2 shorter for the scaling law of , ie. 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 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 rollover size cannot place anymore 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 6 − 10 7 objects per square degree. Detection by serendipitous occultation of star light is currently the only method of detection sensitive to these small KBOs. These exciting experiments have currently placed an upper limit of ∼ 10 8 − 10 9 objects per square degree at this size (Bickerton et al. 2008;Zhang et al. 2008) 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, suggests 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.

Conclusions
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 a divot has been formed, then 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 6 − 10 7 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).

Acknowledgements
A special thanks goes to JJ Kavelaars and Mike Brown, without who's guidance, this work would likely have never come to fruition. I would also like to thank David Trilling for his useful advice on what I should be looking for. Finally I would like to 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.

A. Appendix A
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 Here, d dt n k is the time rate of 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) has demonstrated that, in his constant mass bin formalism, the numerical solutions are typically accelerated compared to the analytic counter parts. It is worth noting that, in the massbatch formalism of Wetherill (1990), the numerical solution lags the analytic one. As we use the constant mass bin formalism, the conclusions about the behaviour of his model are applicable here. We discuss this below. The acceleration decreases with the bin-spacing δ, 5 and varies with each kernel used. Thus, comparison of the numerical and analytic solution 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, ie. A ij = λ = constant, ie. the collision probability doesn't depend on target mass. This kernel produces smooth orderly growth, and the solution to the coagulation equation for this kernel, ie. the number of objects in a bin, n k , is n k no = f 2 (1 − f ) k−1 where f = (1 + η/2) −1 and η = λn o t 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 18 particles 6 , 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 δ.
Another kernel for which an analytic solution exists is for collision probabilities proportional to the sum of colliding masses, ie. A ij = β (i + j). The analytic solution for this kernel is and η = βn o t is the dimensionless time. The analytic and numerical solutions for this kernel are presented in Figures 6 and 7 for n o = 10 18 particles, 128 bins, and δ = 1.08. 7 This kernel creates large objects more rapidly than the constant kernel. 5 As done in the main text, we still refer to δ here in terms of object size. Our δ is the cube of that used in the coagulation equation solution.
6 it was found that for no 10 10 , the results were independent of the choice in no 7 chosen to ensure the largest mass bin is larger than the largest object in the analytic solution.
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 over-abundance 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 δ.
The kernel proportional to the product of the colliding masses, ie. A ij = γ (i * j) produces runaway growth at η = 1 (again η = γn o t 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 re-create the proper evolution. The analytic solution for the non-runaway bodies is n k = no(2k) k−1 k!k η 2 k−1 e −kη , and the mass of the runaway body is m = n o exp − N k=1 k 2 n k no dη ′ (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. Figure 8 is a comparison of the numerical and analytic solutions to the coagulation equation before and after runaway begins using n o = 10 18 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.

Presented in
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 −4 and 5.6 × 10 −4 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 −4 and 6 × 10 −4 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 utilize our model.

B. Appendix B
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.
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.
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 become 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. Miscellaneous Parameters ρ 1.2 g cm −3 v rel 1 km s −1 r min 0.01 cm δ 1.1 Initial Size Distribution Parameters q 1 4.8 q 2 0 − 2 q 3 3.5 r b1 2 km r b2 0.5 km Solid line -500 Myr. a q 2 = 0. b q 2 = 1. c q 2 = 2. The arrows mark the approximate divot center radius r div predicted from Equation 11. NOTE: The radius bin-widths increase proportionally with the bin-centers. This correctly presents the number of objects in a bin, but causes the size distribution slopes to appear −1 shallower than the true slope.    Figure 6, but with the numerical solution presented at a time 12% earlier than the analytic solution.