Computing invariant manifolds for libration point missions

The goal of this lecture is to review several methodologies for the computation of invariant manifolds, having in mind the needs of preliminary mission design of libration point missions. Because of this, the methods reviewed are developed for and applied to the circular, spatial restricted three-body problem (RTBP), although most of them can be applied with few changes, or almost none, to general dynamical systems. The methodology reviewed covers the computation of (families of) fixed points, periodic orbits, and invariant tori, together with the stable and unstable manifolds of all these kinds of invariant objects, and also homoclinic and heteroclinic connections between them. The methods reviewed include purely numerical and semi-analytical ones. No background is assumed except for a graduate level knowledge of calculus, differential equations and basic numerical methods. In particular, the notions from the theory of dynamical systems required for the development of the methods are introduced as needed.


Introduction
In libration point missions, spacecraft are sent to orbits that stay close to the fixed points of the circular, spatial, restricted three-body problem (RTBP) with primaries the Sun and a planet, or a planet and a moon. The RTBP model describes the motion of an infinitesimal particle under the attraction of two massive bodies known as primaries, that are assumed to revolve uniformly in circles around their center of mass. In rotating coordinates, this model has five equilibrium points: three of them, L 1 , L 2 , L 3 , also known as collinear, were discovered by Euler, and two more, L 4 , L 5 , also known as triangular, were discovered by Lagrange. Compared to orbits around the Earth or other planets, orbits around the collinear libration points provide ideal locations for space observation. Among their advantages are the absence of shadow of a celestial body, thus providing a more stable thermal environment, and continuous access to the whole celestial sphere, except for a direction, that is not fixed but rotates with the primaries. Also, the instability of the collinear libration points gives rise to a very rich dynamical structure, that can be exploited not only to search for operational nominal orbits but also to find low-energy passageways between them. These operational orbits could be either of the libration point type or around celestial bodies.
Four examples of libration point missions of different kinds are: SOHO, launched in Dec. 1995, to an Halo orbit around the collinear point L 1 of the Earth-Sun system. Its goal is to provide continuous observations of the Sun, and is still operational.
WMAP, launched in June 2001, to a Lissajous orbit around the collinear point L 2 of the Earth-Sun system. Its goal was to map the temperature fluctuations of the cosmic microwave radiation.
Genesis, launched in Aug. 2001, to an Halo orbit around the collinear point L 1 of the Earth-Sun system. Its goal was to collect solar wind samples and deliver them to Earth in daylight. For this, an additional excursion close to the collinear point L 2 of the Earth-Sun system was necessary.
Artemis, started in Jan. 2009 as an extension of the mission of two of the spacecrafts of the Themis mission, that, using the remaining fuel, were sent from high, eccentric Earth orbits to lunar orbits using L 1 and L 2 Earth-Moon dynamics.
Illustrations of the trajectories of these four missions are shown in Fig. 1. The nominal trajectories of these four missions can be identified among the families of periodic orbits and invariant tori related to the collinear libration points of the RTBP. In the case of SOHO and Genesis, the nominal trajectory is part of the Halo family of periodic orbits. In the case of WMAP, it is part of the Lissajous family of invariant tori. In the case of Artemis, the nominal trajectories would be the final lunar ones, but invariant tori of the L 1 and L 2 Lissajous family play a fundamental role in the transfer from Earth to lunar orbits. The invariant stable (resp. unstable) manifolds of all these periodic orbits and tori can be used to arrive to (resp. depart from) them. In the case of the P1 spacecraft of Artemis, an heteroclinic  connection is closely followed in order to go from the Lissajous torus around L 2 to the Lissajous torus around L 1 . Such connections are obtained as intersections of the stable manifold of the arriving object and the unstable manifold of the departing object. An heteroclinic connection is also outlined by the Genesis mission. The preliminary mission design of these kind of missions is based in being able to compute families of trajectories, in order to be able to select the one that best satisfies the requirements of the mission. The goal of this lecture is to review some of the numerical and semianalytical techniques available in the literature that can be used in order to compute families of periodic orbits and invariant 2D tori of the RTBP as a dynamical system, as well as their invariant stable and unstable manifolds. Some discussion will be done also on the computation of homoclinic and heteroclinic connections. Although preliminary mission design is the main motivation for this lecture, all the techniques that will be described can be used for the numerical computation of invariant objects in other conservative dynamical systems. Many of them can be directly used in or easily adapted to the dissipative case as well.
The methods that will be described can be divided in two classes: numerical and semi-analytical. Semi-analytical methods provide expansions around a base object that must be previously known. They are more convenient than numerical methods for parametric studies of trajectory features, since a single expansion covers a family or many families of trajectories. They have as a drawback that the expansions provide a good approximation of the dynamics in a neighborhood of the base object, but not outside this neighborhood. Numerical methods are able to compute individual objects anywhere in phase space, but parametric studies with them are more tedious, since they require the previous generation of a large database of trajectories obtained by numerical continuation of one or many families of trajectories. This does not mean that parametric studies are not feasible: an example of systematic continuation families of periodic orbits and invariant tori in order to cover large regions of phase space is given in subsection 3.9.
The numerical method chosen for the computation of the 2D tori of the RTBP is based on looking for the Fourier series of a curve in the torus invariant by the time-T flow, being T one of the periods of the quasi-periodic trajectories inside the torus [6,16]. It is a wellestablished method that has proven to be among the most adequate in this context (see [3] for a review of several methods). Since its use requires starting from the normal part of periodic orbits that need to be previously obtained by continuation, this lecture also includes a discussion on the numerical computation of fixed points and periodic orbits, and develops explicit formulation of the linear approximation of their normal dynamics. On the semi-analytical side, this lecture will cover a technique based on the parameterization method [20,21], that produces Taylor expansions of the center manifold of a collinear point and the corresponding reduction of the vector field. In the reduced field, the collinear libration point is no longer unstable, so the trajectories in a neighborhood of it can be obtained by direct numerical integration. An earlier technique known as reduction to the center manifold [13,25] produces the same results. There is another technique, known as the Lindstedt-Poincaré method [25,29], which is still more convenient for parametric studies because it produces expansions of the trajectories instead of the center manifold, at the expense of a slightly smaller neighborhood of validity of the expansions.
The lecture is structured as follows. Section 2 reviews some of the common nomenclature in dynamical systems and, in doing so, introduces the relevant features of the RTBP. After that, section 3 describes numerical techniques for the computation of periodic orbits and invariant tori, whereas section 4 explains how to compute the same objects semi-analytically using the parameterization method. Attention is then focused on the computation of the stable and unstable manifolds of periodic orbits and tori. Section 5 reviews numerical techniques to compute their linear approximation, whereas section 6 explains how to obtain trajectories in these manifolds semi-analytically via the parameterization method. Finally, section 7 addresses the computation and continuation of homoclinic and heteroclinic connections.

Dynamical systems and the RTBP
This section recalls some notions from the theory of dynamical systems and also introduces the circular, spatial Restricted Three-Body Problem (RTBP). Although most readers will probably be familiar with these notions, recalling them will allow us to introduce notations that will be used in the rest of the lecture.

Continuous dynamical systems
The theory of dynamical systems provides an abstract framework for the mathematical study of systems that evolve with time in a deterministic manner. Continuous dynamical systems are those in which time is considered a continuous variable, this is, t ∈ R. They are usually defined in terms of a system of autonomous (time-independent) Ordinary Differential Equations (ODE)         ẋ 1 = X 1 (x 1 , x 2 , . . . , x n ), x 2 = X 2 (x 1 , x 2 , . . . , x n ), . . .
Assuming that this system of ODE can be integrated for all time, for t ∈ R the time-t flow, φ t : R n −→ R n , is defined by the initial value problem It can be thought as a map that "flows" initial conditions along the corresponding trajectories for t time units. The subscript notation for t is in order to stress this fact. It is also common to refer to a continuous dynamical system as "a flow". Given an initial condition x 0 , the corresponding orbit is {φ t (x 0 )} t∈R . A fixed point of a continuous dynamical system is a point whose orbit is itself, that is, φ t (x) = x, ∀t ∈ R. This can only happen if f (x) = 0. An orbit {φ t (x)} t∈R is said to be periodic if there is T > 0 such that φ T (x) = x, φ t (x) = x, for 0 < t < T.
Then T is said to be its period. A set of initial conditions A ⊂ R n is said to be an invariant set if φ t (x) ∈ A ∀t ∈ R, ∀x ∈ A.
A straightforward example is an orbit (in particular, a fixed point or a periodic orbit).
A manifold is a set of points defined (maybe piece-wise) by equations, either implicit or parametric. An invariant manifold is an invariant set that is a manifold (e.g. a torus). We will usually speak of general invariant manifolds as invariant objects and reserve "invariant manifold" to denote stable, unstable or center manifolds associated to an invariant object. Given an invariant set A, its stable set (resp. unstable set) is the set W s (A) (resp. W u (A)) of initial conditions that tend to the object asymptotically through the flow forward (resp. backward) in time. In other words, the set of initial conditions that approach (resp. depart from) A. This is, For several cases in which A is a manifold (e.g. a fixed point, a periodic orbit, an invariant torus), W s (A) (resp. W u (A)) is also a manifold, and is called the stable manifold (resp. unstable manifold) of A.

The circular, spatial Restricted Three-Body Problem
The circular, spatial restricted three-body problem (RTBP) is an example of continuous dynamical system. It can be written as a Hamiltonian system with three degrees of freedom (details on the theory of Hamiltonian systems can be found in e.g. [31]) with Hamiltonian The system of ODE that defines it is, therefore, The RTBP describes the motion of a massless particle ("massless" in the sense that it is considered not to influence gravitationally any other body) under the gravitational attraction of two bodies, called primaries, with masses m 1 > m 2 . The primaries are assumed to revolve uniformly in circles around their common center of mass. The coordinate system used is a synodic one, that rotates with the primaries so that the primary of mass m 1 is fixed at (µ, 0, 0, 0, µ, 0), and the primary of mass m 2 is fixed at (µ − 1, 0, 0, 0, µ − 1, 0). The RTBP depends on the mass parameter µ = m 2 /(m 1 + m 2 ). As it is common with Hamiltonian systems, the (x, y, z) coordinates are called positions, and the (p x , p y , p z ) coordinates are called momenta. The space of positions (this is, 3D physical space) is called configuration space. See Fig. 2. In short, the RTBP can be denoted asẋ = X(x), with (2)

