Two Periodic Models for the Earth-Moon System

This paper discusses two alternative models to the Restricted Three Body Problem (RTBP) for the study of a massless particle in the Earth-Moon system. These models are the Bicircular Problem (BCP) and the Quasi-Bicircular Problem (QBCP). While the RTBP is autonomous, the BCP and the QBCP are periodically time dependent due to the inclusion of the Sun’s gravitational potential. Each of the two alternative models is suitable for certain regions of the phase space. More concretely, we show that the BCP is more adequate to study the dynamics near the triangular points while the QBCP is more adequate for the dynamics near the collinear points.


Introduction
During the last years, the scientific community has increased its interest in the natural motions occurring in the Earth-Moon system. The list of possible applications is vast, for instance: the study of the far side of Moon and the relation with the translunar point L 2 ; the aim to exploit the cis-lunar space and the convenience of using the invariant structures related to L 1 .
We have mentioned a couple which are specifically related to the Lagrangian points but, obviously, the list goes on covering a wide range of interests. Efficient mission designs depend ultimately on the understanding of the natural dynamics. To fulfill this goal, it is advisable to use simplified models. Simple models allow us to understand the underlying mechanisms that lead to interesting phenomena. From the dynamical systems point of view, the comprehension of the invariant structures and their stability of simple models has helped to shed light on difficult problems such as the motion of asteroids through the solar system, station keeping of spacecraft and taking advantage on the natural dynamics to design spacecraft missions. Perhaps the most illustrative example for the purpose of this work is the existence of the Trojan asteroids that can be predicted using the effective stability of the triangular points of the Sun-Jupiter Restricted Three Body Problem. This example is convenient for the purposes of this work because the existence of objects in the triangular points has a counterpart in the Earth-Moon system: the Kordylewsky clouds. We shall come back to this example during this work (Section 4.3)but, for the moment being, we want to stress that the existence of these clouds cannot be established by using the same theoretical mechanisms as the Trojan asteroids [GJL05,PL15,PE17]. In fact, the literature related to the Kordylewsky clouds has been stumbling around the existence or nonexistence of objects in the Earth-Moon triangular points, mostly because of the lack of observations. Therefore it is convenient to analyze whether a simple model is suitable for the problem we want to study.
The Earth-Moon Restricted Three Body Problem (RTBP) is the most used simple model for the motion of a small body in the Earth-Moon system. There is, however, a remarkable number of works that take into account the presence of Sun's gravitation (see, for instance, [GJMS91b,GJMS91a,GJMS93,GLMS01,GLMS85]). The most relevant effect ignored by the Earth-Moon RTBP is the gravitational attraction of Sun. In this respect, a simple model to study the dynamics near the Earth-Moon system needs to take into account Solar gravity. The problem has a natural non-autonomous periodic time dependence formulation. An advantage of the periodic models is that they can be handled by means of a stroboscopic map i. e. the map defined by the evaluation of the flow at the period of the vectorfield. This is crucial because, while the complexity of the system increases, the study of maps (even if they are numerically defined) is, in some aspects, more comfortable than the study of flow. In periodic time dependent systems, the simplest invariant objects, the ones the dynamics is organized from, are the periodic orbits with the same period as the vectorfield. These periodic orbits appear as fixed points of the stroboscopic maps and their robustness is assured by the classical Implicit Function Theorem. We would like to remark that, in quasi-periodic models the simplest invariant objects are invariant tori. The computation and study of these objects is more difficult. The discussion in this paragraph vindicates a closer look to periodic models for the Earth-Moon system. We have selected two among the literature, the Bicircular Problem and the Quasi-Bicircular Problem. Both models include Sun's gravity and can be written as periodic perturbations of the RTBP.
The Bicircular Problem (BCP) is a restricted four body problem [Hua60,CRR64]. There are three primaries and a fourth, massless, test particle. In our case, the three primaries are Earth, Moon and Sun. However, this model has been utilized in other cases [BGMO16]. It is assumed that Earth and Moon move as in the RTBP, that is, along a circular orbit around its common centre of masses. Let us name C EM this barycentre. Name C SEM the centre of masses of the Sun-C EM system. As Moon and Earth move, it is assumed that Sun and C EM move in another circular orbit around C SEM . We refer to [GJMS91c] for a detailed derivation of the equations of motion. The (BCP) is a periodic perturbation of the RTBP that takes into account the direct gravitational effect of a third primary (in our case, Sun) on the particle. This model captures the non-stable character of the triangular points. Henceforth, it is suitable to use it when studying problems related with these locations (see, for instance, [CJ00,Jor00]). A remarkable shortcoming of the BCP is its lack of coherence i. e. the motion assumed for the primaries does not verify Newton's laws. Moreover, the BCP has no translunar dynamical structure. This justifies the seek for a more complex model for the study of, at least, the L 2 point.
The Quasi-Bicircular Problem (QBCP) is a version of the four body problem. It is conceived to be a coherent counterpart of the BCP. This model was introduced by C. Simó, and it has been used in several works, see [And98,AS00,And02] and, more recently, [BMGLD17]. A characteristic of the BCP is the lack of coherence of the bicircular solution assumed for the primaries. However, there exist solutions of the three body problem which are close to bicircular [SM71]. To build the QBCP it is necessary to compute a quasibicircular solution of the three body problem, in this case, for the Earth-Moon-Sun case. There are several ways to do such a thing. In [And98, AS00, And02] the authors build a specific algebraic manipulator and compute directly the Fourier coefficients of the quasibicircular solution. In [Gab03, GJ01, GJR04] the authors use a continuation scheme to compute the desired solution starting from a solution of the two body problem. After that, a Fourier transform is applied to get the Fourier coefficients of the solution. The QBCP is suitable for the study of the collinear points, especially L 1 and L 2 .
With this paper, we aim to provide a general insight about the dynamics of these models for a particle in the Earth-Moon system. We care about (practical) stable motion near the triangular points and, to do so, we use the BCP. We also study invariant manifolds related to the collinear points in the QBCP. We believe that the value of this work is, precisely, giving a wide perspective and help the interested reader to choose a suitable simple model to face a first exploration related to a problem concerning the Earth-Moon system.
The paper is organized as follows: Section 2 is devoted to a brief description of the RTBP. We explain how the phase space near the Lagrangian points is organized referencing some remarkable works and mentioning the techniques used to study the problem. In Section 3 we give some words on the stroboscopic maps and periodic time-dependent Hamiltonian systems. We also explain how to compute high order unstable manifolds re-lated to fixed points using the parameterization method with single and parallel shooting. Section 4 describes how the BCP can be used to study the motion near the triangular points. The advantage of this model with respect to the RTBP is that it captures the unstable character of the triangular points in the real system. The results presented are mainly devoted to stable motion in an extended vicinity of the triangular points. In Section 5 we describe results concerning the QBCP. We focus on the unstable manifolds related to the periodic orbits that replace the collinear points. Finally, Section 6 is devoted to conclusions and Section 7 to technical details.

