Massimiliano Guzzo

Elena Lega

Theory and applications of Fast Lyapunov Indicators for the computation of transit orbits in the three-body problem

In the last decades finite time chaos indicators have been used to compute the phase-portraits of complex dynamics as well as the center, stable and unstable manifolds originating at the partially hyperbolic equilibria, and the Lagrangian Coherent Structures of aperiodic flows. While the definition of most chaos indicators is clearly inspired by the Characteristic Lyapunov Exponent theory, their use is oriented to extract all the information which is contained in the solutions of the variational equations in short time intervals. We here review through examples why the computation of short time chaos indicators is particularly powerful for those systems whose solutions may have an asymptotic behaviour very different from the short-term one, as it can be the case of sequences of close encounters in gravitational systems and the advection of particles in aperiodic flows. The main case study here considered is the computation of transit orbits in the restricted three-body problem.

Introduction

In the last decades finite time chaos indicators have been used to compute the phase-portraits of complex dynamics as well as the center, stable and unstable manifolds originating at the partially hyperbolic equilibria, and the Lagrangian Coherent Structures of aperiodic flows. In the literature we find the definition of several finite time chaos indicators: here we mention the Fast Lyapunov Indicator (FLI) [START_REF] Froeschlé | Fast Lyapunov indicators. Application to asteroidal motion[END_REF], the Finite Time Lyapunov Exponent (FTLE) [START_REF] Tang | Finite time Lyapunov exponent and advectiondiffusion equation[END_REF], the Mean Exponential Growth of Nearby Orbits (MEGNO) [START_REF] Cincotta | Simple tools to study global dynamics in nonaxisymmetric galactic potentials-I[END_REF]. Despite the differences in the definition of these chaos indicators, all of them are obtained from the characteristic Lyapunov exponents theory. For a system of first-order ordinary differential equations:

1 ẋ = F (x), (1)
where x ∈ D ⊆ R n with D open set, and F : D → R n is a smooth vector field on D, the Lyapunov Characteristic Exponent of an initial condition x 0 ∈ D and of an initial tangent vector v 0 ∈ R n \0, is defined by

LCE(x 0 , v 0 ) = lim t→+∞ 1 t ln v(t) v 0 (2)
where v(t) is the solution of the variational equation:

v = ∂F ∂x (x(t)) v (3)
with initial conditions v(0) = v 0 , and x(t) is the solution of (1) with initial conditions x(0) = x 0 . Since the ratio v(t) / v(0) provides an estimate of the amplification of the separation at time t of the solutions with initial conditions x 0 and x 0 + v(0) for any suitably small | |, the Lyapunov exponents (2) characterize the exponential separation of the solutions with initial conditions close to x 0 . With mild hypotheses on the differential equation [START_REF] Anderson | Isolating blocks as computational tools in the circular restricted three-body problem[END_REF] the limit in eq. ( 2) exists for almost all initial conditions x 0 and provides at most n possible different values. Moreover, the largest LCE for the initial condition x 0 is expected for all the vectors v 0 ∈ R n \0, except possibly for those belonging to an hyperplane of R n . As a consequence, a random choice of v 0 provides the largest LCE, and is suitable for the computation of the strongest asymptotic law of separation of solutions with close initial conditions (we refer to [START_REF] Benettin | Tous les Nombres Caracteristiques de Lyapunov sont Effectivement Calculables[END_REF] for the numerical computation of all the LCE). While the LCE provide the asymptotic law of separation, the finite time chaos indicators exploit as much as possible the information which can be extracted from the growth of v(t) on "short" time intervals (see Figure 1). In particular, the computation of finite time chaos indicators provides information which is lost in the limit of t tending to infinity for those systems whose solutions may have an asymptotic behaviour very different from the short-term one, as it can be the case of sequences of close encounters in gravitational systems and the advection of particles in aperiodic flows. A remarkable example is provided by the dynamics of the comets of the Jupiter family, having sequences of close encounters with Jupiter which possibly change the orbital parameters during the short time of the close encounter, and also can expel the comet from the Solar System. In this context the LCE are not significant to analyze the dynamics of the comet, since the most interesting dynamics is occurring in the transient time interval in which the comet is visiting the inner Solar System. A case study has been the comet 67P/Churyumov-Gerasimenko, the target of the recent mission Rosetta, whose past orbit has been analyzed with Fast Lyapunov Indicators [START_REF] Guzzo | A study of the past dynamics of comet 67P/Churyumov-Gerasimenko with Fast Lyapunov Indicators[END_REF][START_REF] Guzzo | Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr[END_REF]. One of the simplest chaos indicators is indeed the so called Fast Lyapunov Indicator (from [START_REF] Froeschlé | Fast Lyapunov indicators. Application to asteroidal motion[END_REF] and subsequent modifications, see [START_REF] Froeschlé | Graphical Evolution of the Arnold Web: From Order to Chaos[END_REF][START_REF] Lega | Theory and Applications of the Fast Lyapunov Indicator (FLI) Method[END_REF]): the FLI of an initial condition x 0 , of an initial tangent vector v 0 , at time T is defined by F LI(x 0 , v 0 ; T ) = max converging to a positive value of the LCE (see eq. ( 2)). While the convergence of the function to the LCE is appreciated after quite a long time t ∼ 10 6 , the finite time chaos indicators are computed using much smaller time intervals.

where v(t) is the solution of the variational equation [START_REF] Cardin | Integrability of the spatial three-body problem near collisions (an announcement)[END_REF] with initial conditions v(0) = v 0 , defined for the solution x(t) of (1) with initial conditions x(0) = x 0 . The closely related quantity:

λ T (x 0 ) = max 0≤t≤T sup v(0) =1 log 10 v(t) v(0)
characterizes the maximum separation of orbits with initial conditions close to x 0 in the time interval [0, T ], and can be used to estimate the reliability of a numerical computation of the solution x(t) in the time interval [0, T ]. In fact, the possible exponential divergence of the solutions with close initial conditions is responsible also of the divergence of the numerically computed solutions of a differential equation from the target solution, determining an amplification of the errors that we introduce with the numerical scheme at each integration step, or simply with the computer floating point representation. In addition, when the initial conditions are known only within an experimental error, we have an error already at time t = 0 which is amplified by the dynamics. As already remarked, the ratio v(t) / v(0) provides an estimate of the amplification of the separation at time t of the solutions with initial conditions x 0 and x 0 + v(0) for any suitably small . Therefore, the quantity λ T (x 0 ), which is obtained by taking the maximum over all possible directions v(0) and over all the times t ∈ [0, T ], provides an estimate of the amplification in the interval [0, T ] of the errors discussed above by a factor 10 λ T (x 0 ) . Therefore, λ T (x 0 ) provides the loss of precision digits due to the numerical computation of the solution in the time interval [0, T ].

For the differential equations whose transient dynamics can be very different with respect to the asymptotic ones, λ T (x 0 ) may be very different from the asymptotic value estimated from the LCE. In figure 2 we report one of these cases, with the transient time evolution of the finite time chaos indicator Λ T (x 0 ) (whose definition is very similar to λ T (x 0 ), see the caption of figure 2 for details) representing the loss of precision digits in the backward numerical computation of the orbit of comet 67P/Churyumov-Gerasimenko. The comet could be observed from Earth for the first time in 1959, after a close encounter with Jupiter which changed the orbit by reducing the perihelion distance to about 1.2 au (see Figure 2, right panel). At this perihelion distance the comet has an important emission of dust and gas, which reflect the light of the Sun and it becomes so bright to be observed from the Earth. At previous epochs, due to larger perihelion distances, the comet was darker and has not been observed from Earth, and therefore the only possibility of reconstructing its past orbit relies on numerical computations. The computation of chaos indicators constrains the time interval of validity of numerical integrations to few centuries: numerical integrations on longer time spans have a statistical interpretation, that can be analyzed with an additional FLI analysis [START_REF] Guzzo | Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr[END_REF]. In fact, the loss of precision digits that we have during a numerical integration is localized in the phase-space at certain hyperbolic dynamical structures, such as, for example, the stable and unstable manifolds of hyperbolic periodic orbits, or the collision manifolds with a massive body.

In the different context of fluid dynamics, finite time chaos indicators have provided original insights to the study of aperiodic flows, defined by non autonomous ODE: ẋ = F (x; t),

