Statistics of close approaches between asteroids and planets: Project SPACEGUARD

ICARUS

88,

292-33s

Statistics

(1990)

of Close Approaches between Asteroids Project SPACEGUARD A. MILANI,“.?

M. CARPINO,S

AND

and Planets:

F. MARZARIS

*lkpar~tnen~ of‘Astronomy. Cornrll Univcrsiry. Itircrc~o. Nc14, Yorh 14853; 1.Dipar~imrnto di Mtrfc~tntrtic~tr, Utliversitd di Piss, Via Buonurrori 2. I-56100 Pi.w. Ilcrliu; $ Ossc~rvcrlorio Astronomicw di Miltrno-Brerrr. Vicr Brcwr 28, I-20121 Milww. Iiulicr; trnd ii Dipurtimento di Fisiccc. Univ~~rsitci di Pudovcr, Podovrr, Iicditr Received

November

6, 1989; revised

April

30, 1990

Within Project SPACEGUARD, numerical integration for -200,000 years of a large number of orbits of planet-crossing asteroids has been used to generate a data base of close approaches to the major planets. These data are used to test the possibility of predicting-for a specific asteroid and for an entire population-the number of close approaches within a given distance, by using statistical theories such as the ones by Kessler (1981, Icarus 48,39-48) and Wetherill(1%7, J. Geophys. Res. 72,2429-2444). The fundamental assumption of the statistical theories is the absence of mean motion resonance locking between the asteroid and the target planet. The orbits belonging to the class Toro-as defined in Milani et al. (1989b, Zcurus 78,212-269)-violate this assumption and indeed the numbers of close approaches in the data base are very different from the statistical predictions. For the other dynamical classes the statistical theories give very good predictions whenever the assumptions on the relative orbital geometry are fulfilled; both Kessler theory and Wetherill theory have difficulties in predicting accurately close approaches between orbits which are either almost coplanar or almost tangent. We propose a modified version of Kessler theory which can cope with these singular cases and give very stable and reliable predictions of the number of close approaches at large to moderate distances. Even the use of the instantaneous orbital elements allows a fairly accurate computation of the total number of close approaches for an entire population; therefore, an accurate computation of the orbits is not necessary. However, these results cannot be extrapolated in an accurate way to predict the frequency of very close approaches, including physical collisions. By means of correlation analysis we have also tested the assumptionused in Monte Carlo studies-that the orbital elements of a planet-crossing orbit change only because of the close approaches; this is essentially true for the semi-major axis of such an orbit, provided it is not in resonance, either with the target planet (Tore class) or with Jupiter (Alinda class). The qualitative dynamical classification proposed in Milani et al. (1989b) is confirmed as a useful tool, although some specific attributions to a class have to be corrected on the basis of objective tests of the definitions of the classes. o MO Academic PM. IW.

been one of the greatest triumphs of human reason, the mechanisms involved in either And then a Plank in Reason broke, very close approaches or physical colliAnd I dropped down, and downsions of two celestial bodies have not been And hir u World, af every plunge, understood as well as it would be necesAnd Finished knowing--thensary. Until quite recently, not only to the Emily Dickinson, from Poem 280, 1861 layman but even to many scientists working While the understanding of the motion of in related fields, a planetary collision apthe celestial bodies as determined by their peared as a powerful symbol of the breaklong-range gravitational interactions has ing down of human reason and of the very 292 I. INTRODUCTION

0019-1035/90

$3.00

Copyright 0 1990 by Academic Press, Inc All rights of reproducfmn in any form reserved.

PROJECT SPACEGUARD

annihilation of conscience, as in the quotation of the New England poet we are using as a header. More than this psychological and philosophical resistance, which has anyway been overcome by the last generation of scientists, the main stumbling block has been the lack of a suitable analytical technique. The main tool of Celestial Mechanics, namely the expansion of both the problem and the solution in series and the use of perturbation theory, is inapplicable in a neighborhood of the singularity of collision. Although it is always possible to represent the orbit of an asteroid undergoing a close approach to a planet by means of some local approximation-e.g. the path is approximately a hyperbolic two-body orbit in the planetocentric frame-the dynamical behavior is globally chaotic, with Lyapounov exponents of the order of unity in 10,000 years. The orbit becomes unpredictable over time spans a few times longer than the inverse of the Lyapounov exponent, much shorter than the expected collisional lifetimes. The attempts to model the occurrence of both close approaches and collisions for planet-crossing orbits have been based upon a statistical approach, starting from the pioneering effort by ijpik (1951). A statistical theory can be constructed by adopting three simplifying assumptions: (1) The orbits between two consecutive close approaches are assumed to be regular and modeled in some simplified way, such that an explicit computation is possible. In the simplest possible approach, the orbits are modeled as keplerian ellipses with constant semimajor axis a, eccentricity e, and inclination I, uniformly precessing with time; more elaborate models exploit secular perturbation theories, obtained by either analytical or semianalytical methods, in which the constants of the motion are “proper” a, e, and Z as opposed to the corresponding instantaneous, or “osculating,” elements. (2) Close approaches and collisions are

293

modeled as purely random events. More specifically, whenever the keplerian ellipses representing the osculating orbits of one asteroid and one planet are close enough to allow for a close approach, the position of the planet and of the asteroid on the orbit (i.e., the anomalies) are assumed to be randomly distributed and uncorrelated; mean motion resonances, which constrain the two anomalies, are not allowed in the model. Moreover, each close approach is seen as a random event uncorrelated with any previous one. (3) The orbital elements a, e, I of the asteroid are assumed to change only as a result of close approaches. Since at deep close approaches these changes can be computed, e.g., by an hyperbolic two-body encounter approximation, the orbital evolution can be described as a random walk. For each present state, the probability (hence the expected waiting time) for a given change in the elements can be explicitly computed; this makes possible the use of Montecarlo simulations. In some more sophisticated theories, other events which can result in changes of the elements, such as the capture in some resonance, are taken into account (Wetherill 1988). These statistical models have played an essential role in allowing a quantitative, although approximate, model of phenomena such as collision and cratering rates on the terrestrial planets, transport of meteorites, asteroids, and extinct comets to Earthcrossing orbits, etc. However, are these models legitimate? That is, are the assumptions such as (l), (2), and (3) good enough approximations of the more complex reality to allow for reliable quantitative estimates? The main difficulty is that the line of reasoning underlying the statistical theories is quite long, and if any of the “planks” breaks down-e.g., if any one of (1), (2), and (3) turns out to be a very poor approximation-then unreliable conclusions would be obtained. The purpose of this paper is to test the assumptions underlying the statistical theo-

294

MILANI,

CARPINO,

ries of collisions, close approaches, and orbital evolution for planet-crossing asteroids. This can be done by using the data base-containing a complete list of close approaches and time series for the orbital elements-generated by numerical integration within Project SPACEGUARD (Milani et al. 1989a, b). Assumption (2) can be tested by using as input the time series of orbital elements for a large enough sample of planet-crossing orbits, by using a statistical theory to compute the expected frequency of close approaches (within a given minimum distance, or better impact parameter), and by comparing the expected number with the actual number of close approaches as recorded in the data base. Assumption (3) can be tested for a given “proper” element by computing the changes it undergoes as a result of each recorded close approach and then by forming the correlation coefficient between these changes due to close approaches and the overall changes. Assumption (I) is a little more difficult to test, because it can take many different forms, depending upon the type of approximate orbit computation used; we can however compute the expected number of close approaches by means of a statistical theory using the instantaneous orbital elements at different times, and see how this number changes with time. For the purpose of this paper, we have used the data on a set of 89 asteroids with either Earth-crossing or at least Earthapproaching orbits computed over a timespan of 200,000 years. This paper is organized as follows: in Section 2 (and Appendix I) we describe the SPACEGUARD data base and the way it has been postprocessed to extract the informations necessary for this work. In Section 3 we test the statistical theory for the frequency of collisions developed by Kessler (1981), which uses the precessing ellipses approximation with constant a, e, and I, and we discuss whether such a theory can be modified to provide also a reliable prediction of the number of close approaches;

AND

MARZARl

the formal derivation for the improved formulas are given in Appendix 2. In Section 4 we test the statistical theory for collision frequency developed by Wetherill (1967), which assumes that information is also available on the secular behavior of some angular variables, especially the argument of perihelion o. In Section 5 we discuss the correlation between changes in the asteroid semimajor axis and the changes resulting from the close approaches; we also explain why we cannot perform the same test on the eccentricity and the inclination. Our conclusions are presented in Section 6. Although our conclusions are complex, we can schematically summarize them as follows. The statistical theories for the frequencies of collisions and close approaches are essentially correct in a statistical sense. That is, every statistical theory can fail in a catastrophic way in its prediction for a single asteroid and/or for a specific short time span. However, if the averages over a large enough population over a long enough time is compared with the experimental results-obtained by numerical integration-the agreement is very good: the cases in which the simplifying assumptions used by the statistical theories fail occur less often than the cases in which they are good approximations; moreover the statistical methods can fail both ways, and a large scale average almost eliminates the discrepancies. We also show that the statistical predictions of these theories are based upon simplified models of the geometry of the two orbits of the asteroid and of the planet, and these geometrical constructions could be improved to remove singularities and uncomputable cases. We present here some of these improvements and they result in much more reliable predictions of close approach frequency. We show that the applicability of any statistical theory depends upon the understanding of the fundamental dynamical property of each class of orbits; we have used here the dynamical classes defined in Milani ef al. 1989b (from now on quoted as

PROJECT SPACEGUARD

Paper I) on the basis of the SPACEGUARD data, and this classification turns out to be an essential tool: e.g., the Geographos dynamical class is the one for which both Kessler and Wetherill theory work, the Toro class is defined by the failure of both theories, while for the Kozai class Kessler theory fails and Wetherill theory works. The relevance of the dynamical classification is further confirmed by the correlation analysis between changes in a and close approaches; e.g., the Geographos members have orbits for which the semimajor axis changes “only because of” the close approaches with Earth and/or Venus, the Alinda and Toro members are those with a changing “only because of” resonances, where of course “only because of” must be rigorously defined in terms of correlation coefficients. However, the classification of dynamical types we are using should be understood as a conceptual tool, that is a hypothesis to be continuously tested, not as a fact. As an example, we introduce in this paper three corrections to our previous classification; we think that this shows that our taxonomy (the class definitions) was meaningful and significant, but the specific attribution of each orbit (or part thereof) to a given class was and is subject to error; the data in the tables are not corrected for this change, which has been made possible by the analysis of the discrepancies occurring therein. 2. THE SPACEGUARD CLOSE APPROACH DATABASE

The data on orbital evolution and close approaches were generated by numerical integration of 410 orbits, starting from the initial conditions as given in an asteroid catalog (Marsden 1986) for the known planetcrossing objects. For some of these asteroids the orbits have later been improved by taking into account recent observations, but our computations use the initial conditions as known in 1986; for this reason, we still use the temporary designation for some recently numbered asteroids. These orbits

295

were computed for the time span between -100,000 and 100,000 years from present, within a model including gravitational perturbations by the seven major planets (Venus to Neptune), by using an experimental parallel and vector supercomputer developed by E. Clementi at the IBM Kingston, NY laboratory. The procedure used to perform the computation is described in Milani et al. (1989a). Given the available computing power, the limit to the size of the experiment was not set by the CPU time, but by the capability to store, handle, postprocess, and understand the output data. Although a significant amount of postprocessing was performed on line, with digital filtering techniques such as those described in Carpino et al. (1987), the amount of output was more than 2 Gigabytes, and to devise ways to analyze it efficiently has been, and still is, a challenge. A first selection of the output data was performed by singling out the orbits which, within the integration time span of 200,000 years, have a perihelion distance which is-at least at some time-not larger than the aphelion distance of the Earth orbit plus 0.1 AU. The number of these Earthapproaching orbits in our sample turns out to be 89. For each of these 89 three data sets have been prepared in the postprocessing phase to be used for the present work: (1) Post-filtered orbital elements. The time series for the orbital elements has been passed through a further digital filtering stage, designed in such a way that all the short periodic perturbations (with a period shorter than 86.2 years) have been removed, and all the perturbations with periods longer than 150 years are preserved without distortions. In this way it is possible to sample the orbital elements only once every 20,000 days = 54.75 years without loss of information on the long term behavior. Thus for each orbit we have 3661 regularly spaced data points. (2) List of minimum distance states at each close approach. The latter are defined as every occurrence of a distance smaller

296

MILANI,

CARPINO.

than 0.1 AU between the asteroid and either Venus, Earth, or Mars, and of a distance smaller than 1 AU from Jupiter. By interpolation, the time at which the minimum distance occurs is found and the position and velocity vectors of both the planet and the asteroid are computed. This allows one to derive not only the minimum distance, but also the relative velocity, hence the impact parameter (allowing correction for the gravitational focusing effects; see Sections 3 and 4). (3) List of orbital changes produced by each close approach. Whenever the asteroid is first recorded to be inside the close approach sphere-of radius 0.1 AU for the inner planets, 1 AU for Jupiter-the orbital elements are computed, and compared with those at the last recorded step inside the same sphere. This also requires interpolation, because the velocities were not stored. Although the procedures used to obtain (1). (2), and (3) are straightforward in almost all the cases, some problems occurred in special cases. Having examined 32 1,560 cases of close approach, it is not surprising to find that some rare combinations of circumstances do occur. Some of these cannot be properly handled by a completely automated postprocessing software, and some turn out to be caused by numerical exceptions and instabilities which cannot be handled properly in a parallel computation. On the other hand, for the purpose of this work we are not really interested in understanding the dynamics of each close approach, but rather in the statistical behavior of large numbers of such events; thus we need only to apply a quality control algorithm to be able to discard the wild cases. For the minimum distance state, this resulted in discarding a very small number of very special events (such as temporary captures as a satellite of Jupiter). For the orbital changes, some orbits have too large a number of “defective” close approaches and have to be removed from the correlation analysis of

AND

MARZARI

Section 5. Details on the computational procedures used in the postprocessing stage, on the numerical instabilities and special event problems, and on the quality control are given in the Appendix 1. In the analysis of the close approach data presented in the following sections we shall often refer to orbits belonging to the dynamical classes Geographos, Toro, Kozai, Alinda, Eros, Oljato, and “Comet.” The definitions are given in Paper I, Section 5 and Table IV. As it is apparent from the column headed “Close approaches” in that table, the distinction between the different classes is especially relevant for the comparison of the actual number of close approaches with the numbers predicted by any statistical theory. As an example, the Toro type orbits are protected from close approaches to the Earth by mean motion resonances; the Kozai types avoid any node crossing with the Earth, and therefore any deep close approach, by means of a protection mechanism resulting from secular perturbations; the Geographos types are the real Earth-crossers and Earthapproachers, together with the very high eccentricity and low inclination Oljato types. The Eros types have a perihelion above 1 AU; the Alinda types raise and lower their perihelion as a result of resonant perturbations by Jupiter. One of the purposes of this paper, apart from checking the validity of the statistical theories for the entire population of near-Earth asteroids, is to confirm-or, if necessary, to correctthe classification performed in an intuitive way (mostly through a visual inspection of the plots) in our previous paper, by means of objective measures of the protection mechanisms. In this paper, the results will be presented as class averages as well as global averages over all the orbits; however, we have also generated tables with the results of the statistical theories and correlation analysis for each orbit, and used these as objective tests for our classification.

PROJECT

SPACEGUARD

One final comment on the SPACEGUARD data base. As we have already pointed out, all the orbits we have examined are chaotic; thus the fact that the initial conditions have been copied from a catalog of real, observed asteroids does not imply that the orbits we have computed are orbits of real objects. The exponential divergence of nearby orbits is so effective in mixing the phase space that after a few tens of thousands years two initial conditions different by a fraction of an arcsecond can evolve into completely unrelated orbits; the times and depths of the close approaches of the two orbits are totally uncorrelated. Thus the use of an asteroid name and/or number has to be understood as a mnemonic code rather than as a reference to the real physical object. For example, the statement “1937 UB (Hermes) has 8015 close approaches to Venus within 0.1 AU, including 4 within 0.001 AU” would be nonsense if it was to refer to the real asteroid. These statements always have to be understood as referring to the object in our sample, whose initial conditions are compatible (within the observational uncertainty) with those of the real objects with that name, and even for those specific initial conditions only statistical statements should be used (e.g., “a moderate inclination, shallow Venus crossing orbit such as that of 1937 UB results in many thousands of close approaches”). We apologize in advance with the readers for sometimes skipping these cautionary formulas to write shorter sentences, but we would like to stress that everything we say has to be understood in the sense of statistical mechanics. 3. COLLISION

PROBABILITIES METHOD

BY

KESSLER

