A&A 408, 1179-1196 (2003)
DOI: 10.1051/0004-6361:20031039

## Resonant returns to close approaches: Analytical theory

G. B. Valsecchi1 - A. Milani2 - G. F. Gronchi3 - S. R. Chesley4

1 - Istituto di Astrofisica Spaziale e Fisica Cosmica, via Fosso del Cavaliere 100, 00133 Roma, Italy
2 - Dipartimento di Matematica, Università di Pisa, via Buonarroti 2, 56127 Pisa, Italy
3 - Dipartimento di Matematica, Università di Pisa, via Buonarroti 2, 56127 Pisa, Italy
4 - Jet Propulsion Laboratory, 4800 Oak Grove Drive, CA-91109 Pasadena, USA

Received 19 October 2001 / Accepted 22 May 2003

Abstract
We extend Öpik's theory of close encounters of a small body (either an asteroid or a comet) by explicitly introducing the nodal distance and a time coordinate. Assuming that the heliocentric motion between consecutive close encounters is Keplerian, or given by an explicit propagator, we can compute the initial conditions for an encounter as functions of the outcomes of a previous one; in this way it is possible to obtain a completely analytical theory of resonant returns. It is found that the initial conditions of a close encounter that lead to a resonant return must lie close to easily computable circles on the b-plane of the first encounter. By further assuming that the nodal distance varies uniformly with time, due to secular perturbations, and considering the derivatives of the coordinates on the b-plane of the second encounter with respect to those on the b-plane of the first encounter, we compute in the latter the location, shape and size of collision keyholes.

Key words: minor planets, asteroids - Solar system: general

### 1 Introduction

The problem dealt with in this paper has deep roots in the history of astronomy. Since the end of the 18th century the study of the motion of a newly discovered small solar system body, comet D/Lexell, made it clear that comets can pass very close to the Earth, and that the gravitational perturbations at close planetary encounters can modify their orbits significantly (Lexell 1778).

Le Verrier (1844, 1848, 1857) re-analyzed the observational record available for this comet with the purpose of studying its motion, and found that it was not possible to determine a unique set of orbital elements at a given epoch: a region in elements space of non negligible size was compatible with the observations. Le Verrier found that the six elements could be expressed as functions of a single free parameter, which he called , whose value could vary between 1.5. He found also that the difference between his own nominal orbital solution and the one given earlier by Clausen was within the allowed range of  and, in fact, for his subsequent studies on the motion of the comet adopted Clausen's orbit. Moreover, Le Verrier determined the range of possible orbits for comet D/Lexell both before the 1767 encounter with Jupiter (lowering the perihelion distance q and making the comet observable from Earth) and after the 1779 encounter (removing the comet from visibility). This second encounter may have been very deep and as a consequence the comet might even have been ejected from the solar system on a hyperbolic orbit.

The papers by Le Verrier on this subject are remarkable, the more so since they pre-date a very recent line of research (Milani 1999; Milani et al. 2000a), for which even the stated purposes are similar. Much of the current work aims at making the recovery of poorly observed near-Earth asteroids feasible, and Le Verrier concluded his first paper on the subject stating that if, in the future, a comet would be discovered whose observations were accountable for by using his elements for comet D/Lexell, with within the prescribed range, then the new comet would just be a new apparition of D/Lexell.

For more than a century the work by Le Verrier was forgotten and the astronomers interested in recovering lost asteroids and comets used a much simpler approach. They assumed that five orbital elements were well known and only the orbital phase, that is the mean anomaly, was subject to uncertainty. Thus the Line Of Variations (LOV), consisting, in this simple approximation, of a segment along the mean anomaly axis in the orbital elements space, replaced the curved line parametrized by Le Verrier's . This approximation is actually not bad when the asteroid/comet has been lost by decades, because the uncertainty along the orbit increases with time as a consequence of the uncertainty in mean motion.

Starting from Milani (1999) a precise way to compute the LOV has been developed and adapted to modern computational techniques, and it is now a standard procedure to compute a large number of multiple solutions for the orbital elements, sampling the LOV, which is, in general, curved. The motivation for this work, apart from the recovery problem, was the need to detect possible Earth impact solutions, the so called Virtual Impactors, in the near future (Milani & Valsecchi 1999; Milani et al. 1999). The assessment of the impact risk from an asteroid observed only during a single apparition cannot be accomplished by the study of one "nominal'' orbit (the solution of the least squares fit to the observations). The nominal orbit, and its neighborhood (which can be studied by linearization techniques), can provide only information on the impacts with probabilities of the order of 1/1000, and given the catastrophic nature of such an impact it is desirable to detect Virtual Impactors with much lower probability. The method currently being used (Milani et al. 2000b,c) scans hundreds of thousands of possible future close approaches, obtained by propagating all multiple solutions into the future, and applies quasi-linear techniques of target plane analysis to the most suspicious looking ones. However, it is clear that a qualitative, geometric understanding of the way these huge lists of close approaches are organized is a requirement to ensure the reliability of these automatic scans, for which human intervention must be exceptional.

The way to gain such an understanding is to go back to a simpler mathematical setting. 50 years ago Öpik began developing a theory of planetary encounters (Öpik 1951, 1976) based on a piecewise 2-body approach. That is, the small body (asteroid, comet, meteoroid) is considered to be in a heliocentric ellipse until the time of the encounter with some planet, then the dynamics is switched to a planetocentric 2-body orbit, which is (in this approximation) always hyperbolic. Then the standard formulas of 2-body scattering are applied to obtain the initial conditions of a new post-encounter heliocentric orbit. Öpik's theory was successfully used to study the statistical properties of the orbital changes resulting from close approaches, and to some extent it is still in use (Valsecchi & Manara 1997; Valsecchi et al. 1997). The main limitation is that Öpik developed the theory only for the case in which the two orbits, of the small body and of the target planet, are actually touching, that is the Minimum Orbital Intersection Distance (MOID) is zero.

