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

Resonant returns to close approaches: Analytical theory

G. B. Valsecchi¹ - A. Milani² - G. F. Gronchi³ - S. R. Chesley⁴

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 $\mu$, whose value could vary between $\pm$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 $\mu$ 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 $\mu$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 $\mu$. 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 XF₁₁ and 1999 AN₁₀. 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 $2\pi$. 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 $\vec{U}$of the small body are (Carusi et al. 1990)

$$\left[\begin{array}{c} U_x \\ U_y \\ U_z \end{array}\right] =... ...e^2)}\cos{i}-1 \\ \pm\sqrt{a(1-e^2)}\sin{i} \end{array}\right]$$

and the planetocentric velocity is

$$U=\sqrt{3-\frac{1}{a}-2\sqrt{a(1-e^2)}\cos{i}}.$$

This can be rewritten as

$$U=\sqrt{3-T}$$

where T is the Tisserand parameter with respect to the planet

$$T=\frac{1}{a}+2\sqrt{a(1-e^2)}\cos{i}.$$

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

$$\left[\begin{array}{c} U_x \\ U_y \\ U_z \end{array}\right] =... ...\phi \\ U\cos\theta \\ U\sin\theta\cos\phi \end{array}\right]$$

and, conversely

$$\left[\begin{array}{c} \cos\theta \\ \tan\phi \end{array}\right] =\left[\begin{array}{c} U_y/U \\ U_x/U_z \end{array}\right].$$

	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. $\vec{U}$ 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 $\vec{U}$ containing the center of the planet. The vector $\vec{b}$ extends from the planet to the intersection of the incoming asymptote with the b-plane; $b=\vert\vec b\vert$ is the impact parameter.

We use a planetocentric coordinate system $(\xi, \eta,\zeta)$ such that $(\xi,\zeta)$ are coordinates on the b-plane and the $\eta$-axis is directed along $\vec{U}$. The $\zeta$-axis is in the direction opposite to the projection on the b-plane of the heliocentric velocity of the planet. The $\xi$-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 $\xi$-axis; this is because this axis is perpendicular, by definition, to both the Y-axis (the direction of motion of the planet), and $\vec{U}$, the planetocentric velocity of the small body. The $\zeta$-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 $\zeta$ coordinate, and not in $\xi$. 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 $\xi$ and the $\zeta$-axis.

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

$$\left[\begin{array}{c} b\sin{\psi} \\ b\cos{\psi} \end{array}\right]= \left[\begin{array}{c} \xi \\ \zeta \end{array}\right].$$

The transformation from the planetocentric reference frame (X, Y, Z) to the b-plane frame $(\xi, \eta,\zeta)$ is accomplished by first rotating by an angle $-\phi$ about Y, then rotating by $-\theta$ about $\xi$ (which is perpendicular to the old Y-axis and to $\vec{U}$). In matrix notation

$$\left[\begin{array}{c} \xi \\ \eta \\ \zeta \end{array}\right... ...(-\phi) \left[\begin{array}{c} X \\ Y \\ Z \end{array}\right].$$ (1)

Similarly, the inverse transformation is accomplished by rotating by $\theta$ about $\xi$, then by $\phi$ about Y. In matrix notation

$$\left[\begin{array}{c} X \\ Y \\ Z \end{array}\right] ={\bf\h... ...left[\begin{array}{c} \xi \\ \eta \\ \zeta \end{array}\right].$$

   
The rotation of $\vec{U}$

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

  

$$\tan\frac{\gamma}{2} \frac{m}{bU^2}=\frac{c}{b}$$ =

$$\cos\gamma \frac{b^2U^4-m^2}{b^2U^4+m^2}=\frac{b^2-c^2}{b^2+c^2}$$ = (2)

$$\sin\gamma \frac{2mbU^2}{b^2U^4+m^2}=\frac{2bc}{b^2+c^2};$$ = (3)

where c=m/U². This quantity plays in the theory the role of a characteristic length; when b=c, the deflection angle $\gamma$ is $\pi/2$. For U=0.5, a typical value for many near-Earth asteroids, $c=1.22\times 10^{-5}$ 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 $\theta'$ and $\phi'$, defining the direction of the post-encounter velocity vector $\vec{U'}$, can be obtained in terms of $\theta$, $\phi$, $\gamma$, $\psi$ by (Carusi et al. 1990)

$$\cos\theta' = \cos\theta\cos\gamma+\sin\theta\sin\gamma\cos\psi$$

$$\tan(\phi-\phi') = \frac{\sin\gamma\sin\psi} {\sin\theta\cos\gamma-\cos\theta\sin\gamma\cos\psi}$$

$$\tan\phi' = \frac{\tan\phi-\tan(\phi-\phi')} {1+\tan\phi\tan(\phi-\phi')}$$

$$\cos\phi' = \frac{1}{\sqrt{1+\tan^2\phi'}}$$

$$\sin\phi' = \cos\phi'\tan\phi'.$$

U is an invariant of the problem and, once it is given, ais a function of $\cos\theta$ only, and does not depend on $\phi$. 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 $\vec{U}$; for fixed U, the magnitude of the sum depends only on the angle between the two vectors, i.e. on $\theta$.

The geometry of the rotation of $\vec{U}$ 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 $\vec{U}$. The vectors $\vec{U}$ and $\vec{U'}$ lie on, respectively, the incoming and outgoing asymptotes of the planetocentric hyperbola; b (the segment PA) is the impact parameter, $\gamma$ 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 $\vec{U'}$.
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 $\theta$ and $\phi$ 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, $\theta$, $\phi$, $\xi$, $\zeta$ and the time t₀ 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

$$\left[\begin{array}{c} X(t) \\ Y(t) \\ Z(t) \end{array}\right... ...a(t-t_0)+Y_0 \\ U\sin\theta\cos\phi(t-t_0) \end{array}\right],$$ (4)

where the planetocentric coordinates of the node at time t₀ are X₀=X(t₀), the signed nodal distance, and Y₀=Y(t₀), 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-t₀, we obtain

$$\left[\begin{array}{c} X \\ Z \end{array}\right] =\left[\begi... ...U_x/U_y)~(Y-Y_0)+X_0 \\ (U_z/U_y)~(Y-Y_0) \end{array}\right];$$

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

$$D_y(w)^2=X^2+Z^2=\frac{U_x^2+U_z^2}{U_y^2}w^2+2\frac{U_x}{U_y}X_0w+X_0^2,$$

and its derivative is

$$\frac{{\rm d}\left(D_y^2\right)}{{\rm d}w}=2\frac{U_x^2+U_z^2}{U_y^2}w+2\frac{U_x}{U_y}X_0.$$

The derivative is zero at

$$w=-\frac{U_xU_y}{U_x^2+U_z^2}X_0,$$

and the minimum value is

$$\min{D_y^2}=X_0^2\left[1-\frac{U_x^2}{U_x^2+U_z^2}\right]=X_0^2\cos^2\phi.$$

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

$$\min{D_y}=\vert X_0\cos\phi\vert.$$

It is possible to define a signed MOID as $X_0\cos\phi$; 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 $\xi$-axis, with $\xi=X_0 \cos \phi$ and the MOID is just $\vert\xi\vert$.

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=t₀ the node is at (X₀,Y₀,0), we can compute the motion of the small body from the second Eq. (4). The distance from the planet is

$$\begin{eqnarray*}D^2 & = & X^2+Y^2+Z^2 \\ & = & U^2t^2+2U(X_0\sin\theta\sin\phi... ...^2t_0^2 -2U(X_0\sin\theta\sin\phi+Y_0\cos\theta)t_0+X_0^2+Y_0^2. \end{eqnarray*}$$

So, we take the derivative with respect to t

$$\frac{{\rm d}(D^2)}{{\rm d}t}=2U^2t-2U(Ut_0-X_0\sin\theta\sin\phi-Y_0\cos\theta)$$

and find its zero

$$t_b=t_0-\frac{X_0\sin\theta\sin\phi+Y_0\cos\theta}{U}\cdot$$

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

$$\left[\begin{array}{c} X \\ Y \\ Z \end{array}\right] =\left[... ...a\sin\phi+Y_0\cos\theta)\sin\theta\cos\phi \end{array}\right],$$ (5)

where the distance is

$$\begin{eqnarray*}D\equiv b & = & \sqrt{X^2+Y^2+Z^2} =\sqrt{X_0^2\cos^2\phi+(X_0\cos\theta\sin\phi-Y_0\sin\theta)^2}, \end{eqnarray*}$$

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

$$\left[\begin{array}{c} \xi \\ \zeta \end{array}\right] =\left... ...hi \\ X_0\cos\theta\sin\phi-Y_0\sin\theta \end{array}\right].$$

