Families of Halo-like invariant tori around L 2 in the Earth-Moon Bicircular Problem

The Bicircular Problem (BCP) is a periodic time dependent perturbation of the Earth-Moon Restricted Three-Body Problem that includes the direct gravitational effect of the Sun. In this paper we use the BCP to study the existence of Halo-like orbits around L 2 in the Earth-Moon system taking into account the perturbation of the Sun. By means of computing families of 2D invariant tori, we show that there are at least two different families of Halo-like quasi-periodic orbits around L 2 .


Introduction
In recent years, major space agencies have shown interest in concepts that involve using the Moon and its neighboring area as candidates to host space assets to support scientific missions or commercial endeavors. An example is the Lunar Gateway, a permanent space This article is part of the topical collection on Toward the Moon and Beyond. Guest Editors: Terry Alfriend, Pini Gurfil and Ryan P. Russell. This work has been supported by the Spanish grant PGC2018-100699-B-I00 (MCIU/AEI/FEDER, UE) and the Catalan Grant 2017 SGR 1374. The project leading to this application has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant Agreement #734557. The authors thank the comments by the reviwers that helped to improve this manuscript. As we have mentioned, in this paper we explore the neighborhood of the translunar point in the BCP. In particular we are interested in the counterparts (in the BCP) of the well-known Halo families (see Breakwell and Brown 1979). In this model, the Halo families are no longer composed by periodic orbits but by quasi-periodic orbits with two basic frequencies, one coming from the Halo orbits of the RTBP, and another from the frequency of the Sun (ω S ). To compute these families we use a combination of a method to approximate invariant curves with multiple shooting, with the continuation method to generate a complete atlas of the dynamical equivalents of the Lyapunov and Halo families near the translunar point. Notice that, due the absence of a natural replacement of L 2 , the properties of some of these families change near the coordinates of the translunar point (which is no longer an equilibrium point in the BCP). In particular, we report the existence of a family of Halo-like orbits that does not come from the original Halo family in the RTBP.
This paper is structured as follows: Sect. 2 contains an analysis of the effect that the Sun, as modeled in the BCP, has on the L 2 point. This serves as a motivation for Sect. 3, where after a brief discussion on two approaches to study the dynamics around the L 2 point, the method of tori continuation is justified as appropriate for the L 2 region and explained. Section 3 also includes the strategy employed to find the different families. Section 4 elaborates on the results obtained from tori continuation, focusing on the Halo-like tori and their stability. The focus in Halo-like orbits is not arbitrary, and it responds to the application these trajectories have for lunar missions. This is also discussed in Sect. 4. Finally, Sect. 5 presents the conclusions and further work. We have added "Appendix 1", an extra section that contains figures of Lyapunov-like orbits. This is done for completeness and to remark the difference between Halo-like and other trajectories.

