The evolution of the orbit distance in the double averaged restricted 3-body problem with crossing singularities

We study the long term evolution of the distance between two Keplerian confocal trajectories in the framework of the averaged restricted 3-body problem. The bodies may represent the Sun, a solar system planet and an asteroid. The secular evolution of the orbital elements of the asteroid is computed by averaging the equations of motion over the mean anomalies of the asteroid and the planet. When an orbit crossing with the planet occurs the averaged equations become singular. However, it is possible to define piecewise differentiable solutions by extending the averaged vector field beyond the singularity from both sides of the orbit crossing set. In this paper we improve the previous results, concerning in particular the singularity extraction technique, and show that the extended vector fields are Lipschitz-continuous. Moreover, we consider the distance between the Keplerian trajectories of the small body and of the planet. Apart from exceptional cases, we can select a sign for this distance so that it becomes an analytic map of the orbital elements near to crossing configurations. We prove that the evolution of the 'signed' distance along the averaged vector field is more regular than that of the elements in a neighborhood of crossing times. A comparison between averaged and non-averaged evolutions and an application of these results are shown using orbits of near-Earth asteroids.


Introduction
The distance between the trajectories of an asteroid (orbiting around the Sun) and our planet gives a first indication in the search for possible Earth impactors. We call it orbit distance and denote it by d min . 1 A necessary condition to have a very close approach or an impact with the Earth is that d min is small. Provided close approaches with the planets are avoided, the perturbations caused by the Earth make the asteroid trajectory change slowly with time. Moreover, the perturbations of the other planets produce small changes in both trajectories. The value of the semimajor axis of both is kept constant up to the first order in the small parameters (the ratio of the mass of each perturbing planet to the mass of the Sun). All these effects are responsible of a variation of d min . We can study the evolution of the asteroid in the framework of the restricted 3-body problem: Sun, planet, asteroid. Then it is easy to include more than one perturbing planet in the model, in fact the potential energy can be written as sum of terms each depending on one planet only.
If the asteroid has a close encounter with some planet, the perturbation of the latter generically produces a change in the semimajor axis of the asteroid. This can be estimated, and depends on the mass of the planet, the unperturbed planetocentric velocity of the small body and the impact parameter, see [18].
The orbits of near-Earth asteroids (NEAs, i.e. with perihelion distance ≤ 1.3 au) 2 are chaotic, with short Lyapounov times (see [19]), at most a few decades. After that period has elapsed, an orbit computed by numerical integration and the true orbit of the asteroid are practically unrelated and we can not make reliable predictions on the position of the asteroid. For this reason the averaging principle is applied to the equations of motion: it gives the average of the possible evolutions, which is useful in a statistical sense. However, the dynamical evolution often forces the trajectory of a NEA to cross that of the Earth. This produces a singularity in the averaged equations, where we take into account every possible position on the trajectories, including the collision configurations.
The problem of averaging on planet crossing orbits has been studied in [8] for planets on circular coplanar orbits and then generalized in [5] including nonzero eccentricities and inclinations of the planets. The work in [8] has been used to define proper elements for NEAs, that are integrals of an approximated problem, see [9]. In this paper we compute the main singular term by developing the distance between two points, one on the orbit of the Earth and the other on that of the asteroid, at its minimum points. This choice improves the results in [8], [5], where a development at the mutual nodes was used, because it avoids the artificial singularity occurring for vanishing mutual inclination of the two orbits. Moreover, we show that the averaged vector field admits two Lipschitz-continuous extensions from both sides of the orbit crossing set (see Theorem 4.2), which is useful for the numerical computation of the solutions.
The orbit distance d min is a singular function of the (osculating) orbital elements when the trajectories of the Earth and the asteroid intersect. However, by suitably choosing a sign for d min we obtain a map, denoted byd min , which is analytic in a neighborhood of most crossing configurations (see [11]).
Here we prove that, near to crossing configurations, the averaged evolution ofd min is more regular than the averaged evolution of the elements, which are piecewise differentiable functions of time.
The paper is organized as follows. Section 2 contains some preliminary results on the orbit distance. In Sections 3, 4, 5 we introduce the averaged equations, present the results on the singularity extraction method and give the definition of the generalized solutions, which go beyond crossing singularities. In Section 6 we prove the regularity of the secular evolution of the orbit distance. Section 7 is devoted to numerical experiments: we describe the algorithm for the computation of the generalized solutions and compare the averaged evolution with the solutions of the full equations of motion. We also show how this theory can be applied to estimate Earth crossing times for NEAs.

