Are Main-Belt Asteroids a Sufficient Source for the Earth-Approaching Asteroids?

ICARUS

127, 33–54 (1997) IS965668

ARTICLE NO.

Are Main-Belt Asteroids a Sufficient Source for the Earth-Approaching Asteroids? I. Predicted vs Observed Orbital Distributions DAVID L. RABINOWITZ Jet Propulsion Laboratory, California Institute of Technology, Mail Stop 183-501, 4800 Oak Grove Drive, Pasadena, California 91109 E-mail: [email protected] Received December 18, 1995; revised December 9, 1996

also the estimated number of main-belt asteroids near the resonances, main-belt asteroids could supply the Earth approachers at a rate sufficient to counter their destruction by planetary collisions or ejection from the Solar System. ¨ pik approximaUsing a Monte Carlo program based on O ¨ tions (Opik 1963, 1976) to model the evolution of mainbelt asteroids and their collisional fragments, Wetherill (1988, hereafter Paper I) predicted a distribution of orbits for kilometer-sized Earth approachers similar to the observed distribution. Wetherill (1991) also showed that encounters with Earth and Venus could perturb Jupitercrossing, short-period comets to orbits similar to kilometersized Earth approachers. However, because most kilometer-sized Earth approachers and main-belt asteroids have similar size-frequencies, reflectance colors, and rotational periods (McFadden et al. 1989; Binzel et al. 1992; Rabinowitz 1996), main-belt asteroids could well be the dominant source. Since the publication of Paper I, there have been major advances in methods of numerical integration which have lead to an improved understanding of the effects of dynamical resonances (Farinella et al. 1994a,b, Gladman et al. 1995, Morbidelli and Moon 1995, Gladman and Burns 1996). The general impact of these studies has been to emphasize the way in which both the n6 and 3 : 1 resonance can accelerate bodies into extremely high-velocity orbits, often on a very short time scale (, 1 myr). The associated high eccentricities lead to orbits with perihelia very near the Sun, and aphelia near Jupiter, either of which is sufficient to remove the body from the inner Solar System. As a consequence, numerically calculated orbits for Earthapproaching asteroids and meteorites stand in marked contrast to the observed steady-state distributions for Earth approachers and to inferred orbits for meteorites. Neither of these results are likely to be incorrect, however. The difference probably results from the relatively small contributions made to the steady state by the rapidly evolving

We predict the orbital distribution of Earth-approaching asteroids with diameter, d, ranging from 10 m to 10 km, assuming they originate as the fragments of main-belt asteroids. The ¨ pik theory to account for the effects of close calculation uses O encounters with the planets and analytic theory to account for the effects of the n6 secular resonance and the 3 : 1 mean motion resonance with Jupiter. We use corrections to the O¨pik formulas to account for long-range encounters with Jupiter, outside the Jovian sphere of influence. We also model the ongoing effects of fragmentation, subsequent to the initial, collisional events that eject parent fragments from the main belt. Biasing our predicted orbital distributions to account for the selection effects of the Spacewatch search, and comparing these biased distributions to the orbits of the Earth approachers detected by Spacewatch, we find that we are able to predict the total number and distribution of orbits for kilometer-sized Earth approachers of moderate inclination, but we do not predict the observed concentration of low-eccentricity orbits among the small Earth approachers (SEAs) with d , p50 m. We therefore conclude that main-belt asteroids are a sufficient source for kilometer-sized Earth approachers, but not for the SEAs. It is plausible that a significant fraction of the SEAs (greater than 5%) are not of main-belt origin.  1997 Academic Press

1. INTRODUCTION

In previous years, various researchers have shown that main-belt asteroids in stable orbits are a probable source for the Earth-approaching asteroids (Williams 1973, Wetherill and Williams 1979, Wisdom 1983, Wetherill 1985, 1988, Greenberg and Nolan 1989, Farinella et al. 1994b). Relatively low-velocity ejecta from those asteroids with orbital motions nearly in resonance with the n6 secular frequency, or with the jovian 3 : 1 commensurability resonance, may evolve from initial orbits inside the main belt to eccentric, Earth-approaching orbits on million-year time scales. Given the estimated efficiency for this process, and 33

0019-1035/97 $25.00 Copyright  1997 by Academic Press All rights of reproduction in any form reserved.

34

DAVID L. RABINOWITZ

orbits, simply as a consequence of their being short lived (for the meteorite orbital distribution, not considered in this paper, the selection effects of atmospheric deceleration and ablation are also of major importance). Because most of the recent numerical work is as yet unpublished, we use a generic method, here, to study the dynamical influences leading to high-eccentricity orbits. This allows us to explore the extent to which these effects can be reconciled with the observed orbital distribution of the Earth approachers. Another purpose of this paper is to address recent observations showing that main-belt asteroids may not be a sufficient source for small Earth approachers (SEAs) with diameter, d 5 5 to 50 m. This population was not observed until the start of the Spacewatch search for Earth-approaching asteroids at the University of Arizona in 1990 (Rabinowitz 1991, Gehrels 1991, Scotti 1994). Using a large-format, charge-coupled device (CCD) and computer software to enhance detection capabilities compared to conventional searches, the Spacewatch search has extended the low end of the detectable size range down from d p 500 m to d p 5 m. The observations reveal a cumulative size distribution for Earth approachers that increases in slope with decreasing d, coincident with an increased fraction of low-eccentricity orbits (Rabinowitz 1993, 1994, Rabinowitz et al. 1993). The colors of the small Earth approachers also show significant deviations compared to the colors of larger Earth approachers and main-belt asteroids (Rabinowitz et al. 1993, Rabinowitz 1996). These observations are not predicted by previous models that assume sources in the main belt, and they may indicate the admixture of inactive cometary fragments, and/or impact ejecta from a local source, such as the Moon, Mars, or Earth. Recently, Bottke et al. (1996) used Monte Carlo calculations to predict the orbits for Earth approachers with d p 50 m, assuming these bodies originate as mainbelt fragments near the 3 : 1 mean motion resonance with Jupiter, and also as impact ejecta from Mars, the Earth– Moon system, Venus, and other local sources. Although their model does not account for the ongoing influence of resonant perturbations, they conclude that the observed SEAs have orbits inconsistent with a main-belt origin, and most similar to ejecta from hypothetical Amor asteroids with Mars-like orbits. They also conclude that the SEAs are unlikely to originate from the stochastic breakup of a single, large Earth approacher. In this paper, we extend the calculation presented in Paper I to cover the size range explored by the Spacewatch observations, d 5 10 m to 10 km. We use results from Bottke et al. (1994, 1996) to account for the dependence of disruption probability on semimajor axis, a, eccentricity, e, and inclination, i. We also determine the size distribution for the Earth approachers as a function of the cumulative size distribution for main-belt asteroids, nmb(d). Our approach is to first develop a model that is consistent with

the Spacewatch observations of larger Earth approachers with d 5 0.1–10.0 km, for which the observational bias is well determined (Rabinowitz 1993, 1994), and the number of observations is large (p70 objects as of June 1995). Agreement between our model and the observations justifies its extension to predict the orbits of the SEAs. The extent to which the predictions and observations differ for SEAs, taking into account both the uncertainties in the model and the effects of observational bias, determines the extent to which sources other than main-belt asteroids are required to supply the population. An important difference between the present work and paper I is that we no longer assume that any asteroidal body injected into the 3 : 1 Kirkwood gap is quickly transferred to a marginally Earth-crossing orbit, in accordance with the results of Wisdom (1983, 1985) and as assumed by Bottke et al. (1996). Also, we now evaluate orbital evolutions in discrete time steps, as is in numerical integrations, rather than using a probabilistic algorithm to evaluate the circumstances only after planetary encounters (Arnold 1965). For marginally planet-crossing bodies, these changes allow us to account for the effects of the phasing of the resonant n6 perturbations with the secular oscillations in the eccentricities of the planet. This is particularly important because it is now known that the close encounters with Earth cannot be expected to efficiently remove bodies from the 3 : 1 resonances before their eccentricities increase to very high values (Farinella et al. 1994a,b). We also recognize that our method for introducing the effects of the resonances, essential to the transfer of asteroidal material to the terrestrial planet regions, is very approximate. Rather than attempting to provide anything like a definitive treatment of this problem, we hope to provide a framework into which more sophisticated algorithms can be introduced as our understanding of these phenomena improves. The present model also allows us to investigate the relationship between the Earth-approacher and main-belt size distributions and to assess the relative contribution of different source regions within the main belt. However, we forestall detailed discussions of these relationships to later papers. Here, we restrict our attention to the case of selfsimilar fragmentation, for which nmb(d) Y d22.5 in the steady state (Dohnanyi 1969, Williams and Wetherill 1994). We also limit our goal to predicting the gross orbital distribution for the Earth approachers as a function of size, with enough detail to determine if the main-belt asteroids are a sufficient source. 2. THE MONTE CARLO CALCULATION

2.1. Overview The Monte Carlo program we use in this paper to determine the steady-state orbital distribution for the Earth approachers is an extension of similar calculations per-

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

formed by Wetherill (1979, 1985, Paper I) and others (Bottke et al. 1996), which are based on the original work ¨ pik (1963). The fundamental asof Arnold (1965) and O sumptions are as follows. Main-belt asteroids with stable orbits experience collisions with other main-belt asteroids that liberate fragments. The collisions impart some of these fragments with velocities that move them into unstable orbits. Owing to close encounters with planets and also to long-range, gravitational perturbations, some of these unstable bodies migrate to orbits approaching Earth. While the orbits for these Earth approachers continue to intersect the main belt, they are subject to ongoing collisions. Each primary body thereby generates a ‘‘retinue’’ of descendant fragments, owing to a sequence of secondary (tertiary, etc.) collisions. These descendant fragments replace the primary bodies on their Earth-approaching migrations. Eventually, all of these Earth-approaching bodies and their progeny are ejected from the Solar System, collide with a planet or the Sun, or are destroyed by collisions with other asteroids. The steady-state population of Earth approachers is composed of all the surviving fragments, generated over the age of the Solar System (4.5 byr) that currently have Earth-approaching orbits. A calculation to determine the orbital distribution of this population must therefore account for the ongoing effects of (i) the collisional ejection of fragments from main-belt asteroids, (ii) their subsequent dynamical evolution, and (iii) their ongoing collisional fragmentation. As in Paper I and similar calculations, the approach we take is to envision the Earth approachers as an ensemble of bodies resulting from the collisional ejection at regular intervals in the past (1 myr) of a fixed number of mainbelt fragments. The present-day orbital distribution of the Earth approachers is given by the survivors of those fragments injected 1 myr ago, 2 myr ago, and so on back to the formation of the Solar System 4.5 byr ago. However, we calculate this distribution by modelling not the past effects, but the future effects of processes (i)–(iii), described above, for a large number of main-belt fragments, and we examine the residual population at 1 myr intervals, or ‘‘episodes’’. The contribution of 1-myr-old bodies is represented by the residual population after 1 myr, 2-myrold bodies by those surviving 2 myr, etc. To further simplify this calculation, we account only for the effects of catastrophic collisions subsequent to the initial event that injects each parent fragment from the main belt. Catastrophic events completely disrupt a body and replace it with a family of smaller fragments, each with a mass less than half the original. We ignore the effects of cratering events, which incrementally reduce the mass of a body, because such events would require too much computer time to model. However, because the disruptive effects of cratering and catastrophic collisions are comparable for a self-similar population (Wetherill 1967, Farinella

35

et al. 1994b), we do not expect the additional effects of cratering to alter our predictions. It is also important to note that this calculation makes no assumptions regarding the nature of the collisional events that initially inject the parent bodies from the main belt. Either catastrophic collisions or cratering events, or a combination of the two processes, could be responsible, but it is of no importance in this calculation. As a final simplification to the calculation, we assume that each catastrophic collision introduces only a negligible divergence in the orbits of the descendant fragments. Each episode in the orbital evolution of a single parent body determines the mean orbit for an entire lineage of descendants. We then assign random orbits with small scatter about the mean when we plot the orbits. Of course, this assumption is not realistic. Owing to the chaotic evolution of planet-crossing orbits, small perturbations lead to widely diverging orbits. However, as long as the number of parent bodies that we model is large enough so that the resulting orbital distribution completely covers the range of orbits possible for Earth approachers, it will accurately represent the steady state. 2.2. Details To implement the calculation summarized above, we use three, sequential computer programs, ‘‘INJECT,’’ ‘‘EVOLVE,’’ and ‘‘FRAGMENT,’’ which we describe in detail in the following subsections, and also in Appendices A and B. First, INJECT models the ejection of parent fragments from the main belt owing to collisional events. Then EVOLVE models the dynamical evolution of these parent fragments, but in the absence of further collisional evolution. As discussed above, we allow each episode in the orbital evolution of a given parent fragment to represent a separate body in the steady state. Finally, FRAGMENT models the collisional evolution for each of these separate bodies. Here, we are able to separate the calculations performed by FRAGMENT and EVOLVE because we allow the orbit for each body in the steady state, as calculated by EVOLVE, to represent the mean orbit of its collisional descendants, as generated by FRAGMENT. Furthermore, in order to compare our predicted results with the observations of the Spacewatch Telescope, we use an additional program, ‘‘OBSERVE,’’ to model the selection effects of the Spacewatch search. OBSERVE biases the steady-state orbits we predict in the same way that the observations of the Spacewatch Telescope are biased. 2.2.1. INJECT. Starting with an orbital distribution, Nmb(a, e, i), characterizing the main-belt asteroids, INJECT randomly chooses orbits for parent, main-belt bodies. For each parent orbit, INJECT then computes the altered orbit for a fragment ejected at a random position along the parent orbit, in a random direction relative to

