Efficient representation of invariant manifolds of periodic orbits in the CRTBP

This paper deals with a methodology for defining and computing an analytical Fourier-Taylor series parameterisation of local invariant manifolds associated to periodic orbits of polynomial vector fields. Following the Parameterisation Method, the functions involved in the series result by solving some linear non autonomous differential equations. Exploiting the Floquet normal form decomposition, the time dependency is removed and the differential problem is rephrased as an algebraic system to be solved for the Fourier coefficients of the unknown periodic functions. The procedure leads to an efficient and fast computational algorithm. Motivated by mission design purposes, the technique is applied in the framework of the Circular Restricted Three Body problem and the parameterisation of local invariant manifolds for several halo orbits is computed and discussed.


1.
Introduction. Invariant sets, as for instance equilibria, periodic orbits, invariant tori and associated stable/unstable manifolds are cornerstones in the study of a dynamical system. They provide insight into the orbit structure of the system and constitute a stylised skeleton around which the complicated dynamics evolves. For example, co-dimension 1 invariant manifolds form tubes which separates regions in the phase space and the intersections of stable and unstable invariant manifolds produce connecting orbits and, in some circumstances, establish the existence of chaotic motion.
In astrodynamics, invariant manifolds associated to equilibria and periodic or quasi-periodic orbits of the Circular Restricted Three Body problem (CRTBP) are considered in the design of energy efficient transfers. In order to reduce the payload of a mission, a target orbit can be reached at cost zero, theoretically at least, by placing the spacecraft on the stable manifold associated to the orbit. Also, a multi body system can be properly decomposed in three-body systems and the hyperbolicity of the invariant manifolds can be exploited to dynamically connect far apart regions in the phase space and to design efficient and complicated itineraries, see for instance [1,7,20,29,42]. The invariant manifold dynamics is also considered in other astronomical problem, as for instance to explain the mass transfer in the solar system [38] or the formation of arms and rings in galaxies [39]. Because of the large number of reasonable three-body systems that can be coupled for potential space mission scenarios, bounded orbits, invariant manifolds and low-energy transfers play a fundamental role in the design of new missions for the deep space exploration.
In this perspective, invariant manifolds are considered as building blocks to be combined to detect portions of trajectories. One strategy to compute the manifolds is to perturb some data of the invariant set (as the periodic orbit) in the direction of the stable/unstable tangent bundle and then to integrate the system. It results a collection of trajectories on the manifold. Any other point on the manifold can be computed by choosing a finer discretisation of the orbit and repeating the process for the new data. This procedure is not efficient and, in some situation, not even practical, see [25].
Several methodologies propose algorithms for globalizing invariant manifolds so to avoid the recomputation of single trajectories as needed [31]. For instance, approximation by collection of geodesics level sets [30], fat trajectories and cells [21] or box covering [14,15,46]. Also, spline interpolation techniques have been presented to approximate the manifold surface, [25,43]. Beside purely numerical techniques, it is worth mentioning that a variety of analytical and semi-analytical methods have been developed to represent the invariant manifolds. Lindstedt-Poincaré series and normal form of the Hamiltonian are employed to compute high order expansion of the centre manifold of the collinear points [28] and invariant manifolds for halo and Lissajous orbits in the circular and elliptic restricted three body problem, see [35,32] and references therein. See also [40] for a short survey.
Whatever technique is used, typically an algorithm for the design and optimisation of a space trajectory requires several evaluations of points on the manifolds and the efficiency of the computation would increase if a continuous and handy representation of the manifold is provided. The demand of an efficient and practical continuous parameterisation of the manifolds is the main motivation of this work.
We present an application of the so-called Parameterisation Method for the computation of a parameterisation of the invariant manifolds associated to periodic orbits in the CRTBP. In a series of paper of Cabré, Fontich, de la Llave [3,4,5], the Parameterisation Method has been presented as a tool to establish existence and regularity of manifolds associated to fixed points and periodic orbits for general dynamical systems. The method defines a functional equation, named invariance equation, for the parameterisation of the manifold and for the dynamics on it. The technique has been applied to a wide range of situations, as for instance for continuous and discrete maps [34,37], in finite and infinite dimension [36], for difference [13] and delay equations [33]. Both analytical and numerical aspects of the Parameterisation Method have been investigated, see also [23,24,26,27]. Moreover, the method is well suited for rigorous computational studies, as shown for instance in [6,45,17,18]. A complete review of the literature is out of the scope of this introduction, we refer to [22] for a more exhaustive discussion.
Following the approach presented in [9], where the Parameterization Method is combined with the Floquet theory, the goal of this paper is to present a fast and efficient algorithm for the computation of the coefficients of a Fourier-Taylor power series parameterisation of the stable and unstable manifolds of periodic orbits in the CRTBP. Since the periodic orbits we consider have one stable and one unstable direction, the manifolds are seen as the image of a function of two variables, say θ, σ. Imaging an invariant manifold as a tube in the phase space, the variable θ parameterises circles on the manifold, while the variable σ parameterises the longitudinal fibers along the tube. The analytical representation of the invariant manifold provided by the parameterisation is well suited and particularly convenient in the process of designing first guess trajectories for space mission. For practical missions, these trajectories necessitate a further optimisation in more complete and realistic system, as N -body problem or JPL ephemeris, therefore any precise analysis of the convergence and aposteriori error evaluation of the parametrisation is beyond the scope of this work.
We hasten to remark that no numerical integration is needed out of the computation of the orbit and of the variational system. Rather, the computation of the high-order coefficients of the parameterisation is performed completely in Fourier space and, thanks to the Floquet decomposition, the system reduces to a straightforward diagonal algebraic problem in the Fourier coefficients.
The paper is organised as follow: in section 2 we present all the theoretical aspects of the method for computing the Fourier-Taylor parameterisation of invariant manifolds for a general vector field. At first, in section 2.1 we recall the Floquet normal form decomposition of the fundamental matrix solution of linear systems of non autonomous ODEs. Then in section 2.2 we quickly review the Parameterisation Method for two dimensional invariant manifolds associated to periodic orbits. The invariance equation and the notion of conjugating chart map are then introduced. In section 2.3 we select a formal Fourier-Taylor series as conjugating map and we solve the invariance equation. That leads us to a sequence of non autonomous linear ODEs (homological equations) to be solved for the coefficients of the parameterisation. The existence of solution is also briefly discussed. Solving for the Fourier coefficients of the wanted periodic functions, the differential problems move into a sequence of algebraic systems. Taking into account the Floquet decomposition, we show how each algebraic system can be diagonalized and iteratively solved. The method is then rewritten in form of a numerical algorithm in section 2.4. Section 3 concerns the dynamical system governing the CRTBP. At first we remind the classical formulation in the rotating reference frame. Then we provide an equivalent polynomial vector field formulation.
In section 4 we apply the technique and we collect the numerical results for the parameterization of invariant manifolds associated to halo orbits. Also we present an application to the problem of computing heteroclinic connections and a optimal transfer from halo to science orbit.
2. Parameterisation of invariant manifold: Theoretical formulation. In this section we describe the theoretical framework of the method for computing Fourier-Taylor parameterisation of invariant manifolds associated to hyperbolic periodic orbits of vector fields.
Following the Parameterisation Method, introduced in [3,4,5], the idea is to conjugate the dynamics on the invariant manifold to a linear dynamics on the parameter space and then to study an invariance equation that characterises chart maps for the local manifold. By inserting a certain formal series ( in our case a Fourier-Taylor series) into the invariance equation, the coefficients of the parameterisation are given recursively as periodic solutions of a sequence of ordinary differential equations with periodic coefficients.
If the eigenvalues of the Monodromy matrix of the periodic orbit satisfy a certain non-resonance condition, in [5] it has been shown that these linear equations admit a solution and an explicit expression is provided. Also, in the quoted papers, the convergence of the series is formally discussed.
To compute the solution of the invariance equation we follow the approach presented in [9]. The unknown periodic functions are expended in Fourier series and the linear differential equations are rephrased as algebraic problems. Exploiting the Floquet theory, the unknown Fourier coefficients can be computed by solving a straightforward equation. That allows to easily compute the coefficients of the parameterisation up to any desired finite order.
In the following we discuss with more details the technique above outlined. We begin the exposition by recalling some fundamental results of the Floquet theory, see [12] for more details.
2.1. Floquet theory. Consider the non autonomous linear system of ordinary differential equationsẋ = A(t)x, where A(t) is a continuous T -periodic n×n matrix valued function. Since the system is linear, the solutions of (1) form a vector space. The principal fundamental matrix solution is the matrix function Φ(t) that solves the initial value problem where I is the n × n identity matrix. From now on we omit the term principal and we refer to Φ(t) as the fundamental matrix solution. The fundamental matrix solution alone determines all the solutions of (1) in the sense that the solution x(t) with initial data x(0) = x 0 is given by The Floquet theorem provides a canonical form for the fundamental matrix solution. Here we recall the complex form.
Theorem 2.1 (Floquet normal form decomposition). Let A(t) be a T -periodic continuous matrix valued function and denote by Φ(t) the fundamental matrix solution of (1). Then there exist a possibly complex matrix B and a possibly complex nonsingular T -periodic matrix function P (t) such that Φ(t) = P (t)e tB for all t ∈ R.
The representation Φ(t) = P (t)e tB is known as the complex Floquet normal form.
The matrix Φ(T ) is called Monodromy matrix and the eigenvalues µ i of the Monodromy matrix are called characteristic multipliers of the corresponding timeperiodic system (1). A complex number λ is called a characteristic exponent if µ is a characteristic multiplier and e λT = µ.
Assume that B is diagonalisable and denote by {λ i , v i } i the eigenvalues and eigenvectors of B. In the next lemma some properties of Floquet's normal form decomposition are collected. Suppose that γ(t) is hyperbolic and denote by W u loc (Γ), W s loc (Γ) the local unstable and stable manifold of Γ, respectively. As it is well known, the dimension of the invariant manifold as like as the linear approximation of the manifold are provided by the variational equation where Df is the Jacobian of the vector field. Denoting by M 0 = Φ(T ) the monodromy matrix, the number of eigenvalues of M 0 with modulo less (larger) than one prescribes the dimension of the stable (unstable) manifold, and the corresponding eigenvectors v i provide the linear approximation of the manifold at γ(0). Throughout this section, we discuss only the stable manifold. Consider the Floquet normal form decomposition Φ(t) = P (t)e Bt and assume that B is diagonalisable. For later purposes, we remind that P (t) and B solve the equationṖ = Df (γ(t))P − P B.
For the purposes of this work, we restrict to the case that B has only one stable eigenvalue λ = λ 1 and hence one stable direction v = v 1 . Indeed, we are interested in the invariant manifolds of periodic orbits in the CRTBP which have one stable and one unstable direction.
For Lemma 2.2 point 3) v is an eigenvector of M 0 and represents the direction at γ(0) tangent to the stable manifold. As shown in [8], the vector v(θ) := P (θ)v is the direction tangent to the manifold at γ(θ), for any θ ∈ [0, T ]. Therefore the function P 1 (θ, σ) = γ(θ) + v(θ)σ parameterises the stable tangent bundle of Γ. The aim of the following discussion is to define a nonlinear correction to P 1 that parameterises the local stable manifold. Remark 1. In the astrodynamical community it is common to refer to the operator Φ(t) as the State Transition Matrix (STM) and the stable direction along the orbit is usually computed asṽ(θ) = Φ(θ)v. Since Φ(θ) = P (θ)e θB and v is an eigenvector of B, it holdsṽ(θ) = e λθ P (θ)v = e λθ v(θ). That is, the two vectorsṽ(θ) and v(θ) are proportional and both can be used to span the tangent space at γ(θ). However, we notice that while v(θ) is a periodic function, that same is not true forṽ(θ). Therefore v(θ) = P (θ)v should be preferred in parameterising the stable bundle, that is geometrically a periodic object.
The variable θ and σ are called the parameters and we refer to the cylinder C T,ν := T T × B ν as the parameter space, where T T = [0, T ] /{0,T } is the circle of length T and B ν := {σ : |σ| < ν}. Consider the linear vector field on C T,ν θ = 1,σ = λσ.
As already mentioned, the purpose is to define a map P : C T,ν → W s loc (Γ) that parameterises the local stable manifold and conjugates the dynamics on the manifold with the linear vector field (5) on the parameter space. Definition 2.3 (Conjugating map for the local stable manifold). A function P: C T,ν → R n is said a conjugating map for a local stable manifold of the periodic orbit Γ if 1. P is a continuous, injective mapping of the cylinder C T,ν which is real analytic on the interior of the cylinder.