The orbit distance
Let (E j , v j ), j = 1, 2 be two sets of orbital elements of two celestial bodies on confocal Keplerian orbits. Here E j describes the trajectory of the orbit and v j is a parameter along the trajectory, e.g. the true anomaly. We denote by E = (E 1 , E 2 ) the two-orbit configuration, moreover we set V = (v 1 , v 2 ). In this paper we consider bounded trajectories only.
Choose a reference frame, with origin in the common focus, and write X j = X j (E j , v j ), j = 1, 2 for the Cartesian coordinates of the two bodies.
For a given two-orbit configuration E, we introduce the Keplerian distance function d, defined by where T 2 is the two-dimensional torus and | · | is the Euclidean norm.
The local minimum points of d can be found by computing all the critical points of d 2 . For this purpose in [6], [7], [12], [2] the authors have used methods of computational algebra, such as resultants and Gröbner's bases, which allow us to compute efficiently all the solutions.
Apart from the case of two concentric coplanar circles, or two overlapping ellipses, the function d 2 has finitely many stationary points. There exist configurations attaining 4 local minima of d 2 : this is thought to be the maximum possible, but a proof is not known yet. A simple computation shows that, for non-overlapping trajectories, the number of crossing points is at most two, see [7].
For each choice of the two-orbit configuration E, d min (E) gives the orbit distance.
The maps d h and d min are singular at crossing configurations, and their derivatives do not exist. We can deal with this singularity and obtain analytic maps in a neighborhood of a crossing configuration E c by properly choosing a sign for these maps. We note that d h , d min can present singularities without orbit crossings. The maps d h can have bifurcation singularities, so that the number of minimum points of d may change. Therefore the maps d h , d min are defined only locally. We say that a configuration E is non-degenerate if all the critical points of the Keplerian distance function are non-degenerate. If E is non-degenerate, there exists a neighborhood W of E ∈ R 10 such that the maps d h , restricted to W, do not have bifurcations. On the other hand, the map d min can lose regularity when two local minima exchange their role as absolute minimum. There are no additional singularities apart from those mentioned above. The behavior of the maps d h , d min has been investigated in [11]. However, a detailed analysis of the occurrence of bifurcations of stationary points and exchange of minima is still lacking.
Here we summarize the procedure to deal with the crossing singularity of d h ; the procedure for d min is the same. We consider the points on the two orbits corresponding to the local minimum points We introduce the vectors tangent to the trajectories E 1 , E 2 at these points and their cross product τ We also define and ∆ h are parallel, see Figure 1. Denoting 00 00 00 11 11 11 Earth orbit Asteroid orbit is an analytic function in a neighborhood of most crossing configurations. Indeed, this smoothing procedure fails at crossing configurations such that τ are parallel. A detailed proof can be found in [11]. Note that, to obtain regularity in a neighborhood of a crossing configuration, we lose continuity at the configurations with τ We shall call (signed) orbit distance the mapd min .