Discrete dynamical systems
Discrete dynamical systems are those in which time is considered as a discrete variable, this is, t ∈ Z. They are are defined by diffeomorphisms (smooth 1-1 maps) We denote by F −1 the inverse map of F , and use superscript notation for the composition of maps: In this way, F n is "the discrete time-n flow". Via this notion, all the previous notions from continuous dynamical systems translate to the discrete case. Given an initial condition, its related orbit is the set A fixed point is an initial condition such that its orbit is itself, If A is a manifold, the stable and unstable sets of A, defined by are usually manifolds. An paradigmatic example of a discrete dynamical system is Chirikov's standard map, that, in one of its formulations is written as Here a is a parameter and x, y ∈ T = R/[0, 2π], that is, we assume that (x, y) and (x,ȳ) are the same point if x −x = j2π, y −ȳ = k2π, for j, k ∈ Z. The standard map is area-preserving. In two dimensions, being area preserving is equivalent to being sympletic, which is the discrete analog of the Hamiltonian formalism (for more details see any textbook in the subject, e.g. [31]). The global dynamics (a phase portrait) of two-dimensional area-preserving maps on compact regions can be swiftly obtained by iteration of the map F . Fig. 3 is obtained by considering the initial conditions {p j := (−π +j2π/100, 0)} 100 j=0 and plotting the points {F k (p j ))} 1000 k=0 , for j = 0, . . . , 100. Several kinds of invariant sets (that are manifolds) can be found, like fixed points, periodic points of different periods and invariant curves. Invariant sets with chaotic dynamics can also be observed.

Orbit generation in a dynamical system
Orbits in discrete dynamical systems can be generated just by iteration of the map, as it has been done in figure 3. In continuous dynamical systems, numerical methods for integration of ODE have to be used. In order to have error control, variable-step methods are preferred over constant-step ones. A popular family of variable-step methods are the Runge-Kutta-Fehlberg (RKF) ones, of which there are some high-order versions like RKF78 [10]. There are many alternatives (see e.g. [41,19]). In the case of a system of ODE given by closed formulae, like the RTBP, a particularly good choice is Taylor's method, of which there are freely available implementations [27,1]. Here we will discuss briefly the black-box usage of a one-step method with step size control.
For a system of n possibly non-autonomous ODE, with x, X(t, x) ∈ R n , denote by φ(t, t 0 , x) its flow from time t 0 to time t, defined by the conditions Given t 0 ∈ R, x 0 ∈ R n , h 0 ∈ R (small), and a tolerance δ, a routine implementing a one-step method with step size control will return (c) h 1 is a recommended step length for the next call.
In the algorithmic descriptions that will follow, we will denote a call to such a routine as In order to implement φ(t 1 , t 0 , x 0 ) for arbitrary t 1 , t 0 , x 0 , it is necessary to write a routine that calls ODEstep many times using as input step h 0 the recommended step h 1 of the previous call, plus a final call with h 0 the step needed to in order to reach the final time t 1 (or more than one such calls, if the step given is reduced by the step size control). In some of the algorithmic descriptions that will follow, a call to such a routine will be denoted as In the systems of non-linear equations that we will solve in order to compute invariant objects, we will need to be able to numerically evaluate both the flow and its differential with respect to initial conditions, that we will denote as D x φ(t, t 0 , x), or simply Dφ t (x) in the autonomous case. It can be found by numerically integrating the system of ODE together with its first variational equations: where x is an n-dimensional vector and A is a n × n matrix. If x(t) and A(t) are solutions of (5) with x(t 0 ) = x 0 and A(t 0 ) = I n the n × n identity matrix, then D x φ(t, t 0 , x 0 ) = A(t). System (5) can be written as a system of n + n 2 ODE aṡ . . , x n ) a k,j , i, j = 1, . . . , n.

Poincaré maps
A way to simplify the study of a continuous dynamical systems is to consider a discrete dynamical system that has essentially the same dynamics. One way to do it is, for a fixed T > 0, to consider the time-T flow (or stroboscopic map), φ T , which is a discrete dynamical system. In this way, for instance, T -periodic orbits are turned into fixed points. Another way to turn a continuous dynamical system into a discrete one is through a Poincaré map. For a continuous dynamical system given byẋ = X(x), let Σ be a hypersurface of R n , and assume it is transversal to the vector field, that is, X(x) is not tangent to Σ for all x ∈ Σ. Let x 0 be such that φ T 0 (x 0 ) ∈ Σ for some T 0 > 0. The implicit function theorem ensures the existence of a neighborhood U x 0 and a map τ : U → R, known as time-return map, such that τ (x 0 ) = T 0 and The map P (x) := φ τ (x) (x) is called Poincaré map corresponding to Σ. If x 0 ∈ Σ and P (Σ ∩ U ) = Σ ∩ U , the restriction of P to Σ ∩ U defines a discrete dynamical system. In going to the starting continuous dynamical system to the discrete one defined by P , periodic orbits are turned into fixed points, and invariant tori are turned into invariant curves. In general, invariant objects lose one dimension, which is an advantage both from the theoretical and the computational point of view.
Orbit generation in this discrete dynamical system requires the numerical evaluation of a Poincaré map, which has as a difficulty that the time-return map τ (x) is unknown. It can be adjusted by Newton iterations once we get close to the section Σ. This is done in Algorithm 1.
If the differential of the Poincaré map, DP (x), is also needed, it can be computed as where Dφ τ (x) (x) is to be understood as D y φ τ (x) (y)| y=x . An expression for Dτ (x) can be obtained by implicit differentiation on g(P (x)) ≡ 0. After substitution in the previous equation, In a routine implementing Algorithm 1 for the evaluation of a Poincaré map, it is useful to keep as an option the integration of the system of ODE defining our continuous dynamical system together with its first variational equations (5), in order to have available Dφ τ (x) (x) to be used in (6).

Numerical computation of periodic orbits and 2D tori
The goal of this section is to review some numerical methods for the numerical computation of the periodic orbits and invariant 2D tori related to collinear libration points. Since the RTBP is Hamiltonian, both periodic orbits and tori are not isolated but embedded in families. Once an invariant object (periodic orbit or torus) has been computed, the remaining objects of its family can be obtained by continuation. The first object of a family is usually computed from the linear approximation of the dynamics around a simpler object (e.g. a torus from the linear dynamics around a periodic orbit, or a periodic orbit from the linear dynamics around a fixed point). This approach can be followed hierarchically in order to do a systematic study of the dynamics around a collinear point.
This section starts recalling the predictor-corrector or pseudo-arclength continuation method as described by standard references (e.g. [2]). After that, subsection 3.2 provides an strategy for the numerical solution of not necessarily square non-linear systems of equations, that simplifies considerably the practical implementation of the methods described later. The subsections that follow (from 3.3 to 3.8) provide methods for the computation of invariant objects and formulation for the linear dynamics around them, necessary to implement a systematic numerical exploration of the dynamics around a collinear libration point. This is actually done in subsection 3.9 for the L 1 collinear point of the Earth-Moon RTBP.

Numerical continuation
A classical way to introduce numerical continuation is as a technique to find a (unknown) solution of a system of non-linear equations G(x) = 0 from a known solution of another system F (x) = 0, that is close to G(x) = 0 in some sense. In order to look for a zero of G, a one-parametric family of intermediate systems H(λ, x) is considered with H(0, x) = F (x) and H(1, x) = G(x). For instance, the convex homotopy between F and G, Then we can try to continue the known solution x 0 of H(0, x) = 0 up to a solution of H(1, x) = 0 with respect to the parameter λ. The algorithm below provides a straightforward approach.
Algorithm 2 Continuation of H(λ, x) = 0 with respect to the parameter λ. input: x 0 ∈ R n such that H(0, x 0 ) = 0, m ∈ N do: x := x 0 ∆λ := 1/m ∀i = 1, . . . , m λ := i∆λ solve H(λ, y) = 0 iteratively for y taking x as starting value x := y output: x Algorithm 2 breaks down if there is a turning point with respect to λ along the continuation curve. An alternative that can cope with this case is the predictor-corrector or pseudo-arclength continuation method (see e.g. [2]). Its basic idea is to consider arclength instead of λ as the continuation parameter. "Pseudo" stands for the fact that the actual parameter is not truly arclength but distance along a line tangent to the continued curve. Define y = (λ, x) ∈ R n+1 . Then H(y) := H(λ, x) = 0 defines implicitly a curve in R n+1 as long as rank DH(y) = n, which is a condition we will assume. The continuation can be done according to the algorithm stated next.

input:
y = (λ, x) ∈ R n+1 such that Π 1 y := λ = 0, H(y) = 0. do: while (Π 1 y < 1) let v ∈ ker DH(y), v 2 = 1, pointing ahead take z := y + γv, for suitable γ (see the comments below) if (Π 1 z < 1) solve H(z) = 0 iteratively for z by a modified Newton's method taking minimum-norm corrections else γ := (1 − Π 1 y)/Π 1 v z := y + γv solve H(z) = 0 by Newton iterations keeping Π 1 z constant y := z output: y A convenient way to control the step length γ is in order to keep constant the number of Newton iterations when solving H(z) = 0 at each continuation step. A simple rule is to assume that this number of iterations is a linear function of the step length chosen: if n old is the number of iterations performed in the last continuation step, γ old is the last step length used and n des is the desired number of Newton iterations, we can take γ = n des n old γ old .
Note that, except for the start and stop criteria, in the pseudo-arclength method there is no distinguished coordinate to be considered a parameter. It can therefore be applied to any system of non-linear equations H(y) = 0, as long as its solution is a curve.

