Long term dynamics for the restricted N-body problem with mean motion resonances and crossing singularities

TL;DR

This paper develops a mathematical framework to analyze the long-term dynamics of asteroids in mean motion resonance with planets, especially addressing crossing singularities, and demonstrates the dynamical protection mechanism against collisions.

Contribution

It extends existing models by allowing the extension of the vector field across crossing singularities and analyzing the differentiability of the orbit distance near crossings.

Findings

Extended the vector field to handle crossing singularities.

Proved differentiability of the orbit distance near crossing times.

Numerical validation with asteroid classes Alinda and Toro.

Abstract

We consider the long term dynamics of the restricted N-body problem, modeling in a statistical sense the motion of an asteroid in the gravitational field of the Sun and the solar system planets. We deal with the case of a mean motion resonance with one planet and assume that the osculating trajectory of the asteroid crosses the one of some planet, possibly different from the resonant one, during the evolution. Such crossings produce singularities in the differential equations for the motion of the asteroid, obtained by standard perturbation theory. In this work we prove that the vector field of these equations can be extended to two locally Lipschitz-continuous vector fields on both sides of a set of crossing conditions. This allows us to define generalized solutions, continuous but not differentiable, going beyond these singularities. Moreover, we prove that the long term evolution of the 'signed' orbit distance (Gronchi and Tommei 2007) between the asteroid and the planet is differentiable in a neighborhood of the crossing times. In case of crossings with the resonant planet we recover the known dynamical protection mechanism against collisions. We conclude with a numerical comparison between the long term and the full evolutions in the case of asteroids belonging to the 'Alinda' and 'Toro' classes (Milani et al. 1989). This work extends the results in (Gronchi and Tardioli 2013) to the relevant case of asteroids in mean motion resonance with a planet.

Equations

$$\ddot{\mathbf{r}}=-{\mathtt{k}}^2\frac{\mathbf{r}}{|\mathbf{r}{|}^3}+{\mathtt{k}}^2\sum _{j=1}^{N-2}{\mu }_j\left(\frac{{\mathbf{r}}_j-\mathbf{r}}{|{\mathbf{r}}_j-\mathbf{r}{|}^3}-\frac{{\mathbf{r}}_j}{|{\mathbf{r}}_j{|}^3}\right),$$

$$\dot{\mathbf{p}}=-\frac{\partial \mathcal{H}}{\partial \mathbf{r}},\dot{\mathbf{r}}=\frac{\partial \mathcal{H}}{\partial \mathbf{p}}=\mathbf{p},$$

$$\mathcal{H}(\mathbf{p},\mathbf{r},t)=\frac{|\mathbf{p}{|}^2}{2}-\frac{{\mathtt{k}}^2}{|\mathbf{r}|}-{\mathtt{k}}^2\sum _{j=1}^{N-2}{\mu }_j\left(\frac{1}{d_j(\mathbf{r},t)}-\frac{\mathbf{r}\cdot {\mathbf{r}}_j(t)}{|{\mathbf{r}}_j(t){|}^3}\right).$$

$$L$$

$$G$$

$$Z$$

$$\mathtt{n}=\frac{{\mathtt{k}}^4}{L^3}.$$

$$\mathcal{H}={\mathcal{H}}_0+ϵ{\mathcal{H}}_1,$$

$${\mathcal{H}}_0=-\frac{{\mathtt{k}}^4}{2L^2},$$

$${\mathcal{H}}_1=\sum _{j=1}^{N-2}{\mathcal{H}}_1^{(j)},{\mathcal{H}}_1^{(j)}=-{\mathtt{k}}^2\frac{{\mu }_j}{{\mu }_5}\left(\frac{1}{d_j}-\frac{\mathbf{r}\cdot {\mathbf{r}}_j}{|{\mathbf{r}}_j{|}^3}\right),$$

$${\mathcal{H}}_1={\mathcal{H}}_1(L,G,Z,\ell ,g,z,t).$$

$${\ell }_j={\mathtt{n}}_{j}t+{\ell }_j(0),g_j={\mathtt{g}}_{j}t+g_j(0),z_j={\mathtt{s}}_{j}t+z_j(0)$$

$$\mathbf{\ell }=(\ell ,{\ell }_1\dots ,{\ell }_N),$$

$${\mathbf{\ell }}_j=(\ell ,{\ell }_j),$$

$$\tilde{\mathcal{H}}=-\frac{{\mathtt{k}}^4}{2L^2}+\sum _{j=1}^{N-2}({\mathtt{n}}_j{L}_j+{\mathtt{g}}_j{G}_j+{\mathtt{s}}_j{Z}_j)+ϵ{\tilde{\mathcal{H}}}_1(L,G,Z,\mathbf{\ell },\mathbf{g},\mathbf{z}),$$

$${\tilde{\mathcal{H}}}_1=\sum _{j=1}^{N-2}{\widetilde{{\mathcal{H}}_1}}^{(j)},{\tilde{\mathcal{H}}}_1^{(j)}=-{\mathtt{k}}^2\frac{{\mu }_j}{{\mu }_5}\left(\frac{1}{{\tilde{d}}_j}-\frac{\mathbf{r}\cdot {\tilde{\mathbf{r}}}_j}{|{\tilde{\mathbf{r}}}_j{|}^3}\right),$$

$${\tilde{\mathbf{r}}}_j={\tilde{\mathbf{r}}}_j({\ell }_j,g_j,z_j),{\tilde{d}}_j=|{\tilde{\mathbf{r}}}_j-\mathbf{r}|.$$

$$h^{\ast }\mathtt{n}+h_5^{\ast }{\mathtt{n}}_5=0\text{mbox}forsome(h^{\ast },h_5^{\ast })\in {\mathbb{Z}}^2.$$

$$\mathbf{\phi }=(\mathbf{\ell },\mathbf{g},\mathbf{z}),{\mathbf{\phi }}_j=({\mathbf{\ell }}_j,{\mathbf{g}}_j,{\mathbf{z}}_j)$$

$$\mathbf{I}=(\mathbf{L},\mathbf{G},\mathbf{Z}),{\mathbf{I}}_j=({\mathbf{L}}_j,{\mathbf{G}}_j,{\mathbf{Z}}_j)$$

$${\Phi }_{\chi }^{ϵ}({\mathbf{I}}^{\prime },{\mathbf{\phi }}^{\prime })=(\mathbf{I},\mathbf{\phi }),$$

$${\mathcal{H}}^{\prime }=\mathcal{H}\circ {\Phi }_{\chi }^{ϵ}=\mathcal{H}+ϵ\{\mathcal{H},\chi \}+O({ϵ}^2)={\mathcal{H}}_0+ϵ({\mathcal{H}}_1+\{{\mathcal{H}}_0,\chi \})+O({ϵ}^2).$$

$${\mathcal{H}}_1+\{{\mathcal{H}}_0,\chi \}=f$$

$${\mathcal{H}}_1=\sum _{j=1}^{N-2}{\mathcal{H}}_1^{(j)}$$

$${\mathcal{H}}_1^{(j)}=\sum _{(h,h_j)\in {\mathbb{Z}}^2}{\hat{\mathcal{H}}}_{(h,h_j)}^{(j)}{e}^{i(h\ell +h_j{\ell }_j)}.$$