with vector field F (x; t) depending explicitly on the time [START_REF] Tang | Finite time Lyapunov exponent and advectiondiffusion equation[END_REF][START_REF] Mezić | A method for visualization of invariant sets of dynamical systems based on the ergodic partition[END_REF][START_REF] Haller | Distinguished material surfaces and coherent structures in threedimensional fluid flows[END_REF][START_REF] Haller | Lagrangian coherent structures from approximate velocity data[END_REF][START_REF] Shadden | Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[END_REF][START_REF] Wiggins | The dynamical systems approach to Lagrangian transport in oceanic flows[END_REF][START_REF] Lekien | Lagrangian coherent structures in n-dimensional systems[END_REF]. For those models, equation (4) typically represents the motion of passive particles advected by a fluid, with F (x; t) representing the velocity field of the fluid (at the point x and time t). For typical problems, such as the transport of particles in the atmosphere or in the ocean, the velocity field F (x; t) is nonperiodic, and therefore it is hard to recover asymptotic results based on the existence of invariant periodic orbits of hyperbolic tori. Nevertheless, rather than using the chaos indicators as a tool to compute the stable and unstable manifolds of invariant periodic orbits, they have been used to define structures, called Lagrangian Coherent Structures ( [START_REF] Mezić | A method for visualization of invariant sets of dynamical systems based on the ergodic partition[END_REF][START_REF] Haller | Distinguished material surfaces and coherent structures in threedimensional fluid flows[END_REF][START_REF] Haller | Lagrangian coherent structures from approximate velocity data[END_REF][START_REF] Shadden | Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[END_REF][START_REF] Wiggins | The dynamical systems approach to Lagrangian transport in oceanic flows[END_REF][START_REF] Lekien | Lagrangian coherent structures in n-dimensional systems[END_REF], see also [START_REF] Haller | Lagrangian Coherent Structures[END_REF] and references therein) where the separation of the solutions with nearby initial conditions is the strongest one during finite time intervals. This Chapter is organized as follows: in Section 2 we describe through examples some basic applications of the fast Lyapunov indicator; in Section 3 we describe a major modification to the FLI that we specifically introduced to improve the detection of stable and unstable manifolds; Section 4 is dedicated to the description of applications of the FLI and modified FLI to the restricted three-body problem; Section 5 specifically focuses on the spatial three-body problem; Section 6 is dedicated to transit orbits. Left panel:

The quantity Λ := Λ T (x 0 ) = max 0≤t≤T max j=1,...,6 w j (t) , where w j (t) are solutions of the variational equations with initial conditions w 1 (0), . . . , w 6 (0) forming a basis of R 6 (see [START_REF] Guzzo | Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr[END_REF] for more details), represents the loss of precision digits in the backward numerical computation of the orbit of comet 67P (general relativity and nongravitational forces included). Please note that in this particularly complex case, the indicator Λ T , slightly different from λ T , has been the preferred choice. Right panel: change of the orbit of comet 67P due to the 1959 close encounter with Jupiter (only one revolution of the comet before and after the close encounter is represented in the picture; Jupiter is symbolically represented in its location at the epoch of the close encounter). The left panel is reprinted from Guzzo M., Lega E., Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr, MNRAS 469, S321-S328, 2017, Figure 1.

The fast Lyapunov indicator in simple model examples 2.1 The standard map

The simplest examples to introduce applications of the FLI are provided by smooth symplectic maps. So, for convenience, we first refer to the dynamics defined by the iteration of a smooth map:

Φ : D -→ D x -→ Φ(x)
with D ⊂ R n . For any initial condition x 0 ∈ D the sequence x k+1 = Φ(x k ) defines the orbit of initial condition x 0 , and for any x 0 ∈ D and any initial vector v 0 ∈ T x 0 D the iterations of the tangent map:

T Φ : T D -→ T D (x, v) -→ Φ(x), ∂Φ ∂x (x)v
defined by the sequences:

x k = Φ(x k-1 ) v k = ∂Φ ∂x (x k-1 )v k-1
provide the time evolution of the tangent vectors, which are used to define the chaos indicators. For example we have:

LCE(x 0 , v 0 ) = lim t→+∞ 1 t ln v t v 0
as well as the finite time indicators:

F LI(x 0 , v 0 ; T ) = max 0≤t≤T log 10 v t v 0 , λ T (x 0 ) = max 0≤t≤T sup v 0 =1 log 10 v t v 0 .
Since the variational equations are linear, for any initial condition x 0 we have a complete knowledge of the tangent dynamics by computing the iterations of the tangent map with initial tangent vectors forming a basis e 1 , . . . , e n ∈ T x 0 D of the tangent space.

Let us consider as model example the standard map:

R × T 1 -→ R × T 1 (I, ϕ) -→ (I , ϕ ) = Φ(I, ϕ)
with T 1 = R/(2πZ), defined by:

ϕ = (ϕ + I)mod(2π) I = I + sin(ϕ + I) (5)
In the top-left panel of Figure 3 we represent the phase-portrait of this map for = 0.6, obtained by computing numerically the orbits of a sample of initial conditions: we appreciate the well known typical phase-portrait characterized by invariant curves and chaotic motions. In the top-right panel we represent the values of the indicator λ T computed on a very refined grid of initial conditions regularly spaced on the phase-space; the values of the chaos indicator are represented using a color scale: on each point (ϕ, I) of the grid, we represent a pixel with a color corresponding to the value of λ T (ϕ, I) (the color scale is reported on the bottom). From this panel we appreciate much more details of the distributions of regular and chaotic motions. On the bottom panel we plot the values of λ T (ϕ, I) versus ϕ, for all the points of the grid used to compute the top-right panel. The distribution of most of the points in the lower part of the panel allows to recognize that most of the orbits in the previous panel are numerically computed with high accuracy: in fact, for almost all the initial conditions we have λ T < 14. Since λ T provides an estimate of the loss of precision digits due to the numerical integration in the time interval [0, T ], the corresponding orbits are computed with high accuracy using a quadruple floating point precision, despite the strong exponential separation of nearby orbits detected for Figure 3: Numerical experiments on the standard map defined in eq. ( 5) for = 0.6. Top-left panel: phase-portrait of the map [START_REF] Ceccaroni | Halo orbits around the collinear points of the restricted three-body problem[END_REF], where ϕ, I are the horizontal and vertical coordinates on the panel respectively. Top-right panel: representation of the indicator λ T computed for T = 100 iterations on a very refined grid of initial conditions regularly spaced on the phase-space; the value of λ T is represented using a color scale: on each point (ϕ, I) of the grid, we represent a pixel with a color corresponding to the value of λ T : λ T = 0 is reported in black; λ T ≥ 6 is reported in yellow; the values intermediate between 0 and 6 are reported with the color scale represented below the panel. Bottom panel: for all the points of the two-dimensional grid used for the top-right panel, we plot the values of λ T versus ϕ. some of the initial conditions. It is also evident that the largest values of λ T are found in the chaotic area generated by the stable and unstable manifolds of the hyperbolic fixed point (ϕ, I) = (0, 0). In the left panels of Figure 4 we report the time evolution of λ T , for increasing values of T , computed for four selected initial conditions: a regular libration orbit (red orbit in the top-right panel), a regular circulation orbit (blue orbit in the top-right panel), a strongly chaotic orbit (violet orbit in the top-right panel) and a weekly chaotic orbit (black orbit in the top-right panel). For all the orbits the value of λ T increases step-wise, but for the two regular orbits the increments accumulate to a logarithmic growth, while for the two chaotic orbits the increments accumulate to a linear law. The asymptotic laws are clearly identified already for T > 50 (top-left panel) for the regular and the strongly chaotic orbits, while it requires T > 1000 for the weakly chaotic one (bottom-left panel). Despite the fact that we need time intervals of T > 1000 to extrapolate reliable values for the slope of the linear law for the weakly chaotic orbit, much smaller time intervals of T ∼ 60 are sufficient to discriminate it from the two regular ones. As a matter of fact, the top-right panel of Figure 3 computed for T = 100 (here reproduced in the bottom-right panel with the four orbits highlighted in black) allowed clearly to identify many weakly chaotic regions related to hyperbolic periodic orbits, including the region containing the selected weakly chaotic orbit. Therefore, the computation of chaos indicators allows one to discriminate regular orbits from chaotic ones in the shortest integration times.

The double gyre

The double gyre defined below has been used since [START_REF] Shadden | Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[END_REF] as a model problem of fluid dynamics to be analyzed with finite time chaos indicators. It was defined by the planar incompressible flow: ẋ = ∂ψ ∂y (x, y; t)

ẏ = - ∂ψ ∂x (x, y; t) (6)
associated to the stream function:

ψ(x, y, t) = -sin(πf (x, t)) sin(πy) where:

f (x, t) = x [1 + (x -2) sin(2πt)]
and is a small parameter. The phase-space variables (x, y)

∈ D = {(x, y) : x ∈ [0, 2] , y ∈ [0, 1]} represent Cartesian variables in a rectangular box D,
where the flow is defined, and the solutions of (6) represent motions in the box of particles advected by the velocity field defined by the stream function ψ. For = 0 the stream function ψ is autonomous and defines an integrable Hamiltonian system on the phase-space D, with phase-portrait represented in the top left panel of Figure 5 (all the panels represent the Poincaré map of the flow defined at the time T = 1). For = 0 the system is not integrable, and we appreciate the appearance of chaotic motions on the phase portraits of the Poincaré map (see the top-right and bottom panels of Figure 5). In particular, for = 0.25, chaotic orbits spread apparently over most of the phase-space. It is interesting to compare the numerically computed phase portraits of the Poincaré map with the computation of a finite time chaos indicator at different integration times T . In [START_REF] Shadden | Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[END_REF] it was shown that short integration times reveal the appearance of a peculiar structure in the phase-space D, identified as a Lagrangian coherent structure, associated to the initial conditions providing the orbits with the strongest separation of the nearby solutions in the time interval [0, T ]. In Figure 6 we report the values of the chaos indicator λ T computed on the phase-space D for the short integration times T = 1, 2, 3, 4 as well as for the longer times T = 10, 100. As the integration time T increases, the structure which is already clearly visible for T = 1, rapidly fills the phasespace. For T = 10 (bottom-left panel) the structure almost fills the phase-space and for T = 100 (bottom-right panel) we do not have enough resolution to appreciate its unfolding in the box D, so that the initial conditions with the largest FLI values appear as a uniform area. While the FLI computation at T = 100 reaches the goal of computing the phase-portrait of the Poincaré map, the intermediate computations reveal the appearance of a structure which could not be appreciated from the computation of the phase-portraits. The structure, which is closely related to the stable manifold of a periodic orbit constrained on the bottom side of the phase-space, appears as the ridge of the finite time chaos indicator, and it has been identified in [START_REF] Shadden | Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows[END_REF] as an example of Lagrangian coherent structure.