Moreover, the basic theory does not consider that the subsequent encounters of the small body with the same planet (or even with another one) are not independent of the occurrence of the previous ones. The idea, now called resonant return, was implicitly contained in the work of Lexell and Le Verrier on Lexell's comet, was in fact used in spacecraft navigation since the 70s, but was first applied to asteroid close approaches only recently (Marsden 1999; Milani et al. 1999). In fact, recent work by our group has shown that the main organizing principle of the huge lists of close approaches obtained by propagating multiple solutions is that of the resonant return (also, to a lesser extent, non-resonant returns play a role).

With this background, it is possible to state the purpose of this paper in a simple way. We extend Öpik's theory of close encounters to near misses, which can occur also for a finite value of the MOID. Having developed this mathematical tool, we use it to describe how each encounter changes the condition for the next encounter, thus spawning a complicated (in principle, fractal) structure of resonant returns. This we achieve with an analytical theory, whose formulas, although long, can be implemented in a software tool allowing us to explore the geometrical structure of the following encounters in the target plane defining the circumstances of the first encounter.

Note that the analytical formulas we provide are not a replacement for the accurate numerical integrations. However, numerical integrations handle one orbit at a time, and we need to have a global view of everything that could happen as a consequence of a given encounter. Each subsequent return that could lead to an impact defines a keyhole on the target plane of the first encounter, such that an orbit through it would indeed collide with the planet. We give an explicit, semi-analytic description of the keyholes for all possible resonant returns. This allows us to draw the LOV footprint on the target plane, and its intersections with the resonant (and non-resonant) return keyholes provide information on all the subsequent encounters that are possible for an asteroid given the available observations. In the present paper the notation is sometimes different, and the geometrical derivations somewhat simplified, with respect to those contained in Valsecchi (2001), where an earlier, less complete version of this theory is presented.

The paper is organized as follows: in Sect. 2 we describe the classical theory of encounters first introduced by Öpik, and in Sect. 3 we extend it to finite nodal distances, introducing also a time coordinate, which we use in Sect. 4 to compute the initial conditions of a second encounter as functions of those of the first. This mapping between the target planes of the first and of the second encounter, together with its derivatives, allows us to compute the keyholes. In Sect. 5 we give some examples of practical application of the theory, in particular to the well known cases of the asteroids 1997 XF11 and 1999 AN10. Since the details of the analytical development may obscure the discussion of the main ideas, we have collected most of them in the appendix.

### 2 Öpik's theory of encounters

Öpik's theory of close encounters (Öpik 1976) consists of modeling the motion of a small body approaching a planet as a planetocentric two-body scattering. The relative velocity of the small body with respect to the planet defines the direction and speed of the incoming asymptote of the planetocentric hyperbolic orbit. This direction and speed are simple functions of the semimajor axis, eccentricity and inclination (a, e, i) of the heliocentric orbit of the small body, in the approximation neglecting terms of the order of the miss distance. The effect of the encounter is then computed as an instantaneous deflection of the velocity vector in the direction of the outgoing asymptote of the planetocentric hyperbolic orbit, ignoring the perturbation due to the Sun and the time it actually takes for the small body to travel along the curved path that "joins'' the two asymptotes. Interestingly, the errors involved in such an approach are smaller for closer approaches, and Öpik's theory is exact only in the limit of the minimum approach distance going to zero.

#### 2.1 The components of the planetocentric velocity

Let us consider a small body encountering a planet that moves on a circular orbit around the Sun. To simplify the formulae, we use a system of units such that the distance of the planet from the Sun is 1, and the period of the planet is . We also assume that both the mass of the Sun and the gravitational constant are equal to 1. We disregard the mass of the planet in the heliocentric orbit of both the planet and the small body, thus the heliocentric velocity of the planet is also 1.

We use a planetocentric reference frame (X, Y, Z) such that the Y-axis coincides with the direction of motion of the planet, and the Sun is on the negative X-axis. In this system, the components of the unperturbed planetocentric velocity vector of the small body are (Carusi et al. 1990)

and the planetocentric velocity is

This can be rewritten as

where T is the Tisserand parameter with respect to the planet

The direction of the incoming asymptote is defined by two angles, and  (see Fig. 1), such that

and, conversely

 Figure 1: Basic geometric setup of Öpik's theory. The planet is at the origin, moving in the direction of Y, with the Sun at unit distance on the negative X-axis. is the planetocentric velocity vector of the small body (see text). Open with DEXTER

#### 2.2 The b-plane frame

We define the b-plane as the plane orthogonal to containing the center of the planet. The vector extends from the planet to the intersection of the incoming asymptote with the b-plane; is the impact parameter.

We use a planetocentric coordinate system such that are coordinates on the b-plane and the -axis is directed along . The -axis is in the direction opposite to the projection on the b-plane of the heliocentric velocity of the planet. The -axis completes the right-handed reference frame, and is perpendicular to the heliocentric velocity of the planet.

In this way, as will also be seen in more detail later, the shortest segment joining the orbit of the planet and that of the small body, corresponding to the MOID, turns out to be directed along the -axis; this is because this axis is perpendicular, by definition, to both the Y-axis (the direction of motion of the planet), and , the planetocentric velocity of the small body. The -axis, then, can be seen as a "time coordinate'', that is, a shift in the time of arrival of the small body at the b-plane will mean a change only in its  coordinate, and not in . In other words, this coordinate system nicely decouples the two factors governing the possibility of a very close encounter, i.e. distance between orbits and encounter timing, mapping them into, respectively, the  and the -axis.

Following Carusi et al. (1990), we define the angle  by

The transformation from the planetocentric reference frame (X, Y, Z) to the b-plane frame is accomplished by first rotating by an angle about Y, then rotating by  about  (which is perpendicular to the old Y-axis and to ). In matrix notation

 (1)