$${\hat{\mathcal{H}}}_{(h,h_j)}^{(j)}={\hat{\mathcal{H}}}_{(h,h_j)}^{(j)}(L,G,Z,{\mathbf{g}}_j,{\mathbf{z}}_j)=\frac{1}{(2\pi {)}^2}{\int }_{{\mathbb{T}}^2}{\mathcal{H}}_1^{(j)}{e}^{-i(h\ell +h_j{\ell }_j)}d\ell d{\ell }_j$$

$$\chi =\sum _{j=1}^{N-2}{\chi }^{(j)},{\chi }^{(j)}={\chi }^{(j)}(L^{\prime },G^{\prime },Z^{\prime },{\mathbf{\ell }}_j^{\prime },{\mathbf{g}}_j^{\prime },{\mathbf{z}}_j^{\prime })$$

$${\hat{\chi }}_{(h,h_j)}^{(j)}={\hat{\chi }}_{(h,h_j)}^{(j)}(L^{\prime },G^{\prime },Z^{\prime },{\mathbf{g}}_j^{\prime },{\mathbf{z}}_j^{\prime })$$

$${\chi }^{(j)}=\sum _{(h,h_j)\in {\mathbb{Z}}^2}{\hat{\chi }}_{(h,h_j)}^{(j)}{e}^{i(h{\ell }^{\prime }+h_j{\ell }_j^{\prime })}.$$

$$\mathcal{H}_{1}+\{\mathcal{H}_{0},\chi\}=\sum_{j=1}^{N-2}\Bigl{(}\mathcal{H}_{1}^{(j)}-\frac{\partial\mathcal{H}_{0}}{\partial{\bf I}}\cdot\frac{\partial\chi^{(j)}}{\partial{\bm{\varphi}}}\Bigr{)},$$

$$\displaystyle\mathcal{H}_{1}^{(j)}-\frac{\partial\mathcal{H}_{0}}{\partial{\bf I}}\cdot\frac{\partial\chi^{(j)}}{\partial{\bm{\varphi}}}=\sum_{(h,h_{j})\in\mathbb{Z}^{2}}\bigl{[}\widehat{\mathcal{H}}_{(h,h_{j})}^{(j)}-i(h\mathtt{n}+h_{j}\mathtt{n}_{j})\widehat{\chi}_{(h,h_{j})}^{(j)}\bigr{]}e^{i(h\ell^{\prime}+h_{j}\ell_{j}^{\prime})}.$$

$$f=\sum _{j=1}^{N-2}{f}_j,$$

Full text

1 Introduction

It is well known that for $N\geq 3$ the $N$-body problem is not integrable, even in the restricted case. In particular, the evolutions of near-Earth asteroids (NEAs) have short Lyapunov times, beyond which the orbit computed by numerical techniques and the true orbit are completely uncorrelated [14]. However, we can obtain statistical information on the long term evolution by considering a normal form of the Hamiltonian of the problem, where we try to filter out the short periodic oscillations. More precisely, we would like to eliminate the dependence on the fast angles from the first order part of the Hamiltonian [1]. Outside mean motion resonances this program can be successfully completed and corresponds to averaging Hamilton’s equations over the mean anomalies of the asteroid and the planets. In case of mean motion resonances, the resonant combination of the mean anomalies is a slow angle and must be retained in the normal form.

In both cases, the elimination of the fast angles is usually obtained through a canonical transformation, in the spirit of classical perturbation theory. However, the intersections between the trajectories of the asteroid and the planets introduce singularities in the standard procedure. Actually, even the coefficients of the Fourier series expansion of the generating function are not defined in a neighborhood of crossings. On the other hand, since the trajectory of a near-Earth asteroid is likely to cross the trajectory of the Earth, we cannot avoid to deal with these problems. Note that the minimal distance between the trajectories of an asteroid and a planet is crucial in the study of possible Earth impactors. Actually, a small value of this quantity, that we denote by $d_{min}$, is a necessary condition for an impact. An orbit crossing singularity occurs whenever $d_{min}=0$.

After the preliminary study by Lidov and Ziglin [8], in the case of orbits uniformly close to a circular one, the problem of averaging over crossing orbits was studied in [5]. Here the authors assumed the orbits of the planets being circular and coplanar, and excluded mean motion resonances and close approaches with them. In [4] the results were extended to the case of non-zero eccentricities and inclinations. In these works, the main singular term is computed through a Taylor expansion centered at the mutual nodes of the osculating orbits. These results were improved in [7], where the main singular term is expanded at the minimum distance points (see Section 4) and where it is proved that the averaged vector field admits two different Lipschitz-continuous extensions in a neighborhood of almost every crossing configuration. The latter property allows us to define a generalized solution, representing the secular evolution of the asteroid, that is continuous but not differentiable at crossings. Moreover, one can suitably choose the sign of $d_{min}$ and obtain a map $\tilde{d}_{min}$ that is differentiable in a neighborhood of almost all crossing configurations [6]. The secular evolution of $\tilde{d}_{min}$ along the generalized solutions turns out to be differentiable in a neighborhood of the singularity.

The basic model considered in these works comes from the averaging principle. Therefore, it is assumed that the dynamics is not affected by mean motion resonances. However, the population of resonant NEAs is not negligible. Moreover, mean motion resonances are considered responsible for a relatively fast change in the orbital elements leading some asteroids to cross the planet trajectories [15]. Hence it is important to extend the analysis to such asteroids, which is the purpose of this paper.

For the resonant case, the averaging process suffers the presence of small divisors. Hence, the dependence on the mean anomalies cannot be completely eliminated, and the terms corresponding to their resonant combination still appear in the resonant normal form, see (7). We observe that in this relation the averaged Hamiltonian considered in [7] is still present. However, a new term appears in the form of a Fourier series, that we truncate to some order $n_{\rm max}$. This term, denoted by $\mathcal{H}_{res}^{n_{\rm max}}$, is singular at orbit crossings and needs to be studied. Another difference with the non-resonant case is that the semimajor axis of the asteroid orbit is not constant, and the number of state variables to consider in the equations is six.

We will prove that, despite these differences, the vector field of the resonant normal form computed outside the singularities admits two different locally Lipschitz-continuous extensions on both sides of a set of crossing conditions, as in [7]. We can also define generalized solutions, continuous but not differentiable, going beyond the crossing singularities and the long term evolution of the map $\tilde{d}_{min}$ along these solutions is differentiable in a neighborhood of crossings.

The analysis of the singularity is performed in two different ways, depending if the crossed planet is the one in mean motion resonance with the asteroid or not. In case of crossings with the resonant planet we show that, in the limit for $n_{\rm max}\to\infty$, we recover the known dynamical protection mechanism against collisions between the asteroid and the planet [9].

The article is organized as follows. In Section 2 we derive the equations of the long term dynamics outside the crossing singularities for a given mean motion resonance. In Section 3 we recall the definition of the signed orbit distance $\tilde{d}_{min}$. The main results are stated and proved in Section 4. In Section 5 we define the generalized solutions and prove the regularity of the evolution of $\tilde{d}_{min}$. In Section 6 we show the relation between the resonant normal form that we use and the averaged Hamiltonian used in the literature, recovering the dynamical mechanism that protects from collisions. We conclude with some numerical examples in Section 7, showing the agreement between the long term evolution and the full evolution in a statistical sense.