So, given $\theta$ and $\phi$, the coordinates $\xi$ and $\zeta$ on the b-plane depend only on X₀ and Y₀; we can invert the relationship and obtain

$$\left[\begin{array}{c} X_0 \\ Y_0 \end{array}\right] =\left[\... ...\cos\theta\tan\phi-\zeta\right)/\sin\theta \end{array}\right].$$ (6)

Note that all the computations of this subsection are linearized about the origin both in the $(\xi,\zeta)$-plane and in the (X₀,Y₀)-plane, that is on the ecliptic. Thus these expressions are applicable neither to very shallow encounters (with large b) nor to encounters with low $U_z/U=\sin\theta \cos\phi$ (e.g., for very low inclination asteroids), for which X₀ and Y₀ can be large despite the impact parameter being small.

Let the planetocentric angular momentum of the small body be $\vec{h}=\vec{b}\times\vec{U}$. At the time t_b corresponding to the minimum unperturbed distance b, we instantaneously rotate the velocity vector, which is parallel to the incoming asymptote, about $\vec{h}$, so as to make it parallel to the other asymptote

$$\vec{U'}\equiv\left[\begin{array}{c} U\sin\theta'\sin\phi' \... ...i' \end{array}\right]={\bf\hat{R}}_{\vec{h}}(\gamma)\vec{U}.$$

