The fossilized size distribution of the main asteroid belt

Icarus 175 (2005) 111–140 www.elsevier.com/locate/icarus The fossilized size distribution of the main asteroid belt William F. Bottke Jr. a,∗ , Daniel D. Durda a , David Nesvorný a , Robert Jedicke b , Alessandro Morbidelli c , David Vokrouhlický d , Hal Levison a a Department of Space Studies, Southwest Research Institute, 1050 Walnut St, Suite 400, Boulder, CO 80302, USA b Institute for Astronomy, University of Hawaii, Honolulu, HI 96822-1897, USA c Observatoire de la Côte d’Azur, B.P. 4229, 06034 Nice Cedex 4, France d Institute of Astronomy, Charles University, V Holešovičkách 2, 180 00 Prague 8, Czech Republic Received 19 May 2004; revised 4 October 2004 Available online 14 March 2005 Abstract Planet formation models suggest the primordial main belt experienced a short but intense period of collisional evolution shortly after the formation of planetary embryos. This period is believed to have lasted until Jupiter reached its full size, when dynamical processes (e.g., sweeping resonances, excitation via planetary embryos) ejected most planetesimals from the main belt zone. The few planetesimals left behind continued to undergo comminution at a reduced rate until the present day. We investigated how this scenario affects the main belt size distribution over Solar System history using a collisional evolution model (CoEM) that accounts for these events. CoEM does not explicitly include results from dynamical models, but instead treats the unknown size of the primordial main belt and the nature/timing of its dynamical depletion using innovative but approximate methods. Model constraints were provided by the observed size frequency distribution of the asteroid belt, the observed population of asteroid families, the cratered surface of differentiated Asteroid (4) Vesta, and the relatively constant crater production rate of the Earth and Moon over the last 3 Gyr. Using CoEM, we solved for both the shape of the initial main belt size distribution after accretion and the asteroid disruption scaling law Q∗D . In contrast to previous efforts, we find our derived Q∗D function is very similar to results produced by numerical hydrocode simulations of asteroid impacts. Our best fit results suggest the asteroid belt experienced as much comminution over its early history as it has since it reached its low-mass state approximately 3.9–4.5 Ga. These results suggest the main belt’s wavy-shaped size-frequency distribution is a “fossil” from this violent early epoch. We find that most diameter D
120 km asteroids are primordial, with their physical properties likely determined during the accretion epoch. Conversely, most smaller asteroids are byproducts of fragmentation events. The observed changes in the asteroid spin rate and lightcurve distributions near D ∼ 100– 120 km are likely to be a byproduct of this difference. Estimates based on our results imply the primordial main belt population (in the form of D < 1000 km bodies) was 150–250 times larger than it is today, in agreement with recent dynamical simulations.  2004 Elsevier Inc. All rights reserved. Keywords: Asteroids; Asteroids, dynamics; Collisional physics; Impact processes; Origin, Solar System 1. Introduction The main asteroid belt is a living relic, with ongoing collisional and dynamical evolution slowly obscuring traces left behind by planet formation processes. Despite this, the main belt retains critical clues that, properly read, can be used to * Corresponding author. Fax: +1-303-546-9687. E-mail address: bottke@boulder.swri.edu (W.F. Bottke). 0019-1035/$ – see front matter  2004 Elsevier Inc. All rights reserved. doi:10.1016/j.icarus.2004.10.026 discern the initial conditions and evolution processes that occurred during the planet formation epoch (e.g., the nature and mass of the solar nebula between Mars and Jupiter, the timing of Jupiter’s formation, the distribution of volatiles in the inner Solar System, the size distribution produced during runaway growth phase of planetary accretion, the scaling laws that control collisional evolution both during and after planetary accretion, the presence of planetary embryos inside Jupiter’s orbit, the migration of the giant planets and whether sweeping resonance ever crossed the main belt, the 112 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 degree of material mixing that occurred between the feeding zones, etc.). In order to interpret the record left behind in the main belt, we need to develop an understanding of how collisional and dynamical evolution has affected this region over the last 4.6 Gyr. A prerequisite for this kind of work is that we develop tools capable of modeling these processes as precisely as possible. At present, there are several different scenarios for modeling the dynamical evolution of the main belt. While these scenarios come in different flavors and have varying degrees of success at matching observational constraints (e.g., see Petit et al., 2002, for a recent review), they all share certain similarities. For example, they all conclude that a massive primordial main belt experienced a short but intense period of comminution during/after the accretion phase. Less then ∼10 Myr later, the primordial main belt quickly lost most of its mass via dynamical processes (e.g., sweeping resonances, excitation via planetary embryos), with the event presumably triggered by the formation and orbital evolution of Jupiter and/or the dispersion of the solar nebula (Petit et al., 2002). We refer to this event in this paper as the “dynamical depletion event,” or DDE for short. The relatively few planetesimals that stayed behind in the main belt region continued to undergo comminution at a reduced rate for billions of years. The net effect of these processes left the main belt in its current state. Although this main belt evolution scenario is considered state of the art, it has yet to be modeled in any modern collision code. A straightforward way to do this would be to combine the current generation of collision codes with the best available dynamical models. The modified code could then be used to track asteroid comminution over the last 4.6 Gyr. The problem with this approach is that it would make use of enough unknown parameters that obtaining unique (or even useful) results would be difficult to impossible. For example, the asteroid disruption and fragmentation routines used in current collision codes contain significant uncertainties (e.g., see Holsapple et al., 2002; Asphaug et al., 2002; and Davis et al., 2002, for recent reviews). We also lack a good understanding of both the conditions that existed in main belt during/after accretion and of the timing/nature of dynamical events that occurred in the main belt over the last 4.6 Gyr. To overcome these problems, we employ in this paper an alternative and more approximate approach that retains the essential aspects of the scenario described above but eliminates several model parameters. This method requires that we make two key assumptions about main belt evolution: (i) comminution among diameter D < 1000 km planetesimals in the main belt zone were very likely dominated by the same collision probabilities and impact velocities found there today (Petit et al., 2001), and (ii) an immense planetesimal population undergoing comminution for a short period of time is equivalent, for our purposes, to a much smaller population undergoing comminution for an extended period of time. As we will show, the application of (i) and (ii) allow us to bypass questions related to the initial size of the main belt population after accretion as well as the timing/nature of the DDE that scattered main belt material. To obtain the best possible results from our collision code, we go to some effort in this paper to derive state of the art model constraints. We do this by taking advantage of recent insights into the main belt population (e.g., Jedicke et al., 2002), asteroid disruption events (Benz and Asphaug, 1999; Durda et al., 2004a, 2004b) and the fragment size distribution produced by real asteroid breakup events (Tanga et al., 1999; Nesvorný et al., 2002a, 2003). We also make use of new estimates for the disruption frequency of asteroid families produced by the breakup of diameter D > 100 km bodies. We find these data are crucial to deriving a unique solution for the asteroid disruption scaling law. An additional way our work differs from other recent efforts is that we explore a wide range of initial conditions. Over the last several decades, nearly all main belt collisional models have limited themselves to initial populations that were more massive in every size range of interest than the observed one. Models of this type, however, may produce results that are inconsistent with the available constraints. For example, the starting conditions used by these models require the elimination of so many asteroids via comminution that the observed main belt cannot be reproduced without the use of disruption scaling laws that are highly discordant with those derived in laboratory and numerical experiments (e.g., Durda et al., 1998). A second example is that these models tend to produce far more asteroid families than those observed today. A third example stems from the fact that these models predict the main belt population should have decayed substantially over the age of the Solar System. If true, the near-Earth object (NEO) population, which is almost entirely replenished by the main belt (Bottke et al., 2000, 2002a, 2002b), should have decayed by a factor of ∼3 or more over the last ∼3 Gyr (Davis et al., 2002). Studies of the lunar and terrestrial cratering record, however, provide no evidence for such a decline; instead, they suggest that the NEO population (and hence the main belt population for D  30 km asteroids) has been relatively constant over this time (e.g., Grieve and Shoemaker, 1994; Shoemaker, 1998). Instead of following this path, our solution has its roots in several pioneering works on main belt evolution (Kuiper et al., 1958; Anders, 1965; Hartmann and Hartmann, 1968). Our best fit collisional model requires us to use an initial main belt population that contains relatively few bodies in the diameter D  120 km range. Accordingly, we argue that the current small body population (D  120 km) is predominantly composed of fragments produced by breakup events among larger asteroids (D
120 km). As we will show, this model produces results that are much more consistent with available constraints than previous efforts (e.g., Davis et al., 2002). The asteroid disruption scaling law derived from our best fit model is also remarkably similar to estimates pro- The fossilized main belt size distribution vided by numerical hydrocode experiments of asteroid collisions (Benz and Asphaug, 1999). Overall, our results lead us to conclude that the main belt size distribution is a “fossil” produced by numerous collisions that occurred early in Solar System history. They also explain why the main belt size distribution has been in steady state for the last ∼3 Gyr. Given the insights provided by this work, we believe we are now ready to attack the main belt evolution problem using increasingly realistic scenarios. A brief outline of our paper is as follows. In Section 2, we discuss some background on the main belt evolution problem and accomplishments (and limitations) of previous efforts. In Section 3, we present our collisional model. In Section 4, we describe our model constraints. In Section 5, we discuss both our approach to the problem and how we determined the nature of our starting population. In Section 6, we show our model results, where we use our collision code to derive the specific shape of the main belt size distribution after accretion ended among the D < 1000 km planetesimals as well as the shape of the scaling relationship controlling asteroid disruption. Some implications of our work are discussed in Section 7. Finally, in Section 8, we list our conclusions. 2. Background In this section, we review several issues related to the collisional modeling work and what insights we have gleaned from previous work. Those wishing to jump to a discussion of the collisional model should go to Section 3. 2.1. A brief history of planet formation and their effects on the main belt region To set the stage for our work, we briefly discuss how the asteroid belt was affected by planet formation. The classical view of planet formation in the inner Solar System, which involves the gradual coalescence of many tiny bodies into rocky planets, can be divided into four stages: (i) the accumulation of dust in the solar nebula into km-sized planetesimals, (ii) runaway growth of the largest planetesimals via gravitational accretion into numerous protoplanets isolated in their feeding zones; (iii) oligarchic growth of protoplanets fed by planetesimals residing between their feeding zones; and (iv) mutual perturbations between Moonto-Mars-sized planetary embryos and Jupiter, causing collisions, mergers, and the dynamical excitation of small body populations not yet accreted by the embryos (e.g., Safronov, 1969; Greenberg et al., 1978; Wetherill and Stewart, 1989; Weidenschilling et al., 1997; Chambers and Wetherill, 1998; Agnor et al., 1999; Weidenschilling, 2000; Petit et al., 2001; Kokubo and Ida, 2002). It is believed that runaway growth occurs over a timescale of 0.01–1.0 Myr while the latter stages required 10–100 Myr (for a recent review, see Weidenschilling and Cuzzi, 2004). 113 One powerful constraint that the main belt provides to planet formation models has to do with the putative large depletion of material in the main belt zone. While the surface density distribution of solid material in our protosolar nebula is currently unknown, estimates can be made using the assumption that the nebula contained just enough material of solar composition to form the planets at their current locations and compositions. These model results suggest that the Solar System’s surface density may have varied as r −1 to r −3/2 between Venus and Neptune, where r is heliocentric distance (Weidenschilling, 1977; Lissauer, 1987). Compared to this prediction, however, the amount of solid material in the main belt zone today is nearly 1000 times lower than our expectations. This is a serious problem when one considers that (i) we have no reason to think that the surface density of the nebula was anything but smooth and (ii) to produce large asteroids on timescales consistent with constraints provided by the meteorite record, the surface density in the main belt had to be at least 100 times higher than currently observed (Wetherill, 1989). To circumvent this mass depletion problem, many groups (e.g., see Petit et al., 2002) now assume that the primordial main belt originally contained M⊕ of material, enough to allow the asteroids to accrete on relatively short timescales An important implication of this assumption is that planetary embryos on the scale of Moon- to Mars-size bodies probably formed in the primordial main belt at the same time (e.g., Wetherill, 1992; Chambers and Wetherill, 1998, 2001; Petit et al., 2001). If so, these bodies may have dynamically excited planetesimals in the main belt region enough to initiate fragmentation. Thus, the primordial asteroid belt may have experienced an early collisional evolution phase where a significant fraction of the total main belt’s mass was still in D < 1000 km bodies. The elimination of bodies from the primordial main belt is likely to have come from a combination of collisions and dynamics, presumably triggered by the formation of Jupiter several Myr after the birth of the Solar System (Wetherill, 1992; Franklin and Lecar, 2000; Morbidelli et al., 2000; Chambers and Wetherill, 2001; Petit et al., 2001; Nagasawa et al., 2002; see review by Petit et al., 2002). This socalled “dynamical depletion event” (DDE) would have left the main belt in a state comparable to its current condition, with a total mass ∼5 × 10−4 M⊕ . The timing and strength of the mechanism that eliminated the mass is constrained by the presence of Vesta’s basaltic crust. If the main belt were massive for too long, Vesta’s crust would have been obliterated by collisions (Davis et al., 1985). We caution that much of this scenario is supposition. Although numerical modeling work makes a good case that planetary embryos once formed in the main belt, no one has yet proven it had to happen, nor have they shown that the main belt had to undergo an early phase of dynamical evolution. On the other hand, models that have attempted to investigate whether most of the main belt’s early mass was either never there (e.g., Kortenkamp et al., 2001) or that it was lost 114 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 by comminution alone (e.g., Durda et al., 1998) have not yet met with success. To probe these questions more deeply, we need to understand these collisional modeling efforts and determine their merits. Hence, we now review previous efforts to simulate the collisional evolution of the main belt. 2.2. Previous efforts to model the collisional evolution of the main belt Some of the earliest attempts to model the collisional evolution of the main belt were made by Anders (1965) (see also Hartmann and Hartmann, 1968). In Anders (1965), it was postulated that the initial main belt size distribution may have had a Gaussian or bell-shape, making most small asteroids in the current population byproducts of a relatively small number of collisional disruption events among D > 60 km asteroids.1 Anders’ conclusions were not based on collisional evolution models per se but rather on his efforts to derive the shape of the primordial main belt size distribution by reconstructing the parent bodies of asteroid families. This work produced several insights that we will invoke later in this paper. The first collisional evolution model was produced by Dohnanyi (1969), who predicted that an asteroid population in collisional equilibrium should eventually evolve to a steady state size distribution with a differential power law slope index of −3.5 (see also Hellyer, 1970, 1971; Paolicchi, 1994; Williams and Wetherill, 1994; Tanaka et al., 1996). To get this result, Dohnanyi’s model assumed fragmentation occurs for a fixed projectile-to-target mass ratio and that the self-gravity of the target was negligible. While his solutions provided important insights into how small body populations might evolve via a collisional cascade, his work did not use a realistic scaling law for disrupting asteroids. These scaling relationships are now commonly defined as the critical impact specific energy, Q∗D , the energy per unit target mass delivered by the projectile required for catastrophic disruption of the target (i.e., such that one-half the mass of the target body escapes) (Durda et al., 1998; Holsapple et al., 2002; Asphaug et al., 2002). The diameter of the projectile ddisrupt capable of disrupting a target asteroid (Dtarget ) is defined as:
 2 1/3 ddisrupt = 2Q∗D /Vimp (1) Dtarget , where Vimp is the impact velocity and the target and projectile are both assumed to have the same density. While there have been many different estimates for Q∗D over the years (e.g., Fig. 6 from Holsapple et al., 2002), several trends remain constant. Small asteroids are considered part of the “strength-scaling” regime, where the fragmentation of the target body is governed by its tensile strength. Large asteroids, on the other hand, are considered part of the “gravityscaling” regime, where fragmentation is controlled by the 1 Note that Anders ideas were based, in part, on Kuiper et al. (1958), who speculated the “bump” seen in observed population was primordial. self-gravity of the target. Values for Q∗D have been estimated using both laboratory experiments (e.g., see recent review by Holsapple et al., 2002) and numerical hydrocode experiments (e.g., see recent review by Asphaug et al., 2002). As summarized by Asphaug et al. (2002), Holsapple et al. (2002), and Davis et al. (2002), laboratory experiments and hydrocode modeling work suggests the transition between the regimes occurs in the range 100 < D < 200 m. The next important advance in collisional modeling was made by Davis et al. (1979, 1985), who created a selfconsistent 1-D collisional evolution code to track the evolution of the main belt from shortly after the runaway growth epoch (see also Davis et al., 1989). Davis et al. were the first to use a disruption scaling law that accounted for size dependent strength among smaller bodies and gravitational binding among larger asteroids. They were also possibly the first group to point out that accurate collisional models had to account for the mass depletion of the main belt, the intact crust of Vesta, the shape of the main belt size distribution, and the population of large families observed in the main belt. Unfortunately, the available data on these topics at the time their paper was written was limited. Despite these handicaps, results from Davis et al. (1979, 1985) provide us with several important insights (or paradoxes) into modeling the collisional evolution of the main belt over 4.6 Gyr. They are: Paradox 1. It is difficult for collisions alone to grind away the size distribution of the main belt predicted by accretion models without eliminating Vesta’s crust or producing size distributions that are inconsistent with the observed main belt size distribution. Paradox 2. Collisional evolution results generated using “realistic” Q∗D disruption laws have difficulty breaking up D > 100 km asteroids, mainly because these objects have significant self-gravity. Because these bodies cannot be eliminated by comminution, collision codes are driven toward initial main belt size distributions that contain roughly the same number of D > 100 km asteroids as those observed today. Results from planetary accretion codes, however, indicate that large main belt asteroids cannot grow over reasonable timescales unless the planetesimal disk in the main belt zone once held far more mass (e.g., Wetherill, 1989). To deal with these paradoxes, Davis et al. (1985) hypothesized that dynamical processes related to planet formation might have eliminated much of the main belt’s mass early in Solar System history. We agree; one of the goals of this paper is to provide evidence that this scenario is viable. Interestingly, Davis et al. (1979, 1985) also tested the Gaussian initial main belt population postulated by Anders (1965). Because their results suggested that these size distributions retain the same basic shape over 4.6 Gyr of comminution, they concluded that Gaussian populations were unlikely to serve as a reasonable starting point for main belt evolution. In hindsight, we believe that Davis et al. may have The fossilized main belt size distribution premature in dismissing the Anders (1965) scenario, in part because their results were based on Q∗D scaling laws that allowed D < 100 km asteroids to disrupt far more easily than current hydrocode models suggest (e.g., Benz and Asphaug, 1999). Nevertheless, the Davis et al. results appear to have been influential enough that all subsequent efforts to model main belt comminution have used initial size distributions with more mass in every size bin of interest than the observed population. We will return to this important issue in Section 5. Since Davis et al. (1979, 1985), several groups (Davis et al., 1989, 1994; Durda, 1993; Campo Bagatin et al., 1994, 2001; Durda and Dermott, 1997; Durda et al., 1998; Marzari et al., 1995, 1999; O’Brien and Greenberg, 2003; Cheng, 2004) have investigated the collisional evolution of the main belt using comparable codes and/or methods. (Some groups have used these codes to investigate small body populations in the outer Solar System as well; e.g., Stern, 1996; Charnoz and Morbidelli, 2003.) The primary motivation of these groups was to explore the effects of different Q∗D laws, different starting main belt populations, and, in some cases, different fragmentation laws (e.g., see Petit and Farinella, 1993; Campo Bagatin et al., 2001). One noteworthy advance produced by these works was the recognition that “bumps” in the main belt size distribution, one near D ∼ 3–4 km and one near D ∼ 100 km (Fig. 1), might be by-products of collisional evolution using “V”-shaped Q∗D laws. To create a wavy size distribution, we need Q∗D function to undergo an abrupt change. As asteroids increase in size, changes from negative Q∗D slopes in the strength regime to positive slopes in the gravity regime make asteroids just beyond the transition point more difficult to disrupt. Because these objects survive longer, an excess number of projectiles is created that is capable of disrupting still larger asteroids. This perturbation launches a wave into the size distribution. For a transition near D = 200 m between strength- and gravity-scaling disruption regimes, a bump is created near D ∼ 3–4 km (Campo Bagatin et al., 1994; Durda et al., 1998; O’Brien and Greenberg, 2003; see Davis et al., 2002, for a recent review). Some groups claim this pattern may also create the bump found near D ∼ 100 km (e.g., Durda et al., 1998). These issues will be discussed further in Section 5. Although all the collisional models mentioned above provide useful insights, we believe none has found a truly satisfying resolution to Paradoxes 1 and 2. To resolve these issues, we constructed a code that incorporates the insights described above into a collisional model. 3. Collisional evolution model To model the comminution in the asteroid belt, we use a modified version of the self-consistent 1-D collisional evolution model described in Durda and Dermott (1997) and Durda et al. (1998). In this paper, this code will be referred to as CoEM, which stands for Collisional Evolution Model. 115 Fig. 1. An representation of the incremental main belt size frequency distribution computed from a parametric representation of the absolute magnitude H distribution (Jedicke et al., 2002). The H bins have been transformed into D bins using Eq. (9) and a geometric albedo pv = 0.092. The dots show the position of our bin. The H data is based on the catalog of known asteroids for H < 12 and results of the Sloan Digital Sky Survey (SDSS) for H > 12 (see Ivezić et al., 2001). The main belt size-frequency distribution is wavy, with “bumps” near D ∼ 100 km and D ∼ 3–4 km. Using this population, the cumulative number of D > 1, 50, and 100 km asteroids is 1.36 × 106 , 680 and 220, respectively. Note that when these data are plotted on a log scale, the apparent slope is shallower by unity than the power law index shown in Eq. (7) (Colwell, 1993). 3.1. Collision model framework To start CoEM, we enter an initial main belt sizefrequency distribution where the population (N ) has been binned between 0.001 km < D < 1000 km in logarithmic intervals d Log D = 0.1. All particles in our bins are assumed to be spherical and are set to the same density. We set the bulk density of each body here to be 2.7 g cm−3 . This value is something of a compromise between the measured densities of several different groups: the average bulk densities of several multi-km C-type asteroids (∼1.3 g cm−3 ), the average bulk density of several multi-km S-type asteroids (∼2.7 g cm−3 ), and the grain densities of several different meteorite classes that may be more representative of the bulk densities of sub-km asteroids (2.2 g cm−3 for CI meteorites, 2.7 g cm−3 for CM meteorites, 3.5 g cm−3 for CV meteorites; 3.5–3.8 g cm−3 for ordinary chondrites) (Britt et al., 2002). Note that moderate changes to this density do not change our results. The characteristic size of the particles in each bin is determined from the total mass and number of particles per bin. CoEM computes the time rate of change in the differential population N per unit volume of space over a size range between diameter D and D + dD (Dohnanyi, 1969; Williams and Wetherill, 1994): ∂N (D, t) = −IDISRUPT + IFRAG − IDYN . ∂t (2) 116 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 Here, IDISRUPT is the net number of bodies that leave the bin between D and D + dD per unit time from catastrophic disruption (see Eq. (3) below). In CoEM, we assume that most of our disruption events are barely-catastrophic (i.e., 50% of the target body’s mass is ejected). This means we neglect both cratering events, which produce much less ejecta over time than barely-catastrophic disruption events (e.g., Dohnanyi, 1969; Williams and Wetherill, 1994) and highlyenergetic catastrophic disruption events, which are unlikely to occur (e.g., Love and Ahrens, 1997). Exceptions to this rule are made for the largest asteroids; see Section 3.5 for details. We define IFRAG as the number of bodies that enter the size interval per unit time that were produced via the fragmentation of larger bodies. Finally, IDYN is the number of bodies lost from the size interval via dynamical processes (e.g., dynamical excitation via planetary embryos, removal from the main belt via a resonance). As we will describe below, we employ a numerical strategy in this paper that allows us to overlook IDYN for now.2 Once we have the components in place, CoEM integrates Eq. (2) over a time range of interest. Output from CoEM includes the evolved size distribution, the time when each asteroid disruption took place, and the total number of asteroids breakups that have occurred in each size bin over time. We now describe various important components of CoEM in greater detail. 3.2. Asteroid collisional lifetimes The component IDISRUPT from Eq. (2) can be defined as: IDISRUPT = N τ (3) where τ is the collisional lifetime of a target body of diameter DT between D and D + dD. Assuming a projectile of diameter ddisrupt can barely disrupt DT , we can define the lifetime of the target body τ as: 1 1 = τ 4 DT N Pi (DT + d ′ )2 dd ′ , (4) ddisrupt where d ′ is the projectile diameter, ddisrupt computed using Eq. (1) and Pi is the “intrinsic collision probability” that projectiles will strike DT . In this model, Pi is defined as the probability that a single member of the impacting population will hit the target over a unit of time (Öpik, 1951; Wetherill, 1967; Greenberg, 1982; Farinella and Davis, 1992; Bottke 2 To avoid adding additional free parameters to CoEM, we purposely ignore the effects of Yarkovsky thermal drift that cause D  30 km asteroids to drift into resonances where they can escape the main belt (e.g., Bottke et al., 2002a, 2002b; Morbidelli and Vokrouhlický, 2003). Preliminary modeling results (O’Brien and Greenberg, 2004) and our own tests indicate that the Yarkovsky effect affects the multikilometer main belt population but not to such a degree that our conclusions are seriously impacted. and Greenberg, 1993).3 Gravitational focusing is generally neglected here because asteroid escape velocities are ∝ m s−1 whereas asteroid impact velocities are ∝ km s−1 (Bottke et al., 1994). For our model to work, we need Pi and Vimp values in both the past and present-day main belt. For the moment, we will concentrate on the current values. To compute Pi and Vimp in the present-day main belt, Bottke et al. (1994) took a representative sample of main belt asteroids (all 682 asteroids with D > 50 km) and calculated the collision probabilities and impact velocities between all possible pairs of asteroids, assuming fixed values of semimajor axis, eccentricity, and inclination (a, e, i). Note that Öpiklike codes like that in Bottke et al. assume the orbits can be integrated over uniform distributions of longitudes of apsides and nodes; this approximation is considered reasonable because secular precession randomizes the orientations of asteroid orbits over ∼104 yr timescales. After all possible orbital intersection positions for each projectile-target pair were evaluated and weighted, they found that main belt objects striking one another have Pi ∼ 2.86×10−18 km−2 yr−1 and Vimp ∼ 5.3 km s−1 . These values have been corroborated by several different groups and methods (Farinella and Davis, 1992; Vedder, 1998; dell’Oro and Paolicchi, 1998; Manley et al., 1998). To get Pi and Vimp in the past, we need to understand how dynamical evolution affected planetesimals shortly after the accretion phase. Depending on how the asteroid belt reached its current state, this could lead to a wide variety of values. To winnow the possibilities, we decided to runs tests on particles evolving in the best available model of main belt dynamical evolution from Petit et al. (2001). We caution, however, that Petit et al. assume that all planetary embryos in the inner Solar System reach their full size more-or-less simultaneously and that the effects of gas drag on planetesimals is negligible. While these calculations currently represent the state of the art, it is possible these approximations are too simplistic. Our results indicate that prior to the formation of Jupiter, planetary embryos excited primordial asteroids to Pi ∼ 2.5 × 10−18 km−2 yr−1 and V ∼ 6 km s−1 . The values are nearly the same as current values. To explain the similarity, we created randomly distributed test bodies within various regions of (a, e, i) space that were centered on the main belt and then computed their Pi and Vimp values. Our results indicate that ensembles of particles in (a, e, i) zones roughly equivalent to the main belt produce Pi and Vimp values comparable to main belt values; big changes in Pi and Vimp only occur when the (a, e, i) zone is stretched to high e, i values. This explains the Petit et al. (2001) simulations, where planetary embryos push planetesimals into a region of space only slightly larger than the observed main belt zone over 10 Myr. 3 An asteroid’s cross-section is usually defined as (π/4)D 2 , but here the T π value is included in Pi . The fossilized main belt size distribution Given these results, we decided it is a reasonable approximation in CoEM to assume that Pi and Vimp values have stayed more or less constant, except for short intervals, since accretion ended among D < 1000 km bodies. This allows us to ignore the timing and precise nature of the DDE on main belt comminution, and it simplifies the computation of Eq. (4). The validity of this approximation hinges on two issues. The first is that the Petit et al. (2001) scenario is a reasonable approximation of reality or that other main belt evolution scenarios allow some dynamical excitation to occur prior to the formation of Jupiter and/or the dynamical depletion of the main belt. The second is that the main belt’s dynamical depletion phase is short enough that high velocity collisions do not dominate the other phases of the main belt’s collisional history (Petit et al., 2001, 2002). If the latter turns out to be false, our model will still produce useful results if high impact velocities fail to make target asteroids significantly easier to disrupt. Note that the record from hydrocode simulations is supportive of this idea; results from Benz and Asphaug (1999) indicate that higher impact velocities can lead to a smaller fraction of the projectile’s kinetic energy coupled to the target body. These issues will be investigated more closely in an upcoming paper. A more detailed discussion of how varying impact velocities in CoEM could affect our results is given in Appendix A. 3.3. Modeling Q∗D in the main belt A serious impediment to obtaining accurate results with CoEM is our lack of knowledge about Q∗D in the main belt. Note that without an accurate Q∗D function, we cannot derive a reliable ddisrupt value from Eq. (1). As described in Holsapple et al. (2002) and Asphaug et al. (2002), there is considerable debate in the catastrophic disruption community over which values of Q∗D are appropriate for particular material properties, impact velocities, and object diameters. Our lack of knowledge about Q∗D is also thwarting progress in planetary accretion codes that need to model the transition period between accretion and fragmentation as accurately as possible. For these reasons, a primary goal of this paper is to test a range Q∗D functions and determine those solutions that provide the best fit to our model constraints.4 To describe Q∗D , we use the formalism described by Durda et al. (1998) and adopt a four-parameter hyperbolic representation. The scaling law is a rotated and translated hyperbola in log Q∗D and log D space that follows the form: Ex 2 + Fxy + Gy 2 + H = 0 (5) 4 We caution that using a single Q∗ function for all asteroid disrupD tions is an oversimplification because asteroids of various taxonomic classes (e.g., S-type, C-type) and internal structures (e.g., monolithic, fractured, shattered, rubble-pile, porous, non-porous; see Richardson et al., 2002) are believed to react differently to impacts. Our Q∗D function should therefore be thought of as a “global average” over all asteroid Q∗D functions. 117 with x = x − x0 , x = log D (km), y = y − y0 , and y = log Q∗D (erg g−1 ). For all Q∗D cases described in this paper, x0 = −0.753 and y0 = 2.10; our rationale for these values is described in Section 6. Like Durda et al. (1998), we assume our Q∗D functions pass through the normalization point Q∗D = 1.5 × 107 erg g−1 and D = 8 cm. Although this value was determined using laboratory impact experiments, differences in target materials and impact velocities suggest other values may be equally valid here. The importance of the slope of Q∗D in the strength regime (ss ) should not be underestimated. As demonstrated by O’Brien and Greenberg (2003), in a collisional steady state scenario, ss determines the power law index of the main belt size distribution for bodies in the strength regime: 7 + ss /3 . ps = (6) 2 + ss /3 Because our attention in this paper is focused on collisional evolution in the gravity regime, however, the normalization point in the strength regime is less critical. O’Brien and Greenberg (2003) point out that when the gravity-scaled component of the population is wavy, it follows the general trend of a power law that is independent of both the population index and the Q∗D law in the strength regime. We have confirmed this numerically; modest changes to the slope of Q∗D in the strength regime do not noticeably affect our results. We will investigate the main belt size distribution in the strength regime in a future paper. 3.4. Modeling the statistical frequency of catastrophic disruption events After ddisrupt is calculated using Eq. (1), CoEM computes the number of objects with ddisrupt < D < DT from the input size distribution. If ddisrupt happens to be smaller than the smallest bin available in CoEM (Dsmall ), which in our runs is set to 0.1 m, the number of projectiles is estimated by extrapolating the shape of the size distribution to values where ddisrupt < Dsmall . With all components in hand, CoEM computes the collisional lifetime τ for each size bin. The timestep for the evolution model is automatically set to be 10 times smaller than the minimum τ value. To remove disrupted bodies from the bins of our sizefrequency distribution, we, like Durda et al. (1998), created two complementary versions of CoEM. In the first version, what Durda et al. called the “smooth” version of his code (or what we call CoEM-SM), the number of bodies removed from each bin is directly taken from analytic expressions using τ , the timestep, and the number of available particles in each bin. Because this number is generally not an integer, it is possible to end up with fractional remnants of individual bodies in the bins (e.g., one bin might be left with 0.7 bodies). The advantages of CoEM-SM are speed and its deterministic nature; a size distribution inserted into CoEMSM with a specific set of input parameters will always produce the same result. As we describe below, these attributes 118 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 can be used to determine certain aspects of the main belt Q∗D function. The main disadvantage of CoEM-SM is that it does a poor job of simulating the stochastic nature of breakups/fragment production in the main belt. For example, at each timestep the code produces a small number of fragments from all non-empty bins. Because this procedure is so different from reality, where a sudden family-forming event can flood the main belt with fragments in an instant, we are unsure whether these results can be realistically compared to observations. For our second version, called the “stochastic” version (or what we call CoEM-ST), breakups are treated as Poisson random events, where integer numbers of particles removed (or not removed) from a size bin within a timestep according to Poisson statistics (Press et al., 1989). Because CoEM-ST mimics breakups in the main belt in a more realistic manner than CoEM-SM, its results can be more directly compared with observations. The main disadvantage of CoEM-ST is that it is not deterministic; a different seed for the random number generator will produce a different outcome. Therefore, results from CoEM-ST need to be treated in a statistical manner; to get a quantitative measure of how good a given set of input parameters reproduces observations, we perform numerous trial runs using different random seeds before comparing our results to observations. 3.5. Asteroid fragmentation To get IDISRUPT for Eq. (3), we need to compute the size distribution of the fragments produced by each catastrophic disruption event. Previously, collisional evolution codes similar to CoEM (e.g., Davis et al., 1985; Durda and Dermott, 1997; Durda et al., 1998; O’Brien and Greenberg, 2003) have distributed ejecta from disruption events into smaller diameter bins using the power-law size distribution: dN = BD p dD. (7) We define dN as the number of fragments between D and D + dD, B is a constant determined so there is only one object as large as the largest remnant, and p, the differential power law index, is determined from the parameter b, the fractional diameter of the largest fragment in terms of the parent body. If one assumes that the total mass of the fragments is equal to the mass of the target asteroid, p is fixed by the expression (Greenberg and Nolan, 1989; O’Brien and Greenberg, 2003):   4 b +4 . p=− 4 (8) b +1 Interestingly, O’Brien and Greenberg (2003) showed analytically that the choice of b does not affect the shape of a power-law size distribution in steady state. Previous work, however, suggests that the entire main belt has yet to achieve collisional equilibrium over 4.6 Gyr of comminution (e.g., Section 2). Thus, the stochastic breakup of D > 100 km asteroids have important short- and long-term effects on the evolution of the main belt. To assess the accuracy of these approximations, we compare the predicted b, p values from Eq. (8) to those produced by asteroid families. We use asteroid families here rather than the results of hydrocode simulations (e.g., Michel et al., 2001, 2003) because we are still learning how to relate the former to the latter. In our discussion below, we concentrate on large family members that are the least likely to have been affected by observational incompleteness and/or collisional evolution. Thus, small family members with shallow p values (e.g., family members with absolute magnitude H > 15; Morbidelli et al., 2003) are purposely excluded from this analysis. Table III from Tanga et al. (1999) shows mLR /mPB , the ratio of the mass of the largest remnant to the mass of the parent body, for 14 prominent asteroid families. Computing b from these values (i.e., b = (mLR /mPB )1/3 ) and examining plots of the fragment size distributions from Tanga et al. (1999), Zappalà et al. (2002), and Morbidelli et al. (2003), we find we can divide asteroid-family-forming breakups into two broad classes: barely-catastrophic disruptions, where half the mass is ejected away at escape velocity (mLR /mPB ∼ 0.5 or b ∼ 0.8), and super-catastrophic disruptions, where more than half the mass is lost (mLR /mPB ≪ 0.5 or b ≪ 0.8). Asteroid families like Adeona, Erigone, and Flora were produced by barely-catastrophic disruption events (Tanga et al., 1999). For these cases, Eq. (8) predicts that b = 0.8 should yield a fragment size distribution with a differential power law index of p = −3.13. The observed families, however, show a very different distribution. To estimate the fragment size distribution of the Adeona and Flora families, we took H < 15 data from each family (see Morbidelli et al., 2003) and combined it with the relationship D ∝ 10−H /5 (e.g., Fowler and Chillemi, 1992; see also Section 4.1). We find the Adeona and Flora families have p  −2 for H < 11 and p = −4.0, −4.6 for 11 < H < 15, respectively. Though Morbidelli et al. (2003) did not investigate Erigone, Fig. 11 from Tanga et al. (1999) shows it shares comparable p values. For all 3 cases, the transition point between the shallow and steep slopes occurs roughly at 1/3 the diameter of the largest remnant. Hence, barely-catastrophic disruption events produce fragment size distributions that are poorly represented by a single p value. Comparable results are found for the large Eunomia (DPB ∼ 284 km) and Hygiea (DPB ∼ 481 km) families, although neither has technically undergone a catastrophic disruption event (mLR /mPB values are 0.73 and 0.61, respectively). Like the three families described above, Eunomia and Hygiea have p  −2 for their largest family members and steeper p values for their smaller objects: the Eunomia family members with 11 < H < 15 have p = −4.2 and Hygiea family members with 12 < H < 15 have p = −4.8. The fossilized main belt size distribution Fig. 2. The fragment size frequency distribution produced by catastrophic disruption events in our code. In this plot, our ejecta was produced by the breakup of a D = 150 km target asteroid. The “Themis-style” fragment distribution is modeled after prominent families like Themis whose largest remnant is <80% the diameter of the parent body (see Tanga et al., 1999). For our code, we assume the largest remnant is 50% the diameter of the parent body, and that the incremental power law index between the largest remnant and fragments 1/60 the diameter of the parent body was −3.5. The p value for fragments smaller than this threshold was −1.5. The “Flora-style” fragment size distribution was modeled after the barely-catastrophic disruption events that produced families like Flora and Adeona. The largest remnant was set to 80% the diameter of the parent body. We chose slope changes to occur at 1/3 and 1/40 the diameter of the parent body. The p values, from the large diameters to the small diameters, are −2.3, −4.0, and −2.0. Once again, no single p value accurately characterizes these fragment distributions. Super-catastrophic events, on the other hand, produce fragment size distributions that are better characterized by a single p value. Examples include families produced by the catastrophic breakup of D > 200 km asteroids such as Themis (DPB ∼ 370 km) and Eos (DPB ∼ 220 km), with mLR /mPB values of 0.31 and 0.11, respectively. Our predicted b values for Themis and Eos using Eq. (8) are 0.68 and 0.48, which corresponds to p = −3.5 and −3.8, respectively. The p values measured from Eos and Themis’s absolute magnitude data between 7.5–8 < H < 15 are −3.5, close to the Eq. (8) values. In both cases, the diameter of the largest remnant is roughly half the size of the parent body. The same trends can be found for Gefion (predicted b = −3.9; observed b(11 < H < 15) = −3.6) and Dora (predicted b = −4.0; observed b(11 < H < 15) = −3.7), but not for Koronis, whose four largest members are comparable in size (predicted b = −4.0; observed b(9 < H < 15) = −3.4). Taken together, these comparisons suggest that fragment size distributions have a strong dependance on the individual nature of each impact event, and that great care must be given when applying any simple metric to this issue. Another potential limitation in using Eqs. (7) and (8) in CoEM has to do with their prediction that b extends from 119 the largest remnant to dust-sized particles. As pointed out by Tanga et al. (1999), geometrical considerations, the fact that observed asteroids have convex shapes, and conservation of volume arguments between parent bodies and their fragments drive fragment size frequency distributions to shallow power law slopes at small fragment sizes. This effect is also observed in laboratory impact experiments, where the quantity of small fragments produced by a high speed collision becomes limited well before the detection limit is reached (e.g., Durda et al., 2002). Thus, one must be careful in stretching the predictions of Eqs. (7) and (8) too far when computing how much mass goes into small bodies. Based on these results and our desire to match the observed fragment size distributions of asteroid families as closely as possible, we developed two fragment size distributions to treat catastrophic breakups in CoEM (Fig. 2). “Themis”-style fragment size distributions (FSDs) were developed for D > 150 km catastrophic disruption events like Eos and Themis, a few of the largest catastrophic disruption events observed in the main belt. Note that it is surprising that so many large families were produced via superrather than a barely-catastrophic disruption events; this may suggest that some of the largest asteroids were pre-shattered before disruption, which would make them easier to disrupt (e.g., Michel et al., 2003). For these distributions, we assume the largest remnant is 50% the diameter of the parent body. The differential power law index p between the largest remnant and fragments 1/60th the diameter of the parent body was set to −3.5. The change in slope roughly describes where the initial Eos and Themis family FSDs are expected to roll over according to geometrical arguments (Tanga et al., 1999). The p value for fragments smaller than 1/60th the diameter of the parent body was set to a very shallow value (p = −1.5). “Flora”-style FSDs, developed for barely-catastrophic disruption events, were chosen for all breakups among D < 150 km bodies. Our logic here was that (i) CoEM explicitly solves for barely-catastrophic disruption cases, (ii) the frequency of super-catastrophic disruptions cannot be easily determined from the observational record, and (iii) barely-catastrophic disruptions are more likely to occur than super-catastrophic disruptions. The values chosen for the Flora FSD (the diameter of the largest remnant was set to 80% the diameter of the parent body; slope changes occur at 1/3 and 1/40 the diameter of the parent body; p values, from the large end to the small end, were −2.3, −4.0, and −2.0) were produced as a compromise between FSD measurements from Adeona, Erigone, and Flora and our goal of producing a distribution comparable to the Themis-style fragments for D > 0.1 km. As a check, we tested how the transition between Florastyle and Themis-style FSDs affected our results (for the best fit cases described in Section 6 using CoEM-ST). We found few observable differences in the evolved main belt size distribution when D < 100 km breakup events with Flora-style FSDs were replaced with Themis-style FSDs. 120 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 More noticeable differences were found when we extended this replacement to 100 km < D < 150 km breakup events. The reason is that Themis-style FSDs decimate the parent body to such a degree that we get significant depletion in the 100 km < D < 150 km size range and an excess of fragments in the D < 100 km range. Because this trend would make it increasingly difficult to reproduce the observed main belt, we believe that Themis-style FSDs are unlikely to dominate breakup events in the 100 km < D < 150 km size range (though they still must play an important role). To resolve this issue, future models will need to include realistic FSD derived by hydrocode models. This completes our basic description of CoEM. In the next two sections, we will discuss our model constraints and our methods for including the insights described in the introduction. 4. Constraints on the collisional evolution of the main belt To obtain useful results with CoEM, we need to have accurate constraints. As we describe below, determining solid constraints was one of the more challenging aspects of our modeling effort, partly because of an ongoing dispute about the true shape and size of the main belt size distribution but also because our understanding of the observed asteroid families has significantly advanced during the writing of this paper. 4.1. Constraint #1: The main belt size frequency distribution The goal of our CoEM runs is to reproduce, as precisely as possible, the current main belt size frequency distribution. To get this, we first need a measure of the main belt absolute magnitude H distribution. Estimates of the main belt H distribution have been produced by several different surveys and instruments: The Palomar–Leiden survey (PLS; van Houten et al., 1970), IRAS (Cellino et al., 1991); Spacewatch (Jedicke and Metcalfe, 1998), the Sloan Digital Sky Survey (SDSS; Ivezić et al., 2001), and Subaru (Nakamura and Yoshida, 2002; Yoshida et al., 2003). As described by Jedicke et al. (2002), SDSS results produced a large number of main belt detections over a small number of images during a controlled observational campaign, enough to produce better statistical results than either the PLS or Spacewatch surveys. For this reason, Jedicke et al. (2002) argued the state of the art for the debiased main belt H distribution is the SDSS H distribution for H > 12 coupled with the set of known main belt asteroids with H < 12 (Table 1; see also Jedicke et al., 2002). The reason for this split is described below. The next step, namely turning the H distribution into a size distribution, is surprisingly problematic. The relationship between asteroid diameter D, absolute magnitude H , Table 1 The observed main belt model parameters H D dN Family dN 3.25 3.25 4.25 4.25 5.25 5.75 6.25 6.75 7.25 7.75 8.25 8.75 9.25 9.75 10.25 10.75 11.25 11.75 12.25 12.75 13.25 13.75 14.25 14.75 15.25 15.75 16.25 16.75 17.25 17.75 18.25 980.9 779.2 618.9 491.6 390.5 310.2 246.4 195.7 155.5 123.5 98.1 77.9 61.9 49.2 39.1 31.0 24.6 19.6 15.6 12.4 9.81 7.79 6.19 4.92 3.91 3.10 2.46 1.96 1.55 1.23 0.98 1.0 0.0 0.0 2.0 1.0 3.0 8.0 17.0 38.0 64.0 91.0 116.0 164.0 185.0 224.0 338.0 554.0 789.7 1548.0 2992.3 5671.8 10463.9 18630.7 31739.6 51398.5 78939.8 115400.3 162026.4 221080.1 296503.1 394278.9 – – – – 1.0 1.0 1.0 5.0 5.0 5.0 – – – – – – – – – – – – – – – – – – – – – Column 1 is the absolute magnitude H . Column 2 is central diameter of the bin assuming geometric albedo pv = 0.092 (see Eq. (9)). Column 3 is the incremental number in each bin (Jedicke et al., 2002). Column 4 is the number of observed asteroid families in each bin. and visual geometric albedo pv , can be written as (e.g., Fowler and Chillemi, 1992): 1329 D(km) = √ 10−H /5 . pv (9) Thus, once we have representative pv values for main belt asteroids in every size bin of interest, Eq. (9) can be used to turn our incremental H distribution into an incremental D distribution. Unfortunately, the distribution of visual albedos from IRAS data among main belt asteroids with D < 50– 100 km are incomplete and potentially biased (Veeder et al., 1989; Tedesco et al., 2002). The best available estimates for asteroid diameters and albedos come from the Supplemental IRAS Minor Planet Survey (SIMPS), who used IRAS data to examine 2228 different multiply observed asteroids (Tedesco et al., 2002). The mean and median albedos for D > 20 km asteroids in this set are 0.10 and 0.060, respectively, while the same values for D < 20 km asteroids are 0.126 and 0.090, respectively. The SIMPS albedo distribution peaks near 0.06 for D < 20 km, but it also has a long tail that stretches to values beyond 0.3. This wide range of The fossilized main belt size distribution possible albedo values makes it problematic to turn our H distribution into a size distribution. Given our unknowns, and after much experimentation, we decided a reasonable compromise strategy was to choose a representative albedo value, pv = 0.092, to transform the Jedicke et al. (2002) H distribution into a size distribution (Fig. 1 and Table 1). In doing so, we, like Morbidelli and Vokrouhlický (2003), assume the power law slopes of the H distribution produced by SDSS/Jedicke et al. (2002) do a good job at reproducing the shape of the real main belt size distribution. The only meaningful change made to the Jedicke et al. distribution was to include the observed asteroids for D > 300 km using the IRAS/color-albedo-derived diameters cited in Farinella and Davis (1992). To check our results for large asteroids, we computed the cumulative number of asteroids with D > 50 and 100 km. Our values, 220 and 680, respectively, are almost exactly those determined by Farinella and Davis (1992) using IRAS/color-albedo data (e.g., Tedesco et al., 2002). For smaller asteroids, we found our main belt size distribution should have N(D > 1 km) ≃ 1.36 × 106 . This value is a good match to several different estimates: Morbidelli and Vokrouhlický (2003) showed using a numerical simulation that ∼1.3 × 106 km-sized main belt asteroids are needed to keep the main belt source regions of the near-Earth object population constantly replenished via the Yarkovsky effect and resonances, while Tedesco and Desert (2002), using observations from the infrared satellite ISO, estimated that the main belt should have (1.2 ± 0.5) × 106 asteroids with D > 1 km. Our results are more discordant, however, with estimates from Ivezić et al. (2001, 2002), who, using SDSS data alone, report there should be 7 × 105 km-sized main belt objects. This result, while within the error bars of Tedesco and Desert (2002), is still lower by nearly a factor of 2 than other estimates; a discussion of this issue can be found in both Jedicke et al. (2002) and Morbidelli and Vokrouhlický (2003). We summarize this discussion below. Ivezić et al. (2001) claim the main belt should have 68,000 bodies with H < 15.5, a value significantly lower than reported by the Spacewatch survey (Jedicke and Metcalfe, 1998). The predicted number of asteroids with H < 14 by Ivezić et al. (2001, 2002) also lies significantly below the number of known asteroids. Jurić et al. (2002) claim this difference may have been produced by systematic errors that permeate the H values reported in the ASTORB catalog (see ftp://ftp.lowell.edu/pub/elgb/astorb.html) and in the Minor Planet Center database. This explanation, however, is unsatisfying: (i) the bright H < 12–13 asteroids are thought to be far less susceptible to systematic observational errors than dimmer objects, yet they suffer from the same discrepancy, and (ii) population estimates made solely from Spacewatch data should not suffer from the same systematic errors as those from other surveys, yet Spacewatch’s results are higher than predictions made from SDSS data. In addition, Morbidelli and Vokrouhlický (2003) report that 121 a main belt population with 7 × 105 km-sized objects would result in a flux of near-Earth objects from the main belt with H < 18 far lower than numerical modeling results would predict (e.g., Bottke et al., 2002a, 2002b). Thus, to overcome this apparent discrepancy yet still take advantage of the large quantity of superb SDSS data, Jedicke et al. (2002) merged the observed distribution with H < 12 with the SDSS absolute magnitude distribution with H > 12. We believe our size distribution estimates, while imperfect, are probably the best we can do until Spitzer space telescope data is used to produce a more complete set of asteroid diameters in the 1 km < D < 30 km range. 4.2. Constraint #2: Asteroid families Asteroid families provide critical constraints for main belt collisional evolution models. These remnants of catastrophic collisions (e.g., Zappalà et al., 2002) are identified by their clustered values of proper semimajor axes a, eccentricities e, and inclinations i (Milani and Knežević, 1994; Bendjoya and Zappalà, 2002; Knežević et al., 2002). Ideally, by cataloging a complete set of main belt families over a given size range, we can determine their disruption frequency, which in turn can be used to constrain Q∗D . Obtaining complete sets of asteroid families, however, is problematic for two main reasons: 1. Our current understanding of planet formation suggests it is probably futile to look for asteroid families formed early in Solar System history. While meteoritical evidence points to the idea that large asteroids did break up during this time (Keil et al., 1994; Bogard, 1995; Scott et al., 2001; Scott, 2002, 2004) dynamical modeling work suggests their families would have likely been dispersed via encounters with planetary embryos (Petit et al., 2001) and/or sweeping resonances produced by the dissipation of the gas component of the disk and/or the dynamical evolution of Jupiter (Franklin and Lecar, 2000; Nagasawa et al., 2002; see review by Petit et al., 2002). Even if such dramatic events did not occur in the main belt, it is unclear whether reliable proper (a, e, i) values could be computed in a system where the giant planets had not yet achieved their final configuration. Hence, we believe the asteroid families observed today had to form after the main belt reached its current dynamical state. It is plausible, perhaps even likely, that this state was not reached until some time after the Late Heavy Bombardment, which occurred (or ended) ∼3.9 Gyr ago (Hartmann et al., 2000; Levison et al., 2001; Gomez et al., 2004). This observational limit would also explain the surprising paucity of observable families in the current main belt (i.e., where are all the families that presumably were produced by Late Heavy Bombardment projectiles?). 2. Like piles of leaves on a windy day, asteroid families dynamically disperse over time. Family members injected 122 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 into dynamical resonances undergo secular changes in both eccentricity and inclination, while the Yarkovsky effect causes all D < 30 km asteroids to slowly drift inward or outward in semimajor axis (Bottke et al., 2001; Nesvorný and Bottke, 2004). In some cases, family members are driven into powerful resonances capable of increasing their eccentricity to planet-crossing orbits; these bodies are lost from the main belt (Farinella and Vokrouhlický, 1999; Bottke et al., 2000; see review by Bottke et al., 2002a, 2002b). Even close encounters between family members and large asteroids like (1) Ceres can produce distinct (a, e, i) changes over time (Nesvorný et al., 2002a; Carruba et al., 2003). Working together with collisions, these processes gradually erase the signatures of catastrophic collisions for all but the largest breakup events. These mechanisms help explain why there is limited observational evidence for small break-up events across the main belt. The precise time needed to spread asteroid families beyond recognition in different regions of the asteroid belt is unknown. Much depends on factors such as the size of the fragments, the characteristics of the break-up event, and the proximity of the family to mean motion and secular resonances. Broadly, we can say that the timescale for family erasure changes from region to region, with the inner main belt (2.1 AU < a < 2.5 AU) much more chaotic and diffusive than the central (2.5 AU < a < 2.8 AU) and outer main belt (2.8 AU < a < 3.2 AU) (Nesvorný et al., 2002b). Determining precise timescales for family dispersion in these regions will require an extensive campaign of numerical integration work and is beyond the scope of this paper. To bypass these problems while still producing useful constraints for CoEM, we concentrate here on large families that cannot have been dispersed since the Late Heavy Bombardment (Bottke et al., 2001; Nesvorný et al., 2003). Previous work has suggested that at least 9 main belt families were produced by the catastrophic (or near catastrophic) disruption of D > 100 km asteroids (Tanga et al., 1999). The size of the parent body in each case was estimated using a geometrical/numerical model that summed the masses of the observed family members and then extrapolated the family’s FSD to smaller sizes. While this technique provides us with many useful insights, it is limited in that it cannot easily constrain the amount of mass hidden below the observational completeness limit (D  10 km). For example, because super-catastrophic disruption events are expected to produce numerous small fragments (Michel et al., 2001, 2003; Durda et al., 2004a, 2004b), Tanga et al. may underestimate the number of observed D > 100 km breakup events, particularly if numerous family members are below the detection limit. The issue is how to constrain the amount of mass hidden below our current detection limits. One way to do this is to compare the observed families to results from numerical hydrocode experiments. For the former, we determined family membership using a cluster-detection algorithm to a set of over 100,000 asteroid proper elements found at the AstDyS node (Milani and Knežević, 1994; Knežević et al., 2002). Note that this method has been successfully used to identify the 5.8 Myr old Karin cluster and many other main belt families (Nesvorný et al., 2002c, 2003, 2005; Nesvorný and Bottke, 2004). For the latter, we computed the fragment size distribution produced by 161 numerical simulations of impacts into 100 km diameter asteroids using a 3-D smooth particle hydrodynamics (SPH) code (Benz and Asphaug, 1999) combined with the N -body code pkdgrav (Richardson et al., 2002; Leinhardt et al., 2000; Leinhardt and Richardson, 2002). The computational details of this modeling effort are described in Durda et al. (2004a, 2004b). Each target body was assumed to be an undamaged basaltic asteroid. The projectiles striking the target body were given impact velocities Vimp between 2.5–7 km s−1 , impact angles θ between 15◦ –75◦ , and diameters Dimp between 10– 46 km. Figure 3 shows a representative example of our results, where we compare the observed members of the Erigone family to a run from Durda et al. (Vimp = 7 km s−1 , θ = 30◦ , and Dimp = 18.5 km). We caution that these impact parameters may not be the true values; because the physical state Fig. 3. A comparison between the observed size distribution of the Erigone family and that of a numerical hydrocode run described in Durda et al. (2004a, 2004b). The Erigone family, plotted as open circles, was determined using the methods described by Nesvorný et al. (2005). The largest member of the family is (163) Erigone, a diameter D ∼ 73 km C-type asteroid. The solid line shows the results of a numerical impact experiment from Durda et al. (2004a, 2004b) where a D = 100 km undamaged basaltic asteroid was struck by a 18.5 km projectile at a velocity of 7 km s−1 and an impact angle of 30◦ . The diameter of the largest remnant from this impact experiment was set to the diameter of (163) Erigone, with the rest of the model FSD scaled accordingly. Note the fit between model and data is satisfactory until D < 4 km, where the observed size distribution bends downward from observational selection effects and the model size distribution bends upward from resolution issues and mass conversation. We estimate that the Erigone parent body was D ∼ 110 km. The fossilized main belt size distribution of Erigone’s parent body prior to the family-forming event is unknown, tests with pre-shattered or rubble-pile target bodies could yield different parameters. For our purposes, however, it may not matter; it is likely that any reasonable fit between Erigone’s observed FSD and a hydrocode-produced FSD will yield a comparable amount of mass hidden below the observational detection limit (and thus the parent body’s true diameter). We find a satisfactory match between model and data for D > 4 km, with the hydrocode results scaled so the largest remnant is the same diameter as that in the Erigone family. The upward bend in the model’s FSD slope for D < 4 km is an artifact produced by the hydrocode’s resolution limit and the fact that the model must explicitly conserve mass. We estimate that the largest remnant left behind from the Erigone family-forming event has 28% of the mass of its original parent body. This suggests the Erigone parent body was D ∼ 110 km, about 20% larger than suggested by Tanga et al. (1999). While we reserve a full description of our results for a future paper, a preliminary analysis of the families identified via the Nesvorný et al. (2005) method suggests there are ∼20 families produced by the breakup of D
100 km asteroids. This value is twice that estimated by Tanga et al. (1999). The number of families produced in each size bin are given in Table 1. The majority of these families have 100 km < D < 200 km. A comparison between our parent body diameters and those described in Tanga et al. (1999) reveals agreement for several parent body diameters (e.g., Adeona, Eunomia, all cratering event-type cases), 20–50% mismatches for several more (e.g., Eos, Erigone, Flora, Koronis, Themis) and >100% mismatches for several super-catastrophic disruption cases (e.g., Dora, Gefion, Maria, Merxia). To use these families in CoEM, we need to know approximately how long they have been collecting in the main belt. Preliminary estimates of the ages of these families, determined by Yarkovsky modeling (Nesvorný et al., 2003; Vokrouhlický et al., in preparation) suggest that few are older than ∼3 Gyr. It is curious that no family is significantly older than this age; after all, the Solar System was formed 4.6 Gyr ago. Given the discussion above, we believe the most plausible solution to this enigma is that the same dynamical instability that produced the Late Heavy Bombardment ∼3.9 Gyr ago (e.g., Hartmann et al., 2000; Levison et al., 2001) also scrambled our ability to compute useful proper (a, e, i) elements beyond this epoch; in effect, this would have produced a “clean slate” in the main belt for the production of new families. Thus, we hypothesize that the largest families were produced over a time period stretching from somewhere between 3.0–3.9 Gyr ago to the present day. Given our uncertainties, we decided to split the difference and assume that the largest families formed over the last ∼3.5 Gyr ago. Note that changing this value to 4.6 Gyr ago would only introduce a ∼20% error into our estimate. 123 4.3. Constraint #3: The intact basaltic crust of (4) Vesta Asteroid (4) Vesta is one of the three largest main belt bodies (D = 529 ± 10 km; Thomas et al., 1997; Standish, 2001; see Britt et al., 2002). It is also the only known differentiated asteroid with an intact internal structure, presumably consisting of a metal core, an ultramafic mantle, and a basaltic crust (Keil, 2002). If Vesta is the ultimate source of the HED meteorites, as many believe, it differentiated and formed its crust ∼6 Myr after the formation of the first solids (i.e., CAIs) (Shukolyukov and Lugmair, 2002). Vesta’s crust is currently intact, making it unlikely that Vesta has gone through a catastrophic breakup-andreassembly episode since its crust formed (e.g., Thomas et al., 1997). Observations from the Hubble Space Telescope indicate that Vesta has a 460 km basin on its surface, which was most likely the result of an impact from a ∼35 km projectile (Marzari et al., 1996; Thomas et al., 1997; Asphaug, 1997). Collisional modeling results suggest this impact event is responsible for the so-called Vestoids, a family of D < 10 km asteroids on Vesta-like orbits whose spectral features strongly resemble eucrites and howardites (Marzari et al., 1996; Burbine et al., 2001). We use this singular crater to set limits on the frequency of large Vesta impacts in both the primordial and present-day main belt. 4.4. Constraint #4: The lunar and terrestrial impactor flux over the last 3 Gyr A fourth constraint comes from the estimated lunar and terrestrial cratering rates over the last ∼3 Gyr. Because most NEOs come from the main belt via the Yarkovsky effect (Bottke et al., 2002a, 2002b), the impactor flux on the Earth and Moon provides information on how the main belt size distribution has changed over time. For example, a steadily declining main belt population over the D < 30 km size range should produce a similarly declining NEO population/impactor flux, while a main belt population in a quasisteady state should produce a constant impactor flux. Data from crater studies indicate the lunar and terrestrial impact fluxes have been relatively constant over 0.5–0.8 to 3 Ga (e.g., Shoemaker, 1998). Proterozoic impact structures in Australia with D > 20 km formed from 0.54–2.6 Ga have a production rate of 3.8 ± 1.9 × 10−15 km−2 yr−1 (Shoemaker and Shoemaker, 1996). This value is in good agreement with a production rate of 3.7 ± 0.4 × 10−15 km−2 yr−1 reported for D > 20 km Eratosthenian craters on the Moon (0.8–3.2 Ga) (McEwen et al., 1997). Shoemaker (1998) claims that the uncertainties in these values are of the order of a factor of 2, such that it has commonly been assumed, by default, that these cratering rates have been constant for the last 3.2 Ga (e.g., Wilhelms et al., 1987). The only reported change comes from a putative factor of 2 increase in the impact flux occurring over the last 120 Ma (e.g., Grieve and Shoemaker, 1994; Neukum and Ivanov, 1994; S. Ward, personal communica- 124 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 tion). Some claim this change could have occurred over the last 400–800 Ma (McEwen et al., 1997), though this is considered controversial (Grier et al., 2001). Additional data on the lunar impactor flux over time comes from Culler et al. (2000), who dated the formation age of 155 lunar spherules found in Apollo 14 soil samples using the 40 Ar/39 Ar isochron technique. These spherules, 100–500 µm in size, are presumably droplets of lunar surface material that were melted and thrown several meters to hundreds of kilometers by an impact. If these spherules come from a variety of different craters across the Moon, their formation ages may reflect the impact history of the Moon. Culler et al. claim that the lunar impactor flux has decreased by a factor of 2–3 over the last ∼3.5 Gyr to a low about 500 to 600 Myr ago, then increased by a factor of 3.7 ± 1.2 over the last 400 Myr. This interpretation, however, is considered by some to be controversial (e.g., Hörz, 2000), with comparable tests applied to spherules from the Apollo 14, 16, 17 sites showing trends that appear to be dominated by local impact events (J. Delano, personal communication; Zellner et al., 2003). For these reasons, we believe these data do not (yet) supersede results from crater studies. The scenario most consistent with these constraints is that the main belt size distribution with D < 30 km reached a quasi-steady state several Gyr ago, with the factor of 2 increase a possible consequence of a recent large-scale disruption event (e.g., the creation of the Flora family ∼500 Ma; Nesvorný et al., 2002a). It is less consistent with main belt population models that show a factor of 3 decline over the last 3 Ga (Durda et al., 1998; Davis et al., 2002). With that said, however, it is unclear whether the accuracy claims made by Shoemaker (1998) are strong enough to rule out this possibility entirely. We leave this issue for future work. 5. Methodology In this section, we describe how CoEM accounts for the DDE as described in Sections 1 and 2. This leads into the selection of our initial conditions. 5.1. Accounting for the putative dynamical depletion event in CoEM To briefly recap our discussion from Sections 1 and 2, many modern models of main belt evolution claim the main belt region once held several Earth-masses of material, enough to accrete large planetesimals and possibly planetary embryos via runaway growth. While planetary embryos presumably would have agglomerated most of this mass, planetesimals smaller than the Asteroid Ceres would still have contained hundreds of times more mass than the current main belt. The observed main belt population, however, contains relatively little mass. To account for this difference, many models now assume that a DDE, possibly triggered by the formation of Jupiter, scattered numerous planetesimals out of the main belt region and drove the main belt to its current state. Our problem is to model this scenario in CoEM while avoiding unknowns like the size of the initial main belt population and the nature/timing of the DDE. To account for the nature of the DDE, we assumed in Section 3.2 that collisions among diameter D < 1000 km planetesimals in the main belt zone over the last 4.6 Gyr have been dominated by the same collision probabilities and impact velocities found there today. This approximation simplifies the nature of the DDE by making Pi and Vimp constants. This leaves us with two unknowns: the mass of the initial population and the length of time that the primordial main belt population experienced comminution. To treat these variables, we invoke the following approximation. Lets assume our simulations start with an massive mass belt whose D < 1000 km members have been dynamically excited enough to depart the accretion phase and enter the fragmentation phase. Assuming Pi and Vimp are constants, we can define the collisional lifetime τ for planetesimals between D and D + dD as: 1 Pi = τ 4 DT (DT + d ′ )2 (Nrem + Ndep ) dd ′ , (10) ddisrupt where DT is the representative target body diameter in the size interval and d ′ is the projectile diameter. Here Nrem and Ndep are defined as two parts of the same initial main belt size distribution: Nrem is the remnant of the initial population that stays behind in the main belt zone, while Ndep is the population that is eventually ejected from the main belt zone via dynamical processes (e.g., sweeping resonances, gravitational interactions with planetary embryos and Jupiter). Asteroids in Nrem can be struck by bodies in Nrem and Ndep , and vice versa. We assume the ratio of Ndep over Nrem at time t = 0 Gyr over all size bins is f . Since both populations uniformly occupy the same volume of space, f will remain nearly constant as t increases (i.e., they would be precisely identical if Poisson statistics did not govern the infrequent breakups of large asteroids). This “self-similarity” concept is just another way of saying that one can divide a single size distribution into smaller versions of itself at any time, with each subpopulation sharing the overall shape. If we now substitute Ndep = f Nrem into Eq. (10), we find that τ becomes f + 1 times shorter than τ where Ndep = 0. Hence, by placing more objects in Ndep , we produce faster collisional evolution among both sub-populations. Conversely, Eq. (10) indicates that when Ndep = 0, Nrem can reach the same evolutionary state if time is increased by a factor f + 1. This means there is a direct tradeoff between the size of the initial main belt population and evolution time. Using this concept, we compensate for the unknown size and departure time of Ndep by extending our evolution time beyond 4.6 Gyr. In the process, the model time in CoEM becomes a “pseudo-time” that must be interpreted in a different The fossilized main belt size distribution way than “real” time. For example, in a real time scenario, the primordial main belt would be Nrem + Ndep just prior to the onset of fragmentation. After undergoing comminution for time t1 , the main belt experiences a DDE and loses Ndep . The Nrem population then undergoes collisional evolution for an additional time t2 = 4.6 Gyr − t1 . In our pseudotime scenario, however, the collisional evolution of Nrem is tracked until it matches our model constraints or until it becomes clear that the trial has failed. If we find a good match with our constraints for values of pseudo-time longer than 4.6 Gyr, it tells us that the primordial main belt must have originally contained more mass than Nrem (i.e., Ndep > 0). We use this example case to further clarify our procedure. Lets assume we performed a CoEM trial where Nrem undergoes comminution for a pseudo-time longer than 4.6 Gyr. If we get our first match to our model constraints at tpseudo = 20 Gyr, it means the main belt required more than 4 times the degree of comminution (20 Gyr/4.6 Gyr) than could have taken place if Nrem were tracked by itself over Solar System history. To account for this extra comminution in real time, we can only infer that Ndep > 0 and that a DDE event took place. The strengths of this technique are that we can ignore the values of t1 and Ndep in our model while allowing CoEM to investigate possible solutions for the Q∗D disruption function and the initial size and shape of Nrem . It also yields information on how much collisional evolution could have ever taken place in the main belt. The weakness is that we cannot solve for t1 and Ndep , which describe the nature and timing of the DDE event. Thus, we cannot tell if Ndep was massive and short-lived or less massive and longer-lived. To do this, we will have to take our Q∗D and Nrem solutions and apply them to a real time model where we include the effects of a realistic DDE. This will the subject of a separate paper in the near future. 5.2. Initial main belt size distribution Using the method described in Section 6.1, we can now explore the nature of the initial size-frequency distribution for Nrem as well as Q∗D . As described in Section 2, this can be challenging work, with many combinations of these two functions capable of reproducing the observed main belt. Our task is to filter out the true pair from the impostors. One powerful method used by Durda et al. (1998) was to fix the shape of the initial population and then use a CoEMSM-like code (Section 3.4) to determine whether particular Q∗D functions could reproduce the observed size distribution. Note that in their study, they assumed the Q∗D function had to have a hyperbolic shape and it had to match values found in laboratory shot experiments. Using initial populations with starting masses 3–6 times the current main belt’s mass, size distributions containing more bodies in every size bin than the current main belt, and no DDE (see Section 2), Durda et al. found that collisions could not reproduce the shape of the main belt size distribution unless the starting population con- 125 tained barely more D
100 km asteroids than the observed population, regardless of the shape of the Q∗D function. Our interpretation of their results is that the D
100 km part of the main belt size distribution must be a by-product of accretion rather than collisional evolution. Despite the claim that Durda et al. (1998) can reproduce the main belt size distribution, we find their best-fit solution for Nrem and Q∗D to be unrealistic. Our computations indicate the Durda et al. (1998) best-fit case yields 44 D > 100 km asteroid disruptions over the last 3.5 Gyr, a factor of 2 higher than our estimates (Section 4.2). Moreover, the shape of their derived Q∗D function is highly discordant when compared to results from hydrocode modeling (e.g., Benz and Asphaug, 1999). The odd shape is needed to eliminate the large excess of main belt material. Finally, the evolution of their initial population would produce a steadily declining NEO population that is likely discordant with the constant impact flux observed on the lunar maria for the last ∼3 Gyr (see Section 1). These problems suggest there may be a more satisfying combination of Nrem and Q∗D than the one deduced by Durda et al. (1998). In order to find that combination, we used CoEM-SM to test a multitude of shapes for the initial main belt size distribution ala the procedure described by Durda et al. (1998). In each case, CoEM-SM searched parameter space for the best possible Q∗D function that would allow Nrem to reproduce the observed population. Our initial results showed good agreement with the conclusions of Durda et al.; Nrem for D
100 km asteroids are likely to mimic the observed population, and breakup events among these bodies do not occur frequently enough to significant modify the shape of the population in this size range. To test the shape of Nrem for D  100 km asteroids, we selected populations both smaller and larger than the observed population. In some cases, we even experimented with size distributions that had multiple elbows. After numerous runs, our results can be summarized as follows. For Nrem (D  100 km) populations set higher than the observed one, our results repeat those of Durda et al. (1998); either we could not reproduce the shape of the observed main belt or we solved for Q∗D functions that created too many families and were highly dissimilar to those reported elsewhere (e.g., Asphaug et al., 2002; Holsapple et al., 2002). Particular problems were found when we used initial Nrem populations with more D = 20–100 km bodies than the observed main belt. In these cases, CoEM-SM had difficulty finding a Q∗D function that could eliminate numerous D = 20–100 km bodies without also decimating the D > 100 km population. Numerous large breakup events lead to too many families, too many fragments, and in some cases, too many bodies capable of creating D ∼ 460 km craters on Vesta. When Nrem (D  100 km) was set lower than the observed main belt, however, our results produced Q∗D functions similar to those described by hydrocode modeling (Benz and Asphaug, 1999). If we assume the number of 126 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 D = 20–100 km bodies does not exceed the observed population, Nrem (D  100 km) is restricted to power-law slopes that are fairly flat. Note that modest variations in the slope do not affect our results because the number of D = 1–10 km bodies in the initial population are 2–3 orders of magnitude lower than those observed. Interestingly, our findings are analogous to the estimated initial conditions predicted by Anders (1965) (see Section 2). They are also comparable to accretion code results where fragmentation has been turned off (e.g., Weidenschilling et al., 1997).5 For example, asteroid belt accretion simulations without fragmentation presented in Wetherill (1989); (see Wetherill and Stewart, 1989 for additional details) show the slope of the cumulative planetesimal size distribution between D ∼ 50–2000 km follows a cumulative power law index between −4.0 and −4.5. This value is similar to the observed slope of those main belt asteroids with D > 100 km (Fig. 1). For D  50 km, the accretion size distribution goes through a change in slope and becomes shallow (cumulative power law index near 1.5). The shallow slope is a byproduct of planetary embryos agglomerating nearly all of the small bodies within their feeding zone. While CoEM-SM provides us with many useful results, we do not use it for our production runs because it does not properly account for the effects of stochastic breakups. Instead, we used insights from these results to generate our initial Nrem populations for the more realistic code CoEMST (Fig. 4). Because very few asteroid families come from the disruption of D > 200 km asteroids, our initial population for D > 200 km uses the observed asteroids, with a few objects added in, to account for the original Eos and Themis parent bodies. For D < 200 km objects, we follow an incremental power law index of −4.5, producing slightly more objects than the observed main belt size distribution. This slope continues until it reaches D = Dx . We treat the location of this size distribution elbow as an unknown and test values Dx = 80, 100, and 120 km. As described above and below, tests with Dx < 80 km are probably not warranted. Finally, for D < Dx , we gave Nrem an incremental slope of −1.2 to limit the initial population of D = 50–100 km asteroids. As described above, modest variations to this value should do affect our results. (As an aside, we note that all these values could be tweaked a bit in order to reach the true initial Nrem , but given available information and our computational limitations, we believe they are a reasonable place to start.) 5 Note that great care must be taken when comparing accretion model results with fragmentation turned on (e.g., Greenberg et al., 1978; Wetherill and Stewart, 1989, which included cases with and without fragmentation; Stern, 1996; Wetherill and Inaba, 2000; Kenyon, 2002; Inaba et al., 2003; Kenyon and Bromley, 2004) to our results because the initial populations, Q∗D functions, and fragmentation laws used by these codes have yet to be tested against the constraints described in Section 4. Fig. 4. Our initial main belt size distribution, which is based on a combination of the observed main belt population and accretion code results. We set the number of D > 200 km asteroids close to those in the observed main belt, assuming that a limited number of these objects ever disrupt. For D < 200 km objects, we follow an incremental power law index of −4.5 until reaching the transition point Dx . We tested Dx = 80, 100, and 120 km in our code. For D < Dx , the size distribution is given a shallow slope (−1.2). Modest changes to this value do not affect our results. Note that this primordial population is only a mathematical convenience; we believe the real population was probably hundreds of times larger than this, with most of the mass eliminated by dynamical processes associated with the formation of Jupiter (e.g., Petit et al., 2002). 6. Model runs In this section, we use two sample CoEM-ST trial cases to demonstrate our methods (Sections 6.1, 6.2) before presenting results from our production runs (Section 6.3). Our best fit results are then applied to the constraints provided by Asteroid (4) Vesta (Section 6.4). 6.1. Demonstration case #1: A good match with constraints To illustrate how CoEM-ST works and how we analyze our results, we describe here a sample trial for a single random seed (Fig. 5). In this case, we use CoEM-ST with the input size distribution described in Section 6.1, where the size distribution changes slope at Dx = 120 km. For our disruption scaling law (Q∗D ), the Eq. (5) parameters are E = 0.861, F = −0.913, G = −0.502, and H = −0.308 (see Fig. 9, where this curve labeled Run 15). Figure 5 shows six snapshots from the evolution of our size distribution, which is tracked for a pseudo-time of 50 Gyr. At time tpseudo = 1.0 Gyr, we see the size distribution has already developed a bump near D ∼ 2 km. The bump is produced by the change in Q∗D slope between the strength and gravity regimes near D ∼ 0.2 km (see Section 2). As our model time increases, Fig. 5 shows a better and better fit to the observed main belt population. To determine the goodness of fit between the binned main belt data and our model, we originally used a χ 2 test (Press et al., 1989). The fossilized main belt size distribution 127 Fig. 6. Tracking the goodness-of-fit for the test case shown in Fig. 5. The 2 , our metric for determining a match between the line represents ψSFD 2 model and observed size distributions. The dotted line indicates ψSFD < 20, our estimated value for a positive match. The open squares represent where 2 the probability associated with χFAM is greater than 30%. This value measures the likelihood that our model population produces the same number of asteroids families as those described in Table 1. We see that the combi2 2 nation of ψSFD and χFAM are satisfied for many values between tpseudo = 9.25–17.25 Gyr. 2 : tribution. We call this more subjective criterion ψSFD  NMODEL (D) − NMB (D) 2 2 = . ψSFD (11) 0.2NMB (D) D Fig. 5. Six snapshots from a representative run where we track the collisional evolution of the main belt size distribution for a pseudo-time of 50 Gyr. This run uses a starting population with Dx = 120 km and a Q∗D function associated with Run 15 (Fig. 9). The bump near D ∼ 120 km is a leftover from accretion, while the bump at smaller sizes is driven by the transition at D ∼ 200 m between strength and gravity-scaling regimes in Q∗D . Our model main belt achieves the same approximate shape as the observed population at tpseudo = 9.25 Gyr. The model closely adheres to the observed population for many Gyr after this time. Eventually, comminution eliminates enough D > 100 km bodies that the model diverges from the observed population. Several attempts using these methods, however, did not lead to satisfying results. Our main problem is that we are trying to fit our model results to an observed population that extends over 1 km < D < 1000 km, with the incremental number of bodies in each size bin varying from >105 on the small end to 1–3 bodies on the large end. We found that standard χ 2 tests were not diagnostic of visually satisfying fits, partly because the observational errors in our size bins are not 1σ (see Section 4.1) but also because our size distribution is logarithmic in nature. Instead, we decided to use a hybrid- or quasi-χ 2 test which, while not following Poisson statistics, does provide us with a useful metric describing when our model provides a reasonable visual fit across the entire span of our size dis- We assume that our model is a good fit if lies within 20% of the observed data across all bins. Note that the 20% value was determined experimentally via numerous comparisons between our model results and data. For reference, this range is slightly smaller than the dot size used to plot the observed main belt size distribution in Fig. 5. Tests indi2 < 20 generally provides a good match becate that ψSFD tween model and data. The 6 snapshots shown in Fig. 5 2 values of 390, 356, 230, 76, 9.0, and 55 for have ψSFD tpseudo = 0, 1.0, 2.0, 5.0, 9.25, and 30.0 Gyr, respectively. 2 Figure 6 shows a plot of ψSFD for all timesteps. The small jumps are byproducts of catastrophic breakup events among large bodies, which flood the main belt with fresh fragments until they are beaten back by collisional evolution. At tpseudo = 9.25 Gyr, our model population reaches 2 < 20 for the first time. This value of tpseudo means that ψSFD in this simulation, Nrem cannot reach the same shape as the observed size distribution over Solar System history unless the main belt once held more mass (i.e., Ndep > 0; see Section 5.1). In addition, we find the main belt size distribution enters into a quasi-steady state close to the observed values after tpseudo = 9.25 Gyr. An apparent balance is achieved between disruption events grinding the small body population down and larger-scale breakup events building the small body population back up. We believe this explains why the lunar and terrestrial crater record shows evidence for a con- 128 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 stant impactor flux over the last 3 Gyr (see Section 4.4). Note that large-scale breakup events do cause perturbations to the small body population, which can cause the population to fluctuate a bit over time. If the breakup event is big enough (e.g., the disruption of Vesta, Pallas, or Ceres), the entire population is modified to such an extent that it will fail the 2 test for an extended time. ψSFD It is plausible that even more comminution took place than this nominal amount. Figure 6 indicates that our model 2 < 20 values between tpseudo = size distribution has ψSFD 9.25–20.5 Gyr. Over this interval, stochastic breakups among large asteroids kept our model population replenished with fragments, enough to allow it to closely adhere to the constraints of the observed population. Figure 5 (last frame) indicates this period only ends when the reservoir of D > 2 envelope. These 80 km bodies falls outside the ψSFD timescales imply we may not be able to tell whether the main belt has only recently entered its current collisional state or if it has been there for an extended time. We know, however, that the lunar impactor flux has been relatively constant over the last 3 Gyr and that most NEOs come from the main belt (Bottke et al., 2002a, 2002b). We infer from this information that the main belt population has not changed very much over the last 3 Gyr. It is necessary, but not sufficient, for our model size dis2 < 20. We also need to determine tribution to reach ψSFD whether our model can match the constraints provided by asteroids families. We do this by having CoEM keep a running Fig. 7. Four snapshots from a second representative run where we track the collisional evolution of the main belt size distribution for a pseudo-time of 50 Gyr. Here we use Dx = 120 km and a Q∗D function associated with 2 < 20 constraint is Run 3 (Fig. 4). The only acceptable match to our ψSFD found at tpseudo = 5 Gyr. The remaining snapshots show timesteps where the observed population was not reproduced. tally of the total number of destroyed bodies produced every timestep. Recall that in Section 4.2, we argued that most observed families formed over the last 3.5 Gyr of Solar System history; families created prior to this time were likely dispersed by large scale dynamical processes associated with the DDE and possibly the Late Heavy Bombardment. Thus, when a CoEM trial reaches t > 3.5 Gyr, we compute the change in the number of destroyed bodies Ndisrupt (D) over a time window of 3.5 Gyr. This value is compared to the observed number of families in the 100 km < D < 400 km size bins (Table 1) using a standard χ 2 test (Press et al., 1989). We assume that the error in the number of known families in each bin is normally distributed, and that the value obtained, 2 , must be better than 1σ (i.e., probability >30%) for χFAM the trial to be considered a positive match. 2 is conThe open squares in Fig. 6 show where χFAM sidered a good fit to our family data. For this run, we find matches for tpseudo = 9.25–13 Gyr and a few near 17 Gyr. If we combine the observed main belt size distribution and 2 2 and χFAM are family constraints together, we find that ψSFD satisfied for 4.25 Gyr. These values will be used in Section 6.3 to compare different sets of runs to one another. 6.2. Demonstration case #2: A bad match with constraints To compare and contrast these results, we show a less successful trial that uses a Q∗D function with E = 0.923, F = −0.659, G = −0.410, and H = −0.246; see Eq. (5). This function is labeled as Run 3 on Fig. 9. In this case, the slope of Q∗D in the gravity regime is shallower than in the previous case, allowing large asteroids to disrupt more easily. Snap2 and shots from this trial are shown in Fig. 7, while the ψSFD 2 2 χFAM results are shown in Fig. 8. Overall, ψSFD < 20 is met 2 once at tpseudo = 5 Gyr while the χFAM criterion is met near tpseudo = 17 and 24 Gyr. The combination of the two produced no matches. An investigation of this trial provides some interesting in2 < 20), we found the sights. For tpseudo = 5 Gyr (where ψSFD number of model families produced by D > 100 km disruption events over a 3.5 Gyr window was slightly more 2 fit at that than half the observed number. The good ψSFD time appears to be a fluke produced by good timing and the stochastic breakup of a D ∼ 250 km asteroid. Timesteps between 5–15 Gyr produce 2–5 times the observed number of families in the D ∼ 120 km bin. This explains the excess number of bodies in the main belt size distribution seen over 2 values (Fig. 8). Reathis time (Fig. 7) and their high ψSFD 2 sonable values for χFAM do not return until the D > 100 km population has become highly depleted for tpsuedo > 15 Gyr. At this point, however, there are so few large asteroids re2 is impossible. maining that matching ψSFD 6.3. Overall results Using the methods described above, we are now ready to explore which combinations of initial main belt size dis- The fossilized main belt size distribution 129 Fig. 8. Tracking the goodness-of-fit for the test case shown in Fig. 7. See 2 Fig. 6 for plot details. In this case, there are no examples where ψSFD and 2 2 χFAM metrics are met simultaneously. The positive matches with χFAM occur after the D > 100 km population has been highly depleted by comminution. tributions and Q∗D functions produce the best fit to our constraints. We begin with the initial size distribution with Dx = 120 km (Fig. 4). The Q∗D functions tested by CoEMST are plotted in Fig. 9. For reference, the Q∗D function described by Benz and Asphaug (1999) for basaltic targets being disrupted by projectiles at V = 3 km s−1 is labeled as Run 13. Each combination was tested using at least 100 separate trials of CoEM using different random seeds, with the maximum pseudo-time set to 50 Gyr. We define this as a “run.” Our results indicate this time is sufficiently long for most individual trials to achieve at least a few posi2 2 and χFAM before comminution drags tive matches in ψSFD the model size distribution irrevocably away from the observed main belt. For the Q∗D functions producing the best matches to our constraints, we ran an additional 100–200 test cases. The results of our runs are summarized in Table 2. Column 1 is the run index. Columns 2 and 3 describe successes and failures, where runs are declared a “success” if they 2 2 < 20 and χFAM better than 1σ (i.e., probabilmatch ψSFD ity >30%) at least once over 50 Gyr. Runs 8–22 all have success rates greater than 70%. Columns 4 and 5 describe the mean and median pseudo-times that each set of runs had 2 2 and χFAM . Except for the their first positive match with ψSFD first few runs, where successes are infrequent and/or short, we see a gradual increase in these values as the slope for Q∗D in the gravity regime steepens. Recall that the steeper the slope, the more difficult it is to break up large asteroids, which in turn means fewer and fewer fragments are available to enhance the small body population. Columns 6 and 7 are the mean and median of the cumulative time that each of the 100 test cases from each run maintains a positive match 2 2 < 20 and χFAM . These values give us some sense with ψSFD of whether the match between model and data is a likely occurrence or a fluke. Column 8 describes the cumulative Fig. 9. The Q∗D functions tested in this paper. The numbers to the right of each curve correspond to their run number, with a run is defined as 100 trials of CoEM started with different random seeds (see Table 2). The bottommost curve is Run 1 and the topmost curve is Run 46. For each run, the size distribution and number of breakup events in each size bin for 100 trials were tracked for 50 Gyr. The colors represent the total time in each 2 2 < 20 and χFAM ; see Table 2, run that the success criteria was met (ψSFD column 8). For reference, the Q∗D function described by Benz and Asphaug (1999) for V = 3 km s−1 impacts into basalt is given by Run 13. The curves representing are best-fit cases are shown in red. The transition point between strength- and gravity-scaling is D ∼ 200 m. All of the plotted functions pass through the normalization point at D = 8 cm and Q∗D = 1.5 × 107 erg g−1 . 2 Fig. 10. A histogram showing the times our goodness-of-fit criteria ψSFD 2 and χFAM were met for 200 test cases of Run 15 (see Table 2 and Fig. 5). The mean for the histogram is tpseudo = 11.8 Gyr, with peaks at 7.375 and 8.625 Gyr. 130 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 Table 2 A statistical description of our results for tests that track the evolution of our initial size distribution with Dx = 120 km (Fig. 4) Run % good 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 6 10 23 37 45 64 69 79 77 81 78 78 83 79 84 79 83 78 80 77 77 72 63 66 56 60 59 56 47 49 33 47 29 35 35 28 29 23 17 22 10 10 12 0 0 0 % bad 94 90 77 63 55 36 31 21 23 19 22 22 17 21 16 21 17 22 20 23 23 28 37 34 44 40 41 44 53 51 67 53 71 65 65 72 71 77 83 78 90 90 88 100 100 100 tfirst 4.5 ± 0.5 4.6 ± 1.0 5.0 ± 1.4 5.0 ± 1.4 5.7 ± 2.1 5.8 ± 2.3 7.0 ± 3.1 6.6 ± 3.3 7.3 ± 3.9 8.2 ± 4.2 8.8 ± 5.0 9.2 ± 5.0 9.7 ± 5.3 9.6 ± 4.7 9.5 ± 5.3 10.1 ± 5.6 10.5 ± 5.5 11.3 ± 6.3 11.6 ± 6.1 11.2 ± 5.8 12.4 ± 7.3 10.7 ± 5.6 11.7 ± 5.2 14.5 ± 8.3 12.3 ± 4.9 13.7 ± 6.3 12.5 ± 7.5 14.4 ± 6.8 14.8 ± 6.9 15.9 ± 7.3 15.7 ± 7.2 15.5 ± 7.1 21.0 ± 8.7 19.6 ± 9.6 18.4 ± 9.8 18.8 ± 10.1 20.8 ± 10.1 17.6 ± 6.9 17.5 ± 8.4 21.3 ± 8.1 18.3 ± 8.8 19.3 ± 10.2 18.2 ± 4.4 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 Median tfirst 4.8 4.0 4.5 4.8 5.0 5.0 6.0 5.2 6.2 7.0 7.5 7.8 8.5 8.5 8.0 8.0 9.5 9.5 10.2 9.8 10.8 9.5 10.8 13.0 12.0 13.0 9.8 13.5 12.0 14.8 14.2 13.5 20.0 15.2 15.5 18.0 20.2 16.5 15.0 21.5 16.5 18.5 18.8 0.0 0.0 0.0 tfit 0.2 ± 0.0 0.4 ± 0.2 0.8 ± 0.8 0.8 ± 0.7 1.0 ± 1.0 1.3 ± 1.0 1.3 ± 1.0 1.8 ± 1.2 1.9 ± 1.4 1.9 ± 1.4 2.1 ± 1.4 2.1 ± 1.4 2.0 ± 1.4 2.1 ± 1.4 2.2 ± 1.5 1.9 ± 1.3 2.2 ± 1.4 2.1 ± 1.6 2.0 ± 1.4 2.0 ± 1.5 1.9 ± 1.2 2.0 ± 1.3 2.2 ± 1.6 1.4 ± 1.2 1.9 ± 1.2 1.6 ± 1.0 1.7 ± 1.1 1.4 ± 1.1 1.4 ± 1.2 1.4 ± 1.1 1.1 ± 0.8 1.0 ± 0.9 1.2 ± 0.8 1.2 ± 0.8 1.5 ± 1.2 1.0 ± 0.8 1.1 ± 0.9 1.1 ± 0.8 1.1 ± 0.7 0.9 ± 0.8 0.7 ± 0.6 0.9 ± 0.6 0.9 ± 0.6 0.0 ± 0.0 0.0 ± 0.0 0.0 ± 0.0 Median tfit 0.2 0.2 0.5 0.5 0.8 1.0 1.0 1.5 1.5 1.8 1.8 2.0 1.8 1.8 1.8 1.5 2.0 1.8 2.0 1.8 1.5 1.5 1.5 1.0 1.8 1.5 1.5 1.0 1.2 1.0 0.8 0.8 1.0 1.0 1.0 0.8 1.0 1.0 1.0 0.5 0.5 1.0 0.8 0.0 0.0 0.0  tfit 1 4 18 30 45 85 90 138 144 153 162 167 169 166 181 147 180 161 159 155 143 140 135 90 107 94 101 77 65 70 34 49 35 43 50 28 32 26 19 19 6 9 10 0 0 0 Every entry is based on at least 100 trials (i.e., 100 CoEM-ST test cases using different random seeds) that were tracked over tpseudo = 50 Gyr; this defines a “run.” Runs 11–18 are based on 200–300 trials. Column 1 is the run index. The slope of Q∗D increases with the run index. Columns 2 and 3 describe successes 2 2 and failures, where our success criteria is ψSFD < 20 and χFAM better than 1σ (i.e., probability >30%) at least once during their 50 Gyr of evolution. Columns 4 and 5 describe the mean and median pseudo-times that each set of runs has their first positive match with our success criteria. Columns 6 and 7 describe the mean and median of the cumulative pseudo-times that test cases from each run maintains a positive match with our success criteria. Column 8 describes the cumulative pseudo-time over all test cases (in Gyr) where our success criteria are satisfied. For Runs 11–18, these values have been averaged to correspond to 100 test cases. 2 pseudo-time over all test cases (in Gyr) where ψSFD and 2 χFAM are satisfied. We use these results to color-code the Q∗D curves from Fig. 9. The runs with cumulative pseudotimes >160 Gyr are plotted in red. The highest values are found for Runs 11–18, with the exception of Run 16 (red curves in Fig. 9). We believe Run 16 is likely a statistical fluke. Because the highest cumulative success times are found for Runs 11–18, we conclude that these Q∗D functions give us our best opportunity of reproducing the observed main belt over the age of the Solar System. For reference, the Eq. (5) parameters for Run 11 are E = 0.871, F = −0.877, G = −0.489, and H = −0.299, while those for Run 18 are E = 0.854, F = −0.940, G = −0.511, and H = −0.314 The fossilized main belt size distribution 131 A visual inspection of trial cases in these runs indicates why these initial Nrem populations are unlikely to reproduce our constraints. The Dx = 80 and 100 km cases create a bump near D ∼ 80–100 km that is difficult to grind away using Q∗D functions with linear slopes in the gravity regime. Though decreasing the gravity regime slope makes it easier to disrupt D > 100 km asteroids, it also tends to produce too many families and too many fragments, the latter which can drive the model size distribution away from the observed main belt. Thus, we conclude that our initial size distribution with Dx = 120 km is more likely to yield the observed main belt than those with Dx = 80 or 100 km. 6.4. Constraints from (4) Vesta Fig. 11. A histogram showing the times when our goodness-of-fit criteria 2 2 ψSFD and χFAM were met for the first time for the test cases from Run 15 (see Table 2 and Fig. 5). The distribution peaks at tpseudo = 6.5 Gyr, with tfirst having mean and median values of 9.5 ± 5.3 Gyr and 8.0 Gyr, respectively. The paucity of values with tpseudo < 4.6 Gyr indicates the main belt had to undergo more collisional evolution over its history than could have been produced by simply running the observed population backward in time 4.6 Gyr. (see also Table 2 and Fig. 9). Figure 10 shows a plot of the cumulative time for all Run 15 test cases where our model constraints are met. Note that the Q∗D function predicted by the hydrocode modeling results of Benz and Asphaug (1999) are represented by Run 13, right in the middle of our best fit results. We take this as confirmation of the validity of our approach and the general accuracy of the results. The median first success pseudo-times for Runs 11–18 are 7.5–9.5 Gyr. A plot showing the distribution of first success pseudo-times for Run 15 is shown in Fig. 11; its mean and median are 9.5 ± 5.3 and 8 Gyr, respectively. Thus, to reach its current state, the main belt needed roughly twice the amount of collisional evolution that a nominal Nrem would have experienced over 4.6 Gyr. This implies that Ndep > 0 and the primordial main belt experienced an intense early phase of collisional evolution. We believe that much of the main belt’s wavy-shaped size-frequency distribution was created during this time, such that its shape can be considered a “fossil” remnant from this violent epoch. Using the same procedure, we have also investigated initial size distributions from Section 6.1 that have Dx = 80 and 100 km (Fig. 4). To cover this parameter space in a reasonable amount of computation time, we set our maximum tpseudo value to 25 Gyr, half the value used above, and we only generated 50 test cases per run. We also selected a somewhat broader range of Q∗D functions than shown in Fig. 9. Interestingly, none of the Dx = 80 and 100 km tests produced results nearly as successful as our Dx = 120 km runs, with no successes for the former and 1–3% successes for the latter. Because these values are so much lower than those obtained by Dx = 120 km runs using the same CoEMST parameters, we believe the successful matches in the Dx = 100 km runs are predominately flukes. Using our best-fit combination of Nrem and Q∗D , we can now address the constraints provided by the intact basaltic crust of (4) Vesta. As described in Section 4.3, Vesta is a D = 529 ± 10 km asteroid with an intact basaltic crust (e.g., Thomas et al., 1997; Standish, 2001; see Britt et al., 2002). It also possesses a 460 km diameter crater, 13 km deep, that completely dominates one hemisphere (Thomas et al., 1997). From Vesta’s morphology, we can infer that this crater was produced after Vesta formed its crust. Thus, this crater provides a more restrictive constraint on the main belt’s collisional history than Vesta’s intact basaltic crust; instead of worrying about the obliteration of the crust, we instead focus on the likelihood that Vesta experienced one (and only one) such event over Solar System history. Using the results from Run 15 and Dx = 120 km described in Section 6.3, we can readily check this scenario. We estimate that Vesta’s intrinsic collision probability is Pi = 2.8 × 10−18 km−2 yr−1 , nearly the same as that assumed for typical main belt asteroids (e.g., Farinella and Davis, 1992; Bottke et al., 1994). The projectile that created the crater was D ∼ 35 km (Marzari et al., 1996; Thomas et al., 1997; Asphaug, 1997). To get the approximate number of available D ∼ 35 km bodies over 4.6 Gyr of main belt history, we interpolate between the central values of bins D = 31 and 39 km for each timestep in our model size distribution. Thus, the average interval between impacts is given by: 1 τimpact = Pi (DVesta + dproj )2 4 × N (31 km  D  40 km). (12) To simulate the timing and likelihood of such large impacts on Vesta, we integrated Eq. (12) in our Run 15 model. We assumed the impacts followed Poisson statistics. We made this computation every 0.25 Gyr in each of Run 15’s 100 test cases, and we reran our results 10 times using different random seeds. The probability distribution obtained from this model indicates that the median number of D ∼ 35 km objects striking Vesta over tpseudo intervals of 0–12, 13–27, and 28–43 Gyr was 0, 1, and 2, while the mean number for 132 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 tpseudo  10, 24, and 41 Gyr was 0.5, 1.5, and 2.5, respectively. These values agree with more simplistic estimates from the current main belt, where ∼280 objects between 31 km < D < 39 km yields an average impact interval of tpseudo = 18 Gyr. Though one must be careful with small number statistics, the singular nature of Vesta’s crater suggests the main belt did not undergo the equivalent of more than ∼27 Gyr of collisional evolution; if it had, the odds are that Vesta would have been struck by at least two D ∼ 35 km impactors. These results are consistent with our predictions from Section 6.3 that the main belt experienced tpseudo ∼7.5– 9.5 Gyr of collisional evolution (Table 2). Assuming these values are reasonable, we estimate the probability that a single 31 km < D < 39 km would have struck Vesta over Solar System history was roughly 42–53%, while the probability one would have struck over the last 3.5 Gyr of real time and produced the Vesta family is ∼19%. 7. Discussion and implications Up to this point, we have only discussed the implications of our results for the asteroid belt in general terms. In this section, we use our results to address number of issues dealing with the asteroid disruption rates, the size of the primordial main belt, and the spin rate distribution of the largest asteroids. We also compare our Q∗D estimates to recent hydrocode modeling work. Fig. 12. The mean disruption rate of asteroids in each bin per Gyr for the 200 test cases of Run 15. The rates were computed by summing the total number of disruption events across the test cases for each timestep, dividing by the total number of runs, and then subtracting these values from one another over 1 Gyr intervals. The mean interval for Karin-size breakups (25 km < D < 35 km) is 15–30 Myr, while those for D > 100 km objects is 0.25 Gyr. Nesvorný and Bottke, 2004) occur ∼30–60 times per Gyr, with a mean interval between breakups of 15–30 Myr. 7.1. Frequency of asteroid disruption 7.2. Estimating the size of the primordial main belt Using Run 15, we have computed the average number of catastrophically disrupted bodies produced by the Nrem population over time. The average number of D > 100 km breakups after tpseudo = 1, 2, 5, 8, 9, and 10 Gyr across the 100 trials of Run 15 are 3, 6, 17, 28, 31, and 35, respectively. Thus, if the main belt has experienced the equivalent of 7.5–9.5 Gyr of comminution over the last 4.6 Gyr, as suggested by our best-fit model results, it appears that 28–35 D > 100 km asteroids in the Nrem population should have disrupted. Note that this estimate excludes the disruption and fragmentation products produced by the Ndep population. For this reason, our CoEM values should be used carefully. The mean disruption rate across our 100 test cases per Gyr is shown in Fig. 12. The rates were computed by summing the total number of disruption events across our test cases for each timestep, dividing by the total number of runs (100), and then subtracting these values over 1 Gyr intervals. We find that for main belt populations that have undergone tpseudo = 10 Gyr of collisional evolution, roughly ∼4 asteroids with D > 100 km disrupt every Gyr. Not surprisingly, this value is consistent with our constraints from Section 4; recall that a metric of CoEM’s success was the disruption rate of observed families. Karin-size breakups (25 km < D < 35 km; Nesvorný et al., 2002c; Our results from Section 6 indicate that a significant amount of collisional evolution occurred while the primordial main belt is massive, and that much of this mass had to be removed dynamically. Though a quantitative study of the collisional and dynamical evolution of the main belt will be investigated in a future paper, we can use our results, together with a simple back-of-the envelope calculation, to crudely estimate the primordial size of the main belt just prior to the onset of fragmentation. The number of breakups events that ever occurred in the main belt can be written as T /τ , where τ is the collisional lifetime of an asteroid with diameter D and T is the time needed for the input size distribution to attain a good match 2 2 and χFAM . If we assume this quantity is a conwith ψSFD stant, and that the main belt was once massive, we get: TPrim TNow T = + , τ τPrim τNow (13) where TPrim and τPrim are the evolution time and collisional lifetime of objects in the primordial main belt before the DDE and TNow and τNow are the values that have existed in the main belt over the last ∼4.5 Gyr. We assume for now that the only main belt dynamical excitation event occurred for t > 4.5 Gyr. By substituting Eq. (12) for τ , we can rewrite The fossilized main belt size distribution the equation as: Nmb T = (xNmb )(TPrim ) + Nmb TNow . (14) Here, we use Nmb as a stand-in for the approximate size of the observed main belt. We set x to be a factor describing the size of the primordial main belt. We assume that the observed main belt has been near its current size for TNow = 4.5 Gyr. Given that our model results predict T = 7.5–9.5 Gyr of collisional evolution in the main belt, we can solve for x provided we know TPrim . Here is where things get tricky. Recall that TPrim represents not only the time the main belt was massive but also the time when asteroids ejected by the DDE were crashing into those that stayed behind in the main belt. For this computation, we decided, somewhat arbitrarily, to set TPrim = 20 Myr; 10 Myr for Jupiter to accrete its gas (Pollack et al., 1996), and another ∼10 Myr for significant numbers of dynamically-excited asteroids to collide with the leftovers in the asteroid belt (e.g., Petit et al., 2001). If these values are reasonable, the primordial main belt for D  1000 km bodies was roughly 150–250 times the size of the current main belt. These values are similar to predictions of DDE models that suggest that the main belt may have lost ∼99.5% of the bodies in its original population (Petit et al., 2001). If the observed asteroid belt is ∼5 × 10−4 M⊕ , these values yield a mass range of 0.075– 0.125 Earth masses for D  1000 km objects. Note that these values are likely to be just a tiny fraction of the entire main belt zone’s mass, most which likely went into Moonto Mars-sized planetary embryos. 7.3. Evidence for primordial bodies from asteroid spin rates and lightcurves Over the last 30 years, numerous groups have examined asteroidal rotation rates and lightcurves to discern clues about the collisional history of the main belt. A short list of work on this topic, much of which is still relevant, includes McAdoo and Burns (1973), Harris and Burns (1979), Tedesco and Zappalà (1980), Farinella et al. (1981), Farinella et al. (1982), Dermott et al. (1984), Dobrovolskis and Burns (1984), Binzel et al. (1989), Davis et al. (1989), Cellino et al. (1990), Farinella et al. (1992), Fulchignoni et al. (1995), Donnison and Wiper (1999), Pravec and Harris (2000), and Pravec et al. (2002). As we describe below, we believe these data provide additional circumstantial evidence supporting our results. A tool used by many of the groups listed above is to plot asteroid spin rates vs. diameter using a “running box” mean method. Fig. 2 from Pravec et al. (2002) shows where this method was applied to the spin rates of 984 asteroids. More specifically for our purposes, Pravec et al. examined nearly 400 asteroids with D > 50 km, a size range where the thermal spin up/down mechanism (YORP) is not expected to significantly modify asteroid spin rates (Rubincam, 2000; Vokrouhlický and Čapek, 2002). A plot of these data reveals 133 a minimum in the mean spin rate distribution for D ∼ 90– 120 km asteroids (spin rate = 1.8 d−1 ). For reference, D ∼ 200 km asteroids have ∼3 d−1 while D ∼ 50 km asteroids have ∼2.5 d−1 . The minimum has also been observed in spin rate distributions that show S, C, and M asteroids separately (Dermott et al., 1984; Binzel et al., 1989). A discontinuity is also observed in plots of running box mean lightcurve amplitudes vs. diameter (Binzel et al., 1989). Asteroids with D > 125 km show mean amplitudes near 0.2, with a linear trend vs. log diameter that stretches from 0.21 for D ∼ 200 km to 0.18 for D ∼ 125 km. Asteroids with D < 125 km, however, quickly jump to larger mean amplitudes, with 0.24 for D ∼ 100 km, 0.28 for D ∼ 70 km, and 0.31–0.32 for D < 40 km. It is unclear if this discontinuity is related to the spin rate discontinuity; we caution that to some unknown extent, lightcurve amplitudes must be affected by an asteroid’s self gravity. Still, we find it a surprising coincidence that both discontinuities occur at D ∼ 100–125 km; the simplest scenario would imply a genetic relationship between these two. Many groups have speculated about the cause of these discontinuities. Dobrovolskis and Burns (1984) suggest the spin rate minimum may be caused by an effect called “angular momentum drain,” where asteroid cratering events preferentially lose ejecta in the same direction as the asteroid’s rotation (and thereby slow their spin rate). Cellino et al. (1990) claimed that a similar effect, called “angular momentum splash,” would occur during marginally disruptive collisions. Numerical hydrocode modeling of asteroid breakups, however, suggest a different story. Love and Ahrens (1997) found that the trajectories of the impact ejecta in their simulations were highly directional (mainly downrange), enough that the signal of angular momentum drain or splash could not be determined from their data. Overall, they found that small erosive collisions have a minimal effect on an object’s spin, while catastrophic disruption events essentially destroy all “memory” of the target body’s initial spin. If true, it seems unlikely that the spin rate and lightcurve amplitude features observed in the data were produced by catastrophic collisions. 6 A different solution was proposed by Tedesco and Zappalà (1980) and Dermott and Murray (1982), who suggested this discontinuity marked the dividing line between primordial asteroids and their collisional products. This hypothesis is consistent with predictions made by several pioneering authors (e.g., Kuiper et al., 1958; Anders, 1965; Hartmann and Hartmann, 1968) and with our model results. It also provides some supporting evidence for our prediction of a significant slope change in initial Nrem at Dx ∼ 120 km (Fig. 4). Finally, it would argue that the spin rates and lightcurve amplitudes for most D > 120 km bodies can be 6 It is possible that variants of these effects were important during accretion when planetesimal impact velocities were only a few m s−1 (Leinhardt et al., 2000). 134 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 used to probe the processes that produced planetesimals during runaway accretion. 7.4. Caveats about our best-fit Q∗D disruption function We are encouraged that our best-fit Q∗D functions in the gravity regime (Runs 11–18; see Table 2 and Fig. 9) are similar to results from hydrocode simulations that study the disruption of undamaged basaltic targets (i.e., our Q∗D function for Run 13 is essentially the same as Benz and Asphaug, 1999, for Vimp = 3 km s−1 ; see also Love and Ahrens, 1996; Holsapple et al., 2002; Asphaug et al., 2002). It is not yet clear, however, how our results fit in with recent smoothedparticle hydrocode (SPH) simulations of pre-shattered asteroids that show that these bodies disrupt more easily than undamaged targets (Michel et al., 2001, 2002, 2003). For this reason, it is useful to describe some of the limitations of current SPH modeling and our CoEM results. While SPH models are the best tools we have for simulating asteroid disruption events, there is still room for improvement. A recent summary of some of their limitations can be found in Holsapple et al. (2002). We point out one additional issue here. Up to now, SPH codes have had limited success including macro- and micro-porosity in their asteroid disruption simulations. This is important when you consider that most S- and C-type asteroids have estimated porosities in the range of 15–20 and 30–65%, respectively (Britt et al., 2002). Because porous materials absorb impact energy and prevent tensile waves from propagating through a target body, they can make asteroids more difficult to disrupt (e.g., Asphaug, 1999; Housen and Holsapple, 2003). In fact, preliminary hydrocode tests suggest that pre-shattered asteroids with 25% macroporosity are easier to disrupt than undamaged targets (W. Benz, personal communication). If these results hold up, it would imply that the best-fit Q∗D functions with shallow gravity regime slopes (e.g., Run 11) are more likely to reflect reality than those with steeper slopes (e.g., Run 18). As for CoEM, our Q∗D estimate, for better or worse, is intrinsically linked to the accuracy of the input and constraints included in CoEM. Here we list some of the factors that could degrade the quality of our results: (i) the observed main belt size distribution may differ from our estimates, (ii) our choice of initial size distribution may need further refinement, (iii) our fragment size distribution may not represent the broad spectrum of cratering/disruption events occurring in the main belt, (iv) the approximations made to model the DDE in CoEM may be overly-simplistic, and (v) the true Q∗D function may not have a simple hyperbolic shape or it may be velocity-dependant (though see Appendix A). To mitigate against these potential problems, we have used goodness-of-fit metrics that allow for some variability 2 ; χ2 in our constraints (e.g., ψSFD FAM ). Still, future modeling work is needed to determine whether our best-fit Q∗D functions are as accurate as they can be. 8. Conclusions In this paper, we created a collisional evolution model (CoEM) capable of tracking main belt comminution from the end of accretion among D < 1000 km bodies to the present day. Our method accounted for the possibility that main belt population was once far more massive than the current population and that it lost the majority of its mass via a dynamical depletion event (DDE). We bypassed questions about the initial size of the main belt population and the timing/nature of the DDE by assuming that: (i) main belt comminution among diameter D < 1000 km bodies has been dominated by the same collision probabilities and impact velocities found in the main belt today, and (ii) a large planetesimal population undergoing comminution for a short period of time is mathematically equivalent to a much smaller population undergoing comminution for an extended period of time. Using CoEM, we tested possible solutions for the asteroid disruption scaling-law (i.e., the critical impact specific energy Q∗D ) and the shape of the initial main belt population. Constraints for CoEM came from the main belt’s size-frequency distribution, the existence of a D = 460 km crater on Asteroid (4) Vesta, the number of asteroid families produced by D > 100 km disruption events over the last 3–4 Gyr, and the relatively constant crater production rate of the Earth and Moon over the last 3 Gyr. These constraints helped drive our results toward a unique solution. We summarize our major findings: • Best fit solution for Q∗D . Our best fit solutions for the Q∗D functions are described by Runs 11–18 (Table 2; Fig. 9). Equation’s (5) parameters for Run 11 are E = 0.871, F = −0.877, G = −0.489, and H = −0.299, while those for Run 18 are E = 0.854, F = −0.940, G = −0.511, and H = −0.314 (see Table 2 and Fig. 9). For reference, the Q∗D function predicted by the hydrocode modeling results of Benz and Asphaug (1999) is represented by Run 13, located in the middle of our best fit results. The positive match gives us increased confidence that the method, constraints, and results described in this paper are accurate. • Shape of the initial main belt size distribution. Using our model’s assumptions, we estimated the shape of the initial main belt population that existed prior to the onset of fragmentation among D < 1000 km bodies. The initial main belt size distribution for D < 1000 km bodies was divided into two components, Nrem and Ndep . We solved for the former, which is a stand-in for the current main belt population. The latter is a hypothetical population that may have contributed impactors to the primordial main belt before it was dynamically eliminated. Until this removal, Ndep would have had the same shape as Nrem . After numerous trials using both CoEM-SM and CoEM-ST, we found the best fit Nrem initial population had nearly the same number of D
120 km asteroids as The fossilized main belt size distribution • • • • the observed main belt (incremental power law index of −4.5) and a much more limited number of D  120 km asteroids. We found the power law slope of the D  120 km asteroid did not change out results, provided it was shallow enough not to exceed observed number of D ∼ 50–100 km asteroids. This shape is relatively consistent with predictions made by several pioneering papers from the 1950 and 1960’s (Kuiper et al., 1958; Anders, 1965; Hartmann and Hartmann, 1968). Degree of collisional evolution in main belt. The first time our best-fit CoEM runs for Nrem matched our constraints was at a median pseudo-time of 7.5–9.5 Gyr. We interpret this to mean that the main belt size distribution could not have attained its current wavy shape without going through a period where it was exposed to many more projectiles than are observed today. Accordingly, these results support the idea that the main belt was once more massive than it is today, with much of that mass lost via a dynamical mechanism rather than comminution (i.e., Ndep > 0). Hence, the wavy main belt size distribution is predominately a “fossil” produced by collisional evolution in the primordial main belt. Our results also suggest that most D
120 km objects have never been disrupted, while many, perhaps most D  120 km asteroids are byproducts of fragmentation events among D
120 km asteroids. Stability of main belt and NEO populations. Our best fit models suggest that once the shape of the main belt size distribution approaches the observed population, it will remain close to those values for several Gyr. Because NEO population is predominately replenished by D < 30 km main belt asteroids via the Yarkovsky effect (Bottke et al., 2000, 2002a, 2002b), this result explains why the impactor flux on the lunar maria has been nearly constant for the last ∼3 Gyr (e.g., Grieve and Shoemaker, 1994). Check of results using large crater on (4) Vesta. The intact basaltic crust of (4) Vesta has superposed on it a D = 460 km crater. This feature was formed by the singular impact of a D ∼ 35 km projectile, making it a better constraint for our model results than Vesta’s intact crust. Using our best-fit model runs, we estimate that two such impactors should have struck Vesta over a pseudo-time of tpsuedo > 27 Gyr. This value suggests the degree of comminution experienced by the main belt had to be less than 6 times the amount it would have received in the current low-mass main belt over the last 4.5 Gyr. These values are consistent with our estimates that the main belt experienced roughly the equivalent of tpseudo ∼ 7.5–9.5 Gyr of collisional evolution. The probability that any D ∼ 35 km impactor struck Vesta over the last 3.5 Gyr of real time and produced the Vesta family is ∼20%. Asteroid disruption frequency. We found that approximately four D > 100 km objects disrupt every Gyr, and that Karin-size breakups (25 km < D < 35 km; 135 Nesvorný et al., 2002c; Nesvorný and Bottke, 2004) occur ∼15–30 times per Gyr. • Constraints from asteroid spin rates and lightcurve data. Supporting evidence for our claim that most D
120 km objects are primordial comes from asteroid spin rates and lightcurve data. A running box mean of asteroid spin rate vs. diameter shows a minima near D ∼ 100–120 km, while a running box mean of lightcurve amplitudes vs. diameter show a discontinuity near D ∼ 125 km. Given our model results, we believe the simplest explanation is that D
120 km bodies are predominantly unaffected by catastrophic collisions, while the D  120 km population is increasingly dominated by collisional fragments as one goes to smaller and smaller sizes. • The estimated size of the primordial main belt. Using our best-fit model runs and a scenario where a massive main belt went through an intense phase of early comminution prior to the formation of Jupiter, we estimate that the main belt population (in the form of D < 1000 km bodies) was once 150–250 times larger than it is today (0.075–0.125M⊕ ). These values are the same as those predicted by numerical models of the dynamical excitation and clearing of the primordial main belt population (Petit et al., 2001). They are also significantly lower than the several Earth masses of material predicted by solar nebula models, suggesting the remaining mass was taken up by the growth of numerous Moon- to Mars-sized planetary embryos in the main belt zone. Acknowledgments We thank Ed Scott, Patrick Michel, Dave O’Brien, and Kleomenis Tsiganis for valuable discussions and input to this study. We also thank Stu Weidenschilling and Sarah Andre for their careful and constructive reviews of this paper. Research funds for William Bottke were provided by NASA’s Origins of Solar Systems Program (Grant NAG510658) and NASA’s Discovery Data Analysis Program (grant NNG04GA75G). Appendix A. Exploring how varying impact velocities affect our results In this paper, we concentrated on CoEM simulations where Vimp = 5.3 km s−1 . A potential problem with this approximation is that main belt impact velocities could have been very different in the past. For example, it is plausible that small asteroids could have fragmented at low velocities while the largest bodies in the main belt were still forming. Alternatively, when the Ndep population was ejected from the main belt, collision velocities between those asteroids and the Nrem population were probably comparable to 136 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 10 km s−1 (e.g., Petit et al., 2001). If either phase was critically important to the collisional evolution of the main belt, the model results and Q∗D function derived here may be inaccurate. Moreover, no one knows precisely how Q∗D varies with velocity. Some hydrocode results suggest that asteroids readily disrupt when struck at Vimp < 1 km s−1 (Benz, 2000). To investigate how our results are affected by varying impact velocities, we used CoEM-SM to set up a series of “toy” simulations where we could compute the best possible Q∗D functions and evolution timescales for a population being bombarded by projectiles with Vimp = 1, 3, 5.3, and 10.6 km s−1 . Our goal was not to generate a realistic simulation, but instead to explore how far Q∗D would change itself in an extreme situation to reproduce the observed main belt. To make it as easy as possible to compare our runs, we set Pi to a clearly unphysical value of 2.86 × 10−18 km−2 yr−1 for all cases (note that in reality, each Vimp value corresponds to its own Pi value). Our initial Nrem population was set to Dx = 120 km (Fig. 4). We found that for each choice of Vimp , CoEM-SM found a Q∗D function that, while very different from one another, always yielded a best fit to our constraints at tpseudo ∼10 Gyr. Fig. 13 shows the derived Q∗D function for each Vimp value. We see that the Q∗D function for Vimp = 3 and 10.6 km s−1 are within a factor of 3 of the reference 5.3 km s−1 run, but that the 1 km s−1 is nearly an order of magnitude lower than the reference run. Fig. 13. Best fit Q∗D functions from our toy CoEM-SM runs where Vimp = 1, 3, 5.3, and 10.6 km s−1 . We used the Dx = 120 km initial size distribution (Fig. 4) and set Pi to the non-physical value of 2.86 × 10−18 km−2 yr−1 for all our velocities. The evolution timescale was tpseudo ∼ 10 Gyr. These runs determine how far Q∗D must vary to produce the observed main belt under our chosen conditions; they were not designed to be realistic. Each simulation produces a comparable number of asteroid disruption events. The reason our evolution timescales stayed constant across our toy simulations was that the observed main belt size distribution can only be reproduced if there are just the right number of large-scale disruption events at particular intervals over time; too many or too few yield size distributions that are inconsistent with observations. These results suggest that if main belt impact velocities were once wildly different than they are today, Q∗D had to be wildly different as well, and in just the right way, in order to limit the number of large-scale breakup events, or that relatively few disruption events occurred during that period. Either way, it suggests our results are probably a reasonable approximation of collisional evolution in the main belt. Another way to understand these results is as follows. If the number of disruption events needed to make the main belt size distribution does not change while Pi is held constant and Vimp varies, a target asteroid struck at a lower Vimp value can only be disrupted by the same sized impactor if Q∗D decreases as well. Figure 13 shows that, to a good approximation, Q∗D in the gravity scaling regime varies as 2 . That is, the ratio of impact energy to that required for Vimp the target’s disruption must remain essentially constant in order to produce the same outcome. A different question to ask is whether real Q∗D functions change in the radical ways suggested by Fig. 13. The limited evidence we have today suggests it is unlikely. Benz and Asphaug (1999) used hydrocode simulations to investigate collisions onto undamaged basaltic target bodies being struck at Vimp = 3 and 5 km s−1 and for undamaged ice target bodies being struck at Vimp = 0.5 and 3 km s−1 . Their results suggest that Q∗D in the gravity regime only varies by a factor of 2 in each case, significantly smaller than the factor of 30 difference between our Vimp = 1 and 5.3 km s−1 cases. If true, these results suggest that an extended period of low velocity impacts among planetesimals during the accretion of planetary embryos would not lead to significant collisional erosion. A caveat to these results comes from Benz (2000), who investigated Q∗D for planetesimal collisions occurring at very low velocities (Vimp = 5 and 40 m s−1 ). He found that for such low velocity impacts, Q∗D in the gravity regime could drop by roughly an order of magnitude over values derived for 3 km s−1 . It is unclear from planet formation codes whether D < 1000 km planetesimals experienced impact velocities in this range for any extended period of time. Nevertheless, because these impact velocities and Q∗D functions would be subject to the same constraints on the number and timing of asteroid breakups as the results described above, we do not believe they modify our conclusions. References Agnor, C.B., Canup, R.M., Levison, H.F., 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142, 219–237. Anders, E., 1965. Fragmentation history of asteroids. Icarus 4, 399–408. The fossilized main belt size distribution Asphaug, E., 1997. Impact origin of the Vesta family. Meteorit. Planet. Sci. 32, 965–980. Asphaug, E., 1999. Survival of the weakest. Nature 402, 127–128. Asphaug, E., Ryan, E.V., Zuber, M.T., 2002. Asteroid interiors. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 463–484. Bendjoya, P., Zappalà, V., 2002. Asteroid family identification. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 613–618. Benz, W., 2000. Low velocity collisions and the growth of planetesimals. In: Benz, W., Kallenbach, R., Lugmair, G.W. (Eds.), From Dust to Terrestrial Planets. Kluwer Academic, Dordrecht, pp. 279–300. Benz, W., Asphaug, E., 1999. Catastrophic disruptions revisited. Icarus 142, 5–20. Binzel, R.P., Farinella, P., Zappalà, V., Cellino, A., 1989. Asteroid rotation rates—distributions and statistics. In: Binzel, R.P., Gehrels, T., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, pp. 416–441. Bogard, D., 1995. Impact ages of meteorites: a synthesis. Meteoritics 30, 244–268. Bottke, W.F., Greenberg, R., 1993. Asteroidal collision probabilities. Geophys. Res. Lett. 20, 879–881. Bottke, W.F., Nolan, M.C., Greenberg, R., Kolvoord, R.A., 1994. Velocity distributions among colliding asteroids. Icarus 107, 255–268. Bottke Jr., W.F., Rubincam, D.P., Burns, J.A., 2000. Dynamical evolution of main belt meteoroids: numerical simulations incorporating planetary perturbations and Yarkovsky thermal forces. Icarus 145, 301–331. Bottke, W.F., Vokrouhlický, D., Broz̆, M., Nesvorný, D., Morbidelli, A., 2001. Dynamical spreading of asteroid families by the Yarkovsky effect. Science 294, 1693–1696. Bottke, W.F., Morbidelli, A., Jedicke, R., Petit, J., Levison, H.F., Michel, P., Metcalfe, T.S., 2002a. Debiased orbital and absolute magnitude distribution of the near-Earth objects. Icarus 156, 399–433. Bottke, W.F., Vokrouhlický, D., Rubincam, D., Brož, M., 2002b. Dynamical evolution of asteroids and meteoroids using the Yarkovsky effect. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 395–408. Britt, D.T., Yeomans, D., Housen, K., Consolmagno, G., 2002. Asteroid density, porosity, and structure. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 485–500. Burbine, T.H., Buchanan, P.C., Binzel, R.P., Bus, S.J., Hiroi, T., Hinrichs, J.L., Meibom, A., McCoy, T.J., 2001. Vesta, Vestoids, and the howardite, eucrite, diogenite group: relationships and the origin of spectral differences. Meteorit. Planet. Sci. 36, 761–781. Carruba, V., Burns, J.A., Bottke, W., Nesvorný, D., 2003. Orbital evolution of the Gefion and Adeona asteroid families: close encounters with massive asteroids and the Yarkovsky effect. Icarus 162, 308–327. Campo Bagatin, A., Cellino, A., Davis, D.R., Farinella, P., Paolicchi, P., 1994. Wavy size distributions for collisional systems with a small-size cutoff. Planet. Space Sci. 42, 1079–1092. Campo Bagatin, A., Petit, J.-M., Farinella, P., 2001. How many rubble piles are in the asteroid belt? Icarus 149, 198–209. Cellino, A., Zappalà, V., Davis, D.R., Farinella, P., Paolicchi, P., 1990. Asteroid collisional evolution I. Angular momentum splash. Despinning asteroids through catastrophic collisions. Icarus 87, 391–402. Cellino, A., Zappalà, V., Farinella, P., 1991. The size distribution of mainbelt asteroids from IRAS data. Mon. Not. R. Astron. Soc. 253, 561–574. Chambers, J.E., Wetherill, G.W., 1998. Making the terrestrial planets: N -body integrations of planetary embryos in three dimensions. Icarus 136, 304–327. Chambers, J.E., Wetherill, G.W., 2001. Planets in the asteroid belt. Meteorit. Planet. Sci. 36, 381–399. Charnoz, S., Morbidelli, A., 2003. Coupling dynamical and collisional evolution of small bodies: an application to the early ejection of planetesimals from the Jupiter–Saturn region. Icarus 166, 141–156. Cheng, A.F., 2004. Collisional evolution of the asteroid belt. Icarus 169, 357–372. 137 Colwell, J.E., 1993. Power-law confusion: you say incremental, I say differential. Lunar Planet. Sci. 24, 325. Abstract. Culler, T.S., Becker, T.A., Muller, R.A., Renne, P.R., 2000. Lunar impact history from 40 Ar/39 Ar dating of glass spherules. Science 287, 1785– 1788. Davis, D.R., Chapman, C.R., Greenberg, R., Weidenschilling, S.J., Harris, A.W., 1979. Collisional evolution of asteroids. Populations, rotations, and velocities. In: Gehrels, T. (Ed.), Asteroids. Univ. of Arizona Press, Tucson, pp. 528–557. Davis, D.R., Chapman, C.R., Weidenschilling, S.J., Greenberg, R., 1985. Collisional history of asteroids: evidence from Vesta and the Hirayama families. Icarus 63, 30–53. Davis, D.R., Weidenschilling, S.J., Farinella, P., Paolicchi, P., Binzel, R.P., 1989. Asteroid collisional history. Effects on sizes and spins. In: Binzel, R.P., Gehrels, T., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, pp. 805–826. Davis, D.R., Ryan, E.V., Farinella, P., 1994. Asteroid collisional evolution: results from current scaling algorithms. Planet. Space Sci. 42, 599–610. Davis, D.R., Durda, D.D., Marzari, F., Campo Bagatin, A., Gil-Hutton, R., 2002. Collisional evolution of small body populations. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 545–558. Dermott, S.F., Murray, C.D., 1982. Asteroid rotation rates depend on diameter and type. Nature 296, 418–421. dell’Oro, A., Paolicchi, P., 1998. Statistical properties of encounters among asteroids: a new, general purpose, formalism. Icarus 136, 328–339. Dermott, S.F., Harris, A.W., Murray, C.D., 1984. Asteroid rotation rates. Icarus 57, 14–34. Dobrovolskis, A.R., Burns, J.A., 1984. Angular momentum drain. A mechanism for despinning asteroids. Icarus 57, 464–476. Dohnanyi, J.W., 1969. Collisional models of asteroids and their debris. J. Geophys. Res. 74, 2531–2554. Donnison, J.R., Wiper, M.P., 1999. Bayesian statistical analysis of asteroid rotation rates. Mon. Not. R. Astron. Soc. 302, 75–80. Durda, D.D., 1993. The collisional evolution of the asteroid belt and its contribution to the zodiacal cloud. PhD thesis. Durda, D.D., Dermott, S.F., 1997. The collisional evolution of the asteroid belt and its contribution to the zodiacal cloud. Icarus 130, 140–164. Durda, D.D., Greenberg, R., Jedicke, R., 1998. Collisional models and scaling laws: a new interpretation of the shape of the main-belt asteroid size distribution. Icarus 135, 431–440. Durda, D.D., Flynn, G.J., Hart, S.D., Asphaug, E., 2002. Impact disruption of three ordinary chondrite meteorites. Lunar Planet. Sci. 33. Abstract 1535. Durda, D.D., Bottke, W.F., Enke, B.L., Merline, W.J., Asphaug, E., Richardson, D.C., Leinhardt, Z.M., 2004a. The formation of asteroid satellites in large impacts: results from numerical simulations. Icarus 170, 243– 257. Durda, D.D., Bottke, W.F., Enke, B.L., Merline, W.J., Asphaug, E., Richardson, D.C., Leinhardt, Z.M., 2004b. Erratum to “The formation of asteroid satellites in large impacts: results from numerical simulations [Icarus 167 (2004) 382–396]. Icarus 170, 242. Farinella, P., Davis, D.R., 1992. Collision rates and impact velocities in the main asteroid belt. Icarus 97, 111–123. Farinella, P., Vokrouhlický, D., 1999. Semimajor axis mobility of asteroidal fragments. Science 283, 1507–1510. Farinella, P., Paolicchi, P., Zappalà, V., 1981. Analysis of the spin rate distribution of asteroids. Astron. Astrophys. 104, 159–165. Farinella, P., Paolicchi, P., Zappalà, V., 1982. The asteroids as outcomes of catastrophic collisions. Icarus 52, 409–433. Farinella, P., Davis, D.R., Paolicchi, P., Cellino, A., Zappalà, V., 1992. Asteroid collisional evolution—an integrated model for the evolution of asteroid rotation rates. Astron. Astrophys. 253, 604–614. Fowler, J.W., Chillemi, J.R., 1992. IRAS asteroid data processing. In: Tedesco, E.F. (Ed.), The IRAS Minor Planet Survey. Phillips Laboratory, Hanscom Air Force Base, MA, pp. 17–43. Tech. Report PL-TR92-2049. 138 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 Franklin, F., Lecar, M., 2000. On the transport of bodies within and from the asteroid belt. Meteorit. Planet. Sci. 35, 331–340. Fulchignoni, M., Barucci, M.A., di Martino, M., Dotto, E., 1995. On the evolution of the asteroid spin. Astron. Astrophys. 299, 929–932. Gomez, R., Tsiganis, T., Morbidelli, A., Levison, H.F., 2004. The late heavy bombardment as a cataclysimic event caused by a delayed planetary migration. Nature. Submitted for publication. Greenberg, R., 1982. Orbital interactions. A new geometrical formalism. Astron. J 87, 184–195. Greenberg, R., Nolan, M.C., 1989. Delivery of asteroids and meteorites to the inner Solar System. In: Binzel, R.P., Gehrels, T., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, pp. 778–801. Greenberg, R., Hartmann, W.K., Chapman, C.R., Wacker, J.F., 1978. Planetesimals to planets. Numerical simulation of collisional evolution. Icarus 35, 1–26. Grier, J.A., McEwen, A.S., Lucey, P.G., Milazzo, M., Strom, R.G., 2001. Optical maturity of ejecta from large rayed lunar craters. J. Geophys. Res. 106, 32847–32862. Grieve, R.A.F., Shoemaker, E.M., 1994. The record of past impacts on Earth. In: Gehrels, T., Matthews, M.S. (Eds.), Hazards Due to Comets and Asteroids. Univ. of Arizona Press, Tucson, pp. 417–462. Hartmann, W.K., Hartmann, A.C., 1968. Asteroid collisions and evolution of asteroidal mass distribution and meteoritic flux. Icarus 8, 361–381. Hartmann, W.K., Ryder, G., Dones, L., Grinspoon, D., 2000. The timedependent intense bombardment of the primordial Earth/Moon system. In: Canup, R., Righter, K. (Eds.), Origin of the Earth and the Moon. Univ. of Arizona Press, Tucson, pp. 805–826. Harris, A.W., Burns, J.A., 1979. Asteroid rotation. I. Tabulation and analysis of rates, pole positions and shapes. Icarus 40, 115–144. Hellyer, B., 1970. The fragmentation of the asteroids. Mon. Not. R. Astron. Soc. 115, 148–383. Hellyer, B., 1971. The fragmentation of the asteroids. II. Numerical calculations. Mon. Not. R. Astron. Soc. 154, 279–284. Hörz, F., 2000. Time-variable cratering rates? Science 288, 2095a. Housen, K.R., Holsapple, K.A., 2003. Impact cratering on porous asteroids. Icarus 163, 102–119. Holsapple, K., Giblin, I., Housen, K., Nakamura, A., Ryan, E., 2002. Asteroid impacts: laboratory experiments and scaling laws. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R.P. (Eds.), Asteroids III. Univ. of Arizona Press, Tucson, pp. 443–462. Inaba, S., Wetherill, G.W., Ikoma, M., 2003. Formation of gas giant planets: core accretion models with fragmentation and planetary envelope. Icarus 166, 46–62. Ivezić, Ž., 32 colleagues, 2001. Solar System objects observed in the Sloan Digital Sky Survey commissioning data. Astron. J. 122, 2749–2784. Ivezić, Ž., Lupton, R.H., Jurić, M., Tabachnik, S., Quinn, T., Gunn, J.E., Knapp, G.R., Rockosi, C.M., Brinkmann, J., 2002. Color confirmation of asteroid families. Astron. J. 124, 2943–2948. Jedicke, R., Metcalfe, T.S., 1998. The orbital and absolute magnitude distributions of main belt asteroids. Icarus 131, 245–260. Jedicke, R., Larsen, J., Spahr, T., 2002. Observational selection effects in asteroid surveys. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 71–88. Jurić, M., 15 colleagues, 2002. Comparison of positions and magnitudes of asteroids observed in the Sloan Digital Sky Survey with those predicted for known asteroids. Astron. J. 124, 1776–1787. Keil, K., 2002. Geological history of Asteroid 4 Vesta: the “smallest terrestrial planet.” In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 573–584. Keil, K., Haack, H., Scott, E.R.D., 1994. Catastrophic fragmentation of asteroids: evidence from meteorites. Planet. Space Sci. 42, 1109–1122. Kenyon, S.J., 2002. Planet formation in the outer Solar System. Publ. Astron. Soc. Pacific 114, 265–283. Kenyon, S.J., Bromley, B.C., 2004. Detecting the dusty debris of terrestrial planet formation. Astrophys. J. 602, L133–L136. Knežević, Z., Lemaître, A., Milani, A., 2002. The determination of asteroid proper elements. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 603–612. Kokubo, E., Ida, S., 2002. Formation of protoplanet systems and diversity of planetary systems. Astrophys. J. 581, 666–680. Kortenkamp, S.J., Wetherill, G.W., Inaba, S., Trilling, D.E., 2001. Asteroid formation with a pre-existing Jupiter. Lunar Planet. Sci. 32. Abstract 1796. Kuiper, G.P., Fugita, Y., Gehrels, T., Groeneveld, I., Kent, J., van Biesbroeck, G., van Houten, C.J., 1958. Survey of asteroids. Astrophys. J. Suppl. Ser. 32 (III), 289–428. Levison, H.F., Dones, L., Chapman, C.R., Stern, S.A., Duncan, M.J., Zahnle, K., 2001. Could the lunar “late heavy bombardment” have been triggered by the formation of Uranus and Neptune? Icarus 151, 286– 306. Leinhardt, Z.M., Richardson, D.C., 2002. N -body simulations of planetesimal evolution: effect of varying impactor mass ratio. Icarus 159, 306– 313. Leinhardt, Z.M., Richardson, D.C., Quinn, T., 2000. Direct N -body simulations of rubble pile collisions. Icarus 146, 133–151. Lissauer, J.J., 1987. Timescales for planetary accretion and the structure of the protoplanetary disk. Icarus 69, 249–265. Love, S.G., Ahrens, T.J., 1996. Catastrophic impacts on gravity dominated asteroids. Icarus 124, 141–155. Love, S.G., Ahrens, T.J., 1997. Origin of asteroid rotation rates in catastrophic impacts. Nature 386, 154–156. McAdoo, D.C., Burns, J.A., 1973. Further evidence for collisions among asteroids. Icarus 18, 285–293. McEwen, A.S., Moore, J.M., Shoemaker, E.M., 1997. The Phanerozoic impact cratering rate: evidence from the farside of the Moon. J. Geophys. Res. 102, 9231–9242. Manley, S.P., Migliorini, F., Bailey, M.E., 1998. An algorithm for determining collision probabilities between small Solar System bodies. Astron. Astrophys. Suppl. Ser. 133, 437–444. Marzari, F., Davis, D., Vanzani, V., 1995. Collisional evolution of asteroid families. Icarus 113, 168–187. Marzari, F., Cellino, A., Davis, D.R., Farinella, P., Zappalà, V., Vanzani, V., 1996. Origin and evolution of the Vesta asteroid family. Astron. Astrophys. 316, 248–262. Marzari, F., Farinella, P., Davis, D.R., 1999. Origin, aging, and death of asteroid families. Icarus 142, 63–77. Michel, P., Benz, W., Tanga, P., Richardson, D., 2001. Collisions and gravitational reaccumulation: a recipe for forming asteroid families and satellites. Meteorit. Planet. Sci. 294, 1696–1700. Michel, P., Tanga, P., Benz, W., Richardson, D.C., 2002. Formation of asteroid families by catastrophic disruption: simulations with fragmentation and gravitational reaccumulation. Icarus 160, 10–23. Michel, P., Benz, W., Richardson, D.C., 2003. Disruption of fragmented parent bodies as the origin of asteroid families. Nature 421, 608–611. Milani, A., Knežević, Z., 1994. Asteroid proper elements and the dynamical structure of the asteroid main belt. Icarus 107, 219–254. Morbidelli, A., Vokrouhlický, D., 2003. The Yarkovsky-driven origin of near-Earth asteroids. Icarus 163, 120–134. Morbidelli, A., Chambers, J., Lunine, J.I., Petit, J.M., Robert, F., Valsecchi, G.B., Cyr, K.E., 2000. Source regions and time scales for the delivery of water to Earth. Meteorit. Planet. Sci. 35, 1309–1320. Morbidelli, A., Nesvorný, D., Bottke, W.F., Michel, P., Vokrouhlický, D., Tanga, P., 2003. The shallow magnitude distribution of asteroid families. Icarus 162, 328–336. Nakamura, T., Yoshida, F., 2002. Statistical method for deriving spatial and size distributions of sub-km main-belt asteroids from their sky motions. Publ. Astron. Soc. Japan 54, 1079–1089. Nagasawa, M., Ida, S., Tanaka, H., 2002. Excitation of orbital inclinations of asteroids during depletion of a protoplanetary disk: dependence on the disk configuration. Icarus 159, 322–327. Nesvorný, D., Bottke, W.F., 2004. Detection of the Yarkovsky effect for main-belt asteroids. Icarus 170, 324–342. Nesvorný, D., Morbidelli, A., Vokrouhlický, D., Bottke, W.F., Brož, M., 2002a. The Flora family: a case of the dynamically dispersed collisional swarm? Icarus 157, 155–172. The fossilized main belt size distribution Nesvorný, D., Ferraz-Mello, S., Holman, M., Morbidelli, A., 2002b. Regular and chaotic dynamics in the mean-motion resonances: implications for the structure and evolution of the asteroid belt. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 379–394. Nesvorný, D., Bottke, W.F., Dones, L., Levison, H.F., 2002c. The recent breakup of an asteroid in the main-belt region. Nature 417, 720–771. Nesvorný, D., Bottke, W.F., Levison, H.F., Dones, L., 2003. Recent origin of the Solar System dust bands. Astrophys. J. 591, 486–497. Nesvorný, D., Jedicke, R., Whiteley, R.J., Ivezić, Ž., 2005. Evidence for asteroid space weathering from the Sloan Digital Sky Survey. Icarus 173, 132–152. Neukum, G., Ivanov, B.A., 1994. Crater size distributions and impact probabilities on Earth from lunar, terrestrial-planet, and asteroid cratering data. In: Gehrels, T., Matthews, M.S. (Eds.), Hazards Due to Comets and Asteroids. Univ. of Arizona Press, Tucson, pp. 359–416. O’Brien, D.P., Greenberg, R., 2003. Steady-state size distributions for collisional populations: analytical solution with size-dependent strength. Icarus 164, 334–345. O’Brien, D.P. Greenberg, R., 2004. The collisional and dynamical evolution of the main belt and NEA size distributions. Icarus. Submitted for publication. Öpik, E.J., 1951. Collision probability with the planets and the distribution of planetary matter. Proc. R. Irish Acad. 54, 165–199. Petit, J., Farinella, P., 1993. Modelling the outcomes of high-velocity impacts between small Solar System bodies. Celest. Mech. Dynam. Astron. 57, 1–21. Petit, J., Morbidelli, A., Chambers, J., 2001. The primordial excitation and clearing of the asteroid belt. Icarus 153, 338–347. Petit, J.-M., Chambers, J., Franklin, F., Nagasawa, M., 2002. Primordial excitation and depletion of the main belt. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 711–723. Paolicchi, P., 1994. Rushing to equilibrium: a simple model for the collisional evolution of asteroids. Planet. Space Sci. 42, 1093–1097. Pollack, J.B., Hubickyj, O., Bodenheimer, P., Lissauer, J.J., Podolak, M., Greenzweig, Y., 1996. Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124, 62–85. Pravec, P., Harris, A.W., 2000. Fast and slow rotation of asteroids. Icarus 148, 12–20. Pravec, P., Harris, A.W., Michalowski, T., 2002. Asteroid rotations. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 113–122. Press, W.H., Flannery, B.P., Teukolsky, S.A., Vetterling, W.T., 1989. Numerical Recipes in Fortran. The Art of Scientific Computing. Cambridge Univ. Press, Cambridge. Richardson, D.C., Leinhardt, Z.M., Melosh, H.J., Bottke, W.F., Asphaug, E., 2002. Gravitational aggregates: evidence and evolution. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 501–515. Rubincam, D.P., 2000. Radiative spin-up and spin-down of small asteroids. Icarus 148, 2–11. Safronov, V.S., 1969. Evolution of Protoplanetary Cloud and Formation of the Earth and Planets. Nauka, Moscow. Scott, E.R.D., 2002. Meteorite evidence for the accretion and collisional evolution of asteroids. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 697– 709. Scott, E.R.D., 2004. Meteoritic constraints on collision rates in the primordial asteroid belt and its origin. Lunar Planet. Sci. 35. Abstract 1990. Scott, E.R.D., Haack, H., Love, S.G., 2001. Formation of mesosiderites by fragmentation and reaccretion of a large differentiated asteroid. Meteorit. Planet. Sci. 36, 869–891. Shoemaker, E.M., 1998. Long term variations in the impact cratering rate on Earth. In: Grady, M., Hutchison, R., McCall, G.J.H., Rothery, D.A. (Eds.), Meteorites: Flux with Time and Impact Effects. Special Publications, vol. 140. Geological Society, London, pp. 7–10. 139 Shoemaker, E.M., Shoemaker, C.S., 1996. The proterozoic impact record in Australia. AGSO. J. Aust. Geol. Geophys. 16, 379–398. Shukolyukov, A., Lugmair, G.W., 2002. Chronology of asteroid accretion and differentiation. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 687–695. Standish, E.M., JPL Interoffice Memorandum 312.F-01-006, April 11, 2001. Stern, S.A., 1996. On the collisional environment, accretion time scales, and architecture of the massive, primordial Kuiper belt. Astron. J. 112, 1203–1210. Tanaka, H., Inaba, S., Nakazawa, K., 1996. Steady-state size distribution for the self-similar collision cascade. Icarus 123, 450–455. Tanga, P., Cellino, A., Michel, P., Zappalà, V., Paolicchi, P., dell’Oro, A., 1999. On the size distribution of asteroid families: the role of geometry. Icarus 141, 65–78. Tedesco, E.F., Desert, F., 2002. The Infrared Space Observatory Deep Asteroid Search. Astron. J. 123, 2070–2082. Tedesco, E.F., Zappalà, V., 1980. Rotational properties of asteroids. Correlations and selection effects. Icarus 43, 33–50. Tedesco, E.F., Noah, P.V., Noah, M., Price, S.D., 2002. The supplemental IRAS Minor Planet Survey. Astron. J. 123, 1056–1085. Thomas, P.C., Binzel, R.P., Gaffey, M.J., Storrs, A.D., Wells, E.N., Zellner, B.H., 1997. Impact excavation on Asteroid 4 Vesta: Hubble Space Telescope results. Science 277, 1492–1495. van Houten, C.J., van Houten-Groeneveld, I., Herget, P., Gehrels, T., 1970. The Palomar–Leiden survey of faint minor planets. Astron. Astrophys. Suppl. Ser. 2, 339. Vedder, J.D., 1998. Main belt asteroid collision probabilities and impact velocities. Icarus 131, 283–290. Veeder, G.J., Tedesco, E.F., Matson, D.L., 1989. Asteroid results from the IRAS survey. In: Binzel, R.P., Gehrels, T.J., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, pp. 282–289. Vokrouhlický, D., Čapek, D., 2002. YORP-induced long-term evolution of the spin state of small asteroids and meteoroids: Rubincam’s approximation. Icarus 159, 449–467. Weidenschilling, S.J., 1977. The distribution of mass in the planetary system and solar nebula. Astrophys. Space Sci. 51, 153–158. Weidenschilling, S.J., 2000. Formation of planetesimals and accretion of the terrestrial planets. Space Sci. Rev. 92, 295–310. Weidenschilling, S.J., Spaute, D., Davis, D.R., Marzari, F., Ohtsuki, K., 1997. Accretional evolution of a planetesimal swarm. Icarus 128, 429– 455. Weidenschilling, S.J., Cuzzi, J.N., 2004. Accretion dynamics and timescales: relation to chondrites. In: Lauretta, D., Leshin, L.A., McSween, H. (Eds.), Meteorites in the Earth Solar System II. Univ. of Arizona Press, Tucson. In press. Wetherill, G.W., 1967. Collisions in the asteroid belt. J. Geophys. Res. 72, 2429–2444. Wetherill, G.W., 1989. Origin of the asteroid belt. In: Binzel, R.P., Gehrels, T., Matthews, M.S. (Eds.), Asteroids II. Univ. of Arizona Press, Tucson, pp. 661–680. Wetherill, G.W., 1992. An alternative model for the formation of the asteroids. Icarus 100, 307–325. Wetherill, G.W., Inaba, S., 2000. Planetary accumulation with a continuous supply of planetesimals. Space Sci. Rev. 92, 311–320. Wetherill, G.W., Stewart, G.R., 1989. Accumulation of a swarm of small planetesimals. Icarus 77, 330–357. Wilhelms D.E., McCauley J.F., Trask N.J., The geologic history of the Moon. USGS Professional Paper 1348. Williams, D.R., Wetherill, G.W., 1994. Size distribution of collisionally evolved asteroidal populations. Analytical solution for self-similar collision cascades. Icarus 107, 117–128. Yoshida, F., Nakamura, T., Watanabe, J., Kinoshita, D., Yamamoto, N., Fuse, T., 2003. Size and spatial distributions of sub-km main-belt asteroids. Publ. Astron. Soc. Japan 55, 701–715. Zappalà, V., Cellino, A., dell’Oro, A., Paolicchi, P., 2002. Physical and dy- 140 W.F. Bottke Jr. et al. / Icarus 175 (2005) 111–140 namical properties of asteroid families. In: Bottke, W.F., Cellino, A., Paolicchi, P., Binzel, R. (Eds.), Asteroids III. Univ. Arizona Press, Tucson, pp. 619–631. Zellner, N.E.B., Spudis, P.D., Delano, J.W., Whittet, D.C.B., Swindle, T.D., 2003. Geochemistry and impact history at the Apollo 16 landing site. Lunar Planet. Sci. 34, 1157–1158.