2 The equations for the long term evolution

We consider the differential equations

[TABLE]

where $\bm{r}$ describes, in heliocentric coordinates, the motion of a massless asteroid under the gravitational attraction of the Sun and $N-2$ planets. The heliocentric motions of the planets $\bm{r}_{j}=\bm{r}_{j}(t)$ are known functions of the time $t$ that never vanish: that is we exclude collisions between a planet and the Sun. Moreover, $\mathtt{k}=\sqrt{\mathcal{G}m_{0}}$ is Gauss’s constant, $\mu_{j}=m_{j}/m_{0}$ with $m_{0}$ the mass of the Sun and $m_{j}$ the mass of the $j$-th planet. Equations (1) can be written in Hamiltonian form as

[TABLE]

with Hamiltonian

[TABLE]

In (2) $d_{j}=|\bm{r}_{j}-\bm{r}|$ stands for the distance between the asteroid and the $j$-th planet. We use Delaunay’s elements $(L,G,Z,l,g,z)$ defined by

[TABLE]

where $a,e,I,\Omega,\omega,t_{0}$ represent semimajor axis, eccentricity, inclination, longitude of the ascending node, argument of perihelion, and epoch of passage at perihelion. For the definition of $\ell$ we use the mean motion

[TABLE]

In these coordinates, the Hamiltonian (2) can be written as

[TABLE]

with $\epsilon=\mu_{5}$,

[TABLE]

and

[TABLE]

and $\bm{r}_{j}=\bm{r}_{j}(t)$. Note that in (3)

[TABLE]

To eliminate the dependence on time in $\mathcal{H}_{1}$ we overextend the phase space. We assume that the planets move on quasi-periodic orbits with three independent frequencies $\mathtt{n}_{j},\mathtt{g}_{j},\mathtt{s}_{j}$.

This is the case considered by Laplace (see for example [11]), where the mean semi-major axis $a_{j}$ is constant and the mean value of the mean anomaly $\ell_{j}$ grows linearly with time, i.e. up to a phase, $\ell_{j}=\mathtt{n}_{j}t$. Here $\mathtt{n}_{j}$ is the mean motion of planet $j$. Moreover, every planet is characterized by two more frequencies $\mathtt{g}_{j},\mathtt{s}_{j}$, describing the slow motions of the other mean orbital elements. We introduce the angles

[TABLE]

and their conjugate variables $L_{j},G_{j},Z_{j}$.

Note that these variables do not correspond to the Delaunay’s elements of planet $j$, since they are functions of the orbital elements of the asteroid and planet $j$. We use the following notation:

[TABLE]

and analogously we define $\bm{L},\bm{G},\bm{Z},\bm{L}_{j},\bm{G}_{j},\bm{Z}_{j}$.

The dynamics in this overextended phase space is determined by the autonomous Hamiltonian

[TABLE]

where

[TABLE]

with

[TABLE]

Here we are assuming that $\bm{r}_{j}$ evolves according to Laplace’s solution for the planetary motions, and we write it as a function of its frequencies, denoted by $\tilde{\bm{r}}_{j}$. Hereafter we shall omit the ’tilde’, to simplify the notation.

The frequencies $\mathtt{g}_{j}$ and $\mathtt{s}_{j}$ are of order $\epsilon$ [11]. In order to study the secular dynamics, we would like to eliminate all the frequencies corresponding to the fast angles $\bm{\ell}$. In case of a mean motion resonance with a planet this is not possible.

In the following we shall assume that there is only one mean motion resonance with a planet and no close approaches occur. To expose our result we shall consider a $|h_{5}^{*}|:|h^{*}|$ mean motion resonance with Jupiter given by

[TABLE]

A mean motion resonance with another planet can be treated in a similar way. We denote by

[TABLE]

the vectors of the angles and by

[TABLE]

the corresponding vectors of the actions.

We use the Lie method [11] to search for a suitable canonical transformation close to the identity, that is we search for a function $\chi=\chi({\bf I}^{\prime},{\bm{\varphi}}^{\prime})$ such that the inverse transformation is

[TABLE]

where $\Phi^{t}_{\chi}$ is the Hamiltonian flow associated to $\chi$. The function $\chi$ is selected so that the transformed Hamiltonian $\mathcal{H}^{\prime}=\mathcal{H}\circ\Phi^{\epsilon}_{\chi}$ depends, at least at first order, on as less fast angular variables as possible. Using a formal expansion in $\epsilon$ we have

[TABLE]

In the resonant case we search for a solution $\chi$ of the equation

[TABLE]

for some function $f=f({\bf I}^{\prime},h^{*}\ell^{\prime}+h_{5}^{*}\ell_{5}^{\prime},\bm{g}^{\prime},\bm{z}^{\prime})$. To solve (5) we restrict to the case where no orbit crossings with the planets occur. We shall see in the next sections how we can deal with the case of crossings.

We develop

[TABLE]

in Fourier’s series of the fast angles:

[TABLE]

Here

[TABLE]

are the Fourier coefficients. We observe that $\widehat{\mathcal{H}}_{(h,h_{j})}^{(j)}$ are defined also in case of orbit crossings, since the integral in (6) converges (see e.g. [7]).

Moreover, we can write $\chi$ as

[TABLE]

and search for the coefficients

[TABLE]

in the Fourier series development

[TABLE]

Inserting these Fourier developments into (5) we obtain

[TABLE]

where

[TABLE]

This expression suggests to choose the function $f$ in (5) in the following form:

[TABLE]

where $f_{5}=f_{5}({\bf I}^{\prime}_{5},h^{*}\ell^{\prime}+h_{5}^{*}\ell_{5}^{\prime},\bm{g}^{\prime}_{5},\bm{z}^{\prime}_{5})$ and $f_{j}=f_{j}({\bf I}^{\prime}_{j},\bm{g}^{\prime}_{j},\bm{z}^{\prime}_{j})$ for $j\neq 5$. This can be accomplished by choosing

[TABLE]

when the denominator does not vanish. Hence, we exclude the case $(h,h_{j})=(0,0)$ and the resonant case $(h,h_{5})=n(h^{*},h_{5}^{*})$ for some $n\in\mathbb{Z}^{*}=\mathbb{Z}\setminus\{0\}$, for which we assume that the corresponding Fourier coefficient of $\chi$ vanishes. With this choice we have

[TABLE]

We truncate the Fourier series to some order $n_{\rm max}$ and consider

[TABLE]

as resonant normal form of the Hamiltonian, where

[TABLE]

and

[TABLE]

with $\Re(z)$ the real part of $z\in\mathbb{C}$, where we used $\overline{\widehat{\mathcal{H}}_{(h,h_{5})}^{(5)}}=\widehat{\mathcal{H}}_{(-h,-h_{5})}^{(5)}$. For simplicity, we shall write $\mathscr{H}$, $\mathcal{H}_{res}$ in place of $\mathscr{H}_{n_{\rm max}}$, $\mathcal{H}_{res}^{n_{\rm max}}$. It is easy to see that, for every $j$,