Restricted Three Body Problem
The (Circular) Restricted Three Body Problem (RTBP) is a simplified model for the motion of a massless particle under the gravitational attraction of two massive bodies, the so-called primaries [Sze67]. The primaries are assumed to revolve along circular orbits around their common centre of masses. It is usual to take units of masses so the sum of the masses of the primaries is equal to one. The units of length are taken so the distance between the primaries is equal to one and the units of time are taken so the period of the revolution of the primaries is equal to 2π. It is also standard to take a rotating frame of reference that fixes the primaries at the horizontal axis. The RTBP is an autonomous Hamiltonian system with three degrees of freedom. The Hamiltonian function writes as: where r 2 P E = (x − µ) 2 + y 2 + z 2 and r 2 P M = (x − µ + 1) 2 + y 2 + z 2 . The parameter µ is called the mass parameter and it is the mass of the smallest primary. In the case of the Earth-Moon system µ ≈ 0.012. It is well known that the RTBP has five equilibrium points. Three of them, the collinear points, are located in the horizontal axis. The other two, the triangular points, are located at the third vertex of an equilateral triangle whose other two vertices are the position of the primaries.
The only integral of motion of the RTBP is its Hamiltonian. In many texts, this integral of motion is presented under a slightly different form as the Jacobi integral. Each surface level of this integral is a five dimensional manifold. If the velocities are set to zero, this defines the so-called Zero Velocity Surface. These surfaces separate the configuration space in different regions. The trajectories of the system cannot cross the boundary between two of these regions. The shape of these regions change with the value of the Jacobi integral. As the Lagrangian equilibrium points are critical points of the Jacobi integral, the topology of these regions change when the energy value crosses the value associated to one of the Lagrangian points (for more details, see [Sze67]).
The three collinear points are of type saddle×centre×centre. This means that, under generic conditions, a 4-dimensional centre manifold emerges from each of these points. These manifolds are tangent to the elliptic eigendirections at the points. There exist, as well, one dimensional stable and unstable directions emerging tangentially to the hyperbolic eigendirections. Moreover: • By the Lyapunov Centre Theorem [MH92], two families of periodic orbits, the Lyapunov families, emanate from the equilibria. One of the families is born tangent to the (z, p z ) plane so it is called vertical family. The other family is contained in the (x, y, p x , p y ) plane and it is called horizontal family. One can parametrize each families by the amplitude of the orbits. The horizontal families related to L 1 and L 3 can be continued up to trajectories which collide with Earth. The horizontal family related to L 2 can be followed up to collisions with Moon. The vertical families end up in bifurcating planar orbits [Mon01,GJMS91c].
• The Lyapunov families can be regarded as the non-linear continuation of the harmonic oscillator given by each elliptic direction of the linearization around the equilibria. When the amplitude tends to zero, the frequency of the family tends to the normal modes of the equilibria.
• As the frequency varies, the horizontal family undergo a 1 : 1 resonance and the Halo [FK73, BB79, GJMS91c, CPS15] family are originated (by means of a pitchfork bifurcation). Secondary families of Halo-type orbits appear by duplication and triplication of the main family [GM01,Bro68].
The centre manifold can be computed by means of normal form techniques [Jor99,JM99b,JM99a], with the parametrization method [HL05a,HL05b,FJ10a,FJ10b] and also numerically [GM01]. The dynamics restricted to the centre manifold can be described by a Hamiltonian with two degrees of freedom. By fixing a level of energy and taking a Poincaré section, one can reduce the problem to the study of a family of area preserving maps. This methodology suffices to observe the phase space during the bifurcation that give rise to the Halo families.