36

DAVID L. RABINOWITZ

the parent body and at fixed ejection velocity, n0 . The size of the fragment is irrelevant since we are interested only in the orbital distribution as a function of n0 . Those fragments that are potential planet crossers, owing to the effects of the n6 resonance or the 3 : 1 mean-motion resonance with Jupiter (see Appendix A) determine an initial ensemble of injected orbits, Iinj . We ignore those fragments that do not meet these criteria, as they do not contribute to the Earth-approaching population. Consistent with the results of Zappala et al. (1984), who analyze the dispersion in the orbital elements of the members of various asteroid families, we adopt a nominal value of n0 5 100 m/sec, with a factor of 2 uncertainty. 2.2.2. EVOLVE. For each member of Iinj , EVOLVE computes an orbital evolution, as described in Appendix A. The program takes into account only dynamical effects and ignores fragmentation. The sizes of the evolving bodies remain irrelevant. Hence, EVOLVE generates an ensemble of state vectors, Idyn 5 huj, where each vector, u, represents orbital elements a, e, and i, and also age, Dt, at 1 myr intervals in the orbital evolution of every member of Iinj . Interpreting each vector to represent a separate body, the resulting ensemble represents the steady-state distribution in a, e, i, and Dt of the present-day Earth approachers, but only when the effects of fragmentation are negligible. To prepare for the calculation to be performed by FRAGMENT, EVOLVE performs one additional function. For each body with age, Dt, EVOLVE computes a new dynamical variable, Dt9, which we call the ‘‘effective main-belt age.’’ For a given Earth approacher, which has acquired a given, cumulative probability of catastrophic disruption since the time of its injection from the main belt, Dt9 is the average interval of time that same body would have to spend on a stable, main-belt orbit in order to accrue the same cumulative probability of disruption. For most Earth approachers, which spend most of their dynamical lifetime on orbits intersecting the main belt, Dt9 . Dt because they collide with main-belt asteroids at higher velocities than main-belt asteroids collide with each other. This exposes the Earth approachers to more frequent disruption by smaller main-belt bodies, which are more numerous. Hence, main-belt asteroids have longer collisional lives than most Earth approachers of the same diameter (Bottke et al. 1994). Section 2.3 describes the method EVOLVE uses to calculate Dt9. 2.2.3. FRAGMENT. For each body represented by Idyn (where we assume that each orbit represents a different body, as discussed above), FRAGMENT assigns an initial diameter, d0 , corresponding to the diameter at the time of injection from the main belt (we assume a maximum value, d0 , 10 km, corresponding to the approximate diameter of the largest observed Earth approachers). The program randomly chooses d0 , weighted by an assumed production

function, ­nprod(d)/­t, describing the cumulative production rate of main-belt asteroids larger than diameter, d. Assuming the main belt is in collisional equilibrium, it follows (Wetherill 1985) that ­nprod(d)/­t is given by ­2nprod(d) 1 ­nmb(d) 5 , ­d ­t tmb(d) ­d

(1)

where tmb(d) is the average collisional lifetime of a mainbelt asteroid with diameter, d. Given our assumptions that nmb(d) Y d22.5 and that fragmentation is self-similar, it follows that tmb(d) Y d 0.5 (see Section 2.3). Hence, Eq. (1) yields ­nprod(d)/­t Y d23.0. After assigning d0 , FRAGMENT goes on to model the effects of randomly occurring, disruptive collisions. The program models a different fate for each body represented by Idyn , but there are three general cases: (a) The parent body suffers no disruptive collisions, and the end result is one body with d 5 d0 . (b) The parent body suffers a disruptive collision and generates a retinue of fragments. These descendant fragments may suffer additional disruptive collisions. The end result is a family of bodies with a range of values for d. (c) The parent body suffers a disruptive collision and generates a retinue of fragments, but all the descendant fragments are too small to count (d , dmin 5 10 m). The end result is a null set. FRAGMENT models fates (a)–(c) as follows. Given the values of d0 , Dt9, and tmb(d0) associated with each member of Idyn , the probability of surviving without disruption to effective age, Dt9, is P 5 e2[Dt9/t mb(d0)]. FRAGMENT therefore randomly chooses fate (a) with weight, P. For each body not receiving fate (a), the program must then model the evolution of a collisional retinue. The first step is to choose a random time since injection, t0 , for an initial disruption to occur. Since the parent body does not survive to effective age Dt9, FRAGMENT must choose t0 in the interval t 5 0 to Dt9 and weight this choice by the cumulative probability of disruption, 1 2 e2[t0 / tmb(d0)], in the interval t 5 0 to t0 . The program models the initial disruption by replacing the parent body with a first-generation family of descendant fragments. For each family, the cumulative number of fragments, nfrag(d0 , d), larger than diameter, d, is a fixed function of d and d0 (discussed below). All descendants share the same orbit as their parent. FRAGMENT then goes on to model the subsequent disruption of each first-generation fragment with d . dmin , and also every second-generation fragment with d . dmin , and so on, until it has modeled the disruption of all latergeneration fragments with d . dmin . This in an inherently recursive procedure, which we have modeled accordingly (see Appendix B). Either of fates (a)–(c) are possible for

37

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

each generation of fragments, and they all share the same orbit as their parent. The superset of bodies, composed of all the undisrupted members of Idyn plus all the descendants with d . dmin that replace the disrupted parents, represents the steady-state size and orbital distribution of the Earth approachers, Iea . This distribution accounts for the effects of both dynamical and collisional evolution. To fix nfrag(d, d9), we again refer to the study of asteroid families by Zappala et al. (1984). They show that the cumulative distribution of masses, n(m), for families with no members larger than 100 km can be described by n(m) 5 [m/m2]12k, where m2 is the mass of the largest fragment, and k is in the range 5/3 to 11/6. These results are generally consistent with laboratory experiments (Fujiwara et al. 1989, Martelli et al. 1994). The equivalent diameter distribution is nfrag(d, d9) 5 [d9/d2(d)]2e, where d2(d) is the size of the largest fragment, and e 5 3(k 2 1). Furthermore, the requirement that the total mass of the fragments equals the mass of the original body implies d2(d) 5 d [(3 2 e)/ e]1/3. We adopt these relations, taking e 5 2.4 as a nominal specification. Note that these relations for nfrag(d, d9) and d2(d) are scale-invariant, which is consistent with our assumption of self-similar fragmentation. 2.2.4. OBSERVE. Starting with Iea , OBSERVE simulates the selection effects of the Spacewatch search, which are determined by the limiting magnitude, limiting angular rates, and other detection criteria of the Spacewatch Telescope. OBSERVE is identical to program ‘‘ESP,’’ described in detail by Rabinowitz (1993, 1994). The program assigns random mean anomalies, longitudes of ascending node, and arguments of perihelion to each member of Iea , so that it can then calculate the apparent motion, visual magnitude, and position in the sky for each body. The ensemble of bodies, Ibias , that would be detectable by the Spacewatch search based on the known detection criteria represents the orbital distribution of the Earth approachers modulated by the effects of detection bias.

belt asteroid capable of disrupting a body of diameter, d, at the mean velocity of collisions, nmb , between asteroids in the main belt. The upper limit, dmax , is the diameter of the largest main-belt asteroid. Parameter Pmb is a geometric quantity that depends only upon the orbital distribution of the main-belt asteroids. For a body with a given orbit, the intrinsic probability of collision, Pi , while passing through the main belt is the probability of collision per unit cross-section of the body when the total number of main-belt asteroids is normalized to one. Hence, Pmb is the average value of Pi for main-belt asteroids. In this paper, we assume self-similar (i.e., scale-invariant) fragmentation. Hence, d1(d) Y d. Given the observation that d1 ! d for hypervelocity impacts (Holsapple 1994) and also our assumption that nmb(d) Y d22.5, Eq. (2) yields tmb(d) Y d 0.5 for d ! dmax . We normalize this relation by fixing tmb(d) 5 3.8 3 10 8 yr at d 5 1.0 km. As shown by Paper I, this normalization is consistent with the collisional lifetimes of meteorites (assuming the self-similar solution for tmb(d) extends to d p 10 cm). This normalization is also consistent to within a factor of 2.5 with recent estimates by Bottke et al. (1994). Following the work of Bottke et al. (1994), we now use Eq. (2) to determine the probability of disruption, tea21(d, a, e, i), of a given Earth approacher with diameter, d, and orbital elements a, e, and i, by collision with a main-belt asteroid. To do this, we replace Pmb with the intrinsic probability of collision, Pi (a, e, i), between an Earth approacher with orbital elements a, e, and i, and a main-belt asteroid. We also replace nmb with the mean impact velocity, nea (a, e, i), for these collisions. Both Pi (a, e, i) and nea (a, e, i) have been mapped by Bottke et al. (1994), given the orbital distribution of the known main-belt asteroids larger than 50 km. Making these replacements, Eq. (2) yields dmax

E

P (a, e, i) 1 5 i tea(d) 4 d95d (d,n 1

(d 1 d9)2

ea (a,e,i))

2.3. Effective Main-Belt Age, Dt9

(3)

As discussed in the previous section, Dt9 governs the collisional evolution of each member of the dynamical ensemble, Idyn . Here, we show how Dt9 depends on the orbital history of each body. We start with an expression given by Wetherill (1967) that relates the average probability of catastrophic collision in the main belt, tmb21(d), to the average, intrinsic probability of collision between mainbelt asteroids, Pmb , and also to nmb(d): dmax

E

1 P 5 mb tmb(d) 4 d95d (d,n 1

mb)

(d 1 d9)2

­nmb(d9) ­d9. ­d9

­nmb(d9) ­d9. ­d9

Equations (2) and (3) now allow us to form an expression for Dt9. Recall that we define Dt9 to be the age at which the average cumulative probability of disruption for a main-belt asteroid with diameter, d, matches the cumulative probability of disruption of a given Earth approacher of the same diameter, but age, Dt. Hence, Dt

E

Dt9 ­t 5 . tmb(d) t50 tea (d, a, e, i)

(4)

(2)

Here, d1(d, nmb), is the diameter of the smallest main-

Here, the integral over time from t 5 0 to Dt takes into account the time-dependence of tea21(d, a, e, i), owing to the time-dependence of the Earth approacher’s orbit.

38

DAVID L. RABINOWITZ

In order to solve Eq. (4) explicitly for Dt9, we assume that the criterion that determines d1(d, n) is a threshold value for the impactor’s kinetic energy relative to the strength of the target body. This assumption is the basis for most theories of catastrophic breakup. It leads to the conclusion that d1(d, n) Y d1(d)n22/3, with d1(d) determined by the size dependence of target strength (Davis et al. 1989). Although recent laboratory investigations indicate that this assumption may fail when the self-gravity of the target body is important (Holsapple 1994), this result does not apply to bodies in the diameter range of Earth approachers (d , 10 km). Given the relation d1(d, n) Y d1(d)n22/3, and also the observations that d1 ! d and dmax @ d, Eqs. (2)–(4) combine to yield Dt9 P

1 5/3 Pmbnmb

Dt

E P (a, e, i)[n i

ea(a,

e, i)] 5/3 ­t.

(5)

t5 0

Note that this expression yields Dt9 p Dt for any body with a stable, main-belt orbit because Pi (a, e, i) p Pmb and nea(a, e, i) p nmb for all t. For most Earth approachers Dt9 . Dt because their orbits intersect the main belt for most of their lives, and because nea(a, e, i) . nmb and Pi(a, e, i) . Pmb (Bottke et al. 1996). Also note that Eq. (5) is independent of d1(d) and dmax . We do not need to parameterize the target strength or the threshold kinetic energy for catastrophic disruption in order to calculate Dt9. To normalize Eq. (5), we assume Pmb 5 2.86 3 10218 km22 yr21 and nmb 5 5.3 km/sec, as calculated by Bottke et al. (1994). 2.4. The Uncertainties for Nmb(a, e, i) and nmb(d) Uncertainties for both nmb(d) and Nmb(a, e, i) arise because observations of main-belt asteroids are complete only for d . p30 km (Cellino et al. 1991). For the range of sizes relevant to this paper (d 5 10 m to 10 km), only a small fraction of the population is known. These observations are also biased toward those bodies that are nearest at opposition (i.e., small a). If it were safe to assume that Nmb(a, e, i) were independent of size, or equivalently, that nmb(d) were independent of a, e, and i, then there would be little uncertainty in either distribution. Both functions would be well determined by the known orbits and diameters of p1000 main-belt asteroids with d . p30 km (absolute magnitude, H , 10.4). Unfortunately, previous efforts to determine nmb(d) by observations of main-belt asteroids reveal a size distribution that varies widely across the main belt and also for different ranges of d. These results are summarized by Fig. 1. Fitting nmb(d) to a power law, d p21, where p is a constant, Fig. 1 shows 1 2 p vs a for different ranges of d as determined by the Palomar–Leiden survey, IRAS observations of main-belt asteroids, and observations by the Galileo spacecraft of small craters on (951) Gaspra.