The probabilities of collisions between two orbiting objects were first derived by Gpik (1951; 1976). His formula was later derived in a more concise way and improved upon by Kessler (Kessler and CourPalais 1978, Kessler 1981), Shoemaker et

297

al. (1979), and Steel and Baggaley (1985). Kessler defines an average spatial density for an orbit with fixed a, e, Z and the angular arguments R, o, 4 (longitude of the node, argument of perihelion and mean anomaly) changing linearly with time; that is, the orbit is described as a keplerian ellipse of fixed shape precessing uniformly, and the average is the average over the angle variables, which coincides with a long term average provided there are no resonances (by the standard averaging theorem; Arnold 1976, Section 51). The resulting density S(r, ~3)at a radius r from the Sun and at an ecliptic latitude p-the density is axially symmetric; thus, the longitude does not appear-is given by (Kessler 1981, formula (21)):

S(r, PI = 1

=

2n3ra d(sin*

Z - sin*p)(r - q)(Q - r)’

(1) provided r is larger than q = a(1 - e) and smaller than Q = a(1 + e), and (PI is smaller than I; S is zero if any of these constraints is not satisfied. Then the most likely number of close approaches per unit time through a cross section A can be computed from S and from the relative velocity V of the two orbits: Or, PI = Sk, P)VA.

(2)

Then formula (2) needs to be averaged over the target orbit; in the approximation of a circular orbit for the target planet, if ZM is the mutual inclination-the inclination of the asteroid orbit with respect to the orbital plane of the planet-and r = a’ is the semimajor axis of the planet: P=

VA 2~r~ua’ sin IM d((~’ - q)
where V can be computed from the known value of the velocity of the asteroid orbit when node crossing occurs, that is, when

298

MILANI.

CARPINO,

r = a’ and /3 = 0: v2

Y

dw(l

YE

-

I

+pa(l

e2>

cos

g2 Y d-1

z

M-

-e’) r?

. (sm’ fM - I) +p

2 1 --a, r

i

)

(4)

with p = GMsuN and p”I = G(MsuN + M’), with M’ the planet mass (actually there are four possible values for the relative velocity vector, all with the same length). The cross section A for a close approach within an impact parameter smaller than D is just rD2, and the relationship between impact parameter and minimum distance &IN is defined by the two-body hyperbolic encounter formula: D2 = (1 + s)

d&.

(3

More complicated formulas have been derived to take into account the eccentricity e’ of the orbit of the planet (Kessler 1981); as we shall see, the main problem with (3) is not in the approximation e’ = 0; thus, we have used this approximation in all the following discussion. To test the accuracy and reliability of Kessler’s formula (3), we use the time series for the orbital elements of the asteroid and of the target for a(t), e(t), Zdt), a’(t), then we compute V(t) from (4) and P(t) for a given D from (3); the integral of P(t) over a time span T gives the expected number Nk of close approaches within the value D of the impact parameter over the same time span. Since the actual number of close approaches can be deduced from the minimum distance status file (by applying the correction (5)), a direct comparison is possible. This comparison has been performed for all the orbits in our sample; however, two orbits became hyperbolic for part of the integration time span; thus, the comparison actually refers to 87.78 orbits, each computed for 201,300 years.

AND

MARZARl

Before discussing the results, we would like to point out the a priori obvious limitations of such a comparison. Kessler’s formula contains D only in the cross section A, and the spatial density S is computed only at the distance a’. That is, the formula is actually obtained as A times a probability per unit area; the latter is obtained as a limit for D + 0, e.g., if a’ is just below q, the expected number of close approaches is zero for every D, even when a’ + D > q and therefore close approaches are indeed possible. Moreover, the formula is singular for a’ = q, a’ = Q and sin ZM = 0, that is it predicts an infinite number of close approaches for every D. This singularity occurs because the average spatial density is actually infinite, thus the limit for D -+ 0 of the close approach probability per unit area is infinite even though the probability for a finite cross section area is always finite. We shall see later how this could be improven upon; as it is, because of this limit process, Kessler’s formula is expected to better predict the occurrence of close approaches at very low D-such as the actual collisions-than the shallow close approaches. On the other hand, the very deep close approaches are rare, and the SPACEGUARD close approach data base may not contain enough of these for a statistically reliable comparison. For example, in the comparison Table I for all the orbits, the column containing the Kessler’s formula expected numbers for D I 0.1 AU has an entry of 152,340 for the Earth, while at D 5 0.002 the expected number is 61. The corresponding numbers of close approaches found from the numerical integrations are 126,199 and 65, respectively. Since the numbers for D % 0.1 are very large, the 18.8% difference is significant; for D C= 0.002 it is not. To show this, we have computed the probability that such a discrepancy could occur by chance, by assuming a Poisson distribution and using the probability function

PROJECT SPACEGUARD TABLE COLLISION

Planet

D

NK

VENUS VENUS VENUS VENUS VENUS

VENUS VENUS

0.100

0.050 0.020 0.010

0.005 0.002 0.001

(4)

124551.2 31550.5 5092.5 1300.7 1278.6 325.2 321.8 52.0 51.6

32516.8 5202.7

12.9

EARTH

0.100

EARTH

0.050 0.020

38084.9 6093.6

34562.8 5758.9

0.010

1523.4

1470.9

EARTH EARTH EARTH

EARTH EARTH

o.cKl5 0.002 0.001

152339.6

380.9 60.9 15.2

MARS

0.100

69210.7

MARS MARS MARS MARS MARS MARS

0.050 0.020

17302.7 2768.4

0.010

692.1

0.005 0.002

173.0 27.7

0.001

6.9

JUPITER

1.000

JUPITER JUPITER JUPITER JUPITER

0.500 0.200

776.9 124.3

0.100

31.1

JUPITER JUPITER

0.050 0.020 0.010

CK

3107.5

7.8 1.2 0.3

129240.7

370.2 60.4 15.3 69478.7 17541.4 2778.7 691.9 172.5

27.5 6.9 8752.0 1065.9 130.2 31.4

7.8 1.2 0.3

nK

(5)

(6)

0.099 o.ooo

130067.3

13.0

(KESSLER)-ALL

NOb.

(3)

(2)

I

PROBABILITIES

N~rn

29424 4666 1076 284 41 13

299

0.100 0.109 0.189 0.135

0.000 0.000 0.000 0.011

0.237 0.068

126199 34094

0.001

0.572

0.188 0.111

0.000 0.000

5828

0.045 0.017 380 0.002 65 -0.065 15 0.015

1498

ORBITS CKVPB

P,(Km)

Pc(obs)

(8)

(9)

(10)

(11)

0.056 0.070 0.087 0.172 0.125 0.230 -0.008

0.000 0.000 0.000 0.000 0.018 0.075 0.526

1.93 1.93 1.93 1.93 1.93 1.93 1.93

1.76 1.81 1.84 1.86

1.70 1.73 1.75 1.58

1.89 1.90 1.91

1.69 1.50 2.31

3.55 3.55 3.55 3.55 3.55 3.55 3.55

2.75 2.98 3.18 3.30 3.35 3.49 3.57

2.80 3.09 3.43 3.61 3.65 3.96 3.92

0.25 0.25 0.25 0.25 0.25 0.25 0.25

0.25 0.25 0.25 0.25 0.25 0.25 0.25

0.23 0.23 0.23 0.22 0.26 0.28 0.20

3.77 3.77 3.77 3.77 3.77 3.77 3.77

10.40 4.76 3.93 4.01 4.08 4.06 4.00

14.76 5.42 3.91 2.92 0.99 0.03 0.00

0.024 0.000 0.014 0.006 -0.012 0.179

65947 0.048 16208 0.065 2598 0.064 614 0.120 178 -0.028 30 -0.080 6 0.143

0.335 0.355 0.462

11623 -1.156 1135 -0.375 131 -0.052 27 0.140 5 0.434 0.217 0 2.000

0.000 -0.282 o.oocl -0.063 0.286 -0.006 0.266 0.151 0.213 0.439 0.647 0.214 0.733 2.000

0.001 0.001

P,(K)

(7)

0.000 0.263 -0.018 0.496 -0.026 0.318 -0.074 0.544 0.018 0.000 0.000

nKm

0.235

0.295 0.291 0.540

0.052 0.000 0.079 0.000 0.067 0.000 0.119

0.001

-0.032 0.319 -0.086 0.344 0.145 0.459 0.000 0.017 0.485 0.248 0.209 0.648 0.741

Note. Comparison of the expected number of close approaches with each planet up to a maximum impact parameter D (column I), computed with Kessler method both in its original formulation (Ns, column 2) and with the modification proposed in the text (N s,,,, column 3) (see Appendix 2) with the corresponding number recorded in the SPACEGUARD integration (Nub,, column 4). Columns 5 and 7 give the relative discrepancy with respect to the observed value, for instance,

” =

(NK

f

NC,,,,)/2

Columns 6 and 8 give the probability of a discrepancy equal to or greater than the observed one, assuming for Nobp a Poisson NK) (see Eq. (6): IIs = &,F(N, NK) (where the distribution F(N,,,,, Columns 9 to I I summation is extended over 0 5 N 5 Nobr if Nob% 5 Ns, over Nabs -< N 5 p otherwise). give the expected number P, of physical collisions per IO9 years for the three methods, obtained by resealing the corresponding values N to the collision cross section rr@: P, = (AT/lo9 year) (D,lD)‘N (where AT = 201,216.B years is the integration timespan); the values of P, are not simply proportional to N because the collision impact parameter D, is also a function of the encounter velocity (see Eq. (7)). These tables have not been corrected for the changes in the dynamical classification introduced in this paper; in particular the header of the second half of Table IV is totally wrong. Total number of orbits = 87.78.

MILANI,

300

CAKPINO,

N

= exp C (log v - log j) - v (6) cj=l 1 for the probability when the expected integral of F gives discrepancy could

that the outcome is N number is V; a suitable the probability that the be as bad, or worse, by

TABLE COLLISION

D

Planet

Class

NK

(1)

(2)

= EROS

No. of orbits

PROBABILITIES N~rn

NOb.

(3)

(4)

AND

MARZARI

chance. From this computation as reported in Table I we can see that the discrepancy between Kessler’s formula and the actual number of close approaches to the Earth becomes statistically not significant for D C= 0.01. The Tables I, 11, and III contain the following information: for a given value of D (between 0.1 and 0.001 AU) and for each

11

(KESSLER)-PROTECTED CK (5)

nK

(6)

BKWZ (7)

CLASSES nKm

E(K)

P,(K”J)

(8)

(9)

(10)

J’c(obs)

(11)

= 20.94

EARTH EARTH EARTH

0.100 0.050 0.020

644.1 161.0 25.8

4506.2 520.1 28.6

4728 235 0

-1.520 -0.374 2.000

0.000 0.000 0.000

-0.048 0.755 2.000

0.001 0.000 0.000

0.09 0.09 0.09

0.87 0.37 0.11

0.86 0.15 0.00

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

20518.4 5129.6 820.7 205.2 51.3

21149.0 5260.2 829.4 205.9 51.4

20024 4498 744 180 52

0.024 0.131 0.098 0.131 -0.014

0.000 0.000 0.004 0.040 0.479

0.055 0.156 0.109 0.134 -0.012

0.000 0.000 0.001 0.036 0.483

0.32 0.32 0.32 0.32 0.32

0.33 0.33 0.32 0.32 0.32

0.31 0.28 0.29 0.28 0.31

Class

= KOZAI

No. of orbits

= 9.17

EA.RTH EARTH EARTH

0.100 0.050 0.020

3 166.5 791.6 126.7

2904.3 748.9 122.9

696 14 0

1.279 1.930 2.000

0.000 0.000 0.000

1.227 1.927 2.000

0.000 0.000 0.000

0.46 0.46 0.46

0.43 0.44 0.45

0.11 0.01 0.00

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

4337.7 1084.4 173.5 43.4 10.8

4248.2 1077.7 172.9 43.0 10.7

5249 1318 214 40 7

-0.190 -0.194 -0.209 0.081 0.431

0.000 0.000 0.001 0.339 0.153

-0.211 -0.201 -0.213 0.073 0.422

0.000 0.000 0.001 0.357 0.161

0.13 0.13 0.13 0.13 0.13

0.13 0.13 0.13 0.13 0.13

0.16 0.16 0.16 0.12 0.08

Class

= TORO

No. of orbits

= 2.29

VENUS VENUS VENUS VENUS VENUS

0.100 0.050 0.020 0.010 0.005

7569.9 1892.5 302.8 75.7 18.9

4978.0 1158.6 216.8 62.4 17.6

4465 697 47 10 2

0.516 0.923 1.463 1.533 1.618

0.000 0.000 0.000 0.000 0.000

0.109 0.497 1.287 1.447 1.591

0.000 o.oou 0.000 0.000 o.oco

6.57 6.57 6.57 6.57 6.57

3.59 3.40 4.25 5.09 5.91

3.26 2.28 1.01 0.71 0.69

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

8411.8 2102.9 336.5 84.1 21.0

7999.3 2216.4 339.2 84.3 21.2

4474 681 28 2 0

0.611 1.022 1.693 1.907 2.000

0.000 0.000 0.000 0.000 0.000

0.565 1.060 1.695 1.907 2.000

0.000 0.000 0.000 0.000 O.OCCl

6.52 6.52 6.52 6.52 6.52

6.20 6.99 6.59 6.54 6.56

3.50 2.43 0.65 0.19 0.00

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

1698.2 424.5 67.9 17.0 4.2

1971.4 447.3 68.2 17.0 4.2

1754 454 72 15 4

-0.032 -0.067 -0.058 0.124 0.060

0.086 0.074 0.326 0.373 0.581

0.117 -0.015 -0.054 0.123 0.058

0.000 0.364 0.340 0.375 0.582

0.27 0.27 0.27 0.27 0.27

0.33 0.29 0.28 0.27 0.27

0.26 0.27 0.27 0.23 0.28

See Note

to Table

I.

Note.

301

PROJECT SPACEGUARD TABLE COLLISION Planet

D

PROBABILITIES

NK

(2) Class

= GEOGRAPHOS

(KESSLER)-NONPROTECTED

NKm

NOb.

(3)

(4)

No. of orbits

III

CK

(5)

nK

(‘3)

CKW8

(7)

CLASSES nKm

PC(K)

Pc(Km)

ebbs)

(8)

(9)

00)

(11)

= 33.38

VENUS VENUS VENUS VENUS VENUS

0.100 0.050 0.020 0.010 0.005

88429.2 22107.3 3537.2 884.3 221.1

93741.8 22445.6 3519.3 874.2 218.3

90095 21915 3528 810 207

-0.019 0.009 0.003 0.088 0.066

0.000 0.098 0.443 0.006 0.181

0.040 0.024 -0.002 0.076 0.053

0.000 0.000 0.438 0.015 0.234

3.60 3.60 3.60 3.60 3.60

3.58 3.54 3.53 3.52 3.52

3.51 3.54 3.66 3.31 3.43

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

102912.7 25728.2 4116.5 1029.1 257.3

91069.0 24187.1 4061.9 1033.0 258.0

96366 27543 4909 1284 321

0.066 -0.068 -0.176 -0.220 -0.220

0.000 0.000 0.000 0.000 0.000

-0.057 -0.130 -0.189 -0.217 -0.217

0.000 0.000 0.000 0.000 0.000

6.34 6.34 6.34 6.34 6.34

5.11 5.74 6.23 6.37 6.37

5.68 6.91 8.08 8.61 8.65

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

25306.1 6326.5 1012.2 253.1 63.3

26264.2 6522.5 1024.1 254.4 63.2

24649 6281 982 239 80

0.026 0.007 0.030 0.057 -0.234

0.000 0.286 0.175 0.198 0.024

0.063 0.038 0.042 0.062 -0.234

0.000 0.001 0.096 0.175 0.023

0.24 0.24 0.24 0.24 0.24

0.25 0.25 0.24 0.24 0.24

0.23 0.24 0.23 0.23 0.31

JUPITER JUPITER

1.000 0.500

0.0 0.0

782.5 0.0

515 0

-2.000 0.000

0.000 0.000

0.412 0.000

0.000 o.ooo

0.00 0.00

2.80 0.00

3.96 0.00

Class

= ALINDA

No. of orbits

= 12.12

VENUS VENUS VENUS VENUS VENUS

0.100 0.050 0.020 0.010 0.005

3259.5 814.9 130.4 32.6 8.2

3105.8 813.9 128.2 32.0 8.0

3133 840 121 30 6

0.040 -0.030 0.075 0.083 0.304

0.013 0.184 0.220 0.366 0.295

-0.009 -0.032 0.058 0.064 0.289

0.309 0.175 0.279 0.407 0.310

0.29 0.29 0.29 0.29 0.29

0.28 0.29 0.28 0.28 0.28

0.29 0.30 0.28 0.28 0.25

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

17543.5 4385.9 701.7 175.4 43.9

6781.7 2009.5 403.5 120.1 32.8

6216 1535 219 57 18

0.954 0.963 1.049 1.019 0.836

0.000 0.000 0.000 0.000 0.000

0.087 0.268 0.593 0.712 0.582

0.000 0.000 0.000 0.000 0.004

4.55 4.55 4.55 4.55 4.55

1.26 1.53 2.06 2.66 3.05

1.00 0.93 0.72 0.73 0.85

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

7676.7 1919.2 307.1 76.8 19.2

7998.3 2002.4 309.7 77.0 19.2

7022 1767 263 60 16

0.089 0.083 0.155 0.245 0.181

0.000 0.000 0.006 0.028 0.278

0.130 0.125 0.163 0.248 0.181

o.ooo 0.000 0.004 0.027 0.278

0.20 0.20 0.20 0.20 0.20

0.21 0.21 0.20 0.20 0.20

0.19 0.19 0.17 0.16 0.17

JUPITER JUPITER JUPITER

1.000 0.500 0.200

403.8 100.9 16.1

1953.1 170.2 15.9

848 24 0

-0.710 1.232 2.000

0.000 0.000 0.000

0.789 1.506 2.000

0.000 0.000 0.000

1.44 1.44 1.44

15.30 3.81 1.43

5.11 1.20 0.00

Note.

See Note

to Table

I.

planet (as specified in column I), both the expected number of close approaches by Kessler’s formula and the actual number found in SPACEGUARD with impact parameter 5 D are given in columns 2 and 4, respectively; the relative error &K = ~(NK N&l(NK + Nabs) is in column 5, and the

Poisson distribution probability of the discrepancy NK - Nabs occurring by chance is in column 6. Column 9 contains the Kessler estimate of the probability of collision per IO9 years per asteroid, or more exactly the expected number of close approaches with an impact parameter less than the collision

302

MILANI,

CARPINO, TABLE

D

Planet

Class

= OLJATO

NK

N~rn

No. of orbits

Kb*

AND

III-Continued. EK

nK

0.100 0.050 0.020 0.010 0.005

26799.5 6699.9 1072.0 268.0 67.0

18505.0 6105.8 1064.8 269.1 67.7

16705 5159 853 196 59

0.464 0.260 0.228 0.310 0.127

0.000 0.000 0.000 0.000

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

17417.5 4354.4 696.7 174.2 43.5

13675.7 4319.4 713.3 175.0 43.6

12057 3640 604 134 32

0.364 0.179 0.143 0.261 0.306

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

8796.8 2199.2 351.9 88.0 22.0

6956.4 2011.2 339.4 85.8 21.5

6432 1685 286 69 15

JUPITER JUPITER JUPITER

1.000 0.500 0.200

0.0 0.0 0.0

1642.9 37.0 0.0

= COMET

CKWI

nKm

P,(K)

s(Km)

Pc(obs)

= 5.21

VENUS VENUS VENUS VENUS VENUS

Class

MARZARI

No. of orbits

0.000 0.000 0.000 0.000

0.181

0.102 0.168 0.221 0.314 0.138

0.000 0.000 0.000 0.001 0.042

0.126 0.171 0.166 0.265 0.306

0.311 0.265 0.207 0.242 0.378

0.000 0.000 0.000 0.021 0.077

2513 25 0

-2.000 -2.000 0.000

0.159

5.19 5.19 5.19 5.19 5.19

3.68 4.80 5.16 5.22 5.24

3.40 4.14 4.27 3.94 4.80

0.000 0.000 0.000 0.001 0.042

4.14 4.14 4.14 4.14 4.14

3.21 4.10 4.24 4.16 4.13

2.87 3.56 3.68 3.35 3.28

0.078 0.177 0.171 0.217 0.357

0.000 0.000 0.002 0.036 0.092

0.46 0.46 0.46 0.46 0.46

0.37 0.42 0.45 0.45 0.45

0.34 0.36 0.38 0.36 0.32

0.000 0.000 0.000

-0.419 0.388 0.000

0.000 0.024 0.000

0.00 0.00 0.00

39.33 2.33 0.00

64.22 1.71 0.00

0.208 0.224 0.329 0.310 0.015

0.000 0.000 0.000 0.045 0.564

0.82 0.82 0.82 0.82 0.82

0.85 0.83 0.84 0.84 0.83

0.71 0.69 0.63 0.66 0.89

0.000 0.000 0.011 0.442 0.113

0.55 0.55 0.55 0.55 0.55

0.57 0.55 0.55 0.55 0.55

0.41 0.44 0.41 0.49 0.92

= 4.66

VENUS VENUS VENUS VENUS VENUS

0.100 0.050 0.020 0.010 0.005

4009.3 1002.3 160.4 40.1 10.0

4162.6 1017.9 163 .O 41.0 10.2

3384 813 117 30 10

0.169 0.209 0.313 0.288 0.002

0.000 0.000 0.000

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

2243.5 560.9 89.7 22.4 5.6

2304.5 561.4 89.4 22.4 5.6

1662 446 68 21 9

0.298 0.228 0.276 0.066 -0.464

0.000 0.000 0.435 0.115

0.324 0.229 0.272 0.062 -0.469

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

876.9 219.2 35.1 8.8 2.2

891.1 220.1 34.9 8.7 2.2

817 205 37 11 4

0.071 0.067 -0.053 -0.226 -0.584

0.022 0.178 0.395 0.267 0.179

0.087 0.071 -0.058 -0.230 -0.589

0.006 0.163 0.385 0.263 0.177

0.05 0.05 0.05 0.05 0.05

0.05 0.05 0.05 0.05 0.05

0.05 0.05 0.05 0.06 0.09

JUPITER JUPITER JUPITER JUPITER JUPITER

1.000 0.500 0.200 0.100 0.050

2703.7 675.9 108.2 27.0 6.8

4372.8 858.6 114.3 27.5 6.8

7747 1086 131 27 5

-0.965 -0.465 -0.191 0.001 0.299

0.000 0.000 0.018 0.548 0.333

-0.557 -0.234 -0.136 0.017 0.307

o.ocm 0.000 0.067 0.516 0.326

67.22 67.22 67.22 67.22 67.22

91.87 77.02 70.31 71.77 73.17

164.38 97.01 73.57 54.95 18.67

value D,: 0: = R’* (I + s)

(7)

a function of the physical radius of the planet R’ and of the encounter relative velocity V; the latter is computed from (4).

0.060 0.580

0.010

This expected number P, must be understood in the sense that the number of such orbits escaping collision after T x lo9 years is most likely to decrease as exp(-PJ). In column 11 the same expected number of collisions is estimated by resealing the number of close encounters within D, as

PROJECT SPACEGUARD

found in SPACEGUARD, by the factor Dt/ D2 derived from (5) and (7), with V the encounter velocity (outside the sphere of influence) corresponding to the velocity VMIN at the minimum distance: V2 = V& 2GM’ld~t~. Columns 3, 7, 8, and 10 are explained below. For larger values of D one would expect Kessler’s formula not to match the SPACEGUARD output too closely, and this is indeed the case, even though the average discrepancy is less than 20% for all the inner planets and all the impact parameters (for which the difference is significant at the 0.06 level). The main causes of this difference are three. (A) There are asteroids with perihelia just above the semimajor axis of a planet; these do have shallow close approaches, for values of D such that a( 1 - e) < a’ + D, while Kessler formula predicts none. The asteroids in the class we have named Eros behave in this way with respect to the Earth; Table II shows the Kessler prediction performed for the asteroids in this class and for the timespan over which they remain in this class (as given in our previous paper, Table XIII); the prediction is completely wrong for D 5 0.1, much less so for lower values. The same phenomenon can occur for orbits with aphelion just below a’. (B) Some asteroids are, at least temporarily, protected from close approaches. There are orbits protected from close approaches to a planet by secular perturbations which lift the perihelion far from the orbital plane whenever the argument of perihelion w is near either 0 or 7~, and the eccentricity is maximum; the asteroids with this behavior with respect to the Earth belong, by definition, to the Kozai class. The Kessler’s formula predictions for the objects in this class are therefore far too high, as can be seen from Table II, and more so for small D because the Kozais have no close approach at all at distances less than 0.045 (see previous paper, Table X). This arises from the fact that the Kessler’s formula averages over all the possible values

303

of the angle w for a constant value of e, while e actually changes by a significant amount in phase with the argument 20. In the extreme cases, w librates and therefore an average over the angle w has nothing to do with a time average (this is the case for the asteroid 3040 Kozai, as well as for other orbits in the solar system such as that of the planet Pluto; Milani et al. 1989~). Another protection mechanism can involve locking in a mean motion resonance; if the position of the asteroid and of the planet along their orbits are not independent, it is possible that the asteroid never passes through the nodes of its orbit when the Earth is nearby. Such an orbit belongs to the Toro class; the number of close approaches is much smaller than predicted by any statistical theory, and for small D the discrepancy increases (Table II). (C) as discussed above, Kessler’s formula had the purpose of estimating the probability of collisions, or anyway close approaches at very small D; for a low inclination and/or shallow crosser, the number of close approaches at impact parameters less than D is not proportional to D2 (see Appendix 2), but can grow linearly with D or even slower than D. This can result in a serious overestimate of the expected number of shallow close approaches. However the problem does not disappear entirely even for very small D. For orbits which belong to the same orbital plane of the planet, and/or for orbits with either perihelion or aphelion exactly at a’, Kessler’s formula is singular. It could be argued that such a singularity occurs only over an hypersurface in the element space, that is in a set of measure zero, and the probability for a random orbit to be exactly on the singularity is zero. However, the orbital elements evolve with time, and the curve a(t), e(t), I,(t) can intersect the singularity hypersurface at a special time L If this occurs, the time integral of the probability per unit time P(t) is an improper integral, which actually converges if the singularity hypersurface is crossed in general position (as it is almost

304

MILANI,

CARPINO, AND MARZARl

always the case), e.g., if dqldt f 0 when q = a’ (see Appendix 2). In these circumstances the expected number of close approaches over a finite time span is finite and well defined, but to numerically compute an improper integral by a sum of values at uniformly spaced times is not the right way to proceed. In other words, the singularity is real in that there is a very short time span over which the probability of collision is very high, and the contribution of such events to the overall collision probability over a long interval is large but finite; the problem is the numerical instability of the collision probability computation, which can result in an unreliable estimate of the number of expected collisions. However, it is not obvious-from this purely analytical argument-whether such an instability could result in an overestimate or rather an underestimate of the collision probability. While the existence of these three problems is apparent from the very analytical expression and derivation of Kessler’s formula, the relative importance of each of these effects could not be predicted in advance. Thus we need the results of the comparison, given in Table I for all the orbits in our sample, in Table 11 for the orbits in the “protected” classes, and in Table III for the unprotected Geographos, Alinda, Oljato, and “comet” classes. As we can see, the general trend is an overestimate of the expected numbers, most sensible for the Earth for rather shallow close approaches. The effect of problems (A) and (B) can be estimated from Table II, while problem (C) shows up in Table III, and in the most striking way for the orbits in the classes Alinda and Oljato, which have often a very low inclination (see previous paper, Tables XI and XV). We need to stress again that the statistical methods should be used as such, that is they should not be used to attempt a prediction of the behavior of a single orbit. However, the contribution of low inclination and/or almost tangent orbits to the probability of collision (and close approaches

with small D) is so important that even a single orbit can change the result in a significant way. By inspecting a larger table containing the result of such a comparison asteroid by asteroid, we found that there is indeed such a case of catastrophic failure of the formula (3); Table IV documents such a catastrophe. The culprit is 1982 DB, whose inclination with respect to the orbital plane of the Earth is often low, and actually goes through zero, with the ascending and descending node exchanged instantaneously, at = +63,600 years (Fig. 1). Moreover, the perihelion distance is around 1 AU; after the time +46,000 years, because of resonant perturbations, the eccentricity decreases and the perihelion distance passes through 1 AU while the inclination is very low, thus giving rise to the worst possible case for formula (3), very close to a double singularity. The actual number of close approaches is very high, but no close approach occurs around the time of the zero inclination; as we shall see, this is a very important clue to some dynamical mechanism which had not been detected. As can be seen in Table I, the difference between the prediction in column 2 and the actual number of close approaches in column 4 is about 26,000 for the Earth and for D 5 0.1; 1982 DB contributes to this discrepancy by about 21,700. From Tables III and IV we see that the very large discrepancy in the Alinda class is almost entirely due to this single orbit. What should we do with this misbehaved asteroid? Should we remove this orbit from our sample of 89 and then proceed to proclaim that the formula (3) agrees very well with the results of the numerical integration? We feel that this would be wrong, not only because it amounts to cheating, but also because in this way we would miss the opportunity to understand a very important dynamical phenomenon. If we were to remove 1982 DB (the orbit consistent with the initial conditions of) from our sample, we would also remove its 7024 close approaches to Earth (within 0.1 AU) from our

305

PROJECT SPACEGUARD TABLE COLLISION Planet

D

NK

PROBABILITIES

IV

(KESSLER)-ASTEROID

1982 DB

NICE

Net..

SIC

l-h

(3)

(4)

(5)

(6)

= GEOGRAPHOS

Timespar

1: from

T = -100608.4yr

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

17245.2 4311.3 609.8 172.4 43.1

5871.1 2507.5 631.1 176.1 44.0

6984 3092 816 230 60

0.847 0.329 -0.168 -0.286 -0.328

0.000 0.000 0.000 0.000 0.009

-0.173 -0.185 -0.255 -0.265 -0.307

0.000 0.000 0.000 0.000 0.013

71.20 71.20 71.20 71.20 71.20

23.30 43.03 65.74 72.25 72.31

30.63 57.46 96.19 106.49 113.32

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

3216.4 804.1 128.7 32.2 8.0

2971.6 901.1 133.6 32.4 8.0

2917 890 135 28 11

0.098 -0.101 -0.048 0.138 -0.311

0.000 0.001 0.299 0.265 0.188

0.019 0.012 -0.011 0.144 -0.311

0.161 0.364 0.462 0.254 0.188

1.51 1.51 1.51 1.51 1.51

1.40 1.70 1.57 1.52 1.51

1.35 1.64 1.55 1.24 1.97

(2)

(1)

Class

Class

= ALINDA

(??!a#?)

Tiespan:

from

T = +460OO.Oyr

&Km

HKm

(7)

(8)

to T = : +46000

p=(K)

(9)

.oyr

to T = +100608.4yr

EARTH EARTH EARTH EARTH EARTH

0.100 0.050 0.020 0.010 0.005

11533.2 2883.3 461.3 115.3 28.8

593.5 347.6 142.5 56.5 17.2

40 0 0 0 0

1.986 2.000 2.000 2.000 2.000

0.000 0.000 0.000 0.000 0.000

1.747 2.000 2.000 2.000 2.000

0.000 0.000 0.000 0.000 0.000

167.77 167.77 167.77 167.77 167.77

9.16 19.53 49.04 79.22 97.68

0.33 0.00 0.00 0.00 0.00

MARS MARS MARS MARS MARS

0.100 0.050 0.020 0.010 0.005

462.8 115.7 18.5 4.6 1.2

519.8 122.6 18.6 4.6 1.1

478 108 17 2 2

-0.032 0.069 0.085 0.793 -0.534

0.232 0.254 0.421 0.160 0.322

0.084 0.127 0.090 0.792 -0.540

0.034 0.100 0.414 0.161 0.319

0.54 0.54 0.54 0.54 0.54

0.60 0.57 0.54 0.54 0.54

0.53 0.47 0.47 0.22 0.89

See Note

to Table

Now.

I.

data base. In other words, such a low inclination orbit does not only give a big contribution to the inaccuracy of the formulas such as (3), it also contributes in a significant way to the real probability of close approaches and collisions; from the opposite point of view, such an orbit can have a very rapid evolution as a result of very frequent and deep close approaches. Thus the question arises whether it is possible to improve the formulas such as (3) in such a way that they are not singular any more. The answer to this question is a qualified yes. It is possible to derive a formula which computes the probability of close approaches with impact parameter 5 D, such that the probability is always finite for a finite D, even for zero inclination and tangent orbits. However, in the singular cases it is still true that the probability per unit

area has an infinite limit for D --, 0; as an example, for zero mutual inclination the probability is essentially proportional to D, not to D2, since the only cross section which matters is a segment and not an area. This implies that the probability of close approaches for comparatively large D is not only nonsingular but can be computed in a stable way; for very small D, e.g., for the physical collisions, a very small change in the orbital elements results in a very large change in probability. In other words, the essential problem of the collisional probability calculations is that in the near-singular cases the probability is extremely sensitive to the value of the orbital elements, which therefore need to be know very accurately and with a very high temporal resolution; this problem cannot be removed by a change in the formula. On the other hand,

306

MILANI,

CARPINO,

AND

MARZARI

FIG. 1. Orbit of asteroid 1982 DB from the SPACEGUARD numerical integration (all the Fourier components having a period shorter than 86.2 years have been removed with a low-pass digital filter): (a) mutual inclination with respect to the Earth; (b) distance from the Earth’s orbit along the line of mutual nodes (dashed line: ascending node; solid line: descending node); the vertical lines at t = 63,600 year correspond to a passage through zero mutual inclination, with exchange between ascending and descending node; (c) “comb plot” of the close approaches with the Earth (each close approach is represented by a vertical segment having height proportional to the logarithm of the reciprocal minimum distance); the zero-inclination crossing at t = 63,600 year does not correspond to any actual close approach, because of the Toro type protection mechanism, due to the 8 : 19 resonance with the Earth. The horizontal lines indicate the limits of the sphere of attraction (dashed line), the sphere of activity (dotted line) and the Hill’s sphere (dashed-dotted line) according to Chebotarev ( 1967).

the computations of the frequency of shallow and moderately deep close approaches can be done in a more stable way, which avoids the 1982 DB type catastrophic failures . Although most details of the derivation are given in the Appendix 2, we give here the geometric idea behind the nonsingular formula. The probability (3) is computed by multiplying the averaged spatial density (1) where the planet is (that is, at Y = a’ and

/3 = 0) by a cross section A and a velocity V. However, if the cross section A is finite, the total probability should rather be the average of S(r, /?) over the volume swept by the cross section A moving with velocity V. Let us simplify this integral by orienting the cross section area A in such a way that it lies in the plane Y, z = r sin p; let us also neglect the dependence of the relative velocity V from the position (Y, p); these are of course approximations, which are some-

PROJECT SPACEGUARD

307

what justified in Appendix 2. Then the nonsingular probability can be approximately computed as

where the area A is defined by (r - u’)~ + r2 sin2 /3 5 D2. Of course if the area A is well inside the region of the r, z space were the average spatial density in nonzero (that is, well inside the trapezoidal region q 5 r % Q, IzI 5 r sin I, see Fig. 2), then the average over the area A will not be very different from the value in the central point Y = a’, z = 0; that is, the result of (8) and the simpler formula (3) will be almost the same. The difference becomes relevant when the area A pokes outside the trapezoid; then the average performed by the integral is very different from the value in the central point; in the singular cases for formula (3), when the center is on one of the sides of the trapezoid, the integral is always convergent (as an improper integral) even if the value at the central point is actually infinite. This is

FIG. 2. Representation of the spatial density function S(r, p) of the Opik-Kessler method, obtained by plotting 8000 positions of the asteroid corresponding to random values of the angle variables (mean longitude, longitudes of the node and the perihelion). In the Kessler method the integral of the function over a disk is approximated by the product of the surface of the disk by the value of the function at the Earth’s position (the solid point at r = 1, p = 0). A nonsingular formula for the close approach probability can be obtained by actually computing the integral (see Appendix 2).

FIG.3. Integral of the spatial density function S(r, Z3) over the close-approach region (r - al)2 + r2 sin2 p 5 D2 (which is proportional to the collision probability) as a function of the impact parameter D, for a lowinclination crosser. The collision probability is proportional to Dz only for small values of the impact parameter (D < r sin I); at larger distances the dependency becomes linear. proven in Appendix 2, which also contains the algorithm to be used in a semianalytical computation of the integral (8). Formula (8) achieves two interesting goals: it is nonsingular and indeed very suitable for numerical computation, and the dependence of the probability upon the limiting value of the impact parameter D significantly deviates from the quadratic law which is assumed by formula (3), at least for low inclination (that is, for sin ZM 5 D/a’) and shallow/almost crossing cases (Fig. 3). To evaluate its performance in a comparison with the SPACEGUARD close approach data base we need to look in Tables I-IV at columns 3 (expected number of close approaches), 7 (relative error with respect to actual number found by numerical integration), 8 (Poisson probability of the discrepancy occurring by chance), 10 (expected number of collisions per lo9 years), which all refer to the nonsingular formula (8). From Table I we find that the overall agreement with the output of the numerical integration is almost too good to be believed (differences smaller than 10% for all

308

MILANI.

CARPINO,

AND

MARZARI

the terrestrial planets and down to the entry (two-body approximation) value a = 1.7788 for D I 0.02). It is worth stressing again AU, thus 1982 DB was classified as Alinda. that in the formulas used in all the Tables IHowever, at CI = 1.7801 AU there is the IV we have entirely neglected the eccennominal value for the 8 : 19 resonance with tricities of the planetary orbits (even for the Earth; this possibility had been disMars!). And indeed we do not believe that missed because of the high degree of the the formulas such as (8) can work for all resonance (the degree is 19 - 8 = 11, and is kind of orbits. Inspection of Tables II and also the degree in e of the lowest degree III shows that there is still a significant and term containing the resonance in the develostensibly systematic underestimate of the opment of the perturbing function). After number of close approaches for the Geo- the discrepancies reported in the discussion graphos class, more or less compensated by above pointed to some protection mechaoverestimates for the other classes. Some nism, the range of values for the semimajor of these overestimates correspond to probaxis was computed more accurately and lem (B), i.e., to protection mechanisms found to be 1.7783 < a < 1.7817; moreover, (such as for the Kozai and the Toro); probthe critical argument 19h - 8h’ - 116 was lem (A), that is shallow crossing, is also found to be in a remarkably stable libration handled much better than with formula (3) (A, A’ are the mean longitudes of the aster(e.g., for the Eros); other discrepancies are oid and the Earth, respectively; (3 is the hard to understand. Even for the ap- longitude of perihelion of the asteroid); proaches to Jupiter at an impact parameter therefore, 1982 DB has to be classified as a 50.5 AU, which hardly qualify as “close” Toro. A very similar behavior is found for approaches, the nonsingular formula per- 1982 XB. forms well, and allows one to detect the However, the virtues of the nonsingular protection mechanism preventing the formula are even more apparent when it is Alinda class from approaching Jupiter. used to tackle another problem. Following As for our most pathological example the schematic presentation of an abstract 1982 DB, the comparison of columns (3) statistical theory we gave in Section 1, we and (4) in Table IV tells an interesting have tested assumption (2), namely that the story. For the time span for which the orbit number of close approaches can be preis a Geographos-that is, no resonances in dicted by a probabilistic computation once mean motion-the nonsingular formula the orbital elements are known. However, gives a fairly good agreement. For the time apart from the cumbersome procedure of span over which the orbit had been classi- performing a numerical integration of all fied as an Alinda, the disagreement is still the orbits, how can we assume we know the serious in relative terms (although the abso- time series for the orbital elements, for orlute numbers in column (3) are not as large bits so plagued by both mean motion and secular resonances that essentially any simas in the singular computation in column plified dynamical model is bound to be un(2)). As we discussed in our previous paper, Section 4, the attribution of each orbit to a reliable? That is, what is the meaning of dynamical class has been obtained by in- assumption (I), the computability of the orspecting one by one and comparing some- bits? A very simple test of the sensitivity of thing like 1,400 plots; such a cumbersome the results-in terms of close approach and the method used procedure cannot be completed without er- collision frequencies-to rors and oversights. After +46,000 years in the computation of the orbits can be performed as follows. Let us assume that we the semimajor axis of 1982 DB remains locked around an average value a = 1.78 use the crudest possible orbital model, AU, with small oscillations. The 5 : 1 reso- namely we only use the osculating orbital elements at a given time t, compute the nance with Jupiter occurs at a nominal

PROJECT

probability per unit time P(r) and then estimate the probability over the entire interval of length T by Z’(t) x T. How much would this estimate change with time? Fig. 4 shows the result of such a computation for all the orbits in our sample; Fig. 4a uses the singular formula (3), while Fig. 4b uses the nonsingular formula (8), computed as described in Appendix 2. On the ordinate axis is the total expected number of close approaches to the Earth within 0.1 AU over

-1 10%

-5s10’ T,me ,Y,

0 s Id” 5:PACEC”ARDlnlegrallon

309

SPACEGUARD

the integration time span of SPACEGUARD; for comparison, the actual value (from Table I, column 4) is 126,199. The most striking feature of Fig. 4a is the sharp peak with a value above 8 x 106; this peak is entirely due to the well known rascal 1982 DB, which happens to be almost in the ecliptic plane and almost tangent to the Earth’s orbit. Therefore the values of the singular prediction are totally unreliable unless the orbital elements are known over a

105 rime ,Y,

SPACECUmD lnlegratlo”

FIG. 4. Estimate of the total number of close approaches to the Earth (within a maximum impact parameter D = 0.1 AU) for all the 89 orbits in our sample, obtained by applying Kessler’s formula to the filtered orbital elements (from numerical integration), as a function of time. Plot a refers to the original formulation of the Kessler method, while plot b shows the results obtained with the improved formula proposed in Appendix 2. The fantastic peak at r = 63,600 year (in plot a) is entirely due to the asteroid 1982 DB. Plot b shows a systematic effect: the expected number of close approaches appears to have values significantly above average for some time around t = 0. This can be explained as the result of observational selection effects (asteroids are more likely to be discovered when their close approach probability with the Earth is higher). This interpretation is confirmed by plotting separately the expected number of close approaches for the faintest (with absolute visual magnitude Vabs 2 18, plot c) and the brightest objects (Vabr < 18, plot d). These plots show not only that the upside-down Vshape of b is entirely due to the asteroids with V abs2 18, but also that the average number of close approaches per asteroid is significantly higher for fainter objects.

310

MILANI,

CARPINO,

long enough time to allow averaging out the almost singular cases. On the contrary, the nonsingular computation of Fig. 4b shows that the expected number, as computed by using the orbital elements at a single time, is never off by more than 20%. It is also interesting to notice that the expected number of close approaches appears to have values significantly above average for some time around t = 0. This could be interpreted as the result of the observational selection effect by which an asteroid is more likely to be discovered when it has frequent close approaches (e.g., because the orbit is almost tangent to the one of the Earth). This interpretation can be confirmed by plotting the expected number of close approaches, computed with the same nonsingular formula, only for the asteroids in our sample with an absolute visual magnitude ~18 (Fig. 4~); the similar plot for asteroids with absolute magnitude ~18 does not show the same signature (Fig. 4d). We conclude that a large fraction of the smallest asteroids are discovered when they have an exceptionally high frequency of close approaches. This result not only shows the advantage of the nonsingular formula, it also confirms that statistical theories can be used to predict the frequency of close approaches, even if the orbit computation is very rough, because the result depends very little upon the orbital evolution, provided there are enough orbits to average over. In other words, we have obtained an empirical confirmation of the ergodic principle, by which an average over enough orbits and an average over a long enough time gives similar results. We would like to give a much weaker statement as far as the predictability of the collision rates is concerned. As it is clear from the derivation, the stability of the nonsingular formula arises from the fact that it computes close approaches up to a finite impact parameter D; the smaller D, the less stable the formula; at the limit for D + 0, formula (8) becomes formula (3) and is again singular. Thus at small D, a figure

AND

MARZARI

such as Fig. 4b shows much larger oscillations; for D 5 0.001 the peak produced by 1982 DB is almost as bad as the one of Fig. 4a. For very close approaches, and more so for physical collisions, the probability may never become infinite, but it can be very high indeed for peculiar orbital configurations such as the one occurring for 1982 DB at = +63,600 years. This also implies that resonant protection phenomena, such as the one experienced by 1982 DB, can be important even if they are not stable for a very long time, if they occur at these critical times. As the convergence of the values of column 10 and column 9 shows, the expected number of 271 collisions per IO9 years for 1982 DB in the Geographos phase does not show any disagreement between the two formulas; moreover, the real value as obtained from the actual close approach count may even be larger. Thus very few exceptional cases can contribute disproportionately to the overall collision probability. To say it even more bluntly, the discovery of a single new low inclination Aten may increase the expected impact rate by a significant amount, possibly larger than the uncertainty of the present calculations. On the other hand a new computation of the orbit of 1982 DB in which, by virtue of the combined effect of the many close approaches to the Earth, this dangerous neighbor gets thrown out of the inner solar system may change the result in the other direction by an equally significant amount. We can conclude this assessment of the statistical formulas of the Kessler type by stressing again that the statistical formulas should not be expected to work outside the rules of statistics; e.g., they can be good enough to predict a large number of shallow close approaches, even for a single asteroid, but they are not likely to predict very accurately the frequency of a very rare event such as a collision. For these, the best estimate can be obtained by averaging over the largest possible data set (both in terms of number of orbits and of time span); the most likely value appears (from Table I,

PROJECT SPACEGUARD

columns 9, 10, and 11) to be -3.5 collisions with the Earth per lo9 years per asteroid. 4. COLLISION PROBABILITIES WETHERILL METHOD

BY

Wetherill (1967) improved the collision probability formulas previously computed by opik (1951) by taking into account the eccentricity e’ of the target body, and also by using an explicit description of the geometry of the orbit crossing as a function of the change in the orbital orientation angles, especially the argument of perihelion w of the asteroid. These formulas can be used in two ways. By assuming that the argument of perihelion is changing uniformly with time, i.e., that averaging with respect to w is equivalent to averaging with respect to time, a long term average probability of collision can be computed from the present value of a, e, I. If on the contrary the orbital elements time series-including ware available from either number integration or some secular perturbation theory, then an explicit computation of the collision probability can be given for each node crossing.

311

Wetherill (1967, formula (15)) gives the probability (or better expected number) of close approaches with an impact parameter 5 D per orbital period of the asteroid as Pw =

r sin (Y’ 4at2 +FF

X

DV (9) d/v2 - (V, cos (Y’ + VT sin CL)* where V is the relative velocity, with components VR in the radial direction and V, in the direction along the orbit of the planet, (Y’ is the angle between position and velocity of the planet, r the distance from the Sun, all computed at the time at which the two orbits touch. However the cumulative probability over a long time span T is not just the product Pw x 27rTln, with n the asteroid mean motion; the probability has the average value Pw only over the time span in which the distance d between the two orbits is smaller than D, and is zero otherwise. For the distance between the two orbits, Wetherill (1967, formula (4)) gives the formula X

tan2 (Y’ sin2 ZM sin2 ff ii d2 = d (1 - cotg (Y’ cos ZM tan (Y)~sin2 (Y’ + tan2 ff sin2 ZM where dN is the distance along the mutual node, CYis the angle between position and velocity, this time for the asteroid, at the touch time, and ZM the mutual inclination. Formula (IO) is derived by the approximation of the two orbits-in a neighborhood of the mutual nodal line-with two straight lines; it is of course a good approximation only provided the mutual inclination is not too small, i.e., provided r sin ZM > D. We have used Wetherill’s formulas for the purpose of testing the assumption (2) of the introduction, namely that close approaches occur at random whenever the geometrical relationship between the two orbits is such that the minimum distance is small enough. To this purpose we use, for

, (10)

each asteroid and for each planet, a complete set of orbital elements (apart from the anomalies) for both orbits at each time t, then compute the distance between the two orbits d(t) by formula (10). In this way we derive an explicit list of times t in which the distance d is smaller than D; the sampling of the time series being frequent, we can transform this list of times into a list of intervals [tr, f2] in which d < D; each of these we call a crossing. However, a real node crossing occurs when the distance d passes through the value 0 and changes sign (the distance dN , thus d, can be considered signed, with a positive value when the asteroid orbit is outside the planet orbit at the mutual node). The simplest case, a deep

MILANI,CARPINO.AND

312

crossing (DC), occurs when the values +D, 0, -D are crossed each once, either in this order or in the order -D, 0, +D; then the time to at which d = 0 can be found, e.g., by interpolation, and the values of (Y, CX’, V, r etc. can be computed at co, giving by (9) a value of PW. The expected number of close approaches (within D) during the deep crossing is then P = (tz - t,) $

Pw

(II)

and the total expected number for one asteroid is the sum of all the values of P over all the deep crossings occurring in the entire time span. These expected numbers can then be compared with the actual numbers of close approaches in the data base derived from the numerical integration, either for each node crossing or for the total over the numerical integration time span. As we did in Section 3, we need to stress that such a comparison has its intrinsic limitations. Wetherill’s formula is meant to compute the probabilities for either collisions or very close approaches, and again we need to use it to estimate the number of close approaches up to somewhat large impact parameters D, otherwise the Poisson statistical significance of the comparison with the SPACEGUARD data would be low. In Table V we have computed Wetherill’s formula for D = 0.1 AU, and in Table VI for D = 0.02 AU. On each line of these tables we give the crossing type (see below) in column 1, the number of crossing of each type in column 2, the expected number of close approaches by Wetherill’s formula, and the number found in SPACEGUARD (columns 3 and 4), the relative error and the Poisson probability of the discrepancy occurring by chance (columns 5 and 6), and finally in columns 7 and 8 the expected number of collisions per lo9 years per asteroid obtained by resealing the numbers in columns 3 and 4 to the collision impact parameter (7). The tables indicate a rough agreement; however, for moderate to large D the predictions based upon formu-

MARZARI

las (9), (IO), and (11) cannot be accurate because they fail to address the possibility of more complicated geometrical configurations. This can be understood by analyzing the different possible types of node crossings (see Fig. 5). If -D < d < + D for tl < t < t2 a generalized crossing occurs, but it is by no means obvious that there is one and only one time to in the interval in which d(to) = 0. If there are two times at which d = 0 the crossing is shallow (SC), if there is no such time it is an almost crossing (AC), and in both cases d(tJ and d(t,) have the same sign; if the number of times at which d = 0 is larger than 2 we have a multiple crossing (MC). Over a finite time span another case can occur: if the time span over which we have data ends (and/or begins) while IdI < D, the type of crossing cannot be classified (unknown crossing (UC)). Wetherill’s formula, in its original formulation, applies only to deep crossings, essentially because it is obtained by a geometrical construction which is correct at the limit for D + 0. Thus, if we use Wetherill’s formula as it is, we need to restrict the comparison with the close approach data base to the close approaches occurring during deep crossings; this is done on the lines marked DC in the tables. Even at D = 0.1 AU, the prediction of the total number of close approaches for all the 87 orbits (excluding the ones going hyperbolic) agrees very well with the SPACEGUARD data. Of course this average agreement hides substantial disagreements for single orbits; the most obvious example is given by the Toro class of asteroids, which avoid close approaches to the Earth-while node crossings are occurring-by means of mean motion resonances. This effect is masked in the total by a compensating underestimate of the number of close approaches for the Geographos class. The class averages and totals, for D = 0.02 AU, are given in table VII (if a transition between classes occurs at a crossing, the data belonging to that crossing are not used; the collision probabilities are not

313

PROJECT SPACEGUARD TABLE

V

COLLISION PROBABILITIES (WETHERILL)-ALL crossing

type

(1) Planet

Kb.

CW&

HW&

P,(Wet)

Pc(4

(2)

(3)

(4)

919

50242.6 11225.9 5024.0 6136.7 72629.2 69081.9 0.0 0.0 141711.1

51200 11150 4961 5046 72357 14464 2597 28351 117769

-0.019 0.007 0.013 0.195 0.004 1.307 0.000 0.000 0.185

0.000 0.239

0.641 0.175

0.189 0.000 0.157 0.000 0.000 0.000 0.000

0.112 1.007 0.955 0.000 0.000 1.961

0.668 0.164 0.077 0.069 0.977 0.277 0.033 0.425 1.712

54568.2 23690.4 6492.1 4501.5 89253.2 38232.5 0.0 0.0 127485.7

53906 24540 6150 3756 88352 11631 3186 22921 126090

0.012 -0.035 0.054 0.181 0.010 1.067 0.000 0.000 0.011

0.002 0.000 0.000 0.000 0.001 0.000 0.000 0.000 0.000

0.990 0.610 0.145 0.117 1.862 0.670 0.000 0.000 2.532

0.986 0.604 0.130 0.078 1.797 0.373 0.057 0.591 2.819

37491.5 6005.1 3978.8 1401.9 48877.3 8902.4 0.0 0.0 57779.8

37960 6109 4054 1534 49657 4186 2237

-0.012 -0.017 -0.019 -0.090 -0.016 0.721 O.OC0 0.000 -0.131

0.132 0.022 0.016 0.005 0.175 0.030 0.000 O.OW 0.205

0.133 0.022 0.016 0.005 0.177 0.015

(5)

03)

(7)

(8)

164 119 70 1272 1269 2541

= EARTH

Deep cross Shallow crow Multiple cross Unknown cross DC+SC+MC+UC Almost cross Not at cross Not at nodes All types Planet

NW,

= VENUS

Deep cross Shallow cross Multiple cross Unknown cross DC+SC+MC+UC Almost cross Not at cross Not at nodes All types Planet

#Cr.

ORBITS

1873 320 181 86 2460 1524 3984

= MARS

Deep cross Shailow cross Multiple cross Unknown cross

DC+SC+MC+UC Almost cros8 Not at CT088 Not at nodes All types

3948

151 131 82 4312 497 4809

0.089

0.115 0.000 0.000 0.000 0.000 0.000 0.000

0.036 0.236

Note. Comparison of the expected number of close approaches computed according to Wetherill’s formula for the different types of crossing (column I) with the results of the SPACEGUARD integration. For each type (for the definition see Section 4 and Fig. 51, column 2 gives the total number of crossings detected in the SPACEGUARD orbital element tiles; column 3 gives the number of close approaches Nwe, obtained by adding the results for each crossing computed with Wetherill’s formula (Eq. (9) and (1 I )); the corresponding number N ,,b, recorded during the integration is indicated in column 4. Columns 5 and 6 give the relative discrepancy eWefwith respect to the observed value and the Poisson probability IIwer of a discrepancy equal or greater than the observed one. Columns 7 and 8 give the expected number P, of physical collisions per IOyyears resealed from Nwet and Nob,. Maximum impact parameter D = 0.1 AU. Total number of orbits = 87.

given, because the partitioning of the time span makes them less significant). The absence of the Kozai class from this table has to be rated as a success of the method: indeed the Kozais have no generalized orbit crossing with the Earth at D = 0.02; the

predicted number of close approaches within this value of D is zero and is in perfeet agreement with the numerical results. This comparison restricted to the case in which the Wetherill’s formula works is not enough; it amounts to saying that the for-

314

MILANI,

CARPINO, TABLE

AND

MARZARI

VI

COLLISION PROBABILITIES (WETHERILL)-ALL crossing

type

(1) Planet

(3)

03)

1514 143 58 11 1726 362 2088

3358.1 291.5 184.3 59.6 3893.4 771.2 0.0 0.0 4664.6

3328 329 218 76 3951 424 133 158

0.009 -0.121 -0.168 -0.242 -0.015 0.581 0.000 0.000 0.000

0.306 0.014 0.006 0.023 0.176 0.000 0.000 0.000 0.488

1.139 0.118 0.096 0.020 1.373 0.327 0.000 0.000 1.700

1.137 0.144 0.123 0.025 1.430 0.234 0.045 0.053 1.761

2779

3950.2

245 76 21 3121 602

864.2 270.1 47.6 5132.2 728.5 0.0 0.0 5860.7

3878 598 295 31 4802 459 191 370 5822

0.018 0.364 -0.088 0.423 0.066 0.454 0.000 0.000 0.007

0.127 0.000 0.062 0.007 0.000 0.000 0.000 0.000 0.309

2.081 0.625 0.171 0.028 2.905 0.368 o.oQo 0.000 3.273

2.088 0.400 0.193 0.013 2.693 0.349 0.115 0.300 3.457

2081.2 233.8 86.9 7.8 2409.6 168.4 0.0 0.0 2578.0

2004 245 78 6 2333 115 77 71 2596

0.038 -0.047 0.108 0.258 0.032 0.377 0.000 0.000 -0.007

0.046 0.220 0.185 0.341 0.060 0.000 0.000 0.000 0.356

0.186 0.022 0.008 0.001 0.217 0.015 o.cKIo 0.000 0.233

0.178 0.023 0.007 0.000 0.209 0.010 0.007 0.006 0.233

3723

= MARS

Deep cross Shallow cross Multiple cross Unknown cross DC+SC+MC+UC Almost cross Not at cross Not at nodes All types Notr.

(2)

= EARTH

Deep cross Shallow mom Multiple cram Unknown cross DC+SC+MC+UC Almost cross Not at cram Not at nodes All types Planet

NV.t

= VENUS

Deep crow Shallow crces Multiple cross Unknown cross DC+SC+MC+UC Almost cram Not at cross Not at nodes All types Planet

#Cr.

ORBITS

See Note

to Table

4394 121 55 16 4586 243 4829

V. Maximum

impact

parameter

mula works when it is geometrically correct, that is the assumption that the positions of the asteroid and the planet on their orbit are uncorrelated is sound, apart from the occurrence of a Toro state. However, as it can be easily appreciated from the values in column 4, the number of close approaches occurring in geometrical configurations of the two orbits different from deep crossing is substantial; a model of the orbital evolution of a shallow crossing asteroid which does not take into account the effect of the close approaches occurring at shallow crossings is wrong, and more so

D = 0.02 AU.

Total

number

of orbits

= X7.

because a shallow crosser has on average many more close approaches than a deep crosser. Thus we decided to generalize Wetherill’s formula to the shallow crossing case, by splitting each SC into two generalized crossings separated by the time at which the distance IdI has a local maximum. Since each of the two has one time at which d = 0, we can apply Wetherill’s formula twice. For a multiple crossing we use the formula without changes, by using the last time at which d = 0 to compute (Y, V, etc. For an unknown crossing we use as to either a time at which d = 0, if such a time

PROJECT SPACEGUARD

FIG. 5. Schematic representation of the different types of crossings: deep crossing (DC), almost crossing (AC), shallow crossing (SC), multiple crossing (MC), and unknown crossing (UC). The classification depends on the behavior of the signed distance between the orbits d (= local minimum distance between the orbits with the sign of the distance along the mutual nodal line) with respect to the crossing zone (here defined as the stripe IdI < 0.1 AU). The solid points mark the times at which the orbital elements are interpolated for the computation of Wetherill’s formula for collision probability. A crossing is classified as “unknown” when it is not completely contained within the integration timespan; the correct classification would require the knowledge of the nodal distance outside the integration interval (dashed lines).

exists, or the time at which IdI is minimum. In all these cases we have introduced some further approximations in Wetherill formula; by reviewing the derivation of the formula we find that the corrective coefficient 7r/4 in formula (7) of Wetherill (1967) is not really justified in all these generalized cases, but the correct coefficient would not be very different. The lines marked SC, MC, UC, and the total DC + SC + MC + UC in the Tables V, VI, and VII have been obtained with this procedure, and again the total shows a very good agreement in the sum over all the orbits, which results from the compensation of different effects (such as the Toro-type protection mechanism). The largest discrepancy in Table VIIapart from the Toros-is in the shallow crossing case for the Alinda class, and is due to the wrong classification of 1982 DB;

315

it would not appear in a further iteration of the classification procedure (that is, including the result of this paper). Even after this generalization, the comparison is not satisfactory in that almost one-third of all the close approaches to the Earth within 0.1 AU, and over a sixth of those within 0.02 AU are not included in the comparison. There are many close approaches occurring at almost crossings, and there are close approaches which do not occur at any crossing at all. The latter are due to two troublesome effects: if the inclination is low, that is, r sin I,,,, < D, a close approach within D can occur at a large angle from the mutual node; moreover, the formula (10) is singular for zero inclination; since the two orbits are approximated by straight lines, the distance between the two is almost always zero when the two are coplanar. Thus, a close approach can occur at such a large angle from the mutual line of nodes that it is impossible to attribute it to either the ascending or the descending node; this type of close approaches is entered in Tables V and VI on the lines marked “Not at nodes.” On the other hand, a close approach may occur at such a distance from the node that formula (10) fails to give a correct estimate of the distance between the orbits at that angle; that is, a close approach within D can occur when the distance between the orbits, as computed from (lo), is larger than D; this case is entered as “Not at crossing.” We do not know any simple way to generalize Wetherill’s formula in such a way that these cases are handled correctly. The only way to compare the data for all the close approaches in our data base with Wetherill’s formula is to use a very crude approximation, namely to use the same formula even for the almost crossings. That is, the values of CY,V, etc used in (9) are computed at the time to at which IdI is minimum, even though the actual crossing d = 0 does not occur at all. In this case the corrective factor 7r/4 of formula (7) of Wetherill(1967) is not correct, and the error introduced can

316

MILANI,

CARPINO,

AND

TABL,E CLOSE Crossing

APPROACHES

type

#Cr.

(1) Clam

Nwet

(2)

(3)

(4)

(5)

1143 41 21 10 1215 126 1341

3130.1 491.2 188.7 28.1 3838.1 363.1 4201.1

3330 520 224 16 4090 396 4486

-0.062 -0.057 -0.171 0.549 -0.064 -0.087 -0.066

0.000 0.093 0.005 0.010 0.000 0.041 0.000

39 4 1 44 44

290.7 104.9 11.4 407.0 407.0

9 0 12 21 21

1.880 2.000 -0.055 I.804 1.804

0.000 0.000 0.463 0.000 0.000

21 21

57.5 57.5

0 0

2.000 2.ooo

0.000 0.000

664 159 42 2 067 363 1230

al.0 222.2 35.8 13.8 352.8 265.0 617.8

66 36 25 6 133 37 170

0.204 1.442 0.356 0.786 0.905 1.510 1.137

0.050 0.000 0.037 0.016 0.000 0.000 0.000

544 21 7 9 581 48 629

370.7 34.8 31.6 5.8 442.8 35.0 477.8

391 32 32 9 464 23 487

-0.053 0.085 -0.014 -0.439 -0.047 0.414 -0.019

0.140 0.355 0.493 0.129 0.151 0.021 0.327

385 17 3 405 43 448

54.2 2.6 0.8 57.5 4.0 61.5

57 2 0 59 1 60

-0.051 0.252 2.000 -0.025 1.200 0.025

0.368 0.524 0.452 0.441 0.092 0.456

(6)

= OLJATO

DC SC MC UC DC+SC+MC+UC AC DC+SC+MC+UC+AC Class

(W~THERILL)

= ALINDA

DC SC MC UC DC+SC+MC+UC AC DC+SC+MC+UC+AC Class

EARTH

= EROS

AC DC+SC+MC+UC+AC Class

THE

= TORO

DC SC MC DC+SC+MC+UC DC+SC+MC+UC+AC Class

VII

= GEOGRAPHOS

DC SC MC UC DC+SC+MC+UC AC DC+SC+MC+UC+AC Class

WITH

MARZARI

= COMET

DC SC MC DC+SC+MC+UC AC DC+SC+MC+UC+AC Notr. See Note AU. Total number

to Table of orbits

V. Maximum = X7.

be large. The only reason why this ruthless procedure can give apparently good results is that the error is an overestimate of the likely number of close approaches, and this

impact

parameter

D = 0.02

large positive error compensates the negative error resulting from the low inclination cases. Thus the line marked “All types” in Tables V and VI shows a surprising agree-

PROJECT SPACEGUARD ment, with relative errors below 1% for Venus, Earth, and Mars at D = 0.02 AU; even at D = 0.1 AU the agreement for the Earth is very good, while for Venus and Mars the weight of the almost crossing and low inclination cases is too large. This type of results is a good illustration of a phenomenon we have encountered many times in this kind of computations: when the uncertainties and approximations made necessary by the difficulty of the problem are both large and numerous, the different sources of error somewhat compensate and leave a surprisingly good overall result. For the Earth at D = 0.1 AU, the error in the AC case is = +26,600 and the error due to the “Not at node” and “Not at crossing” cases is = -26,100, giving an almost perfect compensation. The conclusions of this comparison are quite surprising. We had expected Wetherill’s formula to be somewhat more accurate than Kessler’s formula in predicting the number of close approaches; this is because Wetherill’s formula uses more orbital information, such as the value of w(t). Although we have obtained an apparent good agreement by fiddling with the formula and the conditions under which it is applicable, the fact is that the Wetherill method gives either no information at all or very unreliable information on the behavior of both almost crossing and low inclination asteroids. As in the case of the singular Kessler’s formula, one could argue that the method should only be used at very low values of D, where such problems are less frequent and/or less severe (as can be easily seen by comparing Table V with Table VI); however, this means that the method, when used for small D, gives results critically dependent on very small changes in the orbital elements, and would require both high accuracy and high time resolution in the orbital element time series. Moreover, the use of a statistical theory should not be limited to the computation of the number of very close approaches; it should also provide information on the or-

317

bital evolution. When the orbital evolution of an almost crosser is studied, to assume that no close approaches occur and that the semimajor axis is essentially constant would be a very serious mistake (see Section 5). We have to acknowledge that we have given a rather unfair treatment to Wetherill’s formula with respect to Kessler’s formula, in that we have gone to some effort to improve the latter for both the singular cases. Could Wetherill’s formula be improved by some kind of removal of the singularity such as the one discussed in Appendix 2? This might well be possible, especially for the almost crossing case; some semianalytic formula could be devised to replace the already mentioned factor n/4, which is the result of the correction from the area of the square with side 20 to the cross section rD2, with a suitable integral applicable in the shallow and almost crossing case. However the low inclination case looks more difficult, since the approximation used in formula (IO)-the distance computed by the formula for the distance between straight lines-breaks down completely. This might only mean that we have had neither the time nor the energy to try hard enough to improve Wetherill’s formula, but we are left with the impression that an improved formula would be so complicated that it might not be efficient. On the other hand, when Wetherill’s formula applies-i.e., in the deep crossing and even in the shallow crossing case-it is so accurate and reliable that it can be used for a purpose very different from the one it was invented for. We routinely use Wetherill’s formula, coupled with our fully automated crossing analyzer’s “Toro detector” in a way which is best illustrated by an example. The identification of a Toro requires three steps: if the semimajor axis plot shows oscillations with the periods typical of mean motion resonances (a few hundred years; see e.g., Fig. 6a), a table of resonant values is consulted; if the value of a corresponds to a possible resonance with the

318

MILANI,

CARPINO,

Earth, the nodal distance plot (e.g., Fig. 6b) is consulted to check that node crossings occur. If this is the case, then the “comb plot” (e.g., Fig. 6c) is used to assess whether a protection mechanism is active at a number of node crossings. In most cases this procedure is very effective (see, e.g., the cases of 1865 Toro and 3103 1982 BB, Figs. 6 and 7 in our previous paper); however, there are always less obvious cases which can be overlooked (not to speak of the fact that in doing this work the human eye gets tired very quickly, and can simply fail to see what is visible). The example of Fig. 6, 1980 PA, is such an oversight, and the printout from our Wetherill formula computer code (that is Table VIII) clearly shows a protection mechanism effective over two deep crossings and two shallow crossings, as well as over some almost crossings. 1980 PA is a Toro over the entire time span in which its perihelion is low enough to allow very close approaches and is never a Geographos; it becomes an Eros later, when the eccentricity is lower, still being in the resonance for some time. A further confirmation can be obtained by looking for a critical argument in libration; 3X’ - 8k + 56 plotted in Fig. 6d is such a critical argument. 5. CORRELATION

ANALYSlS

The main difficulty of the dynamics of planet-crossers is due to the fact that such orbits are simultaneously perturbed in two very different ways: long range perturbations, with continuous effects accumulating slowly with time, especially whenever some resonances (either mean motion or secular) are encountered, and occasional short range interactions with an almost impulsive effect. We would like to find an algorithm which allows one to separate the two kinds of perturbations and to give a rigorous meaning to qualitative statements such as “the orbital evolution of a Geographos is determined by the close approaches to Earth and Venus.” The method we are going to discuss in this sec-

AND

MARZARI

tion is correlation analysis; it provides a comparatively simple algorithm and a quantitative measure of the dependence of the orbital evolution upon the immediate effects of close approaches. We shall first discuss the effects on the semimajor axis a of the asteroid orbit. Let us assume that for some span of time there are no close approaches. a(t) will still undergo short periodic perturbations, with periods of the order of the orbital period; these can easily be removed by a digital filter such as the one we have used to generate the postfiltered time series (see Section 2 and Appendix 1). The filtered a(t) will undergo resonant perturbations and/or secular perturbations. The former arise only when the mean motion n is in a ratio with the mean motion IZ’ of some perturbing planet very close to a fraction p/q, with p, q small integers; the effect is roughly inversely proportional to the smafl divisor and the period 2rl(qn - pn’) qn - pn’, must be longer than -90 years for these effects not to be wiped out by the filters. Contrary to widely held but incorrect beliefs there are also secular perturbations on the semimajor axes, and there are known formulas to compute at least their order of magnitude (Milani et al. 1987); however, the secular effects are of the second order in the mass of the perturbing planet, and they are only a minor effect with respect to the large changes occurring in planet-crossing orbits. Thus the filtered semimajor axis a(t) should either remain almost constant or oscillate around a resonant value. When a close approach occurs, a sudden jump to a new value should be apparent from the time series a(t). If a comparatively long time span without close approaches is followed by a deep node crossing with the corresponding shower of close approaches, including some very deep ones, then the semimajor axis plot looks like the example in Fig. 7a. These orbits “regular between node crossings” can be classified as Geographos at a glance, and we do not really need confirmation of the fact that the semi-

PROJECT SPACEGUARD

319

L.802- L.60i; 1.40m w 1.208 2 1.00IO , 0.805 = 0.602 ; 0.405 ; 0.200 0.00-

FIG. 6. Orbit of asteroid 1980 PA from the SPACEGUARD numerical integration (filtered data): (a) semimajor axis; (b) distance from the Earth’s orbit along the line of mutual nodes (dashed line: ascending node; solid line: descending node); (c) close approaches with the Earth; (d) critical argument of the 3 : 8 mean motion resonance with the Earth (given as number of revolutions).

major axis changes essentially only because of close approaches. However, there are cases which cannot be handled in such a simple way. As an example, the orbit in our sample consistent

with the initial conditions of 1954 XA has 16923 approaches to Venus and 14266 to Earth within 0.1 AU; the plot of a(t) (Fig. 8), with a much lower number of data points displayed than the number of

320

MILANI,

CARPINO,

FIG.

AND

MARZARI

h-Conrinued

has the typical fractal appear‘jumps,” ante of the outcome of a random walk. Nevertheless we can show that the semimajor axis of 1954 XA changes only because of the close approaches, by using the list of orbital changes derived from the close appreach data base. Given a chronologically

ordered list of close approaches experienced by a given asteroid, we compute for close approach number k the semimajor axis a; at a time tk “just before” and a: at a time th+“just after”; we assume the entire change to occur in an impulsive approximation at tk = (tk+ + &)/2 and we define the

close

NC t1 0 -100205.338S 0 -49140.7070 1 -72512.6443 1 -60175.8506 0 -49950.0318 0 -47997.0009

not

not

PA

at

at

at

-49020.0043 -47967.1450

-61953.1996

to -100205.33SS -87227.9261 -71510.9409

not

to -92911.4121 -79452.4290 -59301.84SO -43805.6126 -41943.6741

1900

node

nod,

nodes

-66970.5011 -49824.3750 -45660.1404

-70612.6674

-99152.1093 -S5534.7910

t2

crossings

t2 -926Il7.0015 -76940.7025 -56796.4351 -43601.7246 -41195.1000

crossings

=

0

=

I

2337.6665

Dt 1053.2295 3605.9160 1099.7569 1905.26011 125.6560

0

4724.7875 5154.1176 1787.5307 165.6286

Dt 310.9024

0

@=/Y 0.00000 0 .ooooo 0.00000 0.00000 0.00000 0.00000

N-/Y 0.00000 0.00000 0.00000 0.00000 0.00000

3.8473 8.2672 4.9065 4.9050 0.4044 4.1653

NclWet

15.4049 14.1006 3.6863 1.4322

NclWat 0.069

CROSSING(WETHERILL'S PARAMETER = 0.02 AU)

NS 0 0 0 0 0 0

AS 0 0 0 0 0

FORMULA)

2.000 2.000 2.000 2.000 2.000 2.000

RalEr

2.000 2.000 2.000 2.000 2.000

RelEr

0.0003 0.0074 0.0014 0.6674 0.0005

Pois 0.0213

0.0250 0.2388

0.0000

0.0000

Pois 0.4119 *a*

l

**

**

**

*a

**

**

9"

l

l

l

l

l

Note. Example of output of the crossing classification program, for the orbit of 1980 PA and the Earth. Each generalized crossing is classified according to the definitions explained in the text and Fig. 5 (TYPE). For each generalized crossing the table gives the number of times at which the distance between the two orbits passes through zero (NC), the starting time (t,) and the ending time (r2), the time (r,,) of minimum orbital distance, the duration (Dr = t2 - t,), the number of observed close approaches per unit time (NCL/y), the number of close approaches predicted by Wetherill’s formula (NclWet) and the total number of close approaches recorded in the SPACEGUARD integration between t, and t2 (NS), the relative error of the prediction cW = (NW - N&l(Nw + No,,)/2 (RelErr) and the Poisson probability that such an error is occurring by chance (Pois). The asterisks mark the crossings for which the discrepancy is big and statistically significant; however, due to the approximations introduced, the comparison indicates in a reliable way a protection mechanism only for deep crossings (DC).

AC+ AC+

DCDC+

AC+

TYPE UC?+

DESCENDING N. of close

NODE approaches

-61952.5527

-45509.2553 -41960.7286

2 2

0 0

sc+ sc+

AC+ AC+

AC+

t1 -92997.9039 -61673.4900

NODE approaches

16 171 approaches

UC 0

A

TYPE

ASCENDING 1. of close

of

EARTH

A.

VIII

FOREACH IMPACT

TABLE

EXPECTEDNUMBEROFCLOSEAPPROACHES FORORBIT 1980 PA(MAXIMUM

322

MILANI,

CARPINO,

AND

MARZARI

FIG. 7. (a) Semimajor axis of (186.5) Cerberus, an easily identifiable example of a Geographos type orbit; the solid line (above) is the filtered value obtained from the numerical integration, the dashed line (below) is reconstructed by adding up the changes produced by close approaches only (Eq. (12)): this plot is shifted 0.05 AU below the value-for the arbitrary constant ui of (12)-which gives the best agreement with filtered data, in order to avoid superposition of the two lines. The evolution of the semimajor axis is essentially produced by close approaches. On the contrary, the changes in eccentricity (plot b) are mainly produced by long-range secular perturbations.

step function u’(t) = c
(12)

representing the change due to the close approaches only; that is, u’(t) would be essentially equal to u(t) if the semimajor axis were exactly constant between close approaches. The arbitrary constant a6 can be used to obtain for u”(t) the same mean value

as for u(t). Then we compute the correlation coefficient between a(t) and a”(t) over a time span [T, T + AT]: ~(7.3 AT) = 7 IAl J7 (u(t) - ii)(u”(r) - U”)df l“,l (a(t) - ii)? fir j;+yu’(f)

FIG. 8. Semimajor axis of asteroid 1954 XA (SPACEGUAKD

filtered data).

- 2)? fit (13)

PROJECT

where li,z are the averages over the same interval; the integrals are of course approximated by the sum over the sampled time series available. When the correlation (13) is computed over the entire data time span, the behavior of two apparently different cases such as 1865 Cerberus and 1954 XA is shown to be the same: both have a correlation p 2 0.99. By performing such a test on all the orbits for which we have reliable orbital changes data (see Appendix l), the averages for each class are quite distinctive (see Table IX). More information can be obtained by separating the contribution of the different planets; that is, given a subset K of the list of close approaches experienced by a given asteroid (e.g., the ones with Venus only), we can form the partial step function & by restricting the sum (12) to the k in K. However, the partial correlation pK is not the correlation of a(t) with ak(t), because the correlation coefficient is not a linear functional; we define as “partial correlation” pK the partial contribution of the subset K to the integral (13) in such a way that if the total list of close approaches is the disjoint union of two subsets K, H, then

P =

Class (1)

Geographos TOM Aliida KOZai ERM

OF CHANGES Pv

PI3

PM

(4

(3)

(4)

(5)

0.91 0.11 0.08 0.93 0.65

0.17 0.06 -0.02 0.00 0.00

0.73 0.05 0.07 0.11 0.27

0.01 0.00 0.02 0.82 0.38

+

PH:

-ii) -a’,)dt c:,,,(a(r)

AT)

=

T+AT(u(t) - ii)2 dr j-rAT(uc(t) - aC)?dt (14) where the r.m.s. in the denominator are those of u(t) and of the total step function u’(t) rather than that of the partial ai( In this way we can display in Table IX, for each class, the average value of the cot-relation for all close approaches (column 2) and the partial correlations for the close approaches to each planet (columns 3-6), with column 2 equal to the sum of columns 3 to 6. We also give the average value of the total variation and the relative contribution of each planet (columns 7-11); the total variation is defined as vc = c ia: - ai1

(15)

and is given as a measure of the total size of the impulsive changes. Table IX is a good confirmation of the soundness of our taxonomy; that is, the definitions of the dynamical classes we

IN SEMIMAJOR

PM

PK

pK(T,

TABLE CORRELATION

323

SPACEGUARD

1X AXIS

WITH

VE

V,‘IV”

(6)

(7)

(8)

0.00 0.00 0.01 0.00 0.00

3.303 0.896 1.217 0.073 0.325

0.23 0.20 0.23 0.00 0.00

PJ

CLOSE vgv=

APPROACHES VA/V=

V,“IV’

(9)

(10)

(11)

0.74 0.74 0.59 0.54 0.60

0.01 0.05 0.07 0.45 0.39

0.00 0.00 0.09 0.00 0.00

Notr. Average values of the correlations between observed variation in the semimajor variation produced by close approaches only, for orbits belonging to different dynamical averaging the orbits have been weighted proportionally to the time of residence within each class the table shows the total correlation (pm,, column 2) defined by Eq. (13), the ptO, of each planet from Venus to Jupiter (columns 3-6, Eq. (14)), the total variation of axis produced by close approaches (Vc, column 6) according to Eq. (19, and the relative to Vc of each planet (columns 8 to 1 I).

axis and the classes; in the the class. For contribution to the semimajor contribution

324

MILANI.

CARPINO, AND MARZARI

have given in our previous paper gave quantifiable prediction which are, on average, confirmed by this test. The Geographos members have a large p, mostly due to the Earth. A separate computation performed only on the Geographos members which are deep crossers of Venus shows that pv = 0.59, pa = 0.35; this can be understood because close approaches to Venus are more frequent when they are geometrically possible, because synodic periods with Venus are shorter. Toros and Alindas have a low p; that is, the changes in a are mostly due to the resonances. Kozais have again a large p, but this is due to Mars, The situation is somewhat more complicated for the Eros class; some have large p, in many cases due only to Mars; for others, shallow approaches to the Earth are more relevant. Moreover, resonances with the Earth can occur even when the orbit is not Earth-crossing, and when this is the case the changes in a due to the resonance are much larger than the changes due to close approaches to Mars, resulting in a very low correlation. For the reasons explained in Appendix 1, we could not apply this method to most of the Oljato type orbits, and therefore Table IX does not have an entry for this class; however, in the cases in which we can obtain a correlation, it is usually high (p > 0.8 for three orbits out of four). However, Table IX cannot be used to confirm our classification, that is the attribution of a single orbit (or part thereof) to a specific class. To this purpose we need to perform at least two tests. Given the classification as it is, we can compute p for each orbit, and if the classification reports a transition, the computation must be performed separately for the time spans of residence in each class. Then we can check whether, e.g., all the supposed Geographos members have a high p; and of course we find that there are both marginal cases and actual mistakes. As an example, the orbit of 1980 PA for the time span in which it had been classified Geographos has p = -0.06; since

we have already been caught red-handed by Wetherill method, we better acknowledge that 1980 PA is never a Geographos. Two more dubious classifications are: 1917 Cuyo, with p = 0.33 for the time span [-- 100200, +73700] in which it had been classified Geographos (this appears to be a really borderline case, which behaves in a way intermediate between the Geographos, the Kozai, and the Eros class); and 1979 VA, with p = 0.38 (which experiences short spells of resonance with Jupiter). Apart these three bad cases, all the remaining ones have p 2 0.76, and most of the Geographos members we have been able to analyze in this way (28 out of 34) have p z- 0.9. In the same way, all the Alindas we have been able to test have p 5 0.25 with the exception of 1863 Antinous (p = 0.69); all the Toros have p I 0.38; all the Kozais p 2 0.88 but for 1863 Antinous (p = 0.55). A more refined test can be performed by computing the correlation p(T, T + AT) for a fixed AT as a function of T. This is a very powerful method, but unfortunately we have not been able to perform it in a fully automated way, because the choice of the time span AT must be done by exploiting some knowledge of the dynamics of the specific orbit under investigation. Of course if there are no close approaches the correlation is zero, thus we need to make sure that AT is longer than the longest interval between two successive node crossings. On the other hand if AT is too long, a comparatively short-lived transition (such as many Toro states) might not show up. As an example, Fig. 9a shows the correlation as a function of time for 1685 Toro; here AT = 20,000 years, and the transition to the Toro class is perfectly obvious. A study of the partial correlations with Earth and Venus shows more interesting details; when close approaches with the Earth are absent (because of the 5 : 8 resonance) the partial correlation with Venus rises and becomes dominant, until the 5 : 13 resonance with Venus generates another protection; this very strange double resonance state lasts

PROJECT

-100000

-50000

325

SPACEGUARD

0 SPACEGUARD

50000 integration

100000

t LLLLLLLLLLLLLLI,,li, 100000 -50000 0 Tmne (y) SPACEGUARD

50000 mtegratlon

100000

Time

(y)

-

FIG. 9. Correlation of semimajor axis changes with the changes produced by close approaches as a function of time for asteroids (1685) Toro (a) and (1865) Cerberus (b). The total correlation is shown (solid line) along with the contribution from each planet (dotted line: Venus; short-dashed line: Earth; long-dashed line: Mars).

only for a single almost crossing with Venus. A further peculiarity is that this rare state is occurring to the real 1685 Toro just now (see our previous paper, Section 5.2; Janiczek et al. 1972, Williams and Wetherill 1973). Figure 9b shows the pure Geographos type orbit 1865 Cerberus, with AT = 40,000 years; since the period of o is quite long, a similar plot with AT = 20,000 years would show sharp decreases of p which only correspond to intervals between node crossings. The alternation between

Earth and Venus in controlling the orbital evolution of 1865 Cerberus also corresponds to the alternation of the node crossings with the two planets. On the basis of our experience with such plots, we can say that every occurrence of values of p 5 0.6 in a supposed Geographos needs to be investigated; if there are node crossing in that interval, in almost all the cases a short-lived resonance is the cause of the drop in the correlation. So far we have used the correlation

326

MILANI,

CARPINO,

method to confirm and/or correct our own classification, but this is not the only purpose of this Section. We are interested in testing assumption (3) (see Section 1) of an abstract statistical theory of planet-crossing orbits, that is, the possibility of describing the orbital evolution as a random walk controlled by close approaches, in such a way that a Monte Carlo simulation is possible. At this point, we can state that this is indeed possible and is a good approximation (p = 0.9) for the semimajor axis of a Geographos and of a Kozai. The low values for the Alinda and Toro classes, and the mixed results for the Eros class, indicate that a test for resonances with both Jupiter and the Earth should be performed before using such an approximation. Another interesting question is the following: if the evolution of a is determined by the close approaches, which are the most important ones in controlling the orbital evolution, the few very deep or the many more shallow ones? That is, if the changes in a have to be modeled as a random walk, is it better to use a short time step-corresponding to the waiting time for a close approach within a comparatively large impact parameter D-or rather take into account only the possible occurrence of low D approaches? Of course, the possibility of using a long stepsize would be a bonus for Monte Carlo simulations. To investigate this problem, we have again used the partial correlations as defined in Eq. (14), but this time we have partitioned the list of all close approaches not only by planet but also by minimum approach distance 6 = &rN. The interval [O, 0. l] for 6 has been partitioned into 10 equal intervals [0, 0.011, [O.Ol, 0.021, etc., and we compute the partial correlation with a(t) of the step function constructed only with the close approaches in a given subinterval. To be more accurate we could have used the impact parameter instead of the minimum distance, but this does not actually matter because the two are significantly different only below 0.01 AU. Of course this compu-

AND

MARZARI

tation can be performed only as an average over many orbits, such as all the Geographos members; otherwise, the numbers of close approaches would be too small for statistical significance. Table X shows the results averaged over the Geographos, Kozai, and Eros orbits for which we have good orbital change data (in this table, 1980 PA is still counted among the Geographos members although we now know this to be wrong). For each minimum distance bin (column 1) we give the partial correlation, the total variation and the number of close approaches in the bin for Venus (columns 2-4), for Earth (columns 5-7) and for Mars (column 8-10). The result can be described in a qualitative way as follows. If the orbit belongs to a class in which close approaches down to physical collisions are possible, then the correlation is significant only for the smaller distances, while the total variation is almost constant up to quite large distances. If very close approaches are not possible, such as in the Eros and Kozai cases with the Earth, then the most important ones in controlling the orbital evolution are the closest possible, even when they are not very close. A simple model to help understand this experimental result can be devised by assuming that the number of close approaches N is proportional to 6*; thus, the contribution of those between 6 and 6 + d6 is proportional to 6dS. Since in the hyperbolic two-body encounter approximation the change in semimajor axis is proportional to l/6, the contribution to the total variation (1.5) is roughly independent from 6. This can be seen in the total variation data, which are constant up to distances much larger than the radius of the sphere of influence-that is, up to distances at which the two-body approximation used in the above argument breaks down completely. To estimate the contribution to the correlation, that is, to the overall evolution of the variable a, we have to use the random walk formula by which the r.m.s. of the changes in a is proportional to da%%, with da the

PROJECT SPACEGUARD TABLE

X

CORRELATIONOF SEMIMAJORAXISCHANGESASA APPROACH DISTANCE 6k bin (AU)

(1) Class

Pk

VQ.Illls v;

Nk

Pk

(2)

(3)

(4)

(5)

FUNCTIONOF Mars

Earth v;

Nk

Pk

v;

Nk

(6)

(7)

(8)

(9)

(10)

0.149 45 0.48 0.152 157 0.21 0.148 263 0.08 0.143 374 0.04 0.143 500 -0.05 0.141 645 -0.01 0.130 808 -0.03 0.116 902 0.01 0.094 1006 0.01 0.057 1052 0.00

0.379 0.370 0.356 0.317 0.269 0.236 0.200 0.164 0.119 0.059

59 165 267 342 391 458 525 598 652 666

0.01 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.009 0.009 0.009 0.009 0.008 0.008 0.007 0.006 0.005 0.002

10 36 60 85 109 136 154 179 193 203

0 0.56 0 0.27 0 0.01 0 -0.05 3 -0.06 15 -0.01 22 0.06 34 0.00 35 0.04 52 0.02

0.006 0.004 0.004 0.003 0.004 0.004 0.003 0.002 0.002 0.001

4 14 21 26 36 46 46 54 61 68

0.016 0.017 0.015 0.016 0.014 0.014 0.012 0.011 0.009 0.005

11 34 53 76 91 121 137 158 190 193

= GEOGRAPHOS

030-0.01

0.01-0.02 0.024.03 0.03-0.04 0.04-0.05 0.05-0.06 0.06-0.07 0.07-0.08 0.08409 0.09410 Class

327

0.17 0.00 -0.04 0.06 0.02 0.00 0.01 0.01 0.01 0.00

= KOZAI

0.0Oa.01 0.01-0.02 0.02-0.03 0.03-0.04 0.04-0.05 0.05-0.06 0.06-0.07 0.07-0.08 0.08409 0.09-0.10

-

-

-

0.00 0.00 0.00 0.00 -0.05 -0.17 0.12 0.02 0.09 0.09

0.000 o.l.Mo 0.000 0.000 0.001 0.006 0.008 0.011 0.009 0.006

Class=EROS 0.00-0.01 0.01-0.02 0.02-0.03 0.03-0.04 0.04-0.05 0.05-0.06 0.06-0.07 0.07-0.08 0.08409 0.09-0.10

-

-

-

0.00 0.00 0.00 0.01 0.01 0.05 0.06 0.09 -0.01 0.12

0.000 0 0.000 0 0.000 1 0.008 5 0.014 7 0.019 20 0.039 45 0.053 78 0.058 112 0.040 155

0.16 0.20 -0.01 0.06 -0.02 -0.01 0.03 -0.01 -0.03 0.01

Note. Average values (weighted as in Table IX) of the correlations between observed variation in the semimajor axis and the variation produced by close approaches only, for the orbits belonging to the three dynamical classes with higher values of correlation and showing the contribution of close approaches at different minimum distance. For each distance bin (column I) and for each planet the table shows the contribution to the correlation (& and to the total variation (Vi) of the Nk close approaches included in the bin.

average value of the jump in a. As a result, the correlation of a bin of close approaches around distance 6 should decrease as l/V%. Of course this is a rough approximation,

since as we have shown in Section 3 and Appendix 2, the number of close approaches within 6 fails to be proportional to a2 in many cases.

328

MILANI.

CARPINO,

This finding gives an important indication on the procedure to be used for a Monte Carlo experiment: for an orbit which can have very close approaches, such as a Geographos, it is a good approximation to wait for a long time for the next likely close approach within a very small impact parameter D; for an orbit which can have close approaches only with D > Do, such as the Eros and Kozai cases, only the close approaches with D just above Do really matter. All this applies to the semimajor axis. We do not know how to apply the same technique we have used for a to describe the evolution of the other orbital elements. Ideally we would like to describe the long range perturbations by some sort of “proper elements” u*, e”, Z* which, unlike the osculating orbital elements u, e, I, are exactly constant unless a close approach occurs. If such proper elements are computed starting from each output of an orbit computation which includes the effect of the close approaches, they will not be constant, and their change will represent the effect of the close approaches only. Then the correlation coefficient technique could be applied as above. Unfortunately we do not have an algorithm to compute such proper elements. Outside the main mean motion resonances, u” = N is a good enough approximation; e* = e would be a very bad approximation, and therefore one has to resort to a secular theory. The traditional tool of linear secular perturbation theory fails for high eccentricity and/or inclination orbits. Semianalytical secular perturbation theories based upon the technique of averaging do provide proper elements even for high r, I orbits (Kozai 1962, Wilhams 1969, 1979) and can be used inside mean motion resonances (Schubart 1982, Schubart and Bien 1986); however, none of the methods used so far can be applied to crossing orbits. Thus, we have not attempted to compute a correlation coefficient for e.g., the eccentricity; as it is apparent from any example, such as Fig. 7b, the eccentricity undergoes significant

AND

MARZARI

changes even between close approaches, and its overall evolution is determined more by secular perturbations, including secular resonances, than by the close approaches. 6. CONCLUSIONS

The amount of data which can be extracted from a massive numerical integration such as SPACEGUARD is enormous, and we think that the main result of this paper is in having shown that it is possible to devise algorithms to retrieve useful information on complex problems such as close approaches and collision frequency. We have also presented some specific results, which we would like to summarize and comment upon. (I) We have shown that, given the time series of the orbital elements, the number of close approaches is predictable in a statistical sense. This is true with a very good accuracy for comparatively shallow close approaches (with impact parameter D between 0. I and 0.02 AU), and provided the expected number is computed in a way which avoids the singularities of zero inclination and of tangent orbits. We have derived and tested a nonsingular formula which allows one to avoid such a singularity in all cases; the agreement is within 10%. This statistical agreement results from the averaging out of significant discrepancies for individual orbits; many of these discrepancies can be understood in terms of protection mechanisms, such as the Toro and the Kozai states. Vice versa, if the number of close approaches is known from a complete numerical integration, the discrepancies with a statistical theory such as the one by Wetherill provide a reliable indication of the existence of some protection mechanism. (2) The smaller the value of the impact parameter the more difficult and the less reliable is the statistical prediction of the number of close approaches. This results from the existence of near-singular cases, in which a very high frequency of very deep

PROJECT

SPACEGUARD

close approaches occurs for a very limited range of orbital elements, and possibly for a short time. In this respect, numerical integrations and analytical probability computations are similar in that a reliable estimate of the number of very close approaches (down to physical collisions) is possible only by averaging over a large number of orbits and over a long time span. This problem is made more severe by the fact that a small fraction of the orbits contributes a large fraction of the total number of close approaches; e.g., for close approaches to Earth with D 5 0.02 AU, eight orbits contribute more than half of the close approaches found in SPACEGUARD, and of these, two (1982 DB and 1954 XA) contribute almost 30%. Therefore, results presented as an average over 89 orbits are in fact determined mostly by an average over a handful. Moreover, the close approaches are not uniformly distributed in time, but appear in bursts; this occurs both in the numerical data and in the analytical formulas, in the former case because of the sporadic nature of the node crossing events, in the latter because of the singularities. Thus an orbit needs to be computed for a long enough time span, and the orbital evolution needs to be followed with good accuracy and time resolution. (3) The dynamical classes defined in our previous paper are a good tool to improve the predictions of the frequency of close approaches. However, as was pointed out by Wetherill (1979), a classification should not be used in the context of a box model. Our classes define a partition of the phase space into subsets where different rules have to be used to compute the frequency of close approaches (e.g., use Kessler formula for a Geographos, but make sure to use the nonsingular formula for an Eros; take into account e-w coupling for a Kozai, etc.). This does not at all mean that all the orbits in the same class evolve in the same “average” way; e.g., the expected number of collisions with Earth, for different orbits all in the Geographos class, can change by a

329

factor more than 100, as a result of the very different values of a, e, I which can occur within the class. (4) Protection mechanisms based upon the mean motion resonances can be important-for a given orbit and over a given span of time-in decreasing the rate of close approaches, e.g., to the Earth, and therefore in freezing the orbital evolution. However, for a large enough population and for a long enough span of time, the change in the average frequencies is quite small, and barely detectable, given the other uncertainties and sources of error. In other words, in statistical studies and Monte Carlo simulations, the distinction between Geographos and Toro could either be ignored or accounted for by means of a “Toro correction” of a few percent. This approximation might fail to predict the correct frequency of very rare events, such as physical collisions, because we cannot exclude that resonance might occur exactly at the critical times when the probability of collision is very high; however, to predict an average number of moderately deep close approaches an explicit model of the resonances with the Earth might not be necessary. (5) The accuracy of the orbital model used in generating the orbital element time series is very critical in the computation of the collision probabilities and of the frequency of close approaches at small D, it is not critical at all in the computation of the number of shallow close approaches. Surprisingly, for 0.1 2 D 2 0.02 AU, even a model with precessing keplerian ellipseswith constant a, e, Z-is good enough. (6) The approximation by which the changes in the semimajor axes of an Earthcrossing orbit arise only from the close approaches is a good one, provided close approaches do occur, i.e., for a Geographos. The same approximation is totally wrong whenever a mean motion resonance is encountered; this occurs in the Alinda case of a resonance with Jupiter, which changes in a significant way a independently from the

330

MILANI.

CARPINO.

close approaches, and in the Toro case, where the resonances with the Earth are less effective in changing a, but can prevent deep close approaches. For a Mars-crosser the changes in a can be driven either by close approaches to Mars, or by shallow approaches to Earth, or by resonances. However, if the goal is to study the long term evolution of the orbit, it is legitimate to argue that the mean motion resonances with Jupiter are more important, since the resonances with the Earth do not last very long and do not result in dramatic changes of the orbital elements. Having listed what we think we have achieved, we would like to conclude by listing the problems which remain to be solved. (1) We do not know how to compute proper elements for planet-crossing asteroids; therefore, we do not really know how to describe the results of the numerical integrations, apart from showing plots of e, I, etc., and we do not know how to assess the relevance of the close approaches in modifying this evolution. From a Monte Carlo simulation performed by Wetherill (1988) it appears that the synergic effect of close approaches to Mars and of a secular resonance is much larger than the linear superposition of the two effects. In the same way, we suspect that the orbital evolution of the orbits of the Eros and Kozai classes could be much more interesting than one would expect from the effects of the close approaches to Mars alone. It might be possible to use averaging methods to provide proper elements, since to ignore the perturbations of Mars altogether is a reasonably good approximation (Williams and Faulkner 1981) for time spans of the order of a few tens of thousands years. It might in the future become possible to use higher order proper elements (such as KneieviC and Milani 1989, Milani and Kneievid 1990) suitably modified to give at least a first approximation to the secular behavior in the near-Earth region. However, secular perturbations for Earth-crossing orbits require new analytical tools.

AND

MARZARI

(2) We need a better understanding of the resonance phenomena both with Jupiter and with Earth. In particular the problem of discriminating between resonances with Jupiter and resonances with the Earth appears to be more difficult than we had expected. This we have learned from our mistakes in the cases of 1982 DB and 1982 XB. This could be investigated by coupling the output of SPACEGUARD with automated tools of resonance analysis. We are working on this. (3) We do not understand at all the orbits we have called Oljato. This, results also from the fact that they are the ones which are not suitable to be investigated by the methods used within SPACEGUARD, such as parallel computing. On the other hand the Oljato are another example of synergic effect of close approaches to the inner planets and long range perturbation by Jupiter, which result in very large changes in the orbital elements; they are likely to be very interesting to understand the decoupling of the short periodic comets from Jupiter. (4) We have developed very powerful tools to analyze the behavior of an orbit computed numerically; at this stage, the problem appears how to integrate the informations obtained by different tools in a way less cumbersome than the one we have used here, with hundreds of cross-checks between plots and tables (see the 1980 PA story in Sections 4 and 5 as an example). Can these methods be integrated in something of the nature of an expert system? Can this allow an objective and essentially automated qualitative classification of the dynamical behavior of planet-crossing orbits? APPENDIX

1: POSTPROCESSING SPACEGUARD DATA

OF THE

The output data from the SPACEGUARD numerical integration have been postprocessed to generate the data files (l), (2), and (3) described in Section 2. No special problem was found in the generation of (I), the postfiltered element time series; the

331

PROJECT SPACEGUARD

filters used are described in our previous paper, Section 3 and Tables I-II. The theory behind the use of these filters is in Carpino el al. (1987). The close approach output of SPACEGUARD was generated by checking the distance between each asteroid and each planet at each integration step; if the distance was less than 0.1 AU (1 AU for Jupiter) the position of both the planet and the asteroid was stored every 2 integration steps for the inner planets and every 10 steps for Jupiter. In the postprocessing phase, for each close approach the minimum distance among the recorded states was found; the minimum distance state (2) (including relative position and velocity) was obtained by polynomial interpolation of degree 10, with nodes at the 11 data points surrounding the minimum. If less data points were available, the degree of the polynomial was decreased accordingly. All this is well provided there is a single minimum. However, there are rare occurrences of close approaches with a much more complex structure (see Fig. 10). It is not obvious whether we should count them as a single close approach or as many, not to speak of the fact that the interpolation algorithm might fail to correctly identify all the minima. Out of the 321,560 close approaches we have used for the statistical treatment in this paper, 42 turn out to be multiple minima occurring in 19 distinct close approaches; most of these (32 minima in 14 close approaches) are with Jupiter. Moreover, for 8 close approaches our interpolation software either failed or encountered a numerical instability. From the point of view of understanding very complex and extremely rare events, such as temporary satellite captures, each one of these wild cases would be interesting (although such a problem has already been investigated; Carusi and Valsecchi 1981). For the purpose of our statistical investigation, such a small number of cases does not matter. Therefore, we simply ignored the 8 cases of complete failure and used each minimum of the multiple minima cases as a

1983

VA -

Juprter

FIG. IO. Trajectory of the orbit of 1983 VA during a close approach with Jupiter lasting for 13.8 year, starting from f = 81745.1 year; the relative position of the asteroid with respect to Jupiter (in astronomical units) is projected in the ( Y, 2) plane of the inertial (ecliptic) frame. The event is actually a temporary satellite capture; the distance of the asteroid from the planet has four separate local minima, with dMIN = 0.447, 0.322, 0.262. and 0.225 AU.

separate close approach. To give an idea of how little our conclusions depend upon such choices, we can use as an example Table IV: since a close approach of 1982 DB to the Earth (with minimum distance hard to compute, but below 0.003 AU) is between the total failure cases, all the entries in column 4 for the Earth and the Geographos state should be increased by 1. The preparation of the list of the immediate orbital changes due to the close approaches (3) was even more plagued by exceptional cases and numerical problems, essentially because correlation analysis was an afterthought and the output of SPACEGUARD was not designed to provide this information. In particular, the asteroid velocity was not stored in the close approach records, because to compute the asteroid velocity for only one asteroid implies the loss of both vectorization and parallelism. Thus a significantly reduced computational efficiency on a parallel supercomputer such as the one we were using (this is the “load unbalance” problem dis-

332

MILANI.

CARPINO.

cussed in Milani et al. 1989a, Section 6). Our idea to recover the missing data was to use the first five and the last five records of each close approach to interpolate with a polynomial of degree 4 the velocity of the asteroid at the initial and final record; in this way the orbital elements at the beginning and at the end of each close approach could be computed. This procedure did not work well for all the orbits; e.g., for the orbits computed with a stepsize of 2 days, the close approach data were sampled only every 4 days, and with fast approaches to the inner planets lasting only 10 days or even less the amount of data stored was not enough. For the orbits computed with a stepsize of 3 day the number of data points was enough, but for some very high velocity and very deep approaches a form of numerical instability arises. In these rare cases, an oscillation of the numerically integrated solution around the real orbit builds up and peaks a few steps after the closest approach, to slowly damp down after the close approach is over. This kind of instability does not degrade the minimum distance state, but can result in a very large error in the computation of the orbital elements “just after” the close approach; since the oscillation damps down, the elements “after some time” are much better. For some orbits we succeded in detecting the unstable cases by comparing interpolations of different degrees, then by replacing the corrupted elements “just after” with the elements “after some time” extracted from the orbital elements time series. However such a procedure is legitimate only for the cases in which a very large jump in a can easily be discriminated from the background; e.g., the plot of Fig. 7a contains one such patch. Since the purpose of the orbital changes file was to test our classification, and the random walk model for the changes in a, we did not really need a very large sample of orbits; the main point was quality control on the data we were going to use rather than quantity. Thus we have used the

AND

MARZARI

above mentioned test for instability, and discarded 10 more orbits from our sample as unsuitable for a reliable correlation computation because of the instability problem; moreover, we also discarded the Jupitercrossing orbits because the effects of the close approaches to Jupiter have to be investigated in a different setting. Therefore, the data of Section 5 actually refer to only 47 orbits. Unfortunately, the discarded orbits included most of the Oljatos, including the only three cases of orbits classified as Oljato for the entire integration time span. This did not occur by chance. The Oljatos have by definition a very large eccentricity, a low inclination, and are Venus crossing; thus, they have many very fast approaches to Venus. Actually the Oljatos have all the characteristics which make an orbit unsuitable for an investigation using vector and/ or parallel computers (very high eccentricity, very fast close approaches, very large changes in semimajor axis), and their orbits have not been computed reliably enough. Even the assumption that the overall orbital evolution is not dependent upon the numerical errors introduced locally is not necessarily true if there are many numerical unstable close approaches (see Froeschlt and Rickman 1981). The only credit we can claim for SPACEGUARD is that we have pointed out the existence of orbits of this type, and therefore the need to investigate them by much more sophisticated methods, which have to include the use of a numerical integrator with either numerical or analytical stepsize control (see Hahn and Lagerkvist 1988). APPENDIX

2: NONSINGULAR FORMULA

KESSLER’S

The possibility of removing the singularities from Kessler’s formula (l)--(2) has been discussed by Steel and Baggaley (1985); they show that near both singularities (zero mutual inclination and tangent orbits) a suitable average should replace the singular spatial density (Eq. (1)) in formula (2).

333

PROJECT SPACEGUARD

However, their formula does not account for the double singular case in which both ZM = 0 and q = a, and it turns out that there are two dangerous neighbors-1982 DB and 1982 XB-whose orbits are sometimes very close to this double singularity. The integral (8) is always finite, since

imply:

1

(r - q)(Q - r) = u*e*

.

(20)

Thus the appropriate

substitution

is

r = a(1 - e cos u). 27~ I, S(r, PM; A = {(r - ~2’)~+ r* sin2(/3) 5 D2}

(16)

(21)

That is, the angle u is the eccentric unomuly, and the integral (17) becomes

is the probability per unit time of finding the 2 arcsin g(u)du, (22) asteroid in the solid torus defined by (r, /3) E A, A E [0,27r]. This can be understood by remarking that all the singularities of S(r, where the extremes of integration are chop) are of order t; thus, (16) converges sen to ensure that the integral of S is not as an improper integral and can be trans- computed outside the region where the formed into an ordinary integral by a suit- square root in the denominator of S is posiable change of coordinates. To find such a tive: nonsingular computational procedure we cos uj = (U - rj)lue for 4 < rj < Q expand the integral as a double integral in r and z = r sin p: = 0 for = 4 Uj

W=

IAS(r.P)dA

=-

1 r? 2r’U I rt -\/(r

X

=

Uj

Yj

(23

dr - q)(Q

for rj = Q

= 77

- r) x

g(u) = min

Z(T) I

m’(” r sin ZM j/k’

(17)

where the extremes of integration are defined by the geometry of A (see Fig. 2): rl = max(q, a’ - D); t-2 = min(Q, a’ + D)

(18)

z(r) = min(r sin Z,, V/o2 - (r - a’)*). (19)

c 1,

-\/D? - (r - u’)~ 1’ r sin ZM

(24)

The integral (22) is defined in every case (if we set by definition W = 0 for either a’ + DQ,andg(u)= lforsin ZM = 0) and is nonsingular. An efficient computational procedure for (22) should take into account that if the disk A touches the line z = r sin ZM, there is a subinterval [uj, u4] of [ur, u2] on which g(u) = 1. Thus, I 4” 2 arcsin g(u)du = T(Q - u3). (25)

The integral in dz can be simplified by the standard substitution x = zl(r sin ZM) and That is, the low ZM approximation (see Steel gives an arccosine function. The integral in and Baggaley 1985, formula Eq. (32)) can dr has a singularity of the same type; thus, be used in that subinterval; outside [uj, u4] it can be regularized by the substitution of the more complicated formula (22) can be the square root in the denominator with a integrated numerically. This is the algosuitable sine function multiplied by a con- rithm we have used. stant. The appropriate substitution, and The probability per unit time of a close also its geometrical meaning, can be easily approach through the cross section B found from the definitions of 4, Q, which should be the integral of the spatial density

334

MILANI.

CARPINO.

over B multiplied by the relative velocity, where the disk B is orthogonal to the direction of V. Formula (8) assumes that it is possible to replace it by the integral over A, a disk of the same size as B but oriented perpendicular to the velocity of the circular orbit of the planet (in a plane A = const). This could lead to the same kind of errors of the singular formula, e.g., in case q is larger than the maximum of r on B but less than a’ + D. However, in these near singular cases the angle between the velocity of the asteroid and the velocity of the planet is very small. The same occurs near the ZM = 0 singularity: the disk A can sample a different area in the r, z plane than the disk B, but the two areas are almost the same when ZM is almost zero. Another approximation is involved in formula (8), namely the size of the relative velocity V is taken as a constant, thus pulled out of the integral and computed at the point Y = a’, z = 0. A more accurate formula could also account for the change of V(r, z). However, the relative change is not very large in most cases; the only exception is the case with u = u’, where the relative velocity can change significantly in size, but we do not think that statistical theories can be applied to nearTrojan cases (see the example of 2063 Bacthus, an Earth Trojan, which is a Toro, in Table VII of Paper I). There is no doubt that a more accurate estimate of close approach frequencies could be computed by more complicated formulas, and for that matter it would be possible to take into account the eccentricity e’ of the planet with a minor modification of the algorithm we have used. However when the agreement between statistical theory and close approach data is within 10 to 20%, as is the case for moderately close approaches (D 5 0.02 AU), it is doubtful that an improvement could be possible and worthwhile. The area in which improvements of the statistical theories would be worthwhile is in the prediction of the frequency of very close approaches, down to physical colli-

AND

MARZARI

sions. Although we do not know how to do this, we would like to point out that it is not impossible in principle. The probability for unit time for a cross section zero (that is, the probability for unit cross section area at the limit for D + 0) is just given by the function S of Eq. (2) multiplied by the smooth function V(r, z). The probability over a finite time span, over which the orbital elements of the asteroid and the planet are given by the functions u(t), e(t), ZM(t), u’(t), is a time integral of k(t) = S(u’, 0) . V(u’, 0). In almost all the cases, the time derivative of the function f(t) = a(t)(l e(t)) - u’(t) will be different from zero at the times whenf(t) = 0, that is, when k(t) has a singularity. Then the singularity of k(t) is of order i and the integral exists. A similar argument can be made for the ZM = 0 case. This abstract argument shows that for the very dangerous orbits which touch the singular surface q = a’ (or Q = a’, or ZM = 0) to use a time series for the orbital elements is not enough; the (suitably smoothed) time derivatives of the orbital elements should be used. However, how to perform such a computation in a numerically stable way is not obvious to us. ACKNOWLEDGMENTS Project SPACEGUARD has been sponsored by IBM Corporation, by the Italian Space Agency (A.S.I.) and by N.A.T.O. This work was completed while one of the authors (A.M.) was at Cornell University with a visiting professorship funded on NASA grants NAGW 392 and NAGW 310. We thank E. Clementi and his co-workers at IBM Kingston, New York, for their essential assistance in the preparation of the integration program to be run on their experimental parallel supercomputer LCAP-2. We thank A. M. Nobih, G. Hahn, P. Farinella, G. Valsecchi, G. Wetherill, and C. Froeschle for useful discussions, and A. Curioni who helped us in writing the software for statistical analysis. We thank the referees for suggesting many improvements, and in particular J. Williams for proposing the observation selection test of Figs. 4c and 4d. REFERENCES ARNOLD, V. 1976. Les merhodes mathPmatiques mtcanique classique. MIR, Moscow. CAIWINO, M., A. MILANI, AND A. M. NOBILI

de la 1987.

Long-term numerical integrations and synthetic the-

335

PROJECT SPACEGUARD ories for the motion of the outer planets. Astron. Astrophys. 181, 182-194 CARUSI, A., AND G. B. VALSECCHI 1981. Temporary satellite capture of comets by Jupiter. Asrron. Astrophys.

94, 226-228

CHEBOTAREV, G. A. 1967. In Modern Computational its, Vol. 9,

Methods

in Science

Analyiic and and Mathemar-

pp. 266-274. Elsevier, New York. FROESCHL~, C., AND H. RICKMAN 1981. A Monte Carlo investigation of Jovian perturbations on shortperiodic comet orbits. Icarus 46, 400-414. JANICZEK, P. M., P. K. SEIDELMANN, AND R. L. DUNCOMBE 1972. Resonances and encounters in the inner solar system. Astron. J. 77, 764-773 HAHN, G., AND LAGERKVIST, C.-I. 1988. Orbital evolution studies of planet crossing asteroids. Celest. Mech.

43, 285-302

KESSLER, D. J., AND B. G. COUR-PALAIS 1978. Collision frequency of artificial satellites: The creation of a debris belt. J. Geophys. Res. 83, 2637-2646. KESSLER, D. J. 1981. Derivation of the collision probability between orbiting objects: The lifetimes of Jupiter’s outer moons. Icarus 48, 39-48. KNEZEVIC, Z., AND A. MILANI 1989. Asteroid proper elements from an analytical second order theory. In Asteroids ZZ (R. Binzel, T. Gehrels, and M. S. Matthews, Eds.), pp. 1073-1089. Univ. of Arizona Press, Tucson. KOZAI, J. 1962. Secular perturbations of asteroids with high inclination and eccentricity. Astron. J. 67, 591-598 MARSDEN, B. G. 1986. Catalogue of Orbits of Minor Planets, 2nd ed. Minor Planet Center, SAO, Cambridge, MA. MILANI, A., AND Z. KNEZEVIC 1990. Asteroid proper elements. Celest. Mech., in press. MILANI, A., A. M. NOBILI, AND M. CARPINO 1987. Secular variations of the semimajor axes: Theory and experiments. Astron. Asrrophys. 172, 265-279. MILANI, A., M. CARPINO, AND D. LOGAN 1989a. Project Spaceguard. Parallel Computation of PlanetCrossing Asteroid Orbits. Adv. Parallel Computing, in press. MILANI, A., M. CARPINO, G. HAHN, AND A. M. NOBILI 1989b. Dynamics of planet-crossing asteroids: Classes of orbital behaviour. Project Spaceguard. Icarus 78, 212-269.

MILANI, A., A. M. NOBILI, AND M. CARPINO 1989~. Dynamics of Pluto. Icarus 82, 200-217. GPIK, E. J. 1951. Collision probabilities with the planets and the distribution of interplanetary matter. Proc. R. Zr. Acad. 54A, 165-199. GPIK, E. J. 1976. Interplanetary Encounters: CloseRange Gravitational Interactions. Elsevier, New York. SCHUBART, J. 1982. Numerical determination of proper inclination of Hilda-type asteroids. Celesr. Mech. 28, 189-194. SCHUBART, J., AND R. BIEN 1986. On the computation of characteristic orbital elements for the trojan group of asteroids. In Asteroids, Comets, Meteors ZZ (C.-I. Lagerkvist, B. A. Lindblad, H. Lundstedt, and H. Rickman, Eds.). Reprocentralen HSC, Uppsala. SHOEMAKER, E. M., J. G. WILLIAMS, E. F. HELIN, AND R. F. WOLFE 1979. Earth-crossing asteroids: Orbital classes, collision rates with Earth and origin. In Asteroids (T. Gehrels, Ed.), pp. 253-282, Univ. of Arizona Press, Tucson. STEEL, D. I., AND W. J. BAGGALEY 1985. Collisions in the solar system--I. Impacts of Apollo-Amor-Aten asteroids upon the terrestrial planets. Mon. Not. R. Astr.

Sot.

212, 817-836

WETHERILL, G. W. 1967. Collisions in the asteroid belt. J. Geophys. Res. 72, 2429-2444. WETHERILL, G. W. 1979. Steady-state population of Apollo-Amor objects. Zcarus 37, 96-l 12. WETHERILL, G. W. 1988. Where do the Apollo objects come from? Icarus 76, 1-18. WILLIAMS, J. G. 1969. Secular Perturbations in the Solar System. Ph.D. thesis, University of California, Los Angeles. WILLIAMS, J. 1979. Proper elements and family membership of the asteroids. In Asteroids (T. Gehrels Ed.), pp. 1040-1063. Univ. of Arizona Press, Tucson. WILLIAMS, J. G., AND J. FAULKNER 1981. The positions of secular resonance surfaces. Icarus 46, 390399.

WILLIAMS, J. G., AND G. W. WETHERILL 1973. Minor planets and related objects. XIII. Long-term orbital evolution of (1685) Toro. Astron. J. 78, 510-515.