Numerical solution of non-square, non-linear systems of equations
In subsections 3.6 and 3.8, the computation of periodic orbits and invariant 2D tori will be done in terms of solving non-linear systems of equations. In the case of the computation of a single object, the system to be solved will have (locally) unique solution. It is standard practice in this case to require such a system to be square, this is, of the form G(y) = 0 with G : R N −→ R N for some N , and to use Newton's method (see any textbook on numerical analysis, e.g. [41]). In the case of the continuation of a family, the system to be solved will not have unique solution but a curve of solutions. It is standard practice in this case to require such a system to have one more unknown than equations, this is, to be of the form G(y) = 0 with G : R N +1 → R N for some N , and to use Newton's method with some modification to account for non-uniqueness (see e.g. [2,39]). In order to keep the systems of equations of subsections 3.6 and 3.8 simple, it will be convenient not to require them to be either N × N or (N + 1) × N . A way to be able to solve these systems is to consider a modified Newton method y n+1 = y n − (∆y) n in which the linear system to be solved for the correction, DG(y n )(∆y) n = G(y n ), is solved for its minimum-norm, least-squares solution. The minimumnorm, least-squares solution always exists and is unique for any linear system of equations, square or not. Assuming that the starting nonlinear system G(y) = 0 has a solution (perhaps non-unique) and that the initial guess y 0 is close to a solution, this strategy will converge to a nearby solution using minimum-norm corrections.
We discuss briefly how to compute the minimum-norm, least-squares solution of a linear system Ay = b using QR decomposition with column pivoting. 1 Assume that A is an arbitrary m × n matrix with r := rank A ≤ min(m, n). A least squares solution of Ay = b, always exists. If r = n, there is an unique least-squares solution. If r < n, there is a linear subspace of least-squares solutions of dimension d := n − r. Nevertheless, as mentioned previously, is always unique. By applying to A Householder transformations with column pivoting [12], we obtain a decomposition with Q an m × m orthogonal matrix, R 11 an r × r upper-triangular matrix with non-zero diagonal elements, and P a n × n permutation matrix. In order to perform this decomposition, r (or, equivalently, d = n − r) must be known. If we denote where I d is the d × d identity matrix. Finding the minimum-norm element of the previous set is an standard full-rank least-squares problem, that can be solved via a standard (without column pivoting) QR decomposition.
In order to solve the systems of equations of subsections 3.6 and 3.8, it is convenient to write a routine that, for a general general m × n linear system of equations, finds the minimum-norm least-squares solution and, optionally, a basis of the kernel of A, which can be obtained from

Computation of fixed points of flows and maps
For the computation of a fixed point of a flowẋ = X(x), we look for p such that G(p) := X(p) = 0. For the computation of a fixed point of a map x → F (x), we look for p such that G(p) := F (p) − p = 0.
In any case, we can use Newton's method in several variables in order to look for a zero of G.
Algorithm 4 Newton's method in order to find a zero of a function G : R n → R n with tolerance tol, starting from a first guess p 0 , allowing for a maximum of maxit iterations.
input: p 0 , G, tol, maxit do: In the RTBP, it can be analytically seen (see e.g. [42]) that the distance from L j , j = 1, 2, 3 to the closest primary, that will be denoted γ j , is given by the only positive root of the corresponding Euler's quintic equation: Therefore, in this case it is enough to use Newton's method in one dimension. Good guesses are (µ/3) 1/3 for L 1,2 and 1 − (7/12)µ for L 3 .

Linear behavior around fixed points of flows
For a flowẋ = X(x) with a fixed point p, since X(x) = X(p) + DX(p)(x − p) + O( x − p 2 ) and X(p) = 0, its linearization around with A := DX(p). The eigenvalues of A are known as the exponents of the fixed point p.
Assume that λ ∈ R, λ = 0 is an eigenvalue of A, and v is a corresponding eigenvector. Then is a solution of the linearized flow (8). If λ < 0, ϕ(t) → p as t → +∞, so {ϕ(t)} t∈R is a trajectory in the stable manifold of p in the linearized flow. If λ > 0, ϕ(t) → p as t → −∞, so {ϕ(t)} t∈R is a trajectory in the unstable manifold of p in the linearized flow. The stable manifold theorem for flows (see e.g. [18,35]) ensures the existence of a stable (resp. unstable) manifold of the full non-linear flowẋ = X(x) that contains p and is tangent to the linear subspace spanned by the eigenvectors of A corresponding to eigenvalues with strictly negative (resp. positive) real part.
Assume now λ = iω for ω ∈ R, ω = 0, where i denotes the imaginary unit. In this case, −iω is also an eigenvalue, so we can assume that ω > 0. Let v 1 + iv 2 be a corresponding eigenvector, with v 1 , v 2 ∈ R n . Define By using Av 1 = −ωv 2 and Av 2 = ωv 1 (that follows from A(v 1 +iv 2 ) = iω(v 1 + iv 2 )), it is seen that ϕ γ (t) is a solution of the linearized flow. Therefore, by varying γ, ϕ γ (t) provides a one-parametric family of periodic orbits of period 2π/ω of the linearized flow. If the remaining eigenvalues λ j of A satisfy that λ j /(iω) is not an integer, then Lyapunov's center theorem (see e.g. [31,38]) ensures the existence of a one-parametric family of periodic orbits of the non-linear flow with periods that tend to 2π/ω as the periodic orbits collapse to p. These periodic orbits are part of the center manifold of p, which is an invariant manifold tangent to the linear subspace spanned by eigenvectors corresponding to eigenvalues of A with zero real part. The existence of the this manifold is ensured by the center manifold theorem for flows (see e.g. [18,35]). Expressions for trajectories of the linear flow in the case λ = a + iω a, ω ∈ R, a, ω = 0, can be obtained similarly. They will not be necessary in what follows. These trajectories would be close to trajectories of the non-linear flow in the stable, unstable or center manifold according to whether a < 0, a > 0 or a = 0, respectively.

Linear behaviour around fixed points of maps
For a discrete dynamical system given by with A = DF (p). The eigenvalues of A are also called the multipliers of p.
The stable manifold theorem for maps (see e.g. [18,35]) ensures the existence of a stable (resp. unstable) manifold of the full non-linear map F that contains p and is tangent to the linear subspace spanned by the eigenvectors of A corresponding to eigenvalues with modulus strictly smaller (resp. larger) than one.
Assume now λ ∈ C, |λ| = 1, λ = cos ρ + i sin ρ and let v 1 + iv 2 be an associated eigenvector, with v 1 , v 2 ∈ R n . Then is an invariant closed curve of the linearized map. Therefore, by varying γ, ϕ γ (θ) provides a one-parametric family of invariant curves of the linearized map with rotation number ρ. Under number-theoretical hypotheses of ρ and non-degeneracy ones of F , KAM theory (see e.g. [26]) ensures the existence of a Cantorian 2 one-parametric family of invariant curves of the full non-linear map F , with rotation numbers that tend to ρ as the invariant curves collapse to p.

Computation of periodic orbits of flows
The computation of periodic orbits is a classical and well-known subject. There are publicly available software packages, like AUTO-07p [9], that are capable of both computing individual periodic orbits and performing continuation of families. Nevertheless, the simplicity of the methodology that will follow makes its implementation worthwhile, both for computational efficiency and for easier interaction with the methods of computation of invariant tori of section 3.8. We discuss in this section how to compute initial conditions for periodic orbits by solving non-linear systems of equations stated in terms of the flow. The discussion will partially follow [39]. An initial condition for a periodic orbit of a flow can be thought as a fixed point of a discrete dynamical system. In order to turn this idea into a numerical method, consider first a non-autonomous T -periodic system of n ODE,ẋ = X(ωt, x), with ω = 2π/T and X(θ, x) 2π-periodic in θ. 3 Denote its flow by φ(t, t 0 , x 0 ), defined as in (4). Consider the map F (x) := φ(t 0 + T, t 0 , x), with t 0 fixed. An initial condition for a T -periodic orbit of (12) is a fixed point of the discrete dynamical system defined by F , which is found as a zero of G(x) := F (x) − x, as discussed in subsection 3.3. The differential of φ(t 0 + T, t 0 , x) with respect to x is computed by integrating the first variational equations, as discussed in subsection 2.4. Assume now that we have an autonomous system of ODĖ with flow φ t , defined as in subsection 2.1. If we wanted to apply the previous approach, we would look for a fixed point of the discrete dynamical system defined by defines the whole periodic orbit as a manifold, DG(x) is singular at every point of the periodic orbit. We could still use the modified Newton strategy of subsection 3.2, but that would introduce difficulties for continuation. 4 A better strategy is to get rid of the singularity by considering a different discrete dynamical system: a Poincaré map. If Σ is a surface of section transversal to the flow and intersected by the periodic orbit we are looking for, denote as P (x) = φ τ (x) (x) the corresponding Poincaré map, where τ (x) is the time-return map (see subsection 2.5). Then, by looking for a fixed point of P as a zero of G(x) := P (x) − x, we would be looking for an initial condition of the periodic orbit in the Poincaré section Σ, which is locally unique. The previous approach works as long as the periodic orbit we are looking for is isolated, which is usual in generic dynamical systems. But in Hamiltonian systems like the RTBP, periodic orbits are embedded in families. Assume that we are given a Hamiltonian continuous dynamical system with Hamiltonian H(x). The intersections of the periodic orbits of a family with a Poincaré section Σ define locally a curve. On all the points of this curve, DG(x) is singular. A way to get rid of this singularity would be to first reduce our starting dynamical system (13) to an energy manifold {H(x) = h}. Then, an initial condition of a periodic orbit (of energy h) as a fixed point of P would be locally unique. Nevertheless, instead of modifying (13), it is simpler to add an energy equation to the fixed point condition on the Poincaré map. In this way, we would solve for x the non-linear system This system is not square, so a standard approach using Newton's method would not work, but it can be solved by the modified Newton approach of subsection 3.2 with d = 0. In doing this, P and DP can be evaluated as discussed in subsection 2.5.

Practical implementation
An strategy for the computation of periodic orbits still simpler to implement than the one just discussed is to add the Poincaré section as an additional equation. In this way, at the cost of one extra equation, the evaluation of P and DP is avoided. Assuming that the Poincaré section is Σ = {g(x) = 0}, we would solve the (n + 2) × (n + 1) system for (T, x). This can be done using the modified Newton strategy of subsection 3.2 with d = 0, y = (T, x) and An additional advantage of this approach is that the period of the periodic orbit appears explicitly. The system of equations (14) can also be used for the continuation of a family of periodic orbits. This would be done using Algorithm 3 with H := G defined as in (15) but for y = (h, T, x). In an implementation of Algorithm 3, the routine proposed at the end of subsection 3.2 would be able to compute ker DH(y) and to solve H(y) = 0 with minimum-norm corrections.
From the system of equations (14), other systems of interest for the computation and continuation of periodic orbits can be obtained by eliminating equations and unknowns. For instance, if we eliminate the unknown h, keep T constant and eliminate the energy equation H(x) = 0, the resulting system of equations, that is to be solved only for x, would allow us to compute a periodic orbit of a given period. A routine implementing the evaluation of G(y) (as defined in (15)) and DG(y) can also be used in order to solve a system like (16) by giving it the option of eliminating components of G(y) and files and/or columns of DG(y).