FIG. 1. Observed values for the power-law exponent, 1 2 p, assuming nmb(d) Y d p21 for different ranges in a. The asterisk symbol shows the results from the Galileo observation of (951) Gaspra, which constrain p in the range d 5 10–100 m (Belton et al. 1992). Octagons show the results of the Palomar–Leiden survey, which constrain p in the range d 5 1–20 km (van Houten et al. 1970). Squares and triangles show the results from measurements by the IRAS satellite of asteroid albedos, which constrain p in the range d 5 30–100 km and d 5 100–300 km, respectively (Cellino et al. 1991).

Given the wide range of values for 1 2 p, it is apparent that nmb(d) depends on both a and d. It is therefore unlikely that Nmb(a, e, i) is size independent throughout the main belt. Nonetheless, for the range of diameters covered by the observed Earth approachers and for the region of the main belt that is likely to supply most of these bodies (a , 2.6 AU; see Paper I and also Sections 4.1 and 5.1), the observations are consistent with a uniform power law with exponent 1 2 p 5 2.5 6 0.5. This paper therefore uses the orbits of the known main-belt asteroids with d . p30 km to represent Nmb(a, e, i) in the range d 5 0.01 to 10.0 km. The uncertainty in this extrapolation grows with decreasing d in proportion to the uncertainty in nmb(d). Since the uncertainty in 1 2 p is 60.5 for a , 2.6 AU, the

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

39

resulting uncertainty for Nmb(a, e, i) is the factor (d/30 km)20.5. We emphasize, here, that we are not attempting to identify the particular .30 km asteroids that define Nmb(a, e, i) as actual sources for the Earth approachers, in the same way that (4) Vesta is thought to be the source of the basaltic achondrites (Binzel and Xu 1993). Rather, we use these large asteroids that have been sampled completely as indicators of how asteroids bodies are distributed in the main belt. For the size distribution that we have assumed, we can expect that most of the Earth approachers and meteorites will be fragments of much smaller bodies, and that these smaller sources will have orbits distributed in the same was as the .30 km bodies used in the calculation. 3. COMPARISONS WITH NUMERICAL EVOLUTION

Before implementing the complete Monte Carlo calculation described in Section 2, two comparisons were made to determine the validity of program EVOLVE. In Comparison I, EVOLVE computed the evolution of a few known asteroids for which numerically integrated evolutions have been published by Farinella et al. (1994a,b). These are special case asteroids predicted to collide with the sun on a time scale of p106 yr, owing to the effects of secular and mean motion resonances. Starting with the known orbit for each of these test asteroids, EVOLVE computed a series of 10 random evolutions. We then compared each of these predictions to the corresponding numerical integrations. In Comparison II, EVOLVE computed 10 random evolutions for each of 47 known asteroids for which Farinella et al. (1994a) computed both backward and forward integrations. Again, these were special objects chosen by Farinella et al. because of their propensity for solar collisions. Of the 94 orbits integrated by Farinella et al. and of the 470 orbital evolutions calculated by EVOLVE, a comparison was made of the fraction that ended by collision with the sun or by ejection on a hyperbolic orbit within 106 yr. 3.1. Results of Comparison I Figure 2 shows the time evolution of a, e, i, and the longitude of perihelion with respect to Saturn, g 2 gs for asteroid (3551) Verenia. For this asteroid, the n6 resonance is the dominant long-range perturbation. The effects of the resonance are strongest when g 2 gs is slowly varying. The right side of Fig. 2 shows the numerical solution determined by Farinella et al., which includes both backward and forward integrations. The left side of Fig. 2 shows both the first and the sixth evolutions determined by EVOLVE, with the first evolution plotted backward in order to compare with the backward numerical integration. Because of the chaotic nature of both calculations, and of the real orbital evolution of this body, it should not be expected

FIG. 2. Semimajor axis, a, inclination, i, eccentricity, e, and longitude of perihelion with respect to Saturn, g 2 gs , vs time, t, for asteroid (3551) Verenia. The left and right sides of the figure, respectively, show the predictions of program EVOLVE and of Farinella et al. (1994b).

that the comparisons will exactly match one another. What will be seen, however, is that similar patterns of evolution are found in both cases. The numerical integrations show more detail because our Monte Carlo calculations do not model the finer aspects of the planetary perturbations. Comparing the left and right sides of Fig. 2, it is apparent that EVOLVE calculates variations in a, e, i, and g 2 gs that are similar to the numerically predicted variations. In the backward evolutions for Verenia, both EVOLVE and Farinella et al. predict a comparable range of variation for a and also for i. EVOLVE also predicts an increase for e from 0.5 to 0.8 while g 2 gs is resonant (t 5 0 to 20.2 3 106 yr). This is followed by a period of roughly constant e, while g 2 gs is out of resonance (t 5 20.2 3 106 to 20.8 3 106 yr). Finally, e increases from 0.8 to p1.0 (collision with the sun) when g 2 gs returns to resonance (t 5 20.8 3 106 to 21.8 3 106 yr). Farinella et al. predict a similar backward evolution for e and g 2 gs , although the intervals when g 2 gs moves out and then back into resonance are shorter. Hence, Farinella et al. predict collision with the sun at t 5 20.6 3 106 yr. For the forward evolution of Verenia, EVOLVE predicts constant values for a and i, and a slowly varying value for g 2 gs. Meanwhile, e descends from 0.5 to 0.1 and then quickly rises up to 1.0 for a solar collision at t 5 0.8 3 106 yr. Farinella et al. predict a similar forward evolution until t 5 0.8 3 106 yr, when a close encounter with Earth kicks g 2 gs out of

40

DAVID L. RABINOWITZ

each of 47 known asteroids. Panel a shows the number of evolutions, N, that ended either by solar collision or by hyperbolic ejection vs the time, Dt, at the end of the evolution. The solid black distribution represents the number, Ns , that ended by solar collision only. Fig. 3b shows N vs Dt and Ns vs Dt determined by Farinella et al. by forward and backward numerical integrations of the same 47 asteroids. Both EVOLVE and Farinella et al. predict comparable values for the total fraction of bodies surviving less than 106 yr (22 6 2 and 31 6 6%, respectively) and for the relative fraction of bodies going into the sun (50 6 7 and 70 6 20%, respectively). Both predictions also yield distributions for N vs Dt that are peaked near Dt 5 106 yr, and fall off to N p 0 for Dt , 104.5. 3.3. Comparison Implications Comparisons I and II show that EVOLVE accounts for the gross effects of the n6 and 3 : 1 resonances in a reasonable way. For those asteroids that go into the sun on times scales of p106 yr, EVOLVE often predicts the same outcome as numerical integrations. Comparison II shows that EVOLVE also predicts the correct distribution of dynamical lifetimes for Dt , 106 yr, and a realistic fraction of bodies going into the sun. EVOLVE is therefore capable of modelling the dynamical effects of both close planetary encounters and long-range resonant forces in a manner that mimics the influences revealed by numerical integration. The Monte Carlo calculation is more useful, in the present application, because methods of numerical integration are not yet fast enough to perform the required number of evolutions in a manageable amount of time. FIG. 3. The number of orbital evolutions, N, vs the time, Dt, of ejection from the Solar System or collision with the Sun for asteroids with orbits integrated by Farinella et al. (1994a): (a) the predictions of program EVOLVE; and (b) the results of the numerical integration by Farinella et al. The shaded histograms show the number of orbital evolutions, Ns , ending only by solar collision. The vertical axes on the right show percentage of computed orbits.

resonance. Hence, e remains constant beyond that point, and no collision with the sun occurs. In addition to Verenia, we used EVOLVE to compute 10 evolutions for 1988 VP4, 1992 SZ (misidentified as 1991 SZ in Farinella et al. 1994a), and also for a fictitious fragment of (6) Hebe (from Fig. 2 of Farinella et al. 1994b). For the Hebe fragment, the 3 : 1 mean motion resonance is the dominant perturbation. In each of these cases, we were also able to find at least one evolution similar to the numerical evolution. The qualitative agreement was comparable to the comparison in Fig. 2. 3.2. Results of Comparison II Figure 3 shows the results of Comparison II, in which EVOLVE computed 10 evolutions up to t 5 106 yr for

4. MONTE CARLO RESULTS

The Monte Carlo calculation described in Section 2 was performed given the nominal values for Nmb(a, e, i), nmb(d), tmb(d), n0 , and nfrag(d, d9) summarized by Table I, and also given the model parameters described in Appendix A. The following is a description of the results from each phase in the calculation, corresponding to the output of programs INJECT, EVOLVE, FRAGMENT, and OBSERVE. TABLE I Nominal Specification and Range of Uncertainty for Critical Parameters Parameter

Nominal Specification

Uncertainty

nmb(d) tmb(d) Nmb(a, e, i) n0 nfrag(d, d9)

1069 3 (d/30 km)22.5 [3.9 3 108 yr21](d/1 km)0.5 known orbits for d . 30 km 100 m/sec (d9/d2)22.4

(d/30 km)20.5 for d , 30 km factor of 2 at d 5 1 km (d/30 km)20.5 for d , 30 km factor of 2 (d9/d2)22.0 to (d9/d2)22.4

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

4.1. Injection Results Figures 4a and 4b, respectively, show sin(i) vs a and e vs a for ensemble Imb , the known asteroids with d . p30 km (H , 10.4). Here we have used the proper elements calculated by Milani and Knezevic (1994), or else the elements calculated by Lemaitre and Morbidelli (1994), to initialize a, e, and i (see Appendix A). Given this distribution, program INJECT determined the orbits of the injected ensemble, Iinj , for which Figs. 4c and 4d show sin(i) vs a and e vs a, respectively. The total number of injected

FIG. 4. Orbital distribution, Imb , for known asteroids with d . p30 km (a and b) and the orbital distribution, Iinj , of injected fragments (c and d). Dashed lines in (a) and (c) show the assumed location for the n6 resonance for proper eccentricity eo 5 0.1, 0.2, and 0.3 (bottom, middle, and top curves, respectively) and also the assumed range of activity for the 3:1 mean motion resonance with Jupiter. The dashed curve in (b) and (d) shows q 5 Qmars , where Qmars 5 1.7036 AU, the maximum aphelion for Mars. In (a) and (b), dark squares show potential sources that contribute to Iinj , crosses show sources rejected because they are planet crossers (qmin , Qmars), dark squares marked by larger, concentric squares show sources for which eo has been reduced to keep qmin . Qmars , and small dots show potential source bodies too far from the n6 and 3 : 1 resonances to contribute to Iinj .

41

orbits is 1000. The three dashed curves in Figs. 4a and 4c show the assumed location of the n6 resonance for e0 5 0.1, 0.2, and 0.3. A gap appears in each curve in the range a 5 2.48 to 2.52 AU, which is the zone of activity we assume for the 3 : 1 mean motion resonance with Jupiter. Of the 1000 members of Iinj , 599 are outside the 3 : 1 resonance, 401 are inside. The dashed curve in Figs. 4b and 4d shows q 5 Qmars , where Qmars is the maximum aphelion for Mars owing to secular perturbations (1.7036 AU, see Appendix A). As discussed in Section 2.2.1 and Appendix A, Iinj includes only those fragments that are potential planet crossers owing to the orbit change introduced by n0 and to the perturbations of either the n6 resonance or the 3 : 1 meanmotion resonance with Jupiter. These selection criteria restrict the contributing source bodies to the subset represented by dark squares in Figs. 4a and 4b. Those sources represented by small dots are too far from the n6 and 3 : 1 resonances, or too far from the curve marking q 5 Qmars , to produce planet-crossing fragments. On the other hand, those members of Imb represented by crosses are rejected as sources because they are too close to the resonances, or too close to Mars crossing, to be counted as stable orbits for 4.5 byr. They are asteroids (234) Barbara, (329) Svea, (391) Ingeborg, and (1036) Ganymede. We find, however, that including Barbara and Svea has little effect on the outcome of the calculation, as they are long-lived Mars crossers that are only marginally unstable. Ganymede and Ingeborg, on the other hand are deep Mars crossers with dynamical lifetimes that are too short (!1 byr) to serve as steady-state sources. They are probably large fragments ejected relatively recently from stable main-belt asteroids. It is also important to recognize that our assumption that all large parent bodies (d . p30 km) are on stable orbits does not imply stable orbits for their ejecta. In our model, secondary fragments with d , 10 km are continuously injected onto unstable orbits. These secondary bodies also supply the population of smaller Earth approachers by continued fragmentation. Figures 4a and 4b also show eight other sources (represented by dark squares inside larger, concentric squares) that are barely Mars crossing. Given their mean elements as calculated by the above mentioned authors, our model yields a minimum perihelion, qmin , Qmars , for each body. However, if we assign only slightly lower values for e0 (by p0.005), then qmin . Qmars . We have assumed that the sharp drop-off in the number of asteroids of all sizes represents a stability boundary with respect to Mars crossing at minimum asteroid perihelion and maximum Mars aphelion. If so, it is likely that these bodies really are stable and that our model for the stability boundary is not precise enough to indicate otherwise. We have therefore made the minimum adjustments to e0 necessary to have qmin . Qmars for each body.