Averaged equations
Let us consider a restricted 3-body problem with the Sun, the Earth and an asteroid. The motion of the 2-body system Sun-Earth is a known function of time. We denote by X , X ′ ∈ R 3 the heliocentric position of the asteroid and the planet respectively. The equations of motion for the asteroid arë where k is Gauss' constant and µ is a small parameter representing the ratio of the Earth mass to the mass of the Sun. We study the motion using Delaunay's elements Y = (L, G, Z, ℓ, g, z), defined by where (a, e, I, ω, Ω, ℓ) are Keplerian elements, n is the mean motion and t 0 is the time of passage at perihelion. Delaunay's elements of the Earth are denoted by (L ′ , G ′ , Z ′ , ℓ ′ , g ′ , z ′ ). We write E = (E, E ′ ) for the two-orbit configuration, where E, E ′ are Delaunay's elements of the asteroid and the Earth respectively. Using the canonical variables Y, equations (3) can be written in Hamiltonian form asẎ where we use for the symplectic identity matrix of order 2n. The Hamiltonian is the difference of the two-body (asteroid, Sun) Hamiltonian and the perturbing function The function R is the sum of two terms: the first is the direct part of the perturbation, due to the attraction of the Earth and singular at collisions with it. The second is called indirect perturbation, and is due to the attraction of the Sun on the Earth.
We can reduce the number of degrees of freedom of (4) by averaging over the fast angular variables ℓ, ℓ ′ , which are the mean anomalies of the asteroid and the Earth. As a consequence, ℓ becomes a cyclic variable, so that the semimajor axis a is constant in this simplified dynamics. For a full account on averaging methods in Celestial Mechanics see [1].
The averaged equations of motion for the asteroid are given bẏ where Y = (G, Z, g, z) t , Y = (G, Z, g, z) t are some of Delaunay's elements, and is the vector of the averaged partial derivatives of the perturbing function R. Equation (5) corresponds to the scalar equationṡ We can easily include more planets in the model. In this case the perturbing function is sum of terms R i , each depending on the coordinates of the asteroid and the planet i only, with a small parameter µ i , representing the ratio of the mass of planet i to the mass of the Sun.
Note that, if there are mean motion resonances of low order with the planets, then the solutions of the averaged equations (5) may be not representative of the behavior of the corresponding components in the solutions of (4).
Moreover, when the planets are assumed to move on circular coplanar orbits we obtain an integrable problem. In fact the semimajor axis a, the component Z of the angular momentum orthogonal to the invariable plane 3 and the averaged Hamiltonian H are first integrals generically independent and in involution (i.e. with vanishing Poisson's brackets). Taking into account the eccentricity and the inclination of the planets the problem is not integrable any more.
In [14] the secular evolution of high eccentricity and inclination asteroids is studied in a restricted 3-body problem, with Jupiter on a circular orbit. Nevertheless, crossings with the perturbing planet are excluded in that work. In [15] there is a similar secular theory for a satellite of the Earth. The dynamical behavior described in [14], [15] is usually called Lidov-Kozai mechanism in the literature and an explicit solution to the related equations is given in [13].
If no orbit crossing occurs, by the theorem of differentiation under the integral sign the averaged equations of motion (5) are equal to Hamilton's equationsẎ where is the averaged perturbing function. The average of the indirect term of R is zero. When the orbit of the asteroid crosses that of the Earth a singularity appears in (5), corresponding to the existence of a collision for particular values of the mean anomalies. We study this singularity to define generalized solutions of (5) which go beyond planet crossings. Since the semimajor axis of the asteroid is constant in the averaged dynamics, we expect that the generalized solutions can be reliable only if there are no close approaches with the planet in the dynamics of equations (4).

Extraction of the singularity
In the following we denote by E c a non-degenerate crossing configuration with only one crossing point, and we choose the minimum point index h where is Taylor's remainder in the integral form. 4 We introduce the approximated distance where for a vector V = (v1, v2) and a smooth function V → f (V ).

Remark 1.
If the matrix A h is non-degenerate, then it is positive definite because V h is a minimum point, and this property holds in a suitably chosen neighborhood W of E c . At a crossing configuration E = E c the matrix A h is degenerate if and only if the tangent vectors τ h , τ ′ h are parallel (see [11]): First we estimate the remainder function 1/d − 1/δ h . To this aim we need the following: holds for E in W and for every V ∈ T 2 .
Proof. From (8), (9) we obtain the existence of a neighborhood U of (E c , V h (E c )) and a constant C 5 > 0 such that We select U so that no bifurcations of stationary points of d 2 occur and there exists a constant Relation (13) together with (7),(10) yield (11) for some C 1 , C 2 > 0. Moreover, we can find a neighborhood W of E c such that there are no bifurcations of stationary points of d 2 , and the inequalities (12) hold for some Proof. By Lemma 4.1 we can choose two neighborhoods W, V of E c and V h (E c ) respectively such that both (11) and (12) In this case we restrict W to a smaller set (the inner circle), as explained in the proof of Proposition 1. Figure 2). In U \ U Σ we have for a constant C > 0. Using the boundedness of we conclude the proof.
where Σ = {E ∈ W : d h (E) = 0}, can be extended continuously to the whole set W.
Proof. In the following we denote by C j , j = 8 . . . 14 some positive constants. We write ∂ ∂y k and give an estimate for the two terms at the right hand side. We choose a neighborhood U = W × V of (E c , V h (E c )) as in Proposition 1 so that, using (11), (12) and the boundedness of the remainder function, we have and the derivatives are uniformly bounded for E ∈ W since bifurcations do not occur. Hence the relation holds in U 0 , with C 11 = C 8 C 10 . We also have for (E, V ) ∈ U 0 , in fact for each α = (α 1 , α 2 ) with |α| = 3. Using again (11), (12) we obtain From (17), (20) we obtain (14) and the assert of the proposition follows using the boundedness of ∂ ∂y k 1/d , ∂ ∂y k 1/δ h in W × (T 2 \ V).
From (14) in Proposition 2 the average over T 2 of the derivatives of 1/d− 1/δ h in (15) is finite for each E in W, and can be computed by exchanging the integral and differential operators: therefore the average of the remainder function is continuously differentiable in W.
On the other hand, the average over T 2 of the derivatives with respect to y k of 1/δ h are non-convergent integrals for E ∈ Σ: for this reason the averaged vector field in (5) is not defined at orbit crossings. Next we show, exchanging again the integral and differential operators, that the average of these derivatives admit two analytic extensions to the whole W from both sides of the singular set Σ. For this purpose, given a neighborhood W of E c , we set withd h given by (1).
There exists a neighborhood W of E c such that the maps where y k is a component of Delaunay's elements Y , can be extended to two different analytic maps G + h,k , G − h,k , defined in W.
Proof. We choose W as in Proposition 2 and, if necessary, we restrict this neighborhood by requiring that τ = 0 in W, so thatd h | W is analytic. To investigate the behavior close to the singularity, for each E ∈ W, we can restrict the integrals to the domain for a suitable r > 0. By using the coordinate change ξ = A 1/2 h (V − V h ) and then polar coordinates (ρ, θ), defined by (ρ cos θ, ρ sin θ) = ξ, we have (1), and define on W the two analytic maps We observe that G + h,k (resp. G − h,k ) corresponds to the derivative of T 2 1/δ h dℓ dℓ ′ with respect to y k on W + (resp. W − ). Now we state the main result.
Theorem 4.2. The averages over T 2 of the derivatives of R with respect to Delaunay's elements y k can be extended to two Lipschitz-continuous maps ∂R ∂y k ± h on a neighborhood W of E c . These maps, restricted to W + , W − respectively, correspond to ∂R ∂y k . Moreover the following relations hold: with the derivatives ofd h given by (2).
Proof. Using the results of Propositions 2, 3 we define the extended maps by with G ± h,k given by (22). We show that the maps E → ∂R ∂y k ± h (E) are Lipschitz-continuous extensions to W of ∂R ∂y k . The maps G ± h,k are analytic in W, thus we only have to study the integrals T 2 ∂ ∂y k (1/d − 1/δ h ) dℓdℓ ′ . From Proposition 2 we know that these maps are continuous.
We only need to investigate the behavior close to the singularity, therefore we restrict these integrals to the domain D introduced in (21). We prove that the maps are bounded in W \ Σ, otherwise we could not find the analytic extensions G + h,k , G − h,k introduced in Proposition 3. 5 Then we show that the maps are bounded in W \ Σ. Using (7), (10) we write the integrand function in the right hand side of (25) as the sum of and of terms of the following kind: The integrals over D of the terms in (26) are not bounded in W \ Σ, but their sum is bounded and corresponds to (24). In the following we denote by C j , j = 15 . . . 34 some positive constants. Moreover we use the relation 3 and the developments 1 5 Actually we can prove that j,k are bounded in W \ Σ, and is unbounded but cancels out in the difference.
To estimate the integrals of the terms in (29) we make the following decomposition: p For the first term in (29) we obtain where we have used polar coordinates and the inequalities (35). The two integrals at the right hand side of (38) vanish: in fact using Fubini-Tonelli's theorem and passing to polar coordinates (ρ, θ), by relations (37) we obtain The computation for the other integral is analogous. The second term in (29) is estimated in a similar way: and the integral at the right hand side vanishes as well.
We conclude the proof observing that, using (22) and the theorem of differentiation under the integral sign, the derivatives ∂R ∂y k + h , ∂R ∂y k − h , restricted to W + , W − respectively, correspond to ∂R ∂y k , and their difference in W is given by (23).