Let us recall the mathematical meaning of ridge which is used in the analysis of Lagrangian coherent structures with finite time chaos indicators. There are in the literature several definitions, adapting to different dimensions of the phasespace and solving specific issues related to the problem of transport in fluid dynamics. We here provide the simplest definition for a ridge of a function of two variables F : R 2 → R, as a curve Γ having a family of curves γ η (s) transverse to Γ at γ η (0) and such that F (γ η (s)) has a strict maximum at s = 0. Therefore, let us consider the ridge of the chaos indicator λ T computed for the double gyre with T = 10, and represented in the panels of Figure 7 with different zoom levels. To appreciate the structure of the ridge of λ T we need to zoom into the phase-space, so the small black rectangular box highlighted in the top panel is first enlarged in the center-left panel; we still need to zoom in the Representation of λ T computed on a one-dimensional grid of initial conditions having the same value of y; the grid is defined close to the point marked with a cross in the bottom left panel. Bottomright panel: Representation of the two orbits with initial conditions in the onedimensional grid defined for the bottom center panel, and there highlighted by bold points. The two orbits are represented with green color as long as they are so close that cannot be distinguished on the picture (this happens in the time interval [0, 9]); when their separation is sufficiently large one orbit is represented in blue and the other one in red (this happens in the time interval [START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF][START_REF] Guzzo | On the identification of multiple close-encounters in the planar circular restricted three body problem[END_REF]). The separation of the orbits occur within the time interval [START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF][START_REF] Froeschlé | Graphical Evolution of the Arnold Web: From Order to Chaos[END_REF], consistently with the computation of the chaos indicator for T = 10.

panel and therefore a new black rectangular box highlighted in the center-left panel is enlarged in the center-right panel. In all these zoomed panels the set of initial conditions with λ T ≥ 6 appears resolved into many curves, which we are able to associate to a ridge of the indicator λ T with high precision. In fact, by selecting a point close to one of these curves, for example the point highlighted with a cross in the center-right panel, and by computing the indicator on a very refined one-dimensional grid of initial conditions having the same value of y, we clearly find the point on the top of the ridge with an high precision. In the bottom-left panel of the Figure we represent the result of such a computation, and we identify the point on the ridge with a precision of at least 9 digits (the precision can be sharply improved by repeating the computation on smaller refined intervals close to the maximum). For the use of the computed ridge in the analysis of the dynamics, it is important to identify the points which are on opposite sides of the ridge. In fact, let us consider the solutions of the equations of motion of the double gyre with initial conditions identified by the two points highlighted in the bottom-left panel of Figure 7, which have the same value of y and a difference in the x coordinate of about 10 -7 . In the bottom-right panel of the same Figure we report both orbits: the initial conditions are in the bottom-left part of the box, and cannot be distinguished on the picture. For a certain time interval both orbits (represented in green) remain very close, and cannot be distinguished on the picture. The evident separation of the two orbits occurs as they approach an hyperbolic periodic orbit constrained in the bottom side of the box: from this moment we represent one orbit in red and the other in blue. Therefore the ridge of the chaos indicator sharply separates the initial conditions whose orbits will be deviated in opposite directions as a consequence of a close approach with an hyperbolic periodic orbit.

These numerical experiments suggested that the intersection of the stable manifolds of hyperbolic periodic orbits (or hyperbolic tori) with two dimensional surfaces of the phase-space can be computed as the ridges of finite time chaos indicators. In [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF], Guzzo and Lega have shown that this heuristic statement can be transformed into a rigorous one only if strong hypotheses are satisfied by the dynamical system, unless a major modification in the definition of chaos indicators is considered. In the next section we first analyze with more detail additional examples and counterexamples of this heuristic statement, and then we introduce a major modification in the definition of chaos indicators which enlarges the possibility of their use to detect the stable/unstable manifolds.

Modified Chaos Indicators

In [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF][START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] we introduced a major modification to the FLI to specifically improve the detection of stable and unstable manifolds in model systems and the three-body problem; in [START_REF] Guzzo | A study of the past dynamics of comet 67P/Churyumov-Gerasimenko with Fast Lyapunov Indicators[END_REF][START_REF] Guzzo | Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr[END_REF] and [START_REF] Scantamburlo | Short-period effects of the planetary perturbations on the Sun-Earth Lagrangian point L3[END_REF] the method has been adapted to the study of the collisional manifolds of comet 67P and to the detection of libration orbits around the Lagrangian point L3 of the Sun-Earth system respectively, in a model of the solar system which is compatible with the JPL digital ephemerides. In this Section we review through some basic examples the definition and use of these modified FLI.

Stable manifolds and the ridges of chaos indicators: a very basic example

Let us first consider the mechanical pendulum of Hamiltonian: The light blue color identifies the orbits and the time evolution of λ T corresponding to the initial conditions of the grids with the largest values of λ 100 . With evidence, for both grids of initial conditions the maximum value of λ T is obtained for the points which are closer to the separatrix value I(0) = 2. From the time evolution of λ T of all the points in the grids (middle-left and bottom-left panels), we appreciate that the value of the chaos indicator increases stepwise: precisely, the value increases when the orbit approaches the hyperbolic equilibrium; then, the orbit is scattered by the equilibrium and the value of λ T stabilizes; when the orbits returns close to the equilibrium λ T increases again. The values of the increments depend on the distance of the orbit from the separatrix: orbits closer to the separatrix spend longer time intervals close to the hyperbolic equilibrium; in such a time interval the chaos indicator increases with the largest slope. As a consequence, the initial conditions with the largest value of the λ T are closer to the separatrix. The condition max 0≤t≤T in the definition of