Motion near the Triangular points
The Earth-Moon triangular points of the RTBP are linearly stable [Sze67]. KAM theory can be used to establish the existence of a dense set of Lagrangian invariant tori close enough to the equilibria [MS86]. This has important consequences on the nonlinear stability of the triangular points. If we restrict ourselves to the planar case, these KAM tori (of dimension two) act as barriers for the dynamics in a fixed level of energy. Therefore, KAM tori enclose stable motion for initial conditions which are close enough to the triangular points. This argument based on KAM theory falls apart in the spatial case. Indeed, Lagrangian tori have, in that case, dimension three and the phase space, for a fixed level of energy, is five dimensional. There is, in general, no way to avoid Arnold diffusion [Arn64] . However, using normal form techniques, it is possible to derive bounds on the diffusion time [GDF + 89] (these techniques can be extended to the periodic [JV98] and the quasi-periodic case [GJL05]). These, make us think about regions of practical stability i. e. regions in which the motion is non-stable but initial conditions take a long time, maybe longer than the expected age of the solar system, to escape. These theoretical results are valid for a small region near the triangular points and numerical simulations provide evidences of large regions of practical stability [SSST13]. It is natural to look for other invariant structures that play a remarkable role to define the shape of the (numerically computed) region. In this regard, [SSST13] provides numerical evidence on the role of the centre-unstable and centre-stable manifolds of the collinear point L 3 . These manifolds are of dimension five and act as barriers of the dynamics. Obviously, the motion driven by these manifolds escape from the vicinity of the triangular points at some moment, but, again, the required time to do so can be large.

Periodic time-dependent Hamiltonians and stroboscopic maps
The alternatives to the Restricted Three Body Problem presented in this paper, the Bicircular Problem and the Quasi-Bicircular Model are both periodic time dependent Hamiltonian systems that can be seen as perturbations of the RTBP. Because of this periodic time-dependence, the Lagrangian points are no longer equilibria but they are replaced by minimal periodic orbits i.e. periodic orbits with the same period as the perturbation. These periodic orbits are known as the dynamical equivalents of the Lagrangian points. The usual tools to study numerically the RTBP are the combination of fixing suitable energy levels and suitable Poincaré sections. We note that, in time dependent models like the BCP and the QBCP the Hamiltonian is no longer preserved. A standard tool to deal with these periodic time-dependent systems is the so called stroboscopic map: Let U ⊂ R n be an open set, T the period of the vectorfield and ϕ : where ϕ(t 0 , t, x 0 ) stands for the solution which, at time t 0 lies at x 0 evaluated at time t, the flow of the differential equation. We define the stroboscopic map for x ∈ U as f (x) = ϕ(0, T, x). In this work we care about Hamiltonian differential equations. In this case, the stroboscopic map is symplectic.