The conjugacy relation
is satisfied for all θ ∈ T T , σ ∈ B ν , and t ≥ 0.
The above condition implies that if P is a conjugating map then that is, any P(θ, σ) is a point on the stable manifold of Γ. See also Figure 1 for a visualisation of the relation (6). The function P(θ, σ) maps the cylinder in parameter space into the local manifold, in such a way that trajectories of the linear vector field (5) are mapped into trajectories on the manifold.
Let P(θ, σ) be a conjugating map. If we differentiate on both sides of (6) we obtain The following theorem states that the above condition is also sufficient for P(θ, σ) to be a conjugating map, see [5].
Theorem 2.4 (Invariance equation for a conjugating chart map). Suppose that P: C T,ν → R n is a continuous function such that P(θ, 0) = γ(θ), and P(θ, σ) solves the partial differential equation on the interior of C T,ν . Then P is a conjugating chart map for W s loc (Γ) in the sense of Definition 2.3. It follows from (7) that image(P) is a local stable manifold for Γ.
The invariance equation is the main object of our study, since it characterises a conjugating map and hence provides a parameterisation of the local invariant manifold. In the next section we quickly discuss the existence of solutions of equation (10) and we develop a formal power series solution together with a constructive algorithm for computing an approximation.

2.3.
Fourier-Taylor parameterization. We look for a solution of (10) in the form with a α,k ∈ R n for any α, k. We refer to this expansion as a Fourier-Taylor series and we name θ and σ as the Fourier and Taylor parameter respectively. In practice, the unknowns are the periodic functions a α (θ) or, equivalently, the sequence of Fourier coefficients {a α,k } k of any α. We insert (11) into the invariance equation and expand f (P(θ, σ)) as its Taylor series with respect to σ at σ = 0. We obtain where f (α) (P(θ, 0)) is the derivative d α dσ α (f (P(θ, σ))) evaluated at (θ, 0). Now, by matching the powers, we obtain d dθ a 0 (θ) = f (a 0 (θ)) whose solution is clearly the periodic orbit, hence a 0 (t) = γ(t) and d dθ The function a 1 (θ) is the tangent bundle, that we already said to be v(θ) = P (θ)v.
Indeed we now show that P (θ)v solves the above equation. For (4) d dθ since v is an eigenvector of B with eigenvalue λ. Going back to (12), matching the powers σ α with α ≥ 2, we obtain the equations d dθ where we separated the contribution of a α (θ) in 1 α! f (α) (P(θ, 0)), so that R α (θ) only depends on lower order terms. We refer to (13) as the homological equation at level α.
For any α ≥ 2 we iteratively solve the homological equation, hence providing one after the other the coefficients a α (θ) of the parameterisation. Before, let us recall a theorem ensuring sufficient conditions for the existence and uniqueness for T -periodic solution of the homological equation, see [5].
Theorem 2.6. If e αλT is not an eigenvalue of the monodromy matrix Φ(T ), then for any T -periodic function R α (θ) there exists a unique T -periodic solution a α (θ) of equation (13). An expression for the solution is given by where a α (0) solves Remark 2. In the situation we are interested in, the condition that e αλT is not an eigenvalue of the monodromy matrix is always satisfied, since the periodic orbits in the CRTBP have only one stable and one unstable direction. However, in different situations a periodic orbit may have more stable directions and the eigenvalues of B may be rationally dependent. If that happens the eigenvalues are said to be resonant and the construction of the parameterisation here discussed fails. In studying manifolds for stationary solutions, a variation of the method that takes into account resonances is proposed in [44].
We are now in the position to solve for any α the homological equation and, meanwhile, to present the algorithm to construct layer by layer the coefficients of the parameterisation. A possible approach to the problem could be to obtain a α (0) from relation (15) and hence compute a α (θ) from (14). From the computational point of view this approach is not completely straightforward and it returns the aimed solutions a α (θ) only at some grid points of the Fourier parameter θ. Rather, using the Floquet normal form decomposition, a different scheme that allows to directly compute the Fourier coefficients of the function a α (θ) is here presented.
Suppose it is given the Floquet decomposition Φ(t) = P (t)e Bt of the fundamental matrix solution of the variational system (3). Assuming that B is diagonalizable, denote M the matrix such that For any |α| ≥ 2 define the functions w α (θ) : T T → R n by the coordinate transformation Taking into account (4), the homological equation (13) is transformed into the constant coefficient ordinary differential equation Consider the Fourier expansion of w α and Then equation (16) gives The latter equation is diagonalized and solving component- for 1 ≤ j ≤ n. Working back from v α,k we obtain w α,k and the Fourier coefficients of the desired solution a α (θ) as the Fourier coefficients of P (θ)w α (θ),