[TABLE]

being null the average of the indirect perturbation (see [3]). We observe that in the Fourier coefficient $\widehat{\mathcal{H}}_{n(h^{*},h_{5}^{*})}^{(5)}$ the term corresponding to the indirect perturbation does not vanish. We can write

[TABLE]

where

[TABLE]

with $I_{5}^{c,n},I_{5}^{s,n}$ depend on $L,G,Z,\bm{g}_{5},\bm{z}_{5}$.

Moreover, since the new Hamiltonian does not depend on $\ell_{j}$ for $j\neq 5$ we have

[TABLE]

We now introduce the resonant angle $\sigma$ through the canonical transformation

[TABLE]

with

[TABLE]

We chose the matrix $A$ so that $L$ does not depend on $S_{5}$. For this reason we could not use a unimodular matrix. However, this will not affect our analysis.

We shall still denote by

[TABLE]

the resonant normal form of the Hamiltonian in these new variables, with

[TABLE]

Since the Hamiltonian does not depend on $\sigma_{5}$, the value of $S_{5}$ will remain constant and we will treat it as a parameter. Calling $\mathcal{Y}=(S,G,Z,\sigma,g,z)$ we consider the equations for the motion of the asteroid given by

[TABLE]

where

[TABLE]

is the symplectic identity of order $6$. In components, system (9) is written as

[TABLE]

where $\mathcal{H}_{res}$, $\overline{\mathcal{H}}_{1}$ are functions of $(S,G,Z,\sigma,{\bm{g}_{5},\bm{z}_{5}})$ and $(S,G,Z,{\bm{g},\bm{z}})$ respectively. Since $\epsilon C_{j}=-\mathtt{k}^{2}\mu_{j}$, we get

[TABLE]

The derivatives of $\mathcal{H}_{res}$ and $\overline{\mathcal{H}}_{1}$ are not defined at orbit crossings with the planets. In the following sections we shall discuss how we can define generalized solutions of system (9) in case of orbit crossings.

3 The orbit distance

We recall here some facts and notations from [6], [7]. Let $(E,v)$, $(E^{\prime},v^{\prime})$ be two sets of orbital elements, where $E,E^{\prime}$ describe the trajectories of the asteroid and one planet, $v,v^{\prime}$ describe the position of these bodies along them. Denote by $\mu^{\prime}$ the ratio of the mass of this planet to the mass of the Sun. We also introduce the notation $\mathcal{E}=(E,E^{\prime})$ for the two-orbit configuration and $V=(v,v^{\prime})$ for the vector of parameters along the orbits. We denote by $\mathcal{X}=\mathcal{X}(E,v)$ and $\mathcal{X}^{\prime}=\mathcal{X}^{\prime}(E^{\prime},v^{\prime})$ the Cartesian coordinates of the asteroid and the planet respectively. For each given $\mathcal{E}$, $V_{h}(\mathcal{E})$ represents a local minimum point of the function

[TABLE]

We introduce the local minimum maps

[TABLE]

and the orbit distance

[TABLE]

We shall consider non-degenerate configurations $\mathcal{E}$, i.e such that all the critical points of the map $V\mapsto d(\mathcal{E},V)$ are non-degenerate. In this way, we can always choose a neighborhood $\mathcal{W}$ of $\mathcal{E}$ where the maps $d_{h}$ do not have bifurcations. A crossing configuration is a two-orbit configuration $\mathcal{E}_{c}$ such that $d(\mathcal{E}_{c},V_{h}(\mathcal{E}_{c}))=0$ where $V_{h}(\mathcal{E}_{c})$ is the corresponding minimum point. The maps $d_{h}$ and $d_{min}$ are singular at crossing configurations, and their derivatives in general do not exists. Anyway, it is possible to obtain analytic maps in a neighborhood of a crossing configuration $\mathcal{E}_{c}$ by a suitable choice of the sign for these maps. We summarize here the procedure to deal with this singularity for $d_{h}$; the procedure for $d_{min}$ is the same. Let $V_{h}=(v_{h},v_{h}^{\prime})$ be a local minimum point of $d^{2}$ and let $\mathcal{X}_{h}=\mathcal{X}_{h}(E,v_{h})$ and $\mathcal{X}_{h}^{\prime}=\mathcal{X}_{h}^{\prime}(E^{\prime},v_{h}^{\prime})$. We introduce the vectors tangent to the trajectories defined by $E,E^{\prime}$ at these points

[TABLE]

and their cross product $\tau_{h}^{*}=\tau_{h}^{\prime}\times\tau_{h}$. Both vectors $\tau_{h},\tau^{\prime}_{h}$ are orthogonal to $\Delta_{h}=\mathcal{X}^{\prime}_{h}-\mathcal{X}_{h}$, so that $\tau_{h}^{*}$ is parallel to $\Delta_{h}$, see Figure 1.

Denoting by $\hat{\tau}_{h}^{*}$, $\hat{\Delta}_{h}$ the corresponding unit vectors, we consider the local minimal distance with sign

[TABLE]

This map is analytic in a neighborhood of most crossing configurations. Actually, this smoothing procedure fails in case the vectors $\tau_{h},\tau^{\prime}_{h}$ are parallel.

Finally, given a neighborhood $\mathcal{W}$ of $\mathcal{E}_{c}$ without bifurcations of $d_{h}$, we write $\mathcal{W}=\mathcal{W}^{-}\cup\Sigma\cup\mathcal{W}^{+}$, where

[TABLE]

4 Extraction of the singularities

In the following we shall expose a method to investigate the crossing singularities occurring in (9). For simplicity, we shall eventually drop the index 5, referring to Jupiter, and denote simply by a prime the quantities referring to the crossed planet.

Let $\mathcal{E}_{c}$ be a two-orbit crossing configuration and suppose that the trajectories are described by the vector $E=(S,G,Z,g,z)$. In the following we shall write $y_{i}$ for the components of the vector $E$. We choose the mean anomalies as parameters along the trajectory so that $V=(\ell,\ell^{\prime})$. The first step of our analysis is to consider, for each $\mathcal{E}$ in a neighborhood $\mathcal{W}$ of $\mathcal{E}_{c}$, the Taylor expansion of $V\mapsto d(\mathcal{E},V)$ in a neighborhood of $V_{h}=V_{h}(\mathcal{E})$, i.e.

[TABLE]

where ${\cal R}_{3}^{(h)}$ is the remainder in the integral form, and define the approximated distance

[TABLE]

with

[TABLE]

The matrix ${\cal A}_{h}$ is positive definite except for tangent crossings, where it is degenerate. To study the crossing singularities in case of a mean motion resonance with Jupiter we distinguish between the case where the asteroid trajectory crosses the trajectory of another planet and the case where it crosses the trajectory of Jupiter itself. In the first case the crossing singularity appears only in the averaged terms $\frac{\partial\overline{\mathcal{H}}_{1}}{\partial y_{i}}$. In the second case also the derivatives $\frac{\partial I_{5}^{s,n}}{\partial y_{i}}$, $\frac{\partial I_{5}^{c,n}}{\partial y_{i}}$ are affected by this singularity. In both cases the component $\frac{\partial\mathscr{H}}{\partial\sigma}$ is regular.

