A Kustaanheimo-Stiefel regularization of the elliptic restricted three-body problem and the detection of close encounters with fast Lyapunov indicators

We present the Kustaanheimo-Stiefel (KS) regularization of the elliptic restricted three-body problem (ER3BP) at the secondary body $P_2$, and discuss its use to study a category of transits through its Hill's sphere (fast close encounters). Starting from the Hamiltonian representation of the problem using the synodic rotating-pulsating reference frame and the true anomaly of $P_2$ as independent variable, we perform the regularization at the secondary body analogous to the circular case by applying the classical KS transformation and the iso-energetic reduction in an extended 10-dimensional phase-space. Using such regularized Hamiltonian we recover a definition of fast close encounters in the ER3BP for small values of the mass parameter $\mu$ (while we do not require a smallness condition on the eccentricity of the primaries), and we show that for these encounters the solutions of the variational equations are characterized by an exponential growth during the fast transits through the Hill's sphere. Thus, for small $\mu$, we justify the effectiveness of the regularized fast Lyapunov indicators (RFLIs) to detect orbits with multiple fast close encounters. Finally, we provide numerical demonstrations and show the benefits of the regularization in terms of the computational cost.


Introduction
The regularization of the gravitational singularities, appeared at the beginning of the XXth century, has become in the last decades an extremely useful technique to deal with the numerical integration of the N-body problems.Particularly two kinds of regularization techniques are widely known: a geometric one, which basically aims at modifying the equations of motion such that they are defined and regular even on the singularities, and a solution-based one, whose goal consists in an analytic continuation of the original solution through the singular point.In this paper we focus precisely on a celebrated example of the former category, i.e., the Kustaanheimo-Stiefel regularization for a special case of utmost importance, represented by the restricted three-body problem, originally in its circular (CR3BP) and then elliptic (ER3BP) variant.
In a seminal paper [27] Levi-Civita performed a local1 regularization of the planar CR3BP, which relies on the conservation of the so-called Jacobi integral, through the introduction of canonical transformations and a time reparametrization that nowadays are known, after his name, as Levi-Civita (LC) regularization.The issue for the spatial CR3BP was solved by Kustaanheimo and Stiefel in the mid-1990s [22,23].The latter procedure is more complicated than the LC one since it exploits a projection map from a space of four redundant variables to the three-dimensional Cartesian space.Both LC and KS regularizations are iso-energetic, since they exploit the existence of a global first integral, the so-called Jacobi constant.There exist in the literature many uses of the KS transformation to regularize binary collisions in the general 3-body and N-body problems, as well as perturbations of the Kepler problem whose definitions include restricted problems more general than the CR3BP and ER3BP (see [40] and, for example, [1,20,39,9,2,46,24,5] and references therein).For the ER3BP (and for more general restricted problems) a 10-dimensional phase-space is required.In fact, in addition to the 8-dimensional phase-space of the KS variables and an additional variable corresponding to the physical time, another variable corresponding to the energy of the system (variable in time for the ER3BP) is required.While these approaches include the regularizations of the ER3BP, the geometric properties of this system provide specific representations which are more adapted, for example, to the computation of ejection-collision orbits or the dynamics at the Lagrangian points (see, e.g., [42,33]) and, as we consider in this paper, the study of fast close encounters with the secondary body P 2 .
We remark that specific approaches to the regularization of the ER3BP have been developed also using regularizations different from the KS one.Formulations in this regard have been derived in [42,41] and applied in [6,34] for the planar setting; in [44,45,3,28] for the spatial one.
A convenient formulation of the ER3BP uses a synodic rotating-pulsating reference frame and the true anomaly of the secondary body as independent variable [42].Consider an ER3BP defined by the motion of a body P of negligible mass in the gravitation field of two massive bodies P 1 (the primary) and P 2 (the secondary) performing an elliptic Keplerian motion of eccentricity ε ∈ (0, 1).As usual, the simplifying assumptions on the units correspond to setting m 1 = 1 − µ, m 2 = µ for µ ∈ (0, 1/2) as the masses of P 1 , P 2 respectively, while a = 1 and T = 2π are the semi-major axis and the period of the elliptic motion.By denoting with (x, y, z), (p 1 , p 2 , p 3 ) the coordinates of P and their conjugate momenta, and with f the true anomaly of the elliptic motion of P 2 which is used as independent variable, the Hamiltonian reads: where d 1 = ∥P −P 1 ∥, d 2 = ∥P −P 2 ∥.In this paper we first represent the KS regularization at the secondary body of the Hamiltonian (1) by applying the classical KS transformation and the iso-energetic reduction in the extended 10-dimensional phase-space.First, we provide indeed a simple and self-consistent proof on the projection of the solutions of the regularized Hamiltonian to the original solutions, by adapting to the elliptic case a derivation of the KS Hamiltonian of the CR3BP given in [7].Then, we also assess the effectiveness of the regularization for numerical integrations, and the advantage in terms of numerical performances, in a fictitious simple scenario which is nevertheless representative (for the choice of the initial conditions) of realistic close encounters in the Solar System, such as the non-coplanar close encounters with Jupiter, in the Sun-Jupiter ER3BP.Finally, we use the regularized Hamiltonian to extend the definition of fast close encounters to the ER3BP, and justify the effectiveness of the RFLIs to detect orbits with multiple of such fast close encounters.
We call 'close encounter' a transit of a solution (x(f ), y(f ), z(f )) of the ER3BP through the Hill's sphere of P 2 occurring in any finite interval [f 1 , f 2 ] of the true anomaly f (collision solutions are not included).We consider values of the mass parameter µ ≤ 1/10 and we define the Hill's sphere of P 2 by: B(µ Notice that this definition of Hill's sphere is different from the conventional one by a numerical factor: the radius µ with the exclusion of a neighborhood of γ = 0.The condition (2) (or similar ones), which appeared in several studies of close encounters with Levi-Civita regularization (see, e.g., [21,10,13,12,18]) as well as in the heuristic approach known as Öpik's theory (see [32], revised in [43]), is understood by representing the Hamiltonian of the planar CR3BP, Levi-Civita regularized at P 2 , in the form: where (u 1 , u 2 ) are the Levi-Civita coordinates and U = (U 1 , U 2 ) are the momenta conjugate to u = (u 1 , u 2 ) defined as in [27]; ) is regular at u = 0 with Taylor expansion which begins with order 6.If γ > 0, the coefficient of ∥u∥ 2 in the Levi-Civita Hamiltonian (3) is strictly negative, allowing to use methods of hyperbolic dynamics to study the fast close encounters of the CR3BP, as it was done in [21,10,12] with analytic methods, and in [13] to study the effectiveness of RFLIs to detect orbits with multiple close enclounters in the planar CR3BP.We remark that, by considering the higher order corrections to the quadratic approximation of the Hamiltonian, a neighborhood of the limit case γ = 0 should be avoided.Alternative definitions of fast close encounters, as in Öpik's theory, refer to the hyperbolic approximations of the Cartesian solutions which are obtained by considering the Keplerian motion defined by the secondary body P 2 .
Regularized fast Lyapunov indicators have been introduced in [8,25,13] (see also [19], and refeences therein) and have been used in [13,26,17] to detect the several kinds of close encounters with the secondary body P 2 .The ability of FLIs to detect fast close encounters of the planar CR3BP [36,35], i.e., for γ > 0, has been related to the exponential growth of tangent vectors for orbits transiting suitably fast in the Hill's sphere of the planet [13].In fact, the Jacobian matrix X (u, U ) of the Hamiltonian vector field of (3) computed in the limit: where s = s(t) is the Levi-Civita time reparametrization, has eigenvalues ± γ/2, and therefore is hyperbolic when γ > 0.
To develop this idea for the full ER3BP we consider the Hamiltonian K(u, ϕ, U, Φ) of the problem regularized at P 2 , as it will be obtained in Section 2: u = (u 1 , u 2 , u 3 , u 4 ) denote the KS variables; U = (U 1 , U 2 , U 3 , U 4 ) a set of conjugate momenta.The additional conjugate variables ϕ, Φ are needed in the regularization of the elliptic problem: by denoting with s the independent variable (the proper time) of the Hamilton equations of the regularized Hamiltonian, ϕ(s) is the true anomaly of P 2 for the value s of the proper time (see Section 2 for all the details).We obtain for K(u, ϕ, U, Φ) a representation similar to (3): where b(u) = O(∥u∥3 ) ∈ R 4 ; R 6 (u, ϕ) is regular at u = 0, and its Taylor expansion in the vector variable u begins with order 6.There is however a fundamental difference with respect to the circular case, since the coefficient: depends on the variables Φ, ϕ, and therefore is not constant along the solutions (u(s), ϕ(s), U (s), Φ(s)) of the Hamilton equations of K.Moreover, the derivative of Γ s := Γ(ϕ(s), Φ(s)) with respect to the proper time s is proportional to ε, but does not vanish for µ → 0.
Using the representation (4), we notice that the variation of Γ s satisfies: 6 ) , so that, during the transit through the Hill's sphere, where ∥u∥ is smaller than order µ 1 6 (see Section 2), we have: The stability of Γ s up to times of the order of 1/(εµ) provides the opportunity to assess the hyperbolicity of fast close encounters for the ER3BP, for small values of µ.In fact, to establish the hyperbolic character of the close encounter in the ER3BP from the representation (4) we need to establish the stability of the coefficient Γ s for two reasons.First, we notice that the variational matrix X of the Hamiltonian vector field of K has the representation: with Therefore, on one hand we need that motions entering the Hill's sphere with Γ 0 > 0, maintain a value of Γ s > 0 during the transit through the Hill's sphere, so that the matrix X 0 is hyperbolic.On the other hand, we need to provide an upper bound to the elements of the matrix X − X 0 during the transit in the Hill's sphere, to ensure the hyperbolicity of X for suitably small values of µ.In particular, a sufficient upper bound to |U j (s)| is obtained if, for example, during the transit we have Γ s ≤ (3/2)Γ 0 (see Section 3 for all the details).In Section 3 we prove that, if µ satisfies: justify the effectiveness of the RFLIs to detect orbits with multiple fast close encounters also for the ER3BP when the parameter µ is small.We also provide numerical demonstrations of the detection of close encounters with regularized Lyapunov indicators.
The paper is structured as follows.In Section 2 we present a step-by-step construction of the regularization with the final rigorous statement on the projection of the solutions and related proof; Section 3 is dedicated to the discussion of fast close encounters of the ER3BP.Section 4 is dedicated to numerical demonstrations: Subsection 4.1, after a short description of the considered scenarios, deals with numerical explorations in a neighborhood of P 2 and outlines quantitatively the gain as regards the computational effort; Subsection 4.2 reports examples of detection of fast close encounters with RFLIs in the ER3BP.The details about the transformations which are needed to implement numerically the KS regularization are given in the appendix.
Let us now introduce a local regularization on the secondary body2 P 2 .Following [22,23] we introduce the KS space map as a projection from a space of redundant variables u 1 , u 2 , u 3 , u 4 to a space of Cartesian variables q 1 , q 2 , q 3 : where (q 1 , q 2 , q 3 , 0) = A(u)u (10) are related to (x, y, z) by and is a matrix that plays a central role in the KS regularization.In particular, A(u) fulfills the two properties: it is a linear homogeneous function of u 1 , u 2 , u 3 , u 4 and satisfies: where I is the 4-by-4 identity matrix; hence ∥u∥ 2 = d 2 .
In this article we exploit directly such transformation in the elliptic framework by adapting the Hamiltonian derivation of the KS regularization developed in [7] for the CR3BP.We prove that a KS regularization with respect to the secondary body P 2 of the ER3BP is represented by the Hamiltonian: where u, U = (U 1 , U 2 , U 3 , U 4 ) are the KS variables and their conjugate momenta and Φ is an action conjugate to ϕ introduced to make autonomous the problem.The vector b(u) is defined by: Specifically we show that, by adopting a fictitious time s as new independent variable, the solutions (u(s), ϕ(s), U (s), Φ(s)) of Hamilton equations related to K(u, ϕ, U, Φ): that, for s = 0, satisfy: (iii) K(u(0), ϕ(0), U (0), Φ(0)) = 0, ϕ(0) = f 0 project (via π, the translation x → x + 1 − µ and df /ds = ∥u∥ 2 = d 2 ), locally to s = 0, onto solutions (x(f ), y(f ), z(f ), p 1 (f ), p 2 (f ), p 3 (f )) of Hamilton equations from (1).