From RTBP to BCP: the L point case
The most usual formulation of the RTBP is the circular version. In this model, it is assumed that the Earth (E) and the Moon (M) revolve along a circular orbit centered in their common barycentre (B). It is standard to consider a synodic reference frame, that is, a rotating frame that fixes the primaries at the horizontal axis.
In the BCP, the dynamics of the Earth, Moon and Sun (S) are simplified considering that the three bodies orbit in the same plane. Also, it is considered that the Earth and the Moon follow a circular orbit around their barycenter (as in the RTBP), and that B is orbiting around the S-E/M barycenter. Note that this model is not coherent, in the sense that the motion of the three massive bodies is not described by the Newton's equations of motion.
The BCP is usually written as a periodically time-dependent perturbation of the RTBP. In our case, the RTBP is the Earth-Moon system and the perturbing body is the Sun. It is then natural to use the units and reference frame of the Earth-Moon RTBP, so that the Sun is moving around in a circular orbit. For more details see Chapter 3 of Gómez et al. (1993). As in the RTBP, if we consider the momenta p x =ẋ − y, p y =ẏ + x, p z =ż, the BCP admits a Hamiltonian formulation as follows: where r 2 P E = (x − μ) 2 + y 2 + z 2 , r 2 P M = (x − μ + 1) 2 + y 2 + z 2 , r 2 P S = (x − x S ) 2 + (y − y S ) 2 + z 2 , x S = a S cos ϑ, y S = −a S sin ϑ, and ϑ = ω S t with ω S being the frequency of the   Table 1. Note that in this reference system the Sun moves around the origin in a circular motion (see Fig. 1). For the details on the derivation of the BCP equations of motion, the interested reader is referred to Gómez et al. (1993).
An important observation is that the Hamiltonian (1) depends periodically on time. This periodic effect captures the direct gravitational influence of the Sun. Moreover, Hamiltonian (1) can be expressed as a time-periodic perturbation of the RTBP, where X = (x, y, z), P X = ( p x , p y , p z ). Notice that the autonomous part is the Hamiltonian of the RTBP, and H S is the Hamiltonian associated to the perturbation due to the Sun. The consequences on the dynamics of this periodic time-dependency will be explained later. It is well known that L 2 is an equilibrium point of the RTBP. However, as opposed to the RTBP, the BCP is not an autonomous system, and it depends periodically on time so that the L 2 point is not an equilibrium point anymore. In a general setting (not necessarily Hamiltonian), if a periodic perturbation is applied to a differential equation then, under generic conditions of non-degeneracy, an equilibrium point becomes a periodic orbit with the same period as the perturbation. Applying this principle, in this section we explain what are the dynamical consequences that the time-periodic perturbation has on the L 2 point in the context of the BCP.
The approach taken to study the transition of the L 2 point from the RTBP to the BCP is by using a pseudo-arclength continuation scheme with respect to an artificial parameter used to with |ε| ≤ 1. Note that for ε = 0, H 0 = H RTBP , and for ε = 1, H 1 = H BCP . Considering ε = 0, the five Lagrange points (L i , i = 1, . . . , 5) are equilibrium points of the system (2).
When |ε| is small enough, the equilibrium points become periodic orbits around the point L i (now defined only geometrically since they are no longer equilibrium points) with the same period as the perturbation. In our case, the period is equal to the period of the Sun, T = 2π/ω S . More generally, a periodic orbit of the RTBP whose period is p q T , becomes a periodic orbit of period pT once ε is set to be different from zero. This is a consequence of the Implicit Function Theorem.
However, the perturbation due to the gravitational potential of the Sun in the BCP cannot be considered small, i.e., is large enough to produce bifurcations. For example, around the triangular points L 4 and L 5 there is a loss of uniqueness of the periodic orbit, and three periodic orbits appear (see Gómez et al. 1987;Simó et al. 1995). The size of the perturbation also affects the L 2 point.
Besides showing the existence of these periodic orbits, computing their stability is essential to have the full picture. By means of analyzing the spectra of the monodromy matrix, we also show how the linear stability of these new periodic orbits evolves with respect to the continuation parameter ε. To continue numerically the periodic orbits, we use the stroboscopic map, this is, the map obtained from evaluating the flow at the period T . The coordinates of a T -periodic orbit of the system provide fixed points of the stroboscopic map. Moreover, the differential of the map, coincides with the monodromy matrix.
It is important to note that due to the highly unstable nature of the L 2 region, the algorithm to compute periodic orbits has to be implemented using a multiple shooting scheme. This is a pretty standard procedure (see Stoer and Bulirsch 2002;Seydel 2009) and the details can be found in Gómez and Mondelo (2001) for the RTBP. For the present work, the total number of sections used is four. The results of continuing the L 2 point with respect to ε are shown in Fig. 2. The horizontal axis is the x component of the periodic orbit at t = 0, and the vertical axis is the parameter ε applied to the mass of the Sun. Starting form L 2 , and moving to the left the parameter increases until it hits a local maximum, and then decreases to cross the horizontal line and become negative. The point on ε = 0 corresponds to a planar Lyapunov orbit whose period is half the one of the Sun, so we can see it as a closed trajectory traveling twice around L 2 in a single period of the Sun (i.e., a 1:2 resonant planar Lyapunov orbit). Moving from L 2 to the right, although the parameter ε becomes negative (which has no physical sense), it decreases until it hits a turning point, and then increases to become positive and reach ε = 1, that is, the BCP. Again, in this case, the crossing point with ε = 0 corresponds to the previous 1:2 resonant planar Lyapunov orbit. As it has been mentioned before, for ε = 0, the 1:2 resonant orbit consists on a single loop traveling twice around L 2 in one period of time. When ε > 0, this 1:2 resonant orbit becomes a periodic orbit that nearly travels the same trajectory twice before closing the loop. This behavior is maintained until ε = 1, see Fig. 3, where the L 2 point and the Moon are added for reference.
As a summary of the previous discussion, the value of ε is not small enough and there is no natural dynamical substitute of the L 2 point in the BCP. In other words, there is no direct connection between L 2 and a periodic orbit in the BCP. Figure 2 also contains the details on the linear stability of the periodic orbits computed. As a technical remark, it is not a good idea to multiply the differentials of the flow between consecutive points since this produces a matrix with very large entries, and therefore loosing accuracy in the determination of the eigenvalues of small modulus (like the stable and elliptic directions). So, it is better to use a specific technique (see, for instance Gonzalez and Mireles James 2016). This particular approach is the natural way to proceed when dealing with multiple shooting: the spectrum of the orbit is recovered from the spectrum of the Jacobian matrix related to the fixed point map associated to the multiple shooting and it is suitable to compute (stable and unstable) invariant manifolds related to highly hyperbolic periodic orbits and invariant tori. We do not provide more details as the behavior of the stable and unstable manifolds of the orbits are out of the scope of the paper.
It is observed that the periodic orbits alternate between the types saddle×center×center (green regions) and saddle×saddle×center (red regions). Starting from L 2 , the linear stability Due to the Hamiltonian nature of the system, the other three eigenvalues are λ −1 i , i = 1, 2, 3. Also, note that due to the non-autonomous character of the BCP Hamiltonian, there is no double eigenvalue 1 is of the type saddle×center×center. Moving to the left, ε increases and the periodic orbits keep this linear stability type until they hit the local maximum. In this turning point, the linear stability becomes of the type saddle×saddle×center until another bifurcation point at resonant 1:2 planar Lyapunov (ε = 0). At this point, ε becomes negative, and the resulting periodic orbits are of the type saddle×center×center. A similar pattern but with different sign for ε is observed when moving to the right of the L 2 point. In this scenario, ε decreases and maintains the same linear stability type as L 2 until they hit a local minimum. As before, in this turning point, the linear stability becomes of the type saddle×saddle×center until ε = 0, where there is yet another bifurcation. At this bifurcation, the linear stability becomes of the type saddle×center×center. Finally, this resonant planar Lyapunov orbit is continued until the last bifurcation point. This is a pitchfork bifurcation, and it is where the 1:2 resonant (with the Sun) Halo orbit in the RTBP ends (this is shown in Andreu (1998), where a bifurcation diagram and a suitable analysis are provided ). The implication is that the 1:2 resonant Halo orbit in the RTBP does not reach the BCP. As we will see in Sect. 4, this is not the case for all Halo orbits, and there is a dense set of Halo orbits that survive the perturbation of the Sun as modeled in the BCP.
After this point, the stability of the periodic orbits is of the type saddle×saddle×center until ε = 1. The eigenvalues λ i , i = 1, . . . 6 of the monodromy matrix associated to the periodic orbit in the BCP are captured in Table 2. Notice that, the final orbit has not the same stability character as L 2 . It is important to note that the nature of the perturbation shapes the dynamics around an equilibrium point. Moreover, as it is seen in Fig. 3, the size of the final orbit is much larger than the one expected for a dynamical equivalent of L 2 . The comparison between the BCP and QBCP illustrates this phenomena. In the QBCP, L 2 is replaced by a periodic orbit that is small in the sense that its maximal distance to L 2 is of the order of 10 −6 , and it has the same stability type of the L 2 point. See (Andreu 2002;Jorba-Cuscó et al. 2018) and references therein for the details.