Invariant Structures:
The simplest invariant objects of the original system are periodic orbits with the same period as the vectorfield. These appear as fixed point of the stroboscopic map. The monodromy matrix associated to these periodic orbits is the differential of the stroboscopic map. The eigenvalues of this matrix determine the linear behaviour around the fixed points and, under generic conditions, give some insight about the local non-lineal dynamics. In the symplectic case, under generic conditions, each elliptic direction give rise to a family of invariant curves which can be parametrized by the frequency ( [JV97b]). This frequency approaches to the normal mode responsible form the birth of the family at the fixed point. Along the hyperbolic directions, unstable (stable) manifolds depart (arrive) from the fixed points. These invariant objects are crucial to understand the dynamics of the system.
In this paper we focus on the information that can be extracted from the computation of invariant curves and high order approximations of unstable invariant manifolds. While it is quite common, in the literature related astrodynamics, to find discussions on the computation of invariant tori of maps [CJ00, JO04, Jor01] it is not so usual for high order approximations of stable/unstable manifolds.

High order approximation of unstable manifolds using the parametrization method
Let U ⊂ R n be an open set and assume that we are given a map f : U → R n induced by the evaluation at time T of a flow of some ordinary differential equation (stroboscopic map). The following discussion can be done for any Poincaré map as well. Here we assume that the section is temporal for simplicity of the exposition and because it is the natural section to chose in a periodically perturbed autonomous system. Let us suppose also thatx ∈ U is a fixed point i.e. f (x) =x. Obviouslyx is an initial condition for a T -periodic orbit of the original flow. The linearized normal behaviour around the fixed point is given by the eigenvalues of the differential of the map evaluated at the point. Assume spec Df = {λ, λ 2 , . . . , λ n } with |λ| > 1. Under generic conditions, we know that there exist a 1-dimensional unstable invariant manifold related to the fixed point. That is, there exist an open interval I ⊂ R and a map x : I → U such that x(0) =x and Equation (2) is known as the invariance equation of the invariant manifold. The parametrization method [CFdlL03a, CFdlL03b, CFdlL05, Sim90] is a powerful tool to, both, prove the existence of the manifold and compute it. The idea is to expand the parameterization of the manifold in Taylor series at s = 0 and solve equation (2) order by order. This makes sense in the case when both the map and the manifold are analytic. This assumption is fulfilled by the applications we are interested in. Hence, the goal is to compute a semi-analytic approximation of a parametrization of the invariant manifold. Let us name We are interested in numerically compute the coefficients a j for j = 0, . . . , d for a given degree d. This is achieved by an online scheme in which we solve equation (2) order by order: • Order 0 is given by the coordinates of the fixed point.
• Order 1 is given by the eigenvector related to λ.
• For k > 0, order k + 1 is given by the solution of the following linear system: Here, b k+1 is the k + 1-th term of the evaluation of the manifold up to degree k by the map f , that is: where the subindices (·) ≤k+1 denote the truncation of the power expansion of (·) at order k + 1.
Notice that it is mandatory to have a method to compose power expansions of the manifold with the map itself. Sometimes, when, the map is explicit, one is able to find a suitable recurrence expression to compute the terms of this composition. If a recurrence is not available one can compute higher order terms by automatic differentiation 1 . In the case we are interested, the map is not given explicitly but comes from a numerical integration of a differential equation. Here, the only reasonable strategy seems to use Jet Transport. This technique is based in the idea of, instead of integrating a single point, integrate a function given by its expansion in Taylor series. That is, one transports the jet, the set of derivatives of the function at a given point, up to a given order. It is straightforward to construct an integrator of jets from an integrator of numbers. It is only a matter of replacing the standard floating arithmetic by a polynomial arithmetic [FGJ + 18].
There is another obstacle that can appear when dealing with Poincaré maps: if the orbit is very unstable, the hyperbolic direction may lead to a huge error propagation. This problem can be avoided by using parallel shooting. The idea behind parallel shooting is to enlarge the dimension of the system in order to decrease the time of integration. Let us denote by ϕ t f t 0 (x) the solution of the differential equation with initial condition (t 0 , x) evaluated at time t f . Fix k ∈ N, the number of sections, and set h = T /k.
1 In principle, one could also use numerical differentiation but it is a bad approach in terms of efficiency and accuracy Figure 2: The Bicircular model.
1. y is a fixed point of F if and only ifx is a fixed point of f .
3. The projection to the first coordinate of the invariant manifold of F related to ζ coincides with the invariant manifold of f related to ζ k = λ.