Lagrangian formulation in the rotating-pulsating frame
Let L(x, y, z, x ′ , y ′ , z ′ , f ) be the Lagrangian of the spatial ER3BP in the rotating-pulsating frame with explicit dependence on the true anomaly3 f (the superscript denotes the derivative with respect to that): where, explicitly, The origin of the coordinate axes is now moved to one of the two singular positions and thus, as mentioned above in Section 2, we choose P 2 (x 2 , y 2 , z 2 ): Then ( 16) becomes: where the addenda q ′ × (0, 0, 1) have been dropped because they do not contribute to the Lagrange equations.

The space of redundant variables
By following [7] and applying their argument to the elliptic problem, we apply the projection map defined by (10) to the previous Lagrangian L, and we compute the function L (u, u ′ , f ) exploiting the relationship: is the bilinear form appearing in the usual KS regularization.We obtain: for b(u) expressed as in (15).The first task consists in proving the specific invariance of Lagrange equations under the transformation at issue.In practice, the solutions of Lagrange equations for L (u, u ′ , f ), which we write using the operator notation: have to be compared with the solutions of Lagrange equations for L(q, q ′ , f ), denoted by: With the following statement it turns out that this requirement is fulfilled as soon as the solution u(f ) ̸ = 0 for all f ∈ T.
) is a solution of Lagrange equations associated to L(q, q ′ , f ) as soon as u(f ) ̸ = 0.
Proof.For any smooth curve u(f ), reminding that: one gets from the chain rule: As a consequence we have: is for any f in the kernel of the matrix ∂π ∂u (u(f )) T .We claim that the kernel of T if and only if its components satisfy the system: , which admits the unique solution α = β = η = 0 as long as at least one of the components of u is different from zero.This implies