We obtain the following results.

Theorem 1**.**

Let $\mathcal{E}_{c}$ be a non-degenerate crossing configuration with a planet (including Jupiter). Then, there exists a neighborhood $\mathcal{W}$ of $\mathcal{E}_{c}$ such that for each $i=1,\ldots,5$ we can define two maps

[TABLE]

that are Lipschitz-continuous extensions of the maps

[TABLE]

Moreover, the following relation holds in $\mathcal{W}$:

[TABLE]

Proof.

We can show this result by following the same steps as in [7, Theorem 4.2], replacing $R$ by $-\epsilon\overline{\mathcal{H}}_{1}$.

∎

Theorem 2**.**

Let ${\sf h}=(h^{*},h^{*}_{5})$ and $\mathcal{E}_{c}$ be a non-degenerate crossing configuration with Jupiter. Then, there exists a neighborhood $\mathcal{W}$ of $\mathcal{E}_{c}$ such that, for every $n>0$ and for each $i=1,\ldots,5$, we can define four maps

[TABLE]

that are Lipschitz-continuous extensions of the maps

[TABLE]

respectively. Moreover, the following relations hold in $\mathcal{W}$:

[TABLE]

Before giving a proof of Theorem 2 we state some consequences of both theorems. We define the following locally Lipschitz-continuous maps, extending the vector field of Hamilton’s equations (9) in a neighborhood of the crossing singularity,

[TABLE]

where we use the definition above in case of crossings with Jupiter, and the one below for crossings with other planets. Here $\mathcal{H}_{0}$, $\mathcal{H}_{res}$ are defined as in (8), and

[TABLE]

Moreover, we consider the map

[TABLE]

Corollary 1**.**

If $\mathcal{E}_{c}$ corresponds to a crossing configuration with a planet different from Jupiter, then the following relation holds in $\mathcal{W}$:

[TABLE]

Corollary 2**.**

If $\mathcal{E}_{c}$ corresponds to a crossing configuration with Jupiter, then the following relation holds in $\mathcal{W}$

[TABLE]

We recall that, for each $N\in\mathbb{N}$ and $x\neq 2h\pi$, with $h\in\mathbb{Z}$, we have

[TABLE]

where

[TABLE]

is the Dirichlet kernel (see [13]).

Remark 1**.**

With the notation above we have

[TABLE]

that for $n_{\rm max}\to\infty$ converges in the sense of distributions to the Dirac delta $\delta_{\sigma_{c}}$ centered in $\sigma_{c}:={\sf h}\cdot V_{h}$.

Remark 2**.**

The component $\frac{\partial\mathscr{H}}{\partial\sigma}$ is locally Lipschitz-continuous.

4.1 Proof of Theorem 2

We shall prove the result only for the maps (12), the proof for (13) being similar. Since we assume that Jupiter cannot collide with the Sun, the term ${\bf r}_{5}$ will never vanish, so that we study only the derivatives

[TABLE]

for a fixed value of $n\in\mathbb{N}$. We shall refer to some estimates and results proved in [7]. For the reader’s convenience we collect them in Appendix A. Moreover, we shall denote by $C_{k}$, $k=1,\dots,12$, some positive constants independent on $\mathcal{E}$.

Let $\mathcal{E}_{c}$ be a non-degenerate crossing configuration. Let us choose two neighborhoods $\mathcal{W}$ of $\mathcal{E}_{c}$ and ${\cal U}$ of $(\mathcal{E}_{c},V_{h}(\mathcal{E}_{c}))$, as in Lemma 1 in the Appendix. To investigate the crossing singularity we can restrict the integral above to the set

[TABLE]

for some $r>0$. We first note that

[TABLE]

and prove that the first three addenda have a continuous extension to $\mathcal{W}$. From the estimate (36) the map

[TABLE]

admits a continuous extension to $\mathcal{W}$. We now prove that also the map

[TABLE]

admits a continuous extension to $\mathcal{W}$. Indeed we note that

[TABLE]

By (27), (37) the first addendum in the r.h.s. of (16) is summable. For the second, by (29) we get

[TABLE]

From the estimate

[TABLE]

we can conclude using (30).

The existence of a continuous extension to $\mathcal{W}$ of the maps

[TABLE]

comes from (27).

The last term cannot be extended with continuity at crossings. Using Lemma 3 we define the two maps

[TABLE]

that are continuous extensions to $\mathcal{W}$ of the restrictions of $\frac{\partial}{\partial y_{i}}\int_{{\cal D}}\frac{1}{\delta_{h}}dV$ to $\mathcal{W}^{\pm}$ respectively. Then we set

[TABLE]

To conclude the proof we just need to prove that these maps are Lipschitz-continuous. We establish the result by proving that the function

[TABLE]

is uniformly bounded in $\mathcal{W}\setminus\Sigma$. Let us consider the Taylor expansion

[TABLE]

where

[TABLE]

is the remainder in integral form, so that in ${\cal U}$ we have

[TABLE]

for some $C>0$. Using the approximated distance $\delta_{h}$ defined in (11) we can write $F(\mathcal{E})$ as sum of four terms:

[TABLE]

where

[TABLE]

We prove that each term $F_{i}$ is bounded by a constant independent on $\mathcal{E}$. The boundedness of $F_{1}$ comes trivially from (28). From the relation

[TABLE]

and the estimates (26),(29),(31) we obtain

[TABLE]

Then (17) and (30) yield the boundedness of $F_{4}$:

[TABLE]

To show the boundedness of $F_{2}$ we just need to prove that

[TABLE]

so that

[TABLE]

Using $d^{2}=\delta_{h}^{2}+{\cal R}_{3}^{(h)}$ we get

[TABLE]

We prove that each of the four terms in the previous sum satisfies an estimate like (18). For the second term we use estimates (31),(32), for the third (29),(33), and for the last (34). To estimate the first term we note that

[TABLE]

and use

[TABLE]

We can conclude using (26),(29),(33),(35).

Now we show the boundedness of $F_{3}$. We write

[TABLE]

and study the two integrals in the r.h.s. separately. To estimate the first we use (11) and get

[TABLE]

so that

[TABLE]

Then we use the change of variables $\xi=\mathcal{A}_{h}^{1/2}(V-V_{h})$ and polar coordinates $(\rho,\theta)$ defined by $\xi=\rho(\cos\theta,\sin\theta)$. We distinguish between terms with even and odd degree in $(V-V_{h})$. First we consider the ones with even degree. The term of degree $2$ is estimated as follows

[TABLE]

while for the term of degree $4$ we note that

[TABLE]

for some functions $b_{\gamma}$, $c_{\gamma}$, uniformly bounded in ${\cal W}\setminus\Sigma$, and for $\gamma=(\gamma_{1},\gamma_{2})\in(\mathbb{N}\cup\{0\})^{2}$. The terms with odd degree in $(V-V_{h})$ vanish, as can be shown by similar computations, using

[TABLE]

with $\gamma_{1}+\gamma_{2}$ odd. To estimate the second integral in (19) we proceed in a similar way, using