Approach to study the vicinity of L 2
To study the dynamics of a dynamical system, a typical approach is to look for invariant objects and analyze their stability. Typically, the analysis starts looking for equilibrium points, then periodic orbits, 2D tori, and so on. Section 2 covered the analysis of the periodic orbit with the same period as the effect of the Sun. However, this does not provide the full picture and gives little insight on the dynamics around the L 2 point.
One approach to get the full picture of the dynamics is to do a reduction to the center manifold of the periodic orbit. This approach consists in a series of changes of variable to decouple the saddles from the centers. This decoupling allows to reduce the dimension of the system, and to focus only on the invariant objects that live in the center manifold. This technique has been proven very successful to characterize the dynamics around the collinear points in the RTBP (Jorba and Masdemont 1999) for different mass parameters; around the L 1 point the BCP  and L 2 in the QBCP (Andreu 2002;Le Bihan et al. 2017); or around the L 1 and L 2 points in the Sun-Earth RTBP for solar sails (Farrés and Jorba 2010). Note that the systems that can be studied with this technique are very broad: the reference (Jorba and Masdemont 1999) deals with autonomous Hamiltonians, the references Andreu 2002;Le Bihan et al. 2017) with Hamiltonians that depend periodically on the time, and the reference (Farrés and Jorba 2010) with general Ordinary Differential Equations. The interested reader is referred to Carr (1981); Sijbrand (1985); Vanderbauwhede (1989) for a more general treatise on the center manifold and its applications. The main advantage of this method is that it provides a comprehensive picture of orbits staying in a neighborhood of an invariant object and its bifurcations. The disadvantages are that, due to the construction of the center manifold, the neighborhood where it is valid may be very small due to the presence of small divisors.
An alternative to the center manifold is to directly compute the families of invariant objects that shape the phase space of the dynamical system (equilibrium points, periodic orbits, 2D tori and so on). A key advantage of this approach is that it can be applied far away from L 2 . Also, in some cases (equilibrium points, periodic orbits, and 2D tori) there are techniques to compute the stability of each member of the family. The main limitation of this approach is that computing tori of dimension higher that 2 is very expensive computationally (Jorba and Olmedo 2009) and, sometimes, cumbersome. Examples in the context of the BCP can be found in Castellà (2003), where families of 3D tori around the triangular points were computed.
In addition to that, the continuation of these objects involves some level of trial and error, and once the continuation process starts, a lot of fine tuning due to the presence of resonances is needed. Finally, and as opposed to the center manifold approach, this method provides an incomplete picture unless all relevant invariant objects are computed.
Note that the latter approach assumes the existence of families of invariant objects. The assumption deserves some explanation. In the context of the BCP, the existence of invariant tori is inherited from the RTBP. It is well know that around the collinear equilibrium points of the RTBP there are families of periodic orbits (planar and vertical Lyapunov, and Halo orbits) and quasi-periodic orbits (quasi-halos and Lissajous). See (Jorba and Masdemont 1999;Gómez and Mondelo 2001) for details. Under generic conditions of non-resonance and non-degeneracy, adding a small enough periodic (or quasi-periodic) time-dependent perturbation to RTBP, causes the existing invariant objects to inherit the frequencies of the perturbation. It is important to mention that the families of invariant objects become Cantorian because only those frequencies satisfying a suitable non-resonance condition survive. As a consequence, the families of objects are Cantorian, not continuous. The details on the proofs that back these statements can be found in Jorba and Villanueva (1997). Finally, an example of this phenomena in the context of the RTBP and the BCP can be found in .
The study of the L 2 region in the BCP was initially approached using the reduction to the center manifold. Actually, the code used to generate the results in  for L 1 in the BCP, was initially developed to study the neighborhood of L 2 . However, the radius of convergence of the computed center manifold was very small, and it was concluded that this approach was not suitable for L 2 . Hence, it was decided to compute families of 2D tori around the L 2 region along with their stability following the methods described (Castellà and Jorba 2000;Jorba 2001), respectively. The next subsections outline the numerical methods to compute invariant tori, their stability, and the continuation strategy.