The Bicircular Problem and the triangular points
The BCP is a perturbation of the RTBP. It is usual to take the units and the synodic coordinates of the Earth-Moon RTBP (see Figure 2). The BCP is not coherent, that is, the trajectories followed by the primaries do not obey Newton's laws. This is not an inconvenient since the model has been shown to be useful to describe the dynamics near the triangular points [SGJM95]. As a dynamical system, the BCP is a Hamiltonian system with three and a half degrees of freedom, i.e., a non-autonomous periodically time dependent with three degrees of freedom. The Hamiltonian function, written in the RTBP coordinates and units, is given by Here µ, r P E and r P M denote the same quantities as in (1). Moreover, m S denotes the mass of Sun, a S the averaged semi-major axis of Sun, θ = ω S t, ω S is the frequency of Sun in this system of reference, T S = 2π ω S is its period and finally, r 2 P S = (x − a S cos θ) 2 + (y − a S sin θ) 2 + z 2 . Notice that this Hamiltonian can be splitted in two parts: The time dependent part contains two terms, the Coriolis effect due to the rotating frame of coordinates and Sun's gravitational potential. The Taylor expansion of the potential starts with 1 a S 1 + x cos θ − y sin θ a S .
Therefore, the Hamiltonian, if we truncate the Sun's potential at linear order is written as So, the Coriolis term and the truncated Sun's potential cancel out and the dynamics is the one of the RTBP. This is to say that the contribution due to Sun's potential starts at order two, that is, the BCP is a periodic time dependent perturbation with size O( m S a 3 S ) ≈ 0.0056. Anyhow, it is large enough to produce remarkable changes on the dynamics, especially near the triangular points. In Figures 3 (left) and 6 (left) we display continuations from the RTBP to the BCP. The vertical axis in these plots represent an artificial parameter ε which multiplies the mass of Sun. Therefore, when ε = 0 the model corresponds to the RTBP and, when ε = 1 the model corresponds to the BCP. We shall comment these figures in detail in the next sections.

Dynamical equivalents of the Triangular points
First of all let us mention that, due to a symmetry, the dynamics near L 4 is the same as the dynamics near L 5 (in fact, this symmetry maps orbits in the region y > 0 to orbits in the region y < 0). A feature of the BCP to be stressed is that the region around the geometrically defined triangular points is unstable. The influence of Sun's potential is enough to produce a bifurcation in the periodic orbit that replaces L 4 (i.e., L 5 ). It is well known [SGJM95] that each triangular point is replaced by three periodic orbits with the same period as Sun. One small and unstable (the actual replacement of L 4 ) and two which are stable. We have named these orbits P O1, P O2 and P O3. See Figure 3 (left) a continuation diagram from the RTBP to the BCP. The two additional periodic orbits are produced by an imperfect pitchfork bifurcation (i.e., a pitchfork bifurcation broken due to a loss of symmetry). One may ask which is the model that displays the perfect bifurcation and which is the broken symmetry. To address this question we take a look at the order two of the Taylor expansion of Sun's gravitational potential. We have We have named T (x, y, θ) = −x cos θ+y sin θ. We would like to stress again that H 2 S is the first contributing non-autonomous term in the model due to the cancellation produced by the Coriolis acceleration. This term is invariant under the symmetry (x, y, x, θ) → (x, −y, x, −θ). The order three of the expansion is given by: Here ρ 2 = x 2 + y 2 . The polynomial in T is no longer even. This breaks the symmetry and, hence, the pitchfork bifurcation. The non-autonomous model that displays the perfect bifurcation is: Figure 3 (left, curve in red) shows the continuation diagram from the RTBP to this simplified version of the BCP. Due to the symmetry, periodic orbits P O2 and P O3 only differ on the phase on the orbit.

Phase space of the stroboscopic map near the triangular points
The three periodic orbits appear as fixed points of the stroboscopic map. We recall that their linear normal behaviour is of type saddle×centre×centre for P O1 and totally elliptic for P O2 and P O3. There are several ways to justify that, from the elliptic directions of each fixed point, there is a family of invariant curves whose frequency tends to the normal modes of the fixed points [JV97b,JV97a]. Therefore, we have a family of invariant curves for each elliptic direction, that is, two for P O1 (HF 1 in the horizonal plane and V F 1 in the vertical direction), three for P O2 (HF 2F 1 and HF 2F 2 are horizontal, and V F 2 is vertical) and three for P O3 (HF 3F 1 and HF 3F 2 horizontal, and V F 3 vertical). The remaining eigendirection of P O1 is hyperbolic. There exist stable and unstable one-dimensional invariant manifolds associated to these hyperbolic directions. Initial conditions near the triangular points shadow the unstable manifold which wonder some time around the periodic orbits P O2 and P O3 and finally abandon the vicinity of the triangular points. These manifolds are of special interest if one plans to put or take out objects near L 4 and L 5 . The stable and unstable manifolds related to P O1 can be computed up to high order directly on the stroboscopic map (Section 3.1 and [FGJ + 18]). In Figure 5 we observe a projection of the phase portrait of the map. The three points displayed with crosses correspond to P O1 (in the middle), P O2 and P O3. It is displayed as well, semianalytical approximations of  the stable and unstable invariant manifolds. We have used an approximation of order 64. The width curve are the pieces given by the parameterization. The thin curve correspond to some iterations of these pieces. It can be observed, also, some invariant curves growing from P O2 and P O3. These invariant curves are totally elliptic near the fixed points.