The modified Lagrangian
The second matter to tackle regards the Legendre transform (necessary to deduce in Subsection 2.3 the corresponding transformed Hamiltonian and then proceed with the development): where is an ad hoc permutation matrix coming from the bilinear form term (l(u, u ′ ) = Ωu • u ′ ), which is not invertible with respect to the generalized velocities, because the Hessian matrix is identically singular, indeed det H u ′ = 0. To overcome the degeneracy we proceed as in [7]: it is profitable to change the Lagrangian just by adding two times the square of the bilinear form (so that −2l 2 vanishes): such artifice precisely allows to restore the invertibility, thereby: is the modified Lagrangian and in fact, introducing the KS momenta is non-degenerate (thus invertible) in u ′ for u ̸ = 0.

Rotational invariance of the modified Lagrangian
The sum of the quadratic expression 2l 2 (u, u ′ ) of course alters L (u, u ′ , f ) and again one has to make sure that such action is legitimized under appropriate conditions (until now u(f ) ̸ = 0 always).Let then the investigation begin by realizing the well known remarkable symmetry property of the KS transformation.
Proposition 2. The modified Lagrangian L(u, u ′ , f ) is invariant under the one-parameter family of transformations involving the redundant coordinates: where S θ ∈ SO( 4) is the four-dimensional rotation matrix whose orbits define the fibers of the projection π, i.e., π(S θ u) = π(u) for all θ ∈ T.More precisely: Proof.This is a property of the KS transformation, which is the same for the CR3BP and for the ER3BP.Thus, for the proof we refer to [7], Section 2.
This fact implies that there exists, by Noether's theorem, a conserved quantity: which is an autonomous first integral for the Lagrangian L. For convenience the final constant of motion is given by: If the bilinear form is cleverly zeroed out at f = 0 (by proper initial conditions), it will keep taking zero value for further f (since we only consider time intervals such that ∥u∥ ̸ = 0), so the extra factor 2l 2 would become a vanishing contribution to the Lagrange equations.According to such idea, the bilinear form assumes the meaning of constraint to be respected along the motion and the final claim, whose proof in the proposition below resolves completely the issue, is that Lagrange equations associated to L have the same solutions of the Lagrange equations associated to L .Proposition 3. If u(f ) is a solution of the Lagrange equations of L(u, u ′ , f ) with initial data u(0), u ′ (0) satisfying u(0) ̸ = 0 and l(u(0), u ′ (0)) = 0, then it is also a solution of the Lagrange equations of L (u, u ′ , f ) as long as u(f ) ̸ = 0.