42

DAVID L. RABINOWITZ

We also considered the alternate assumption that our stability criteria are correct. Instead, the large asteroids (d . 30 km) that we observe today with qmin , Qmars may be the remnants of a primordial, unstable population, most members of which have been lost. We examined this assumption by repeating our Monte Carlo calculations, and leaving out our alterations to the mean elements of the parent bodies. As a result, we find that the fragments of the parents with qmin , Qmars dominate the steady-state population of the Earth approachers, and yield a heavy concentration of Amors with a , 2.0 AU. Accounting for observational selection, and then comparing to the distribution observed by Spacewatch (see Section 5.4), we also find that this heavy concentration is not observed. If this interpretation is correct, then this represents a discrepancy between theory and observation. In order to judge the sensitivity of Iinj to n0 , we repeated the selection of Iinj by running program INJECT with n0 5 50 and 150 m/sec, the range of uncertainty identified in Section 2.2.1. We do not plot the results because they are similar to the results already plotted in Fig. 4. The primary effects of increasing n0 are to increase the fraction, Finj , of source bodies from Imb that contribute to Iinj and to increase the relative fraction of contributing sources, K, that inject fragments into the 3 : 1 resonance. For n0 5 50, 100, and 150 m/sec, Finj 5 2.5 6 0.5, 4.3 6 0.6, and 6.6 6 0.8%, and K 5 26 6 10, 48 6 10, and 52 6 9%, respectively. 4.2. Evolution Results Figures 5a and 5b show q vs a for the set of orbits, Idyn , generated by program EVOLVE given Iinj as starting orbits. The program computed 2000 evolutions, two for each of the 1000 members of Iinj . Recording each evolving orbit every 106 yr produced a total number of orbits, Ndyn 5 5553, with q , 1.3 AU. Those members of Iinj initially perturbed by the n6 secular resonance (i.e., outside the 3 : 1 jovian resonance) produced the 3746 orbits plotted in Fig. 5a, while those members inside the 3:1 resonance produced the 1807 orbits plotted Fig. 5b. Note that not all of the 2000 evolutions determined by EVOLVE are represented by Figs. 5a or 5b because some bodies never become Earth approachers or become Earth approachers for a short time compared to the recording interval of 106 yr. Hence, these bodies are ejected from the Solar System, collide with a planet, or go into the sun before EVOLVE has a chance to record their Earth-approaching orbits. Of all 2000 evolutions, 43.1% end as hyperbolic orbits, 8.1% end as long-period (.1000 yr) orbits, 3.8% end by collision with a planet, and 19.2% go into the sun. These respective outcomes occur at median ages, kDt9l 5 5.0, 5.0, 200, and 100 myr (note, however, that the age distributions are bimodal: 48% of the bodies going into the sun do so in

FIG. 5. Perihelion, q, versus semimajor axis, a, for steady-state ensembles predicted by EVOLVE (Idyn) and FRAGMENT (Iea): (a) slowtrack members of Idyn , originating outside the 3 : 1 mean motion resonance with Jupiter; (b) fast-track members of Idyn , originating inside the 3 : 1 resonance; (c) members of Iea with d 5 1.0–1.26 km; and (d) members of Iea with d 5 10.0–12.6 m.

less than 3.5 myr). The relative fraction of bodies colliding with planets Mercury, Venus, Earth, Mars, and Jupiter are 2.6, 28.9, 32.6, 6.6, and 28.9%, respectively. No impacts with the planets more distant than Jupiter occurred. Figure 6a shows n(Dt), the cumulative fraction of orbits from Idyn with ages less than Dt. The distribution ranges from 107 to 109 yr, in rough agreement with the dynamical ages determined by Wetherill in Paper I. Figure 6a also shows (dashed lines) the cumulative distribution of effective main-belt ages, n(Dt9). These ages are shifted to higher values relative n(Dt) for the reasons discussed in Sections 2.2.2 and 2.3. Also note that n(Dt) and n(Dt9) have inflection points near Dt 5 2.0 3 108 and Dt9 5 6.0 3 108 yr, respectively. This inflection occurs because bodies originating from sources inside and outside the 3 : 1 resonance achieve Earth-approaching orbits on different time scales. Figure

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

43

4.4. Observation Results Figures 7a and 7c, respectively, show q vs a and sin(i) vs a for the biased distribution of orbits, Ibias , determined by program OBSERVE for bodies with d 5 0.050 2 10.0 km given the predictions for Iea plotted in Fig. 5c. Figures 8a and 8c show similar predictions for bodies with d 5 0.005 2 0.050 km, given the orbits plotted in Fig. 5d. Here, OBSERVE assumes the same detection criteria described in Rabinowitz (1994), except that the search area with the upgraded CCD (Tek2Ke) is extended to 2673 deg2 to account for further observations in the period October 1994 to April 1995. Each point in Figs. 7a, 7c, 8a, and 8c has been randomly selected from Iea , after weighting each member of Iea by its relative probability of detection. Also, random jitter has been added to each point (as in Figs. 5c and 5d) to reveal multiple representation of single orbits. The total number of points in each figure has been selected

FIG. 6. Cumulative age distributions for Idyn : (a) for all members, and (b) separately for fast and slow-track members. Solids lines show n(Dt). Dashed lines show n(Dt9). See Section 4.2 for details.

6b shows n(Dt) and n(Dt9) for each of these separate populations. The respective values for kDtlm and kDt9lm are 1.3 3 107 and 3.2 3 107 yr inside the 3 : 1 resonance, and 1.3 3 109 and 1.6 3 109 yr outside the resonance. 4.3. Fragmentation Results Figures 5c and 5d show the expected orbits, Iea , for Earth approachers with d 5 1.0–1.26 km and d 5 10.0–12.6 m, respectively. These distributions were generated by program FRAGMENT, starting with Idyn as represented in Figs. 5a and 5b. Here, we have chosen an equal, but arbitrary, number of points to represent the distribution for each d interval. We do this by sampling the orbits at random, after weighting them in proportion to their relative frequency. To make evident the relative frequency of overlapping orbits, we have also added random jitter to each point with rms scatter Da 5 0.05 AU an Dq 5 0.02 AU (see Section 2.1). Note that Figs. 5a–5d include the predicted orbits for bodies with aphelion Q , 0.987 AU (Earth’s perihelion), even though these are not normally defined as Earth approachers.

FIG. 7. Perihelion, q, and sin of inclination, i, vs semimajor axis, a, for Earth approachers with H , 24 (d . p50 m): (a) and (c) show predicted orbits with observational bias simulated by program OBSERVE; (b) and (d) show orbits observed by Spacewatch from September 1990 to April 1995.

44

DAVID L. RABINOWITZ

arbitrarily. If the assumptions of the Monte Carlo calculation are realistic, then these biased predictions should match the observed distributions which are plotted to the right of each prediction (Figs. 7b, 7d, 8b, and 8d). The observations cover the range 15 , H , 24 (p50 m , d , p5.0 km) in Fig. 7 and 24 # H , 30 (p5.0 m , d p50 m) in Fig. 8, and were recorded by Spacewatch from September 1990 to April 1995. To further compare Ibias to the observed distributions of orbits, Fig. 9a shows the number of large Earth approachers (d . p50 m) with q . 1.017 AU (Amors) vs a, and Fig. 9b shows the number of small Earth approachers (d , p50 m) with q . 0.9 AU vs a. Both figures show the observed and predicted distributions, where we have separately normalized Ibias in each figure so that the total number of all Earth approachers (Amors, Apollos, and Atens)

FIG. 9. The numbers of Earth approachers, N, vs semimajor axis, a, observed by Spacewatch (solid line) and predicted by OBSERVE (dashed line) for (a) d . p50 m and (b) d , p50 m (see Section 4.4 for details).

matches the observed number (70 for d . p50 m and 28 for d , p50 m). 5. DISCUSSION

5.1. Injection

FIG. 8. Perihelion, q, and sin of inclination, i, vs semimajor axis, a, for Earth approachers with H . 24 (d , p50 m): (a) and (c) show predicted orbits with observational bias simulated by program OBSERVE; (b) and (d) show orbits observed by Spacewatch from September 1990 to April 1995. The dashed line in (a) and (b) shows aphelion, Q 5 1.0 AU.

Comparing Imb to Iinj in Fig. 4, it is apparent that Iinj is dominated by the fragments of main-belt asteroids close to the n6 and 3 : 1 resonances and by the fragments of main-belt asteroids with high e. These are fringe asteroids, represented by dark squares in Figs. 4a and 4b, that require only small perturbations to evolve from stable, main-belt orbits to unstable Mars-crossing orbits. Farinella et al. (1993) independently predict a similar source distribution as a function of i and a. They perform a similar analysis, but use a more accurate theory to determine the boundaries of

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

the n6 and 3 : 1 resonances. They also assume a distribution of ejection velocities, ­N(n)/­n Y n23.25 for n . 100m/sec, and ­N(n)/­n 5 0 for n , 100 m/sec. Our agreement with their results supports the validity of our relatively simple injection model. In order to make a more quantitative test, we used INJECT to predict the yield of fragments into the 3 : 1 resonance from all the asteroids numbered less than 1000 that were also considered by Farinella et al. (1993). Given the same proper elements, and also the same ­N(n)/­n, we predict approximately the same yields (to within a factor of 2) for 33 of 36 cases. We made additional comparisons with the predictions of Morbidelli et al. (1994) to test our predicted yields of Mars-crossing fragments, injected by asteroids close to the n6 resonance. Given the same sources and the same mean elements, our predictions agree to within a factor of 2 in 17 of 20 cases. Assuming that the ensemble from which INJECT has selected source bodies (Imb) is a good representation of the orbits of smaller main-belt asteroids (5 m , d , 30 km), it follows that Iinj is a plausible representation of the orbits for the main-belt sources of the Earth approachers. Certain details in the distribution may be incorrect, such as the specific asteroids from Imb that INJECT has selected as parent bodies and the fact that the . 30 km asteroids are not perfect indicators of the actual sources, most of which are smaller. Depending upon the errors in the mean elements that we have used to define Imb and also upon the precision with which we have modeled the location of the resonances, those asteroids which appear to be closest to the instability boundary will vary. For example, by choosing to define Imb using osculation elements instead of proper elements, we obtain a different set of source bodies. This is not an appropriate procedure because it selects those bodies that are most near to Mars-crossing at the present epoch, but not necessarily on a time scale of 4.5 byr. Nonetheless, the resulting orbits concentrate in the nearly same regions of a, e, i space as the source bodies plotted in Figs. 4a and 4b. We therefore believe that the source population determined by INJECT is representative of the range of populations we can expect, given our uncertainties in the stability boundaries and mean elements for main-belt asteroids. Since our intention is not to identify actual parent bodies, but to determine a distribution of injected orbits that is realistic overall, our model is adequate for these purposes. From the above discussion, it follows that the most crucial uncertainties affecting Iinj are the uncertain dependence of Imb on d, and the uncertainty for n0 . As discussed in Section 4.1, changes in n0 affect the fraction, Finj , of main-belt bodies that contribute in Iinj , but do little to alter the shape of the resulting distribution. Hence, uncertainty in n0 mainly influences the normalization and shape for the predicted size distribution, which we discuss in a later

45