Similarly, the inverse transformation is accomplished by rotating by  about , then by  about Y. In matrix notation

#### The rotation of

As a consequence of the encounter with the planet, is rotated into , aligned with the outgoing asymptote, without changing the length: U=U'. The deflection angle  between the two vectors is a function of U, the mass of the planet m, and the impact parameter b according to

 = = (2) = (3)

where c=m/U2. This quantity plays in the theory the role of a characteristic length; when b=c, the deflection angle  is . For U=0.5, a typical value for many near-Earth asteroids,  AU, i.e. about 0.29 Earth radii; this means that, unless U is very low, large deflections at close Earth encounters cannot be obtained.

The angles and , defining the direction of the post-encounter velocity vector , can be obtained in terms of , , , by (Carusi et al. 1990)

U is an invariant of the problem and, once it is given, ais a function of only, and does not depend on . In fact, in the geometric setup just described we have fixed the heliocentric distance of the small body, thus fixing its potential energy; to obtain its total energy - that means to obtain a - we need to compute its kinetic energy, i.e. its heliocentric velocity. The latter is the vectorial sum of the heliocentric velocity of the planet, whose components are (0,1,0), and of ; for fixed U, the magnitude of the sum depends only on the angle between the two vectors, i.e. on .

The geometry of the rotation of  is illustrated in Figs. 2 and 3; note that c, that is also an invariant of the problem, due to the conservation of U, has a simple geometric meaning.

 Figure 2: The deflection of . The vectors and lie on, respectively, the incoming and outgoing asymptotes of the planetocentric hyperbola; b (the segment PA) is the impact parameter, is the deflection angle, and the length of the segments AB and BC is c. Open with DEXTER

 Figure 3: Same as Fig. 1, but with also the post-encounter velocity . Open with DEXTER

Öpik's theory of close encounters works as long as the region of space in which the encounter takes place is "small,'' so that the interaction can be thought of as taking place in a point. This assumption breaks down as the Tisserand parameter approaches 3, i.e., when the encounters take place at low planetocentric velocity, and therefore the assumption of a point-like interaction does not hold any more. The theory is inapplicable for a Tisserand parameter exceeding 3.

In principle, this theory could be extended to the case of an elliptic orbit of the planet, preserving most of the formulae, provided the angles  and  are defined with respect to the velocity of the planet at the time of the encounter, which would generally not be orthogonal to the Sun-planet direction.

### 3 Extension of Öpik's theory to near misses

In its original formulation, Öpik's theory of close encounters does not use a complete set of state variables. In a complete formulation, the six orbital elements of the small body have to be transformed into a set also containing six coordinates. The six coordinates used in the extension of Öpik's theory presented in this paper are: U, , , and the time t0 of the crossing of the ecliptic plane by the small body.

In Appendix A.1 we provide formulae to compute the transformation from the orbital elements to this encounter-related coordinate set. In this way the theory can be used to compute actual cases of close encounters with non-zero miss distances, and its results can be compared with those coming from accurate numerical modeling of the motion. However, the formulas we are providing are linearized in the actual miss distance, and therefore cannot be used for very shallow encounters.

In the (X, Y, Z) frame, the motion along the incoming asymptote is

 (4)

where the planetocentric coordinates of the node at time t0 are X0=X(t0), the signed nodal distance, and Y0=Y(t0), which measures how early or late the planet is for the encounter.

The simple formulae of uniform rectilinear motion allow one to decompose the search for the minimum encounter distance into two steps. First, we look for the minimum not constrained by encounter time, to find the MOID. Second, taking into account the time, we find the actual minimum encounter distance for a generic initial condition.

#### 3.1 The local Minimum Orbital Intersection Distance

The Minimum Orbital Intersection Distance (MOID) is the minimum distance between the orbit of the small body and that of the planet. In the case of encounters for which Öpik's theory is applicable (Greenberg et al. 1988) there is the possibility of two local minima of the distance between the orbits, one close to each node, and we call each of these minima a local MOID.

We can derive an expression for the local MOID assuming that the small body travels on the incoming asymptote. From the first Eq. (4), eliminating t-t0, we obtain

setting w=Y-Y0, the square of the distance from the Y-axis (i.e. from the direction of motion of the planet) is

and its derivative is

The derivative is zero at

and the minimum value is

This means that the ratio between the local MOID, close to the node under examination, and X0 is

It is possible to define a signed MOID as ; this expression was the starting point of Bonanno (2000) for his study on the uncertainty of the MOID.

Thus a small body having an encounter at its local MOID crosses the b-plane on the -axis, with and the MOID is just .

#### 3.2 The encounter

In general, the close approach does not take place at the minimum possible distance, i.e. at the MOID. To compute the actual minimum distance as a function of the initial conditions we have to consider that in the (X, Y, Z)-frame, and in the rectilinear motion approximation, the node of the orbit of the small body moves backwards along the Y-axis with speed -1. Thus, if at t=t0 the node is at (X0,Y0,0), we can compute the motion of the small body from the second Eq. (4). The distance from the planet is

So, we take the derivative with respect to t

and find its zero

Then, the minimum approach distance occurs when the small body is at

 (5)

where the distance is

and the coordinates on the b-plane are computed by (1)

So, given and , the coordinates  and  on the b-plane depend only on X0 and Y0; we can invert the relationship and obtain

 (6)

Note that all the computations of this subsection are linearized about the origin both in the -plane and in the (X0,Y0)-plane, that is on the ecliptic. Thus these expressions are applicable neither to very shallow encounters (with large b) nor to encounters with low (e.g., for very low inclination asteroids), for which X0 and Y0 can be large despite the impact parameter being small.

Let the planetocentric angular momentum of the small body be . At the time tb corresponding to the minimum unperturbed distance b, we instantaneously rotate the velocity vector, which is parallel to the incoming asymptote, about , so as to make it parallel to the other asymptote