h(I, ϕ) = I 2 2 -cos ϕ (7
λ T (I 0 , ϕ 0 ) = max 0≤t≤T log 10 sup (v I 0 ,v ϕ 0 ) =1 (v ϕ t , v I t ) (8)
is important for this discussion. Let us in fact consider the following chaos indicator defined without the condition max 0≤t≤T :

λT (I 0 , ϕ 0 ) = log 10 sup

(v I 0 ,v ϕ 0 ) =1 (v ϕ T , v I T ) . (9)
The values of this chaos indicator are reported in Figure 9: on libration and circulation orbits, the value of the chaos indicator (9) increases when the orbit approaches the hyperbolic equilibrium and decreases when the orbit goes far away from it. While the neat growth is well captured by the indicator defined in [START_REF] Cincotta | Simple tools to study global dynamics in nonaxisymmetric galactic potentials-I[END_REF], the large oscillations of the indicator defined in (9) make the final result strongly dependent on the choice of the integration time, and as we appreciate in the left panel of Figure 9 the correspondence between the values of the chaos indicator at some time T and the distance from the separatrix is partially lost for λT .

Stable manifolds and the ridges of chaos indicators: counter-examples

Despite the many examples given in the literature where the localization of stable/unstable manifolds with chaos indicators is successful, there are cases where it does not work at all. We here present a trivial example and a more subtle case, introduced in [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF], which motivate the need of a major modification Trivial examples are obtained from the consideration that for all the linear equations, for example the hyperbolic linear saddle:

ẋ1 = x 2 ẋ2 = x 1
the solution to the variational equations (and so the values of chaos indicators) are the same for all the initial conditions x(0). Therefore, the possibility to detect the stable/unstable manifolds with chaos indicators must be due to the non-linear properties of the systems. This is well understood from the following two maps:

• Map (A): ϕ = ϕ + I I = I + 3 2 sin(ϕ ) (-cos ϕ +2) 2 (10) • Map (B): ϕ = ϕ + I I = I + 3 2 sin(ϕ ) (cos ϕ +2) 2 (11)
having both an hyperbolic saddle fixed point at (ϕ, I) = (0, 0). Therefore, for both maps the origin has one-dimensional stable and unstable manifolds. Let us investigate the detection of these stable manifolds with chaos indicators. In Figure 10 we report the values of the chaos indicator λ T computed for the map (A) (left panels) and the map (B) (right panels) on a refined grid of regularlyspaced initial conditions, and different integration times (T = 10 for the top panels, T = 15 for the middle panels and T = 25 for the bottom panels). For the map (A), for increasing integration times T we appreciate the appearance of a ridge (the yellow curve in the left panels) which is easily identified with the stable manifold of the origin. For the map (B), the computation of the chaos indicator (right panels) does not highlight any curve which can be identified with the stable manifold of the origin. In [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF] the difference in the behavior of chaos indicators for the two maps has been identified in the different variation of the norm of the tangent map along the unstable manifold of the origin: for the map (A) this norm has a local maximum at the origin, for the map (B) this norm has a local minimum. The first case has been associated to a success in the determination of the stable manifold with chaos indicators, the second case has been associated to a failure; we refer to [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF] for a more definite and quantitative formulation.

Modified FLI

Since it is not practical to verify if the condition on the tangent map is satisfied or not for a generic system, the problem has been overcome by a major modification in the definition of the chaos indicator obtained by including in the computation of the chaos indicators of an initial condition x(0) only the increments occurring when the orbit x(t) is in a suitable neighbourhood B of a target hyperbolic set (see [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF]).

For example, in [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF] we have defined the modified FLI:

mFLI u (x(0), v(0), T ) = max t∈[0,T ] t 0 u(x(s)) v(s) • v(s) v(s) 2 ds (12)
where u(x) is a test function for the set B, called window function hereafter.

Remark. The integral definition of the mFLI is required to weight the growth of the chaos indicator accordingly to the distance of x(t) from the set B. In fact, since from standard calculus we have:

d ds ln v(s) = v(s) • v(s) v(s) 2 ,
we first notice that we can rewrite the traditional FLI in integral form:

max t∈[0,T ] ln v(t) v(0) = max t∈[0,T ] (ln v(t) -ln v(0) ) = max t∈[0,T ] t 0 v(s) • v(s) v(s) 2 ds,
so that for small ∆t we identify the contribution ∆F LI = v(t)• v(t) v(t) 2 ∆t as the growth of the chaos indicator in the time interval [t, t+∆t]. Let us now consider a test function u(x) for the set B: this means that u(x) = 1 for all x ∈ B, while outside B the function decreases rapidly to u(x) = 0. In the modified FLI defined in Eq. ( 12) the growth ∆F LI is weighted by the function u(x(t)), so that if x(t) ∈ B then u(x(t)) = 1, and the growth is fully accounted; if x(t) is suitably far from B then u(x(t)) = 0, and in the interval [t, t + ∆t] there is no contribution to the integral.

When we use functions u(x) which are test functions of a subset B of the phase-space, the mFLI is higher when computed on initial conditions of solutions which spend most of the time interval [0, T ] inside the set B and close to an hyperbolic set of B. When the set B is a suitable neighbourhood of an hyperbolic invariant set γ, for example an hyperbolic or partially hyperbolic periodic orbit or invariant torus, the mFLI is higher when computed on points which are on the stable manifold of γ. For example, when γ is an equilibrium or a periodic orbit we compute its stable/unstable manifolds using the test function:

u(x) =    1 if |x -γ| ≤ r 2 1 2 [cos(( |x-γ| r -1 2 )π) + 1] if r 2 < |x -γ| ≤ 3r 2 0 if |x -γ| > 3r 2 (13)
with a suitable parameter r > 0. Instead, to compute the stable/unstable manifolds of center manifolds (as for the spatial three-body problem [START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF]) the test function is constructed using a preliminary analytic computation of the Hamiltonian normal form obtained following the methods of reduction to the center manifold [START_REF] Gómez | Dynamics and Mission Design Near Libration Point Orbits[END_REF][START_REF] Jorba | Dynamics in the center manifold of the restricted three-body problem[END_REF][START_REF] Celletti | Lissajous and Halo Orbits in the Restricted Three-Body Problem[END_REF][START_REF] Ceccaroni | Halo orbits around the collinear points of the restricted three-body problem[END_REF][START_REF] Pucacco | Structure of the centre manifold of the L 1 , L 2 collinear libration points in the restricted three-body problem[END_REF]. In addition, one also may consider mFLIs by including in the definition a maximum computed on all the vectors v(0) satisfying v(0) = 1, as in the definition of λ T (x 0 ) (see Section 1). We remark that for an effective definition of the mFLI we do not need to localize in advance the target invariant hyperbolic set with high precision: low order perturbation theories provide good definitions of the set B and of the function u(x).

In Figure 11 we compare the values of λ T with the values of a modified FLI for the map (B): while the largest values of λ T identify only a large neighbourhood of the stable manifold (the yellow area in the left panel), the ridge of the modified FLI sharply identifies the stable manifold of the origin in the right panel.

The efficiency of the modified FLI to detect the stable manifolds has been proved in a quite general framework, we refer to [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF] for all the technical details. We remark that, if mild hypotheses are satisfied by the dynamical system, the mFLI has a small quasi-constant plateau at the top of the ridge (whose amplitude shrinks by extending T ), and its value decreases by increasing the distance. In a suitable range for the distance, the decrement that we have between two points at distances 10 -N1 , 10 -N2 is proportional to |N 1 -N 2 |.

FLI and the restricted three-body problem

The circular restricted three-body problem (CR3BP hereafter) is defined by the motion of a massless body P in the gravity fields of two massive bodies P 1 and P 2 , the primary and secondary body respectively, which rotate uniformly around their common center of mass. In a rotating frame the problem has five equilibria, the so called Lagrangian points L 1 , . . . , L 5 , which are the only known simple solutions of the equations of motion of P :

     ẍ = 2 ẏ + x -(1 -µ) x+µ r 3 1 -µ x-1+µ r 3 2 ÿ = -2 ẋ + y -(1 -µ) y r 3 1 -µ y r 3 2 z = -(1 -µ) z r 3 1 -µ z r 3 2 , ( 14
)
where the units of masses, lengths and time have been chosen so that the masses of P 1 and P 2 are 1 -µ and µ, their coordinates are (-µ, 0, 0) and (1 -µ, 0, 0) and their revolution period is 2π; r 1 = (x + µ) 2 + y 2 + z 2 and r 2 = (x -1 + µ) 2 + y 2 + z 2 denote the distances of P from P 1 , P 2 . In a series of papers [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF][START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] we developed modified chaos indicators which are suitable to detect the stable and unstable manifolds of the center manifolds originating at the Lagrange equilibria L 1 and L 2 , which are known to play a relevant role in the definition of transit motions for the CR3BP.

Let us first review some basic properties of the dynamics of the CR3BP. Since equations ( 14) have an integral of motion, the Jacobi constant:

C(x, y, z, ẋ, ẏ, ż) = x 2 + y 2 + 2 1 -µ r 1 + 2 µ r 2 -ẋ2 -ẏ2 -ż2 , (15)
for any given value C of C we consider its level set M C in the phase-space (x, y, z, ẋ, ẏ, ż), so that we have (x, y, z, ẋ, ẏ, ż) ∈ M C along the motions as well as (x, y, z) ∈ ΠM C where:

ΠM C = {(x, y, z) ∈ R 3 \0 : x 2 + y 2 + 2 1 -µ r 1 + 2 µ r 2 ≥ C}. (16)
The boundary B C of ΠM C separates the so called realm of possible motions ΠM C from the realm of forbidden motions. The Lagrangian equilibria are critical points for the Jacobi constant; the values C 1 , C 2 , ..., C 5 of C at the Lagrange equilibria L 1 , . . . , L 5 correspond to topological changes of the set B C . If C > C 1 , B C disconnects the space into a neighbourhood of P 1 , a neighbourhood of P 2 , and a region external to both bodies P 1 , P 2 ; therefore no transits between these regions are possible. The neighbourhoods of P 1 and P 2 are connected if C < C 1 , so transits between neighbourhoods of P 1 and neighbourhoods of P 2 are possible. These neighbourhoods are connected also to the external region if C < C 2 . Therefore, values of C slightly smaller than C 2 are the larger values of C for which the realm of possible motions connects a region around P 1 , a region around P 2 , and a region external to the binary system P 1 , P 2 .

In Figure 12 we represent an example of intersection of the set ΠM C with the planes z = 0 (left panel) and y = 0 (right panel) for a value C smaller than C 2 (the shaded areas are in the realm of forbidden motions). As it happens for the cases represented in the figure, the connection between these regions is realized through two bottlenecks of B C , where we find the Lagrange points L 1 and L 2 . The transit of motions through the bottlenecks is guided by the stable-unstable manifolds of the center manifolds W c 1 , W c 2 originating at the equilibria L 1 ,L 2 , which are partially hyperbolic. The intersections of these stable and unstable manifolds with M C have been identified as separatrices for the transit of motions through the bottlenecks of B C connecting the region of internal and the region of external motions, see [START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF][START_REF] Simó | Dynamical systems methods for space missions on a vicinity of collinear libration points[END_REF][START_REF] Koon | Dynamical Systems, the three body problem and space mission design[END_REF][START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF] (for the planar three-body problem) and [START_REF] Gomez | Connecting orbits and invariant manifolds in the spatial restricted three-body problem[END_REF][START_REF] Anderson | Isolating blocks as computational tools in the circular restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] (for the spatial three-body problem). The computation of the stableunstable manifolds is therefore essential to compute the different kind of transit orbits. There is a rich literature about the computation of these manifolds and the related dynamics (see for example [START_REF] Simó | Dynamical systems methods for space missions on a vicinity of collinear libration points[END_REF][START_REF] Gómez | Dynamics and Mission Design Near Libration Point Orbits[END_REF][START_REF] Gomez | Connecting orbits and invariant manifolds in the spatial restricted three-body problem[END_REF][START_REF] Koon | Dynamical Systems, the three body problem and space mission design[END_REF][START_REF] Masdemont | High Order Expansions of Invariant Manifolds of Libration Point Orbits with Applications to Mission Design[END_REF][START_REF] Celletti | Lissajous and Halo Orbits in the Restricted Three-Body Problem[END_REF][START_REF] Ceccaroni | Halo orbits around the collinear points of the restricted three-body problem[END_REF][START_REF] Pucacco | Structure of the centre manifold of the L 1 , L 2 collinear libration points in the restricted three-body problem[END_REF][START_REF] Zanzottera | Intersecting invariant manifolds in spatial restricted three-body problems: Design and optimization of Earth-to-halo transfers in the Sun-Earth-Moon scenario[END_REF][START_REF] Paez | A study of temporary captures and collisions in the Circular Restricted Three-Body Problem with normalizations of the Levi-Civita Hamiltonian[END_REF]). We shall focus on the method of computation formulated with chaos indicators.