Generalized solutions
We show that generically we can uniquely extend the solutions of (5) beyond the crossing singularity d min = 0. This is obtained by patching together classical solutions defined in the domains W + with solutions defined in W − , or vice versa.
Let a > 0 be a value for the semimajor axis of the asteroid and Y : I → R 4 be a continuous function defined in an open interval I ⊂ R, representing a possible evolution of the asteroid orbital elements Y = (G, Z, g, z). We introduce with where k is Gauss' constant and E ′ is a known function of time representing the evolution of the Earth. 7 Let T (Y ) be the set of times t c ∈ I such that d min (E(t c )) = 0, and assume that each t c is isolated, so that we can represent the set In the case of one perturbing planet E ′ (t) is constant and represents the trajectory of a solution of the 2-body problem. If we consider more than one perturbing planet then E ′ (t) changes with time due to the planetary perturbations.
as disjoint union of open intervals I j , with N a countable (possibly finite) set.
Definition 5.1. We say that Y is a generalized solution of (5) if its restriction to each I j , j ∈ N is a classical solution of (5) and, for each t c ∈ T (Y ), there exist finite values of Choose Y 0 ∈ R 4 and a time t 0 such that d min . We show how we can construct a generalized solution of the Cauchy probleṁ Let Y (t) be the maximal classical solution of (41), defined in the maximal interval J. Assume that t c = sup J < +∞, and lim t→t − c E(t) = E c , with E c a non-degenerate crossing configuration such that d min (E c ) = d h (E c ) = 0 for some h. Let W, W ± be chosen as in Theorem 4.2. Suppose that there exists τ ∈ (t 0 , t c ) such that E(t) ∈ W + for t ∈ (τ, t c ). Let Y τ = Y (τ ). By Theorem 4.2 there existsẎ c ∈ R 4 such that In fact relation (42) is fulfilled by the solution of the Cauchy problem 8 which corresponds to the solution of (41) in the interval (τ, t c ) and is defined also at the crossing time t c . Let us denote by Y c its value for t = t c . Using again Theorem 4.2 we can extend Y (t) beyond the crossing singularity by considering the new probleṁ The solution of (44) fulfils The vector field in (44) corresponds to −J 2 ∇ Y R on W − , thus we can continue the solution outside W and this procedure can be repeated at almost every crossing singularities. Indeed, the generalized solution is unique provided the evolution t → E(t) is not tangent to the orbit crossing set Σ. Moreover, if det A h = 0 the extraction of the singularity, described in Section 4, cannot be performed.
In case E(t) ∈ W − for t ∈ (τ, t c ) the previous discussion still holds if we exchange ( In this case (45) becomes ) .