We also rotate, again instantaneously, the position vector of the small body from  to the one corresponding to the minimum unperturbed distance on the new asymptote

Thus, lies in the post-encounter b-plane normal to . The post-encounter (i.e., rotated) coordinates in the (X, Y, Z) frame can be obtained by applying the appropriate rotations

The new coordinates of the crossing of the ecliptic can be obtained considering that

the time of crossing is then

the Y-coordinate of the crossing is

and the X-coordinate of the crossing is

In the above procedure, we have assumed that the action of the planet on the velocity vector of the small body is instantaneous; one may wonder whether the error involved is large. In fact, as it is clear from Fig. 2, the time actually taken for the small body to go from true planetocentric anomaly to is about 2c/U for large values of b. For small values of b, the length of the arc of hyperbola goes to zero with b. Thus, the approximation involved in neglecting the time actually taken from true anomaly  to  is good enough for the theory of subsequent encounters discussed in the following sections.

#### 3.3 The post-encounter local MOID

We can apply the appropriate rotations by and  to the post-encounter coordinates X'(tb), Y'(tb), Z'(tb), to obtain the coordinates in the post-encounter b-plane as functions of the pre-encounter  and , as in Eq. (1); it can be checked that , as is clear from the geometry of the rotation . The complete formulas are given in Appendix A.2.

Let us discuss here the explicit expression for the new local MOID, i.e. for the  coordinate in the post-encounter b-plane

proceeding as described just before, the result is:

Note that, if , the new local MOID cannot be zero unless the pre-encounter local MOID is already zero. In other words, only initial conditions on the -axis end up on the -axis of the post-encounter b-plane.

Öpik's theory is not applicable for the exactly tangent and coplanar orbits, and its validity becomes questionable for very low values of (Greenberg et al. 1988). This is unfortunate since a small, but non negligible, number of asteroids undergoing encounters with the Earth has . In these cases the interaction may become extremely complex, possibly including temporary satellite captures, as shown by numerical integrations (Carusi et al. 1981): an analytical theory able to handle these cases would be much more complex than the one presented in this paper.

### 4 Resonant returns and keyholes

The orbital period of the planet is , and that of the small body after the encounter is . If the two periods are commensurable, that is a'3/2=k/h with h and k integers, then after h periods of the asteroid k periods of the planet have elapsed, and both the planet and the small body will be back again in the same position of the previous encounter. Such a subsequent encounter is called a resonant return.

Also if the ratio of the period is not exactly k/h, but is close, a subsequent encounter can take place, but the planet will be earlier or later for the encounter than it was at the previous one. The new encounter conditions can be computed as follows.

#### 4.1 Post-encounter propagation

We are now going to compute the post-encounter propagation by using Keplerian heliocentric motion, in order to be able to derive simple analytical expressions. This approximation turns out to be sufficiently accurate to discuss the timing, that is the coordinate, of the next encounter, but not to determine the MOID, that is the  coordinate. This is further discussed in Sect. 4.4.

If the small body follows, after the first encounter, a Keplerian orbit, it will be again at the same node at time

In the planetocentric (X, Y, Z)-frame the node of the orbit of the small body moves backwards along the Y-axis with speed -1(actually, it revolves backwards with period ). At time t'0the Y-component of its distance from the planet was Y'0. We now compute its displacement along the Y-axis between t'0 and t''0, that is . To avoid a discontinuity near , we use the formula

The planetocentric distance of the small body when it is at the node of the resonant return is then

After the propagation, the initial conditions for the new encounter can be computed from:
• the components of the planetocentric velocity U''x=U'x, U''y=U'y, U''z=U'z (since U''=U', , );
• the time of passage at the node t''0;
• the nodal distance X''0=X'0;
• the distance of the node from the planet .
Then, in the X-Y-Z frame, the minimum approach distance is at

For the coordinates corresponding to the minimum approach distance, we have the same expression as Eq. (5), with X'', Y'', Z'', X''0, Y''0, , instead of X, Y, Z, X0, Y0, , , so that the coordinates on the b-plane are

#### 4.2 Solving for a given final semimajor axis

In order to exploit the simple geometry in the b-plane of the solutions of this problem, we discuss it in the framework of classical Öpik's theory disregarding the corrections due to the fact that the heliocentric distance of the small body at close encounter is not 1 (see Appendix A.1).

A given resonance corresponds to certain specific values of a', i.e. of , that we denote with a'0 and . In fact, if we constrain the post encounter orbit in such a way that the ratio of the periods is k/h, then

In the formula of Sect. 2.3 for the post-encounter

the deflection angle  is a function of c and b given by Eqs. (2) and (3). Thus, for given U, , and , this is an equation in the pre-encounter b-plane giving the locus of points leading to a given resonant return. By solving for and using we obtain

 (7)

Replacing b2 with and rearranging terms we obtain

 (8)

This is the equation of a circle centred on the -axis (Valsecchi et al. 2000); if R is the radius of such a circle, and D the value of the -coordinate of its center, its equation is

thus the circle is centred in (0,D) with

Note that the radius of the circle is zero for and . For both D and R tend to infinity:

and the semimajor axis does not change if the encounter takes place on the straight line . If a'>a then and D>0. The half plane contains the circles with D>0, that is the perturbation from the encounter increases the orbital energy because the planet is pulling from ahead, thus a'>a.

In general, the circle intersects the -axis at the values

which represent the extremal values that b can take for a given a'; for one of the two intersections tends to infinity, and the other to . The circle intersects the -axis at

and the maximum value of for which a given  is accessible is R. The maximum value of a' accessible for a given U is for , and is obtained for

and the minimum value of a' is for , and is obtained for

in both cases we must have , that is this happens for zero local MOID.

To take now into account at first order the non zero planetocentric distance of the small body, we compute from  and  the pre-encounter heliocentric distance :