Regions of effective stability near the triangular points
As we have observed, the triangular points are replaced by three periodic orbits, one of them unstable. This is the reason why the BCP is an interesting model [GJMS91c]. Indeed, the unstable manifold of the triangular periodic orbit takes initial conditions away from the vicinity of the triangular points. However, we can pursue on the seek for regions of (effective) stability out of the plane of motion of the primaries. The mechanism that suggest the existence of regions of effective stability is the stickiness of normally elliptic low dimensional invariant tori. See [JV98,JV97a]. As we discussed before, there are families of invariant tori emanating from the periodic orbits P O2 and P O3. These families are elliptic close enough to the fixed points. This results on two small regions of effective stability in the plane of motion of the primaries related to the totally elliptic orbits.
We put our attention on the vertical families (one for each orbit) of invariant tori. In Figure 3, Right, we display how these families vary when they grow out of the plane of motion of the primaries. We observe that the three families display a broken pitchfork symmetry, analogous to the one of the periodic orbits. The linear normal behaviour of the tori is the same as the periodic orbit near the plane. As a consequence of the pitchfork bifurcation, at some some distance of the plane, the surviving family is totally elliptic. Therefore, the tori are sticky and regions of effective stability are to be expected. We label these families V F 1, V F 2 and V F 3 after the corresponding periodic orbits. The families V F 1 and V F 2 are connected, as it happens with P O1 and P O2. On the other hand V F 3 reaches high amplitudes in the (z, p z ) plane. It is known that, skipping resonances, the three families have the same stability as the corresponding periodic orbits. Therefore the tori of V F 1 have hyperbolic directions while the ones of V F 2 are normally elliptic (except for small intervals of instability produced by resonances involving internal and normal frequencies). Recall that both families are connected and the change of stability takes place at a turning point. The tori of V F 3 are normally elliptic up to very high values (again, except for resonances).
Normally elliptic lower dimensional tori induce regions of effective stability. Numerical estimations of the shape and the size of these regions show that, in the case of V F 2, the regions are small and narrow while in the case of V F 3 large regions exist for sufficiently high values of the vertical amplitude. In Figure 5 we show two stability regions out of the horizontal plane. These regions seem to persist in the real model for time spans of 1000 years [Jor00,HXSW15]. The effect of Solar Radiation Pressure on the effective stability regions of the BCP is discussed in [JCFJ15].  Let us explain how Figure 5 is obtained. Each of the vertical family of invariant tori can be identified by the value of the coordinate p z when z = 0 and p z > 0. Denote by a(p z ) ∈ R 6 the coordinates of the point that identifies a torus. We have to select a set of initial conditions near a(p z ) and integrate them for a long time span. Let us be more precise on how to select the initial conditions. We use a two dimensional grid, the coordinates z, p x , p y and p z will be fixed by the corresponding values of a(p z ). To adapt to the shape of the regions, we use a polar-like grid, centered at Earth: where h r and h α are used to control the density of the grid. The computation goes as follows. Take a point of the grid and integrate the vector field 15000 Moon revolutions. At each integration step, we test if there is a collision with Earth or Moon. If there is a collision, or the coordinate y becomes negative, we stop the integration (Recall that we are interested in the points that remain close to L 4 ). We have used h r = 0.001 and h α = 0.0002. The difference on the sizes of these small quantities is aimed to produce a nearly squared grid.

A weakness of the BCP
The translunar point is one of the most interesting locations at the Earth-Moon system. The reason is that L 2 seems suitable to observe the far side at Moon. Taking into account that, a natural criticism to the BCP is that it does not have a dynamical replacement of L 2 . In Figure 6 (left) it is displayed a continuation of L 2 as a periodic orbit from the RTBP to the BCP. The periodic orbits are identified by their coordinates as fixed points of the Stroboscopic map. These orbits have been computed by means of the parallel shooting method. Again, the vertical axis is an additional parameter, ε, multiplying the mass of Sun. The point L 2 is the middle crossing of the characteristic curve with the homotopy level corresponding to the RTBP, at the bottom. The other two points of the RTBP correspond to the same 1 : 2 resonant planar Lyapunov orbit but with different initial times. We observe that the continuation of L 2 reaches a turning point and it never reaches the homotopy level of the BCP. The result is that the translunar dynamical structure is lost in the BCP. This suggest that a more complex model needs to be used to analyze the natural behaviour near the translunar point. The resonant Lyapunov orbit can be continued to the BCP. The result is a large orbit that remains away from the translunar point, see Figure 6 (right).