The case study that we review here, which is particularly important for the dynamics of comets, is represented by the Sun-Jupiter mass ratio and by values of the Jacobi constant slightly smaller than2 C 2 .

Let us first consider the planar case, which is often studied on the section:

Σ P (C) = {(x, 0, 0, ẋ, ẏ, 0) : C(x, 0, 0, ẋ, ẏ, 0) = C} (17)
parameterized by the variables x, ẋ. A first snapshot of the dynamics is provided by the computation of the traditional FLI in a very refined grid of initial conditions in Σ P (C), such as the ones represented in Figure 13. Since the orbits with such initial conditions may have close encounters with the secondary body, it is more convenient to define the FLI using the equations of motion of the planar CR3BP regularized with respect to the body P 2 [START_REF] Celletti | Some results on the global dynamics of the regularized restricted three-body problem with dissipation[END_REF][START_REF] Guzzo | On the identification of multiple close-encounters in the planar circular restricted three body problem[END_REF][START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF][START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF]. A convenient regularization of the planar problem is the Levi-Civita regularization [START_REF] Levi-Civita | Sur la régularisation qualitative du probléme restreint des trois corps[END_REF]. In Figure 13 we report the computation of the FLI for a specific value of C, on a grid of initial conditions of Σ P (C) with final integration time T equal to half (left panel) and twice (right panel) of Jupiter's period. Even using such short integration times we appreciate the strong divergence of solutions with initial condition in this window of the phase-space. From these pictures we have only a vague idea of the structures forming the chaotic region, which we know are correlated to the stable manifolds originating at the Lyapunov orbits of L 1 , L 2 , but to compute them out of all the chaotic initial conditions we need to use modified chaos indicators.

Modified FLI for the planar CR3BP

In the planar CR3BP the center manifolds W c i originating at L 1 or L 2 are two-dimensional manifolds, and close to the Lagrangian point are foliated by a family of periodic orbits, called planar Lyapunov orbits. Therefore, for values C of the Jacobi constant slightly smaller than C i the intersection of the center manifold W c i with the level set M C contains a planar Lyapunov orbit, that we denote by L i (C). The orbits L i (C) are libration orbits around the equilibrium L i , are located in the bottleneck formed by the boundary B C at L i , and have stable and unstable manifolds, hereafter denoted W s Li , W u Li , which extend on both sides of the bottleneck. Orbits transiting from one side to the other side of the bottleneck can be constructed by shadowing the two branches of the manifolds W s Li , W u Li on the left and on the right of L i (C) (see [START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF]). Modified FLI have been specifically designed in [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF][START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] to detect the intersection of the manifolds W s Li , W u Li with two-dimensional sections of the three-dimensional level set M C .

Definition of the mFLI. Following [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF], we define the mFLI using the equations of motions of the planar CR3BP formulated using the Levi-Civita regularization with respect to the secondary body P 2 . The Levi-Civita regularization is defined by the space transformation

x -(1 -µ) = u 2 1 -u 2 2 y = 2u 1 u 2 ( 18
)
complemented by the time transformation:

dt = r 2 ds, ( 19
)
where t is the physical time and s is called the proper time. The equations of motion in the variables u 1 , u 2 , and proper time s are (see for example [START_REF] Szebehely | Theory of orbits[END_REF]):

u 1 = 1 4 [(a + b)u 1 + cu 2 ] u 2 = 1 4 [(a -b)u 2 + cu 1 ] ( 20
)
where the primed derivatives denote derivatives with respect to the proper time s, and a, b, c are functions of u, u which are regular as u tends to 0 (this condition corresponds to a collision with P 2 ). After rewriting equations [START_REF] Jorba | Dynamics in the center manifold of the restricted three-body problem[END_REF] as first order differential equations:

       u 1 = v 1 u 2 = v 2 v 1 = 1 4 [(a + b)u 1 + cu 2 ] v 2 = 1 4 [(a -b)u 2 + cu 1 ] (21)
and denoting them by: with ζ = (u 1 , u 2 , v 1 , v 2 ), the solutions w(s) of the variational equations of ( 22):

ζ = F (ζ) (22)
w = ∂F ∂ζ (ζ)w, ( 23
)
are used to define the mFLI. Precisely, by denoting with

(x, y, ẋ, ẏ) = W(u 1 , u 2 , v 1 , v 2 )
the transformation from the Levi-Civita to the Cartesian variables:

       x -(1 -µ) = u 2 1 -u 2 2 y = 2u 1 u 2 ẋ = 2 r2 (u 1 v 1 -u 2 v 2 ) ẏ = 2 r2 (u 1 v 2 + u 2 v 1 ), (24)
for any initial condition ξ 0 = (x 0 , y 0 , ẋ0 , ẏ0 ) we consider a solution ζ(s) of ( 22) with initial condition ζ 0 such that ξ 0 = W(ζ 0 ), and the transformation s := σ(ξ 0 , t) between any physical time t and the proper time s for this solution.

Then, for a given physical time T , we define the modified FLI by