paper (Part II). Of greater significance to Iea is the unknown dependence of Imb on d. Hence, in Section 5.6 we qualify our conclusions in terms of these uncertainties. 5.2. Dynamical Evolutions In earlier papers (Wetherill 1979, 1985, Paper I), it was shown that the dynamical routes from the main belt to Earth-crossing orbits could be fast or slow. Asteroid ejecta injected into the 3 : 1 jovian resonance take the fast track, achieving Earth-crossing orbits in p106 yr. Bodies injected into Mars-grazing orbits at the inner edge of the main belt take the slow track. Occasional encounters between Mars and these inner-edge bodies randomly perturb their orbits until they are jostled closer to the n6 resonance. Proximity to the resonance increases the eccentricities of the perturbed orbits, thereby increasing their chances for Mars encounters. Over a time scale of p109 yr, both Mars encounters and n6 perturbations act together (‘‘synergistically’’) to force bodies at the inner edge to Earth crossing. As expected, our injected bodies follow similarly fast and slow tracks. As shown by Fig. 6b, the dynamical ages of bodies initially injected into the 3 : 1 resonance are p100 times shorter than the ages of fragments initially injected into Mars-crossing orbits outside the 3 : 1 resonance. The short-lived bodies achieve Earth-crossing orbits in p106 yr, and many (23%) go on to hit the sun. However, the steady-state distribution is dominated by those bodies that are removed from the 3 : 1 resonance by close encounters with Earth and Venus and are able to survive as Earth approachers for 107 to 108 yr. The slow-track bodies, on the other hand, generally reside in shallow, Mars-crossing orbits for more than 1 byr before they approach Earth. Smaller members of this population that are removed by collisions are replaced by collision fragments of the larger slow-track bodies. Once the orbits of the slow-track bodies approach Earth, however, encounters with Earth and Venus rapidly drive their dynamical evolution. As with the fast-track fragments, the steady-state population is composed of bodies surviving as Earth approachers for 107 to 108 yr. Despite the differences in dynamical ages for slow- and fast-track Earth approachers, Figs. 5a and 5b show that their orbital migrations are similar. Both distributions have concentrations of Amors (orbits with 1.017 , q , 1.3 AU) either in the range a 5 2.0–2.4 AU (slow-track) or a 5 2.4–2.6 AU (fast-track), in agreement with the predictions of Paper I and Bottke et al. (1996). The range in a of these Amors reveals the main-belt location of their parent bodies. The concentrations appear because the Amors have longer dynamical lifetimes than Apollos and Atens (bodies with q , 1.017 AU), and because Amors are efficiently resupplied by main-belt fragments. Apollos and Atens of both slow- and fast-track origin, on the other

46

DAVID L. RABINOWITZ

hand, are dispersed over a broad range in a and q. Owing to the scrambling effect of close encounters with Earth and Venus, their orbital distribution is not specific to their point of origin. Although Figs. 5a and 5b are in overall agreement with the previous predictions of Wetherill and Bottke et al., there are some minor, but important differences. Unlike Paper I, our present model does not predict a strong concentration of Earth approachers with a p 2.05 AU, extending from q 5 0.2 to 1.3 AU. This was an artifact of the earlier n6 calculations, corrected in a later paper (Wetherill 1991). Our present calculation also predicts a new feature. Fragments following both the slow and the fast tracks yield concentrations of low-e Amors with a 5 1.4–1.6 AU and q . 1.1 AU, separated from the high-e Amors by a gap in the range a 5 1.5–1.8 AU. This is a weak feature in Figs. 5a and 5b, but it persists when we take fragmentation and observational selection into account. These low-e Amors have no special place of origin and suffer multiple encounters with Earth and Venus prior to their arrival. As their route is circuitous, their supply rate is low. Of the 2000 evolutions computed by EVOLVE, only 4 produce orbits that fall within the low-e concentration. The feature persists in the steady state only because it is long lived. 5.3. Fragmentation Comparing the dynamical distributions in Figs. 5a and 5b with the orbital distributions we predict for d p 1.0 km and d p 0.01 km (Figs. 5c and 5d, respectively), it is apparent that fragmentation is an important process. Neither Fig. 5c nor Fig. 5d is equivalent to Idyn (the superposition of Figs. 5a and 5b), as we would expect if fragmentation were negligible. Instead, both distributions are dominated by the fast-track orbits in Fig. 5b. This dominance is most conspicuous for the smaller bodies, where the concentration of Amors with a 5 2.4–2.6 AU is a prominent feature (85% of the small body population derives from fast-track parents). For d p 1.0 km, a smaller fraction are of fasttrack origin (64%), but this is still twice the fraction comprising Idyn (33%). To understand how fragmentation increases the portion of fast-track Earth approachers in the steady state, consider the separate size distributions we should expect for the slow- and fast-track bodies. The effective main-belt ages for the fast-track bodies have a median value (kDt9lm p 3.2 3 107 yr) which is much less than the average collisional lifetime (tmb(d) p 4 3 108 yr) for kilometer-sized asteroids in the main belt. Hence, most kilometer-sized Earth approachers of fast-track origin are lost by hyperbolic ejection or by solar and planetary collisions before they have a chance to catastrophically fragment. Their steady-state

distribution of sizes remains proportional to the production distribution in the main belt, ­nprod(d)/­t p d23.0. Slow-track bodies, on the other hand, live long enough (kDt9lm p 1.6 3 109 yr) to suffer several disruptive collisions before ejection or collision. As discussed by Wetherill (1985) and by Greenberg and Nolan (1989), ongoing collisions allow the size distribution to relax from d23.0 toward the equilibrium distribution, d22.5. Since the fast-track Earth approachers have a steeper size distribution, they dominate the steady-state ensemble. A similar argument holds for Earth approachers with d 5 10 m. The small, slow-track bodies have a relatively flat, equilibrium distribution proportional to d22.5 because tmb(d) ! kDt9lm . The small, fast-track bodies have a steeper distribution, nearly proportional d23.0, because tmb(d) and kDt9lm are comparable. In Part II of this series, we present numerical results to support these arguments. Since the effects of fragmentation are small for fasttrack Earth approachers with d . 1.0 km (see Section 5.2), we can use the results of the calculation performed by EVOLVE to predict their steady-state number, Nfast . We start with the total number of fast-track, Earth-approaching orbits plotted in Fig. 5b, Nb 5 1807. This is the steady-state number we would expect in the absence of fragmentation if the cumulative production rate in the main belt were ­nprod /­t 5 [KFinj]21 3 802/106 yr21 (two evolutions for each of the 401 members of Iinj injected into the 3 : 1 resonance, each recorded once every 106 yr; the factor [KFinj]21 accounts for production in the main belt that does not contribute fragments to the fast-track ensemble). Given Eq. (1) and also the discussion in Section 2.2.3, a realistic value for ­nprod /­t at d 5 1.0 km is (2.5/3.0) 3 nmb(d 5 1 km)/tmb(d 5 1 km) p 0.83 3 5.3 3 106/3.8 3 108 yr 5 1.2 3 1022 yr21. Here we have extrapolated nmb(d 5 1 km) from the observed number, nmb(d) 5 1069 at d 5 30 km, assuming nmb(d) Y d22.5. Scaling Nb to account for this corrected production rate yields Nfast 5 21 Nb 3 1.2 3 1022 yr21/[KF inj 3 802/106 yr]. Hence, for our nominal calculation with n0 5 100 m/sec, Finj 5 4.3%, and K 5 48%, we obtain Nfast p 558. Finally, dividing this result by the predicted fraction of km-sized Earth approachers of fast-track origin yields total number, Nea 5 558/0.638 p 875 for d . 1.0 km. The above prediction is in marginal agreement with the number of Earth-crossing asteroids with d . 1.0 km (1500 to within a factor p2) estimated by Rabinowitz et al. (1994). However, since we expect that the number of Earth approachers is p40% greater than the number of Earth crossers (given current discovery statistics) our predicted value for Nea is low. If the assumptions of our nominal calculation were rigid, this would imply that main-belt asteroids are not a sufficient source. As we discuss in Section 2.2.1, however, our uncertainty for n0 covers the range 50 to 150

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

m/sec. By repeating the entire Monte Carlo calculation with n0 5 150 m/sec instead of n0 5 100 m/sec, we obtain Nea 5 1773, consistent with Rabinowitz et al. This change has negligible influence on the resulting orbit distributions (see Section 5.6). Alternatively, it is conceivable that a lower value for n0 is appropriate, but our nominal estimate for nmb(d) is low. As we discuss in Section 2.4, and summarize in Table I, the uncertainty for nmb(d) is a factor of 1.7 at d 5 10 km. Hence, our uncertainties for both n0 and nmb(d) accommodate our difference with Rabinowitz et al. 5.4. Comparison with Spacewatch Observations of Kilometer-sized Earth Approachers Accounting for the observational bias of the Spacewatch search, the predicted distribution in q vs a for kilometersized Earth approachers (Fig. 7a) is a good match to the distribution observed by Spacewatch (Fig. 7b). Both show distinct concentrations of Amors with a p 2.0–2.6 AU. As discussed above (Section 5.2), this is an intrinsic feature expected for Earth approachers originating from the main belt. We also anticipate the appearance in the observations of a secondary concentration of Amors with a p 1.3–1.5 AU and q p 1.0–1.15 AU, nearly coincident with a similarly prominent feature in Ibias . Finally, Fig. 9a shows that the relative numbers of observed Amors in these two concentrations are comparable to our expectation. The most significant discrepancy is for Amors with a 5 2.45–2.55 AU, where only p1/3 the expected number are observed. The distribution of Apollos and Atens in Fig. 7a also agrees with Fig. 7b. Both distributions show roughly the same scatter and cover roughly the same domain in q vs a. The fraction of observed Earth approachers that are Apollos and Atens (57 6 9% sampling error) closely matches our prediction (54%). Both the observed and the predicted distributions also show an apparent decrease in the concentration of Apollos as a increases. This trend is intrinsic to Iea (Fig. 5c). However, because the trend is more conspicuous in Ibias than Iea, selection effects are partially responsible. Comparing our prediction for sin(i) vs a (Fig. 7c) to the observations (Fig. 7d), there is again agreement, but it is only partial. Neglecting the observed orbits with sin(i) . 0.4, both the observed and the predicted distributions scatter within a triangular region narrowing from sin(i) p 0.0–0.4 at a 5 1.0 AU, to sin(i) p 0.0–0.1 at a 5 2.8 AU. In our model, this is the outcome of repeated encounters with Earth and Venus. Each two-body encounter preserves the approaching asteroid’s pseudo-energy in the rotating reference frame of the planet (Jacobi’s constant 5 2/a 1 Ïq(1 1 e) cos(i) when referred to Earth). For bodies that encounter Earth at their perihelion (q p 1.0 AU), any decrease in a requires a corresponding decrease in e, or

47

else an increase in sin(i). Encounter’s with q near Venus have a similar effect. Hence, the Earth approachers with a p 1.0 AU, having started with a . 2.0 AU, tend to have large values of sin(i). For the observed Earth approachers with sin(i) , p0.4, the above scenario may be valid. Although we predict values for sin(i) that are low by a factor p2 for bodies with a . 2.0 AU, this discrepancy does not challenge the theory. The parents of the observed Earth approachers could have higher i than we have assumed. Our model clearly fails, however, to explain the appearance of so many bodies with sin(i) . 0.4 (20% of the observations), extending to values as large as sin(i) p 0.8. This disagreement might indicate that we have neglected an important source for the Earth approachers (such as high-i comets, or perhaps Hungaria asteroids). Secular variations in i (which we ignore in our model) also might have an important influence (see Section 5.6). Despite our failure to predict the orbital distribution of kilometer-sized Earth approachers with large sin(i), our successful prediction of q vs a and of sin(i) vs a for low-i bodies offers strong support for our model, and for the hypothesis that main-belt fragments could be a dominant, if not a sufficient source. We cannot rule out a substantial contribution from short-period comets, as modeled by Wetherill (1991), since these objects achieve orbits similar to Earth approachers. However, if the main belt is the dominant source, we can now trust our model to predict the resulting orbits for the SEAs. 5.5. Comparison with Spacewatch Observations of SEAs Accounting for the observational bias of the SEAs, the values we predict for q vs a (Ibias in Fig. 8a) are not a good match to the observations (Fig. 8b). Compared to our prediction, too many of the observations are clustered near q p 1.0 AU and a 5 1.0 to 1.4 AU. This discrepancy appears more clearly in Fig. 9b. Of the observed SEAs with q . 0.9 AU, 5 have a , 1.2 AU and none have a . 2.2 AU, while the predicted numbers are less than 1 and p3, respectively. Note that OBSERVE predicts an observational bias toward the detection of SEAs with either q p 1.0 AU or Q p 1.0 AU. The concentration of such orbits are enhanced in Ibias relative to Iea (Fig. 5d). However, these selection effects do not significantly enhance the detection of low-e SEAs with both q p 1.0 AU and a p 1.0 AU. Furthermore, repeating the Monte Carlo calculation with n0 5 50 and 150 m/sec, or with an enhanced contribution from sources at the inner edge of the main belt, does not remove the discrepancy with observation (see Section 5.6). Hence, the observed concentration of low-e SEAs is both an unpredicted and an intrinsic feature of their orbital distribution. The observed deficit of Atens

48

DAVID L. RABINOWITZ