The regularized Hamiltonian
The corresponding singular Hamiltonian enters now by performing the Legendre transform: where is the inverse of U with respect to u ′ ; more explicitly we have: and the bilinear equality l(u, u ′ ) = 0 straightforwardly translates in the Hamiltonian formalism as l(u, U ) = 0 for u ̸ = 0, because With all this in hands it is useful to work with an autonomous extension of the transformed Hamiltonian K .So we append one more degree of freedom to form the extended phase space T * ((R 4 \ C ) × T), where is the collision set in KS coordinates, with the extra couple of variables (ϕ, Φ) ∈ T×R and standard symplectic form 4 i=1 du i ∧dU i +dϕ∧dΦ, in such a way to build the autonomous transformed Hamiltonian: and consider the solutions u(f ), ϕ(f ), U (f ), Φ(f ) of the Hamilton equations of ( 35) such that, for given initial value f 0 of the true anomaly, satisfy: At this point we perform a rescaling similar to the one in the Levi-Civita regularization, and define the regularized Hamiltonian: Remark 1.
-For ε = 0, the action Φ is a constant of motion and the Hamiltonian ( 36) is identical to the KS Hamiltonian of the CR3BP, as represented in [7] with Φ = −E.
-K(u, ϕ, U, Φ) is invariant under the same one-parameter family of transformations defined by ( 28) and ( 29), hence J (u, g(u, U )) = l(u, U ) is a first integral also for the Hamilton equations of K(u, ϕ, U, Φ).