[TABLE]

Remark 3**.**

If $\mathcal{E}_{c}$ is an orbit configuration with two crossings, assuming that $d_{h}(\mathcal{E}_{c})=0$ for $h=1,2$, we can extract the singularity by considering the approximated distances $\delta_{1},\delta_{2}$ and considering $1/d$ as sum of the three terms $(1/d-1/\delta_{1}-1/\delta_{2})$, $1/\delta_{1}$, $1/\delta_{2}$.

5 Generalized solutions and evolution of the orbit distance

Following [7, Sections 5-6] we can construct generalized solutions by patching classical solutions defined in the domain $\mathcal{W}^{+}$ with classical solutions defined on $\mathcal{W}^{-}$ and vice-versa. Let $(E(t),\sigma(t))$, with $E(t)=(S(t),G(t),Z(t),g(t),z(t))$, represent the evolution of the asteroid according to (9). In a similar way we denote by $E^{\prime}(t)$ a known function of time representing the evolution of the trajectory of the planet. Setting $\mathcal{E}(t)=(E(t),E^{\prime}(t))$ we let $T(\mathcal{Y})$ be the set of times $t_{c}$ such that $d_{min}(\mathcal{E}(t_{c}))=0$ and suppose that it has no accumulation points.

We say that $\mathcal{Y}(t)$ is a generalized solution of (9) if it is a classical solution for $t\notin T(\mathcal{Y})$ and for each $t_{c}\in T(\mathcal{Y})$ there exist finite values of

[TABLE]

In order to construct a generalized solution we consider a solution $\mathcal{Y}(t)$ of the Cauchy problem given by (9) with a non crossing initial condition $\mathcal{Y}(t_{0})$. Suppose that it is defined on a maximal interval $J$ such that $\sup J=t_{c}\in T(\mathcal{Y})$ and that $\mathcal{Y}(t)\in\mathcal{W}^{+}$ as $t\to t_{c}$ . Suppose that the crossing is occurring with a planet different from Jupiter (resp. Jupiter itself). Applying Theorem 1-(a) (resp. Theorems 1-(a) and 2-(a)) we have that there exists

[TABLE]

and the solution can be extended beyond $t_{c}$ considering the Cauchy problem

[TABLE]

for some $\tau\to t_{c}$, so that we call $\mathcal{Y}(t_{c})=\mathcal{Y}_{c}$. Using again Theorem 1-(a) (resp. Theorems 1-(a) and 2-(a)), we can extend the solution beyond the singularity considering the new Cauchy problem

[TABLE]

whose solution fulfills, from Corollary 1 (resp. Corollary 2)

[TABLE]

Note that the evolution of the orbital elements according to a generalized solution is continuous but not differentiable in a neighborhood of a crossing singularity. More precisely, the evolution of the elements $(G,Z,\sigma,g,z)$ is only Lipschitz-continuous while the evolution of $S$ is $C^{1}$, since $\frac{\partial\mathscr{H}}{\partial\sigma}$ is continuous also at orbit crossings.

Once a generalized solution $\mathcal{Y}(t)=(E(t),\sigma(t))$ is defined, we can consider the evolution of the distance $\tilde{d}_{h}(\mathcal{E}(t))$. Let us define

[TABLE]

and suppose that it is defined in an interval containing a crossing time $t_{c}$ corresponding to a non-degenerate crossing configuration. We have the following

Proposition 1**.**

Let $\mathcal{Y}(t)$ be a generalized solution of (9) and $\mathcal{E}(t)$ be defined as above. Suppose that $t_{c}$ is a crossing time such that $\mathcal{E}_{c}=\mathcal{E}(t_{c})$ is a non-degenerate crossing configuration. Then there exists an open interval $I\ni t_{c}$ such that $\bar{d}_{h}\in C^{1}(I,\mathbb{R})$.

Proof.

We choose the interval $I$ such that $\mathcal{E}(I)\in\mathcal{W}$ with $\mathcal{W}$ defined in Theorem 1 (resp. 2) and suppose that $\mathcal{E}(t)\in\mathcal{W}^{+}$ for $t<t_{c}$ and $\mathcal{E}(t)\in\mathcal{W}^{-}$ for $t>t_{c}$. We can compute, for $t\neq t_{c}$,

[TABLE]

The second addendum is continuous while for the first we need to distinguish between crossing a planet different from Jupiter (the resonant planet) and crossing Jupiter itself. In the first case, we apply Corollary 1 and obtain

[TABLE]

where $\{,\}$ are the Poisson brackets.

In the second case, we apply Corollary 2 and get

[TABLE]

∎

6 Dynamical protection from collisions

In case of crossings with the resonant planet, the resonance protects the asteroid from close encounters with that planet (see [9]). This protection mechanisms is usually derived by a perturbative approach different from ours. Here we describe how this mechanism can be recovered from the normal form (8) in the limit for $n_{\rm max}\to\infty$.

Let us consider, for simplicity, a restricted 3-body problem Sun-planet-asteroid, where the asteroid is in a mean motion resonance with the planet, given by

[TABLE]

and their trajectories cross each other during the evolution. In the following we take a Hamiltonian containing only the direct part of the perturbation, the indirect part being regular. Therefore we set

[TABLE]

where $d$ is the distance between the asteroid and the planet. We consider the following procedures:

(I) Through a unimodular transformation $\Psi$ of the fast variables $V=(\ell,\ell^{\prime})$ we pass to new variables $(\sigma,\tau)$, with

[TABLE]

whose evolution occurs on different time scales: $\sigma$ has a long-term evolution, $\tau$ has a fast evolution. More precisely we have

[TABLE]

where $W=(\sigma,\tau)^{T}$ and $\mathcal{U}$ is a constant unimodular matrix whose first raw is $(h,h^{\prime})$. The transformation $\Psi$ can be extended to a canonical transformation (here denoted again by $\Psi$) by defining the corresponding actions as $(S,T)=\mathcal{U}^{-T}(L,L^{\prime})$ and leaving the other variables unchanged. Then, we average over the fast variable $\tau$ and get the Hamiltonian

[TABLE]

Here $X$ is the vector of the other variables, evolving on a secular time scale. This procedure is used e.g. in [9].

(II) As in Section 2, we consider the resonant normal form obtained by eliminating all the non resonant harmonics from the Fourier series of the Hamiltonian. For each integer $N$ we take the partial Fourier sums

[TABLE]

where

[TABLE]

and

[TABLE]

in which we denote by ${\cal V}$ the vector $(\ell,\ell^{\prime})$ when the latter are integration variables. We formally define

[TABLE]

Note that

[TABLE]

where $D_{N}(x)$ is the Dirichlet kernel. We introduce the functions

[TABLE]

Indeed both $\mathscr{K}_{N}$ and $\mathscr{K}_{\infty}$ do not depend on $\tau$. The Hamiltonian $\mathscr{K}_{N}$ corresponds to the resonant normal form in (8). However, here we used a unimodular matrix ${\cal U}$ in the canonical transformation.

Moreover, we observe that the Hamiltonian $\overline{\cal K}$ defined in (21) can be written as a pointwise limit for $N\to\infty$ of the partial Fourier sums