Computation of highly unstable invariant tori
The method used to compute invariant tori is based on Castellà and Jorba (2000). The general statement of the problem is the following: assume there exists a quasi-periodic orbit with two basic frequencies ω 1 , ω 2 ∈ R such that ω 1 /ω 2 ∈ R Q. This means that there exists a map : T 2 → R 6 (the parameterization of the torus) such that the function W : R → R 6 defined by W (t) = (ω 1 t, ω 2 t) is a trajectory of the system.
In the scenario of the BCP, one of the frequencies is equal to the frequency of the Sun (ω S ) so, from now on, ω 2 = ω S . Now, let us define the stroboscopic map F as the flow of the BCP, φ BCP , at time T = 2π/ω S . Note that now the closed curve θ ∈ T → W (θ, 0) ∈ R 6 is invariant by F, Thus, knowing that one of the fundamental frequencies of the motion is ω S , the problem of computing a torus is reduced to finding a functionŴ : T 1 → R n that satisfies Eq. (3) for a given ρ (note that to know ρ is equivalent to know ω 1 ). Such functionŴ is called an invariant curve with rotation number ρ. Obviously, ifX is an invariant curve with rotation number ρ, it satisfies that From a practical point of view, the approach is to find a zero of G. A convenient way to approximate an invariant curve is to use its (truncated) Fourier series, Hence, the goal is to compute the Fourier coefficients a i , b i , i = 0, . . . , N such that they define a periodic function X which is a zero of (4). This leads to (2N +1)n unknowns. Hence, at least the same number of equations is required to solve for all a i , b i , i = 0, . . . , N . To this end, (4) is discretized by using an equispaced grid of values of θ such that This provides the number of equations needed to solve for the Fourier coefficients Finally, an extra equation specifying a value for the Fourier coefficients at θ = 0 is required to resolve the ambiguity in the Fourier coefficients due to the fact that the map F is autonomous (see Castellà and Jorba 2000 for further details). This system of equations is solved by means of a standard Newton's method using least squares to account for the fact that we have more equations than unknowns.
In the same fashion as for the continuation of periodic orbits described in Sect. 2, the use of multiple shooting is required to mitigate the error growth due to the instability of the L 2 region (see Duarte 2020 for a discussion for the Sun-Jupiter L 1,2 ). We recall from Table 2 that the largest eigenvalue of the monodromy matrix of the periodic orbit around L 2 found in the BCP is of order of 10 6 . The following paragraphs illustrate how this is approached. Let us start with the following definition.
Definition 1 Let g 1 , . . . , g r diffeomorphisms of some subset of R n into itself, let W be the parametric representation of a closed curve of R n , θ ∈ T and let ρ ∈ T. Then, W is called an r-invariant curve for g 1 , . . . , g r with rotation number ρ if Remark 1 It is easy to check that if W is an r-invariant curve then, for any α ∈ R, W (θ + α) is also a r -invariant curve. This implies that there are different sets of Fourier coefficients representing the same r -invariant curve.
Given a r -invariant curve W 0 approximated by a truncated Fourier series (5), the goal is to compute its (2N + 1)n coefficients a i , b i , i = 0, . . . , N . The invariance condition for a r -invariant curve reads As a result, to find W 0 , it is also required to solve for W i , i = 1, . . . , r − 1. This is, there are a total of (2N + 1)nr unknowns corresponding to all the r -invariant curves. Now, we use the grid (6) to discretized each of the equations in (7), and the following set of equations is obtained, An extra equation specifying, for instance, a value of a coordinate at θ = 0 is required to resolve the ambiguity in the Fourier coefficients (see Remark 1). The system of Eqs. (8) is solved by means of a standard iterative Newton's method using least squares to account for the extra equation. The iteration process is stopped when the norm of the function becomes smaller than a prescribed tolerance (typically, a value of the order of 10 −6 is good enough for plots, but for the computation of the stability we have used 10 −10 ). Note that this method ends up computing r curves. This multiple shooting approach is useful to compute invariant curves for very unstable systems. In the case of interest, the L 2 region in the BCP, the maps g j , j = 1, . . . , r are defined as follows: if p denotes a point in the phase space, then where φ( p; t 1 , t 2 ) denotes the flow from time t 1 to time t 2 , and we recall that T is the period of the Sun. In this work, we use r = 4. Note that the convergence of the Newton's method does not guarantee that the solution is a good representation of the torus. Remember that we have computed the torus based upon a truncated Fourier series (5). To estimate the error of the actual representation, the invariance condition is checked on a finer mesh. If the error in the verification of this condition is larger than a prescribed threshold, then more Fourier coefficients are added in the representation (5), and the process starts again (see Castellà and Jorba 2000). Hence, the value of N is not fixed throughout the continuation process, and the estimate of the error of the representation is the main driver on how it evolves. Besides the error of the representation, the shape of the invariant curve is also a factor on how many Fourier nodes are needed in the representation (5) to capture the features of the curve. Typically, invariant curves far from a circle require a higher value of N . For this work, low values of N = 5 are typically used to start the process, with a maximum of N = 252.

Linear stability
To compute the stability of an invariant object is as important as the invariant object itself. The methods in this section are based on the results in Jorba (2001), that here we have adapted to a multiple shooting scheme. The following paragraphs provide an overview of the method to compute the stability of invariant curves, and the modification to work with unstable systems.
Let us assume that W is an invariant curve satisfying condition (3). To study the dynamical behavior close to the curve, we consider a small displacement h ∈ R n with respect to W . Then, Hence, using that F(W (θ )) = W (θ + ρ) and discarding the second order term, we have that the following dynamical system describes the linear normal behavior around the invariant curve,h where A(θ ) = D x F(W (θ )) and h ∈ R n . Let C(T 1 , C n ) be the set of continuous functions between T 1 and C n . If ψ ∈ C(T 1 , R n ), we define the operator T ρ : C(T 1 , C n ) → C(T 1 , C n ) as T ρ (ψ(θ )) = ψ(θ + ρ), θ ∈ [0, 2π). In Jorba (2001) it is shown that: • The stability analysis of an invariant curve of (3) is reduced to the following generalized eigenvalue problem, • If the Poincaré map is autonomous, then 1 is an eigenvalue of (10) with eigenfunction x , where x denotes the invariant curve and the differentiation with respect to θ . • Eigenvalues with norm 1 correspond to elliptic directions, and eigenvalues with norm different from 1 correspond to hyperbolic directions.
From a practical point of view, the goal is to solve a discrete version of (10). Details about how to deal with this problem numerically can be found in Jorba (2001), and will not be repeated here. In the following paragraphs we focus on how to adapt these methods to a multiple shooting scheme. Let us assume that we have computed a r -invariant curve using a multiple shooting scheme with r sections, and that we want to know its stability. Using the same argument as before to construct the linearized dynamical system (9) and the generalized eigenvalue problem (10), the stability of the r -invariant curves is reduced to the analysis of the following generalized eigenvalue problem: ⎡ where g k = g(W k (θ )), k = 1, . . . , r , Dg k is the differential evaluated on W (θ ) and T ρ denotes the operator T ρ : ψ(θ) → ψ(θ + ρ). In a more compact way, this eigenvalue problem can be expressed as This generalized eigenvalue problem is solved identically as the case r = 1. The comments in Jorba (2001) apply also this formulation of the problem. Note that in a simple shooting technique we compute the invariant curve for the map g r • · · · • g 1 . In the same way, the stability for a single shooting invariant curve is given by the eigenvalue problem The relation between the eigenvalues obtained when using single shooting with the ones obtained with multiple shooting is given by the next proposition.