and the post-encounter heliocentric distance :

so that , and are:

Note that is of order c.

Given a, e, i, , we have for and , the values of U and computed for finite (see Appendix A.1):

where U and are the values computed from a, e, i ignoring the first-order correction due to , and , . We have also

For the post-encounter semimajor axis , the value of is:

where is the value computed from a'0 ignoring the first-order correction due to , and .

By replacing in Eq. (7) with and  with

and expanding to first order in and we obtain

Thus, the equation incorporating the first order correction in  is essentially Eq. (8), quadratic in , , c, that defines the b-plane circles, plus terms that are of the third order in the same variables; this essentially follows from the fact that is of the order of c. Since , , c are small numbers, in practice the effect of the third order terms is to distort the shape of the circles to some degree, without altering the overall geometry.

#### 4.3 Mapping to the b-plane of the next encounter

An important advantage of an analytical, albeit approximate, theory of planetary close encounters is that it gives us a key to understand how a given region on the b-plane of an encounter is mapped on the b-plane of the next encounter.

To this end, let us examine the derivatives of the coordinates , on the pre-encounter b-plane of the second encounter. They are functions of the pre-encounter coordinates , on the b-plane of the first encounter (see Appendix A.4).

In particular, let us consider the case in which the encounters are not too close, that is encounters for which . This is not too much of a limiting choice: e.g., let us consider the well-known encounters with the Earth of 1997 XF11 in 2028 and of 1999 AN10 in 2027. For 1997 XF11, , and the local MOID is 0.000 19 AU, so that in the worst case , , while for 1999 AN10, , and the local MOID is 0.000 25 AU, so that in the worst case , . Thus, for both encounters the approximation neglecting c2/b2 is justified, and we can use the approximate expressions given in Appendix A.5. Also neglecting the terms in is justified because of the small value of c for the terrestrial planets, with the exception of unusually slow encounters.

The Jacobian matrix is computed in two steps, as described in the appendix. The first one corresponds to the change from to , that is to the first encounter. The second one corresponds to the change from to , the Keplerian propagation between the encounters; in this step we need to take into account the dependence of , the time delay at the second encounter, on , that is, on a'. On the contrary, the Keplerian propagation does not affect the MOID, that is . Hence the partial derivatives take the form

Considering first the partial derivatives of , we have

The principal parts for of the derivatives appearing in this expression are easily obtained from the formulas of Appendix A.5:

and taking into account that (the deflection is small), this implies

The other partial derivative is

the principal parts of the derivatives are

so that

For the partial derivatives of a similar argument applies. In the derivatives

proceeding term by term the principal parts for are (Appendix A.5)

Thus, as far as the first encounter is concerned, the matrix , has the following structure

 (9)

The important consequence of (9) is that, in the approximations used, the encounter is described by a nearly-area-preserving operator (called  in Appendix A.2).

The partial derivatives of are (disregarding terms in )

 (10)

 (11)

In each of them, the first term comes from the Keplerian propagation, while the second comes from the first encounter, as computed above. The terms describing the Keplerian propagation grow linearly with time due to the presence of h, the number of heliocentric revolutions made by the planet between the first encounter and the next one, and in general can become very large. The divergence of nearby trajectories is expressed by these terms in the majority of cases, that is, for a dynamical evolution dominated by small deflection encounters ( ), and excluding tangential encounters ( ).

The divergence of nearby orbits, having different a', is linear in time. However, sequences of encounters result in multiplicative accumulation of the divergence from each encounter, and thus lead to exponential divergence and chaos, with the maximum Lyapounov exponent proportional to encounter frequency.

As the expressions for and show, the increase of the separation of initially nearby particles on the b-plane of the next encounter contains the factor , defined in Eqs. (10) and (11), given by (disregarding terms in ):

it is clear that, for a'>0, the sign of s is determined by the sign of . Moreover, unless , in which case the use of our theory is questionable, s=0 (for a'>0) if

that happens for

 (12)

If the post-encounter orbit fulfills this condition, then

#### 4.4 Pre-images on the b-plane and keyholes

The term keyhole was introduced by Chodas (1999) to denote the small regions of the b-plane of a specific close encounter such that, if the asteroid passes through one of them, it will hit the planet at a subsequent return. That is, a keyhole is simply one of the possible pre-images of the Earth's cross section on the b-plane. The term keyhole may also be used to indicate a region on the b-plane leading not necessarily to collision, but to a very deep encounter. Thus, a keyhole is linked to a specific value of the post-encounter semimajor axis a', i.e., to the value allowing the occurrence of the next encounter at the given date.

We have seen in this section how to solve for a given value of a', and have discussed the structure of the b-plane circles corresponding to a given a'. However, the algorithm we use includes a Keplerian propagation between encounters, thus only concerns the timing of the next encounter, i.e., the value of the -coordinate, but leaves unchanged the MOID of the next encounter, i.e., the -coordinate.

In fact, the MOID is bound to vary between encounters for two main reasons: on a long time scale, secular perturbations (Gronchi & Milani 2001) make it slowly evolve through the so-called Kozai cycle, or -cycle, while on a shorter time scale significant quasi-periodic variations are caused by planetary perturbations and, for planets with massive satellites, by the displacement of the planet with respect to the center of mass of the planet-satellite system. Figure 4 shows the situation for asteroid 1997 XF11from a precise numerical integration. Large short-period variations superimposed on the secular trend are evident, and these make an analytical modeling of the variation of the MOID problematic.

For the purpose of obtaining the size and shape of an impact keyhole we can, however, just model the secular variation of the MOID as a linear term affecting

 (13)