[TABLE]

Let $\sigma_{c}={\sf h}\cdot V_{h}$. If $d_{h}=0$ then $\sigma_{c}$ is the value of $\sigma$ allowing a collision, occurring for $V=V_{h}$. Assume that ${\cal E}_{c}$ is a non-degenerate crossing configuration, i.e. $d_{h}=0$ and ${\cal A}_{h}$ is positive definite. We use $Y=Y({\cal E})$ to denote the variables different from $\sigma$ and we set $Y_{c}=Y({\cal E}_{c})$.

Proposition 2**.**

The following properties hold.

If ${\cal E}\neq{\cal E}_{c}$, then for each $\sigma$ we have

• i)

$\mathscr{K}_{N}(\sigma;Y)=\overline{\cal K}_{N}(\sigma;Y),\quad\forall N$**

• ii)

$\mathscr{K}_{\infty}(\sigma;Y)=\overline{\cal K}(\sigma;Y)$.

Moreover, these functions are differentiable with continuity with respect to $Y$. 2. 2.

For ${\cal E}={\cal E}_{c}$ we have

• i)

$\mathscr{K}_{N}(\sigma;Y_{c})=\overline{\cal K}_{N}(\sigma;Y_{c}),\qquad\forall N,\ \forall\sigma$

• ii)

$\mathscr{K}_{\infty}(\sigma;Y_{c})=\overline{\cal K}(\sigma;Y_{c}),\quad\qquad\forall\sigma\neq\sigma_{c}$

• iii)

$\lim_{\sigma\to\sigma_{c}}\mathscr{K}_{\infty}(\sigma;Y_{c})=\lim_{\sigma\to\sigma_{c}}\overline{\cal K}(\sigma;Y_{c})=+\infty$. 3. 3.

If ${\cal E}={\cal E}_{c}$ and $\sigma\neq\sigma_{c}$ then, denoting by $y_{j}$ a generic component of $Y$,

• i)

the derivatives $\displaystyle\frac{\partial\mathscr{K}_{\infty}}{\partial y_{j}}(\sigma;Y_{c})=\displaystyle\frac{\partial\overline{\cal K}}{\partial y_{j}}(\sigma;Y_{c})$ exist and are continuous;

• ii)

the derivatives $\displaystyle\frac{\partial\mathscr{K}_{N}}{\partial y_{j}}(\sigma;Y_{c})=\displaystyle\frac{\partial\overline{\cal K}_{N}}{\partial y_{j}}(\sigma;Y_{c})$ generically do not exist. 4. 4.

For each $N$ and for each value of $\sigma$ there exist the limits

[TABLE]

from both sides of the crossing configuration set $\Sigma$. These limits are generically different and their difference converges in the sense of distributions, for $N\to\infty$, to the Dirac delta relative to $\sigma_{c}$, multiplied by the factor

[TABLE]

Remark 4**.**

If ${\cal E}={\cal E}_{c}$, procedure (I) gives a well defined vector field, provided that $\sigma\neq\sigma_{c}$. On the other hand, with procedure (II) it does not make sense to consider

[TABLE]

However, for each $N$ we can extend the vector field of $\mathscr{K}_{N}$ in two different ways on $\Sigma$, and the difference between the two extensions has a very weak behavior for $N\to\infty$: it tends to a Dirac delta in the sense of distribution, being the singularity of the delta just at $\sigma=\sigma_{c}$.

Proof of Proposition 2.

1. For every $N$, by applying the change of variables $V\to\Psi(V)$ and Fubini-Tonelli’s theorem we obtain

[TABLE]

that proves i). Point ii) comes from the fact that, for ${\cal E}\neq{\cal E}_{c}$, $\overline{\cal K}(\sigma;Y)$ is a smooth function of $\sigma$ and the corresponding Fourier series converge pointwise for every $\sigma$. Hence we can pass to the limit as $N\to\infty$ in the previous equality.

The differentiability comes from the fact that the distance function ${\cal H}=1/d$ is bounded for ${\cal E}\neq{\cal E}_{c}$.

2. To prove i), we can repeat the argument used in (22). Indeed, the double integral is finite also for ${\cal E}={\cal E}_{c}$ and we can apply Fubini-Tonelli’s theorem.

To prove ii), we recall that the Fourier series of an $L^{1}$ function converges pointwise at every point of differentiability [13]. Therefore, for every $\sigma\neq\sigma_{c}$, $\overline{\cal K}_{N}(\sigma;Y_{c})\to\overline{\cal K}(\sigma;Y_{c})$ for $N\to\infty$. Hence, using i) and passing to the limit for $N\to\infty$ in $\mathscr{K}_{N}$ we get the result.

To prove iii) we just need to prove that one of the two limits diverges. From Fatou’s lemma

[TABLE]

We can prove that the integral in (23) diverges by a singularity extraction technique. Let us write

[TABLE]

The first term in the r.h.s. of (24) is bounded, while the integral of the second diverges because

[TABLE]

where

[TABLE]

with ${\cal U}$ the unimodular matrix defined in (20), and

[TABLE]

The number $b_{22}$ in (25) is defined by

[TABLE]

and is strictly positive because ${\cal B}_{h}$ is positive definite, being ${\cal E}_{c}$ non-degenerate (and therefore ${\cal A}_{h}$ positive definite).

3. Estimate (25), decomposition (24), and the theorem of differentiation under the integral sign yield the existence and continuity of the derivatives $\displaystyle\frac{\partial\overline{\cal K}}{\partial y_{j}}$, that is i). Point ii) is a consequence of property 4.

4. This follows from Theorem 2 and Corollary 2.

$\square$

7 Numerical experiments

We compare the long term evolution coming from system (9) with the full evolution of equation (1), corresponding to the classical restricted $N$-body problem.

To get the evolution of the planets, we compute a planetary ephemerides database for a time span of 2000 yrs, starting at 57600 MJD with a time step of 0.5 years. The computation is performed using the FORTRAN program orbit9 included in the OrbFit free software111http://adams.dm.unipi.it/orbmaint/orbfit. The planetary evolution at the desired time is obtained from this database by linear interpolation.

Inspired by the classification in [10] we consider two paradigmatic cases, representing the two crossing behaviors discussed in the previous sections. The first case is asteroid (887) Alinda, that is considered in the gravitational field of 5 planets, from Venus to Saturn. This asteroid is in $3:1$ mean motion resonance with Jupiter and we will consider its crossings with the orbit of Mars. The second case deals with the ’Toro’ class: we consider a fictitious asteroid that we call 1685a under the influence of 3 planets: the Earth, Mars and Jupiter. This asteroid crosses the orbit of the Earth, and is in the $5:8$ mean motion resonance with it.

We use the same algorithm as in [7] to compute the solution of system (9). This is a Runge-Kutta-Gauss method evaluating the vector field at intermediate points of the time step. The time step is reduced when the trajectory of the asteroid is close to a planet crossing, in order to get exactly the crossing condition. By Theorems 1-2 we can find two locally Lipschitz-continuous extensions of the vector field from both sides of the singular set $\Sigma$. The difference between the two extended fields is given by Corollary 1 for asteroid 887 (Alinda) and by Corollary 2 for asteroid 1685a. In both cases, we compute the intermediate values of the extended vector field just after the crossing, and then we correct them using Corollary 1 or Corollary 2. We use these corrected values as an approximation of the vector field at the intermediate point of the solution, see Figure 2. This algorithm avoids the computation of the vector field at the singular points, which could be affected by numerical instability.