The Quasi-Bicircular Problem and the collinear points
The quasi-bicircular solution of the Earth-Moon-Sun system is planar i.e. the three bodies move in the same plane. After the quasi-bicircular solution is computed one can write the equations of motion of the test particle, prescribing the quasi-bicircular solution as motion for the primaries. It is usual to compute the quasi-bicircular solution in the Jacobi frame, however, if one has the purpose of describing the dynamics in the Earth-Moon vicinity, it is suitable to use the frame of coordinates corresponding to the Earth-Moon RTBP.
To do so, one has to perform three different transformations. First, one has to use a translation to move the origin from the global barycentre to Earth's and Moon's centre of masses. Second, one has to use a rotating (synodic) frame to keep Earth and Moon fixed on the horizontal axis. Third, the unit of length is scaled so the distance between Earth The resulting model is a Hamiltonian system with three and a half degrees of freedom. The Hamiltonian function can be written as H = 1 2 α 1 (p 2 x +p 2 y +p 2 z )+α 2 (p x x+p y y+p z z)+α 3 (p x y−p y x)+α 4 x+α 5 y−α 6 1 − µ r pe + µ r pm + m S r ps , (4) where, r 2 pe = (x − µ) 2 + y 2 + z 2 , r 2 pm = (x − µ + 1) 2 + y 2 + z 2 , r 2 ps = (x − α 7 ) 2 + (y − α 8 ) 2 + z 2 , and for i = 1, . . . , 8 α i : T → R are periodic functions. That is, Here, θ = ω S t and ω S is the frequency of Sun. Moreover, α i is odd for i = 1, 3, 4, 6, 7 and even for i = 2, 5, 8. Obviously one can only have a numerical approximation of these functions. In this case, we take advantage on the computations done in [And98] and take the same values for the Fourier coefficients of the periodic functions α i 's. To end, and taking into account the properties of the functions α i 's, it is easy to see that the Hamiltonian function (4) has the symmetry (θ, x, y, z,ẋ,ẏ,ż) → (−θ, x, −y, z, −ẋ,ẏ, −ż), The meaning of these periodic functions is the following: 1. (α 7 , α 8 , 0) is the position of Sun in the plane of motion of the primaries.
2. α 1 , α 2 , α 3 and α 6 capture the fact that the distance between Earth and Moon is not constant.
3. α 4 and α 5 take into account the Coriolis effect due to the rotating frame of reference.

Dynamical equivalents of the collinear points
In this section we give some words about the minimal periodic orbits that replace the collinear points in the QBCP. In Figure 7 we display the dynamical equivalents, from left to right, of L 1 , L 2 and L 3 . We observe that the orbits replacing L 1 and L 2 are small, their maximal distance to the corresponding equilibrium point is of order O(10 −6 ). As the original equilibrium points, the linear normal behaviour of these orbits is of type saddle×centre×centre. In Table 1 we display the eigenvalues of each orbit. We notice that the unstable direction of L 1 (of order 10 8 ) and the unstable direction of L 2 (of order 10 6 ) are large and this implies huge propagation of error near these orbits. On the other hand, the dynamical equivalent of L 3 has a very weak unstable direction, at least compared to the other two.  Table 1: Eigenvalues of the three dynamical equivalents of L 1 , L 2 and L 3 . We only put three for each orbit. The rest are given by their inverses due to the symplectic character of the stroboscopic map.