computing the time derivative of  either from a suitable secular theory for crossing orbits (Gronchi & Milani 2001) or using a value deduced from a numerical integration. The result could then be corrected to take into account the short periodic terms, possibly by taking into account the output of a numerical integration such as the one of Fig. 4. Without such correction, the theory would reliably predict very close encounters (e.g., within a few thousandths of an AU), but not collisions. That is, the theory may be used to predict keyholes in the looser definition, i.e., the points on the b-plane leading to very deep encounters but not necessarily to impacts.
 Figure 4: Time variation of the nodal distance and of the MOID with respect to the Earth for 1997 XF11. Open with DEXTER

The actual keyhole computation proceeds as follows. We want to compute the pre-image of the point of coordinates , , lying on the b-plane of the next encounter, that takes place h revolutions of the small body after the current encounter, on the b-plane of the latter. We start by computing the images of two points with coordinates , and , , say , and , ; furthermore, we choose  and  such that . Note that  is in general slightly different from , but by a very small amount, since in the approximation (13) is a slowly varying function of . We then check whether ; if not, we choose another pair of values for and , until the condition is verified. In practice, the segment parallel to the -axis needs to straddle one of the resonant return circles. At this point, we find the pre-image of a point with , by using, for instance, regula falsi iterations; let us call the coordinates of this point , and the coordinates of its image . If , where  is the radius of the Earth  augmented by the gravitational focusing

then to a good degree of precision the pre-image of the point of coordinates , is the point of coordinates , . The accuracy of this algorithm can be easily checked and improved, if necessary, by an iterative procedure based again on the regula falsi.

The basis of the above procedure is that, while the "horizontal'' distance on the b-plane (i.e., along ) remains essentially unchanged between the two encounters, the "vertical'' one, along , is stretched by a large factor that, as seen in the preceding subsection, comes almost entirely from the propagation between the encounters. Therefore, the "pre-image'' of the Earth on the b-plane of the encounter preceding the collision must resemble a thickened arclet closely following the shape of the circle corresponding to the suitable orbital period; the smallness of the impact keyholes is mainly due to the non-area-preserving nature of the propagation between encounters, i.e. to the large values that can be reached by and .

### 5 Examples and applications

#### 5.1 Resonant returns

To test the theory described in the previous Sections, we apply it to the two already mentioned cases of the encounters with the Earth of 1997 XF11 in 2028 and of 1999 AN10 in 2027; let us start from the latter.

The encounters of a rather large number of "virtual'' 1999 AN10's (fictitious asteroids with orbits compatible with the observational record available in March 1999) were analyzed by Milani et al. (1999), using a realistic gravitational model; around the date of closest approach (7 August 2027), the virtual asteroids have spread into a very thin and long wire extending over a good fraction of an AU. Actually, all the virtual 1999 AN10's at that time occupy a very small region of space, but the small differences in a have accumulated into a substantial spreading in the mean anomaly, similar to what happens in the along-track dispersion of meteoroids in streams.

We model the wire as a very large number of particles all sharing the same values of a, e, i - such that U=0.884, , - the same MOID of 0.000246 AU, so that , and differing only in the time of closest approach; this means that the wire, in the b-plane, is just a segment parallel to the -axis that extends to .

 Figure 5: The outcomes, computed with the extended Öpik theory, of the August 2027 encounter with the Earth of asteroid 1999AN10. Left: final states in the a-e plane (a circle marks the pre-encounter orbit); right: final states in the plane  (difference in time from closest approach; for encounter at the MOID) vs. P (post-encounter orbital period). The left panel is equivalent to Fig. 2 of Milani et al. (2000c), the right panel to Fig. 1 of Milani et al. (1999), in order to show how the extended Öpik theory reproduces the important features of the encounter. For the meaning of the inclined lines, see text. Open with DEXTER

Figure 5 shows the post-encounter parameters of the virtual asteroids; its left panel shows the situation in the a-eplane, and has to be compared with Fig. 2 of Milani et al. (2000c), while the right panel shows the plane (difference in time from closest approach) vs. P (post-encounter orbital period), and has to be compared with Fig. 1 of Milani et al. (1999); in this panel we have also traced the same inclined lines of the corresponding figure of Milani et al. (1999). These lines show the conditions to encounter the Earth at in, from right to left, 2040, 2038, 2036, 2034, 2039 and 2032/2037 (the two leftmost lines, crossing on the axis). The encounters correspond to resonant returns due to, respectively, the 7/13, 6/11, 5/9, 4/7, 7/12 and 3/5 mean motion resonances; the 3/5 resonance appears twice, with the two curves on the left crossing each other, corresponding to returns in 2037 (the less inclined one) and 2032.

The comparison of the two panels of the figure with the corresponding figures of Milani et al. (1999, 2000c) shows that the basic phenomenology of the behaviour of the virtual asteroids wire is well reproduced. Careful inspection of the figures reveals that some details are different. According to Fig. 5, some virtual asteroids can reach the 3/5 resonance, but this possibility was excluded by the computations of Milani et al. This is due, of course, to the inherent approximations taken in the analytic theory, which complements, but cannot fully replace the much more realistic and detailed numerical computations.

#### 5.2 Keyholes

Figure 6 refers again to the August 2027 encounter of 1999AN10. In it, the dotted lines show the circles corresponding to the mean motion resonances 7/13, 10/17, and 11/19, leading to returns in 2040, 2044, and 2046, respectively. Close to each circle, the full lines show the keyholes leading to encounters within 4 Earth radii at each resonant return.