mFLI u (ξ 0 , w 0 , T ) = max t∈[0,T ] σ(ξ0,t) 0 u(W(ζ(s)) w(s) • ẇ(s) w(s) 2 ds ( 25
)
where the window function u(ξ) is defined as in eq. ( 13), by identifying γ with the Lyapunov orbit L i (C).

Grid computations of the mFLI: heteroclinic orbits. As a demonstration of the use of the new indicator mFLI in the computation of the stable and unstable manifolds of the Lyapunov orbits, in [START_REF] Guzzo | Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Laypunov Indicators[END_REF] we computed heteroclinic orbits in W u L1(C) ∩W s L2(C) , i.e. orbits which converge in the past to the Lyapunov orbit L 1 (C) and in the future to the Lyapunov orbit L 2 (C). First, we obtained representations of the intersections of both manifolds W u L1(C) , W s L2(C) with the two-dimensional section Σ P (C) by computing two mFLI, which we denote by mFLI 1 , mFLI 2 , on the same refined grid of initial conditions of Σ P (C), using the same initial tangent vector w 0 , the same parameter r = 10 -3 in the definition of the window function u, and:

• a negative time -T 1 < 0, and γ = L 1 (C), in the definition of the window function u for the indicator mFLI 1 ;

• a positive time T 2 > 0, and γ = L 2 (C), in the definition of the window function u for the indicator mFLI 2 .

Therefore, the ridges of mFLI 1 identify the manifold W u L1(C) while the ridges of mFLI 2 identify the manifold W s L2(C) . Finally, to appreciate both manifolds W u L1(C) , W s L2(C) on the same picture (and therefore to identify their intersections) we represent a weighted average of the two indicators:

p mFLI 1 + mFLI 2 (p + 1) , (26)
for a convenient choice of the weight parameter p > 0. The results are represented in Figure 14: the stable and unstable manifolds appear on the two panels, obtained for different times, as the ridges of the indicator [START_REF] Haller | Lagrangian Coherent Structures[END_REF]. We clearly appreciate different lobes of both manifolds, as well as many intersection points, providing initial conditions of heteroclinic orbits. The number of lobes and heteroclinic intersections increases very rapidly for increasing times T 1 . The computation of the indicator [START_REF] Haller | Lagrangian Coherent Structures[END_REF] has been repeated by zooming close to a selected intersection point of the manifolds W u L1(C) , W s L2(C) , in order to detect the initial condition providing the heteroclinic orbit with all the desired precision which is necessary for a numerical computation of the orbit (in Figure 15, the heteroclinic point is obtained with more than 15 digits of precision). In Figure 15, bottom-right panel, we report the orbit corresponding to the initial conditions in the intersection of the two manifolds.

Extraction of the ridges of the mFLI

The precision of the computation of the ridges of the mFLI presented in the previous subsections is limited by the stepsize α of the two-dimensional grids of initial conditions defined on the section Σ C (P ). If one aims to compute long parts of the stable and unstable manifolds with an high precision (for example the precision α < 10 -15 obtained in the computation of Figure 15), it is convenient to follow a specific strategy indicated, for example, in [START_REF] Garth | Efficient Computation and Visualization of Coherent Structures in Fluid Flow Applications[END_REF][START_REF] Lipinski | A 3D fast algorithm for computing Lagrangian coherent structures via ridge tracking[END_REF]. In Figure 16 we sketch the main ideas of the strategies that we have implemented in the papers [START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF][START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] to compute ordered sequences of points π 0 , π 1 , π 2 , . . . on a ridge of the mFLI on the two dimensional sections of M C . Let us consider different connected sets of points with positive value of mFLI on the same two dimensional section Σ of M C . For each connected set, we select a point π0 which is a candidate to be close to a ridge of the mFLI, and we compute again the mFLI on a much more refined one-dimensional grid in Σ centered in π0 (see figure 16). The maximum value of the mFLI on this grid provides a much better estimate of the point π 0 on the ridge of the mFLI. The second point π 1 on the ridge is obtained by computing the mFLI on another one-dimensional grid of points of Σ obtained by changing only one coordinate of π 0 . The point with the maximum mFLI on this second grid provides π 1 . The third point π 2 is computed from a one-dimensional grid of N points defined on an arc whose vertex is in π 1 + (π1-π0)

|π1-π0| , with some small step-size α (see figure 16). The point with maximum mFLI on the arc is denoted by π 2 . The procedure is then iterated, providing us an ordered sample of points π 0 , π 1 , π 2 , . . . on a ridge of the mFLI on the two dimensional section Σ. The values of the parameters , N, α are adjusted at each step, as well as the total integration time.

Examples of ridges computed with this method in [START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF] are reported in Figure 17, where we represent six different ridges of the mFLI on the planar section Σ P (C): two of them, which we denote by µ 1 , µ 2 , are in the realm of motions around the primary body P 1 ; a ridge denoted λ 1 intersects L 1 (C) and extends on both sides of this Lyapunov orbit; three ridges, denoted λ 2 , λ 3 , λ 4 , are in the realm of motions around P 2 ; the red curves represent the projection on the section Σ P (C) of sample orbits connecting two ridges. Since the target orbit γ in the definition of the mFLI is the Lyapunov orbit L 1 (C), and the indicator is computed using positive times T , the computed ridges are approximation of the stable manifold W s L1 . Therefore, the points with initial conditions on the ridges λ 1 , λ 2 , λ 3 , λ 4 , µ 1 , µ 2 have been numerically integrated, to reconstruct (a part of) this stable manifold.

π 0 π 0 ∼ π 1 π 2 π 3 π 4 ε ε ε
Figure 16: Sketch of the algorithm to compute a sequence of points π 0 , π 1 , π 2 , . . . on a ridge of the mFLI. The blue circle represents the starting point, that is a point of the section Σ P (C) with a high value of the modified FLI, obtained from a computation on a 2-dimensional grid of initial conditions. Then, we refine the position of the starting point by re-computing the mFLI on a much more refined 1-dimensional grid marked by red circles centered on the blue one. The black circle corresponds to the point of this grid with the maximum value of the mFLI. Then, we consider a second 1-dimensional grid centered at a point obtained from the previous black circle, and shifting one coordinate of a small quantity . We mark again with a black circle the point of this grid with the largest value of the mFLI. The third point π 2 is computed from a 1-dimensional grid of N points chosen on an arc whose vertex is in π 1 + (π1-π0)

|π1-π0| , with some small step-size. The maximum mFLI on the arc grid provides the point π 2 . The procedure is now iterated, providing us an ordered sample of points on a ridge of the mFLI. Figure reprinted from: Physica D, vol. 325, Lega E. and Guzzo M., Three-dimensional representations of the tube manifolds of the planar restricted three-body problem, 41-52, Fig. 3, Copyright 2016, with permission from Elsevier.

Reconstruction of the stable manifolds

The stable manifold of the selected Lyapunov orbit L 1 (C) has been computed in [START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF] from the forward and backward numerical computations of the orbits with initial conditions on the curves λ j and µ j introduced in the previous Subsection.

Let us briefly analyze how these orbits connect the different curves µ 1 , µ 2 , λ 1 , λ 2 , λ 3 , λ 4 . The lobe λ 1 (represented with a bold curve in the bottom-right panel of Figure 17) is different from all the others, since the orbits with initial conditions on λ 1 converge directly to L 1 without performing circulations around the primary or secondary body. Instead λ 4 , λ 3 , λ 2 and λ 1 are connected by orbits which perform a given number of half circulations around P 2 , in the sequence indicated above (see the red curves in the panels of Figure 17); µ 2 , µ 1 , λ 1 are connected by orbits which perform half circulations around P 1 , in the sequence indicated above (see the red curves in the panels of Figure 17). In Figure 18 we represent the projection on the xy plane of orbits with initial conditions on λ 2 , numerically computed forward (red curves, we appreciate that all of these curves converge to L 1 (C)) and backward (blue curves) in time. In Figure 19 we represent the projection of sample orbits with initial conditions in λ 2 in the three-dimensional space x, y, ẋ: the pink surface is obtained from the forward integration and the violet surface is obtained from the backward integration. From the representation it appears clearly how the computed orbits connect the curve λ 2 to the curves λ 1 , λ 3 . As we extend the integration time of the orbits with initial conditions on the curves λ j we appreciate the complicate development of the stable manifold of the Lyapunov orbit; for additional pictures (also for the manifolds related to L 2 (C)) we refer to [START_REF] Lega | Three-dimensional representations of the tube manifolds of the planar restricted three-body problem[END_REF]. We here review how the heteroclinic orbits, which are in the intersections of the stable manifolds of one Lyapunov orbit, for example L 2 (C), and of the unstable manifold of the other Lyapunov orbit, for example L 1 (C), provide an idea of the complicate structure of the stable and unstable manifolds.

In Figure 20, top-left panel, we represent on the x, ẋ plane some connected components of the intersection of both stable and unstable manifolds of the Lyapunov orbits L 1 (C) and L 2 (C) with the section Σ P (C), computed as the ridges of the mFLI. The black curves in the figure highlight two pieces of the stable manifold of L 1 (C) and of the unstable manifold of L 2 (C), which intersect transversely. The intersection points belong therefore to an heteroclinic orbit connecting all the neighbourhoods of the two Lyapunov orbits. On the top-right and the bottom panels we report the three-dimensional part of the stable and unstable manifolds obtained by numerically computing the orbits with initial conditions on the black curves (only some parts of these strips are represented, in order to obtain a better visualization).

5 Modified FLI and the spatial CR3BP

In the spatial CR3BP the center manifolds W c 1 , W c 2 originating at the Lagrangian points L 1 , L 2 intersect the set M C in surfaces diffeomorphic to threedimensional spheres (provided that C < C i and |C -C i | is suitably small), which we denote by CC 1 , CC 2 ; their stable and unstable manifolds are therefore four dimensional hypertubes. In [START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] we provided examples of different computations and visualizations of these manifolds obtained from the computations of mFLI for the spatial CR3BP. Let us recall the main differences of the mFLI method introduced to study the spatial with respect to the planar CR3BP.

Figure 19: Projection on the x, y, ẋ space of a set of orbits obtained from backward numerical integrations of initial conditions of λ 2 (the xy plane is the horizontal one in the pictures, the ẋ axis is the vertical direction). The pink surface is obtained from the forward integration of orbits with initial condition on λ 2 ; the violet surface is obtained from the backward integration of the same initial conditions; the sets λ 2 , λ 3 , λ 4 as well as the Lyapunov orbit L 1 (C) (the perspective is such that the positive direction of the x axis points to the left) are highlighted as the border of blue, pink, cyan, green flat surfaces. The transparent mesh is used to represent the border of the level set M C in the x, y, ẋ space.

Definition of the mFLI

The definition of the mFLI for the spatial CR3BP uses the Kustaanheimo-Stiefel regularization [START_REF] Kustaanheimo | Spinor regularisation of the Kepler motion[END_REF][START_REF] Kustaanheimo | Perturbation theory of Kepler motion based on spinor regularization[END_REF] at P 2 , KS hereafter, which is the spatial extension of the Levi-Civita regularization (see also [START_REF] Cardin | Integrability of the spatial three-body problem near collisions (an announcement)[END_REF][START_REF] Cardin | Integrability of close encounters in the spatial restricted three-body problem[END_REF]). The KS regularization is defined by the projection map π : (u 1 , u 2 , u 3 , u 4 ) → (x, y, z):

   x = 1 -µ + u 2 1 -u 2 2 -u 2 3 + u 2 4 y = 2(u 1 u 2 -u 3 u 4 ) z = 2(u 1 u 3 + u 2 u 4 )
and by the time transformation dt = r 2 ds where t is the physical time and s is called the proper time. The KS transformation regularizes the equations of motion of the spatial CR3BP for any fixed value C of the Jacobi constant, associating to the equation of motions [START_REF] Gómez | Dynamics and Mission Design Near Libration Point Orbits[END_REF] represented in the Cartesian variables, the equation of motions which we denote by:

u = F (u, u ; C) (27)
represented in the regularized u = (u 1 , u 2 , u 3 , u 4 ) variables; the function F (u, u ; C) depends parametrically on C and is regular at u = 0 (corresponding to the singularity P 2 ). In addition, eq. ( 27) is complemented by the non-holonomic relation:

u 1 u 4 -u 2 u 3 + u 3 u 2 -u 4 u 1 = 0.
In [START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF], we defined the chaos indicators for the spatial CR3BP from the solutions w(s) ∈ R 8 of the variational equations of eq. ( 27): In fact, since in the spatial problem we deal with three-dimensional subsets of the center manifolds CC 1 , CC 2 , it is not computationally convenient to define a neighbourhood of the sets CC 1 , CC 2 from a sample of their points. Rather, we prefer to construct the window functions φ(ζ) using the Hamiltonian reductions to the center manifolds [START_REF] Gómez | Dynamics and Mission Design Near Libration Point Orbits[END_REF][START_REF] Jorba | Dynamics in the center manifold of the restricted three-body problem[END_REF][START_REF] Celletti | Lissajous and Halo Orbits in the Restricted Three-Body Problem[END_REF][START_REF] Ceccaroni | Halo orbits around the collinear points of the restricted three-body problem[END_REF][START_REF] Pucacco | Structure of the centre manifold of the L 1 , L 2 collinear libration points in the restricted three-body problem[END_REF]. The reduction to the center manifold define approximations of the center manifolds from the construction of phasespace coordinates γ, η, such that the approximated manifold is obtained as the level sets η j = 0. Therefore, the numerical computation of j |η j | defines a convenient distance from the manifold, which is used to define the function φ(ζ).

w = [F(u, u )]w , w ∈ R 8 . (28)

Sections of the stable manifolds

In [START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF] we computed the ridges of the mFLI defined in (29) on the two dimensional vertical sections: Σ V (C) = {(x, 0, z, 0, ẏ, 0) : C(x, 0, z, 0, ẏ, 0) = C}, using two window functions φ 1 (ζ), φ 2 (ζ) constructed to detect the stable/unstable manifolds of CC 1 and CC 2 respectively. As for the planar case, for a fixed value of C we obtained a family of curves characterized by a given number of revolutions around P 2 before converging to CC 1 or CC 2 ; with evidence, each ridge of the vertical section intersects the corresponding ridge of the planar section Σ P (C). Figure 21, top panel, represents a sample of ridges of the indicator mFLI defined in [START_REF] Lekien | Lagrangian coherent structures in n-dimensional systems[END_REF], computed in the vertical section Σ V (C), representing connected components of the intersection of both stable and unstable manifolds of CC 1 (black curves) and CC 2 (blue curves) with the section; a magnification is represented in the bottom-left panel, where we highlight in blue two connected components of the ridge λ 4 , characterized by initial conditions whose orbits converge to CC 1 after two revolutions around P 2 . The bottom-right panel reports in the x, ẋ, z space the ridge λ 4 computed in the vertical section Σ V (C), as well as the ridge λ 4 previously computed in the planar section Σ P (C) (see Figure 17), and we appreciate their intersections on the line z = 0, ẋ = 0. The development of the stable and unstable manifolds of CC 1 and CC 2 can be obtained from the numerical integrations of the solutions with initial conditions on the computed ridges of the mFLI. In Figure 22 we represent the projection of the orbits associated to solutions with initial conditions in the lobe λ 2 of Σ V (C) on the x, y plane (top-left panel), on the x, z plane (top-right panel) and on the x, y, z space (bottom panel). For additional computations and representations we refer to the paper [START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF].

Transit orbits

The stable and unstable manifolds of CC 1 , CC 2 (as well as of the Lyapunov orbits L 1 (C), L 2 (C) for the planar problem) are separatrices for the transit of motions through the bottlenecks in the realm of admissible motions which open at L 1 , L 2 for values of the Jacobi constant slightly smaller than C 1 or C 2 respectively.

For typical applications to Celestial Mechanics, such as the dynamics of comets and the space mission design, we are interested in the computation of orbits performing transits occurring within given time intervals. Also, it is important to quantify the numerical errors that we have in the computation of the stable/unstable manifolds. Let us analyze in this Section these two issues with some more detail.

Since the computation of the stable and unstable manifolds is always affected by ( possibly very small) numerical errors, the orbits with initial conditions ξ on the numerically computed stable manifold of CC 1 will approach the set CC i but then (due to numerical errors) may transit on the other side of the bottleneck opened at CC i , or may bounce back, depending on which side of the stable manifold the initial condition ξ actually is. As a matter of fact, the meaning of the numerical computation of each ridge of mFLI is the determination of an annulus σ in the two dimensional section Σ P (C) or Σ V (C) delimited by an inner curve and an outer curve, whose points are at distance from the ridge larger than the numerical error, and have a well defined transit property. Let us consider, for example, the ridge λ 2 in the planar section Σ P (C) and the very thin annulus σ ⊆ Σ P (C) containing λ 2 and identified as indicated above. In Figure 23 we represent the projection on the xy plane of orbits with initial conditions in the inner border (left panel) and in the outer border (right panel) of the annulus. All the orbits with initial conditions in the inner border approach the Lyapunov orbit L 1 (C) and then transit to the left of the orbit, while all the orbits with initial conditions in the outer border approach the Lyapunov orbit L 1 (C) and then bounce back. In Figure 24 we represent with colors the transit properties, within a fixed integration time T , of the orbits with initial conditions in the vertical Section Σ V (C): violet represents initial conditions whose orbit transits towards P 1 within the time T ; orange represents initial conditions whose orbit transits towards the external region within the time T (see the caption of Figure 24 for more details).

The result provides a clear correlation between the transit properties with the computation of some of the ridges λ j Σ V (C) previously provided by the computation of the mFLI defined for the spatial problem. The orbits which do not transit in the time interval [0, T ] (yellow color in the picture) may of course transit in longer time intervals, so that by increasing the integration time T additional ridges have to be considered to define the transit properties in the section Σ V (C).

Figure 24: For any initial condition on a refined grid of 10000 × 2000 initial conditions of Σ V (C) we computed their orbits until they reach a value for the x coordinate smaller than 0.933 or larger than 1.067, within the maximum integration time T = 10. In the former case, we represent the initial condition with a violet pixel, in the latter we represent the initial condition with an orange pixel. If none of the conditions is met, we represent the initial condition with a yellow pixel.

Figure 1 :

 1 Figure 1: Typical time evolution of 1 t ln v(t) v 0

  Figure 2:Left panel:The quantity Λ := Λ T (x 0 ) = max 0≤t≤T max j=1,...,6 w j (t) , where w j (t) are solutions of the variational equations with initial conditions w 1 (0), . . . , w 6 (0) forming a basis of R 6 (see[START_REF] Guzzo | Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr[END_REF] for more details), represents the loss of precision digits in the backward numerical computation of the orbit of comet 67P (general relativity and nongravitational forces included). Please note that in this particularly complex case, the indicator Λ T , slightly different from λ T , has been the preferred choice. Right panel: change of the orbit of comet 67P due to the 1959 close encounter with Jupiter (only one revolution of the comet before and after the close encounter is represented in the picture; Jupiter is symbolically represented in its location at the epoch of the close encounter). The left panel is reprinted from Guzzo M., Lega E., Scenarios for the dynamics of comet 67P/Churyumov-Gerasimenko over the past 500 kyr, MNRAS 469, S321-S328, 2017, Figure1.

Figure 4 :

 4 Figure 4: Left panels: Time evolution of λ T for four selected orbits of the standard map defined by = 0.6 (top-left panel for T ∈ [0, 100], bottom-left panel for T ∈ [0, 10000]). Top-right panel: the four orbits whose FLI evolution is represented on the left panels; the color of each orbit is associated to a curve on the left panels. Bottom-right panel: the four orbits are represented in black on the color map representation of the indicator λ T .

Figure 5 :

 5 Figure 5: Phase portraits of the Poincaré map of the double gyre defined by the flow of (6) at time T = 1 for = 0 (top-left panel), = 0.05 (top-right panel), = 0.1 (bottom-left panel), = 0.25 (bottom-right panel).

Figure 6 :

 6 Figure 6: Color representation of the chaos indicator λ T computed for the double gyre defined by = 0.25 on a refined grid of regularly-spaced initial conditions on the xy plane. The integration times are T = 1 (top-left panel), T = 2 (top-right panel); T = 3 (center-left panel), T = 4 (center-right panel); T = 10 (bottom-left panel), T = 100 (bottom-right panel).

Figure 7 :

 7 Figure 7: Top and center panels: Color representation of the chaos indicator λ T computed for the double gyre defined by = 0.25 on a refined grid of regularly-spaced initial conditions, and different zoom levels. The integration time is T = 10. Bottom-left panel:Representation of λ T computed on a one-dimensional grid of initial conditions having the same value of y; the grid is defined close to the point marked with a cross in the bottom left panel. Bottomright panel: Representation of the two orbits with initial conditions in the onedimensional grid defined for the bottom center panel, and there highlighted by bold points. The two orbits are represented with green color as long as they are so close that cannot be distinguished on the picture (this happens in the time interval [0, 9]); when their separation is sufficiently large one orbit is represented in blue and the other one in red (this happens in the time interval[START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF][START_REF] Guzzo | On the identification of multiple close-encounters in the planar circular restricted three body problem[END_REF]). The separation of the orbits occur within the time interval[START_REF] Conley | Low Energy Transit Orbits in the Restricted Three-Body Problems[END_REF][START_REF] Froeschlé | Graphical Evolution of the Arnold Web: From Order to Chaos[END_REF], consistently with the computation of the chaos indicator for T = 10.

  ) and represent the chaos indicator λ T computed for T = 100, for a grid of initial conditions defined by ϕ(0) = 0 and I(0) ∈[-4, 4]. The intersection of the separatrices with the segment ϕ(0) = 0, I(0) ∈ [-4, 4] appears clearly from the highest values of λ T represented in the top panel of Figure8. Let us zoom on the one-dimensional grid of initial conditions close to a maximum of λ T : we first consider 25 initial conditions with ϕ(0) = 0, I(0) ∈ [1.95, 2.05], and represent their orbits (middle-right panel) as well as the related time evolution of λ T (middle-left panel). Then, we consider a more refined grid of 100 initial conditions with ϕ(0) = 0, I(0) ∈ [1.95, 2.05], and represent their orbits (bottomright panel) as well as the related time evolution of λ T (bottom-left panel).