Evolution of the orbit distance
We prove that the secular evolution ofd min is more regular than that of the orbital elements in a neighborhood of a planet crossing time. We introduce the secular evolution of the distancesd h and of the orbit distanced min : Assume these maps are defined in an open interval containing a crossing time t c , and suppose E c = E(t c ) is a non-degenerate crossing configuration at time t c , as in Section 4. In the following we shall discuss only the case ofd h . The same result holds ford min , taking care of the possible exchange of role of two local minima d h , d k as absolute minimum.
Proposition 4. Let Y (t) be a generalized solution of (41) and E(t) as in (39), (40). Assume t c ∈ T (Y ) is a crossing time and E c = E(t c ) is a nondegenerate crossing configuration with only one crossing point. Then there exists an interval (t a , t b ), t a < t c < t b such that d h ∈ C 1 ((t a , t b ); R).
Proof. Let the interval (t a , t b ) be such that E((t a , t b )) ⊂ W , where W is chosen as in Theorem 4.2. We can assume that E(t) ∈ W + for t ∈ (t a , t c ), E(t) ∈ W − for t ∈ (t c , t b ) (the proof for the opposite case is similar). For Here ∇ E , ∇ Y , ∇ E ′ denote the vectors of partial derivatives with respect to E, Y, E ′ respectively. The derivativeĖ ′ (t) exists also for t = t c . On the other hand, by Theorem 4.2, the restrictions of ∇ Y R(E(t)) to t < t c and t > t c admit two different continuous extensions to t c . By (23), sinced h (E(t c )) = 0, we have We describe the algorithm to compute the solutions of the averaged equations (5) beyond the singularity, where R is now the sum of the perturbing functions R i , i = 1 . . . 5, each related to a different planet. We use a Runge-Kutta-Gauss (RKG) method to perform the integration: it evaluates the averaged vector field only at intermediate points of the integration time interval. When the asteroid trajectory is close enough to an orbit crossing, then the time step is decreased to reach the crossing condition exactly.
From Theorem 4.2 we can find two Lipschitz-continuous extensions of the averaged vector field from both sides of the singular set Σ.
To compute the solution beyond the singularity we use the explicit formula (23) giving the difference between the two extended vector fields, either of which corresponds to the averaged vector field on different sides of Σ. We compute the intermediate values of the extended vector field just after crossing, then we correct these values by (23) and use them as approximations of the averaged vector field in (5) at the intermediate points of the solutions, see Figure 3. This RKG algorithm avoids the computation of the extended vector field at the singular points, which may be affected by numerical instability.
A difficulty in the application of this scheme is to estimate the size of a suitable neighborhood W of the crossing configuration E c fulfilling the conditions given in Section 4.