The reason why there are two keyholes for each resonant circle is that, to have an encounter within a very small distance, say  , at a specified date, two conditions must be met, since both  and , the b-plane coordinates at the next encounter, must be small. The vicinity to the resonant circle takes care of , and Eq. (13) gives the condition on  at the previous encounter in order to have small; the latter condition forces the keyholes not to be symmetric with respect to the -axis. Then, the keyholes are near the intersections of a strip, of width approximately equal to , parallel to the -axis, with the given resonant circle; in most cases there are either two intersections or none.

 Figure 6: The dotted lines show the circles, in the b-plane of the August 2027 encounter with the Earth of 1999 AN10, corresponding to the mean motion resonances 7/13 (upper circle), 10/17 (smaller of the two lower circles), 11/19; these resonances lead to returns in 2040, 2044, 2046. The units for the b-plane coordinates  and  are Earth radii, and the Earth is shown at the origin with a circle drawn to scale. Note that the radius of the Earth on the b-plane is larger than 1, in this case by about 8%, due to gravitational focusing. The lines on, or close to, each circle, denote b-plane coordinates resulting in encounters within 4 Earth radii at the resonant returns, assuming a linear drift in the MOID by 0.35 Earth radii per year. To guide the eye, a vertical segment is traced at , corresponding to the actual MOID of about 0.000246 AU. Open with DEXTER

The vertical segment in Fig. 6 represents the wire of virtual asteroids at the encounter. As it can be easily seen, the ranges in  spanned by the keyholes leading to encounters within 4 Earth radii in 2040, 2044 and 2046 are so large that these encounters could take place even if the MOID at the epoch of the 2027 encounter differed by several Earth radii from that of the real 1999 AN10; this is not true only in the case of the 2040 encounter, corresponding to the smallest circle of Fig. 6, since in that case an increase of the MOID of about 1 Earth radius is already sufficient to prevent the intersection of the wire with the resonant circle.

Let us now see how the extension of Öpik's theory described in this paper allows us to trace the contours of collision keyholes. To this purpose, let us examine the October 2028 encounter of 1997 XF11; for this asteroid, U=0.459, and . The pre-images of the Earth for a collision at a resonant return in 2040, are shown in Fig. 7; there are again two keyholes, whose shape has been computed disregarding in and the smaller terms of first order in , (see Appendix A.4).

 Figure 7: Keyholes, in the b-plane of the October 2028 encounter with the Earth of 1997 XF11, for collision at the resonant return in 2040; units for the coordinates are Earth radii. Upper left: keyhole nearer to the Earth, with the latter shown to scale for comparison (the gravitational focusing is here about 30%); upper right: keyhole farther from the Earth; lower left and lower right: enlargements of, respectively, the near and far keyholes shown in the corresponding upper panels. The resonant circle corresponding to the exact resonance is drawn with dots. Open with DEXTER

Both keyholes span, in , roughly the diameter of the Earth augmented by the gravitational focusing; it is in  that a great compression is noticeable. This compression, so to speak, is necessary to compensate for the divergence between nearby trajectories due to the Keplerian propagation (see Sect. 4.3). The divergence is by a factor between 2 200 and 2 600 for the keyhole nearer to the Earth, and by a factor of about 120 for the other one, whose total area is thus much larger. As already remarked, the result of the large compression in  is that the shape of a collision keyhole closely follows the resonant circle to which it is associated, thus looking like an arclet whose thickness is determined by the derivative .

In Fig. 7 it is easily visible that the keyhole farthest from the Earth is somewhat displaced from the resonant circle. This is because the exact resonance corresponds to the circle and the return does not take place at the resonance, but at a slightly different value of a to compensate for the non zero value of  at the first encounter.

Note also that this figure has been traced ignoring the correction due to heliocentric distance of the encounter not being exactly 1 (i.e., due to ); however, we have checked that the terms of the first order in  would displace the position of the resonant return with respect to the circle of Eq. (8) by about 1% of the distance to the Earth.

Finally, a numerical integration, done with the software described in Milani et al. (1999, 2000c), shows that nearby trajectories actually diverge by a factor 134 between 2028 and 2040: as we have just seen, the extended Öpik gives for the corresponding derivative ( ) a value of about 120; this agreement, to within about 10%, is quite satisfactory, given the approximations involved.

#### 5.3 Interrupted returns

When the MOID (i.e. ) of the orbit of a small body encountering a planet is nearly equal the to radius R of the circle corresponding to a given return, there are interesting consequences if also the value of  is close to the coordinate D of the center of the circle. When this is the case, we must distinguish between two possibilities, , and ; to discuss them, we will again make use of the same type of "wire'' approximation used to discuss the previous examples.

In the first case, , we have that the portion of the wire for which is close to the resonance that leads to the return of interest, but remains always on the same side of the resonance itself; as a consequence, in the b-plane of the next encounter the wire cannot straddle the -axis, but has to bend and leave the vicinity of the planet from the same side from which it has approached it. This is the case of interrupted returns discussed by Milani et al. (2000cFig. 7).

In the second case, when , we have interesting consequences for the portion of the wire for which . In this case, we can have a keyhole rather different from those discussed previously, and Fig. 8 illustrates an example obtained by suitably modifying the initial conditions leading to the two collision keyholes of 1997 XF11 discussed previously.

 Figure 8: Keyhole, in the b-plane of the October 2028 encounter with the Earth of 1997 XF11, for collision at the resonant return in 2040; units for the coordinates are Earth radii. The initial conditions have been modified with respect to Fig. 7 in order to force the two keyholes of that figure to coalesce into one. Left: overall view; right: enlargement. Open with DEXTER

On one hand, there is just one keyhole, as the wire is not crossing the circle, but nearly tangent to it; on the other hand, in this keyhole the divergence of orbits along  is much smaller than in the other cases.

To understand this, let us consider Eq. (11), giving ; in it, the term passes through zero at the tangency of the resonant circle, implying that also the largest term in passes through zero, while the remaining term, , is close to 1; this means that, somewhere in the vicinity, the entire derivative is null, or very small, with the consequence that the keyhole can extend substantially in . This is clearly visible in Fig. 8, where the keyhole spans almost 40 Earth radii in , and is in contrast with "normal'' keyholes, like those of Fig. 7 that, along , are thinner than the b-plane image of the Earth by a large factor (of the order of that, as we have seen, in those cases is rather large).