Multiple shooting
As we will see later (e.g. in figure 4), for many periodic orbits of the neighborhood of the collinear points of the Earth-Moon RTBP, the maximum absolute value of the eigenvalues of Dφ T (x) can be larger than 2000. This means that, after numerical integration for T time units, any error in the initial condition can be amplified by this factor. Even with exact data, the local truncation error of the first step of numerical integration could be amplified by this factor. 5 Then, for example, if the tolerance of numerical integration is set to 10 −14 , we cannot expect an error smaller than 10 −11 . Because of this, initial conditions for Newton's method need to be very accurate in order to obtain convergence, and continuation steps become very small. We can reduce these amplification factors by making use of multiple shooting. Multiple shooting is classically introduced as a way to overcome dynamical instability in the solution of boundary value problems (see e.g. [41]). As a general idea, the multiple shooting strategy can be thought as introducing intermediate objects as unknowns in order to reduce integration time. In our case, we would need to consider points x 0 := x, x 1 , . . . , x m−1 along the periodic orbit and add the corresponding matching equations to the system to be solved. In this way, system (14) would become In order to compute a single periodic orbit, the unknowns to consider would be (T, x 0 , . . . , x m−1 ). In order to continue a family of periodic orbits, the unknowns would be (h, T, x 0 , . . . , x m−1 ). As commented before, other systems of interest can be obtained from this one by eliminating equations and unknowns. By using multiple shooting with m points, the amplification factors are typically reduced to the m-th root of the starting ones, at the cost of multiplying by m the dimension of the system of non-linear equations to be solved.

Linear behaviour around a periodic orbit of a Hamiltonian autonomous system
An initial condition x 0 of a T -periodic orbit is also fixed point of φ T . In the case of a Hamiltonian autonomous systemẋ = X(x), this fact by itself was not enough in order to find x 0 numerically, but it is useful to study the linear behavior of the flow around x 0 . Let M := Dφ T (x 0 ) be the monodromy matrix of our periodic orbit. Because of the autonomous character of our system and the fact that it has a first integral (the Hamiltonian), M has 1 as a double eigenvalue (for a proof see, e.g., [31]). Moreover, M is a symplectic matrix (see e.g. also [31]), which implies that, if λ is an eigenvalue of M , then 1/λ is also eigenvalue. Now assume that our system is, as the RTBP, of three degrees of freedom, this is, x ∈ R 6 . Then the eigenvalues of M are In the remaining discussion, we will assume that |λ i | ≤ |λ −1 i |. The linear behaviour around a periodic orbit in our 3-degrees-offreedom Hamiltonian system is better studied in terms of Hénon's stability parameters [22], that are defined as A calculation shows that From the discussion of subsection 3.5, if s i ∈ R, |s i | > 2 (hyperbolic case), the stable (resp. unstable) manifold of x 0 as fixed point of φ T is tangent to the λ i (resp. λ −1 i ) eigendirection. This means that the periodic orbit has a stable (resp. unstable) manifold, and its section through the λ i , λ −1 i eigenplane is tangent to the λ i (resp. λ −1 i ) eigendirection. If s i ∈ R, |s i | ≤ 2 (elliptic case), assume λ i = cos ρ + i sin ρ and that v is an eigenvector of eigenvalue λ i . As we have seen, there is a continuous, one-parametric family of closed curves invariant by the linearization of φ T around x 0 in the {x 0 + α 1 Re v + α 2 Im v} α 1 ,α 2 ∈R plane, with rotation number ρ. According to the discussion in subsection 3.5, the full non-linear flow φ T possesses a Cantorian family of invariant curves around x 0 , with limiting rotation number ρ.
When transported by the flow, these invariant curves generate twodimensional invariant tori. Rotation numbers of the form ρ = 2πn/m give rise to bifurcated families by multiplication of the period by m (further details on this kind of bifurcations can be found in [37]). The particular values ρ = 0 and ρ = π, which correspond to s i = 2 and s i = −2, respectively, are known as the parabolic case.
Note that, if a stability parameter s i satisfies |s i | < 2 on a range of energies, since for each energy in this range a one-parametric family of tori is born, across energies this family of tori becomes two-parametric.

Computation of 2D invariant tori
This subsection is devoted to the computation of 2D invariant tori. The method discussed, first introduced in [6], consists in looking for a curve inside the torus invariant by the time-T flow, where T is one of the periods of the torus. The formulation will be made explicit for an autonomous Hamiltonian system with three degrees of freedom (as the RTBP), but it can be modified to account for systems with a different number of degrees of freedom, non-autonomous 6 or not Hamiltonian ones.

Looking for a parameterization of an invariant curve
According to KAM theory (see, e.g., [26]), a 2D torus born around the collinear points of the RTBP can be parameterized by a function ψ(θ 1 , θ 2 ), 2π-periodic in θ 1 , θ 2 , satisfying an invariance equation of the form where (ω 1 , ω 2 ) is the vector of frequencies of the torus. Looking for a torus can be reduced to looking for an invariant curve inside it by observing that ϕ(ξ) := ψ(ξ, 0) parameterizes a curve invariant by φ 2π/ω 2 , and satisfies for ρ = 2πω 1 /ω 2 and ∆ = 2π/ω 2 . Once we have ϕ, we can recover ψ by A calculation shows that, if ϕ is a 2π-periodic function satisfying (20), thenψ as defined by (21) is 2π-periodic in each variable and satisfies the invariance equation (19) for ω 1 := ρ/∆, ω 2 := 2π/∆. In order to turn (20) into a finite system of non-linear equations, we can take ϕ as a truncated Fourier series, (20) at as many values of ξ as the number of Fourier coefficients needed. This is, we will look for ϕ defined as in (22) satisfying The fact that our dynamical system is autonomous gives rise to two indeterminacies: An invariant curve inside a torus is not unique: if ϕ(ξ) satisfies (20) or (23), then φ t ϕ(ξ) also does, for any t ∈ R.
The first indeterminacy can be eliminated by prescribing the value of a coordinate of A 0 (the value chosen must be valid for the torus we are looking for). The second indeterminacy can be eliminated by prescribing a coordinate of A 1 to be zero: if we denote A 1 = (A 1,1 , . . . , A 1,6 ), B 1 = (B 1,1 , . . . , B 1,6 ), and assume that j ∈ {1, . . . , 6} is such that (A 1,j , B 1,j ) = (0, 0), since we can always choose ξ 0 such that A 1,j cos ξ 0 − B 1,j sin ξ 0 = 0. With the two indeterminacies removed in this way, there is a one-to-one correspondence between (approximate) Fourier coefficients of parameterizations of invariant curves ϕ solution of (23) and invariant 2D tori of our dynamical system.

The system of equations
By solving system (23) with its two indeterminacies removed, we could compute an invariant curve ϕ of a torus with "longitudinal period" ∆ and rotation number ρ, that via (21) would correspond to a torus with frequencies ω 1 = ρ/∆, ω 2 = 2π/∆. We could also use this system in order to do continuation with respect to ∆ and/or ρ and, in this way, obtain the corresponding 2-parametric family. Nevertheless, we make two more considerations before stating the final system of equations that we will solve: We want to be able to prescribe values for the energy, so we will add an extra equation for this.
We want to overcome the effects of instability, so we will use multiple shooting, as we did in subsection 3.6.2.
We will, therefore, look for ϕ 0 , . . . , (except for a coordinate of A 0 0 and another one of A 0 1 , according to the previous subsection) with h, ∆, ρ ∈ R, A l j , B l j ∈ R 6 and In order to compute a single torus, we can solve system (24) keeping constant, in addition to the coordinates given by the considerations of the previous subsection, two parameters among h, ρ, T . This will fix a torus within its two-parametric family. In order to continue this torus via the pseudo-arclength method, only one of the parameters h, ρ, T must be keep fixed. Two interesting cases are: To fix ρ to a number with good Diophantine properties. For instance, a noble number (a number number with continued fraction expansion equal to one from a point on).
To fix h, in order to follow an iso-energetic family of tori. In this case, care must be taken because the family is not continuous but Cantorian: the pseudo-arclength method will work as long as the gaps due to resonances are small.
Note that, both in the computation of a single torus and in the continuation of a one-parametric subfamily, we end up with a system of non-linear equations with more equations than unknowns. Namely, in the first case the system is 1 . This is not a problem as long as we use the modified Newton method of subsection 3.2. Note that, when solving the linear system for the Newton correction, d must be set to zero in the case of the computation of a single torus, whereas it must be set to one in the case of continuation of a one-parametric subfamily.
Once a torus has been computed (either individually of by continuation), an estimate of its error can be obtained by evaluating the invariance equation in a refinement of the discretization in ξ used for its computation. In this way, we can use the estimate forξ j = j2π/M and M 1 + 2N f . The value of this estimate can be used to choose the number of Fourier coefficients N f . When doing continuation, N f can be increased or decreased in order to keep this error estimate within a prescribed interval. Observe that large values of N f will give rise to large systems of equations. The time needed for their solution, which requires O((6m(1 + 2N f )) 3 ) operations, will overcome the time needed for numerical integration and become the computational bottleneck of the procedure.

Starting from the central part of a periodic orbit
According to the discussion of subsection 3.7, a family of periodic orbits with an elliptic stability parameter in a range of energies gives rise to a 2-parametric family of invariant tori. Here we will develop formulae from the linear flow around the backbone periodic orbit in order to obtain initial conditions to start the continuation of such a family of tori using system (24).
For an arbitrary function G, let us denote the linearization of G around y 0 as L y 0 G (y) = G(y 0 ) + DG(y 0 )(y − y 0 ). Let x 0 be an initial condition of a T -periodic orbit, with stability parameters s i = λ i + λ −1 i satisfying |s i | ≤ 2 and λ i = cos ν + i sin ν. If we define F := φ T , then x 0 is a fixed point of F , and equation (11) provides an expression for an invariant curve of the linearized flow. In this expression, ξ can be substituted by ξ − ξ 0 , and then we have that which is the linearized-flow version of equation (20). Therefore, as initial seed to get a torus around the o.p., we can take The parameter γ should be chosen small enough for equation (27) to be a good approximation of equation (20). All the computations of subsection 3.9 have been done with either γ = 10 −3 or γ = 10 −4 . The free parameter ξ 0 can be chosen in order to make zero a coordinate of A 0 1 , and in this way avoid the second indeterminacy discussed in subsection 3.8.1. An additional problem when computing a first torus around a periodic orbit is that the periodic orbit has a large basin of attraction and is also a (singular) solution of system (24). A way to prevent falling back to it during the Newton iterations is to keep constant a coordinate of A 0 1 or B 0 1 . A good choice is B 0 1,j , for j such that A 0 1,j is being kept equal to zero in order to prevent the second indeterminacy of subsection 3.8.1.
When we obtain a first invariant curve around a periodic orbit in this way we will say that we are "starting longitudinally to the periodic orbit", because we obtain a tiny invariant curve around x 0 for which, in the evaluation of the flow in (24), numerical integration in order to come back to it is "along the periodic orbit". It will be convenient later to be able to obtain a first invariant curve not tiny but approximately of the same size of the periodic orbit and close to it. We will call this second strategy "starting transversally to the periodic orbit".
In order to develop formulae for this second case, we first globalize the invariant curveφ of the linearized flow to a whole 2D torus bȳ A calculation shows thatψ satisfies the linearized-flow equivalent of the invariant equation (19), namely The invariant curve we are looking for will be close toψ(0, θ 2 ). In order to find it, and since ν can be substituted by ±ν + j2π in all the previous expressions, we can take as initial seed and A l k , B l k coming from a Discrete Fourier transform (DFT) of {ψ(0, j 2π N )} N −1 j=0 . Some notation for the DFT and its relation with Fourier coefficients is developed in subsection 5.2.