Initial condition and continuation of invariant tori
As mentioned at the beginning of this section, obtaining a first invariant curve is one of the main challenges. We will use as starting point for a family of invariant tori a periodic orbit which in the Poincaré map is a fixed point of center×saddle type. Hence, we can use as first approximation the linearization of the Poincaré section around this fixed point. The initial frequency of this invariant curve is set to be ρ = ω L + ρ, where ω L is the frequency of one the elliptic directions of the periodic orbit and ρ is a small increment. The sign of ρ is positive or negative depending on whether the frequency increases or decreases when moving away from the periodic orbit along the selected elliptic direction. Then, with this initial approximation, the Newton method is applied as described in Sect. 3.1.
Hence, for now on, let us assume that a torus as expressed in (5) is known. The strategy employed here to continue a family of tori is to parameterize the family with respect to the rotation number. To find a new torus of the family the rotation number is slightly increased (or decreased, depending on which direction the family wants to be continued) as it was done to find the first torus, and then the Newton method is applied to solve for the new torus as described in Sect. 3.1. In this sense, by modifying the rotation number we are using the current torus a seed for the Newton process. This is done until three tori are computed. After the third torus, the initial condition for the next tori of the family and the rotation number are obtained by interpolating the coefficients and the rotation numbers of the previous three tori, and extrapolating them to the new one by an increment ds. This provides a good enough initial guess to find the torus in a few iterations of the Newton method.
The rotation number can be regarded as a variable or extrapolated from the previous values of the family. Both strategies were implemented for this study. Considering the rotation number a variable, in our experience, does not provide any significant benefit and requires an extra equation. The results appearing in this paper have been obtained by using the extrapolation approach but this choice is just a matter of preference.
In order to keep the number of iterations low, the extrapolation step ds needs to be adjustable. The strategy followed is to double the extrapolation step if the number of iterations is less than 6, and divide it by two if it is greater.
The process of continuing tori is not absent of challenges. Hence, we consider relevant to address the main issue found during the continuation: the sensitivity to resonances. As mentioned at the beginning of this section, the family is Cantorian. This means that it has empty interior and positive Lebesgue measure ( Jorba and Villanueva 1997). The gaps in this family are due to resonances and, typically, they are small. Hence, the continuation process jumps over them. However, there are some instances where these gaps are too big and the continuation process has difficulties to continue. In this scenario, in order to restart the process, a new initial guess for the Newton method is required. Two strategies were employed to deal with this issue. The first strategy was to increase the stepsize of the continuation parameter and check if the process jumped over the gap. This involved some trial and error, but worked in instances were the gap was small enough. The second strategy was to stop close enough to a resonance, and then transition from the BCP to the RBTP by decreasing the mass of the Sun. Once in the RTBP, the torus is a periodic orbit that can be easily continued until it crosses the resonance, and then go back to the BCP by increasing the mass of the Sun. Sometimes it is not necessary to reach the RTBP when decreasing the mass of the Sun, it is enough to lower its mass (this reduces the size of the gap) to continue the torus through the resonance and then to increase the mass to be again in the BCP. We consider that the new torus obtained belongs to the same family if a reasonable condition is satisfied: the torus obtained by continuing from the RTBP is close enough to the last torus of the family.
The next section describes how these techniques were applied to find families of invariant tori in the BCP.