To produce the comparison, we consider 64 possible initial conditions for system (1) corresponding to the same initial condition of system (9). For asteroid 887 (Alinda) these are produced by shifting the mean anomalies in the following way. Let $\bar{\ell}_{j}$ and $\bar{\ell}$ be the mean anomalies of planet $j$ and the asteroid, at the initial epoch 57600 MJD. For each planet, we consider the 64 values $\ell^{(k)}_{j}=\bar{\ell}_{j}+k\pi/64$ with $k=0,\dots,63$. For every $k$, we compute the initial value of the mean anomaly $\ell^{(k)}=\bar{\ell}+l^{(k)}$ of the asteroid such that

[TABLE]

The integration of this 64 different initial conditions is performed with the program orbit9. Then we consider the arithmetic mean of the 5 Keplerian elements $a,e,I,\Omega,\omega$ and the critical angle $\sigma=h^{*}_{5}\ell_{5}+h^{*}\ell$ over these evolutions and compare them with the corresponding elements coming from system (9), in which we choose $n_{\rm max}=3$. Figure 3 summarizes the results: the solid line corresponds to the solution of (9) while the dashed line corresponds to the arithmetic mean of the full numerical integrations. The shaded region represents the standard deviation from the arithmetic mean. The correspondence between the solutions is good. The Mars crossing singularity occurs around $t=3786\>yr$.

For asteroid 1685a we proceed in the same way, with the Earth playing the role of Jupiter. For the long term evolution we used $n_{\rm max}=3,15$. In Figure 4 we show the results. Using $n_{\rm max}=15$ we see that the result improves very much. The Earth crossing singularity occurs around $t=2281\>yr$. In this test the value of $\sigma_{c}$ at crossing results to be about $348$ degrees, which is quite different from all the values of $\sigma$ in Figure 4. We cannot really appreciate the effect of the singularity in the evolution since we obtain very small values of the components $\mathrm{Diff}_{h}(\frac{\partial\mathscr{H}}{\partial y_{i}})$.

8 Conclusions

We studied the long term dynamics of an asteroid under the gravitational influence of the Sun and the solar system planets, assuming that a mean motion resonance between the asteroid and one of the planets occurs. We focused on the case of planet crossing asteroids and considered a resonant normal form $\mathscr{H}_{n_{\rm max}}$, see (7),(8). The analysis is performed separately for crossings with the resonant planet or with another one. In both cases, we could define generalized solutions of the differential equations for the long term dynamics, going beyond the singularity. These solutions are continuous but in general not differentiable. We also proved that generically, in a neighborhood of a crossing time, the evolution of the signed orbit distance along the generalized solutions is more regular that the long term evolution of the orbital elements. In case of crossings with the resonant planet, we recovered the protection mechanism against collisions in the limit $n_{\rm max}\to\infty$. This implies that, if the resonant angle $\sigma$ is different from the critical value $\sigma_{c}$ at the crossing times $t_{c}$ (see Sections 5,6) also deep close encounters are avoided, which makes the results of this theory more reliable. Indeed, close encounters can still occur with a planet not involved in the resonance, and this represent a critical case. Actually, in this case the semimajor axis usually suffer a drastic change [12], pushing the asteroid outside the considered resonance. By means of numerical experiments, in some relevant cases, we showed that the model seems to approximate well the full evolution in a statistical sense. We plan to make numerical tests on a large scale, to study different dynamical behaviors of the population of NEAs.

This work extends the results in [7] to the resonant case and gives a unified view of the orbit crossing singularity in case of mean motions resonances with one planet: indeed, comparing the results in Corollaries 1,2 we see how the discontinuity in the derivatives, represented by $\mathrm{Diff}_{h}\frac{\partial\mathscr{H}}{\partial y_{i}}$, vanishes in a weak sense (i.e. in the sense of distributions) for $n_{\rm max}\to\infty$, if $\sigma\neq\sigma_{c}$. Moreover, the resonant normal form introduced in (8) can easily be extended to include more than one resonance, also with different planets, by considering all the harmonics associated to the corresponding resonant module (see [11, Chap.2]).

Appendix A Appendix

From the definition of the approximate distance $\delta_{h}$, we have that

[TABLE]

We summarize below some relevant estimates and results from [7]. In the following, we shall denote by $c_{i}$, $i=1,\dots,14$, some positive constants independent on $\mathcal{E}$. We first recall some Lemmas.

Lemma 1**.**

There exist positive constants $c_{1}$, $c_{2}$ and a neighborhood ${\cal U}$ of $(\mathcal{E}_{c},V_{h}(\mathcal{E}_{c}))$ such that

[TABLE]

holds for $(\cal E,V)$ in ${\cal U}$. Moreover, there exist positive constants $c_{3}$, $c_{4}$ and a neighborhood ${\cal W}$ of $\mathcal{E}_{c}$ such that

[TABLE]

holds for $\mathcal{E}$ in ${\cal W}$ and for every $V\in\mathbb{T}^{2}$.

Lemma 2**.**

Using the coordinate change $\xi={\cal A}_{h}^{1/2}(V-V_{h})$ and then polar coordinates $(\rho,\theta)$, defined by $(\rho\cos\theta,\rho\sin\theta)=\xi$, we have

[TABLE]

with ${\cal B}=\{\xi\in\mathbb{R}^{2}:|\xi|\leq r\}$. The term $-2\pi d_{h}/\sqrt{\det{\cal A}_{h}}$ is not differentiable at ${\cal E}={\cal E}_{c}\in\Sigma$.

Lemma 3**.**

The maps

[TABLE]

can be extended to two different analytic maps $\mathcal{G}_{h}^{+}$, $\mathcal{G}_{h}^{-}$ such that, in $\mathcal{W}$,

[TABLE]

Moreover the following estimates hold, with $\mathcal{U}_{\Sigma}=\{(\mathcal{E},V_{h}(\mathcal{E})):\mathcal{E}\in\Sigma\}$:

• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]
• [TABLE]

Acknowledgments. We are grateful to Alessandro Morbidelli, whose comments induced us to investigate better the relation between this work and the other results present in the literature, as explained in Section 6. The authors acknowledge the support by the Marie Curie Initial Training Network Stardust, FP7-PEOPLE-2012-ITN, Grant Agreement 317185. S.M. also acknowledges financial support from the Spanish Ministry of Economy and Competitiveness, through the ”Severo Ochoa” Programme for Centres of Excellence in R&D (SEV-2015-0554), the project ”Geometric and numerical analysis of dynamical systems and applications to mathematical physics” (MTM2016-76072-P), and the ”Juan de la Cierva-Formación” Programme (FJCI-2015-24917). G.F.G. has been partially supported by the University of Pisa via grant PRA-2017 ‘Sistemi dinamici in analisi, geometria, logica e meccanica celeste’.