Numerical exploration of the dynamics around the L 1 point of the Earth-Moon RTBP
The goal of this subsection is to implement the hierarchical approach mentioned at the beginning of this section in order to systematically compute families of periodic orbits and tori around a collinear libration point of the RTBP. This can be also seen as numerically growing the center manifold of the collinear libration point. The numerical results shown, which are a subset of the ones in [16], will be for the L 1 point and the Earth-Moon mass ratio. In all the computations of this subsection, the flow of the RTBP and its differential with respect to initial conditions have been evaluated according to the discussion of subsection 2.4, using as one-step method with step size control for numerical integration a Runge-Kutta-Fehlberg one of orders 7 and 8 [10] with tolerance 10 −14 . The value used for the Earth-Moon mass ratio is µ = 1.215 0585 6096 2404 · 10 −2 , as obtained from the DE406 JPL ephemeris file [40].

Periodic orbits
The linear behavior around the L 1 point for the Earth-Moon mass ratio is of the type center×center×saddle [42]. This is, if we denote byẋ = X(x) the vector field of the RTBP, as in equations (1), (2), we have for ω p , ω v , λ > 0. As discussed in subsection 3.4, Lyapunov's center theorem ensures that each center gives rise to a family of periodic orbits. In the expression for Spec DX(L 1 ) above, ω p (resp. ω v ) can be chosen in such a way that the eigenplane corresponding to the eigenvalues ±iω p (resp. ±iω v ) is contained in {z = p z = 0} (resp. {x = p x = y = p y = 0}). Because of this, the family of periodic orbits related to the ±iω p (resp. ±iω v ) eigenvalues is known as the planar (resp. vertical) Lyapunov family. Initial guesses to start the numerical continuation (see subsection 3.1) of these families can be obtained from (9). When doing Newton iterations on system (17) to find the first periodic orbit, a convenient way to avoid falling back to the L 1 point (which is a singular solution of system (17) with a large basin of attraction) is to keep constant a coordinate of x 0 . A convenient way to represent a family of periodic orbits that has been obtained by numerical continuation is by plotting the period and the stability parameters (18) of its orbits with respect to energy. The period vs. energy curve is known as characteristic curve. Figure 4 represents the characteristic curve and stability parameters of the vertical Lyapunov family. This family starts at energy −1.59417 (the one of L 1 ), has a bifurcation at energy −1.49590, that will be commented later, and ends at a large planar orbit with energy 0.41391, that surrounds the Earth and the collinear points L 1 , L 3 . Plots of sample orbits of this family and all the other families of periodic orbits that we will consider can be found in [32].
In figure 5 we have represented the characteristic curve and stability parameters of the planar Lyapunov family. This family starts at energy −1.59417 (the one of L 1 ), has several bifurcations and ends at a collision with the Earth. 7 According to [23], the only possible kinds of bifurcation from the planar Lyapunov family to a family of three-dimensional orbits are the ones sketched in figure 6. Types A and B correspond to a stability parameter crossing 2, whereas types C and D correspond to a stability parameter crossing −2 (and thus are period-doubling bifurcations). In cases A, B not one but two families of periodic orbits bifurcate from the planar family. The two bifurcated families are symmetric with respect to {z = 0}. Assuming that the Poincaré section used in the continuation of the planar Lyapunov family is {y = 0} 8 , an initial condition for one of such bifurcated orbits can be obtained by doing a small displacement in the z coordinate for types A, C, D, and in the p z coordinate for type B. The displaced coordinate can be kept constant during Newton iterations on system (17) in order to avoid falling back to an orbit in the planar Lyapunov family. The bifurcations found for the planar Lyapunov family, together with its classification according to [23], are given in  large range of energies halo orbits have complex (non-real) stability parameters; figure 8 zooms figure 7 in order to show the transition from real to complex stability parameters and vice-versa. In figure  8 left, it is also shown how the small stability parameter goes across 2 cos(2π/3) once, at energy −1.52944, and across −2 = 2 cos(2π/2) twice, at energies −1.51081, −1.51033. The first case gives rise to two period-triplicated bifurcated families of periodic orbits, one with elliptic-hyperbolic normal behaviour and the other one with ellipticelliptic normal behaviour. The second case gives rise to a periodduplicated family of periodic orbits with elliptic-elliptic normal behaviour. The third case gives rise to another period-duplicated family of periodic orbits but with elliptic-hyperbolic normal behaviour. These three bifurcations take place for each of two symmetric halo families. As discussed in subsection 3.7, there are many more bifurcations, but these three will play a role in the computations of invariant tori of the next subsection. The actual initial conditions used to find orbits of these families have been found by shooting from invariant tori nearby. The second bifurcation of the planar Lyapunov family gives rise to two families, symmetric with respect to z = 0, that can be thought as a two-lane bridge that connects the planar Lyapunov family with the vertical one at its bifurcation at energy −1.49590. Some orbits of this family are shown in figure 9. Table 1 still reflects a third bifurcation of the planar family that we do not follow, because it takes place at an energy outside the range of energies that will be reached by the continuation of invariant tori of the next subsection.  figure 7 showing the transitions to and from complex saddle.

Invariant tori
The first families of tori around the libration point L 1 that we will compute will be the ones of constant rotation number ρ starting longitudinally from the vertical Lyapunov family of periodic orbits. The range of values of ρ to be considered is thus provided by the values of ν > 0 such that, according to subsection 3.8.3, 2 cos ν is one of the stability parameters of the base vertical Lyapunov orbit. Therefore, initiating the continuation of each constant ρ family requires to find a initial condition of a vertical periodic orbit corresponding to a specific value of ν. This initial condition is obtained by continuation of system (17). Since ν is not a continuation variable, it must be considered a function of a continuation variable, for instance ν = ν(h). The value of h providing a prescribed value of ν(h) can be found by a numerical one-dimensional zero-finding method. A good choice is Brent's (see e.g. [36]), since it has fast, global convergence and does not require computing derivatives. If we represent the value of ν with respect to energy along the vertical Lyapunov family of periodic orbits for the range of energies in which they have central part (see figure 4), we obtain the curve labeled β in figure 10. This curve goes from the point P 2 , that corresponds to the collinear point L 1 , to the point P 3 , that corresponds the bifurcation of the vertical Lyapunov family at energy −1.4959 (see figure 4). Our first continuation of families of tori, with constant rotation number, has been done by choosing an approximately equally spaced grid of noble values of ρ (in order to stay away from resonances, as discussed in subsection 3.8.2), ranging from the ordinate value of the point P 2 of figure 10 to the maximum value of ν along the β curve, and starting longitudinally from the leftmost planar Lyapunov periodic orbit of the β curve with stability parameter 2 cos ρ. Each of the obtained families of tori, with constant ρ, that would be seen in figure 10 as a horizontal line, collapses at a vertical Lyapunov orbit of higher energy, as the shape of the β curve suggests. With this first continuation we cover the region in figure 10 delimited by the curves α, β, γ with ρ ≥ ρ(P 2 ). This continuation of families of invariant tori, and all the remaining continuations that we will describe, have been done by solving system (24) with m = 2, a tolerance of 10 −11 for Newton iterations, and with continuation step size control with n des = 4 in (7). In the continuation of each family, the number of harmonics N f of the Fourier expansions (25) has been determined in order to keep the estimate (26) under 10 −10 . In addition to this, an upper limit of N f = 100 has been set. When this limit is reached, the error estimate (26) is allowed to grow up to 10 −8 and, when this happens, the continuation is stopped. This has never happened in this first exploration. In order to cover the region within the curves α, β, γ with ρ < ρ(P 2 ), a possibility would be to start from the β curve and go downwards. This means to perform continuation of families of tori with h constant. If h is close to the energy of L 1 , the iso-energetic family of tori obtained should end by collapsing to a planar Lyapunov orbit, because this is what happens linearly. The actual tori of such a continuation are shown in figure 11, for h = −1.59. Although the tori do collapse to a planar orbit, the corresponding invariant curves ϕ 0 obtained by solving system (24) do not collapse to a point but tend to the whole ending planar Lyapunov orbit. The limiting value of ρ is numerically checked to be where ν is such that 2 cos ν is a stability parameter of the ending planar Lyapunov periodic orbits. Therefore, according to (28), the same invariant curves within the tori of figure 11 could be obtained by starting transversally from this ending planar Lyapunov orbit. The α curve of figure 10 is obtained by plotting expression (31) as a function of h, with ν such that 2 cos ν is a stability parameter of the planar Lyapunov orbit of energy h. The point with label P 1 in this curve corresponds to the bifurcation of the planar Lyapunov family of periodic orbits to the halo families (see figure 5 and table 1). The family of tori of figure 11 would be seen in figure 10 as a vertical line with h = −1.59 that goes from the curve β to the curve α. In order to complete the computation of invariant tori within the curves α, β, γ, and in order to avoid "jumping over resonances", we go back to the constant ρ continuation strategy. From the discussion in the last paragraph, the remaining tori within the α, β, γ curves can be computed by starting transversally from the family of planar Lyapunov periodic orbits, for an approximately equally spaced grid of noble values of ρ of the form (31), for the range of values of ν that produced the α curve. When doing so, some of the corresponding constant-ρ families of invariant tori with largest ρ value have reached a vertical Lyapunov periodic orbit of higher energy. The remaining ones have stopped due to the N f = 100 computational limit. For each value of ρ in which this has happened, we have also continued for decreasing energies the family with constant ρ starting longitudinally from the rightmost vertical Lyapunov periodic orbit of the β curve with this ρ value. In this way, we have covered with invariant tori all the region within the α, β, γ curves of figure 10 except for the one labeled as > 100. By allowing for N f > 100, some of the tori of this last region could be computed. Many of them, however, simply do not exist, because, as we will see later, as ρ goes to zero for fixed energies larger than the one of the point P 1 in figure 10, we approach homoclinic connections of periodic orbits. The α, β, γ curves of figure 10 delimit a set of tori that can be considered a single family, since all of these tori can be reached by numerical continuation starting from L 1 . Close to L 1 , the tori of this family are the ones given by KAM theory. Trajectories in them are known as Lissajous trajectories by the astrodynamics community. We will thus denote this family as the Lissajous family of invariant tori. Tori in this family can be considered to have "natural" frequencies ω v (T, ρ), ω p (T, ρ), obtained by continuation from the frequencies ω v , ω p of the collinear point L 1 in equation (30). An application of Lyapunov's center theorem shows that  Invariant tori around the elliptic-hyperbolic period-duplicated halo-type family of periodic orbits, in an energy range analogous to the previous one.
Invariant tori around planar Lyapunov orbits, in a short energy range starting at the bifurcation of the two-lane bridge joining it with the vertical one, in order to complete the Poincaré sections of figure 14.
Except for the last one, these families are represented in figure 12 in h-ρ plots analogous to figure 10. Contrary to the Lissajous family of invariant tori, none of these new families has been described completely. The numerical continuations have been stopped when the N f = 100 computational limit has been reached. How these families further evolve is an open question.