Resonant orbits of low order
As the QBCP is a T S -periodic system, the simplest invariant objects are T S -periodic orbits. We already have mentioned that the equilibrium points are replaced by these periodic orbits of minimal period. Periodic orbits of the RTBP whose frequency is resonant with the one of Sun also persist as T S -periodic orbits in the QBCP. The Lyapunov and Halo families of periodic orbits related to the equilibrium points L 1 and L 2 , are a source for these kind of resonant orbits. In contrast with the families related with L 3 , the families of the two first libration points are nourished with low order resonant orbits. A relation of low order resonant periodic orbits of the RTBP can be found in [And98].
In [GM01] the authors show the ranges for the admissible periods for each family. The families related to L 3 are of relatively large period and there are not many periodic orbits whose frequency are in low order rational relation with the frequency of Sun. There is, however, a 1 : 1 resonant periodic orbit near the end of the vertical family. This orbit is enormous in size and cannot be considered in the vicinity of L 3 . In Table 2 details of the continuations of low order resonant orbits from the RTBP to the QBCP are given: The first column corresponds to resonant periodic orbits of the RTBP. The label in this first column consist in three numbers that encode each orbit. The first is a zero and indicates that the orbit belongs to the RTBP (this is intended to distinguish them from the orbits in the last column corresponding to the QBCP). The second number refers to the libration point related to each orbit (all of them belong to Lyapunov and Halo families related to L 1 and L 2 ). The third number is just an enumeration. The second column indicates the order of the resonance. We stress that the influence of Sun is relevant enough to produce bifurcating orbits in each of the continuations. The third column shows how many orbits bifurcate from the original ones when they are continued to the QBCP. Finally the last column contains the labels of the resulting orbits in the QBCP. Table 2 can be found originally in [And98]. We have added the order of the resonance and the color code to indicate the linear normal behaviour of each orbit. Labels in blue stand for orbits of type saddle×centre×centre. Labels in green denote linear character of the kind saddle×saddle×centre. Names in cyan denote totally hyperbolic orbits. The color yellow denotes totally elliptic orbits. The continuation for the orbits in red do not reach the homotopy level of the QBCP and, therefore, are not considered.

High order approximation of the unstable manifolds of the collinear periodic orbits
This section is devoted to the results of implementing the algorithm explained in Section 3.1 to the dynamical equivalents of the collinear points. Figure 8 shows pieces of the stable (dashed) and unstable (solid) manifolds related to the three collinear periodic orbits (the other branches can be obtained by symmetry). From left to right: the one related to L 1 , the one related to L 2 and the one related to L 3 . We would like to remark that these pieces are obtained directly from the evaluation of the approximation (of order 64) of the manifolds. The error is controlled by checking that the contribution of the last term of the approximating polynomial is small. It is also checked that all the points  in the pieces verify the invariance equation with hight accuracy. We can observe that these approximations already give large excursions far away from the collinear points. Especially in the case of L 3 , where the piece of the manifold passes very close to the triangular points. The axes of Figure 8 show the x and y values. These pieces can be mapped through the stroboscopic map to obtain larger pieces of the manifolds if it is necessary. The point of giving high order approximations of the manifold is that, just a fewer number of iterates are necessary. For the computation of the manifold related to L 3 , a simple shooting method has been used. Indeed, the instability associated to this libration orbit is very weak. For the computation of the manifold related to L 2 , multiple shooting is required. We have used two sections. For the computation of the manifold related to L 1 , the most unstable one, we have used a single shooting strategy but with an extended precision arithmetic of 128 bits. This last approach makes the program far slower but very simple to code. In Figure 9 (left) we show the resonant orbit 2G of Table  2. We display also (right) the stable (dashed) and unstable (manifolds).

Conclusions
We have presented two alternatives to the RTBP for the study of the motion of a test particle in the Earth-Moon system. Both models, the BCP and the QBCP, depend periodically on the time. We use the so-called stroboscopic map to study the minimal periodic orbits of the systems and the invariant manifolds related to them. The BCP is as useful model for the study of the triangular points. The simplicity of the vectorfield is a strong point, especially in problems related to effective stability where massive integrations are mandatory. We have also stressed its weakness: it is not suitable to understand the dynamics around the collinear points. The BCP is useless to describe the vicinity of the translunar point. We have used the parametrization method to obtain high order approximations of the unstable manifolds related to the minimal periodic orbits that replace the collinear points in the QBCP. This is helpful to design long excursions between the two primaries and the collinear points. The main novelty is that we have computed the manifolds directly on the stroboscopic map. The QBCP is a complicated model with a numerically computed vectorfield. This makes it a bad candidate (in front of the BCP) to be the model used to face the problems involving massive simulations related to the triangular points.
We would like to stress that the BCP should be used to face problems related to the triangular points. Especially if this problems involve large time integrations to seek for regions of practical stability. The QBCP should be used when dealing with problems involving the collinear points.

Technical details
All the computations appearing in the Figures of this paper, also the ones which appear in the literature, have been performed by the authors. The integrations for the RTBP, the BCP and the QBCP have used a Taylor method with variable order and stepsize. The demanded accuracy for the standard double precision has been 10 −16 . The computations in multiple accuracy have been done using the library mpfr. The LAPACK library has also been used for some computations related to linear algebra. The rest of the programs have been written by the authors in C and C++ languages from the scratch. Table 3 -5.00e-02 0.