Figure 8 :

 8 Figure 8: Top-panel: Values of the chaos indicator λ T computed for the mechanical pendulum (7) for T = 100, in a grid of initial conditions defined by ϕ(0) = 0 and I(0) ∈ [-4, 4]. Middle-left panel: Representation of the time evolution of λ T for a grid of 25 initial conditions close to the maximum. Middleright panel: Orbits corresponding to the 25 initial conditions selected for the previous panel. Bottom-left panel: Representation of the time evolution of λ T for a grid of 100 initial conditions close to the maximum. Bottom-right panel: Orbits corresponding to the 100 initial conditions selected for the previous panel. In the middle and bottom panels the light blue color identifies the orbits and the time evolution of λ T corresponding to the initial conditions of the grids with the largest value of λ 100 .

Figure 9 :

 9 Figure 9: Left-panel: Values of the chaos indicator λT defined in (9) computed for the mechanical pendulum (7) for T = 100, in a grid of 100 initial conditions defined by ϕ(0) = 0 and I(0) close to the maximum. Right-panel: Representation of the time evolution of λT for the 100 initial conditions selected for the previous panel.

Figure 10 :

 10 Figure 10: Color representation of the chaos indicator λ T computed for the map (A) (left panels) and the map (B) (right panels) defined in (10) and (11) on a refined grid of regularly-spaced initial conditions, and different integration times (T = 10 for the top panels, T = 15 for the middle panels and T = 25 for the bottom panels).