Iso-energetic Poincaré sections
Since the center manifold of L 1 , W c (L 1 ), is four-dimensional, its restriction to an energy value, W c (L 1 ) ∩ {H = h}, would be threedimensional, and a Poincaré section in this restriction, W c (L 1 )∩{H = h} ∩ Σ, would be 2-dimensional. Following [13,25], it is convenient to visualize W c (L 1 ) by a sequence of iso-energetic Poincaré sections. This is done in figures 13 and 14, using the Poincaré section Σ := {z = 0, p z > 0}. In order to be able to produce these figures, in the continuation of each constant ρ family of tori of the previous subsection, the tori of the energies of the plots of figure 13 have been obtained by doing Newton iterations keeping h constant, starting from pseudo-arclength predictions from nearby tori (see algorithm 3). All the plots of figure 13 have a similar structure. The exterior curve in each plot is a Lyapunov planar orbit of the energy level cor-responding to the plot. Strictly speaking, the Poincaré section is not valid for this orbit, so it should not have been plotted. Nevertheless, it is useful to use it as boundary of W c (L 1 ) ∩ Σ at the energy of the plot. The closed curves inside the region bounded by the Lyapunov planar orbit are the intersections with Σ of the invariant tori computed in the previous subsection. These intersections are computed through algorithm 1, starting from the invariant curve ϕ 0 (see system (24)) of each torus.
In all the plots there is a fixed point on the x axis associated to the vertical Lyapunov orbit. This point is not represented, but outlined by the smallest blue curves. For small energy values, the whole picture is formed by invariant curves surrounding this fixed point. They are associated to the intersections with Σ of Lissajous-type trajectories around the vertical periodic orbit, whose evolution from the vertical Lyapunov periodic orbits to the planar one is similar to the one displayed in figure 11. At the energy levels in which halo orbits have bifurcated from the planar Lyapunov family, there appear two additional fixed points, again not represented but outlined by the smallest violet invariant curves. Increasing the values of the energy, the halo family undergoes the two bifurcations mentioned in subsection 3.9.1, by period triplication and duplication. Within the bifurcated families there are some with central part, which are surrounded by invariant tori, also computed in the previous subsection, whose Poincaré sections provide here the red invariant curves. These invariant curves give rise to the "island chain" structure typical of two-dimensional area-preserving maps (compare with figure 3). To display more clearly this behaviour, the last plot of figure 13 displays a magnification of the bifurcated periodic orbits and its surrounding invariant tori.
The region between the tori around the vertical Lyapunov orbit and the tori around the halo orbits is not empty, as it appears in the above figures. It should contain, at least, the traces on the surface of section of the invariant manifolds of the Lyapunov planar orbit. These manifolds act as separatrices between both kinds of tori. The same thing happens between the islands of the bifurcated halo-type orbits and the tori around halo orbits. In this case, the region between both kinds of tori is filled with the traces of the invariant manifolds of the bifurcated hyperbolic halo-type orbits. In all these boundary regions, the motion should have chaotic behaviour. The numerical methods of this section are not able to capture this chaotic motion, but the semi-analytical methods of the next section can capture it.
The plot corresponding to energy −1.507, shown in figure 14, has more structure. For this energy level, the two-lane bridge between the planar and vertical Lyapunov families of periodic orbits has al-ready bifurcated, so the planar family has gained central part, and its periodic orbits are again surrounded by invariant tori. The {z = 0} sections of these tori are the outermost curves that appear in figure  14 (in this case, the planar Lyapunov periodic orbit, that surrounds all these curves, is not represented). In the figure, the two diamond points are the fixed points corresponding to the intersections of the two orbits of the bridge with the surface of section. The invariant manifolds of these bifurcated periodic orbits are the ones that must act as separatrices between the tori around the halo orbits and the tori around the vertical Lyapunov orbit of this energy.

Semi-analytical computation of invariant objects using the parameterization method
The parameterization method is an approach to the study of invariant manifolds, whose general idea is to seek for parameterizations of invariant manifolds as solution of invariance equations, that are simplified through changes of variables that exploit geometrical properties. It is a strong point of this approach that "theoretical" and "numerical" are two aspects of the same philosophy. On the one hand, the proofs are constructive and can be turned into algorithms. On the other hand, these algorithms, when implemented with rigorous numerics based on interval arithmetic, can be turned into computer assisted proofs. Since its introduction in [5], it has been used by many authors. A recent review, that also has some original developments, can be found in [21]. Here we will be concerned with the use of the parameterization method for the (non-rigorous) computation of Taylor expansions of invariant manifolds around fixed points of flows. It will be applied to the computation of the center manifold of the collinear points L 1 , L 2 , of the Earth-Moon RTBP. In this way, this variant of the parameterization method can be seen as a semi-analytical technique for the computation of the invariant objects inside the center manifold of the collinear points of the RTBP. An earlier technique, known as reduction to the center manifold ( [13,25]), produces essentially the same results. The parameterization method has some advantages in computational speed, generality (the implementation is independent of the dynamical system under study, the RTBP in our case) and flexibility, since the coordinates of the manifold can be adapted to the dynamics, as we will see in subsection 4.4.2.
The discussion will follow chapter 2 of [21] except for some notational changes, additional computations and plots. The software package in http://www.maia.ub.edu/dsg/param/ includes a C routine that computes expansions of invariant manifolds of fixed points of flows as described below.

The method
Assume we are given a continuous, n-dimensional dynamical systeṁ x = X(x) with a fixed point p at which the differential of the vector field is diagonalizable. We would like to compute a d-dimensional manifold that contains the fixed point and is tangent to a d-dimensional eigenspace of the differential of the vector field. By a change of variables of the form x = p + P y, our original system can be turned intoẏ = Y (y), y = (y 1 , . . . , y n ), with DY (0) = diag(λ 1 , . . . , λ n ), λ i ∈ C, in such a way that the eigenspace of interest is {y ∈ R n : y d+1 = · · · = y n = 0}. Then our goal is to compute an expansion of a d-dimensional manifold that contains the origin, is invariant by the flow, and is tangent to the y 1 , . . . , y d coordinates.
To do this, we look for W : C d −→ C n , parameterization of the manifold, and for f : C d −→ C d , the vector field reduced to the manifold. In this way, if we denote by s ∈ C d the parameters describing the manifold, then the differential equations in parameter space are ṡ = f (s). From the parameterization of the manifold W (s) in the y variables, a parameterization of the manifold in the original x variables can be recovered asW In order to find W , f we need to solve the invariance equation: Assume that W , f are expanded as power series in s, with W k n-vector and f k d-vector of homogeneous polynomials of degree k in s = (s 1 , . . . , s d ),

Now assume that
f <k := f 1 + · · · + f k−1 are known. If we restrict equation (33) to its terms of order k, we obtain the order-k cohomological equation. By putting all the unknown terms in the left-hand side and all the known terms in the right-hand one, we obtain as right-hand side where [ ] k stands for "terms of order k". The evaluation of the second term in the previous expression involves products of homogeneous polynomials. The first term, which consists in plugging the known part of W into the vector field and obtaining the terms of degree k, is computationally more costly. High efficiency is achieved through the use of automatic differentiation, as will be discussed below.
The expression for the left-hand side of the order-k cohomological equation depends on the component. The whole order-k cohomological equation reads whereλ := (λ 1 , . . . , λ d ), λ , m := λ 1 m 1 + . . . λ d m d . The manifold can be computed as long as (36) can be solved, this is, there are no m ∈ N d , i ∈ {d + 1, . . . , n} such that λ i = λ , m , which would be a cross-resonance. The solution of (35) can be done in several ways, that give rise to different styles of parameterization: The graph style, that consists in taking W i k,m = 0, f i k,m = R i k,m , as to obtain W 1 (s) = s 1 ,. . . , W d (s) = s d , so that, in y coordinates, the manifold is the graph of the function (W d+1 , . . . , W d ).
The normal form style, in which the expansion of f is taken as simple as possible: When λ i = λ , m for i ∈ {1, . . . , d}, one speaks of an internal resonance.
The following mixed style, that, given sets of indexes I 1 , . . . , I N ⊂ {1, . . . , n}, turns the sets {s i = 0, i ∈ I l }, l = 1, . . . , N , into invariant submanifolds: if ∃l: i ∈ I l and m j = 0 ∀j ∈ I l , This mixed style allows adapting the parameterization to the dynamics, as will be shown in the examples that follow.
Note that, as a whole, the order-k cohomological equation is linear and diagonal: each unknown monomial of the left-hand side is computed as a constant times the corresponding monomial of the right-hand side. All the computational effort goes in the evaluation of R k .