Dynamics around the L 2 point
This section is devoted to understand the main one-parametric families of two-dimensional invariant tori that exist nearby L 2 in the BCP. While the perturbation due to the Sun is large enough to modify the geometrical structure of the phase space near the translunar point (as it is discussed in Sect. 2), the situation is similar in an extended neighborhood. Let us recall that, once the periodic time dependent perturbation is turned on (ε = 0), the periodic orbits of the planar and vertical Lyapunov families gain, generically, the frequency of the perturbation and become families of two dimensional invariant tori (this is illustrated in Fig. 4, for a particular Halo orbit). As we have mentioned already, the periodic orbits whose period is a Fig. 4 Transition from Halo orbit with energy -1.5244988379312372 in the RTBP (green) to a torus in the BCP (red). The torus is the dynamical equivalent (in the BCP) to the periodic orbit in the RTBP and has rotation number 1.3800185497627542 rational multiple of the frequency of the Sun do not become two dimensional invariant tori but, remain being periodic. Analogously, the two-dimensional invariant tori contained in the center manifold of L 2 become, with the perturbation due to the Sun, three-dimensional tori, gaining the frequency of the perturbation. Again, if one of the frequencies of a two dimensional torus is a rational multiple of the one of the Sun, these torus remains two-periodic. Summarizing, there are two main mechanisms to obtain two-dimensional invariant tori near the translunar point in the BCP: from generic perturbations of periodic orbits or from resonant perturbations of two-dimensional invariant tori. The previous claims on the behaviour of periodic and quasi-periodic motions under periodic perturbations are theoretically supported in Jorba and Villanueva (1997). In this section we make an effort to identify the main families (the ones coming from periodic orbits and the ones coming from low order resonances) and understand how these are related to each other.
As explained in Sect. 3.3, the continuation process requires an initial torus. This initial torus usually is computed from a periodic orbit. In the context of the BCP, two kind of initial periodic orbits were used to find and continue families of invariant tori. The first kind of periodic orbits are Halo orbits from the RTBP. The approach is to pick a Halo orbit in RTBP, and then continue it with respect to ε until it reaches the BCP (ε = 1). As explained before, the original periodic orbit is casted to a two-dimensional invariant torus once ε > 0. Figure 4 shows an example, in different projections, of how a Halo orbit in the RTBP becomes a quasi-periodic orbit in the BCP with two frequencies: the intrinsic one corresponding to the Halo orbit, and the one acquired due to the Sun's perturbation. The two-dimensional tori continued from the Halo orbit belongs to a family in the BCP. Therefore, it can be continued with respect to the non-trivial frequency. Family of tori, however, have gaps corresponding to resonances (this will be discussed later with more detail). The gaps hinder the numerical continuation of the family. A solution to overcome this difficulty is to pick up several Halo orbits from the RTBP, continue them with respect to ε, and then, at the BCP, with respect to the non-trivial frequency. In this way, a complete description of the family can be obtained.
The other periodic orbit used was the one found by continuing the L 2 point from the RTBP to the BCP, i.e., the orbit described in Sect. 2. The latter generates a family of planar quasiperiodic orbits that can be considered the quasi-periodic planar Lyapunov family counterpart of the periodic ones in the RTBP.
The result of computing and continuing families of 2D invariant tori is showed in Fig. 5. The horizontal axis represents the x component of the corresponding invariant curve of the stroboscopic map when θ = 0. The vertical axis is the rotation number. Several resonances have been identified in Fig. 5 to illustrate the argument made in Sect. 3.3 about the gaps in the family. A total of six families were found. Two of these families are planar Lypaunov-type quasi-periodic orbits (families H1 and H2 in Fig. 5), and four have a vertical component. Out of these four, two are Halo-like quasi-periodic orbits (for the moment being, we refer to them as Halo families of Type I and Type II, see Fig. 5 Moon (families V1 and V2 in Fig. 5). In Fig. 6 we show magnifications of several parts of Fig. 5 corresponding to remarkable events and to the different tori plotted in this paper. The Type I Halo family is obtained by continuing some Halo orbits from the RTBP to the BCP, and then continuing them in the BCP. See Table 3 in "Appendix 1" for details on the RTBP Halo orbits used. This family is discussed in greater detail in Sect. 4.1.
The H1 family originates from the periodic orbit obtained by continuing the L 2 point from the RTBP to the BCP (see Fig. 3). The stability of family H1 was analyzed, and most of the tori are hyperbolic. There is always an eigenvalue equal to 1 with multiplicity two, plus one real eigenvalue of the order of 10 6 (and its inverse), and another pair that evolves in a way that the family undergoes two bifurcations. This is illustrated in Fig. 7, where the last pair of eigenvalues are plotted. The horizontal axis corresponds to the x component of the invariant curve, and the vertical axis the absolute value of the eigenvalue. Figure 7 shows that there are two bifurcations where the absolute value of the eigenvalues is equal to one. In these cases, there are two small intervals that contain partially elliptic tori; this is, that the eigenvalues are complex with norm equal to one. These small intervals are zoomed in Fig. 8. The top row of Fig. 8 shows the absolute value of the eigenvalues, and the bottom row the arguments. Note that a similar phenomena appears in the RTBP, where the planar Lyapunov family undergoes a bifurcation that gives rise to the well-know family of Halo orbits. The same happens in the BCP for these two bifurcations. Each one of this families can be continued along the vertical component. Looking at the family H1 in Fig. 5, from left to right, the first bifurcation gives rise to the Type II Halo family, while the second one to the V1 family. During the continuation of the V1 family, it was found that some small resonances needed to be avoided. The strategy of going back to the RTBP by decreasing the mass of the Sun, continuing the resulting object there until the resonance is passed, and going back to the BCP was employed. After returning to the BCP, it was noticed that the resulting torus did not belong to the V1 family, but to a new one labeled as V2. This torus was continued, both increasing and decreasing the rotation number. Eventually, the V2 branch met a planar quasi-periodic Lyapunov orbit of a new family, called H2. Again, this family was continued, hence completing the picture represented in Fig. 5. A complete study of the H and V families is left for another work, although some examples are provided in the "Appendix". The next subsection elaborates on the Type I and Type II Halo-like families, the focus of this paper.

The Type I and Type II Halo-like families
Let us begin showing some representative examples of the members of these two families. The first example of Type I Halo family is shown in Fig. 4. The green curve represents a Halo orbit in the RTBP and, in red, we display the torus of the BCP with two frequencies: the one of the Halo orbit and the frequency of the Sun (ω S ). This 2D-torus is seen as an invariant curve of the stroboscopic map with rotation number ρ = 1.380018549762754. Another example is shown in Fig. 9. In this case, the rotation number (in the stroboscopic map) is ρ = 2.675226847819367. This torus is close to the resonance value of ρ = 6π/7 ≈ 2.6927937 . . .. The effect of being close to a resonance is illustrated in Fig. 10, a torus with rotation number (in the stroboscopic map) ρ = 2.692464347819371. Figure 11 shows a torus of the Type II family. This particular example has rotation number (in the stroboscopic map) ρ = 3.116137168026786. The projection on the x − z plane shows that the orbit is a Halo-like in the sense that when observed from the Earth, the orbit circles around the L 2 Similarly to the case in the Type I Halo family, near a resonance we observe the same phenomena, and the orbit becomes more dense around the periodic orbit corresponding to that resonance. Figure 12 provides and example with rotation number equal to ρ = 3.130357871578353, close the resonance ρ = π.
As it has been mentioned before, the Type I family of Halo-like orbits appears when adding the Sun effect to the family of Halo orbits of the RTBP: the (non-resonant) Halo orbits add the frequency of the Sun to its own frequency and become a quasi-periodic orbit with two basic frequencies. Therefore Type I family is to be understood as the dynamical equivalent in the BCP to the classical Halo family of the RTBP.
To better understand the Type II family, we study its original counterpart in the RTBP. To do so, we continue some orbits of Type II family to the RTBP by decreasing ε down to zero. As an example, in Fig. 13 (top row) we have plotted two RTBP orbits that come from the continuation of the Type II orbits with rotation numbers ρ = 0.739476685309787 and ρ = 0.858771705123796. It is clear from Fig. 13 (top row) that the ones obtained are not periodic but quasi-periodic Halo orbits, however, we can provide more details on the role these quasi-periodic Halo orbits play in the neighborhood of L 2 . To do so, we have to perform a reduction to the center manifold related to L 2 in the RTBP. This technique is rather complicated and is completely out of the scope of the present paper to discuss it. To provide a sufficient understanding on the procedure we will say that it consists of a change of variables (a partial normal form) that uncouples, up to high order, the hyperbolic and elliptic directions related to L 2 . Then, we can restrict the Hamiltonian to the elliptic directions to obtain a two degrees of freedom Hamiltonian system that contains the stable motion in the neighborhood of L 2 . This allows to study the periodic and quasi-periodic motions near L 2 without being affected by the high instability of the region. This technique is applied in Jorba and Masdemont (1999) to study the neighborhood of the Earth-Moon collinear points. In Jorba (1999), a deep discussion on the method together with practical details on the implementation is presented.
Once the reduced Hamiltonian is obtained, one can study it as any other system with two degrees of freedom. In particular, one can use a Poincaré section to reduce the dimension of the system and the Hamiltonian energy to slice the resulting phase space, getting therefore, a family of (two-dimensional) Area Preserving Maps (APM). Notice that, by Hamiltonian energy, we mean the one of the reduced Hamiltonian, not the one of the RTBP. Again, this extra reduction is used in Jorba and Masdemont (1999) to visualize the stable motion around the collinear points. Moving the parameter of the APM (the normalized energy), the authors are able to show how the stable motion around the L 2 is organized when the bifurcation that gives rise to the Halo family takes place.
In Fig. 13 (bottom row), we display the phase portraits of the Hamiltonian of the RTBP restricted to the center manifold of L 2 for two different energy values. The axis in this figure, (q 2 , p 2 ), are coordinates build specifically to represent the center manifold, see (Jorba and Masdemont 1999) for details. One can observe typical features of area preserving maps: Fixed points and invariant curves. The boundaries of the regions displayed in Fig. 13 (bottom row) correspond to planar Lyapunov orbits (when sent back to the RTBP) while the fixed point at the origin correspond to vertical Lyapunov orbits. The two extra fixed points correspond to Halo orbits of the two (symmetric) families of Halo orbits: North (the fixed point on the left) and South (the one on the right). These fixed points corresponding to Halo orbits are surrounded by invariant curves that, when transformed back to the coordinates of the RTBP, correspond to quasi-periodic Halo orbits. Recall that, in Fig. 13 (top row) we display samples of Type II orbits continued back to the RTBP. By means of the change of coordinates, we have sent initial data of each orbit to the center manifold coordinates. Then, we have plotted a Poincaré map for the level of energy of each orbit and we have marked the initial data of each orbit in the map with a big dot (with the same color used to plot the orbits). The results are shown in Fig. 13, down. The axes of this figure represent suitable normalized coordinates to display the center manifold. To define this set of coordinates is out of the scope of this paper, see (Jorba and Masdemont 1999) for a deep discussion.
This shows that the Type II orbits come from quasi-halo orbits of the RTBP that have one of its two frequencies in resonance with the frequency of the Sun. In this way, the effect of the Sun does not add a new frequency and the quasi-halo is continued into the BCP as a quasi-periodic orbit with two basic frequencies that we refer as Type II.

Stability
To fully characterize these orbits, we study their stability. The stability of the tori, in this case, can be characterized by six paired eigenvalues. As the stroboscopic map is autonomous and symplectic, all the tori have a trivial pair of eigenvalues (equal to the unity). The remaining eigendirections may be elliptic or hyperbolic. In the families we study in this section, the tori have a very unstable/stable pair, i.e., a very large eigenvalue and its inverse. In the case of Type I family, the largest eigenvalue ranges from 2300 to 318600 (approximately) while, for the Type II family, the largest eigenvalue ranges from 23400 to 794260 (approximately). Therefore, the Type II family is less stable. The remarkable changes in the spectrum of the family take place in the remaining pair, in which we focus on.
Using the method described in Sect. 3.1, the stability of all the tori computed for each one of the families is obtained. For the Type I Halo family, they mostly behave like their counterparts in the RTBP, the Halo orbits. Due to the Hamiltonian structure, there is always the eigenvalue 1 with multiplicity two. For each tori of this family there is a large real eigenvalue (and its inverse), and, for almost each tori, a complex eigenvalue (and its inverse)  with modulus 1. The absolute value of the latter pair of eigenvalues is shown in Fig. 14 (left) with respect to the x component of the invariant curve at θ = 0. It is observed that most of these pairs of eigenvalues have modulus 1 with the exception of some isolated zones. However, the main takeaway is that most of the tori are partially elliptic with one saddle. On the other hand, the Type II Halo family has a different stability type. In this case, and as in the case of the Type I Halo family, there is always the eigenvalue 1 with multiplicity two. There is also a large real eigenvalue (and its inverse). The other pair, however, is also real and positive. Figure 14 (right) shows the evolution of this eigenvalue with respect to the x component of the invariant curve at θ = 0. Hence, the Type II Halo family has two saddles. We note that the largest eigenvalues of the Type I and Type II families are of the same order of magnitude.

Applications
The existence of two Halo-like families illustrates a resonance between the direct effect of the Sun's gravity, as modeled in the BCP, with a quasi-Halo orbit of the RTBP. We emphasize the dependency on how the effect of the Sun is accounted for because, for example, the QBCP also models the direct effect of the Sun's gravity but, as of today, only the quasi-periodic counterparts of the Halo orbits (Type I family) have been computed (see Andreu 1998). The existence of Type II Halo-like orbits provides mission analysts with new potential candidates to meet the requirements for missions to the vicinity of the Moon. Finally, let us comment a bit more on the Type II Halo family, that shares some topological features with Type I Halo family. Note that there are representatives members of each family that are not blocked by the Moon, making them useful for missions to the neighborhood of the Moon that require continuous line-of-sight with the Earth. Figure 15 shows the projection on the x = 0 plane corresponding to how these orbits would be seen from an observer in the Earth. The projections in Fig. 15 correspond to the same orbits shown in Figs. 11 and 12. In these figures, the center of the Moon is at the origin, and it has been plotted a circle with the approximated radius of the Moon, and another one with a circle twice the radius of the Moon. In both cases it is observed that there is continuous line-of-sight between the Earth and the orbit.

Conclusions and further work
In this paper we have explored the dynamics of a massless particle around L 2 in the Bicircular Problem. In this model, the L 2 point is only defined geometrically because it is not an equilibrium point. By means of a continuation scheme with respect to the Sun's effect, we showed that there is no natural dynamic replacement of the L 2 point in the BCP.
Continuation of invariant tori families was the technique adopted as an alternative to the reduction to the center manifold, given that the size of the domain of validity of the expansions is too small to provide information about the dynamics in a reasonable neighborhood of the L 2 region. Following this approach, we have identified a total of six families. Out of these six families, two of them were planar quasi-periodic orbits, and the other four had vertical component. Two out of these four non-planar families were Halo-like orbits. The first Halo-like family, called Type I Halo, were obtained by sampling and continuing their RTBP counterparts from the RTBP to the BCP. The second family, Type II Halo, was found by analyzing the bifurcation of the planar family H1.
The stability of the Halo-like orbits was computed. One family, the Type I Halo family, can be seen as the natural continuation of the classical Halo orbits of the RTBP and share the same stability type of their RTBP counterparts (saddle×center×center). On the other hand, the Type II Halo family comes from a quasi-Halo orbit of the RTBP which has one frequency in resonance with the frequency of the Sun and their stability is of the type saddle×saddle×center. The shape and location of these two families makes them suitable to space missions, as they allow permanent communication with the Earth (see Fig. 15). Let us remark that, given its stability character, Type I Halo family probably more adequate for practical purposes. There is, however, more work to be done. The future work is summarized in the following paragraph.
As first step, more work needs to be done to study the other families. We also believe that there are some more families of Halo-like quasi-periodic orbits with two basic frequencies, and more work is needed to find them. Although with no obvious application to space missions, they still have academic interest. In parallel, and focusing on mission analysis and to the potential applications of the Halo-like families to mission design, we need study these families in a real ephemeris model and to develop station-keeping strategies. Second, how to transfer from the vicinity of the Earth to one of these orbits would also be of interest to the mission designer . To that effect, the use of the stable/unstable manifold would be very useful, as it has been proven very successfully in other contexts. This also would help to get a deeper understanding of the dynamical skeleton of the BCP around the L 2 point.

Conflict of interest
The authors declare that they have no conflict of interest.

Appendix
In this section we provide additional information concerning the paper. In Table 3 we provide the initial conditions (identified by the energy) of the Halo family of the RTBP that have been used to produce Fig. 5.
We also provide some examples of tori from the other families found are given (see Fig. 5). They are provided here to illustrate the richness of the Sun-Earth-Moon BCP, and to evidence that the vertical families V1 and V2 are not Halo-like. A complete study of their stability properties and how they transition from the RTBP to the BCP is under work, and no details are provided here.
The planar tori from the families H1 and H2 are very similar, and two examples of each one are shown in Fig. 16. The representative of the family H1 (left) has rotation number ρ = 0.522687812628674. The rotation number of the representative of the family H2 (right) is ρ = 0.258684108104417.
More interesting are the families V1 and V2 with a vertical component. Different projections of a representative of the family V1 with rotation number ρ = 0.651014628070470 are  Fig. 6d) shown in Fig. 17. The projection onto the plane x = 0 (bottom-left image) shows that this orbit falls behind the Moon. Finally, and example of the family V2 is illustrated in Fig. 18. This torus has rotation number ρ = 0.585297052915989. It also falls behind the Moon. However, the projection onto the plane x = 0 (bottom-left image) show that has different symmetry than the representative of the V1 family.