We also rotate, again instantaneously, the position vector of the small body from $\vec{b}$ to the one corresponding to the minimum unperturbed distance on the new asymptote $\vec{b'}$

$$\vec{b'}\equiv\left[\begin{array}{c} \xi' \\ 0 \\ \zeta' \end{array}\right]={\bf\hat{R}}_{\vec{h}}(\gamma)\vec{b}.$$

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

$$\left[\begin{array}{c} X'(t_b) \\ Y'(t_b) \\ Z'(t_b) \end{arr... ...left[\begin{array}{c} \xi \\ \eta \\ \zeta \end{array}\right].$$

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

$$\left[\begin{array}{c} X'(t_b) \\ Y'(t_b) \\ Z'(t_b) \end{arr... ... U'_y(t_b-t'_0)+Y'(t'_0) \\ U'_z(t_b-t'_0) \end{array}\right];$$

the time of crossing is then

$$t'_0=t_b-\frac{Z'(t_b)}{U'_z},$$

the Y-coordinate of the crossing is

$$Y'(t'_0)\equiv Y'_0=Y'(t_b)-U'_y\left(t_b-t'_0\right),$$

and the X-coordinate of the crossing is

$$X'(t'_0)\equiv X'_0=X'(t_b)-U'_x\left(t_b-t'_0\right).$$

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 $f_g=-\gamma/2$ to $f_g=+\gamma/2$ 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 $-\gamma/2$ to $+\gamma/2$ 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 $\theta'$ and $\phi'$ to the post-encounter coordinates X'(t_b), Y'(t_b), Z'(t_b), to obtain the coordinates in the post-encounter b-plane as functions of the pre-encounter $\xi$ and $\zeta$, as in Eq. (1); it can be checked that $\eta'=0$, as is clear from the geometry of the rotation ${\bf\hat{R}}_{\vec{h}}(\gamma)$. 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 $\xi'$ coordinate in the post-encounter b-plane

$$D'_y\equiv\xi'=X'_0\cos\phi';$$

proceeding as described just before, the result is:

$$\xi'=\frac{\left(b^2+c^2\right)\xi\sin\theta} {\sqrt{\left[\l... ...2\right)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2}}\cdot$$

Note that, if $\sin\theta \neq 0$, the new local MOID cannot be zero unless the pre-encounter local MOID $D_y\equiv\xi$ is already zero. In other words, only initial conditions on the $\zeta$-axis end up on the $\zeta$-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 $\sin\theta$ (Greenberg et al. 1988). This is unfortunate since a small, but non negligible, number of asteroids undergoing encounters with the Earth has $\sin\theta\approx0$. 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 $2\pi$, and that of the small body after the encounter is $2\pi a'^{3/2}$. 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 $\zeta$coordinate, of the next encounter, but not to determine the MOID, that is the $\xi$ 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

$$t''_0=t'_0+h\cdot2\pi a'^{3/2}.$$

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 $2\pi$). At time t'₀the Y-component of its distance from the planet was Y'₀. We now compute its displacement along the Y-axis between t'₀ and t''₀, that is $\delta Y_0 + 2\pi k = t_0''-t_0'$. To avoid a discontinuity near $\delta Y_0=0$, we use the formula

$$\delta Y_0 = -\left[{\rm mod}\left(t''_0-t'_0+\pi,2\pi\right)-\pi\right].$$

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

$$D(t''_0)=\sqrt{{X'_0}^2+(Y'_0+\delta Y_0)^2}.$$

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', $\theta''=\theta'$, $\phi''=\phi'$);
• the time of passage at the node t''₀;
• the nodal distance X''₀=X'₀;
• the distance of the node from the planet $Y''_0=Y'_0+\delta Y_0$.

Then, in the X-Y-Z frame, the minimum approach distance is at

$$t''_b=t''_0-\frac{X''_0\sin\theta'\sin\phi'+Y''_0\cos\theta'}{U}\cdot$$

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

$$\left[\begin{array}{c} \xi'' \\ \zeta'' \end{array}\right] =\... ...X''_0\cos\theta'\sin\phi'-Y''_0\sin\theta' \end{array}\right].$$

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 $\theta'$, that we denote with a'₀ and $\theta'_0$. In fact, if we constrain the post encounter orbit in such a way that the ratio of the periods is k/h, then

$$a'_0 = \sqrt[3]{\frac{k^2}{h^2}}$$

$$\cos\theta'_0 = \frac{1-U^2-1/a'_0}{2U}\cdot$$

In the formula of Sect. 2.3 for the post-encounter $\theta'$

$$\cos\theta'_0=\cos\theta\cos\gamma+\sin\theta\sin\gamma\cos\psi,$$

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

$$\zeta =\frac{\left(b^2+c^2\right)\cos\theta'_0-\left(b^2-c^2\right)\cos\theta} {2c\sin\theta}\cdot$$ (7)

Replacing b² with $\xi^2+\zeta^2$ and rearranging terms we obtain

$$\xi^2 + \zeta^2-\frac{2c\sin\theta}{\cos\theta'_0-\cos\theta}... ...ac{c^2(\cos\theta'_0+\cos\theta)}{\cos\theta'_0-\cos\theta}=0.$$ (8)

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

$$\xi^2 + \zeta^2-2D\zeta+D^2=R^2;$$

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

$$D=\frac{c\sin\theta}{\cos\theta'_0-\cos\theta}$$

and has radius

$$\begin{eqnarray*}R & = & \left\vert\frac{c\sin\theta'_0}{\cos\theta'_0-\cos\theta}\right\vert \cdot \end{eqnarray*}$$

Note that the radius of the circle is zero for $\theta'_0=0$ and $\theta'_0=\pi$. For $\theta'_0\rightarrow\theta$ both D and R tend to infinity:

$$\begin{eqnarray*}D & = & \frac{c\sin\theta}{\cos(\theta+\delta\theta)-\cos\theta... ...\frac{c}{\delta\theta}-\frac{c\cos\theta}{\sin\theta}\right\vert \end{eqnarray*}$$

and the semimajor axis does not change if the encounter takes place on the straight line $\zeta=c \cot\theta$. If a'>a then $\cos\theta_0'>\cos\theta$ and D>0. The half plane $\zeta>c \cot\theta$ 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 $\zeta$-axis at the values

$$\zeta=D\pm R=\frac{c(\sin\theta\pm\sin\theta'_0)}{\cos\theta'_0-\cos\theta},$$

which represent the extremal values that b can take for a given a'; for $\theta'_0\rightarrow\theta$ one of the two intersections tends to infinity, and the other to $c \cot\theta$. The circle intersects the $\xi$-axis at

$$\xi=\pm c\sqrt{\frac{\cos\theta+\cos\theta'_0}{\cos\theta-\cos\theta'_0}},$$

and the maximum value of $\vert\xi\vert$ for which a given $\theta'_0$ is accessible is R. The maximum value of a' accessible for a given U is for $\theta'_0=0$, and is obtained for

$$\zeta=\frac{c\sin\theta}{1-\cos\theta},$$

and the minimum value of a' is for $\theta'_0=\pi$, and is obtained for

$$\zeta=-\frac{c\sin\theta}{1+\cos\theta};$$

in both cases we must have $\xi=0$, 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 $\xi$ and $\zeta$ the pre-encounter heliocentric distance $r=1+\epsilon$:

$$r^2 = 1+2X+X^2+Y^2+Z^2 =1+b^2+2(\zeta\cos\theta\sin\phi+\xi\cos\phi)$$

$$r = 1+\zeta\cos\theta\sin\phi+\xi\cos\phi,$$

and the post-encounter heliocentric distance $r'=1+\epsilon'$:

$$r'^2 = 1+2X'+X'^2+Y'^2+Z'^2 =1+b^2+\frac{4b^2c\sin\theta\sin\phi+2(b^2-c^2) (\zeta\cos\theta\sin\phi+\xi\cos\phi)}{b^2+c^2}$$

$$r' = 1+\frac{2b^2c\sin\theta\sin\phi+(b^2-c^2) (\zeta\cos\theta\sin\phi+\xi\cos\phi)}{b^2+c^2},$$

so that $\epsilon$, $\epsilon'$ and $\epsilon'-\epsilon$ are:

$$\epsilon = \zeta\cos\theta\sin\phi+\xi\cos\phi$$

$$\epsilon' = \frac{2b^2c\sin\theta\sin\phi+\left(b^2-c^2\right) \left(\zeta\cos\theta\sin\phi+\xi\cos\phi\right)}{b^2+c^2}$$

$$\epsilon'-\epsilon = \frac{2c[b^2\sin\theta\sin\phi -c\left(\zeta\cos\theta\sin\phi+\xi\cos\phi\right)]}{b^2+c^2}\cdot$$

Note that $\epsilon'-\epsilon$ is of order c.

Given a, e, i, $\epsilon$, we have for $U_\epsilon$ and $\cos\theta_\epsilon$, the values of U and $\cos\theta$ computed for finite $\epsilon$ (see Appendix A.1):

$$U_\epsilon = U-\epsilon\frac{2U^3\cos\theta+4U^2\cos^2\theta +U^2+1}{U^2+2U\cos\theta+1} =U-\epsilon f$$

$$\cos\theta_\epsilon = \cos\theta\left[1+\epsilon\frac{4U^2\c... ...2-1}{U\left(U^2+2U\cos\theta+1\right)} =\cos\theta+\epsilon g,$$

where U and $\cos\theta$ are the values computed from a, e, i ignoring the first-order correction due to $\epsilon$, and $f=f(U,\cos\theta)$, $g=g(U,\cos\theta)$. We have also

$$\begin{eqnarray*}\sin\theta_\epsilon & = & \sin\theta-\epsilon g\cos\theta. \end{eqnarray*}$$

For the post-encounter semimajor axis $a'_0=\sqrt[3]{k^2/h^2}$, the value of $\cos\theta'_{0\epsilon'}$ is:

$$\begin{eqnarray*}\cos\theta'_{0\epsilon'} & = & \cos\theta'_0+\epsilon'g', \end{eqnarray*}$$

where $\cos\theta'_0$ is the value computed from a'₀ ignoring the first-order correction due to $\epsilon'$, and $g'=g(U,\cos\theta'_0)$.

By replacing in Eq. (7) $\theta$ with $\theta_\epsilon$ and $\theta'_0$ with $\theta'_{0\epsilon'}$

$$\left(b^2+c^2\right)\cos\theta'_{0\epsilon'}-\left(b^2-c^2\right)\cos\theta_\epsilon -2c\zeta\sin\theta_\epsilon = 0$$

and expanding to first order in $\epsilon$ and $\epsilon'$ we obtain

$$\begin{eqnarray*}\xi^2+\zeta^2 &\!\!-\!\!&\frac{2c\zeta\sin\theta}{\cos\theta'_0... ...sin\phi+\xi\cos\phi)g\cos\theta} {\cos\theta'_0-\cos\theta} = 0. \end{eqnarray*}$$

Thus, the equation incorporating the first order correction in $\epsilon$ is essentially Eq. (8), quadratic in $\xi$, $\zeta$, 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 $\epsilon'-\epsilon$ is of the order of c. Since $\xi$, $\zeta$, 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 $\xi''$, $\zeta''$ on the pre-encounter b-plane of the second encounter. They are functions of the pre-encounter coordinates $\xi$, $\zeta$ 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 $c^2\ll b^2$. This is not too much of a limiting choice: e.g., let us consider the well-known encounters with the Earth of 1997 XF₁₁ in 2028 and of 1999 AN₁₀ in 2027. For 1997 XF₁₁, $c\approx1.4\times 10^{-5}$, and the local MOID is 0.000 19 AU, so that in the worst case $c/b\approx0.076$, $c^2/b^2\approx0.006$, while for 1999 AN₁₀, $c\approx3.9\times 10^{-6}$, and the local MOID is 0.000 25 AU, so that in the worst case $c/b\approx0.016$, $c^2/b^2\approx0.000~3$. Thus, for both encounters the approximation neglecting c²/b² is justified, and we can use the approximate expressions given in Appendix A.5. Also neglecting the terms in $\epsilon'-\epsilon$ is justified because of the small value of c for the terrestrial planets, with the exception of unusually slow encounters.

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

$$\begin{eqnarray*}\frac{\partial\xi''}{\partial\xi} & = & \frac{\partial\xi'}{\p... ...heta'}{\partial\zeta} +\frac{\partial\zeta'}{\partial\zeta}\cdot \end{eqnarray*}$$

Considering first the partial derivatives of $\xi'$, we have

$$\frac{\partial\xi'}{\partial\xi}=\frac{\partial X'_0}{\partial\xi}\cos\phi' +X'_0\frac{\partial\cos\phi'}{\partial\xi}\cdot$$

The principal parts for $c/b\to0$ of the derivatives appearing in this expression are easily obtained from the formulas of Appendix A.5:

$$\frac{\partial X'_0}{\partial\xi} = \frac{\sin\theta +\mathc... ...os\phi +\mathcal{O}(c/b)} =\frac{1}{\cos\phi}+\mathcal{O}(c/b)$$

$$X_0'~\frac{\partial\cos\phi'}{\partial\xi} = \mathcal{O}(c/b),$$

and taking into account that $\cos\phi'=\cos\phi+\mathcal{O}(c/b)$ (the deflection is small), this implies

$$\frac{\partial\xi'}{\partial\xi}=1+\mathcal{O}(c/b).$$

The other partial derivative is

$$\frac{\partial\xi'}{\partial\zeta}=\frac{\partial X'_0} {\partial\zeta}\cos\phi'+X'_0\frac{\partial\cos\phi'}{\partial\zeta},$$

the principal parts of the derivatives are

$$\frac{\partial X'_0}{\partial\zeta} = \frac{\mathcal{O}(c/b)} {\sin\theta \cos\phi +\mathcal{O}(c/b)} =\mathcal{O}(c/b)$$

$$X_0'~\frac{\partial\cos\phi'}{\partial\zeta} = \mathcal{O}(c/b),$$

so that

$$\frac{\partial\xi'}{\partial\zeta}=\mathcal{O}(c/b).$$

For the partial derivatives of $\zeta'$ a similar argument applies. In the derivatives

$$\frac{\partial\zeta'}{\partial\xi} = \frac{\partial X'_0}{\p... ...ial\xi}\sin\theta' +Y'_0\frac{\partial\sin\theta'}{\partial\xi}$$

$$\frac{\partial\zeta'}{\partial\zeta} = \frac{\partial X'_0}{... ...ta}\sin\theta' +Y'_0\frac{\partial\sin\theta'}{\partial\zeta},$$

proceeding term by term the principal parts for $c/b\to0$ are (Appendix A.5)

$$\frac{\partial\zeta'}{\partial\xi} = \frac{\cos\theta'\sin\p... ...theta \cos\phi}\sin\theta' +\mathcal{O}(c/b) +\mathcal{O}(c/b)$$

$$\frac{\partial\zeta'}{\partial\zeta} = \mathcal{O}(c/b)+\mat... ...}+\mathcal{O}(c/b)+\mathcal{O}(c/b) \\ = 1+\mathcal{O}(c/b).$$

Thus, as far as the first encounter is concerned, the matrix $\partial(\xi',\zeta')/\partial(\xi,\zeta)$, has the following structure

$$\left[\begin{array}{cc} \frac{\partial\xi'}{\partial\xi} & \... .../b) \\ \mathcal{O}(1) & 1+\mathcal{O}(c/b)\end{array}\right].$$ (9)

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

The partial derivatives of $\zeta''$ are (disregarding terms in $\epsilon'$)

$$\frac{\partial\zeta''}{\partial\xi} = h\cdot\frac{2\pi a'^{5... ...l\cos\theta'}{\partial\xi} +\frac{\partial\zeta'}{\partial\xi}$$ (10)

$$\frac{\partial\zeta''}{\partial\zeta} = h\cdot\frac{2\pi a'^... ...ta'}{\partial\zeta} +\frac{\partial\zeta'}{\partial\zeta}\cdot$$ (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 ( $c^2\ll b^2$), and excluding tangential encounters ( $\sin\theta\not\approx0$).

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 $\partial\zeta''/\partial\xi$ and $\partial\zeta''/\partial\zeta$ show, the increase of the separation of initially nearby particles on the b-plane of the next encounter contains the factor $s(U', \theta')$, defined in Eqs. (10) and (11), given by (disregarding terms in $\epsilon'$):

$$s(U',\theta')= \frac{2\pi a'^{5/2}\left[U'\cos^2\theta'+\cos\... ...right]} {\left(1-U'^2-2U'\cos\theta'\right)^{5/2}\sin\theta'};$$

it is clear that, for a'>0, the sign of s is determined by the sign of $[U'\cos^2\theta'+\cos\theta'(1-U'^2)-3U']/\sin\theta'$. Moreover, unless $\sin\theta'=0$, in which case the use of our theory is questionable, s=0 (for a'>0) if

$$U'\cos^2\theta'+\cos\theta'(1-U'^2)-3U'=0,$$

that happens for

$$\cos\theta'=\frac{U'^2-1\pm\sqrt{U^4+10U^2+1}}{2U'}\cdot$$ (12)

If the post-encounter orbit fulfills this condition, then

$$\frac{\partial\zeta''}{\partial\xi} = \frac{\partial\zeta'}{\partial\xi}$$

$$\frac{\partial\zeta''}{\partial\zeta} = \frac{\partial\zeta'}{\partial\zeta}\cdot$$

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 $\zeta$-coordinate, but leaves unchanged the MOID of the next encounter, i.e., the $\xi$-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 $\omega$-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 XF₁₁from 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 $\xi''$

$$\xi''=\xi'+\frac{\rm d\xi}{{\rm d}t}\left(t''_0-t'_0\right),$$ (13)

computing the time derivative of $\xi$ 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 $\xi''$, $\zeta''=0$, 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 $\xi_1$, $\zeta_1$ and $\xi_2$, $\zeta_2$, say $\xi_1''$, $\zeta''_1$ and $\xi_2''$, $\zeta''_2$; furthermore, we choose $\xi_1$ and $\xi_2$ such that $\xi_2=\xi_1$. Note that $\xi_1''$ is in general slightly different from $\xi_2''$, but by a very small amount, since in the approximation (13) $\xi''$ is a slowly varying function of $\zeta$. We then check whether $\zeta''_1\zeta''_2<0$; if not, we choose another pair of values for $\zeta_1$ and $\zeta_2$, until the condition is verified. In practice, the segment parallel to the $\zeta$-axis needs to straddle one of the resonant return circles. At this point, we find the pre-image of a point with $\zeta''=0$, by using, for instance, regula falsi iterations; let us call the coordinates of this point $(\xi_0, \zeta_0)$, and the coordinates of its image $(\xi''_0, 0)$. If $\vert\xi''_0\vert<b_\oplus$, where $b_\oplus$ is the radius of the Earth $r_\oplus$ augmented by the gravitational focusing

$$b_\oplus=r_\oplus\; \sqrt{1+\frac{2c}{r_\oplus}}$$

then to a good degree of precision the pre-image of the point of coordinates $\xi''_0$, $\sqrt{b_\oplus^2-{\xi''_0}^2}$ is the point of coordinates $\xi_0$, $\sqrt{b_\oplus^2-{\xi''_0}^2}/(\partial\zeta''/\partial\zeta)$. 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 $\xi$) remains essentially unchanged between the two encounters, the "vertical'' one, along $\zeta$, 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 $\partial \zeta''/\partial\xi'$ and $\partial \zeta''/\partial\zeta'$.

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 XF₁₁ in 2028 and of 1999 AN₁₀ in 2027; let us start from the latter.

The encounters of a rather large number of "virtual'' 1999 AN₁₀'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 AN₁₀'s at that time occupy a very small region of $(a,e,i,\omega,\Omega)$ 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, $\theta=105\fdg3$, $\phi=41\fdg3$ - the same MOID of 0.000246 AU, so that $\xi=0.000246$, 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 $\zeta$-axis that extends to $\zeta=\pm0.21$.

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 $\Delta t$ (difference in time from closest approach; $\Delta t=0$ 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 $\Delta t$(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 $\Delta t=0$ in, from right to left, 2040, 2038, 2036, 2034, 2039 and 2032/2037 (the two leftmost lines, crossing on the $\Delta t=0$ 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 1999AN₁₀. 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  $b_{\rm max}$, at a specified date, two conditions must be met, since both $\xi''$ and $\zeta''$, the b-plane coordinates at the next encounter, must be small. The vicinity to the resonant circle takes care of $\zeta''$, and Eq. (13) gives the condition on $\xi$ at the previous encounter in order to have $\xi''$ small; the latter condition forces the keyholes not to be symmetric with respect to the $\zeta$-axis. Then, the keyholes are near the intersections of a strip, of width approximately equal to $b_{\rm max}$, parallel to the $\zeta$-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 $\xi$ and $\zeta$ 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 $\xi =5.8$, 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 $\xi$ 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 AN₁₀; 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 XF₁₁; for this asteroid, U=0.459, $\theta=84\fdg0$ and $\phi=99\fdg5$. 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 $\partial\zeta''/\partial\xi$ and $\partial\zeta''/\partial\zeta$ the smaller terms of first order in $\xi'$, $\zeta'$ (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 $\xi$, roughly the diameter of the Earth augmented by the gravitational focusing; it is in $\zeta$ 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 $\zeta$ 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 $\partial\zeta''/\partial\zeta$.

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 $\zeta$ 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 $\epsilon\neq0$); however, we have checked that the terms of the first order in $\epsilon$ 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 ( $\partial\zeta''/\partial\zeta$) 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. $\xi$) 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 $\zeta$ is close to the coordinate D of the center of the circle. When this is the case, we must distinguish between two possibilities, $R\la\xi$, and $R\ga\xi$; 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, $R\la\xi$, we have that the portion of the wire for which $D\approx\zeta$ 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 $\xi''$-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 $R\ga\xi$, we have interesting consequences for the portion of the wire for which $D\approx\zeta$. 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 XF₁₁ 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 $\zeta$ is much smaller than in the other cases.

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

Concerning $\xi$, the situation is essentially governed by the value of $\partial\zeta''/\partial\xi$ that, as is easily visible from Eq. (10), can be rather large, so that very small changes in $\xi$ produce large changes $\zeta$, 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

 
• Bonanno, C. 2000, A&A, 360, 411 In the text NASA ADS  
• Carusi, A., Kresák, L., & Valsecchi, G. B. 1981, A&A, 99, 262 In the text NASA ADS  
• Carusi, A., Valsecchi, G. B., & Greenberg, R. 1990, CeMDA, 49, 111 In the text NASA ADS  
• Chodas, P. W. 1999, BAAS, 31, 1117 In the text  
• Greenberg, R., Carusi, A., & Valsecchi, G. B. 1988, Icarus, 75, 1 In the text NASA ADS  
• Gronchi, G. F., & Milani, A. 2001, Icarus, 152, 58 In the text NASA ADS  
• Le Verrier, U.-J. 1844, Compt. Rend. Acad. Sci., 19, 982 In the text  
• Le Verrier, U.-J. 1848, Compt. Rend. Acad. Sci., 26, 465 In the text  
• Le Verrier, U.-J. 1857, Ann. Obs. Paris, 3, 203 In the text NASA ADS  
• Lexell, A. J. 1778, Acta Acad. Sci. Imper. Petropol., pars posterior, 12, published in 1781 In the text  
• Marsden, B. G. 1999, J. Brit. Interpl. Soc., 52, 195 In the text  
• Milani, A. 1999, Icarus, 137, 269 In the text NASA ADS  
• Milani, A., & Valsecchi, G. B. 1999, Icarus 140, 408 In the text NASA ADS  
• Milani, A., Chesley, S. R., & Valsecchi, G. B. 1999, A&A, 346, L65 In the text NASA ADS  
• Milani, A., La Spina, A., Sansaturio, M. E., & Chesley, S. R. 2000a, Icarus, 144, 39 In the text NASA ADS  
• Milani, A., Chesley, S. R., Boattini, A., & Valsecchi, G. B. 2000b, Icarus, 145, 12 In the text NASA ADS  
• Milani, A., Chesley, S. R., & Valsecchi, G. B. 2000c, Planet. Space Sci., 48, 945 In the text NASA ADS  
• Öpik, E. J. 1951, Proc. R. Irish Acad., 54, 165 In the text  
• Öpik, E. J. 1976, Interplanetary Encounters (New York: Elsevier) In the text  
• Valsecchi, G. B. 2002, in Singularities in gravitational systems. Applications to chaotic transport in the solar system, ed. Cl. Froeschlé, & D. Benest (Heidelberg, Germany: Springer Verlag), 145 In the text  
• Valsecchi, G. B., Manara, A. 1997, A&A, 323, 986 In the text NASA ADS  
• Valsecchi, G. B., Froeschlé, Cl., & Gonczi, R. 1997, Planet. Space Sci., 45, 1561 In the text NASA ADS  
• Valsecchi, G. B., Milani, A., Gronchi, G. F., & Chesley, S. R. 2000, CeMDA, 78, 83 In the text NASA ADS

  
Online Material

Appendix A: Computational details

  
A.1 Pre-encounter state vector

The pre-encounter state vector has components (U, $\theta$, $\phi$, $\xi$, $\zeta$, t₀). To compute these components from the heliocentric elements (a, e, i, $\omega$, $\Omega$, f) we have to consider where the encounter of interest takes place.

Let us call $\lambda_{\rm p}$ the longitude of the planet at time t₀; then, if at that time $f=-\omega$, the encounter takes place at the ascending node, in the post-perihelion branch of the orbit if $\pi<\omega<2\pi$, and in the pre-perihelion branch otherwise. If at time t₀ we have $f=\pi-\omega$, the encounter takes place at the descending node, in the post-perihelion branch of the orbit if $0<\omega<\pi$, 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 ${\cal O}(X_0^2+Y_0^2)$, at the ascending node we have

$$\xi = \cos\phi\left[\frac{a(1-e^2)}{1+e\cos\omega}-1\right]$$

$$\zeta = \xi\cos\theta\tan\phi-\sin\theta\left(1+\frac{\xi}{\cos\phi} \right)\tan(\Omega-\lambda_p),$$

with

$$0<\phi<\frac{\pi}{2} \quad {\rm if} \quad 0<f<\pi,~ \pi<\omega<2\pi$$

$$\frac{3\pi}{2}<\phi<2\pi {\rm if} \pi<f<2\pi, 0<\omega<\pi,$$

whereas at the descending node we have

$$\xi = \cos\phi\left[\frac{a(1-e^2)}{1-e\cos\omega}-1\right]$$

$$\zeta = \xi\cos\theta\tan\phi-\sin\theta\left(1+\frac{\xi}{\cos\phi} \right)\tan(\Omega-\lambda_p-\pi),$$

with

$$\frac{\pi}{2}<\phi<\pi {\rm if} 0<f<\pi, ~0<\omega<\pi$$

$$\pi<\phi<\frac{3\pi}{2} {\rm if} \pi<f<2\pi,~ \pi<\omega<2\pi.$$

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, $\theta$ and $\phi$:

$$U = \sqrt{3-\frac{1}{a}-2\sqrt{a\left(1-e^2\right)}\cos{i}}$$

$$\theta = \arccos\frac{\sqrt{a(1-e^2)}\cos{i}-1} {\sqrt{3-1/a-2\sqrt{a\left(1-e^2\right)}\cos{i}}}$$

$$\tan\phi = \frac{\pm\sqrt{2-1/a-a\left(1-e^2\right)}}{\pm\sqrt{a\left(1-e^2\right)}\sin{i}},$$

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

Since

$$\sqrt{a\left(1-e^2\right)}\cos{i}=\frac{3-1/a-U^2}{2},$$

a simple expression for $\theta$ is in terms of U, a

$$\cos\theta=\frac{1-U^2-1/a}{2U}\cdot$$

Finally, a, e, i can be obtained from U, $\theta$, $\phi$

$$\begin{eqnarray*}a & = & \frac{1}{1-U^2-2U\cos\theta} \\ e & = & U\sqrt{(U+2\co... ...\ i & = & \arctan\frac{U\sin\theta\cos\phi}{1+U\cos\theta}\cdot \end{eqnarray*}$$

The non-zero planetocentric distance of the small body affects the computation of U, $\theta$, $\phi$ 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 $r=1+\epsilon$, an explicit computation for $U_\epsilon$, $\theta_\epsilon$ and $\phi_\epsilon$, i.e. the values of U, $\theta$ and $\phi$ corresponding to such a situation, gives, at the first order in $\epsilon$,

$$U_\epsilon = U-\epsilon\cdot\frac{2U^3\cos\theta +4U^2\cos^2\theta+U^2+1}{U^2+2U\cos\theta+1}$$

$$\cos\theta_\epsilon = \cos\theta+\epsilon\cdot\frac{U^3(2\co... ...\cos\theta+1)+U\cos\theta+\cos\theta-1} {U^3+2U^2\cos\theta+U}$$

$$\sin\phi_\epsilon = \sin\phi-\epsilon\cdot\frac{U^3(1-\sin^2... ...\phi)+2\cos\theta} {U\sin^2\theta\sin\phi(U^2+2U\cos\theta+1)}$$

$$\cos\phi_\epsilon = \cos\phi\cdot\left(1+\epsilon\cdot\frac{U^2+2U\cos\theta-1} {U^2+2U\cos\theta+1}\right),$$

where U, $\theta$, $\phi$ are the values computed for $\epsilon=0$, and

$$\epsilon=\zeta\cos\theta\sin\phi+\xi\cos\phi.$$

On the other hand, if we have a generic vector $\vec{U}$ applied in $r=1+\epsilon$, we obtain $a_\epsilon$, $e_\epsilon$, $i_\epsilon$ from the components of $\vec{U}$ with the following expressions, again at first order in $\epsilon$:

$$\frac{1}{a_\epsilon} \approx \frac{1}{a}-2\epsilon$$

$$a_\epsilon \approx \frac{a}{1-2a\epsilon}$$

$$\approx a\left(1+2a\epsilon\right)~~~~{\rm for} \quad a\ll\frac{1}{\epsilon}$$

$$e_\epsilon = e +\epsilon\cdot U(2U-U\sin^2\theta\sin^2\phi+4\cos\theta)$$

$$\tan{i_\epsilon} = \tan{i} \cdot\left(1-\frac{\epsilon}{1+U_y}\right),$$

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

   
A.2 Post-encounter state vector

The close encounter can be seen as an operator $\vec{E}$ that maps the pre-encounter state vector $\vec{V}$, with components $(U, \theta, \phi, \xi, \zeta, t_0)$, into the post-encounter one $\vec{V'}$, with components $(U', \theta', \phi', \xi', \zeta', t'_0)$

$$\vec{V'}=\vec{E}\vec{V}.$$

The components of the post-encounter state vector, as functions of the pre-encounter state vector components (with $b=\sqrt{\xi^2+\zeta^2}$), are

U' = U

$$\theta' = \arccos\frac{\left(b^2-c^2\right)\cos\theta+2c\zeta... ...ht)\sin\theta-2c\zeta\cos\theta\right]^2 +4c^2\xi^2}}{b^2+c^2}$$

$$\phi' = \arccos\frac{\left[\left(b^2-c^2\right)\sin\theta-2c\... ...^2-c^2\right)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2}}$$

$$\xi' = \frac{\left(b^2+c^2\right)\xi\sin\theta} {\sqrt{\left[... ...^2-c^2\right)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2}}$$

$$\zeta' = \frac{\left(b^2-c^2\right)\zeta\sin\theta-2b^2c\cos\... ...^2-c^2\right)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2}}$$

$$t'_0 = t_0+\frac{2c\left[\xi\sin\phi(2\zeta\cos\theta-\xi\tan... ... -2c\zeta\cos\theta\right]\cos\phi+2c\xi\sin\phi\right\}}\cdot$$

A.3 Propagation to the next encounter

A purely Keplerian propagation can be seen as an operator $\vec{P}$ that maps the post-encounter state vector $\vec{V'}$ in the pre-next-encounter one $\vec{V''}$, with components (U'', $\theta''$, $\phi''$, $\xi''$, $\zeta''$, t''₀)

$$\vec{V''}=\vec{P}\vec{V'}.$$

The transformation is given by:

U'' = U'

$$\theta'' = \theta'$$

$$\phi'' = \phi'$$

$$\xi'' = \xi'$$

$$\zeta'' = \zeta' -\left({\rm mod}\left[h\cdot2\pi a'^{3/2}+\pi,2\pi\right]-\pi\right)\sin\theta'$$

$$t''_0 = t'_0+h\cdot2\pi a'^{3/2},$$

where the key parameter is the post-encounter semimajor axis a'. Note that the only components of $\vec{V'}$ that differ from those of $\vec{V''}$ are the $\zeta$-component on the post encounter b-plane $\zeta'$ and the time of nodal passage t'₀.

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

$$a' = \frac{b^2+c^2}{\left(b^2+c^2\right)\left(1-U^2\right) -2... ...[\left(b^2-c^2\right)\cos\theta+2c\zeta\sin\theta\right]}\cdot$$

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

$$\frac{1}{a'_{\epsilon'}} = \frac{1}{a'}-2\epsilon'$$

with

$$\epsilon' = \frac{2b^2c\sin\theta\sin\phi+\left(b^2-c^2\right... ...eft(\zeta\cos\theta\sin\phi+\xi\cos\phi\right)}{b^2+c^2}\cdot$$

   
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

$$\vec{V''}=\vec{P}\vec{E}\vec{V}.$$

The matrix of derivatives is then

$$\begin{eqnarray*}\frac{\partial(U'',\theta'',\phi'',\xi'',\zeta'',t''_0)} {\part... ...\phi',\xi',\zeta',t'_0)} {\partial(U,\theta,\phi,\xi,\zeta,t_0)} \end{eqnarray*}$$

where the propagation derivatives matrix is regarded as a function of the pre-first encounter variables (U, $\theta$, $\phi$, $\xi$, $\zeta$, t₀).

For encounters computed with the extended Öpik theory,

$$\frac{\partial(U',\theta',\phi',\xi',\zeta',t'_0)} {\partial(... ...& \frac{\partial t'_0}{\partial\zeta} & 1 \end{array}\right],$$

and, for a purely Keplerian heliocentric propagation between encounters,

$$\frac{\partial(U'',\theta'',\phi'',\xi'',\zeta'',t''_0)} {\pa... ... \frac{\partial t''_0}{\partial\zeta'} & 1 \end{array}\right],$$

so that the composite matrix has the following structure

$$\frac{\partial\left(U'',\theta'',\phi'',\xi'',\zeta'',t''_0\r... ... \frac{\partial t''_0}{\partial\zeta} & 1 \end{array}\right].$$

We are interested, in particular, in the submatrix ${\partial\left(\xi'',\zeta''\right)}/{\partial(\xi,\zeta)}$, 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., $\xi''=\xi'$, so that the first row of the matrix is

$$\begin{eqnarray*}\frac{\partial\xi''}{\partial\xi} & = & \frac{\partial\xi'}{\pa... ...zeta}\cos\phi' +X'_0\frac{\partial\cos\phi'}{\partial\zeta}\cdot \end{eqnarray*}$$

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

$$\begin{eqnarray*}\frac{\partial\zeta''}{\partial\xi} & = & \frac{\partial\zeta'... ...os\theta'}{\partial\zeta} +\frac{\partial\zeta'}{\partial\zeta}, \end{eqnarray*}$$

with the derivatives of $\zeta'$ with respect to $\xi$ and $\zeta$ given by

$$\frac{\partial\zeta'}{\partial\xi} = \frac{\partial X'_0}{\p... ...ial\xi}\sin\theta' +Y'_0\frac{\partial\sin\theta'}{\partial\xi}$$

$$\frac{\partial\zeta'}{\partial\zeta} = \frac{\partial X'_0}{... ...sin\theta' +Y'_0\frac{\partial\sin\theta'}{\partial\zeta}\cdot$$

Note that the second and the third term in the derivatives of the Keplerian propagation are of the first order in $\xi'$, $\zeta'$, 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

$$\frac{\partial X'_0}{\partial\xi} = \frac{\partial X'(t_b)}{... ...frac{Z'(t_b)}{U'^2_z} \frac{\partial U'_z}{\partial\xi}\right)$$

$$\frac{\partial X'_0}{\partial\zeta} = \frac{\partial X'(t_b)... ...ac{Z'(t_b)}{U'^2_z} \frac{\partial U'_z}{\partial\zeta}\right)$$

$$\frac{\partial Y'_0}{\partial\xi} = \frac{\partial Y'(t_b)}{... ...frac{Z'(t_b)}{U'^2_z} \frac{\partial U'_z}{\partial\xi}\right)$$

$$\frac{\partial Y'_0}{\partial\zeta} = \frac{\partial Y'(t_b)... ...ac{Z'(t_b)}{U'^2_z} \frac{\partial U'_z}{\partial\zeta}\right)$$

$$\frac{\partial U'_x}{\partial\xi} = U\left(\frac{\partial\si... ...\phi' +\sin\theta'\frac{\partial\sin\phi'}{\partial\xi}\right)$$

$$\frac{\partial U'_x}{\partial\zeta} = U\left(\frac{\partial\... ...hi' +\sin\theta'\frac{\partial\sin\phi'}{\partial\zeta}\right)$$

$$\frac{\partial U'_y}{\partial\xi} = U\frac{\partial\cos\theta'}{\partial\xi}$$

$$\frac{\partial U'_y}{\partial\zeta} = U\frac{\partial\cos\theta'}{\partial\zeta}$$

$$\frac{\partial U'_z}{\partial\xi} = U\left(\frac{\partial\si... ...\phi' +\sin\theta'\frac{\partial\cos\phi'}{\partial\xi}\right)$$

$$\frac{\partial U'_z}{\partial\zeta} = U\left(\frac{\partial\... ...i' +\sin\theta'\frac{\partial\cos\phi'}{\partial\zeta}\right)$$

$$\frac{\partial X'(t_b)}{\partial\xi} = \frac{4c^2\xi\sin\phi... ...\right)+\left(b^4-c^4\right)\cos\phi} {\left(b^2+c^2\right)^2}$$

$$\frac{\partial X'(t_b)}{\partial\zeta} = \frac{4c^2\zeta\lef... ...+4c^2\zeta^2\right)\cos\theta\sin\phi}{\left(b^2+c^2\right)^2}$$

$$\frac{\partial Y'(t_b)}{\partial\xi} = \frac{4c^2\xi\left(c\cos\theta-\zeta\sin\theta\right)}{\left(b^2+c^2\right)^2}$$

$$\frac{\partial Y'(t_b)}{\partial\zeta} = \frac{4c^2\zeta\lef... ...ight)+\left(c^4-b^4\right)\sin\theta} {\left(b^2+c^2\right)^2}$$

$$\frac{\partial Z'(t_b)}{\partial\xi} = \frac{4c^2\xi\cos\phi... ...eft(b^4-c^4+2c^2\xi^2\right)\sin\phi} {\left(b^2+c^2\right)^2}$$

$$\frac{\partial Z'(t_b)}{\partial\zeta} = \frac{4c^2\zeta\lef... ...4c^2\zeta^2\right)\cos\theta\cos\phi}{\left(b^2+c^2\right)^2}$$

$$\frac{\partial t_b}{\partial\xi} = -\frac{\tan\phi}{U\sin\theta}$$

$$\frac{\partial t_b}{\partial\zeta} = \frac{1}{U\tan\theta}$$

$$\frac{\partial t'_0}{\partial\xi} = \frac{\partial t_b}{\par... ...al\xi}+\frac{Z'(t_b)}{U'^2_z}\frac{\partial U'_z}{\partial\xi}$$

$$\frac{\partial t'_0}{\partial\zeta} = \frac{\partial t_b}{\p... ...eta}+\frac{Z'(t_b)}{U'^2_z}\frac{\partial U'_z}{\partial\zeta}$$

$$\frac{\partial a'}{\partial\xi} = \frac{8Uc\xi\left(c\cos\th... ...t(b^2-c^2\right)\cos\theta+2c\zeta\sin\theta\right]\right\}^2}$$

$$\frac{\partial a'}{\partial\zeta} = \frac{8Uc\zeta(c\cos\the... ...t(b^2-c^2\right)\cos\theta+2c\zeta\sin\theta\right]\right\}^2}$$

$$\frac{\partial\epsilon}{\partial\xi} = \cos\phi$$

$$\frac{\partial\epsilon}{\partial\zeta} = \cos\theta\sin\phi$$

$$\frac{\partial\epsilon'}{\partial\xi} = \frac{4c^2\xi\left(c... ...n\right)+\left(b^4-c^4\right)\cos\phi}{\left(b^2+c^2\right)^2}$$

$$\frac{\partial\epsilon'}{\partial\zeta} = \frac{4c^2\zeta\le... ...ft(b^4-c^4\right)\cos\theta\sin\phi} {\left(b^2+c^2\right)^2}$$

$$\frac{\partial a'_{\epsilon'}}{\partial\xi} = \frac{\partial... ...{\partial\xi}+4\epsilon'a'\cdot\frac{\partial a'}{\partial\xi}$$

$$\frac{\partial a'_{\epsilon'}}{\partial\zeta} = \frac{\parti... ...rtial\zeta}+4\epsilon'a'\cdot\frac{\partial a'}{\partial\zeta}$$

$$\frac{\partial t''_0}{\partial\xi} = \frac{\partial t'_0}{\p... ...qrt{a'_{\epsilon'}}\frac{\partial a'_{\epsilon'}}{\partial\xi}$$

$$\frac{\partial t''_0}{\partial\zeta} = \frac{\partial t'_0}{... ...{a'_{\epsilon'}}\frac{\partial a'_{\epsilon'}}{\partial\zeta}$$

$$\frac{\partial\cos\theta'}{\partial\xi} = \frac{4c\xi\left(c\cos\theta-\zeta\sin\theta\right)}{\left(b^2+c^2\right)^2}$$

$$\frac{\partial\cos\theta'}{\partial\zeta} = \frac{2c\left[\l... ...\theta-\zeta\sin\theta\right)\right]} {\left(b^2+c^2\right)^2}$$

$$\frac{\partial\sin\theta'}{\partial\xi} = -\frac{\left(b^2-c... ...ht]^2+4c^2\xi^2}} \cdot\frac{\partial\cos\theta'}{\partial\xi}$$

$$\frac{\partial\sin\theta'}{\partial\zeta} = -\frac{\left(b^2... ...]^2+4c^2\xi^2}} \cdot\frac{\partial\cos\theta'}{\partial\zeta}$$

$$\begin{eqnarray*}\frac{\partial\cos\phi'}{\partial\xi} & =& \frac{2\left(\xi\si... ...ight)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2\right\}^3}} \end{eqnarray*}$$

$$\begin{eqnarray*}\frac{\partial\cos\phi'}{\partial\zeta} & = & \frac{2\cos\phi\... ...ight)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2\right\}^3}} \end{eqnarray*}$$

$$\begin{eqnarray*}\frac{\partial\sin\phi'}{\partial\xi} & = & \frac{2\left(\xi\s... ...ight)\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2\right\}^3}} \end{eqnarray*}$$

$$\begin{eqnarray*}\frac{\partial\sin\phi'}{\partial\zeta} &= & \frac{2\sin\phi\l... ...\sin\theta-2c\zeta\cos\theta\right]^2+4c^2\xi^2\right\}^3}}\cdot \end{eqnarray*}$$

   
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/U² is very small. This happens because for these planets $m\la3\times 10^{-6}$, and typical values of U are around 0.5. On the other hand, for an Earth-encountering asteroid with $U\approx0.5$, a rather typical value, the minimum value of b that avoids collision is

$$b=\sqrt{r_\oplus^2+2cr_\oplus}\approx5\times 10^{-5},$$

so that for most encounters the condition $c^2\ll b^2$ is verified (the only exceptions are extremely deep encounters with very low U.)

If we can assume $c^2\ll b^2$, most expressions of this paper can be approximated by much simpler formulas neglecting $\mathcal{O}(c^2/b^2)$terms. We have collected these approximate formulas below.

The post-encounter state vector:

U' = U

$$\theta' \approx \arccos\left(\cos\theta+\frac{2c\zeta}{b^2}\... ...ox\arcsin\left(\sin\theta-\frac{2c\zeta}{b^2}\cos\theta\right)$$

$$\phi' \approx \arccos\left(\cos\phi+\frac{2c\xi}{b^2\sin\the... ...rcsin\left(\sin\phi-\frac{2c\xi}{b^2\sin\theta}\cos\phi\right)$$

$$\xi' \approx \xi\left(1+\frac{2c\zeta}{b^2}\frac{\cos\theta}{\sin\theta}\right)$$

$$\zeta' \approx \zeta\left(1+\frac{2c\zeta}{b^2}\frac{\cos\theta}{\sin\theta}\right) -\frac{2c}{\tan\theta}$$

$$t'_0 \approx t_0 +\frac{2c\left[\xi\sin\phi\left(2\zeta\cos\t... ...^2\theta+\zeta^2\right)\right]}{Ub^2\sin^2\theta\cos\phi}\cdot$$

Note that, in this approximation, unless $\sin\theta$ is very small, $\xi'/\xi=1+\mathcal{O}(c/b)$, 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

$$\epsilon'\approx\epsilon+2c\sin\theta\sin\phi.$$

Post-encounter semimajor axis a':

$$a'\approx a\left(1+\frac{4aUc\zeta}{b^2}\sin\theta\right).$$

Partial derivatives:

$$\frac{\partial X'(t_b)}{\partial\xi} \approx \cos\phi$$

$$\frac{\partial X'(t_b)}{\partial\zeta} \approx \cos\theta\sin\phi$$

$$\frac{\partial Y'(t_b)}{\partial\xi} \approx 0$$

$$\frac{\partial Y'(t_b)}{\partial\zeta} \approx -\sin\theta$$

$$\frac{\partial Z'(t_b)}{\partial\xi} \approx -\sin\phi$$

$$\frac{\partial Z'(t_b)}{\partial\zeta} \approx \cos\theta\cos\phi$$

$$\frac{\partial a'}{\partial\xi} \approx -\frac{8Ua^2c\xi\zeta}{b^4}\sin\theta$$

$$\frac{\partial a'}{\partial\zeta} \approx \frac{4Ua^2c\left(\xi^2-\zeta^2\right)}{b^4}\sin\theta$$

$$\frac{\partial\cos\theta'}{\partial\xi} \approx -\frac{4c\xi\zeta}{b^4}\sin\theta$$

$$\frac{\partial\cos\theta'}{\partial\zeta} \approx \frac{2c\left(\xi^2-\zeta^2\right)}{b^4}\sin\theta$$

$$\frac{\partial\sin\theta'}{\partial\xi} \approx \frac{4c\xi\zeta}{b^4}\cos\theta$$

$$\frac{\partial\sin\theta'}{\partial\zeta} \approx \frac{2c\left(\zeta^2-\xi^2\right)}{b^4}\cos\theta$$

$$\frac{\partial\cos\phi'}{\partial\xi} \approx \frac{2c\left(\zeta^2-\xi^2\right)}{b^4}\frac{\sin\phi}{\sin\theta}$$

$$\frac{\partial\cos\phi'}{\partial\zeta} \approx -\frac{4c\xi\zeta}{b^4}\frac{\sin\phi}{\sin\theta}$$

$$\frac{\partial\sin\phi'}{\partial\xi} \approx \frac{2c\left(\xi^2-\zeta^2\right)}{b^4}\frac{\cos\phi}{\sin\theta}$$

$$\frac{\partial\sin\phi'}{\partial\zeta} \approx \frac{4c\xi\zeta}{b^4}\frac{\cos\phi}{\sin\theta}$$

$$\frac{\partial X'_0}{\partial\xi} \approx \frac{1}{\cos\phi}... ...ta\cos\phi-2\xi\zeta\sin\phi\right]} {b^4\sin\theta\cos^2\phi}$$

$$\frac{\partial X'_0}{\partial\zeta} \approx \frac{2c\xi\left... ...ta\cos\phi+2\xi\zeta\sin\phi\right]} {b^4\sin\theta\cos^2\phi}$$

$$\frac{\partial Y'_0}{\partial\xi} \approx \frac{\cos\theta\s... ...\sin^2\phi\right)\cos\theta\right]}{b^4\sin^2\theta\cos^2\phi}$$

$$\frac{\partial Y'_0}{\partial\zeta} \approx -\frac{1}{\sin\t... ...2\phi\right)\cos\theta\right]}{b^4\sin^2\theta\cos^2\phi}\cdot$$

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 $c\approx0.01$; considering that b could have been as low as 0.0003 (Le Verrier 1857), we have that for the closest approach $b^2\ll c^2$, a case in which the following simplified expressions hold.

Post-encounter state vector:

U' = U

$$\theta' \approx \arccos\left(-\cos\theta+\frac{2\zeta}{c}\si... ...pprox\arcsin\left(\sin\theta+\frac{2\zeta}{c}\cos\theta\right)$$

$$\phi' \approx \arccos\left(-\cos\phi+\frac{2\xi}{c}\frac{\si... ...ft(-\sin\phi -\frac{2\xi}{c}\frac{\cos\phi}{\sin\theta}\right)$$

$$\xi' \approx \xi\left(1-\frac{2\zeta}{c}\frac{\cos\theta}{\sin\theta}\right)$$

$$\zeta' \approx -\zeta\left(1-\frac{2\zeta}{c}\frac{\cos\theta}{\sin\theta}\right)$$

$$t'_0 \approx t_0 -\frac{2\left[\xi(2\zeta\cos\theta-\xi\tan\p... ...^2\theta+\zeta^2)\cos\phi\right]}{Uc\sin^2\theta\cos\phi}\cdot$$

Note that, also in this approximation, unless $\sin\theta$ is very small, $\xi'/\xi=1+\mathcal{O}(b/c)$, that is the local MOID is nearly conserved. Note also that

$$\epsilon'\approx-\epsilon.$$

Post-encounter semimajor axis a':

$$a'\approx a\left(1+\frac{4aU\zeta}{c}\sin\theta\right).$$

Partial derivatives:

$$\frac{\partial X'(t_b)}{\partial\xi} \approx -\cos\phi+\frac{4\xi}{c}\sin\theta\sin\phi$$

$$\frac{\partial X'(t_b)}{\partial\zeta} \approx -\cos\theta\sin\phi+\frac{4\zeta}{c}\sin\theta\sin\phi$$

$$\frac{\partial Y'(t_b)}{\partial\xi} \approx \frac{4\xi}{c}\cos\theta$$

$$\frac{\partial Y'(t_b)}{\partial\zeta} \approx \sin\theta+\frac{4\zeta}{c}\cos\theta$$

$$\frac{\partial Z'(t_b)}{\partial\xi} \approx \sin\phi+\frac{4\xi}{c}\sin\theta\cos\phi$$

$$\frac{\partial Z'(t_b)}{\partial\zeta} \approx -\cos\theta\cos\phi+\frac{4\zeta}{c}\sin\theta\cos\phi$$

$$\frac{\partial a'}{\partial\xi} \approx \frac{8U\xi}{c^2}\frac{\cos\theta}{\left(1-U^2+2U\cos\theta\right)^2}$$

$$\frac{\partial a'}{\partial\zeta} \approx \frac{4U}{c^2}\fra... ...a}{c^2}\frac{\sin^2\theta}{\left(1-U^2+2U\cos\theta\right)^3}$$

$$\frac{\partial\cos\theta'}{\partial\xi} \approx \frac{4\xi}{c^2}\cos\theta$$

$$\frac{\partial\cos\theta'}{\partial\zeta} \approx \frac{2\left(c\sin\theta+2\zeta\cos\theta\right)}{c^2}$$

$$\frac{\partial\sin\theta'}{\partial\xi} \approx \frac{4\xi}{c^2}\frac{\cos^2\theta}{\sin\theta}$$

$$\frac{\partial\sin\theta'}{\partial\zeta} \approx \frac{2\left(c\cos\theta-2\zeta\sin\theta\right)}{c^2}$$

$$\frac{\partial\cos\phi'}{\partial\xi} \approx \frac{2\left(c... ...2\xi\cos\phi-2\zeta\cos\theta\sin\phi\right)} {c^2\sin^2\theta}$$

$$\frac{\partial\cos\phi'}{\partial\zeta} \approx -\frac{4\xi\cos\theta\sin\phi}{c^2\sin^2\theta}$$

$$\frac{\partial\sin\phi'}{\partial\xi} \approx -\frac{2\left(... ...2\xi\sin\phi-2\zeta\cos\theta\cos\phi\right)} {c^2\sin^2\theta}$$

$$\frac{\partial\sin\phi'}{\partial\zeta} \approx \frac{4\xi\cos\theta\cos\phi}{c^2\sin^2\theta}$$

$$\frac{\partial X'_0}{\partial\xi} \approx -\frac{1}{\cos\phi... ...a\cos\theta\cos\phi-2\xi\sin\phi\right)}{c\sin\theta\cos^2\phi}$$

$$\frac{\partial X'_0}{\partial\zeta} \approx \frac{2\xi}{c}\frac{\cos\theta}{\sin\theta\cos\phi}$$

$$\frac{\partial Y'_0}{\partial\xi} \approx -\frac{\cos\theta\... ...phi\right)\cos\theta\sin^2\phi\right]}{c\sin^2\theta\cos^2\phi}$$

$$\frac{\partial Y'_0}{\partial\zeta} \approx \frac{1}{\sin\th... ...a\cos\phi\right)\sin\phi\right]}{c\sin^2\theta\cos^2\phi}\cdot$$