2.4.
Algorithm for the computation of the parameterisation. In this section we summarise the construction of the parameterisation theoretically described before and we provide the algorithm that will be implemented in numerical experiments. Clearly, we can only compute a finite order parameterisation First of all, we notice that we used the operator P −1 (θ). Since P (θ) is non-singular, the inverse P (θ) −1 is well defined for any θ. Also, P (θ) −1 solves a differential problem similar to (4). Indeed, from the Floquet normal form decomposition Φ(θ) = P (θ)e Bθ it follows P −1 (θ)Φ(θ) = e Bθ . Then, differentiating on both sides, we obtain ( Thus, P −1 is solution of the system Taking the transpose on both sides of the differential equation and defining S(θ) = that is exactly the same system as the one for P (θ) with the data H(θ) in place of Df (γ(θ)) and W in place of B.
The recipe to construct the parameterisation of the local invariant manifold is the following.
1. Compute the T -periodic orbit γ(t) of the vector field equationẋ = f (x). 2. Integrate the variational equation (3) for one period to obtain Φ(T ). Compute the matrix B as B = log(Φ(T )) T . 3. Compute P (t), t ∈ [0, T ], by integrating the equation (4) for one period and compute P −1 (t) from (18). Choose m, the number of Fourier coefficients for each term in the Fourier-Taylor series. 4. From the time discretisation of γ(t), P (t), P −1 (t), compute the Fourier coef- Choose N , the order of the finite series. Choose λ among the eigenvalues of B, and v the associated eigenvector. 5. Compute a 1,k = P k v, for any |k| ≤ m. 6. For any α ranging from 2 to N : a) compute R α,k and A α,k with |k| ≤ m, respectively the Fourier coefficients of R α (θ) and P −1 (θ)R α (θ). b) compute a α,k for |k| ≤ m as showed in the last part of the previous section. We now recall some fundamental issues concerning the Fourier sequence space, useful when implementing the algorithm.
Remark 3. Given two periodic functions, f (t) and g(t) and the associated sequences of Fourier coefficients F = {f k }, G = {g k }, the sequence of Fourier coefficients of the product f (t)g(t) is given by the convolution product F * G, where This operation extends to any number of factors: From the computational point of view, the convolution product can be easily and quickly performed by means of FFT and its inverse. Moreover, as showed in [18], the Fourier Transform can be used to control the errors of the convolution products. From Remark 3 one realises that the coefficients A α,k , at step (vi)-a) in the algorithm can be easily computed once the sequence of Fourier coefficients {R α,k } k is known. As outlined in Remark 4, if f (x) is a polynomial vector field then the functions R α (θ) are given as products of a i (θ), with i < α. Therefore the Fourier coefficients of R α are rapidly computed as convolutions of the known sequences {a i,k } k .
To conclude this section, we note that after performing one time the points (i) (ii) (iii) of the algorithm one obtains different approximations for the parameterisation of both the stable and unstable local manifold, by changing the parameter m and N and also the eigenvalues λ and v. Steps (iv)-(vi) of the algorithm only involve algebraic manipulation of the coefficients and can be quickly executed. It is also evident that the required ingredients for the execution of the algorithms are the eigenvalue λ with associated eigenvector v and the Fourier sequences of γ(t), P (t) and P −1 (t). Different approaches, for instance by integration, fixed point argument or degree theory, can be followed to compute the periodic orbit and the Floquet decomposition. Any deeper discussion in this direction is out of the scope of the present work.
3. Dynamical model: The circular restricted three-body problem. In the Circular Restricted Three-Body problem (CRTBP), a massless particle, say p, moves under the gravitational influences of two massive bodies, called primaries, with mass m 1 ≥ m 2 . Due to their mutual gravitational attraction, the primaries are revolving with constant angular velocity in circular orbits around the centre of mass of the system.
By introducing the velocity vector v = (v x , v y , v z ), system (19) is equivalent to the first order vector field equatioṅ The following are well known facts about the CRTBP, [19].
• The CRTBP admits an integral of motion, called Jacobi constant The family of 5-dimensional energy manifolds is a foliation of the six-dimensional phase space. • The system has five equilibrium points L k , called Lagrangian points. Three of them, L 1 , L 2 , L 3 are on the x-axis and represent collinear configurations, whereas L 4 , L 5 correspond to equilateral configurations of the masses. All of them are placed on the z = 0 plane. • Linear stability analysis shows that the collinear Lagrangian points exhibit a center×center×saddle behaviour. From each equilibrium point bifurcate families of periodic planar and vertical Lyapunov orbits and families of threedimensional periodic halo orbits. The monodromy matrix associated to an halo orbit has two eigenvalues |µ u | > 1, |µ s | < 1 related to expanding and contracting directions, µ u µ s = 1; two eigenvalues µ 3 = µ 4 = 1, related to the direction tangent to the orbit and to the variation of energy; two complex eigenvalues µ 5 and µ 6 of modulus 1. Therefore, each of these periodic orbits has a two dimensional stable and unstable invariant manifold, i.e. at any point of the orbit there is a one dimensional attracting and expanding direction. The same dynamical behaviour holds for the Lyapunov orbits.
3.1. Polynomial vector field formulation. As we outlined in Remark 4 and subsequent discussion, the proposed algorithm for the computation of the parameterisation of the invariant manifolds is well suited when the nonlinearity of the vector field is polynomial. Clearly, system (20) is not polynomial because of the division by r 3 1 and r 3 2 . To overcome this hurdle, we now transform system (20) into an equivalent polynomial system. Let us introduce two more functions: Differentiation with respect to t giveṡ Inserting into (20), we obtain a new vector field equatioṅ If we impose the constraint solutions of equation (21) subjected to (22) are equivalent to solutions of system (20). Denote by X(t) = (x, y, x, v x , v y , v z , p, q)(t) the vector of unknowns and F (X) the right hand side of (21). Suppose now γ(t) = (x, y, x, v x , v y , v z , p, q)(t) is a Tperiodic solution ofẊ = F (X). The projection on the first three components of γ(t) gives the periodic orbit for the CRTBP system (19). Note that F (X) is a polynomial vector field. We use a α = {a α,k } k to denote the sequence of vectors of Fourier coefficients and similarly a (j) In our algorithm for the parameterisation of the manifolds we first have to compute the linearised system. The Jaobian of F (X) is explicitly given as . Then, in the computation of the high order terms a α with α ≥ 2 we need to compute R α . We remind that R α only involves lower order terms, that is a α with α < α. According to Remark 3 and Remark 4, in our situation the sequence R α is given by 4. Numerical results. For any numerical computation we report the mass ratio µ, the vector of initial conditions of the periodic orbit [x(0), z(0), v y (0)], ( having y(0) = v x (0) = v z (0) = 0 and p(0), q(0) as in (22)), the period T and λ u , λ s the unstable and stable eigenvalues of the matrix B. Since the stable/unstable eigenvalues of the monodromy matrix satisfies µ s µ u = 1 and µ s,u = e λs,uT , it follows that the λ u > 0 and λ s = −λ u . Also, we provide the Fourier parameter m and the Taylor order N used in computing the parameterisation P m,N (θ, σ).
For the numerical integration of the orbit and of the variational system we used a variable order Adams-Bashforth-Moulton method with relative and absolute tolerance of 10 −14 , 10 −16 respectively.  In Figure 2 we plot the image of the parameterization for both the local stable and unstable manifold. In Appendix we report the first few Fourier coefficients of the functions a α (θ) with α = 0, . . . , 5.
Example 2. Figure 3 shows   We present a quick analysis of the output of the computation. Recall the notation: a α,k denotes the k-th Fourier coefficient of the function a α (θ) and a α = {a α,k } k denotes the sequence of Fourier coefficients for a given α. For α = 1, . . . , N we compute the norm |a α,k | of the vector a α,k for any k = 0, . . . , m. We expect that |a α,k | decreases to zero as k increases. That indicates that the Fourier series converge and a α (θ) is smooth.
In the left of Figure 5 we plot the values of log(|a α,k |). Different lines correspond to different α. We note that the lines are almost linear, showing that the Fourier coefficients are exponentially decreasing and hence providing evidence that the functions a α (θ) are analytic.
Also, we measure the norm of each a α as a α = |k|≤m |a α,k |. Roughly speaking, this norm corresponds to the L 1 norm of the function a α (θ). The parameterisation is supposed to converge for σ small enough and, being a power series, the domain of convergence depends on how fast the functions a α (θ) decay.
In the right of Figure 5 the values of a α are plotted for α = 1, . . . , N , in logarithmic scale. Again we notice that a α decays exponentially hence showing insight about the formal convergence of the series.   and compute that parameterisation for both the stable and unstable local manifolds, whose image are plot in Figure 6.
We want to measure the quality of the parameterisation, that is to check whether the image of the computed parameterisation agrees with the manifold. For this we compute part of the unstable manifold by integrating (21) starting from some data on the linear approximation of the manifold and we plot the manifold so obtained on the top of the image of the parameterisation, see Figure 7 (left). We note that the two pictures overlap almost perfectly, in particular close to the periodic orbit one can not even distinguish the two images. Then, we define two planes Π 1 : −x + 0.15y − 0.9z + 0.987 = 0 and Π 2 : −x + 0.15y − 0.9z + 0.9858 = 0, see Figure 7 (right), and we compute the Poincaré sections of both the computed manifolds. The red points in Figure 8 denote the intersections of the fibers of the manifold obtained by integration with the planes, while the blu line is the Poincaré section of the image of the parameterisation. At Π 1 , the two Poincaré sections coincide proving that the parameterisation method provides reliable results. On Π 2 the two sections are still quite close, but the blu line is open. That occurs Figure 6. Image of the parameterisation of the stable and unstable local manifold associated to the orbit of Exemple 6.
because the image of the parameterisation starts to fold and curl and part of the image does not hit the plane. Hence we can say at the Π 2 the parameterisation is not trustworthy anymore and we are out of the domain of convergence of the parameterisation series. To conclude, Table 1 reports the time required to compute the higher order terms of the parameterisation for different choices of m and N in Matlab 2014a (on a MacBook Pro, Intel Core i7, 2.8 GHz), that is to perform points 5 and 6 of the algorithm described in section 2.4.