with q , 0.5 AU among the SEAs is similarly unpredicted and intrinsic. While our model has trouble predicting q vs a for the SEAs, it has much better success predicting their distribution in sin(i) vs a. Figures 8c and 8d show that both the predicted and the observed orbits fill a triangular region, as we discuss for kilometer-sized Earth approachers in Section 5.4. As for the kilometer-sized Earth approachers, the observed SEAs with a .2.0 AU also have values of sin(i) about twice as high as predicted. Again, however, this may imply that main-belt sources for these bodies have higher i than we have assumed. Unlike the kilometer-sized Earth approachers, however, there are no unanticipated bodies with sin(i) . 0.4. The most conspicuous discrepancy is a peculiar cluster of points appearing in the predicted distribution (near a 5 0.7 AU and sin(i) 5 0.025), but not appearing in the observations. This cluster is an artifact of program OBSERVE. The program has picked out one orbit with sin(i) p 0 and Q p 1.0 AU, which is a very favorable orbit for detection. The resulting cluster of points represents repeated detections of this one object. Given the success with which our Monte Carlo calculation predicts the orbits of kilometer-sized Earth approachers, our failure to predict q vs a for the SEAs leads us to conclude that a main-belt origin is awkward for these small bodies. Even when the five SEAs with q . 0.9 AU and a , 1.2 AU are neglected (consisting of 1991 TT, 1991 VG, 1992 DU, 1992 JD, and 1992 YD3), the remaining distribution is weighted toward values of a that are lower than we expect, and a deficit of Atens persists. These residual discrepancies might then be explained by sampling error in the sparse observations, uncertainties in the model, or perhaps by the influence of dynamical and physical causes that we have not accounted for. For example, perhaps the relative contribution of fast-track Earth approachers is smaller than we have assumed. By correcting our model we would then predict fewer Earth approachers with a p 2.5 AU, in better agreement with the observed numbers of such bodies for both large and small d. However, because the fast track also yields a greater concentration of low-e Earth approachers than the slow-track, decreasing the fast-track contribution leads to fewer predicted SEAs with low e, in poorer agreement with observation (see Section 5.6). There are other possible explanations for our failure to predict the SEA orbits. Solar radiation might lower the albedo of the Atens with low q, so that they are harder to detect than the SEAs with q p 1.0 AU. However, there are no physical observations supporting this hypothesis. Perhaps our theory of secular resonances is not accurate enough to predict the SEA orbits. But then we must find a reason that the theory works for kilometer-sized Earth approachers with sin(i) , p0.4. We are left with the conclusion that our model is reasonable, but fails to predict the

FIG. 10. Perihelion, q, versus semimajor axis, a, predicted for d 5 1.0–1.26 km (top) and d 5 10.0–12.6 m (bottom). From left to right: n0 5 50 km/sec, n0 5 100 m/sec, n0 5 150 m/sec, and Nmb(a, e, i) defined by the known asteroids with d . 10 km.

orbits of the SEAs because a significant fraction are not of main-belt origin. As discussed in Rabinowitz (1994), the p5 SEAs with the low-e orbits that cause the most problem for the theory of main-belt origin comprise p20% of the observed population. However, they represent a smaller fraction of the intrinsic population when observational bias is corrected. Rabinowitz assumes that the total population consists of a near-Earth belt of SEAs with q 5 0.9–1.1 AU and Q , 1.4 AU, and that this belt is superimposed on a background population with the same orbital distribution as kilometersized Earth approachers. This near-Earth belt then comprises 5 6 3% of the total population of SEAs. Given the above discussion, however, it is clear that this is a lower limit. A less conservative estimate might include all of the SEAs with q . 0.9 AU, since their distribution of orbits is skewed to low values of a compared to the model. This orbital range covers 60 6 20% of the total population modeled by Rabinowitz. Hence, p1/2 of the SEAs may not originate from the main belt. 5.6. Uncertainties in the Orbital Distributions In the foregoing discussion, we point out that the most significant source of uncertainty affecting the shape of our predicted orbital distributions is Nmb(a, e, i). These predictions are relatively insensitive to uncertainty in the remaining parameters that control the calculation (see Table I). To show this, we have altered n0 and Nmb(a, e, i) within their ranges of uncertainty and then used INJECT and EVOLVE to recompute the steady-state distribution of orbits for Earth approachers with d p 1.0 km and with d p 0.01 km. Figures 10a–10d and Figs. 10e–10h show the

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

results for large and small bodies, respectively. Note that Figs. 10b and 10f are the same as Figs. 5c and 5d, except that we use fewer points to represent the distributions. For the other figures, we preserve our nominal assumptions, except that n0 5 50 m/sec for Figs. 10a and 10e, n0 5 150 m/sec for Figs. 10c and 10g, and Nmb(a, e, i) is determined by the known asteroids with d . p10 km (instead of d . p30 km) for Figs. 10d and 10h. Using a smaller size cutoff, here, to define Nmb(a, e, i) artificially biases Iinj toward sources at the inner edge of the main belt because the smaller asteroids have not been completely sampled and are preferentially detected nearer to Earth. The distributions plotted in Figs. 10a–10c and 10e–10g verify that changing n0 has little effect. Figures 10a and 10c each show the same two concentrations of kilometersized Amors that appear in Fig. 10b (from a 5 1.4–1.6 AU and from a 5 2.0–2.6 AU). There is some variation in the population density, but it is not significant. Figs. 10a–10c also reveal similar distributions of kilometer-sized Apollo and Atens. There is slightly more disparity between Figs. 10e–10g. For example, a clump of bodies with low e appears only in Fig. 10e (from a 5 1.2 to 1.6 AU and from q 5 1.0 to 1.2 AU). However, this and other variations are not enough to fix the disparity with the orbits of the observed SEAs. Figures 10e–10g all predict too many Atens with q , 0.5 AU and also too many Apollos and Amors with a . 2.0 AU. Comparing Figs. 10d and 10h to Figs. 10b and 10f, it is apparent that our alteration to Nmb(a, e, i) has a strong effect. For both d p 1.0 km and d p 0.01 km, the densest concentrations of bodies now appear among the Amors with a 5 1.4–1.7 AU and q 5 1.1–1.3 AU. Also, low-e bodies with q p 1.0 AU are underabundant compared to our nominal prediction. Using OBSERVE to bias these predictions, neither one makes a good match to the observations. The lesson is that we must be careful to properly define Nmb(a, e, i) and that the our nominal definition is a relatively good one. Another possible source of error is program OBSERVE, which we use to model the detection biases. However, because OBSERVE has been shown to be a reliable program, we do not believe it is significant source of error. As described in Rabinowitz (1993, 1994), OBSERVE accurately predicts the distribution of discovery variables (angular rate, distance at discovery, apparent magnitude, etc.) of the Spacewatch observations. More recently, Bottke et al. (1996) confirm that selection biases predicted by OBSERVE match the biases they independently determine by calculating the ratio of Earth-impact probability to Earth-encounter velocity. As discussed by Shoemaker et al. (1990), this ratio is an alternative measure of observational bias. Thus, it is unlikely that errors in program OBSERVE would account for the differences between our predictions and the observed orbits for SEAs.

49

As a final source of error, we must acknowledge that there are additional resonances that our model does not take into account. Morbidelli et al. (1995) show that mean motion resonances with Jupiter exterior to the 3 : 1 resonance (i.e., 5 : 2, 7 : 3, and 2 : 1) can deliver main-belt fragments to the inner Solar System. However, the paucity of observed Amors with a . 2.5 AU suggests that these resonances play a minor role. Our model ignores the n16 secular resonance, which pumps up the inclinations of asteroids with a p 2.0 AU (Scholl et al. 1989). We also assume that the influence of secular resonances are small for bodies with a , 2.0 AU. Finally, there are Kozai resonances that can protect bodies from encounters with the inner planets (Michel and Thomas 1996). Such resonances are important in special cases, such as lunar and Martian ejecta (Gladman et al. 1995, Gladman and Burns 1996). The question remains, however, whether or not any of these resonances are important influences for Earth approachers of mainbelt origin. 6. CONCLUSIONS

We have used a Monte Carlo program to predict the orbital distribution of the Earth approachers ranging in diameter from 10 m to 10 km assuming they originate as main-belt fragments. Considering the probable oversimplifications of the model, it is able to account rather well for the diverse effects of close planetary encounters, resonant perturbations, collisional fragmentation, and observational selection. A comparison of our dynamical evolutions with numerically computed predictions shows that we can account for the effects of the n6 and 3 : 1 resonances in a way that mimics their computed effects. Taking these resonances into account, our model selects parent bodies for the Earth approachers that have a distribution of orbits similar to the parent bodies determined by an independent, numerical analysis. We also predict ejection yields from such bodies that are comparable to independently predicted values. For these reasons, we believe our model is reasonably realistic. Assuming a range of injection velocities for main-belt fragments that span the values determined from observations of asteroid families, and assuming collisional lifetimes consistent with self-similar fragmentation in the main belt, our model successfully predicts the orbital distribution for the kilometer-sized Earth approachers. Both the observations and our predictions show a major concentration of Amors with a in the range 2.0 to 2.6 AU and a secondary concentration of Amors with a 5 1.3 to 1.5 AU. The relative number of observed Earth approachers that are Atens and Apollos also agrees with our prediction. Finally, given our uncertainties in the total number of main belt supplying the Earth approachers, and also in the ejection velocities of collisional fragments, our calculation yields a total abun-

50

DAVID L. RABINOWITZ