Concerning , the situation is essentially governed by the value of that, as is easily visible from Eq. (10), can be rather large, so that very small changes in  produce large changes , i.e. in the timing of the next encounter, leading to the very small thickness visible in Fig. 8.

In these nearly tangent cases, especially if the condition given by Eq. (12) is fulfilled, we may have large keyholes, which represent a new, unsuspected phenomenon. The assessment of their importance for collisions in the solar system goes beyond the scope of the present paper, and is the topic of continuing research.

Acknowledgements
The authors are grateful to A. Carusi and P. Chodas for useful discussions on the subject of this paper, and are particularly indebted to the referee A. Morbidelli, whose constructive criticism led to a significant improvement of the paper. The work described here was completed at the Istituto di Astrofisica Spaziale e Fisica Cosmica of CNR and at the University of Pisa, under a contract with the Italian Space Agency (ASI), and at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.

## References

### Appendix A: Computational details

#### A.1 Pre-encounter state vector

The pre-encounter state vector has components (U, , , , , t0). To compute these components from the heliocentric elements (a, e, i, , , f) we have to consider where the encounter of interest takes place.

Let us call the longitude of the planet at time t0; then, if at that time , the encounter takes place at the ascending node, in the post-perihelion branch of the orbit if , and in the pre-perihelion branch otherwise. If at time t0 we have , the encounter takes place at the descending node, in the post-perihelion branch of the orbit if , and in the pre-perihelion branch otherwise.

Thus, taking (6) into account, and neglecting terms of second order in the miss distance at the node, that is , at the ascending node we have

with

whereas at the descending node we have

with

For the remaining components, let us see first what Öpik's theory, which assumes that the position of the small body coincides with that of the planet, gives; in this case, we have for U, and :

with the quadrant of determined from the previous considerations. The node-crossing time t0 can be computed by standard two-body formulas, including Kepler's equation.

Since

a simple expression for is in terms of U, a

Finally, a, e, i can be obtained from U, ,

The non-zero planetocentric distance of the small body affects the computation of U, , from a, e, i only if the heliocentric distance r of the small body is not 1: if r=1, Öpik's classical expressions are valid. If , an explicit computation for , and , i.e. the values of U, and corresponding to such a situation, gives, at the first order in ,

where U, , are the values computed for , and

On the other hand, if we have a generic vector  applied in , we obtain , , from the components of  with the following expressions, again at first order in :

where a, e, i are the values computed for ; as shown by the formulae relating a to  , care must be taken in the case of comets, where can easily be non negligible with respect to 1/a. For  and the formulas given above are correct to first order in  if the angles  and  are computed to the same order.

#### A.2 Post-encounter state vector

The close encounter can be seen as an operator that maps the pre-encounter state vector , with components , into the post-encounter one , with components

The components of the post-encounter state vector, as functions of the pre-encounter state vector components (with ), are

U' = U

#### A.3 Propagation to the next encounter

A purely Keplerian propagation can be seen as an operator that maps the post-encounter state vector in the pre-next-encounter one , with components (U'', , , , , t''0)

The transformation is given by:

U'' = U'

where the key parameter is the post-encounter semimajor axis a'. Note that the only components of  that differ from those of  are the -component on the post encounter b-plane  and the time of nodal passage t'0.

Classical Öpik's theory, which assumes that the position of the small body coincides with that of the planet, gives

For finite miss distances the first order correction to Öpik's expression would be, as we have seen before:

with

#### A.4 Derivatives

To write the derivatives, let us consider the total evolution from before the first encounter to before the second as the composition of the first encounter and of the Keplerian heliocentric propagation, as seen in the previous subsections

The matrix of derivatives is then

where the propagation derivatives matrix is regarded as a function of the pre-first encounter variables (U, , , , , t0).

For encounters computed with the extended Öpik theory,

and, for a purely Keplerian heliocentric propagation between encounters,

so that the composite matrix has the following structure

We are interested, in particular, in the submatrix , giving the derivatives of the b-plane coordinates of the second encounter with respect to those of the first encounter. The Keplerian propagation does not affect the MOID, i.e., , so that the first row of the matrix is

On the other hand, the Keplerian propagation affects  only through a', and U is an invariant of the motion. Thus the derivatives of  with respect to  and  have the following structure

with the derivatives of with respect to and  given by

Note that the second and the third term in the derivatives of the Keplerian propagation are of the first order in , , so that they are, in many situations of practical interest, much smaller that the first term, and can be ignored.

The individual derivatives entering the above expressions are

#### A.5 Approximate equations for small deflection

When dealing with encounters of NEAs with the terrestrial planets, we have that in many cases of interest c=m/U2 is very small. This happens because for these planets , and typical values of U are around 0.5. On the other hand, for an Earth-encountering asteroid with , a rather typical value, the minimum value of b that avoids collision is

so that for most encounters the condition is verified (the only exceptions are extremely deep encounters with very low U.)

If we can assume , most expressions of this paper can be approximated by much simpler formulas neglecting terms. We have collected these approximate formulas below.

The post-encounter state vector:

U' = U

Note that, in this approximation, unless is very small, , that is the local MOID is nearly conserved. This corresponds to the intuitive idea that an encounter can hardly change the MOID, because a significant deflection only occurs near the MOID point, and the post-encounter orbit has to pass from the encounter point. Note also that

Post-encounter semimajor axis a':

Partial derivatives:

#### A.6 Approximate equations for large deflection

Encounters of interplanetary objects with the giant planets are characterized by values of c larger than in the cases seen before, because of the large masses of these planets. As an example, the 1779 encounter of D/Lexell with Jupiter was characterized by ; considering that b could have been as low as 0.0003 (Le Verrier 1857), we have that for the closest approach , a case in which the following simplified expressions hold.

Post-encounter state vector:

U' = U

Note that, also in this approximation, unless is very small, , that is the local MOID is nearly conserved. Note also that

Post-encounter semimajor axis a':

Partial derivatives: