# Detection of close encounters and resonances in three-body problems through Levi-Civita regularization

E. Lega, M. Guzzo, Cl. Froeschlé

► To cite this version:

E. Lega, M. Guzzo, Cl. Froeschlé. Detection of close encounters and resonances in three-body problems through Levi-Civita regularization. Monthly Notices of the Royal Astronomical Society, 2011, pp.1353. 10.1111/j.1365-2966.2011.19467.x. [hal-00632663]

HAL Id: hal-00632663

https://hal.science/hal-00632663v1

Submitted on 17 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Detection of close encounters and resonances in three-body problems through Levi-Civita regularization

E. Lega,<fnref n="1" />* M. Guzzo<fnref n="2" />* and C. Froeschlé<fnref n="1" />*

Accepted 2011 July 18. Received 2011 May 25; in original form 2011 April 20

## ABSTRACT

The detection and identification of resonances is a key ingredient in the studies of stability of dynamical systems, with important applications for our Solar as well as for extrasolar systems. In this paper we study the detection of resonances and close encounters in the three-body problem using the fast Lyapunov indicator method. In order to investigate the close encounters, we needed to adapt the method to the model and its singularities. Our technical improvement lies in computation of the solutions of the singular variational equations by measuring the divergence of close initial conditions. We have used the Levi-Civita regularization for the integration of the equations of motion. As an application of the method, we show that it provides a correct detection of the tube manifolds related to the Lagrangian point $L_1$ of the Sun–Jupiter system.

**Key words:** chaos – methods: analytical – methods: numerical – space vehicles – celestial mechanics.

## 1 INTRODUCTION

Since the pioneering work of Hénon and Heiles (Hénon & Heiles 1964) wherein they studied the phase space with the method of Poincaré surface of section, global phase-space studies of dynamical systems have become a standard approach in many different problems. In the field of astronomy, global studies for the puzzling problem of the long-term stability of our Solar system can be found in, for example Nesvorny & Morbidelli (1998), Murray & Holman (1999), Robutel & Laskar (2001), Morbidelli (2002), Guzzo (2005, 2006), Robutel & Gabern (2006) and Hayes (2008); a picture of the global dynamics in galactic potential is provided by, for example Papaphilippou & Laskar (1998), Voglis, Tsoutsis & Efthymiopoulos (2006) and Namouni, Guzzo & Lega (2008); and studies of global dynamics related to cometary motion, close encounters and space mission design can be found in, for example Villac (2008), Koon et al. (2001), Perozzi & Ferraz Mello (2010), Gomez et al. (2004), Valsecchi (2005) and Koon et al. (2008). Actually, many of these studies have been motivated by the celebrated KAM (Kolmogorov 1954; Arnold 1963; Moser 1958) and Nekhoroshev (Nekhoroshev 1977) theorems, providing fundamental information about the long-term stability of a Hamiltonian system from the global knowledge of phase space, specifically from the distribution of resonances.