dance of Earth approachers larger than p1.0 km that is consistent with estimates derived from the observed number. It is thus likely that main-belt asteroids are a dominant source for kilometer-sized Earth approachers. We require no other sources to explain the observed orbits, although a significant contribution from short period comets is possible. Extending our model to predict the orbits for SEAs, our predictions are relatively unsuccessful. Our model fails to predict a distinct concentration of observed SEAs with q p 1.0 AU and a , 1.4 AU, and we predict an abundance of Atens among the SEAs that is large compared to the observed number. Varying our critical assumptions within their respective ranges of uncertainty does not remove these discrepancies. Thus, a significant fraction of the SEAs, no less than 5%, but perhaps as large as 50%, may not originate from the main belt. Further calculations will be required to conclusively test the hypothesis that impact ejecta from Mars, low-e Amors, or some other local source, is the supplier of these exceptional bodies. APPENDIX A This appendix reviews the critical algorithms used by the Monte Carlo calculation (described in Section 2) to account for the effects of (a) close planetary encounters with planets; (b) distant encounters with Jupiter outside its sphere of influence; (c) secular variation of the planet orbits; (d) the n6 secular resonance; and (e) the 3 : 1 mean motion resonance with Jupiter. Most of these algorithms have been described in detail ¨ pik 1963, 1976). elsewhere (Paper I, Wetherill 1979, 1985, Arnold 1965, O Our implementation, here, is only slightly changed from Paper I, mainly to account for the effects of the n6 and 3 : 1 resonances in a more realistic way.

Having calculated P at a given time step, EVOLVE then randomly chooses whether or not a planetary encounter occurs, weighting the decision by P. If an encounter occurs, EVOLVE chooses the planet, the encounter geometry, and r0 , as described by Arnold (1965). If a longrange encounter with Jupiter is chosen, then EVOLVE chooses the geometry and r0 as described by Wetherill (1979). A collision occurs if r0 , rc , and the evolution ends. Otherwise, EVOLVE determines the body’s ¨ pik’s formulas. The time step is incremented, postencounter orbit using O and the procedure repeats until (a) the body’s age exceeds the Solar System’s age, Dt . 4.5 3 109 yr, (b) the body collides with a planet, (c) the body is ejected on a hyperbolic orbit, (d) the body’s orbital period exceeds 1000 yr, or (e) q falls below a minimum value, qmin 5 0.05 AU (i.e., solar collision).

A.2. Secular Variations of the Orbits of the Planets To model the secular variations of e and i for the major planets, i 5 1 to 8, EVOLVE uses periodic functions, Ei(t) and Ii(t), chosen to approximate the variations computed by Cohen et al. (1973) using the secular theory of Brouwer and van Woerkom (1950). We also use the maximum and minimum values for e computed for each planet by Quinn et al. (1991) to fix the corresponding extremes for Ei(t). These functions determine the osculating values of e and i for every planet at each time step during the course of a small body’s evolution (see Section A.1). Both Ei(t) and Ii(t) are of the form Fa(t) or Fb(t), where Fa(t) 5 F0 1 F1 sin(g1t 1 u1) sin(g2t 1 u2), and Fb(t) 5 F0 1 F1 sin(g1t 1 u1) 1 F2 sin(g2t 1 u2), and F0 , F1 , F2 , g1 , g2 , u1 , and u2 are constants (see Table AI).

A.1. Close Encounters Assuming that close encounters between a planet and a small body are effectively two-body interactions, and also that the times between the encounters are long compared to the precession rates of the two ¨ pik determined an expression for the mean time between enbodies, O counters as a function of the impact parameter, r0 , and the orbital elements of the small body and the planet. He also derived formulas to compute the change in the small body’s orbit at each encounter in terms of the geometry of approach. Given these formulas, Arnold, and later Wetherill, developed computer programs to model the orbital evolution of small bodies owing to repeated, random approaches within each plan¨ pik. et’s sphere of action, Ra , as defined by O ¨ pik’s equations in a manner similar to ArProgram EVOLVE uses O nold and Wetherill, but modifies the procedure to account for events occurring between the times of planetary encounters (e.g., fragmentation and long-range perturbations). Hence, the program determines the fate of a given body at regular time steps spaced by dt 5 103 yr, which is much less than the average time between encounters. At each time step, EVOLVE determines the likelihood, P, of an encounter with any of the eight major planets (Mercury to Neptune) within time, dt, and within distance r0 5 10 times the gravitational capture radius, rc , for each planet. The program uses the same algorithms used by Arnold (1965) and Wetherill (1979), including a procedure developed by Arnold to approximate the effects of distant encounters (10rc , r0 , Ra) and a procedure developed by Wetherill to account for long-range encounters interior to Jupiter (within 0.65 AU of Jupiter’s perihelion).

A.3. Effects of the n6 Resonance Wetherill (1979) describes a method for accounting for the effects of the n6 secular resonance that is only slightly changed in the present application. For a given body with osculating elements a, e, i, longitude of perihelion, g, and with proper elements a0 , e0 , i0 , and g0 , the effect of the resonance is to produce periodic variations in e of amplitude, e1 and mean value, e0 . Secular perturbations to variables a and i are ignored. The time dependence for e is determined by the time-dependent phase, g 2 n6 , where n6 is the phase of the secular resonance (nearly equal to the longitude of perihelion for Saturn). Variables e1 and e0 are constant in the absence of close encounters with the planets. If both e0 and i0 are small, then e1 , e, e0 , g, g0 , and n6 are related by e 20 5 e 2 1 e 21 2 2ee1 cos(g 2 n6)

(A.1)

e 2 5 e 20 1 e 21 1 2e0e1 cos(g0 2 n6).

(A.2)

and

Note that Eq. (A.2) corrects a typographical error appearing in Wetherill (1979). This earlier paper also used the following simple formula, based

51

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

TABLE AI Constants Defining Secular Perturbations for the Planetsa Planet

Function

F0

F1

F2

Mercury

E1(t) I1(t) E2(t) I2(t) E3(t) I3(t) E4(t) I4(t) E5(t) I5(t) E6(t) I6(t) E7(t) I7(t) E8(t) I8(t)

0.190 7.159 0.0345 1.691 0.0290 1.462 0.0600 2.927 0.0435 0.365 0.0480 0.860 0.0375 1.000 0.0092 0.680

0.0515 2.378 0.0345 1.691 0.0290 1.462 0.0270 2.427 0.0155 0.185 0.0350 0.240 0.0315 0.180 0.0066 0.125

— 0.25 — — — — 0.031 0.50 0.002 — 0.002 — 0.0 — 0.0 0.075

Venus Earth Mars Jupiter Saturn Uranus Neptune

a

2f/g1 310 6 yr

2f/g2 310 6 yr

u1 (rad)

u2 (rad)

Form

3.5 0.96 4.0 1.9 4.0 1.9 0.093 1.2 0.055 3.8 0.055 3.8 0.82 3.8 0.55 0.58

0.85 0.099 0.1 0.1 0.1 0.1 2.0 0.16 0.82 0.049 0.82 0.049 0.0 0.44 0.0 1.9

0.000 4.712 0.000 4.712 0.000 4.712 0.000 0.000 0.000 0.000 3.142 0.000 3.142 1.571 0.000 0.000

0.000 0.000 3.142 0.000 0.000 3.142 0.000 0.000 0.000 3.142 0.000 0.000 0.000 0.000 0.000 0.000

a b a a a a b b b a b a b a b b

For function Ii(t), the units for F0 , F1 , and F2 are degrees. Forms a and b correspond to functions Fa and Fb , respectively, as defined in Section A.2.

on semi-analytic perturbation theory, to determine e1 as a function of a0 : e1 5 6(0.023)[ua0 2 ares u]20.598,

(A.3)

where ares 5 2.08 AU. Furthermore, e1 . 0 if a , ares, e1 , 0 if a . ares, and the maximum value of ue1u is limited to emax 5 0.999. Thus, for orbits with a0 p ares , the absolute value of the n6 amplitude rises sharply to ue1u 5 emax . The present application differs from Wetherill (1979) in the way that EVOLVE determines g 2 n6 and e1 , and in the treatment of the n6 resonance for asteroids with high i and a . 2.5 AU. Instead of choosing random values for g 2 n6 at random encounter times, EVOLVE increments g 2 n6 at every time step, and at a rate, ­(g 2 n6)/­t, that depends on a0 . Based on the results of Yoshikawa (1987), who provides an analytic estimate for ­(g 2 n6)/­t vs a0 for i0 , p208, EVOLVE assumes ­(g 2 n6)/­t 5 218.080 yr

21

for a0 , 1.2 AU,

­(g 2 n6)/­t 5 220.550 yr AU21 3 (2.08 AU 2 a0) for 1.2 , a0 , 2.08 AU, 21

and ­(g 2 n6)/­t 5 35.200 yr21 AU21 3 (a0 2 2.08 AU) for 1.2 , a0 , 2.08 AU. Here we linearly extrapolate ­(g 2 n6)/­t for 1.2 AU , a , 2.0 AU from Yoshikawa’s results for a p 2.0 AU. We also assume that ­(g 2 n6)/­t is constant for a , 1.2 AU. For a , p1.5 AU, however, these assumptions are unimportant because ue1u is small (,0.03). Unlike Wetherill (1979), EVOLVE also determines ares as a function of e0 and i0 . Interior to the 3 : 1 mean motion resonance with Jupiter (a0 , 2.48 AU), the program assumes 8/3 ares 5 2.08 AU 1 14.5[sin(i*) 0 2 0.275(e* 0 2 0.1)] ,

(A.4)

where i*0 5 i0 if sin(i0) $ 0.02, and e*0 5 e0 if 0.5 $ e0 $ 0.1. Otherwise,

i*0 5 0 if sin(i0) , 0.02, e*0 5 0.1 if eo , 0.1, and e*0 5 0.5 if e0 . 0.5. This expression is based on the analytic theory of Knezevic et al. (1991). Using it to determine ares , EVOLVE then uses Eq. (A.3) to determine e1 . Exterior to the 3 : 1 resonance (a0 . 2.52 AU), however, EVOLVE assumes that the location of the n6 resonance is independent of a0 , and centered at ires , where sin(ires) 5 0.317 1 0.608(e*0 2 0.1).

(A.5)

Furthermore, Eq. (A.3) is superseded for a0 . 2.52 AU by e1 5 6[1.64 3 1029/(sin[i0] 2 sin[ires])4 1 9.0 3 1024]1/2,

(A.6)

where e1 , 0 for sin(i0) , sin(ires) and e1 . 0 for sin(i0) $ sin(ires). As in Eq. (A.3), EVOLVE restricts the range of the n6 amplitude to ue1u , emax , with emax 5 0.999. Given the above assumptions, EVOLVE then uses the following procedure to account for the effects of the n6 resonance. At t 5 0, EVOLVE initializes a0 5 a and i0 5 i, where a and i are the initial osculating elements for a given body. As in Wetherill (1979), we ignore secular perturbations to elements a and i. Their osculating and proper values are always equal, and remain constant until a close encounter with a planet. At t 5 0, EVOLVE also sets e0 5 e, but only if Eqs. (A.3)–(A.6) yield a value for e1 that also satisfies Eq. (A.1) for some value of g 2 n6 . If so, Eq. (A.1) initializes g 2 n6 . Otherwise, EVOLVE chooses g 2 n6 5 0 or f radians, and solves Eqs. (A.1) and (A.3)–(A.6) simultaneously. The resulting solution initializes e1 and e0 in terms of a0 , i0 , e, and g 2 n6 . If g 2 n6 librates (i.e., if e1 . e0 , forcing g 2 n6 to oscillate within a restricted range) then EVOLVE chooses a random, initial sign for ­(g 2 n6)/­t. At each future time step, EVOLVE then increments g 2 n6 by d (g 2 n6) 5 dt 3 [­(g 2 n6)/­t] and uses Eq. (A.1) to determine the corresponding value for e in terms of e0 and e1 . If g 2 n6 lies outside the libration range, then the sign of ­(g 2 n6)/­t is reversed, and the phase increment is reflected to bring g 2 n6 back into its restricted range.

52

DAVID L. RABINOWITZ

The above procedures continues for each time step until a planetary encounter occurs. If the body survives the encounter, and the post-encounter value for a is outside the range for the 3 : 1 resonance (a , 2.48 AU or a . 2.52 AU), then a0 , i0 , e0 , e1 , and ­(g 2 n6)/­t are re-initialized, keeping g 2 n6 at the pre-encounter value, plus a shift determined by the planetary perturbation. However, if Eqs. (A.1)–(A.6) do not yield a consistent solution for g 2 n6 , then EVOLVE sets g 2 n6 5 0 or f and reinitializes e1 and e0 as at t 5 0. If, on the other hand, the post-encounter value for a is inside the range for the 3 : 1 resonance, then EVOLVE ignores the effects of the n6 resonance until another encounter brings a outside the 3 : 1 resonance. At that point, EVOLVE again reinitializes the relevant n6 variables as at t 5 0.

A.4. Effects of the 3 : 1 Mean Motion Resonance with Jupiter To account for the effects of the 3 : 1 mean motion resonance, EVOLVE employs the following procedure. For a given body at t 5 0, if a is in the range 2.48 to 2.52 AU, EVOLVE assigns an initial value for a base eccentricity, ebase 5 e. At each future time step, and until a planetary encounter moves the body outside the resonance range, EVOLVE sets e 5 ebase 1 Desin , where Desin is a sinusoidal component with amplitude, 0.15, and period, 5 3 10 4 yr. At each time step, EVOLVE also alters ebase by adding two random components: a diffusive increment, Dediff , and a jumping increment, Dejump . Component Dediff causes ebase to random walk by 60.15 with a diffusion time of 4.0 3 10 6 yr. The jumping component, on the other hand, causes ebase to make jumps uniformly distributed in amplitude from 20.6 to 0.6, at random times, t, exponentially distributed with a 1/e time of 3.3 3 10 6 yr. The long-term variation of e is thus the combination of sinusoidal, diffusive, and jumping components. These perturbations are reinitialized, as at t 5 0, each time planetary perturbations move a into the range for the 3 : 1 resonance.

A.5. Criteria for Planet Crossing When program INJECT initially determines the orbits of injected fragments of main-belt asteroids, Iinj , the program admits only those orbits that are potential planet crossers. Such an orbit must meet at least one of the following criteria: (i) a within the range for the 3 : 1 resonance (2.48 , a , 2.52 AU); (ii) a minimum perihelion less than the maximum aphelion of Mars; and (iii) a maximum aphelion no less than 0.65 AU smaller than the minimum perihelion of Jupiter. To evaluate criteria (ii) and (iii), INJECT initializes a0 , i0 , e0 , e1 , and g 2 n6 for the given body, as described in Section A.3, above. The resulting upper limit for e, as determined by Eq. (A.1), determines the minimum perihelion and maximum aphelion. The maximum aphelion for Mars and the minimum perihelion for Jupiter, on the other hand, are determined by the secular variations of their orbits, as described in Section A.2. An additional criterion that INJECT uses to limit admission into Iinj is that the orbit of a given fragment’s parent asteroid must not already satisfy criteria (i), (ii), or (iii), above. This criterion ensures that the sources for the Earth approachers are represented by stable bodies that yield unstable ejecta only for perturbed orbits. The subsequent calculation can also assume that the present distribution of source orbits has been stable for most of the age of the Solar System. Sections 4.1 and 5.1 discuss the implications of this selection procedure.

APPENDIX B This appendix describes the recursive algorithm that program FRAGMENT uses to model the collisional cascade of a given fragment into a set of descendant fragments, each with a known diameter, d. Only fragments with d . dmin are counted, so the resulting set will be a empty if none of the predicted fragments meet this criterion. Section 2.2.3 already

describes the parameters that determine the time, t0 , of the first disruptive collision and the size distribution of the first generation of fragments. These same parameters constrain the iterative disruption of each subsequent generation of fragments, as follows. (1) For each logarithmic interval in diameter, log(d) to log(d) 1 0.1, starting at log(dmin), and ending with log(d2) 2 0.1, where d2 is the size of the largest fragment (see Section 2.2.3), FRAGMENT determines the number of fragments, Dn0(d), that will survive to the end of the effective main-belt lifetime, Dt9, of the parent orbit. This is given by Dn0(d) 5 N0 Dnfrag(d0 , d) f0 , where Dnfrag(d0 , d) 5 nfrag(d0 , d)2 nfrag(d0 , 10 0.1 d), f0 5 e2[(Dt92t0)/tmb(d)], and N0 is a normalizing constant equal to unity for the first generation of descendants. The number of first-generation fragments that will not survive is Dn0(d) 5 N0 Dnfrag(d0 , d)(1 2 f0). To account for non-integer values of either Dn0(d) or Dn0(d), FRAGMENT randomly chooses between the two nearest integer values, weighting the two choices inversely by their difference from Dn0(d) or Dn0(d). (2) For those fragments in each log(d) interval that do not survive, FRAGMENT chooses a single, effective main-belt time, t1 , for their mutual disruption. Although a more realistic calculation would choose an independent time of disruption for each non-surviving body, this would require a large increase in computing time. FRAGMENT therefore sets t1 equal to the time when half the non-surviving, first-generation bodies would normally be destroyed by independent, random disruptions in the main belt: t1 5 t0 2 tmb(d) log[(1 1 fo)/2]. (3) If Dn0(d) . 0 for a given log(d) interval, except the lowest one including dmin , then FRAGMENT repeats steps (1), (2), and (3) with d0 5 d, t0 5 t1 , and N0 5 Dn0(d). This is a recursive operation because FRAGMENT will keep repeating steps (1) and (2) each time it returns to step (3). On each iteration, the program thereby handles a subsequent generation of fragments. The iterations stop when there are no more fragments to disrupt (N0 5 0), or all fragments are too small to break into descendant fragments larger than dmin . Once FRAGMENT has implemented the above algorithm, it will have determined Dn0(d) for every generation of descendants. Summing over all generations yields the total number of surviving fragments as a function of size for a given source fragment.

ACKNOWLEDGMENTS I thank George Wetherill (Carnegie Institution of Washington) for his expert and ingenuous help, without which this paper could not have been written. I also thank the members of the Spacewatch team (University of Arizona) for their ongoing discoveries of Earth approachers, large and small, P. Farinella (University of Pisa) for constructive comments, M. Holman (University of Toronto) for numerical investigations, Bill Bottke (University of Arizona) for the use of his collision probability maps, and B. Gladman (Cornell University) for extensive comments. This work was supported by NASA.

SOURCES FOR EARTH-APPROACHING ASTEROIDS?

REFERENCES ARNOLD, J. R. 1965. The origin of meteorites as small bodies. II. The model. Astrophys. J. 141, 1536–1547. BELTON, M. J. S., J. VEVERKA, P. THOMAS, P. HELFENSTEIN, D. SIMONELLI, C. CHAPMAN, M. E. DAVIES, R. GREELEY, R. GREENBERG, J. HEAD, S. MURCHIE, K. KLAASEN, T. V. JOHNSON, A. MCEWEN, D. MORRISON, G. NEUKUM, F. FANALE, C. ANGER, M. CARR, AND C. PICHER 1992. Galileo encounter with 951 Gaspra: First pictures of an asteroid. Science 257, 1647–1652. BINZEL, R. P., AND S. XU 1993. Chips off of asteroid (4) Vesta: Evidence for the parent body of the basaltic achondrite meteorites. Science 260, 186–191. BINZEL, R. P., S. XU, S. BUS, AND E. BOWELL 1992. Origins for the nearEarth asteroids. Science 257, 779–782. BOTTKE, W. F., M. C. NOLAN, R. GREENBERG, AND R. A. KOLVORD 1994. Collisional lifetimes and impact statistics for near-Earth asteroids. In Hazards Due to Comets and Asteroids (T. Gehrels, Ed.), pp. 337–358. Univ. of Arizona Press, Tucson. BOTTKE, W. F., M. C. NOLAN, H. J. MELOSH, A. M. VICKERY, AND R. GREENBERG 1996. Origin of the Spacewatch small Earth-approaching asteroids. Icarus 122, 406–427. BROUWER, D., AND A. VAN WOERKOM 1950. Astron. Papers Am. Ephem. Naut. Alman. 13, part II. CELLINO, A., V. ZAPPALA, AND P. FARINELLA 1991. The size distribution of main-belt asteroids from IRAS data. Mon. Not. R. Astron. Soc. 253, 561–574. COHEN, C. J., E. C. HUBBARD, AND C. OESTERWINTER 1973. Planetary elements for 10,000,000 years. Celest. Mech. 7, 438–448. DAVIS, D. R., S. J. WEIDENSCHILLING, P. FARINELLA, P. PAOLICCHI, AND R. P. BINZEL 1989. Asteroidal collisional history: Effects on sizes and spins. In Asteroids II (R. Binzel, T. Gehrels, and M. Matthews, Eds.), pp. 805–826. Univ. of Arizona Press, Tucson. DOHNANYI, J. S. 1969. Collisional model of asteroids and their debris. J. Geophys. Res. 74, 2531–2554. FARINELLA, P., CH. FROESCHLE´ , CL. FROESCHLE´ , R. GONCZI, G. HAHN, A. MORBIDELLI, AND G. B. VALSECCHI 1994a. Asteroids falling into the Sun. Nature 371, 315–317. FARINELLA, P., CL.FROESCHLE´ , AND R. GONCZI 1994b. Meteorite delivery and transport. In Asteroids, Comets, and Meteors 1993 (A. Milani, M. Di Martino, and A. Cellino, Eds.), pp. 205–222. Kluwer Academic, Dordrecht. FARINELLA, P., R. GONCZI, CH. FROESCHLE´ , AND C. FROESCHLE´ 1993. The injection of asteroid fragments into resonances. Icarus 101, 174– 187. FUJIWARA, A., P. CERRONI, D. R. DAVIS, E. RYAN, M. DI MARTINO, K. HOLSAPPLE, AND K. HOUSEN 1989. Experiments and scaling laws on catastrophic collisions. In Asteroids II (R. Binzel, T. Gehrels, and M. Matthews, Eds.), pp. 240–268. Univ. of Arizona Press, Tucson. GEHRELS, T. 1991. Scanning with charge-coupled devices. Space Sci. Rev. 58, 347–375. GLADMAN, B., AND J. BURNS 1996. The delivery of martian and lunar meteorites to Earth. Bull. Am. Astron. Soc. 28, 1054. GLADMAN, B., J. BURNS, M. DUNCAN, AND H. LEVISON 1995. The dynamical evolution of lunar impact ejecta. Icarus 118, 302–321. GREENBERG, R., AND M. NOLAN 1989. Delivery of asteroids and meteorites to the inner Solar System. In Asteroids II (R. Binzel, T. Gehrels, and M. Matthews, Eds.), pp. 778–804. Univ. of Arizona Press, Tucson. HOLSAPPLE, K. A. 1994. Catastrophic disruption and cratering of Solar System bodies: A review and new results. Planet. Space Sci. 42, 1067– 1078.

53

KNEZEVIC, Z., A. MILANI, P. FARINELLA, CH. FROESCHLE´ , AND CL. FROESCHLE´ 1991. Secular resonances from 2 to 50 AU. Icarus 93, 316–330. LEMAITRE, A., AND A. MORBIDELLI 1994. Calculation of proper elements for high inclination asteroidal orbits. Celest. Mech. Dynam. Astron. 60, 29–56. MARTELLI, G., E. V. RYAN, A. M. NAKAMURA, AND I. GIBLIN 1994. Catastrophic disruption experiments: Recent results. Planet. Space Sci. 42, 1013–1026. MCFADDEN, L. A., D. J. THOLEN, and G. J. VEEDER 1989. Physical properties of Aten, Apollo, and Amor asteroids. In Asteroids II (R. Binzel, T. Gehrels, and M. Matthews, Eds.), pp. 442–486. Univ. of Arizona Press, Tucson. MICHEL, P., AND F. THOMAS 1996. The Kozai resonance for near-Earth asteroids with semimajor axes smaller than 2 AU. Astron. Astrophys. 307, 310. MILANI, A., AND Z. KNEZEVIC 1994. Asteroid proper elements and the dynamical structure of the asteroid belt. Icarus 107, 219–254. MORBIDELLI, A., AND M. MOONS 1995. Numerical evidences on the chaotic nature of the 3 : 1 resonance. Icarus 155, 60–65. MORBIDELLI, A., R. GONZCI, CH. FROESCHLE, AND P. FARINELLA 1994. Delivery of meteorites through the nu(sub6) secular resonance. Astron. Astrophys. 282, 955–979. MORBIDELLI, A., V. ZAPPALA, M. MOONS, A. CELLINO, AND R. GONCZI 1995. Asteroid families close to mean motion resonances: Dynamical effects and physical implications. Icarus 188, 132–154. ¨ PIK, E. J. 1963. The stray bodies in the Solar System. I. Survival of O cometary nuclei and the asteroids. Adv. Astron. Astrophys. 2, 219– 262. ¨ PIK, E. J. 1976. Interplanetary Encounters. Elsevier, Amsterdam. O QUINN, T. R., S. TREMAINE, AND M. DUNCAN 1991. A three million year integration of the Earth’s orbit. Astron. J. 101, 2287–2365. RABINOWITZ, D. L. 1991. Detection of Earth-approaching asteroids in near real time. Astron. J. 101, 1518–1529. RABINOWITZ, D. L. 1993. The size distribution of the Earth-approaching asteroids. Astrophys. J. 407, 412–427. RABINOWITZ, D. L. 1994. The size and shape of the near-Earth asteroid belt. Icarus 111, 364–377. RABINOWITZ, D. L. 1996. Observations constraining the origins of Earthapproaching asteroids. In Completing the Inventory of the Solar System (T. Rettig and J. Hahn, Eds.), pp. 13–28. Astron. Soc. of the Pacific Conf. Ser. 107. RABINOWITZ, D., E. BOWELL, E. M. SHOEMAKER, AND K. MUINONEN 1994. The population of Earth-crossing asteroids. In Hazards Due to Comets and Asteroids (T. Gehrels and M. Matthews, Eds.), pp. 285–312. Univ. of Arizona Press, Tucson. RABINOWITZ, D. L., T. GEHRELS, J. V. SCOTTI, R. S. MCMILLAN, M. L. PERRY, W. WISNIEWSKI, S. M. LARSON, E. S. HOWELL, AND B. E. A. MUELLER 1993. Evidence for a near-Earth asteroid belt. Nature 363, 704–706. SCHOLL, H., CH. FROESCHLE, H. KINOSHITA, M. YOSHIKAWA, AND J. G. WILLIAMS 1989. Secular resonances. In Asteroids II (R. Binzel, T. Gehrels, and M. Matthews, Eds.), pp. 845–861. Univ. of Arizona Press, Tucson. SCOTTI, J. V. 1994. Computer aided near Earth object detection. In Asteroids, Comets, and Meteors 1993 (A. Milani, M. Di Martino, and A. Cellino, Eds.), pp. 17–30. Kluwer Academic, Dordrecht. SHOEMAKER, E. M., R. F. WOLFE, AND C. S. SHOEMAKER 1990. Asteroid and comet flux in the neighborhood of Earth. In Global Catastrophes

54

DAVID L. RABINOWITZ

in Earth History (V. L. Sharpton and P. D. Ward), pp. 155–170. Geol. Soc. of America, Boulder. VAN HOUTEN, C. J., I. VAN HOUTEN-GROENVELD, P. HERGET, AND T. GEHRELS 1970. The Palomar-Leiden survey of faint minor planets. Astron. Astrophys. Suppl. 2, 339–448. WETHERILL, G. W. 1967. Collisions in the asteroids belt. J. Geophys. Res. 72, 2429–2444. WETHERILL, G. W. 1979. Steady state populations of Apollo-Amor objects. Icarus 37, 96–112. WETHERILL, G. W. 1985. Asteroidal sources of ordinary chondrites. Meteoritics 20, 1–23. WETHERILL, G. W. 1988. Where do the Apollo objects come from? Icarus 76, 1–18. WETHERILL, G. W. 1991. End products of cometary evolution: Cometary origin of Earth-crossing bodies of asteroidal appearance. In Comets in the Post-Halley Era Vol. 1 (R. L. Newburn, Jr., et al., Eds.), pp. 537–556. Kluwer Academic, Dordrecht.

WETHERILL, G. W., AND J. G. WILLIAMS 1979. Origin of differentiated meteorites. In Origin and Distribution of the Elements (L. H. Ahrens, Ed.), pp. 19–31. Pergamon, Oxford. WILLIAMS, J. G. 1973. Meteorites from the asteroid belt? Eos 54, 233. [Abstract] WILLIAMS, D. R., AND G. W. WETHERILL 1994. Size distributions of collisionally evolved asteroidal populations: Analytical solution for selfsimilar collision cascades. Icarus 107, 117–128. WISDOM, J. 1983. Chaotic behavior and the origin of the 3/1 Kirkwood gap. Icarus 56, 51–74. WISDOM, J. 1985. Meteorites may follow a chaotic route to Earth. Nature 315, 731–733. YOSHIKAWA, M. 1987. A simple analytical model for the secular resonance n6 in the asteroids belt. Celest. Mech. 40, 233–272. ZAPPALA, V., P. FARINELLA, Z. KNEZEVIC, AND P. PAOLICCHI 1984. Collisional origin of the asteroid families: Mass and velocity distributions. Icarus 59, 261–285.