4.1.
Application. Invariant manifolds structure are surely of primary interest for mission design purposes. In many situations, part of the design strategy consists in finding the insertion point on the invariant manifold that minimizes a certain objective functional. As a paradigm, if one is interested in approaching a periodic halo orbit, a suitable insertion point on the stable manifold is seek. Similarly, a transfer from a periodic orbit to a science orbit can be seen as an optimal trajectory on the unstable manifold. In both cases the objective functional is a function of the points on the manifold. If the manifold is given only in terms of a discretisation, it might be difficult if not impossible to find the best trajectory. Rather, having the parameterisation to represent the manifold, the objective functional results to be a smooth function on the two parameters θ, σ and semi analytical minimisation techniques may be applied to find (local) optimality. The following examples are meant to illustrate the potentialities and possible applications of the parameterisation of the manifolds, not to add any new scientific result.
Consider the problem of computing an heteroclinic connection between two periodic orbits, say γ 1 , γ 2 . Denoting by ϕ(x, t) the flow associated to the CRTBP, the problem consists in finding a point x 0 on the unstable manifold of γ 1 and a point x f on the stable manifold of γ 2 such that ϕ(x 0 , T ) = x f for a suitable T . Having the parametrisation of the unstable manifold of γ 1 , P 1 m,N (θ 1 , σ 1 ), and of the stable manifold of γ 2 , P 2 m,N (θ 2 , σ 2 ), the problem can be rephrased as follows. Fix σ 1 =σ 1 , σ 2 =σ 2 and look for (θ 1 , θ 2 , T ) such that  The next application consists in finding a trajectory joining an halo orbit in the Jupiter-Europa CRTBP (µ = 2.5280176e − 05) to a science orbit around Europa, as considered in [2]. Denoted by γ an halo orbit around L 2 , we aim at connecting γ to a circular orbit around Europa at altitude r. A single impulsive manoeuvre is allowed to insert the spacecraft into the circular orbit once the required altitude is reached. The design strategy is as follows: the spacecraft leaves the halo along a trajectory on the unstable manifold and if a point at altitude r around Europa is reached, the necessary ∆V manoeuvre to circularise the orbit is computed. The optimal trajectory is the one that minimises the required ∆V . Denote by (r, v)(t) the position and velocity vectors in the synodical reference frame. At the moment when r(t) reaches the sphere of radius r around Europa convert the state (r, v) into the Europa-centric inertial frame (R, V ) = (r The required ∆V is given as where v = (µ/r) 1/2 is the modulus of the circular orbit, v 1 = |V | and φ is the angle between V and the tangent to the sphere of radius r. Let P m,N (θ, σ) be the parameterisation of the unstable manifold of γ. Fix a value of σ =σ. Then each trajectory on the unstable manifold is uniquely determined by θ. It results that the required ∆V is a function of θ and a numerical scheme to can be employed to compute (local) minimiser of ∆V . Thanks to the continuous parameterisation P m,N (θ,σ) of the section of the unstable bundle, it is not necessary to compute any trajectory on the unstable manifold close to the periodic orbit at each stage of the minimisation process. In our computation we consider the halo orbit with initial state x = [1.0175667, 0, 0.0222484, 0, −0.0306484, 0] and the target science orbit at altitude r = 200 km. The minimisation scheme returns a minimal ∆V = 0.43 km/s. (In this computation we have considered only those trajectories on the unstable manifold that intersect the sphere of radius r around Europa at the first approach to the moon.)

5.
Conclusions and future plan. In this paper a method to compute the Fourier-Taylor parameterisation of the local invariant manifolds associated to periodic orbits of the CRTBP is presented. Taking advantage from the Floquet normal form decomposition of the fundamental matrix solution, the time dependency in the homological equation is removed, and the functions involved in the power series parameterisation are found in terms of their Fourier coefficients. Some algebraic manipulation allows to diagonals the system in the Fourier space.
The output of the algorithm consists in sequences of Fourier coefficients a α,k for 0 ≤ α ≤ N and |k| ≤ m. Combining the Fourier coefficients, a function of two parameters is produced, After rephrasing the CRTBP system as a polynomial vector field, the computational algorithm has been applied to several halo orbits. Since no numerical integration is needed out of the computation of the orbit and of the variational system, the higher order terms of the parameterisation are quickly approximated. The image of the parameterisation results to be reliable, at least close enough to the orbit, even with few nonlinear corrections, i.e. with N = 4 or N = 5.
Since the function P m,N (θ, σ) is continuous, indeed smooth, and parameterises the invariant manifold, it can be useful for mission design purposes. Besides the couple of applications above presented, one can take advantage from the parameterisation anytime the objective functional is a function of the points on the manifold. Having an analytical representation the manifold, the objective functional results to be a smooth function on the two parameters θ, σ and semi analytical minimisation techniques may be applied. For instance, the parameterisation of invariant manifolds for several periodic orbits can be considered at the same time and combined in the design of multi-moon optimal trajectories [16].
Lastly, we remark one more time that the parameterisation method, in particular in the form of a Fourier-Taylor parameterisation, is well suited for rigorous computational analysis. The aposteriori error evaluation together with a contraction argument on the N -tail space may lead to a computer assisted enclosure of the local invariant manifold, see [10,11]. Thus, a projected boundary value problem can be set up to prove existence of connecting orbits.