Comparison with the solutions of the full equations
We performed some tests to compare the solutions of the averaged equations (5) with the corresponding components of the solutions of the full equations (4). Here we show two tests with the asteroids 1979 XB and 1620 (Geographos). We considered the system composed by an asteroid and 5 planets, from Venus to Saturn. We selected the 8 values kπ/4, with k = 0 . . . 7, for the initial mean anomaly of the asteroid and the same for the planets. Using the program orbit9, we performed the integration of the system with these 64 different initial conditions (i.e. we selected equal  of the four equinoctial 10 orbital elements h, k, p, q of the asteroid over these evolutions, and compared them with the results of the secular evolution. In Figures 4, 5, we show the results: the crosses indicate the secular evolution, the continuous curve is the mean of full numerical one and the gray region represents the standard deviation from the mean. The correspondence between the solutions is good. During the evolution the distance between the asteroid and the Earth for some initial conditions attains values of the order of 10 −4 au for 1620 (Geographos), and 10 −3 au for 1979 XB. In Figure 5 the Earth crossing singularity is particularly evident near the epoch 3000 AD.
Some numerical tests of the theory introduced in [8], with the planets on circular coplanar orbits, can be found in [10].

An estimate of planet crossing times
The results of Section 6 can be used to estimate the epoch in which the orbit of a near-Earth asteroid will cross that of the Earth. We are interested in particular to study the behavior of those asteroids whose orbits will cross the Earth in the next few centuries, so that they must have a small value of d min already at the present epoch. We can consider, for example, the set of potentially hazardous asteroids (PHAs), which have d min ≤ 0.05 au and absolute magnitude H mag ≤ 22, i.e. they are also large.
In Figure 6 we show 3 different evolutions of the signed orbit distancẽ d min for the PHA 1979 XB. Here we draw the full numerical (solid line), secular (dashed) and secular linearized (dotted) evolution ofd min . By Propo- sition 4 the linearization of the secular evolution d min (t) can give a good approximation also in a neighborhood of a crossing time.
We propose a method to compute an interval J of possible crossing times. We sample the line of variation (LOV), introduced in [16], which is a sort of 'spine' of the confidence region (see also [17]), and compute the signed orbit distanced min for each virtual asteroid (VA) of the sample. Then we compute the time derivative ofd min for each VA and extrapolate the crossing times by a linear approximation of the evolution. We set J = [t 1 , t 2 ], with t 1 , t 2 the minimum and maximum crossing times obtained (see Figure 7). In the computation of J we take into account a band centered at the Earth crossing line d min = 0: in this test the width of the considered band is 2 × 10 −3 au. We describe a method to assign a probability of occurrence of crossings in a given time interval, which is related to the algorithm described above. For each value of the LOV parameter s we have a VA at a time t, so that we can computed min (t). Thus, using the scheme of Figure 7 we can define a map T from the LOV parameter line to the time line. The map T gives the crossing times, using the linearized secular dynamics, for the VAs on the LOV that correspond to the selected values of the parameter s. Moreover, we have a probability density function p(s) on the LOV. Therefore, given an interval I in the time line, we can consider the set U I = T −1 (I) and define the probability of having a crossing in the time interval I as P (I) = U I p(s) ds .
Finally, in Figure 8 we show the corresponding interval J ′ obtained by computing the secular evolution (without linearization) of the orbit distance for each VA of 1979 XB. The sizes of J and J ′ are almost equal, but the left extremum of J ′ is ∼ 10 years before.

Conclusions and future work
We have studied the double averaged restricted 3-body problem in case of orbit crossing singularities, improving and completing the results in [8], [5]. This problem is of interest to study the dynamics of near-Earth asteroids from a statistical point of view, going beyond the Lyapounov times of their orbits. We have also proved that generically, in a neighborhood of a crossing time, the secular evolution of the (signed) orbit distance is more regular than the averaged evolution of the orbital elements. The solutions of this averaged problem have been computed by a numerical method and then compared with the solutions of the full equations in a few test cases. The results were good enough; however, we expect that the averaging technique fails in case of mean motion resonances or close encounters with a planet. We plan to perform numerical experiments with a large sample of near-Earth asteroids showing different behaviors: this will be useful to understand the applicability of the averaging technique to the whole set of NEAs.