Efficiency considerations
Once the R k term is computed, the solution of the order-k cohomological equation with any of the styles previously mentioned is very fast. Assuming that we have explicit formulae for the vector field, as is the case in the RTBP, the evaluation of R k as given in (34) depends on both being able to perform sums and products of dense 10 multivariate polynomials and being able to compose truncated (multivariate) power series into elementary functions such as sine or square root. An strategy for an efficient implementation of the product of homogeneous polynomials is to represent them recursively with respect to the number of variables. A d-variate homogeneous polynomial of degree k can be represented as a linear combination of (d − 1)variate polynomials of degrees k, k − 1, . . . , 0: for s = (s 1 , . . . , s d ), s = (s 1 , . . . , s d−1 ), The memory representation can be made to mimic this recursive definition. The use of this strategy avoids the need for hash tables and reduces the product of homogeneous polynomials to dot products of vectors of coefficients.
With respect to the composition of truncated Taylor expansions into elementary functions, an efficient strategy is the use of a form of automatic differentiation 11 based on the notion of radial derivative. The radial derivative of a function f : R n → R is defined as On an homogeneous polynomial of degree k, it satisfies It also satisfies a form of chain rule: for a function ϕ : Now, if ϕ satisfies a differential equation, the previous two properties can be used to deduce a recurrence that relates the series expansions of f and ϕ • f . For instance, for Using this recurrence, p k can be computed from f 1 , . . . , f k−1 and p 0 , . . . , p k−1 . This is, the terms of order < k of ϕ • f are also needed.
Because of this, in order to proceed order by order in the computation W , f , we need to store the power series expansions of all the intermediate operations in the evaluation of [F (W <k (s))] k that involve the composition of a power series with an elementary function. The software package in http://www.maia.ub.edu/dsg/param/ includes a C library for the manipulation of multivariate, truncated power series that implements all these ideas.

Error estimation
Once we have computed

Using the graph style
The first example will be the computation of W c (L 1 ) using the graph style. Denote the vector field of the RTBP in Hamiltonian form aṡ x = X(x), and denote the eigenvalues of DX(L 1 ) as in equation (30). Denote as P a matrix having as columns eigenvectors of eigenvalues iω p , −iω p , iω v , −iω v , λ, −λ, in this order. For this example, we apply the procedure of subsection 4.1 with n = 6, d = 4 to Y (y) := P −1 X(L 1 + P y) , using the graph style. In this way, we obtain a parameterization of the 4D center manifold of L 1 as with W i (s) = s i , i = 1, 2, 3, 4. Expansions of W have been computed for several orders.  Table 2: For several orders, computing times (in seconds) of the expansions of W c (L 1 ) for the Earth-Moon RTBP using the graph style, on a Mac with Intel Core Duo @ 2.16GHz. Figure 15 shows the {s 4 = 0} Poincaré sections of several trajectories at fixed energies. Note that each point in these plots uniquely determines a trajectory: s 3 is computed from s 1 , s 2 and the (fixed) value of the energy. The Poincaré sections in figure 15 are analogous to the ones computed in [13,25]. Since through the parameterization (37) points with s 4 = 0 go to points with z = 0, the Poincaré sections in figure 15 are also analogous to the ones in figure 13. Note that they are obtained in completely different ways: here by direct numerical integration ofṡ = f (s); there by computing individually every torus represented. Figure 13 can reach higher energies because of the numerical approach. Here, the use of the expansions is limited to their domain of validity. An estimation of this domain is shown next. Here, on the other hand, the numerical integration ofṡ = f (s) allows us to capture all the dynamics in the center manifold at each energy level. In the numerical approach of figure 13, we can only display the objects that we individually compute.   10, 20, 30. Compared to the expansions around L 1 , the domain of validity is smaller, but the precision is about the same for the same energies. This is coherent with the fact that the energy of L 2 is larger than the one of L 1 .  an error that is quadratic in the distance to the base object, which is adequate for many applications, including preliminary mission design. Approximations of higher order can be obtained through semianalytical techniques, including the parameterization method, as will be discussed in section 6. The Lindstedt-Poincaré method [29] is another semi-analytical alternative.

Invariant manifolds of periodic orbits
Let x 0 be an initial condition of a periodic orbit of period T , this is, φ T (x 0 ) = x 0 . A parameterization of the periodic orbit as an invariant manifold is given by the 2π-periodic function ϕ : [0, 2π] −→ R 6 defined as ϕ(θ) := φ θ 2π T (x 0 ). Let Λ ∈ R, |Λ| = 1 be an eigenvalue of the monodromy matrix of the periodic orbit with eigenvector v, this is An eigenvalue Λ with |Λ| > 1 (resp. |Λ| < 1) would correspond to an unstable (resp. stable) manifold. For brevity, let us assume for the rest of the discussion that Λ > 0; a comment will be made on the case Λ < 0 at the end. Therefore, Λ > 1 (resp. Λ < 1) would correspond to an unstable (resp. stable) manifold. A parameterization of a set of vectors tangent to the unstable (resp. stable) manifold, also know as unstable bundle (resp. stable bundle), is given by the 2π-periodic function v(θ) := Λ − θ 2π Dφ θ 2π T (x 0 )v. By combining the two previous expressions, we can obtain a parameterization of the linear approximation of the unstable (resp. stable) manifold:ψ (θ, ξ) := ϕ(θ) + ξv(θ).
It satisfies the approximate invariance equation for ω = 2π T , λ = ω ln Λ 2π . It can thus be evaluated for small ξ only. Nevertheless,ψ can be used to globalize the manifold by numerical integration while still providing a cylinder-like parameterization: for ξ not necessarily small, we can take an integer m > 0 (resp. m < 0) such that Λ −m ξ is small and computē  In the case Λ < 0, all the previous discussion is valid if we substitute T by 2T and Λ by Λ 2 . In this way, v is 2π-periodic and the expressions forψ,Ψ still provide cylinder-like parameterizations.