Projection of the solutions of the regularized Hamiltonian
Let us prove that the solutions of the Hamilton equations of the regularized Hamiltonian (36) project on the Hamilton solutions of the Hamiltonian (1) of the ER3BP.Similarly to the classic LC and KS techniques we need an independent variable redefinition, which for the ER3BP is obtained by introducing the fictitious true anomaly s such that: whose inverse is precisely ∂K/∂Φ.Thereby we state our result.
Proof.In light of what already derived in the previous subsections, we only need to prove the equivalence between the solutions associated to the transformed K and the regularized K. Given the initial conditions u 0 , U 0 , f 0 , let us consider the solution ( u(s), ϕ(s), U (s), Φ(s)) of the Hamilton equations of K with: and s in a neighborhood of s = 0 such that ∥ u(s)∥ > 0; in particular we have: for all s.Next, consider: which is invertible (since in the neighborhood of s = 0 we have ∥ u(s)∥ > 0), and (u(f ), ϕ(f ), U (f ), Φ(f )) defined by: We claim that (u(f ), ϕ(f ), U (f ), Φ(f )) are the solutions of the Hamilton equations of K with initial conditions (u In fact, we have: as well as: and: where to obtain the second equality we used K ( u(s), ϕ(s), U (s), Φ(s)) = 0. Finally, we also have: Remark 2. The reason for the success of this Hamiltonian regularization is, as for the spatial CR3BP, the possibility to exploit the symmetry presented in Proposition 2 in the framework of symplectic reductions (see [30,31,37,47]) in the 10-dimensional phase-space (u, ϕ, U, Φ).

Fast close encounters in the ER3BP
The regularized Hamiltonian ( 14) has the representation: where has Taylor expansion with respect to the variables u starting at order 6.In fact, the Taylor expansion of 1/ ∥π + (1, 0, 0)∥ + π 1 − 1 with respect to π 1 , π 2 , π 3 starts at order 2, and the π j are quadratic functions of the u i .We provide upper bounds to the remainder R 6 and its derivatives in the ball: projecting to the ball in the Cartesian space: We will refer to both spheres ( 41) and ( 42) as the Hill's sphere (notice that the terminology is different from the conventional one by a numerical factor in the radius of the sphere).Below, we denote by µ 0 , by c 1 , c 2 , . . .> 0 and by γ 1 , γ 2 , . . .> 0 constants which are independent of µ and ε, as well as on the parameter Γ 0 which will be later introduced to characterize each close encounter.
(ii) We consider solutions q(f ) of the ER3BP transiting through the Hill's sphere of the secondary body , and denote by (u(s), ϕ(s), U (s), Φ(s)) a solution of the Hamilton equations of the regularized Hamiltonian K satisfying the hypotheses of Subsection 2.4 such that π(u(s)) = q(f (s)) for all s ∈ [s 1 , s 2 ] (s 1 , s 2 corresponding to f 1 , f 2 ).If µ ≤ µ 0 , there exists a constant c 5 > 0 such that, during the transit in the Hill's sphere, the variation of the parameter Γ s introduced in Section 1 satisfies: The proof is reported in Subsection 3.2.
(v) Finally, we prove that if µ satisfies the inequality: with suitable c > 0, the motion exits from the Hill's sphere within a proper time interval satisfying (48).
(vi) Therefore, we have proved that if µ satisfies a smallness condition, the matrix X 0 is hyperbolic during the transit in the Hill's sphere.Moreover, if the Hill's radius is suitably small (the smallness condition depending possibly on Γ 0 ), which is obtained by introducing another smallness condition on µ, then also the matrix X is hyperbolic during the transit in the Hill's sphere.Indeed, in virtue of ( 6), (7) and the above bounds, the eigenvalues of X differ from those of X 0 by terms of order µ 1/3 .43),( 44),( 45),( 46) 44) and ( 45) follow immediately from the fact that the vector b(u) is cubic in the u i .Indeed:
To prove inequalities (43) and ( 46) we first replace in R 6 the only non-polynomial term with its Taylor expansion of order 1 in the π 1 , π 2 , π 3 with remainder ∆ 4 (of order 2 in the π j , of order 4 in the u i ): and we obtain: Since for µ ≤ µ 0 < 1/10 we have 3 ) as well as: we obtain: 2 is a polynomial in π 1 , π 2 , π 3 with terms of order 2, 3, 4, there exists a constant γ 1 such that: Inequality ( 43) follows easily.We proceed by estimating the derivatives of ∆ 4 : Since the π h are quadratic function of the u j , with |π h | ≤ ∥u∥ 2 , and since 1 − d 6 1 is a polynomial containing terms of least order in the π i equal to 1, there exists a positive constant γ * such that: and a positive constant γ 2 such that: ) Inequality (46) follows straightforwardly.