Many of these results have been obtained thanks to the development of tools essentially based on the one hand on Fourier analysis (such as Laskar's frequency analysis, Laskar 1990), and on the other

hand on the developments in the Lyapunov exponents theory such as fast Lyapunov indicator (FLI hereafter, Froeschlé, Lega & Gonczi 1997; Guzzo, Lega & Froeschlé 2002) [alternative methods are the Mean Exponential Growth of Nearby Orbits (MEGNO), helicity angles, local Lyapunov characteristic numbers, etc.; the reader can refer to Morbidelli (2002) for an exhaustive review].

The implementation of methods that use the Lyapunov exponents theory requires numerical computation of the solutions of the variational equations, which is very difficult to achieve in the case of close encounters with the massive bodies. In fact, the power of singular terms in the equation of motions increases when computing the variational equations (see Section 2.2).

In this paper we propose a numerical computation of the solutions of the variational equations that works independently of the close encounter with the primary or with the secondary body. The computation is performed by implementing in the three-body problem the method of measure of the divergence of close initial conditions (Benettin, Galgani & Strelcyn 1976), in combination with the Levi-Civita (LC hereafter) regularization with respect to the primary or secondary body. In addition we introduce a technical improvement to the method originally proposed in Benettin et al. (1976), which improves the order of precision of the computation.

In order to avoid the problem of switching between the primary and secondary body, global transformations would provide a simultaneous regularization of the two singularities. One example is provided by the Birkhoff transformation (Birkhoff 1950) suited to the planar circular restricted three-body problem, and other global transformations can be found in Szebehely (1967). The extension to the spatial case is reviewed in Della Penna (2002) and

<footnote n="1">¹Université de Nice Sophia Antipolis, CNRS UMR 6202, Observatoire de la Côte d'Azur, Bv. de l'Observatoire, BP 4229, 06304 Nice cedex 4, France</footnote>

<footnote n="2">²Università degli Studi di Padova, Dipartimento di Matematica Pura ed Applicata via Trieste 63, 35121 Padova, Italy</footnote>

<footnote>*E-mail: elena@oca.eu (EL); guzzo@math.unipd.it (MG); claude@oca.eu (CF)</footnote>
references therein. In addition to the problem of their complexity, global transformations cannot be extended to the N-body problem. Other different regularizations can be found, for example in Froeschlé (1970) and Aarseth (1974), for integrating orbits of the restricted three-body problem in the spatial case, and later generalized to the N-body problem by Heggie (1974). The reader can find a review in Aarseth (2001).

The use of the regularizing transformation in both the equations of motion and the variational equations is quite complicated in the three-body problem.

First we illustrate the method with a test computation of the largest Lyapunov exponent of individual orbits, then we provide a significant application of the FLI method for the detection of the so-called tube manifolds of the Lagrangian points $L_1$ and $L_2$ for the Sun-Jupiter system. We recall that recent applications of the FLI method (Villac 2008; Guzzo, Lega & Froeschlé 2009; Guzzo 2010; Lega, Guzzo & Froeschlé 2010) provided a new way of computation of the so-called stable and unstable manifolds related to the resonances. Investigation of these manifolds is, since Poincaré, the key point for understanding chaos and diffusion in dynamical systems. In the case of the restricted circular three-body problem (RC3BP hereafter), the computation of a piece of the unstable manifolds associated with the Lyapunov orbits of the Lagrangian points $L_1$ and $L_2$, also known as tube manifolds, allows us to show, for example, the role played by heteroclinic intersections on cometary orbits (Koon et al. 2001). A renewed interest in this is related also to the spacecraft trajectory design (see Simó 1995 and references therein; Gomez et al. 2004; Koon et al. 2008; Villac 2008 for an application of FLI).

The paper is organized as follows. In Section 2 we discuss the problem of numerical integration of the variational equations of the restricted three-body problem; we show that the problem is not trivial even in the regularized setting and we propose a method of computation based on the divergence of nearby orbits. Examples of computation of the largest Lyapunov exponent and of the tube manifolds are provided in Section 3. Conclusions are given in Section 4. We discuss in Appendix A the LC regularization in the three-body problem.

## 2 THE RESTRICTED THREE-BODY PROBLEM AND ITS VARIATIONAL EQUATIONS

RC3BP describes the motion of a massless body P perturbed gravitationally by two massive bodies $P_1$ and $P_2$ (called primary and secondary body, respectively). In the rotating frame $xOy$, the equations of motion of P are

$$ \begin{cases} \ddot{x} = 2\dot{y} + x - (1 - \mu)\frac{x+\mu}{r_1^3} - \mu\frac{x-1+\mu}{r_2^3} \\ \ddot{y} = -2\dot{x} + y - (1 - \mu)\frac{y}{r_1^3} - \mu\frac{y}{r_2^3}, \end{cases} \quad (1) $$

where the units of mass, length and time have been chosen so that the masses of $P_1$ and $P_2$ are $1 - \mu$ and $\mu$ ($\mu \le 1/2$), respectively; their coordinates are $(-\mu, 0)$ and $(1 - \mu, 0)$; and their revolution period is $2\pi$. We define $r_1^2 = (x + \mu)^2 + y^2$ and $r_2^2 = (x - 1 + \mu)^2 + y^2$.

As it is well known, the problem has an integral of motion C, the so-called Jacobi integral, defined by

$$ C = x^2 + y^2 + 2\frac{1-\mu}{r_1} + 2\frac{\mu}{r_2} - \dot{x}^2 - \dot{y}^2. \quad (2) $$

In order to introduce the variational equations of (1), we first rewrite it as a first-order system of equations:

$$ \left\{ \begin{aligned} \dot{x} &= v_x \\ \dot{y} &= v_y \\ \dot{v}_x &= 2v_y + x - (1 - \mu) \frac{x+\mu}{r_1^3} - \mu \frac{x-1+\mu}{r_2^3} \\ \dot{v}_y &= -2v_x + y - (1 - \mu) \frac{y}{r_1^3} - \mu \frac{y}{r_2^3}, \end{aligned} \right. \quad (3) $$

and then we introduce its compact form:

$$ \dot{\xi} = F(\xi), \quad (4) $$

where $\xi = (x, y, v_x, v_y)$. The variational equations of (4) are therefore

$$ \left\{ \begin{aligned} \dot{\xi} &= F(\xi) \\ \dot{w} &= \frac{\partial F}{\partial \xi} w, \end{aligned} \right. \quad (5) $$

where $w \in \mathbb{R}^4$ represents the tangent vector. As it is evident, the power of singular terms of the vector field $F$ increases when computing the Jacobian matrix $\frac{\partial F}{\partial \xi}$, and the numerical integration of (5) becomes very difficult to achieve in the case of close encounters.

To solve the problem related to the singularities in the equations of motion (4), different theories of regularization of the Kepler problem have been developed in the last century. However, these regularizations do not resolve the worst singularities of equations (5). We show this fact in Section 2.2, which can also be understood with the simpler one-dimensional regularization.

Since the numerical integration of the variational equations (5) is a hard task, we use the method of computation based on the divergence of nearby orbits, introduced in Benettin & Galgani (1979), as an alternative method of computation of the largest Lyapunov exponent.

### 2.1 The method of divergence of nearby orbits

The first efficient numerical methods for the computations of the largest Lyapunov exponent were introduced in the seventies (Benettin et al. 1976; Benettin & Galgani 1979). In Benettin et al. (1976) the computation was implemented by means of the divergence of nearby orbits. Precisely, let us consider a differential equation

$$ \dot{\xi} = F(\xi), $$

and let us denote by $\phi(t, \xi)$ its flow. For any fixed time $t$, the solution $\xi(t)$, $w(t)$ of the variational equations

$$ \left\{ \begin{aligned} \dot{\xi} &= F(\xi) \\ \dot{w} &= \frac{\partial F}{\partial \xi} w \end{aligned} \right. \quad (6) $$

satisfies

$$ w(t) = \phi(t, \xi(0) + w(0)) - \phi(t, \xi(0)) + \mathcal{O}(\|w(0)\|^2). \quad (7) $$

Therefore, as the norm of $w(0)$ becomes small, the vector $w(t)$ is well approximated by the divergence of orbits with nearby initial conditions $\xi(0) + w(0)$ and $\xi(0)$.

For numerical integrations, for a fixed small time $\tau$, the flow is approximated by an integrator map $\psi(\xi)$, and the evolution of the tangent vector $w(i\tau)$ is obtained by iterating the map:

$$ \left\{ \begin{aligned} \xi((i+1)\tau) &= \psi(\xi(i\tau)) \\ w((i+1)t) &= \psi(\xi(i\tau) + w(i\tau)) - \psi(\xi(i\tau)). \end{aligned} \right. \quad (8) $$

Since the tangent dynamics is linear with respect to the tangent vectors, it is always possible to implement the computation with
normalized vectors,

$$ \tilde{w}(0) = s_0 \frac{w(0)}{\|w(0)\|}, $$

with $s_0$ being suitably small. One chooses $s_0$ so that the second-order terms in (7) are very small, ideally smaller than the required precision of the computation. Since the dynamics may amplify the norm of the vectors, one typically needs to normalize several times the vector during the computation in order to achieve the desired precision. This is the essence of the method described in Benetin et al. (1976). The original method was abandoned afterwards; instead it was modified by using a direct computation of the solutions of the variational equations, with normalizations. We recover in this paper the original method, which we find more suitable for the problem of close encounters, with some modifications that improve the numerical precision. In fact, since the second-order term in equation (7) is symmetric with respect to the initial vector $w(0)$, we gain an order in the precision of the computation by using the divergence of three nearby orbits and by approximating $w((i+1)t)$ with

$$ \left\{ \begin{array}{l} \xi((i+1)\tau) = \psi(\xi(i\tau)) \\ w((i+1)t) = \frac{1}{2}(\psi(\xi(i\tau) + w(i\tau)) - \psi(\xi(i\tau) - w(i\tau))). \end{array} \right. \quad (9) $$

Equation (7) now improves to

$$ w(t) = \frac{1}{2}(\phi(t, \xi(0) + w(0)) - \phi(t, \xi(0) - w(0))) + \mathcal{O}(\|w(0)\|^3). \quad (10) $$

We list below the steps used for the implementation of the method in the computations done in this paper.

(i) Let us consider an initial condition $\xi(0)$, $w(0)$. We first renormalize $w(0)$ to $\tilde{w}(0) = s_0 w(0)/\|w(0)\|$, with suitable $s_0$.

(ii) We iterate the algorithm (9) until the norm $\tilde{w}(t)$ is bigger than a suitably chosen $s_1$ (for example, one may set $s_1 = 10s_0$), or until $|t|$ is bigger than a fixed amount of time $T$. Therefore, we renormalize the vector $w(t)$ as in step (i).

The algorithm produces at any time $t$ a vector $w_t$ which is related to the tangent vector $w(t)$ through the relation

$$ w(t) = \frac{w_t}{s_0^N}, $$

where $N$ is the number of normalizations performed during the integration in the time interval $[0, t]$. We also chose the initial vector tangent to the surface of constant Jacobi integral $C$. The numerical integrator $\psi$ is a fourth-order Runge-Kutta applied to the equations of the three-body problem. Such equations are written using the LC transformation with respect either to the primary or to the secondary, depending on the relative position of the body (see Appendix A).

## 2.2 The variational equations and regularizing transformations

In this subsection we show that the regularizing transformations such as the LC and the Kustaanheimo-Stiefel (Kustaanheimo & Stiefel 1965) do not resolve the singularities in the variational equations due to the divergence of the gravitational potential at collisions. We consider here for simplicity the one-dimensional case, corresponding to bodies moving under their gravitational attraction on a straight line. The same kind of problems are encountered in the planar case with the LC transformation and in the spatial case with the Kustaanheimo-Stiefel transformation.

Let us consider the Lagrangian system

$$ L(x, \dot{x}) = \frac{1}{2} \dot{x}^2 + \frac{1}{x} \quad (11) $$

with $x \in \mathbb{R} \setminus \{0\}$. The variational equations are

$$ \left\{ \begin{array}{l} \dot{x} = v \\ \dot{v} = -\frac{1}{x^2} \\ \dot{w}_x = w_v \\ \dot{w}_v = \frac{2}{x^3} w_x, \end{array} \right. \quad (12) $$

where $w = (w_x, w_v)$ represents the tangent vector.

We now perform the one-dimensional version of LC transformation, and write the equations of motion and the variational equations in the regularized variables. Setting $x = u^2$, and correspondingly $\dot{x} = 2u\dot{u}$, the Lagrangian transforms to

$$ \tilde{L}(u, \dot{u}) = 2u^2 \dot{u}^2 + \frac{1}{u^2}. \quad (13) $$

The regularizing transformation is complemented by the introduction of a fictitious time $s$:

$$ dt = u^2 ds. \quad (14) $$

As a result of these transformations, the motions on any surface of constant energy $E$ satisfy the non-singular equation

$$ u'' = -\frac{1}{2}Eu, $$

$$ \text{where } u'' = \frac{d^2}{ds^2} u. $$

We now check if the transformation regularizes also the variational equations (12). We first perform the change of variables:

$$ (u, p) \rightarrow (x, v) = \left( u^2, \frac{p}{2u} \right), $$

which implies the transformation of the tangent vectors:

$$ (w_u, w_p) \rightarrow (w_x, w_v) = \left( 2u w_u, \frac{w_p}{2u} - \frac{p}{2u^2} w_u \right). $$

The above transformation conjugates, on any surface of constant energy $E$, the variational equations (12) to the following equations:

$$ \left\{ \begin{array}{l} \dot{w}_u = \frac{w_p}{4u^2} - \frac{p}{2u^3} w_u \\ \dot{w}_p = -\frac{2E}{u^3} w_u. \end{array} \right. \quad (15) $$

After the introduction of the fictitious time (14), the above equations become

$$ \left\{ \begin{array}{l} w'_{u} = \frac{w_{p}}{4} - \frac{p}{2u} w_{u} \\ w'_{p} = -2E w_{u}. \end{array} \right. \quad (16) $$

The first equation is still singular because of the term $1/u$. Therefore, the LC-like transformation has not resolved the singularity of the variational equations.

The variational equations of the regularized equation of motions,

$$ \left\{ \begin{array}{l} w'_{u} = \frac{w_{p}}{4} \\ w'_{p} = -2E w_{u}, \end{array} \right. \quad (17) $$

are different from (16), and do not correspond to the variational equations (12) transformed with the regularizing transformation. The reason is that the LC transformation includes a fictitious time change which is a function of the spatial coordinates.

## 3 APPLICATIONS TO THE RESTRICTED CIRCULAR THREE-BODY PROBLEM

In this section we report on the applications of the method of computation explained in Section 2 to RC3BP for the Sun-Jupiter system
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Computation of the largest Lyapunov exponent for a regular (left-hand panel) and a chaotic orbit (right-hand panel) with initial conditions respectively: $x = 0.99$, $\dot{x} = y = 0$, $\dot{y} = \dot{y}_C$ and $x = -1.9$, $\dot{x} = y = 0$, $\dot{y} = \dot{y}_C$. $\dot{y}_C$ is the value of the initial velocity $\dot{y}$ as obtained from the Jacobi constant $C = 3.03$. The top panels show the projection of the two orbits on the $(x, y)$ plane (for a short time). The bottom panels show the computation of the largest Lyapunov indicator, i.e. of the log $||w(t)||$ divided by $t$. In the case of the regular orbit the largest Lyapunov indicator converges to zero. For the other orbit, which performs several transitions between the regions which are inner and outer with respect to Jupiter, the evolution of the Lyapunov indicator is compatible with a chaotic orbit.</figcaption></figure>

mass ratio and the Jacoby constant $C = 3.03$. As a test of the method we first compute the largest Lyapunov indicator, whose limit for $t \to \infty$ is the largest Lyapunov exponent, for orbits having close encounters with Jupiter, regular and chaotic ones. The results are presented in Fig. 1: the top panels present the projection of the orbits in the $x, y$ plane (for a short time interval) and the bottom panels present the evolution of the largest Lyapunov indicator on a logarithmic scale. In the case of the regular orbit, the largest Lyapunov indicator converges to zero as expected. We remark that this computation represents a test of the precision of the algorithm because the integrator switches the regularization from the primary to the secondary many times during the numerical integration. The second orbit performs several transitions between regions that are inner and outer with respect to Jupiter's orbit, and the evolution of the Lyapunov indicator is compatible with a chaotic orbit.

As a further example we provide in Fig. 2 the computation of log $||w(t)||$ as a function of time for a regular orbit having close encounters with Jupiter. It appears clearly that approaching the singularity the tangent vector increases sharply and decreases when

the orbit moves away from Jupiter. The growth of the tangent vector is compatible with a logarithmic growth, typical of regular orbits.

We provide a significant application of the FLI method for the detection of the so-called tube manifolds of the Lagrangian points $L_1$ and $L_2$ for the Sun-Jupiter system. FLI is defined as usual (Guzzo et al. 2002):

$$ \text{FLI}(\xi(0), w(0), t) = \sup_{\tau \le t} \log ||w(\tau)||, \qquad (18) $$

where $w(t)$ is the tangent vector at time $t$ for an orbit with initial conditions $(x(0), w(0))$. As we have shown in several of our previous papers, FLI distinguishes the dynamical behaviour of the orbit with initial condition $x(0)$ in a relatively short time $t$ (Froeschlé, Guzzo & Lega 2000; Guzzo, Lega & Froeschlé 2002).

Therefore, the computation of FLI on two-dimensional grids of initial conditions provides a picture of the phase space (Froeschlé et al. 2000). Recent applications of the FLI method (Villac 2008; Guzzo et al. 2009; Guzzo 2010; Lega et al. 2010) have also provided
<figure><img src="image_3.png" /><figcaption>**Figure 4.** Computation of FLI for a grid of 1000 × 1000 initial conditions regularly spaced in *x*(0), *y*(0). The other initial conditions are *y*(0) = 0 and *y'(0) = y'_c*; *y'_c* is obtained from the Jacobi constant *C* = 3.03. The integration time is *t* = 50. Therefore, we expect to detect additional lobes of the hyperbolic manifold with respect to the computation presented in Fig. 3. Again, for each initial condition we present the value of FLI using the colour scale shown below the figure: the largest values of FLI, denoted by light grey, correspond to an hyperbolic manifold. Actually, with respect to the shorter time computation of Fig. 3, we detect additional lobes belonging to a longer piece of manifold.</figcaption></figure>

4 we show the result of computation of FLI on a grid of 1000 × 1000 initial conditions on the same portion of the phase space as in fig. 7 of Vela-Arevalo & Marsden (2004). Precisely, the initial conditions are regularly spaced on *x*, *x*' with −2.2 ≤ *x*(0) ≤ −1.1 and −0.3 ≤ *x'(0)* ≤ 0.3; *y*(0) = 0 and *y'(0)* is obtained from the Jacobi constant *C* = 3.03. The integration time is *t* = 15 in Fig. 3 and *t* = 50 in Fig. 4. The value of FLI is indicated in the panels using the colour scale shown below them: the largest values of FLI, denoted by light grey,<fnref n="1" /> correspond to the hyperbolic manifolds. In Fig. 3, the points of maximum FLI evidently correspond to the intersection of the tube manifold with the section (see top left of fig. 7 of Vela-Arevalo & Marsden 2004 for comparison). Fig. 4 shows the same computation done on a longer time *t* = 50, revealing additional lobes belonging to a longer piece of the manifold.

## 4 CONCLUSION

In this paper we have proposed a method for the numerical computation of the variational equations of RC3BP, which is particularly suitable for the computation of dynamical indicators. Dynamical indicators provide a detailed knowledge of the phase-space structure, such as resonances and invariant manifolds. The method enriches the possibility of the numerical investigation of the three-body dynamics, as for example the computation of tube manifolds. This result is relevant for many astronomical studies of theoretical interest; it may also be relevant for applications such as mission design.

## ACKNOWLEDGMENTS

MG has been supported by the project CPDA092941/09 of the University of Padova. Part of the computations have been done on the ‘Mesocentre SIGAMM’ machine, hosted by the Observatoire de la Cote d’Azur.

## REFERENCES

* Aarseth S. J., 1974, Celest. Mech. Dynamical Astron., 10, 185

* Aarseth S. J., 2001, in Steves B., Maciejewski A., eds, The Restless Universe. Applications for Gravitational N-body Dynamics to Planetary, Stellar and Galactic Systems. SUSSP Publ. Inst. Phys. Publishing, Bristol, p. 93

* Arnold V. I., 1963, Russian Math. Surv., 18, 9

* Benettin G., Galgani L., 1979, Lyapunov Characteristic Exponents and Stochasticity. Intrinsic Stochasticity in Plasmas. Les Editions de Physique, Orsay, France, p. 93

* Benettin G., Galgani L., Strelcyn J. M., 1976, Phys. Rev. A, 14, 2338

* Birkhoff G. D., 1950, in Collected Math. Papers Vol. 1, The Restricted Problem of Three Bodies. Am. Math. Soc., New York

* Della Penna M. G., 2002, in Benest D., Foreschlé C., eds, Lect. Notes Phys. Vol. 590, Singularities in Gravitational Systems. Springer, Berlin, p. 49

* Froeschlé C., 1970, A&A, 4, 115

* Froeschlé C., Lega E., Gonczi R., 1997, Celest. Mech. Dynamical Astron., 67, 41

* Froeschlé C., Guzzo M., Lega E., 2000, Sci, 289, 2108

* Gómez G., Koon W. S., Lo M. W., Marsden J. E., Masdemont J., Ross S. D., 2004, Nonlinearity, 17, 1571

* Guzzo M., 2005, Icarus, 174, 273

* Guzzo M., 2006, Icarus, 181, 475

* Guzzo M., 2010, in Perozzi E., Mello F., eds, Space Manifold Dynamics: Novel Spaceways for Science and Exploration. Springer, Berlin, p. 97

<footnote n="1">¹ The colour versions of all the figures can be found in the electronic version of the paper; the light grey here corresponds to yellow there and the darker grey here corresponds to red–violet there.</footnote>
Guzzo M., Lega E., Froeschlé C., 2002, Phys. D, 163, 1

Guzzo M., Lega E., Froeschlé C., 2009, Phys. D, 182, 1797

Hayes W., 2008, MNRAS, 386, 295

Heggie D. C., 1974, Celest. Mech. Dynamical Astron., 10, 217

Hénon M., Heiles C., 1964, AJ, 69, 73

Kolmogorov A. N., 1954, Doklady Akademii Nauk SSSR, 98, 527

Koon W. S., Lo M. W., Marsden J. E., Ross S. D., 2001, Celest. Mech. Dynamical Astron., 81, 27

Koon W. S., Lo M. W., Marsden J. E., Ross S. D., 2008, Dynamical Systems, the Three Body Problem and Space Mission Design. Marsden Books, Wellington, New Zealand

Kustaanheimo P., Stiefel E., 1965, J. für reine angewandte Math., 218, 204

Laskar J., 1990, Icarus, 88, 266

Lega E., Guzzo M., Froeschlé C., 2010, Celest. Mech. Dynamical Astron., 107, 115

Levi-Civita T., 1906, Acta Math., 30, 305

Morbidelli A., 2002, Modern Celestial Mechanics: Aspects of Solar System Dynamics. Taylor & Francis, Abingdon, UK

Moser J., 1958, Communications Pure Applied Math., 11, 81

Murray N., Holman M., 1999, Sci, 283, 1877

Namouni F., Guzzo M., Lega E., 2008, A&A, 489, 1363

Nekhoroshev N. N., 1977, Russian Math. Surv., 32, 1

Nesvorny D., Morbidelli A., 1998, AJ, 116, 3029

Papaphilippou Y., Laskar J., 1998, A&A, 329, 451

Robutel P., Gabern F., 2006, MNRAS, 372, 1463

Robutel P., Laskar J., 2001, Icarus, 152, 4

Simó C., 1995, in C. Simó, ed., NATO Advanced Sci. Inst. Ser. C Vol. 533 Math. Phys. Sci., Hamiltonian Systems with Three or more Degrees of Freedom. Kluwer, Dordrecht, p. 223

Perozzi E., Ferraz-Mello S., eds, 2010, Space Manifold Dynamics: Novel Spaceways for Science and Exploration. Springer, Berlin

Szebehely V., 1967, Theory of Orbits. Academic Press, New York

Valsecchi G. B., 2005, Comptes Rendus – Phys., 6, 337

Vela Arevalo L., Marsden J., 2004, Class. Quantum Grav., 21, 361

Villac B. F., 2008, Celest. Mech. Dynamical Astron., 102, 29

Voglis N., Tsoutsis P., Efthymiopoulos C., 2006, MNRAS, 373, 280

## APPENDIX A: THE LEVI-CIVITA REGULARIZATION AND THE THREE-BODY PROBLEM

Since the potential function is infinite at collisions, the numerical investigation of close encounters with both the massive bodies needs special attention. To this aim, various theories of regularization of the Kepler problem have been developed in the last century. In this paper we use two Levi-Civita regularizations (Levi-Civita 1906) with respect to the primary or secondary body, in order to treat both the singularities.

The Levi-Civita transformation is

$$ \begin{pmatrix} x - x_0 \\ y \end{pmatrix} = \begin{pmatrix} u_1 & -u_2 \\ u_2 & u_1 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix} = \mathbf{A} \begin{pmatrix} u_1 \\ u_2 \end{pmatrix}, \quad (A1) $$

where $x_0$ is a parameter that can be set equal to the abscissa of $P_1$ or $P_2$.

The transformation will eliminate only one of the two singularities. We will set $x_0$ depending on whether the third body is closer to the primary or to the secondary. Precisely, the choice is made by comparing the forces, i.e. we regularize with respect to $P_1$ if

$$ \frac{(1 - \mu)}{r_1^2} > \frac{\mu}{r_2^2}, $$

otherwise we regularize with respect to $P_2$.

The LC transformation is completed by the introduction of a fictitious, or regularized, time $s$ defined by

$$ dt = r \, ds, \quad (A2) $$

where $r = u_1^2 + u_2^2 = \sqrt{(x - x_0)^2 + y^2}$. The equations of motion in the regularized coordinates are (see, for example Szebehely 1967)

$$ \left\{ \begin{aligned} u_1'' &= \frac{1}{4}[(a+b)u_1 + cu_2] \\ u_2'' &= \frac{1}{4}[(a-b)u_2 + cu_1] \end{aligned} \right. \quad (A3) $$

with

$$ \left\{ \begin{aligned} a &= \frac{2m}{R} - C + x^2 + y^2 \\ b &= 4y' + 2rx - \frac{2mr((x-x_0)\pm1)}{R^3} \\ c &= 2ry - 4x' - \frac{2mry}{R^3}. \end{aligned} \right. \quad (A4) $$

In (A4) we use $R = r_2$, the sign plus, and $m = \mu$ if we regularize with respect to $P_1$ and $R = r_1$, the sign minus, and $m = 1 - \mu$ if we regularize with respect to $P_2$.

Let us remark that the denominators appearing in (A4) are well separated from zero, therefore the singularity associated either with $1/r_1$ or with $1/r_2$ is removed.

As usual the integration of second-order differential equations is done by splitting each equation into two first-order differential equations. We have therefore a set of four first-order differential equations. Starting with an initial condition $(x, y)$, we obtain the initial values of $(u_1, u_2)$ as follows:

$$ \left\{ \begin{aligned} x &< x_0 \\ u_2 &= \sqrt{\frac{r-(x-x_0)}{2}} & \quad u_1 &= \sqrt{\frac{r+(x-x_0)}{2}} \\ u_1 &= \frac{y}{2u_2} & \quad u_2 &= \frac{y}{2u_1}. \end{aligned} \right. \quad (A5) $$

We have distinguished the cases $x < x_0$ and $x \ge x_0$ in order to avoid zero denominators in the computations of respectively $u_1$ and $u_2$. Concerning the velocities, the relation between the physical and regularized velocities is

$$ \begin{aligned} x' &= 2(u_1 u_1' - u_2 u_2'), \\ y' &= 2(u_1 u_2' + u_2 u_1'), \end{aligned} \quad (A6) $$

where $'$ denotes a derivative with respect to $s$. The first derivatives of the coordinates $u$ with respect to $s$ as a function of $x'$ are

$$ \begin{pmatrix} u_1' \\ u_2' \end{pmatrix} = \frac{1}{2r} \mathbf{A}^T \begin{pmatrix} x' \\ y' \end{pmatrix}, \quad (A7) $$

where $\mathbf{A}^T$ is the transpose of the matrix $\mathbf{A}$. Actually one can write the previous expression with respect to $x'$ or to $\dot{x}$ taking into account that

$$ \begin{pmatrix} x' \\ y' \end{pmatrix} = r \begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix}. \quad (A8) $$

Equations (A1) and (A6), (A8) provide the position and velocities in the physical space as a function of the regularized coordinates while equations (A5) and (A7), (A8) provide the regularized position and velocities as a function of the physical coordinates.

This paper has been typeset from a TeX/LA TeX file prepared by the author.