Figure 11 :

 11 Figure 11: Color representation of the chaos indicator λ T (left panel) and the modified FLI (right panel) computed for the map (B) for T = 20. The modified indicator is defined by the window function u(x) as in (13) with r = 0.05.

Figure 12 :

 12 Figure 12: Representation of examples of intersection of the sets ΠM C defined in eq. (16) with the planes z = 0 (left panel) and y = 0 (right panel) for values of C smaller than C 2 (the shaded areas are in the realm of forbidden motions).

Figure 13 :

 13 Figure 13: Color representation of the (non modified) FLI computed on a grid of initial conditions on the Section Σ P (C) with C = 3.0368... and ẏ > 0; the integration time T is half Jupiter's period (left panel) and twice Jupiter's period (right panel). The left panel is reprinted from: Physica D, vol. 325, Lega E. and Guzzo M., Three-dimensional representations of the tube manifolds of the planar restricted three-body problem, 41-52, Fig. 2, Copyright 2016, with permission from Elsevier.

Figure 14 :

 14 Figure 14: mFLI computation of the unstable manifold of L 1 (C) and the stable manifold of L 2 (C) for the Sun-Jupiter systems and C = 3.0368.... The panels represent the values of the modified FLI (26) computed on the section Σ P (C) (initial conditions with ẏ > 0): on each point of the panel we represent its mFLI value with a color scale, so that both the unstable manifold of L 1 (C) and the stable manifold of L 2 (C) appear as the yellow curves. The left panel is obtained for T 1 = T 2 = 5, the right panel with T 1 = T 2 = 100. On the left-panel we also identify with a black square and with a black circle a selected set of curves belonging to the unstable manifold of L 1 (C) or the stable manifold of L 2 (C) respectively. The left panel is reprinted from: Guzzo M., Lega E., Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Lyapunov Indicators, SIAM J. APPL. MATH., Vol. 74, No. 4, pp. 1058-1086, 2014, Fig. 5; Copyright 2014 Society for Industrial and Applied Mathematics; Reprinted with permission. The figure is here represented with some graphic filters in order to appreciate more clearly the manifolds.

Figure 15 :

 15 Figure 15: On the top panel we report the representation of the indicator (26) computed with T 1 = 5, and select the heteroclinic point in the black box to compute with higher precision. In the bottom-left panel we represent again the indicator (26) computed on an extremely small box (of side 10 -15 ) around the heteroclinic point selected in the top panel (the values of x, ẋ corresponding to the heteroclinic point are denoted with x h , ẋh in the labels). In such a way we identify the heteroclinic point with a precision of more than 16 digits. On the bottom-right panel we report the projection on the plane (x, y) of the orbit with initial conditions corresponding to the numerically computed heteroclinic point. The panels are reprinted from: Guzzo M., Lega E., Evolution of the tangent vectors and localization of the stable and unstable manifolds of hyperbolic orbits by Fast Lyapunov Indicators, SIAM J. APPL. MATH., Vol. 74, No. 4, pp. 1058-1086, 2014, Fig. 4 and Fig. 6; Copyright 2014 Society for Industrial and Applied Mathematics; Reprinted with permission.

Figure 17 :

 17 Figure 17: The panels represent with different zoom levels six ridges denoted µ 1 , µ 2 (on the left of L 1 ) and λ 1 , λ 2 , λ 3 , λ 4 of the mFLI (25) computed for the Sun-Jupiter system on the planar section Σ P (C), for C = 3.0368.... The target orbit γ in the definition of the mFLI is the Lyapunov orbit L 1 (C). The red curves represent the projection on the section Σ P (C) of sample orbits connecting two ridges; the arrows indicate the positive time-direction on the orbits. Top and bottom left panels reprinted from: Physica D, vol. 325, Lega E. and Guzzo M., Three-dimensional representations of the tube manifolds of the planar restricted three-body problem, 41-52, Fig. 8, Copyright 2016, with permission from Elsevier. Bottom right panel reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 2, Copyright 2018, with permission from Elsevier.

Figure 18 :

 18 Figure 18: Projection on the xy plane of orbits with initial conditions on λ 2 , numerically computed forward (red curves) and backward (blue curves) in time. The bold red curves represent the Lyapunov orbits L 1 (C), L 2 (C); the bold blue curve connect the endpoints of the blue curves; the shaded area represents the realm of forbidden motions.Figure reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 3, Copyright 2018, with permission from Elsevier.

Figure 20 :

 20 Figure 20: Top-left panel: representation on the x, ẋ plane of some connected components of the intersection of both stable and unstable manifolds of the Lyapunov orbits L 1 (C) (the red curves) and L 2 (C) (the blue curves) with the section Σ P (C), computed as the ridges of the mFLI. The bold curves highlight initial conditions which are numerically integrated to provide the top-right and the bottom panels. Top-right panel: the violet surface represents λ 2 ; the blue surface represents the forward time evolution of the part of W u L2(C) highlighted in bold black in the top-left panel; the yellow surface represents the backward time evolution of the part of W u L2(C) highlighted in bold black in the top-left panel; the light-blue surface represents the forward time evolution of the part of W s L1(C) highlighted in bold black in the top-left panel; the transparent mesh represents the boundary of the admissible region in the space x, y, ẋ. Bottom panel: we represent with pink color also the backward time evolution of the part of W s L1(C) . Figure reprinted from: Physica D, vol. 325, Lega E. and Guzzo M., Three-dimensional representations of the tube manifolds of the planar restricted three-body problem, 41-52, Fig. 14 and 15, Copyright 2016, with permission from Elsevier.

Figure 21 :

 21 Figure 21: Top panel: Ridges of the mFLI defined in (29), computed in the vertical section Σ V (C), representing connected components of the intersection of both stable and unstable manifolds of CC 1 and CC 2 with the section. Bottom panels: The ridges in the red box are magnified in the bottom-left panel. In the bottom-right panel we represent in the x, ẋ, z space the ridge λ 4 , computed for the planar problem (red curve) and the two ridges λ 4 computed for the spatial problem (blue curves). We appreciate the intersection of the spatial curves with the planar one along the line z, ẋ = 0. Panels reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 6 and 8, Copyright 2018, with permission from Elsevier.

Figure 22 :

 22 Figure 22: Representation of the projection of the orbits associated to solutions with initial conditions in the lobe λ 2 of Figure 21 on the x, y plane (top-left panel), on the x, z plane (top-right panel) and on the x, y, z space (bottom panel). Figures reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 11, Copyright 2018, with permission from Elsevier.

Figure 23 :

 23 Figure 23: Projection on the xy plane of orbits with initial conditions in the inner border (left panel) and in the outer border (right panel) of the annulus σ ⊆ Σ P (C) containing the ridge λ 2 of the planar CR3BP. Panels reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 13, Copyright 2018, with permission from Elsevier.

  Figure24: For any initial condition on a refined grid of 10000 × 2000 initial conditions of Σ V (C) we computed their orbits until they reach a value for the x coordinate smaller than 0.933 or larger than 1.067, within the maximum integration time T = 10. In the former case, we represent the initial condition with a violet pixel, in the latter we represent the initial condition with an orange pixel. If none of the conditions is met, we represent the initial condition with a yellow pixel. Figure reprinted from: Physica D, vol. 373, Guzzo M. and Lega E., Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem, 38-58, Fig. 14, Copyright 2018, with permission from Elsevier.

The analog definition of Lyapunov characteristic exponents for the dynamics defined by the iterations of a map will be given in Section

2.1.

For the Sun-Jupiter mass ratio we have C 1 =

3.0387..., C 2 = 3.0374...; for the numerical experiments we use C = 3.0368.... We remark that the values of the Jacobi constant C 1 , C 2 are very close, and in particular small changes of C close to C 1 , C 2 determine important changes in the shape of the realm of possible motions.

The definition depends on the local inversion map used to choose ζ(0): the indicators (29) are therefore locally defined; as a matter of fact it is sufficient to use a couple of indicators mFLI ± φ (ξ, w(0), T ) to cover the phase-space, see[START_REF] Guzzo | Geometric chaos indicators and computations of the spherical hypertube manifolds of the spatial circular restricted three-body problem[END_REF].

Acknowledgements. M.G. acknowledges the project MIUR-PRIN 20178CJA2B "New frontiers of Celestial Mechanics: theory and applications".