The advantage of regularization: a numerical test
In order to assess the effectiveness of KS regularization of the ER3BP near the singularity at P 2 for numerical integrations we consider a fictitious simple scenario which is nevertheless representative (for the choice of the initial conditions) of realistic close encounters in the Solar System, such as the non-coplanar close encounters of comets with Jupiter, in the Sun-Jupiter ER3BP (here identified by the values µ = 9.536433730801362 • 10 −4 , ε = 0.0489).Precisely we consider the case, critical for the numerical integrations, of fast close encounters, where here we mean 'fast' both in the sense of Section 3 (see the caption of Fig. 1) and that the close encounter does not produce a temporary capture through the Lagrangian points L 1 , L 2 .We emphasize that fast close encounters are observed for celestial bodies in the Solar System (see for example [15,16] where the dynamics of comet 67P Churyumov-Gerasimenko, target of the recent Rosetta mission, is discussed in detail), and are used in Astrodynamics to accelerate spacecrafts.
In this subsection we consider a model example numerically integrated using a quadruple precision floating-point format with an explicit fixed step numerical integrator of the Runge-Kutta family (Luther's method [29], RK6 for brevity) to analyze possible gains in the use of the KS regularization.The choice of using explicitly a fixed step integrator is motivated by the need to avoid any interference of a variable step strategy with the regularization, which automatically performs the reduction of the step size by adopting a fictitiuos independent variable.We also remark that, even if the RK6 integrator is not symplectic, it does not produce a relevant energy loss since fast close encounters occurr in small time intervals.
Therefore, we choose orbits with initial conditions characterized by a high initial energy , where J is the f -dependent Jacobi "integral"4 : having, in the vicinity of the secondary body, an important deflection of the trajectory with respect to the solutions of the Kepler problem defined by the Sun.Moreover, we consider orbits which are non-planar with Jupiter's orbit.An efficient way to visualize the outcome of the regularization on fast close encounters is to consider the initial conditions already at their minimum distance from Jupiter, with E > E 4 , with inclination different from 0, then to numerically integrate the orbit by first running a backward integration up to a sufficiently large distance (at least d 2 = ∥P − P 2 ∥ > 1 = ∥P 2 − P 1 ∥), and then to switch to a forward integration lasting exactly twice the number of iterations of the previous operation, so that the upshot produces almost equal-length branches before and after the encounter (blue and red lines in top panels of Fig. 1).In such a way, we appreciate the whole dynamics with the deviation caused by the planet.More details on the choice of the initial conditions are given in the caption of Fig. 1; the physical parameters are drawn from the NASA Planetary Fact Sheet. 5The details of implementation of the numerical integration, including also the formulas for the computation of the initial condition in (u, ϕ, U, Φ) variables for any given Cartesian initial condition (x, y, z, p 1 , p 2 , p 3 ), are provided in Appendix A.2.In Figs. 1 and 2 we illustrate the particle's orbit: we set initial values so that the planetary flyby is of hyperbolic gravity-assist type and the heliocentric paths before and after the encounter are almost Keplerian ellipses.
In Table 1 we compare the numerical integration of the close encounter represented in Figs. 1 and 2 using both the Cartesian and the KS regularized equations of motion on a smaller true anomaly range, for different values of the integration steps.Precisely, using as initial condition the point of closest approach, we integrate the equations of motion backward up to f − ≈ −0.5, and then forward up to f + ≈ 0.5, since f + − f − ≈ 1 represents  approximately the interval of the true anomaly in which the close encounter takes place.We compare both the conservation of the extended Hamiltonians H = H + Φ, K and the value of ∥r∥, r = (x, y, z), at the end of the two numerical integrations.Since the numerical integration of the Cartesian equations of motion is carried out with a fixed step ∆f , while the numerical integration of the regularized equations of motion is carried out with a fixed step ∆s, in order to compare the outputs of the two numerical integrations at the same values of f we add a further iteration to the Cartesian numerical integration to reach the final value of f obtained with the KS numerical integration.We note the sharp advantage of the KS algorithm in terms of computational efficiency, providing a sharp reduction of the number of iterations needed to maintain a high level of accuracy in the representation of the trajectory around P 2 .We finally argue that, to maximize the computational efficiency in a long simulation, one should implement as usual a switching tool which works with the H-equations out of a ball centered at P 2 , say the body's Hill's sphere, and the K-equations once inside.
As conclusive remark, since K and l are smooth functions in a neighborhood of P 2 , we get that the invariances given by K(u, ϕ, U, Φ) = 0, l(u, U ) = 0 are numerically satisfied with high precision, and eventually the machine precision is quickly reached after lowering ∆s by about a factor of a hundredth (in a convergence profile, the rate of decay reflects the accuracy of the method, i.e., O(∆s  For each case the step size is eventually adapted in order to evaluate the singular and regularized solution at same corresponding times (see text).The asterisk in the last Cartesian experiment indicates failure of the adopted numerical method (explicit RK6) to compute the orbit accurately.Top panel: norm of the solutions (rounded to 16 significant digits) and total number of iterations.Bottom panel: degree of conservation of the singular extended Hamiltonians (rounded to 11 significant digits), where H(0) = K (0) = 0.

Detection of multiple close encounters with RFLIs
Regularized fast Lyapunov indicators, introduced in [8, 25,13] (see also [19]), are defined for the ER3BP by: RFLI(ξ 0 , w 0 ; f 0 , F ) = max where: • ξ 0 = (r 0 , r ′ 0 ) is the initial position and initial velocity vector given in Cartesian coordinates; f is the true anomaly of the elliptic motion and f 0 , F its initial and final values; s(f ) denotes the proper time expressed as a function of f ; • w(s) is the solution of the variational equations of Hamiltonian K obtained from the variational matrix ( 6) with initial conditions w(s(f 0 )) = w 0 , computed for an orbit with initial conditions u(f 0 ), ϕ(f 0 ), U (f 0 ), Φ(f 0 ) provided by a local inversion of the KS transformation (cf.( 77)-(79) in Appendix A.2; by changing the local inversion map, a change in the initial tangent vector should be applied accordingly, see [17]; however, in the experiments below we do not need to apply any local inversion map for the choice of the initial tangent vector and it suffices to fix the same vector in R 8 for all the integrations, see [17] for more details).
According to [14,26,16,17], we also consider the modified indicator mFLI which is suitable to detect the fast close encounters of the ER3BP: where: • q(f ) = r(f ) − (1 − µ, 0, 0) represents the solution of the equations of motion with initial conditions ξ 0 ; • f (s) is the true anomaly of the elliptic motion expressed as a function of s; • χ(q) is a function depending on a parameter λ > 0: where ∥q∥ = d 2 is the Cartesian distance ∥P − P 2 ∥ and λ is a parameter, that for close encounters is conveniently set as from 1 to 2 Hill's radii r h of P 2 .In the following experiments, we aim to detect multiple fast close encounters for small values of µ < 1/10 using the indicators RFLI and mFLI above.We consider the Sun-Earth spatial ER3BP (µ = 3.00347 • 10 −6 , ε = 0.0167).
We start from a planar reference orbit characterized by a remarkably fast close encounter with the Earth, i.e., with Γ 0 > 0 sufficiently large entering the Hill's sphere; then we suitably vary initial conditions (x, x ′ ) in a neighborhood of those of the reference orbit to explore the nearby phase space.We take initial conditions corresponding to a 2:3 mean-motion resonant orbit in the Keplerian approximation intersecting the Earth's orbit, with pericenter interior to the Earth's trajectory and aligned with the planet's one on the X-axis of the inertial frame.We synchronize times in order to have an exact collision in the future at the first intersection on the Keplerian ellipses by starting with P and P 2 located well apart, say, for convenience, at a relative distance expressed in the inertial frame.The resulting orbit is integrated in Fig. 3, complemented by Γ s during the Hill's sphere crossing.
At this stage, we compute the indicator (70) on a grid of initial conditions (x, x ′ ) equally spaced in the two coordinates, while keeping fixed y, z, y ′ , z ′ and setting Φ = −H(x, y, z, p 1 , p 2 , p 3 , f 0 ) to perform the regularization.We slightly move away from the reference orbit until the mFLI portrays fast close encounter structures in the phase space and we therein increase the resolution.Additionally, for a better visualization and to include multiple close encounter orbits, we extend the integration interval from f 0 up to F = f c + 15π + 2, with f 0 , f c given in the caption of Fig. 3.The outcome reproducing fast close encounter loci is reported in the top panel of Fig. 4.
The spatial continuation of these loci can be inspected by integrating orbits with nonzero initial z ′ (or z) starting from the planar structures.In the bottom panel of Fig. 4 we compute such continuations in the (x, z ′ ) section for positive and negative z ′ upon setting, respectively, x ′ = x ′ max and x ′ = x ′ min for all the grid points, where x ′ ∈ [x ′ min , x ′ max ] in the planar chart.The lobes depicted in Fig. 4 correspond to single or multiple close encounter orbits, as well as to deeper or less deep encounters.We explore the nature of these lobes by considering the sample orbits (indicated in the same figure) ℓ M , ℓ L , ℓ U when z ′ = 0, ℓ + M , ℓ + U when z ′ > 0 and ℓ − M , ℓ − L when z ′ < 0. We plot in Fig. 5 RFLI(s), mFLI(s), Γ s for s ∈ [0, s(F )] and d 2 (s) in a neighborhood containing the Hill's sphere for each selected orbit.
In the top left plot, we observe the strong divergence of orbits with close initial conditions for all the orbits pinpointed by a growing RFLI.In the top right plot, we can distinguish between orbits facing either 1, 2 or 3 close encounters with the Earth by looking at their jumps in the mFLI.Qualitatively, we can better determine the strength of the approaches by comparing their minimum distance d 2 in the bottom right plot: we find that the orbits on the middle lobe (label 'M ') have a single strong close encounter (small min d 2 and high mFLI jump); the orbits on the lower lobe (label 'L') have a double close encounter on the plane, the first weak (high min d 2 and low mFLI jump) and the second strong, whilst a single medium-strong close encounter going into the space; the orbits on the upper lobe (label 'U ') have a triple close encounter on the plane, the first weak, the second and third medium, whilst a double close encounter going into the space, the first medium and the second weak.According to this analysis, the development of the dynamics in the third spatial dimension appears to reduce the chance for fast close encounters, which seems intuitive given the extra degree of freedom coming into play.Nonetheless, it is worth mentioning that an intrinsic limitation is due to the finite resolution of the portraits in Fig. 4: more refined grids reveal further distinct branches inside the visible lobes with possible different dynamical features in terms of  close encounter detection.Finally, in the bottom plot of Fig. 5, we verify the hyperbolicity of each encounter by looking at its variation of Γ s in the Hill's sphere (almost null, see Fig. 3), which fulfills the properties of Section 3. The occurrences, position and length of each colored line provide information on the number, moment and duration in proper time of the corresponding encounters along every orbit, where by 'duration' we mean the time spent in the Hill's sphere.
The last result that we report is Fig. 6.On the same grid (x, x ′ ) of the central zoom-in of top panel of Fig. 4, we compute the Tisserand parameter at final value f = F , which is subject to a more significant variation when close approaches occur.Comparing to Fig. 4, we can notice that only fewer and less detailed structures are replicated, like the three main lobes identified above, corresponding to the deepest encounters.This further supports the effectiveness of the indicator (70).

Figure 2 :
Figure 2: Projections on the coordinate planes of the right bottom panel of Fig. 1 on equal axis aspect ratio.

Figure 4 :
Figure 4: mFLI χ charts of the Sun-Earth ER3BP over 1000×1000 regularly spaced initial conditions close to those of Fig. 3.Here λ = r H in χ and w 0 ∈ S 7 .Top panel: planar section with three magnifications of the close encounter lobes emerging diagonally in the figure.In the central magnification we discern three main lobes (middle, lower, upper), on which sample orbits ℓ M , ℓ L , ℓ U are taken, respectively.Bottom panel: spatial sections merged to corresponding above planar magnifications of the main lobes.The effect is a continuation of the close encounter structures along the z ′ -axis.Sample orbits ℓ + M , ℓ + U are taken on the positive z ′ extensions of the main lobes (left), while sample orbits ℓ − M , ℓ − L

Table 1 :
Comparison between Cartesian and KS integration for four consecutively increased step sizes in a small neighborhood of f = s = 0 (close approach).accordingto choices of s such that the true anomaly interval is almost symmetric with respect to the origin and sample nodes are multiple integers of every ∆s considered.