Invariant manifolds of 2D tori
We follow the discussion of [24] with a slightly modified computational strategy. Assume that ϕ parameterizes an invariant curve inside a 2D torus, as in subsection 3.8.1, We want to find Λ ∈ R, |Λ| = 1 and u : R → R 6 , 2π-periodic, s.t.
this is, an invariant bundle associated to the eigenvalue Λ. It will be an unstable (resp. stable) invariant bundle if |Λ| > 1 (resp. |Λ| < 1), that will be tangent to the unstable (resp. stable) manifold of the torus on the invariant curve parameterized by ϕ.
Equation (41) can be compactly written as Assuming that u is expanded as a truncated Fourier series, the eigenvalue problem (42) can be discretized and thus converted in a finitedimensional matrix-vector eigenvalue problem by approximating the Fourier coefficients of Cu by their Discrete Fourier Transform (DFT). We use the following notation for the DFT: for N even, given real data {f j } N −1 j=0 , we denote with δ 0 = δ N 2 = 1, δ k = 2 for k = 1, . . . , N 2 − 1. If the data comes from the regular sampling of a 2π-periodic function, this is, f j = f (θ j ) for θ j = j2π/N and f is 2π-periodic, where the approximation is an equality if θ = θ j , 0 ≤ j ≤ N − 1. In this way, the DFT coefficients provide an approximation of the Fourier coefficients (for a bound on the difference, see e.g. [11,17]).
If we denote then, for a suitable (finite-dimensional) matrix C, The columns of C can be found as the DFT coefficients of the operator C applied to the canonical basis elements in x space, this is, to the functions w k , w k cos(θ), w k sin(θ), w k cos(2θ), w k sin(2θ), etc., being w k ∈ R 6 the k-th element of the canonical basis, k = 1, . . . , 6. Since all these functions can be written in terms of complex exponentials of the form e ikθ , the coefficients of the C matrix can be computed from which, after a few calculations, is found to be for j = 1, . . . , 6 and k, m = 0, . . . , N/2. Since Dφ ∆ (ϕ(θ l − ρ))w j is a 6-vector for each j, the computation of all the needed values of the previous expression is reduced to 36 DFT, which, by using FFT, are computed in O(N log N ) operations each. Some knowledge on the structure of the spectrum of the invariant bundle we are looking for is necessary in order to choose the right eigenvalues of the C matrix of (43). The eigenvalues of C appear grouped in circles. Since the tori we are looking for are reducible, there are as many circles as eigenvalues of the reduced matrix (which can be considered analogous to the monodromy matrix of a periodic orbit). Assuming that (41) has a solution, from the fact that the RTBP is a Hamiltonian system, apart from unit circles there will be a circle containing Λ and another circle containing Λ −1 . These are the ones we are interested in. The corresponding eigenvectors provide the Fourier coefficients of the invariant bundles we are looking for. More details on this discussion and some additional considerations on the accuracy of the computed eigenvalues can be found in [24]. Now, from an invariant stable or unstable bundle u(θ), tangent to the stable or unstable manifold of the torus on the invariant curve ϕ(θ), we can obtain the invariant bundle tangent to the stable or unstable manifold of the torus on the whole torus through This expression assumes Λ > 0. If this is not the case, ∆ needs to be changed to 2∆, so equations (40) and (41) are satisfied with ρ substituted by 2ρ and Λ by Λ 2 . Defined as above, the v function is 2π-periodic in θ 1 , θ 2 and satisfies where ψ is the parameterization of the 2D torus defined in equation (21), ω 1 = ρ/∆ and ω 2 = 2π/∆. From the parameterization of the stable or unstable bundle defined on the whole torus, we can write a parameterization of the linear approximation of the stable or unstable manifold of the torus as ψ(θ 1 , θ 2 , ξ) = ψ(θ 1 , θ 2 ) + ξv(θ 1 , θ 2 ), which is 2π-periodic in θ 1 , θ 2 and satisfies the approximate invariance equation φ t ψ (θ 1 , θ 2 , ξ) =ψ(θ 1 + tω 1 , θ 2 + tω 2 , e tλ ξ) + O(ξ 2 ), for ω 1 = ρ/∆, ω 2 = 2π/∆, λ = ω 2 ln Λ/(2π), and thus equation (44) can be evaluated for small ξ only. For ξ not necessarily small, we can consider an integer m such that Λ −m ξ is small (m > 0 for the unstable manifold, m < 0 for the stable manifold) and computē 6 Semi-analytical computation of stable and unstable manifolds using the parameterization method We have seen how the parameterization method can be used as a semi-analytical technique in order to find the periodic orbits and tori in the center manifold of a collinear libration point. Without any modification, it can also be used to find the invariant stable and unstable manifolds of these trajectories. All the unstable manifolds of the invariant objects of W c (L 1 ) are contained in the center-unstable manifold of L 1 , W cu (L 1 ), which is an invariant manifold tangent to the directions given by the eigenvectors with eigenvalues where we have recovered the notation of equation (30). All the stable manifolds of the invariant objects of W c (L 1 ) are contained in the center-stable manifold of L 1 , W cs (L 1 ), which is the invariant manifold tangent to the directions given by the eigenvectors with eigenvalues The parameterization method does not need any modification to compute W cu (L 1 ) or W cs (L 1 ) instead of W c (L 1 ).
As an example, we can apply the procedure described in section 4.1 with the same choice of P as in subsection 4.4, n = d = 6 and choosing a mixed style parameterization with the sets of indexes I 1 , . . . , I 14 defined by table 3. In this way, we obtain a reparameterization of a whole neighborhood of L 1 that is completely adapted to the dynamics. Table 3 is the recipe to choose initial conditions on the different kind of objects. For example, points of the formW (0, 0, 0, 0, s 5 , 0) are in the unstable manifold of L 1 because of the use of I 1 , whereas points of the formW (0, 0, s 3 , s 4 , 0, s 6 ) are in the stable manifold of the vertical Lyapunov family of periodic orbits because of the use of I 10 . Table 4 shows the computing times of the expansions for several orders. These times are now larger than the ones of section 4.4 because the truncated power series have 6 variables instead of 4.
As before, it is necessary to determine a neighborhood of validity of the expansions. This has been done in figure 21, by an exploration similar to the one done in section 4.4.2, but now taking initial conditions with s 5 , s 6 = 0 in evaluation of the e O estimate, and also integrating both forward and backward in time, in order to test both the stable and the unstable manifold. The maximum e O of the trajectories of each energy tested are represented by a point in figure 21. The pairs of green-violet curves correspond, from left to right, to orders 10, 15, 20, 25, 30. The full details of the exploration can be found in [21].
A sample application of the use of these expansions is the generation of what are known as transit and non-transit trajectories [7,8]. With the choice of the eigenvectors corresponding to ±λ shown schematically in figure 22, orbits with s 5 s 6 > 0 are transit in the sense that go from the Earth to the Moon or vice-versa. Orbits with l I l submanifold described by s i = 0, i ∈ I l  1 {1, 2, 3 Table 4: For several orders, computing times of the expansions of the mixed style reparameterization of the neighborhood of L 1 of the Earth-Moon RTBP.  s 5 s 6 < 0, however, are non-transit in the sense that after departing from a primary they "bounce back" to it. Figure 23 shows some trajectories used in the evaluation of the error estimate of figure 21, which are all transit because they were chosen with s 5 = s 6 > 0. For clarity, the trajectories are not integrated as in the evaluation of the error estimate, but forward in time up to the first cut with x = µ − 1 + R M , where R M is radius of the Moon in dimensionless units (red trajectories), and backward in time up to the second cut with y = 0 after the first passage behind the Earth (blue trajectories). Looking at each blue curve followed by the red one as a single trajectory, the plots show that all of them are Earth-Moon transit.

Computation of homoclinic and heteroclinic connections
An homoclinic connection of a object (with itself) is a trajectory that tends to the object both forward and backward in time. An heteroclinic connection of a departing object and an arrival object is a trajectory that tends to the departing object backward in time and to the arrival object forward in time. From the dynamical systems point of view, these connections play a fundamental role in the global organization of the dynamics. In the RTBP, they also provide low-energy transfers between objects [15] and resonance transitions [28,14]. Using Conley-McGehee tubes [7,30] inside Hill's regions, they allow to prescribe itineraries between the interior and exterior regions of a moon-planet system, as in [28,14,34].

Computing individual connections
Consider ψ u (θ, ξ) a parameterization of an approximation of the unstable manifold of a departure object, and ψ s (θ, ξ) a parameterization of an approximation the stable manifold of an arrival object. These approximations can be the linear ones, or of higher order. Let Σ = {g(x) = 0} be a Poincaré section intersected by the manifolds, and consider two associated Poincaré maps: P + Σ , computed integrating forward in time, and P − Σ , integrating backward in time. This is, where the functions τ + (x), τ − (x) are time-return maps with τ + (x) > 0, τ − (x) < 0 defined implicitly by the conditions g φ τ + (x) (x) = g φ τ − (x) (x) = 0.
The intersections of homoclinic (if the departure and arrival objects are the same) or heteroclinic (in the case of different departure and arrival objects) connections with the section Σ are given by the zeros of the function F (θ u , θ s ) = P + Σ (ψ u (θ u , ξ)) − P − Σ (ψ s (θ s , ξ)).
In this function, ξ is a fixed parameter, that needs to be taken small if ψ u , ψ s are linear approximations, or not necessarily, if ψ u , ψ s are approximations of higher order. The θ u , θ s parameters are vectors of phases of the same dimension of the connecting objects (scalars for periodic orbits, 2-vectors for 2D tori).
In the case of periodic orbits, their stable and unstable manifolds are locally diffeomorphic to 2D cylinders. As long as this remains true when globalizing their manifolds, the computation of connections is reduced to intersect 2D tubes, which can be visualized without much difficulty. Their visualization is particularly simple if the orbit is planar: the planar RTBP has 2 degrees of freedom, and therefore a Poincaré section of fixed energy is 2D. Figure 24 shows the manifold tubes of a planar Lyapunov orbit around L 1 of the Earth-Moon RTBP, and also their intersection with Σ := {x = µ − 1}. The two points of intersection of the two curves coming from the sections of the manifold tubes with Σ (figure 24 right) give rise to two homoclinic connections. Initial conditions in order to find zeros of the function F of equation (46) via Newton iterations can be obtained from this plot. Care must be taken with the number of cuts of the manifold that define the Poincaré maps: according to figure 24 left, P + Σ is defined as the second cut with Σ, whereas P − Σ is defined as the first cut. x, y in the left plot, p x , p y in the right one.
In the cases in which the sections of the manifold tubes with Σ are not easy to visualize, other approaches need to be followed. As an example, consider searching for heteroclinic connections between: a Lissajous torus around L 1 of the Earth-Moon RTBP with energyh := −1.58 and rotation numberρ := 0.2800082, and a Lissajous torus around L 2 with the same energy and rotation numberρ := 0.1700025.
Denote as Ψ u (θ 1 , θ 2 , ξ) (resp. Ψ s (θ 1 , θ 2 , ξ)) a parameterization of the linear approximation of the unstable (resp. stable) manifold of the departing (resp. arrival) torus. Denote also as P + Σ , P − Σ the Poincaré sections defined in equation (45) after the needed number of cuts with the section. Then, in order to look for connections, we can plot in terms of θ 1 , θ 2 the function

Continuation of connections
Since the RTBP is a Hamiltonian system, periodic orbits and tori are not isolated but part of families. As a consequence, the connections between them are part of families too. If we want to compute several connections along a family, it is a tedious procedure to compute them individually as described before. The process of computing homoclinic or heteroclinic connections along families can be automated by the use of continuation on Equation (46), by letting ψ u (θ u , ξ) and ψ s (θ s , ξ) evolve freely along the families of departing and arrival objects. The actual way to do it depends on the way that ψ u , ψ s have been obtained, that can be semi-analytical or numerical. In the following we will focus in the numerical approach.
Assume we wanted to numerically compute a family of homoclinic connections of periodic orbits of the RTBP by continuation. Let Σ 1 = {g 1 (x) = 0} be a Poincaré section for the initial conditions of the periodic orbit, and Σ 2 = {g 2 (x) = 0} a Poincaré section used to match the invariant manifolds of the periodic orbit. Assume these Poincaré sections are valid along the portion of the family we want to compute. We need to consider as unknown everything necessary to determine a periodic orbit of the family and its homoclinic connection: the value of the energy, h, the initial condition of the periodic orbit, x 0 , the eigenvalue of the monodromy matrix related to the unstable (resp. stable) manifold, v u (resp. v s ), the departing (resp. arriving) phase on the linear approximation of the unstable (resp. stable) manifold, θ u (resp. θ s ), and, finally, the time of flight from the linear approximation of the unstable (resp. stable) manifold to the surface of section in which the manifolds are intersected, T u (resp. T s ). The system of equations needs to impose all the conditions for h, x, T , Λ u , v u , Λ s , v s , θ u , T u , θ s , T s to determine a periodic orbit and an homoclinic connection of it. It would thus be with, according to (39), Dφ θ 2π T (x)v j , for j = u, s. Note that the system (47) includes a normalization condition on v u , v s , in order to make them to be locally unique. Also observe that, since we use the linear approximation of the manifolds, ξ is a parameter that must be kept fixed at a small value, e.g. 10 −6 . An actual implementation requires multiple shooting, both in the periodic orbit and in the connection. Additional details can be found in [4]. Figure 27 displays a homoclinic connection (in violet) of a large planar Lyapunov orbit (in blue) around L 1 of the Earth-Moon RTBP that has been reached by such a continuation procedure. In order to aid visualization, all the perigees, apogees, periselenes, and aposelenes have been numbered as their appear along the connection.
Dφ ∆ u (ϕ u (θ))v u (θ) − Λ u v u (θ + ρ u ) = 0, Dφ ∆ s (ϕ s (θ))v s (θ) − Λ s v s (θ + ρ s ) = 0, g φ ∆ u * (Ψ u (θ u 1 , θ u 2 )) = 0, g φ ∆ s * (Ψ s (θ s 1 , θ s 2 )) = 0, φ ∆ u * (Ψ u (θ u 1 , θ u 2 )) − φ ∆ s * (Ψ s (θ s 1 , θ s 2 )) = 0 (48) for as many discrete values of θ as Fourier coefficients needed to be determined in the corresponding equation, and with for i = u, s. Note that a Taylor expansion of the previous expression around ϕ i θ i 1 − (θ i 2 /(2π))ρ i up to first order in ξ i turns it into an expression analogous to (44) except for an error O((ξ i ) 2 ), which is already the error of the linear approximation of the manifold. Compared to (44), expression (49) has as an advantage the fact that it does not contain the differential of the flow. The comments made for system (47) also apply here: system (48) also includes a normalization condition for the invariant bundles v u , v s to be locally unique, ξ i is a parameter that must be kept fixed at a small value (e.g. 10 −6 ), and an actual implementation requires multiple shooting, both in the tori and the connection. Additional details can be found in [33].