Numerical study of the geometry of the phase space of the Augmented Hill Three-Body problem

The Augmented Hill Three-Body problem is an extension of the classical Hill problem that, among other applications, has been used to model the motion of a solar sail around an asteroid. This model is a 3 degrees of freedom (3DoF) Hamiltonian system that depends on four parameters. This paper describes the bounded motions (periodic orbits and invariant tori) in an extended neighbourhood of some of the equilibrium points of the model. An interesting feature is the existence of equilibrium points with a 1:1 resonance, whose neighbourhood we also describe. The main tools used are the computation of periodic orbits (including their stability and bifurcations), the reduction of the Hamiltonian to centre manifolds at equilibria, and the numerical approximation of invariant tori. It is remarkable how the combination of these techniques allows the description of the dynamics of a 3DoF Hamiltonian system.


Introduction
In the recent years, a lot of attention has been devoted to the study of asteroids and comets. On one side, they carry information about the past of the Solar System. On the other hand, we are starting to look at them as potential resources of raw minerals and volatiles, to be used in-situ for space exploration, or to be sent back to Earth. Moreover, the so-called Near-Earth Objects (NEOs) are also of high interest due to the risk of collisions with the Earth. Solar sailing is a novel way of navigating an unmanned spacecraft in the Solar System, and it could allow multi-rendezvous missions to visit several asteroids (Dachwald et al. 2014). It is remarkable that, as the gravitational field of an asteroid is small, solar radiation pressure becomes relevant when moving around them (Villac et al. 2012). For this reason, the use of the solar radiation pressure to navigate around asteroids can give rise to new mission concepts.
The three dominant effects for a small solar sail moving near an asteroid are the gravitational attraction of the asteroid, the gravitational field from the Sun and the solar radiation pressure (SRP). The first effect also includes the shape and rotation of the asteroid, which are relevant within a few radii from its surface (this is known as the gravity regime). If the distance to the asteroid is larger, it is enough to consider it as a point mass. This is the case in which the solar gravitation and SRP are relevant and must be considered. Moreover, as the distance to the Sun is much larger than the distance to the asteroid, Sun's gravity and SRP can be taken as uniform forces. For this reason, we use as a model for the dynamics the Augmented Hill 3-Body problem, AH3BP (Morrow et al. 2001;Giancotti and Funase 2012;Farrés et al. 2014b), which is the classical Hill problem (Hill 1878;Szebehely 1967) with an additional term that accounts for the effect of the solar sail. This term depends on the sail's efficiency, defined by two parameters (β, ξ r ), and orientation, parameterised by two angles (α, δ). The parameters (β, ξ r ) are the sail lightness number and the reflectivity coefficient, and the angles (α, δ) define the orientation of the sail. It is true that asteroids follow eccentric orbits around the Sun, but in most cases their eccentricity is small. For this reason and because a good understanding of the no eccentricity case is still needed we have neglected this effect. A preliminary study on the dynamics of a solar sail in the vicinity of Vesta using the elliptic H3BP can be found in Farrés and Jorba (2012). The Augmented Hill model considered here is an extension of the so-called Photogravitational Hill problem, in which the solar radiation pressure only acts in the opposite direction of the Sun, see Markellos et al. (2000), Kanavos et al. (2002), Papadakis (2006), Liu and Villac (2010), Giancotti et al. (2014) and Yárnoz et al. (2015). Preliminary studies on the dynamics of a solar sail in the Hill model are Morrow et al. (2001), Morrow et al. (2002), Giancotti and Funase (2012) and Farrés et al. (2014a, b).
The goal of this paper is to provide a numerical description of bounded motions near equilibrium points of the AH3BP that could be useful for a space mission around an asteroid. We also want to emphasise how both classical and more recent numerical methodology can be used to have a complete description of a 3 degrees of freedom Hamiltonian system. For the sake of completeness, we have included a section where we briefly describe the methods that we have used.
It is well known (Hill 1878;Szebehely 1967) that the classical Hill problem has two symmetric equilibrium points at (±3 −1/3 , 0, 0). The linear dynamics around these points is of the type centre × centre × saddle, which gives rise to families of periodic and quasiperiodic orbits that are part of the centre manifold. By adding the effect of a given solar sail (i.e. (β, ξ r ) are fixed), the equilibrium points are displaced and give rise to a 2-parametric family of equilibria, parameterised by the two angles (α, δ) that define the orientation of the sail. As the effect of the sail is small, these new points inherit the linear stability of the unperturbed equilibria, together with its surrounding dynamics. It is remarkable that, for a set of parameter values of practical interest, the two linear frequencies of the equilibrium point coincide.
In this paper, we perform a numerical description of the nonlinear dynamics around the equilibrium points corresponding to a representative set of sail orientations, that has been chosen in order to show its global effect on the dynamics. This description is always done in three stages. We begin by performing numerical continuation of the families of periodic orbits that are born at the equilibrium point, together with their main bifurcations, that provide the underlying skeleton that organises all the dynamics. After that, a part of the dynamics surrounding this skeleton is found through the use of two techniques: reduction to the centre manifold and numerical continuation of quasi-periodic motion. Of these two techniques, the first one, which is actually semi-numerical, provides a (truncated) asymptotic expansion of the centre manifold of the equilibrium point that enables the computation of trajectories by direct numerical integration. Due to its asymptotic character, this expansion is accurate in a neighbourhood of moderate size of the equilibrium point. Through the second technique, quasi-periodic motion (which, through reduction to the centre manifold, is found to fill most of the phase space near the equilibrium point) can be continued further away.
This paper is structured as follows. Section 2 describes in detail the AH3BP. Section 3 summarises the computational methodology used throughout the paper, in order to provide context and notation to the numerical results. This is presented in a separate section to simplify the reading. Section 4 describes the families of equilibrium points obtained by varying the sail parameters. This section motivates the choice of the actual sail parameters for the numerical study of periodic and quasi-periodic orbits performed in Sect. 5. In this last section, we also perform a partial numerical unfolding of the 1:1 resonance.

The AH3BP model for a non-perfectly reflecting solar sail
The Augmented Hill 3-Body problem (AH3BP) models the motion of an infinitesimal particle (sailcraft) that interacts with a small mass (asteroid) and is perturbed by a distant large body (Sun), and also takes into account the effect of solar radiation pressure (SRP). As the Sun-sail distance is large compared to the asteroid-sail distance, the effect of the Sun is included as a uniform gravity field and a uniform SRP.
We consider a rotating reference frame centred on the asteroid defined such that: the x direction points from the Sun to the asteroid; the z direction is aligned with the Sun-asteroid angular velocity; and the y direction completes a positive coordinate reference system. We use normalised units such that L = (μ sb /μ sun ) 1/3 R is the unit of distance, and T = 1/ω is the unit of time. In the previous expressions μ sun = Gm sun represents the gravitational parameter of the Sun, μ sb = Gm sb the gravitational parameter of the asteroid, R the Sun-asteroid mean distance, and ω = μ sun /R 3 the angular velocity of the asteroid around the Sun. Using these units, the equations of motion for a solar sail in the vicinity of an asteroid (Liu and Villac 2010;Villac et al. 2012) arë where Farrés et al. Here (X, Y, Z ) denotes the position of the sailcraft in the rotating frame, r = √ X 2 + Y 2 + Z 2 is the distance between the asteroid and sailcraft, and a sail = (a X , a Y , a Z ) is the acceleration due to the solar sail.
To model the acceleration given by the solar sail we take into account the effect of both reflection and absorption of the photons by sail material (McInnes 1999;Dachwald et al. 2005). The force due to reflection is given by F r = 2P A r s , n 2 n while the one due to absorption is F a = P A r s , n r s , where P is the SRP magnitude at a distance R from the Sun, A is the solar sail area, r s is the SRP direction and n is the normal direction to the surface of the sail. We denote by ξ a the absorption coefficient and by ξ r the reflectivity coefficient, that must satisfy ξ a + ξ r = 1. Hence, the acceleration of a solar sail is given by We note that ξ r = 1 corresponds to a perfectly reflecting solar sail, and ξ r = 0 to a solar panel, where the only effect is the absorption of the photons by the panels. We measure the effectiveness of a solar sail in terms of the sail lightness number β = 2P A/m. As we use normalised units of distance and time, the parameter β is re-scaled by a factor (μ sb ω 4 ) −1/3 , where ω = μ sun /R 3 (Liu and Villac 2010). We call the normalised sail lightness numberβ = β(μ sb ω 4 ) −1/3 . Following Broschart et al. (2014), we have that 8502 when A is given in m 2 and m in kg. Notice that β depends on both: the area-to-mass ratio of the satellite (σ = (A/m)) and the gravitational parameter of the small body (μ sb ). Hence, the sameβ can represent two different missions scenarios. For further details on values ofβ depending on the sail performance and the asteroids mass see Farrés et al. (2014a, b).
We parameterise the normal direction to the surface of the sail (n) using two angles α, δ, which represent the horizontal and vertical displacement of n with respect to r s . In the rotating reference frame r s = (1, 0, 0) and n = (cos α cos δ, sin α cos δ, sin δ). Notice that the sail normal direction cannot point towards the Sun, hence α, δ ∈ [−π/2, π/2].
With these definitions, a sail = (a X , a Y , a Z ) in Eq.
Moreover, the AH3BP as given by Eq.
(1) admits an integral of motion, also known as the Jacobi constant, If we introduce the classical momenta P X =Ẋ − Y , P Y =Ẏ + X and P Z =Ż , the system can be written in Hamiltonian form as The Hamiltonian satisfies H = J c /2. Throughout the paper, we focus on the particular case ofβ = 5 and ξ r = 0.85, that corresponds to a sailcraft with similar characteristics to NanoSail-D2 orbiting close to Ceres (Farrés et al. 2014b). Our goal is to describe the dynamics of the system and how this one varies when we change the sail orientation. This knowledge is very useful for practical mission applications when defining a nominal orbit and its control (Farrés and Jorba 2008a).

Computational methodology
As it has been explained in the introduction, the goal of this paper is to provide a numerical description, as complete as possible, of bounded motions around equilibrium points of the AH3BP that could be useful for a space mission around an asteroid. Once an equilibrium point has been chosen, our strategy is, as also mentioned in the introduction, to perform numerical continuation of families of periodic orbits, reduction to the centre manifold, and numerical continuation of invariant tori. The corresponding numerical and semi-numerical methodologies are summarised in Sects. 3.1, 3.2, and 3.3, respectively. References to the literature with full details are also provided.

Computation of periodic orbits
It is well known that, if φ t (x) is the flow associated with an autonomous ODE, a periodic orbit corresponds to a fixed point of the Poincaré map: is transversal to the flow). Hence, finding a periodic orbit is equivalent to finding a zero of G(x) = P(x)− x, which can be solved using the classical Newton method (Simó 1990). Due to the instability of the system close to equilibrium points, we use multipleshooting (Stoer and Bulirsch 2002) to avoid problems with the convergence of Newton's method and increase the accuracy.
As the system is Hamiltonian, periodic orbits are organised in one-parametric families. Hence, all the orbits in the family are solution of G(x) = 0. In order to have unicity we must add an extra condition. We have choosen to fix the period of the orbit (τ (x) = T ), although we could also fix the energy of the system (H (x 0 ) = h). We continue the corresponding families via a predictor-corrector method, using the pseudo-arc-length as continuation parameter. In this way, we avoid having problems with possible turning points (Simó 1990).
We also study the stability of the periodic orbits, as changes in the stability can lead to possible bifurcations (as in happens with the halo orbits in the classical Hill problem). To study the stability of a periodic orbit, one must check the spectrum of the monodromy matrix M(T ) = Dφ T (x 0 ) (also known as state transition matrix). As the system is Hamiltonian, we know that the eigenvalues appear in pairs. Hence, the spectrum is: {1, 1, λ 1 , λ −1 1 , λ 2 , λ −1 2 }. Following Hénon (1965), we define the stability parameters as s i = λ i + λ −1 i for i = 1, 2 (Gómez and Mondelo 2001;Gómez et al. 2005). Notice that if |s i | > 2 we have an unstable direction and if |s i | < 2 we have a linearly stable one. We use these parameters as indicators to detect possible bifurcations whenever |s i | = 2.

Computation of centre manifolds
A centre manifold of an equilibrium point is an invariant manifold tangent to the centre part of the linearised vector field at the equilibrium point. Here, we focus on the dynamics near the L 2 point of the Hamiltonian (6) and, in this case, the centre manifold is a 4D invariant manifold which is normally hyperbolic due to the saddle directions. All the periodic and quasi-periodic orbits are contained in this manifold and, as we will see, the computation of this manifold allows the visualisation of the dynamics in an extended neighbourhood of the equilibrium the point. The centre manifold can be obtained by means of a partial normalisation (up to a finite order) of the power expansion of the Hamiltonian at the equilibrium point (Jorba 1999;Jorba and Masdemont 1999). Another option is to compute a power expansion of a parameterisation of the manifold. For more details on this approach, see Farrés and Jorba (2010a) and Haro et al. (2013). In both cases, the validity of the approximation is restricted to the domain of convergence of the expansions. Classical results on the existence and (lack of) uniqueness of these manifolds can be found in Carr (1981), Sijbrand (1985) and Vanderbauwhede (1989).
In this paper, we will use a partial normalisation on the Hamiltonian of the AH3BP. A very preliminary version of the centre manifold computations in this paper can be found in Farrés and Jorba (2016b). The main idea is to use the Lie series method: if G(q, p) is a Hamiltonian system, then the functionĤ defined bŷ is the result of applying a canonical change to H . This change is the time one flow corresponding to the Hamiltonian G. G is usually called the generating function of the transformation to obtainĤ . As the change is canonical,Ĥ is the Hamiltonian of H in the new reference system, see Giorgilli et al. (1989) and Simó (1989) and references therein for more details. This scheme is applied on a power series expansion of H at the equilibrium point, and G is also expressed as a power series, see Jorba (1999), Jorba and Masdemont (1999), Farrés and Jorba (2010a), Celletti et al. (2015) and Ceccaroni et al. (2016b) for more details. As we do not want a complete normal form, we will only suppress some specific monomials during the process. To explain which monomials are kept, let us introduce some notation. Assume that the origin of coordinates is at the equilibrium point, which means that the power expansion of the Hamiltonian takes the form where H j (x, y) is an homogeneous polynomial of degree j. Moreover, assume that a suitable linear transformation has been applied so that H 2 has the form The goal is to suppress all the monomials in H 3 , H 4 , ..., H N such that the exponent of x 1 is different from the exponent of y 1 (the reason will be clear later). This normalising scheme is performed up to order N , for N between 16 and 32, which gives a very good balance between accuracy and the size of the domain where the power expansions are valid. In general, normalising to all orders is a divergent process, which means that the domain of convergence of the expansions shrinks to zero when the normalising order goes to infinity. In this particular case, as this process has no small divisors (see Jorba and Masdemont 1999), the reduction of the domain is very slow with the normalising order, giving rise to a quite large domain of validity of the expansion for the orders considered. Let us denote as (q, p) the new variables. After neglecting the remainder, the Hamiltonian has the form H = H (q 1 p 1 , q 2 , p 2 , q 3 , p 3 ), where we note that it does not depend of q 1 , p 1 but on their product q 1 p 1 . The key point is that the quantity I 1 = q 1 p 1 is a first integral of the Hamiltonian. Choosing I 1 = 0, we obtain a Hamiltonian system, with two degrees of freedom, of the form where q = (q h , q v ) ∈ R 2 , p = ( p h , p v ) ∈ R 2 and H cm 3 , . . ., H cm N denote the monomials of degrees between 3 and N . The values ω h and ω v are the frequencies of the linear oscillations at the point, that give rise to the horizontal and vertical families of Lyapunov orbits, respectively. It is important to note that, in all cases considered, the frequencies ω h , ω v are positive. This implies that, for h > 0 small enough, the level surfaces H cm = h are bounded and can be visualised as (deformed) spheres in R 4 . Moreover, this also implies that the periodic and quasi-periodic orbits near the point have a value of the energy larger than the energy level of the equilibrium point. In fact, for each energy level, there are two Lyapunov periodic orbits (related to the horizontal and vertical linear oscillations) inside this level.
As a final comment, we remark that this technique does not allow to describe the dynamics far from the equilibrium point. This is due to the use of power expansions and, therefore, they only converge up to the first singularity. For the initial expansion at the equilibrium point, the first singularity corresponds to the position of the asteroid. Then, this radius of convergence is slightly reduced by each normalising transformation and, therefore, the validity of the final Hamiltonian is restricted to a ball centred at the equilibrium point, and radius slightly below the distance to the asteroid.

Computation of invariant tori
The method used in this paper for the computation and continuation of invariant tori is based on Castellà and Jorba (2000). It consists on looking for closed curves, approximated as Fourier series, that are invariant by the flow time one of the periods of the torus. We look for the Fourier coefficients of a curve that satisfies the invariance equation From the parameterisation of such a curve, a parameterisation of a 2D torus that contains it can be recovered as This torus has frequencies ω 1 = ρ/T , ω 2 = 2π/T . Eq. (8) is a functional equation that is turned into a finite-dimensional, nonlinear system of equations by considering a grid of 1 + 2N f values of θ equally spaced in [0, 2π]. In this equation, the unknowns are ρ, T and the Fourier coefficients The parameters T, ρ can be used to describe the families of tori that will be computed. Nevertheless, in order to relate these families to other invariant objects computed with other methodologies, an energy equation is added to the discretised version of (8), in order to have the value of the Hamiltonian h as a parameter. This final discretised system allows for both the computation of individual tori and the continuation of families. For continuation, we use a standard pseudo-arc-length procedure with step size control (Allgower and Georg 1990). For this procedure to work, we take several values of one of the parameters h or ρ, and perform numerical continuation of the one-parametric families of tori corresponding to these fixed values of the parameter.
For a fixed energy level in which there is a periodic orbit with central part, meaning that one of its stability parameters s i (see Sect. 3.1) is such that |s i | < 2, the linearised flow around the periodic orbit possesses a one-parametric family of invariant tori with constant rotation number ρ = arccos(s i /2). For the full, nonlinear flow, this family persists as a Cantorian family of tori in which ρ varies (Jorba and Villanueva 1997). This family could be obtained numerically by numerical continuation of our discretised system, taking as initial approximation a torus of the linear flow close to the periodic orbit. Since the families of periodic orbits have central part in a range of energy levels, the families of invariant tori around them are 2-parametric. In order to avoid problems related to the non-existence of tori for non-Diophantine frequencies, we do not continue invariant tori for constant energies varying ρ, but for constant ρ while increasing the energy. The values of ρ chosen are always noble numbers (continued fraction expansion equal to one from a coefficient on up to double precision). Explicit formulae and additional technical details can be found in Gómez and Mondelo (2001).

Equilibrium points of the AH3BP
The classical Hill model has two equilibrium points L 1 and L 2 , symmetrically located around the origin at (±3 −1/3 , 0, 0) (Scheeres 2012), and we can "artificially" displace them with the effect of a solar sail. In this section, we describe how the different parameters of the sail (β, α, δ, ξ r ) affect the position of the equilibrium points. This will allow us to understand the effect that the different parameters of the system have on the position of the equilibrium points and their stability. For further details see Morrow et al. (2001), Giancotti and Funase (2012) and Farrés et al. (2014a, b).
The equilibrium points of Eq. (1) are given by: Hence, for a given solar sail performance (i.e.β and ξ r fixed), there is a 2D surface of equilibrium points parameterised by the angles α and δ. In order to compute these surfaces of equilibria, we find the zeros of F(X, Y, Z , α, δ) : U ∈ R 5 → R 3 . Note that for α = δ = ±π/2 (no solar radiation pressure effect, since the sail is parallel to the Sun direction) we already know that the equilibrium points are (±3 −1/3 , 0, 0). The other solutions are found using classical continuation methods, with pseudo-arc-length as continuation parameter (Simó 1990). Since 2D surfaces are hard to compute, we always fix a relation between the two sail orientations in order to slice the surface of equilibria. Namely, we fix α = 0 and vary δ ∈ [−π/2, π/2] and vice-versa. When we include the SRP (β = 0) and the sail is perpendicular to the Sun-line (α = δ = 0), asβ increases L 1 and L 2 move towards the Sun along the x-axis (Scheeres 2012). For example, forβ = 5 and ξ r = 1, we have SL 1 = (−1.7727361696, 0, 0) and SL 2 = (0.40146718344, 0, 0). If we change ξ r we are essentially changing the performance of the solar sail to β * =β(1 + ξ r )/2 moving the equilibria away from the Sun. When we fixβ and Family of equilibrium points forβ = 5, ξ r = 0.85. Left α ∈ [−π/2, π/2] and δ = 0 (all points are on the Z = 0 plane); right α = 0 and δ ∈ [−π/2, π/2] (all points are on the Y = 0 plane). The arrows represent n at each equilibrium point

Fig. 2
Family of equilibrium points forβ = 5, ξ r = 0.0, 0.30, 0.75, 1.0. Left α ∈ [−π/2, π/2] and δ = 0 (all points are on the Z = 0 plane), right α = 0 and δ ∈ [−π/2, π/2] (all points are on the Y = 0 plane). The arrows represent n at each equilibrium point change the sail orientation, we displace the equilibrium points in other directions: changes in α imply shifting their position to one side or the other form the Sun-asteroid line, while changes in δ imply displacing their position above or below the ecliptic plane. In this case, changes in the reflectivity parameter ξ r will affect the maximum displacement distance of the equilibrium points.
In Fig. 1, we plot the position of the equilibrium points related to L 1 and L 2 forβ = 5 ξ r = 0.85 for different sail orientations. Here, Fig. 1 (left), corresponds to α = 0, δ = 0 and all points are in the Z = 0 plane, whereas in Fig. 1 (right), corresponds to α = 0, δ = 0 and all the points are in the Y = 0 plane. In both plots, the arrows represent the direction of the acceleration given by the solar sail at each equilibrium point. We can see how, for a solar sail almost perpendicular to the Sun-sail line, the equilibria related to L 2 are closer to the asteroid than those related to L 1 . Hence, in the case of monitoring an asteroid we must place a solar sail close to the displaced L 2 , as L 1 is out of the Hill region. In Fig. 2, we plot the position of equilibrium points related to L 1 and L 2 forβ = 5 and ξ r = 0.0, 0.30, 0.75, 1.0, where we can see the effect of varying the reflectivity parameter of the sail. As before, Fig. 2 (left), corresponds to α = 0, δ = 0 and Fig. 2 (right) to α = 0, δ = 0. Notice that for ξ r = 0 (i.e. only absorption) the acceleration given by the sail is only on the r s direction and we cannot displace equilibria away from the x-axis. As we increase ξ r we increase this displacement reaching its maximum for ξ r = 1, a perfectly reflecting sail.
When we look at the stability of the equilibrium points we have that the points on the L 1 family have different stability properties, where the linear dynamics can be saddle × centre × centre, saddle × saddle × centre, centre × centre × centre or even complex-saddle × centre. The points on the L 2 family, however, are all of type saddle × centre × centre. In this paper, we will focus on the nonlinear dynamics around the points of type saddle × centre × centre, where the eigenvalues are ±λ, ±iω h and ±iω v .
In Fig. 3, we plot the variation of λ and (ω h − ω v ) for the family of equilibria related to L 2 . Recall that the equilibrium points in Fig. 1 left are parameterised by α and those in Fig. 1 right by δ. As we can see, the equilibrium points are more unstable when α, δ ≈ 0, because their position is closer to the asteroid. On the other hand, notice that for the equilibrium points on the Z = 0 plane ( Fig. 3 right), the two centre oscillations undergo a 1:1 resonance (ω h = ω v ). This has an implications on the periodic and quasi-periodic motion around them that we will explain in more detail in the following sections.
In Fig. 4, we plot the variation of λ and ω h − ω v for the equilibrium points related to L 1 that are of type saddle × centre × centre. Notice that now for α, δ ≈ 0 the equilibrium points are less unstable than those with α, δ = ±π/2 (i.e. no solar sail). This is expected as the equilibria are further away from the asteroid. Moreover, we can see that the equilibria in the ecliptic plane (Fig. 4 right) also have values ofβ for which we have 1:1 resonances between the two centre oscillations.
In the case of a mission scenario in which we want to monitor a small body from a spacecraft with a solar sail, we are interested on equilibrium points that are relatively close to the asteroid and placed on the day side. This brings to a dilemma, as the equilibria that are the closest to the surface of the asteroid are the ones belonging to the L 2 family, which is placed on the dark side. If we consider the equilibria that are on the asteroid's day side (belonging to the L 1 family), these points quickly go away from the asteroid's vicinity asβ increases. Actually for realistic values for a solar sail (Farrés et al. 2014a) the displaced L 1 equilibria are out of the Hill sphere of influence and the AH3BP is no longer a good approximation to describe the dynamics around them. On the other hand, previous works (Broschart et al. 2014) show that we can find, around the displaced L 2 equilibria, periodic and quasi-periodic orbits that spend more than half of the time on the day side, which would be useful in a practical mission. For these reasons, we only focus on the dynamics around the equilibrium points belonging to the L 2 family.

Periodic and quasi-periodic dynamics around equilibria
This section provides results on the computation of periodic and quasi-periodic orbits around the L 2 displaced equilibrium points of the AH3BP for different sail orientations. As we have seen in the previous section, the linear dynamics for the displaced equilibrium points in the L 2 family is always centre × centre × saddle. This means that, under generic conditions, two families of periodic orbits emanate from each equilibrium point [Lyapunov's Centre Theorem (Siegel and Moser 1971;Meyer and Hall 1992)]. Each of these families of periodic orbits is related to one of the two pairs of complex eigenvalues (±iω h , ±iω v ) and the period of the orbits tends to 2π/ω h , 2π/ω v when we approach the equilibrium point. Furthermore, the coupling between the two frequencies gives rise to a Cantorian family of invariant tori, or Lissajous-type orbits (Jorba and Villanueva 1997).
We have divided this section in four subsections. In the first one, we describe the nonlinear dynamics for α = 0, δ = 0, i.e. when the solar sail is perpendicular to the Sun-sail line. In the second, we focus on the case when α = 0 and δ = 0 where all the equilibrium points belong to the Z = 0 plane. In the third, we focus on α = 0 and δ = 0 where the equilibrium points are on the Y = 0 plane. Finally, we perform a partial numerical unfolding of the neighbourhood of the 1:1 resonance when we vary the sail orientation. For all four cases, we have computed the different families of periodic orbits that emanate from the equilibrium points and their bifurcations, as well as the invariant tori around them, using both the centre manifold reduction (for energies close to the equilibrium point) and the computation of the invariant tori (for energies far from equilibria).

Dynamics for
This section provides results on the continuation of families of periodic orbits and reduction to the centre manifold when the solar sail is perpendicular, i.e. α = δ = 0. Although this case has already been discussed in detail (Farrés et al. 2014a, b), we use this section to describe the main invariant objects around the equilibria. It will also serve us as a reference to see how changes in the sail orientation affect the dynamics.

Families of periodic orbits
If we orient the solar sail perpendicular to the Sun-sail line, we are essentially changing the gravitational pull of the Sun on the probe. As it happens in the RTBP (Farrés and Jorba 2010b), the nonlinear dynamics close to the displaced equilibrium point is very similar to   (1) around the equilibrium point, SL 2 = (x * , 0, 0), the motion on the Z = 0 plane is decoupled from the rest. Hence one of the two families that emanates from the equilibrium point is completely contained in the plane. We denote it as planar Lyapunov family. The other family oscillates above and below the plane. We denote it as vertical Lyapunov family. As it happens in the classical Hill problem (i.e. no sail) (Gómez et al. 2005), when we move along the family of planar orbits, at some point we find a 1:1 resonance between the frequency of the orbit and its normal frequency and the halo family of periodic orbits appears.
On the left-hand side of Fig. 5, we show the bifurcation diagram for the continuation of these families of periodic orbits. The two axis represent the X and Z coordinates of periodic orbits on the Poincaré section 1 = {Y = 0,Ẏ > 0}, where each point corresponds to one periodic orbit and the colour refers to the stability of the orbits: green stands for centre × saddle, blue for saddle × saddle and purple for centre × centre. As we can see, the planar and vertical periodic orbits, born from the equilibrium point, are both centre × saddle. If we move along the planar family, at some point the stability changes, a pitchfork bifurcation takes place, and two families of halo orbits appear. The orbits on the planar family become saddle × saddle while the orbits in halo family are centre × saddle. On the other hand, the vertical family does not experience any bifurcation: the orbits in this family start by increasing their altitude (maximum Z displacement) up to Z ≈ 0.4 and then decrease until the orbit is almost planar.
In the other two plots in Fig. 5, we show the evolution of the stability parameters s 1 , s 2 as a function of the energy H . Here, each colour represents the family: the planar family is purple, the vertical family is green and the halo family is blue. Notice how the vertical family has one indicator greater than 2 and the other less than 2 throughout all the family. The planar family starts the same way but one of the indicators increases eventually crossing the value 2. At that moment, the two symmetric families halo orbits are born, both with the same stability properties. At some point, the halo orbits undergo an additional change of stability from centre × saddle to centre × centre.
In Fig. 6, we show the XYZ projections of the periodic orbits in these three families. From left to right, we have the planar, vertical and halo family. Notice how the orbits in the planar family start being almost circular and, as we move along the family, they elongate following a boomerang shape with the asteroid close to the central part. On the other hand, the vertical family follows the typical eight shape, with a larger vertical oscillation than horizontal. These orbits become almost planar at the end of the family. Finally, the two families of halo orbits are symmetric to each other. Towards the end of the family, as these orbits approach the asteroid, they are almost circular and become smaller.

Centre manifold
To have a more global description of the dynamics we use the Hamiltonian reduced to the centre manifold. As it has been discussed before, the linearised dynamics at the equilibrium point can be described as centre × centre × saddle. The saddle is responsible for the instability of the neighbourhood of the equilibrium point: almost all orbits inside this neighbourhood go away after a short integration time. It is also known that, due to the Hamiltonian character of the system, each linear centre (elliptic) direction is tangent to the Lyapunov family of periodic orbits, as described in the previous section. The combination of the two linear oscillations gives rise to a family of quasi-periodic motions, sometimes known as Lissajous orbits, with two basic frequencies. Near the equilibrium point, these quasi-periodic motions can be seen as the nonlinear coupling of the families of Lyapunov periodic orbits (Jorba and Villanueva 1997).
After computing the power expansion (up to degree 32) of the Hamiltonian restricted to the centre manifold, we have applied the procedure described in Sect. 3.2. We denote as (q h , p h ) the position-momentum pair of coordinates of the Hamiltonian reduced to the centre manifold, that are related to the horizontal oscillations. Let us select the Poincaré section q v = 0, which is transversal to the vertical oscillations. We start by selecting an energy level H cm = h * . We select a mesh of values (q h , p h ) and we numerically compute the (unique) value of p v such that H cm (q h , p h , 0, p v ) = h * (if such a value does not exist, it implies that there is no orbit in this energy level crossing the section q v = 0 at the given point (q h , p h )). Then, the successive intersections of these trajectories with the Poincaré section q v = 0 are plotted.
In Fig. 7 left, we show this section for energy H cm = 0.4 (which corresponds to H = −4.51907174 in the initial synodical coordinates). The fixed point in the middle of the plot corresponds to the vertical Lyapunov periodic orbit in this energy level, and the invariant curves around this point represent quasi-periodic motions (also called Lissajous orbits). The outer boundary of this plot is the planar Lyapunov orbit. If we increase the energy up to H cm = 0.8 and H cm = 2 (which corresponds to H = −4.45085751 and H = −4.24621480 for the initial Hamiltonian) we obtain Fig. 7 (centre and right). Note the two halo orbits and how they are surrounded by quasi-periodic orbits, the so-called quasi-halos. We can relate the phase portraits of Fig. 7 to the evolution of the families of periodic orbits in Fig. 5 (right), where we can check that for H = −4.51907174 there are no halo orbits, while for H = −4.45085751 there are. Observe that the Poincaré section {q v = 0} is not a plane in synodic coordinates. In summary, for α = δ = 0, the dynamics near the equilibrium point is similar to the one of the RTBP.

Dynamics for α = 0, δ = 0
This section provides results on the continuation of periodic orbits, reduction to the centre manifold and numerical continuation of invariant tori for α = 0, δ = 0. There is also a discussion on how the 1:1 resonance that occurs for α = α crit and δ = 0 affects the dynamics.

Families of periodic orbits
We recall that if we select sail orientations with α = 0, δ = 0 we displace the equilibrium points left or right with respect to the Sun-asteroid line. The same happens with the families of periodic orbits. It is easy to check that, in this case, the Z = 0 plane is invariant. For small variations of α, we have the same qualitative behaviour as for α = 0, but as we increase the sail orientation the families of orbits will experience some changes. It is found that, for α ≈ 0.24, the planar Lyapunov family experiences a second change of stability going from saddle × saddle back to centre × saddle (the first was produced by the halo orbits). This is a pitchfork bifurcation that produces two new families of periodic orbits that are born at the planar Lyapunov family and die at the vertical Lyapunov family (again in a pitchfork bifurcation). We call these orbits Sideway orbits (as we will see from their shape later on). The Sideway orbits already appear in the classical Hill problem (Gómez et al. 2005) (no sail) and in the classical RTBP (Gómez and Mondelo 2001). For a perpendicular solar sail, the location of this second bifurcation point varies (Ceccaroni et al. 2016a). In the case of L 2 , as the sail lightness number, β, increases the second bifurcation point moves towards the primary and disappears at some point.
As we mentioned in Sect. 4, for some sail performance parameters (β, ξ r ), there is a value of α where we have a 1:1 resonance between the planar and vertical oscillations of the linearised system at the equilibrium point. For β = 5, ξ r = 0.85 and δ = 0, this happens for α = α crit ≈ 0.50781958553993878. We have observed that as we get closer to the 1:1 resonance (α < α crit ) the halo and Sideway orbits bifurcate closer to the equilibria, and eventually collapse at the equilibrium for α = α crit . After crossing the 1:1 resonance (α > α crit ), the halo orbits bifurcate from the vertical Lyapunov family instead of the planar one (the two families have changed their role with respect to the halo orbits) and the Sideway orbits disappear.
On the left-hand side of Fig. 8, we plot the bifurcation diagram for α = 0.1 (top) and α = 0.26 (bottom). Again we have the X Z projection of the orbits on the Poincaré section   The other two plots on Fig. 8 show the variation of the stability parameters as a function of the energy H . There we can see the different stability changes and bifurcations. Notice how for α = 0.26 we have that both the planar and vertical families change their stability close to H = 0.6 and the family of Sideway orbits connecting them appears.
In Figs. 9, 10 and 11 we have the XY Z projections of the planar Lyapunov, vertical Lyapunov and halo families of periodic orbits, respectively, for different sail orientations α = 0.26, 0.44 and 0.52. As we can see, taking α = 0 does not affect much to the shape of the orbits, it just displaces them with respect to the Sun-asteroid line. In Fig. 12, we find the two family of Sideway periodic orbits for α = 0.26, 0.44 and 0.50, where we can appreciate how they change from the shape of a planar Lyapunov orbit into the shape of a vertical Lyapunov by gaining vertical amplitude and loosing horizontal one. We also observe how these orbits become smaller as α increases and they bifurcate closer to the equilibrium point.

Centre manifold
As it has been done in the previous section, to show a more global picture of the dynamics we use reduction to the centre manifold. For instance, if we consider the sail orientation α = 0.26 and δ = 0 (see Fig. 8 bottom for the periodic orbit diagram) we see that, near the new equilibrium point, the dynamics is very similar to the previous case (α = δ = 0), and as we increase the energy level, the halo orbits appear. Figure 13 shows the dynamics in the centre manifold for several values of the Hamiltonian, where we can see that there is a similar dynamics as for α = δ = 0, the main difference being that the orbits are shifted in the plane. It is noticeable that, when the energy grows, the quasi-periodic motions around halo orbits occupy a very large portion of phase space.

Invariant tori
The iso-energetic Poincaré sections obtained through reduction to the centre manifold can be extended to higher energy levels through direct numerical continuation of families of invariant tori, as described in Sect The tori computed from items (a) and (b) are of the same family, and are known as Lissajous tori. The tori from item (c) are known as quasi-halo.
The invariant tori obtained are represented in Fig. 14, in terms of the h, ρ parameters described in Sect. 3.3. The points in this figure do not correspond to the tori obtained by numerical continuation, but to tori that have been refined from them to a grid of equally spaced energy values, in order to be able to produce the Poincaré sections that will be shown later. The ranges of central part of the families of periodic orbits used to start the continuation of tori are represented as curves in Fig. 14. For all these curves except the one representing the planar family in item (a) above, ρ is chosen so that the central eigenvalues of the monodromy matrix are λ j , λ −1 j with λ j = e iρ . The relation between ρ and the stability parameter is then For the curve representing the planar family, however, the relation between the stability parameter and ρ cannot be the one of Eq. (10), because then the curves representing the planar and vertical families would not meet at the energy of the equilibrium point. Lyapunov's Centre Theorem (Siegel and Moser 1971;Meyer and Hall 1992) instead of Eq. (10). From all these choices, it follows that, for all the tori represented in Fig. 14, the relation between ρ and the frequencies of the tori are where, of all the possible ways to choose the frequencies of a torus 1 , ω h and ω v have been chosen as to approach the horizontal and vertical frequencies, respectively, of the equilibrium point when the torus is close to it. For both the tori represented in Fig. 14a and the ones represented in Fig. 14b, the continuations have been stopped when either the tori have collapsed to a periodic orbit or a computational limit of N f = 100 (see Eq. 7) and a maximum error of 10 −8 has been reached. Except for the families that reach this computational limit, in Fig. 14a all the continuations finish by collapsing to a vertical Lyapunov periodic orbits. Because of this, the region delimited by the curves given by the planar and vertical families represents the whole Lissajous family of invariant tori. In Fig. 14b, however, all the continuations finish because of the previous computational limit. It is an open question which is the dynamical limit of these families. Figure 15 shows iso-energetic Poincaré sections of the invariant tori represented in Fig. 14. Each curve in any of the plots of Fig. 15 is the intersection with {Z = 0,Ż > 0} of an invariant torus of the corresponding energy level represented by a point in either Fig. 14 a or b. These Poincaré sections are analogous to the ones of Fig. 13, except for the fact that Fig. 13 is in centre manifold coordinates, whereas Fig. 15 is in original (synodic) coordinates. We can observe how, as the energy level increases, quasi-halo tori occupy most of configuration space and Lissajous tori are gradually squeezed by them.    A remarkable characteristic of the AH3BP is that the bifurcation point where the halo orbits are born approaches the equilibrium point as α increases. This means that the energy level for which the halo orbits are born approaches the energy level of the point. In Fig. 16 we show the variation of the stability parameters with respect to the energy for α = 0.3, 0.4 and 0.44 (before α crit ). There we have zoomed the stability diagram in order to appreciate how both the halo and the Sideway orbits bifurcate closer to the origin as α increases. When α = α crit ≈ 0.507819585540, the halo orbits are born from the equilibrium point and, for values of α slightly larger than α crit , the halo orbits are born from the vertical family of periodic orbits. On Fig. 17 (left) we find the bifurcation diagram on the Poincaré section for α = 0.50 (top) and α = 0.52 (bottom), before and after the 1:1 resonance. As we can see, before the 1:1 resonance (α = 0.50) the halo orbits bifurcate from the planar family of orbits and after they bifurcate from the vertical family. Moreover, as we have already mentioned, the Sideway orbits disappear for α > α crit , and therefore have already disappeared for α = 0.52. This is better appreciated in the zoomed stability parameter plots of Fig. 17 right. The evolution across the 1:1 resonance of the quasi-periodic motion that surrounds the periodic one can be observed by either numerical continuation of invariant tori or reduction to the centre manifold. The first technique is the only choice when the bifurcations of the planar Lyapunov family giving rise to the halo and Sideway orbits takes place at energies far from the one of the equilibrium point. As α approaches α crit , these two bifurcations approach the equilibrium point, and the transition can be observed by reduction to the centre manifold.
As an example of the first case, we have performed numerical continuation of invariant tori for α = 0.40, δ = 0. The corresponding energy-rotation number diagrams are shown in Fig. 18. In addition to the families of tori described for α = 0.26, now we have also continued invariant tori starting from the central part of the planar family of periodic orbits for a small range of energies after its second bifurcation, in which the Sideway family of periodic orbits is born (see Fig. 16, centre). This new family of tori, that we denote as quasi-planar, is represented in Fig. 18 c. Corresponding iso-energetic {Z = 0,Ż > 0} Poincaré sections are shown in Fig.19. The different energy levels of Fig. 19 can be related to Fig. 18 as follows. In Fig. 19a only Lissajous tori are present, surrounding the planar periodic orbit of its energy level. In Fig. 19b, quasi-halo tori appear, surrounding the halo orbits of its energy level. In Fig. 19c, quasi-halo tori grow and squeeze the Lissajous family. In Fig. 19d, we can see both quasi-planar tori and the two Sideway periodic orbits of this energy level, whose location is outlined by the quasi-planar and the Lissajous family of tori. The Lissajous family becomes smaller in Fig. 19e, to completely disappear in Fig. 19f, once the vertical family has lost its central part. The location of the vertical periodic orbit in this last plot is outlined both by quasi-halo and quasi-planar tori.
This transition can also be seen in Figs. 20 and 21, obtained from the reduction to the centre manifold and using suitable energy levels for α = 0.50 and α = 0.52, respectively. As we have just mentioned, all these changes appear close to the equilibrium points, hence the centre manifold is able to capture all the dynamics. For this reason, we have not performed an explicit computation of invariant tori to describe this phenomena, as it is a more tedious work and would not provide any extra information. Figure 20 shows the Poincaré section for α = 0.50 and H cm = 0.45, 0.49 and 0.52, respectively (left to right). In all the plots, the equilibrium point in the middle corresponds to the vertical Lyapunov orbit for this energy level and the two equilibrium points close to p h = 0 correspond to the two halo orbits. The motion is bounded by the planar Lyapunov orbit. In the middle plot (H cm = 0.49), we have two more equilibrium point close to the vertical axis (q h = 0) that are unstable (notice the saddle point behaviour). They correspond to the two Sideway orbits (see Fig. 17). As we increase the energy, these two points approach the vertical Lyapunov orbit (middle equilibrium point) and we have a pitchfork bifurcation. For H cm = 0.52 (right plot), we can see that now the vertical Lyapunov orbit is unstable. This behaviour is analogous to the one observed in Fig. 19.
In Figure 21, we show the same analysis for α = 0.52, having considered H cm = 0.01, 0.02 and 0.20 (respectively). Again, in all the plots the equilibrium point in the middle corresponds to the vertical Lyapunov orbits. The two other equilibrium points correspond to the halo orbits, and the motion is bounded by the planar Lyapunov orbit. Notice that for H cm = 0.01 the vertical Lyapunov orbit is stable and for H cm = 0.02 this one is unstable (saddle motion around the equilibrium point). We can also see how the two halo orbits are born from the vertical family after the pitchfork bifurcation. This behaviour remains as we increase the energy, showing no sign of the Sideway orbits. We can relate these picture with the structure of the periodic orbit families that we find in Fig. 17.
In Sect. 5.5, we will describe what happens in the neighbourhood of the 1:1 resonance for δ = 0. A more complete description will be the goal of a future work.

Schematic road-map for the periodic orbits behaviour
As we have seen, the behaviour and configuration of periodic orbits as we move one of the angles defining sail orientation is very rich. This section summarises graphically the different type of orbits that we have for a fixed set of values of the sail parameter and describe what happens when we move α and cross the 1:1 resonance.
For α < α crit , there are four different type of families that play a role on the dynamics close to the equilibrium point. There is the classical planar Lyapunov and vertical Lyapunov orbits that are born at the equilibrium point. Two symmetric halo orbits and two symmetric Sideway orbits [also know as axial-symmetric orbits (Gómez et al. 2005)]. The halo and Sideway orbits appear at two different Pitchfork bifurcation in the Planar Lyapunov family. The first Pitchfork bifurcation gives rise to the halo orbits, which evolve and undergo their own bifurcations, and the second Pitchfork bifurcation gives rise to the Sideway orbits, which connect with the Vertical Lyapunov orbits. In Fig. 22, we have a schematic representation of this bifurcation diagram, which is also observed in the classical Hill and RTBP problem for certain mass parameters (Gómez et al. 2005;Gómez and Mondelo 2001;Ceccaroni et al. 2016a).
As we increase α, the Pitchfork bifurcation that gives rise to the halo family comes closer to the equilibrium point. Moreover, the two Pitchfork bifurcation producing the Sideway family (one bifurcation in the planar family and another one in the vertical family) also come closer to the point so that the full family disappears (merging into the equilibrium point) when α reaches α crit . For α = α crit there is a 1:1 resonance between the planar and vertical frequencies (ω h , ω v ). For this value of α, four families of orbits are born from the equilibrium point (one planar and one vertical Lyapunov and the two halo). For α > α crit the Planar Lyapunov family does not experience any Pitchfork bifurcations, but now the halo orbits are born from a Pitchfork bifurcation of the Vertical Lyapunov family, and there are no Sideways orbits. A schematic representation of this transition can be seen in Fig. 23. This section analyses the periodic and quasi-periodic dynamics around the displaced L 2 equilibrium points for α = 0, δ = 0. We provide results on the continuation of families of periodic orbits and on the reduction to the centre manifold.

Families of periodic orbits
We recall that, if we consider the sail orientations with α = 0, δ = 0, we displace the equilibrium point above and below the horizontal plane Z = 0. As we will see here, the effect of moving the sail up and down with respect to the Sun-sail line has a more drastic effect on the distribution of the periodic orbits. The reason is that the symmetry responsible for the existence of the several pitchfork bifurcations seen before disappears when δ = 0. In particular, the horizontal plane Z =Ż = 0 is no longer invariant. Hence, from now on the planar family of periodic orbits is no longer contained in a plane. Now we have two families of periodic orbits that are born from the equilibrium points, where the horizontal family starts almost planar but it quickly gains inclination looking as halo-type orbits. On the other hand, the vertical family of periodic orbits shows a similar behaviour as for δ = 0. In both cases, the planar and vertical family do not experience any 1:1 bifurcations.
If we look at this configuration, we could think that orbits that used to be there for δ = 0 now have disappeared. This is not true, as there are still planar orbits and the other halo branch, but they are related in a different way as we can see in Fig. 24. What happens is that the pitchfork bifurcation that we have for δ = 0 is now broken, due to the symmetry breaking induced by δ = 0, and now we have a curve with no change in stability and a saddle-node bifurcation. The same phenomena is observed in the RTBP when we include the effect of a solar sail (Farrés and Jorba 2010b). We can see the phenomena explained above: we see how the two families of periodic orbits that are born from the equilibrium point (rightmost point in the plots) evolve, one of them recalling the behaviour of a halo orbit and the other the behaviour of the vertical orbit. There we also see that we still have the second branch of the family of halo orbits, but they are no longer symmetric one with respect to the other. Finally, we also see that, close to the asteroid, we still have planar periodic orbits, that are related to one of the halo branches through a saddle-node bifurcation. Here, we clearly see how the pitchfork bifurcation has been broken. Also, the gap between the saddle-node bifurcation and the equilibrium points increases as we increase δ as expected.
The other plots in Fig. 25 show the variation of the stability with respect to the energy H of the system. Here, each colour corresponds to one type of orbits: purple for the horizontal family, green for the vertical family and blue for the remaining branch that contains planar and halo-type orbits. Here, we can also appreciate that the halo-type orbits are no longer symmetric and that planar and halo orbits are related in a different way. We can also see that as δ increases, the stability parameters of the planar and vertical orbits separate from s i = 2, suggesting that there will not be a second bifurcation giving place to Sideway orbits.
In Fig. 26, we have the XY Z projection of the vertical family of periodic orbits for δ = 0.02, 0.26 and 0.52, respectively. There we can see that the non-symmetric character of the system has also an impact on the shape of the orbits. The vertical family for δ = 0 has a symmetric eight shape, but now one loop is larger than the other. In Fig. 27, we have the XY Z projections of the horizontal family of orbits emanating from equilibria. Notice that their shape recalls a halo orbit rather than a planar one. In Fig. 28, we have the other branch of the family, where we have halo and planar orbits that are now connected. We note that, if we plotted both families together, we would see that we essentially have the same type of orbits as for δ = 0, slightly shifted above the ecliptic plane and related to each other in a different way due to the symmetry breaking that we have already mentioned.

Centre manifold
To have a more complete description, let us give a look at the dynamics in the centre manifold for these cases. We start with the values α = 0 and δ = 0.02, Fig. 29 Fig. 7 for the case α = 0, δ = 0, where we have a equilibrium point in the middle corresponding to the vertical Lyapunov orbit. We must mention that, as δ = 0, the planar orbit is not contained in the {q v = 0} plane (in fact, it intersects transversely this plane) and, hence, it is no longer bounding the motion in the centre manifold. When we look at larger values of the energy, some bifurcations appear. For instance, in Fig. 29 (middle) we can see in the centre the vertical Lyapunov orbit, and in the left-hand side a small island corresponding to a halo orbit. This halo orbit is the one that appears when we follow the planar family that is born from the equilibrium point (see Fig. 25).
If we keep on increasing the energy value we see how a saddle-node bifurcation (coming from the symmetry breaking of the pitchfork bifurcation that gives rise to the halo orbits) takes place and the planar and the other halo orbit appear. We can see them in Fig. 29 (right) for H cm = 1.20. Here, the elliptic point in the centre corresponds to the vertical Lyapunov orbit, the two other elliptic points (close to p h = 0) correspond to the two halo orbits (we can appreciate the asymmetry between these two orbits). Finally, the unstable equilibrium point (saddle behaviour around it) that is at the right hand side of the section corresponds to the planar orbit that, as we have mentioned before, intersects the Poincaré section.

The neighbourhood of the 1:1 resonance
In this section, we focus on the dynamics in a neighbourhood of the displaced L 2 equilibrium point for α = α crit , δ = 0. As we have mentioned before, for these values of the parameters Here, we perform a small numerical study of the dynamics near this bifurcation using the natural parameters α and δ. These two parameters are not enough to fully unfold this bifurcation. A complete unfolding is at present work in progress. We have already described the dynamics near the equilibrium point when δ = 0 and α goes through α crit , showing that the bifurcation that gives rise to the halo orbits goes from the planar to the vertical Lyapunov families. Besides, the family of (Sideway) periodic orbits that connects the two families of Lyapunov periodic orbits collapses to the equilibrium point when α reaches α crit and disappears when α ≥ α crit , as shown in Fig. 30, which is a magnification of the central region of the plots in Fig. 20. The main effect of the parameter δ is to break the symmetry responsible for these pitchfork bifurcations so that they all become saddle-node bifurcations. To show it with more detail, we select a small value for δ (for instance, 10 −4 ) and describe the neighbourhood of the equilibrium point when α approaches α crit . Figure 31 shows the effect of δ = 0 on the two pitchfork bifurcations displayed in Fig. 30. The first row of Fig. 31 refers to the transition from the left to the centre plot in Fig. 30, and the second row corresponds to the transition from the centre to the right plot in Fig. 30. Let us start by describing the first row of Fig. 31, where the saddle node appears near the boundary of the Poincaré section: the planar Lyapunov orbit is hyperbolic and it moves up (from the bottom part of the figure) while a saddle-node bifurcation appears in the upper part of the figure. As a result, there is an elliptic orbit at the centre of the plot and two hyperbolic orbits below and above it as seen in the last plot of the first row in Fig. 31 (compare this plot with the centre plot in Fig. 30). In the second row of Fig. 31 the initial central elliptic periodic orbit moves upwards to meet an hyperbolic orbit and then disappears (this is a saddle-node bifurcation), while the lower hyperbolic orbit moves to the centre of the image (compare the last plot of this second row with the last plot in Fig. 30). We note that the use of other parameters to unfold the 1:1 resonance at α = α crit , δ = 0 could result in a different unfolding. In fact, more than two parameters are needed to fully unfold this bifurcation (which, as mentioned before, is work in progress).

A note on possible mission applications
In the mission scenario of a solar sail close to an asteroid, there are several constraints that must be taken into account in order to select a target orbit. For instance: orbiting close to the asteroid to allow detailed scientific observations; keeping the solar sail almost perpendicular to the SRP to magnify the sail's performance; and being able to observe the day side of the asteroid. As seen in Sect. 4, we must focus on the periodic and quasi-periodic orbits around the displaced L 2 equilibria to remain close to the asteroid. The main drawback of this region is that most of the periodic and quasi-periodic motion remains in the shadow side of the asteroid. Nevertheless, there exist periodic orbits, specially on the planar and vertical families, which spend part of their orbital period on the day side of the asteroid, and are interesting for an observational mission (Farrés et al. 2014a). Similar behaviour has been observed in the terminator and quasi-terminator orbits (Broschart et al. 2013(Broschart et al. , 2014. For a solar sail perpendicular to the SRP (α = δ = 0), the planar and vertical Lyapunov orbits that are at the end of the family (by end we refer to those orbits with a larger energy level) spend more than half of their orbital period on the day side of the asteroid (Farrés et al. 2014a). As we have seen in Sect. 5 by changing the sail orientation we can displace these orbits: changes in α induce an extra force on the XY direction and the orbits are displaced to one side or the other of the Sun-asteroid line; and changes in δ induce an extra force on the Z direction and displaces the orbits above and below the asteroid's orbital plane. In both cases, we will change the symmetric shape of the orbits. The two plots on the left of Fig. 32 show, for a given planar and vertical Lyapunov orbits, their variation with respect to changes in α (top) and changes in δ (bottom). The orbits in purple represent α, δ > 0 and those in green are for α, δ < 0. We can see that these orbits spend most of their period on the day side of the asteroid (i.e. (X < 0)). These orbits are linearly unstable and an active control law must be implemented in order to remain close to them.
As we move along the two halo families, the orbits come close to the asteroid but remain always on the dark side of the asteroid. Again, changing the sail orientation we can displace their position, having orbits that spend part of their orbital period on the day side of the asteroid. On the right hand side of Fig. 32 we see the variation of one of the two halo orbits for changes in α (top) and changes in δ (bottom). The other halo orbit experiences similar effects. We recall that the orbits at the end of the halo family (also known as Terminator orbits) are linearly stable and no station keeping is required. Being able to displace the periodic orbits with a solar sail is very interesting, mainly because it increases the visibility properties of the periodic orbits and their potential for observational missions. Moreover, an active control law can be used to drift along these families of displaced orbits following the ideas by Jorba (2008b, 2016a).

Conclusions
We have provided a fairly complete description of the bounded dynamics around the equilibria of the Augmented Hill 3-Body Problem, which models the motion of a spacecraft around a point-mass asteroid, equipped with a non-perfectly reflecting solar sail, under uniform SRP and uniform gravity from the Sun. For each fixed set of parameter values, this model has two equilibrium points equivalent to the ones of the classical Hill problem. By numerical continuation of the fixed point equations, we have analysed how these equilibrium points are displaced as the parameters are varied. Practical considerations concerning asteroid missions have led to the selection of one of the equilibrium points (for each set of parameter values) and a range of sail orientations, for which a description of their effect on the surrounding periodic and quasi-periodic dynamics has been made. This has been done through a variety of numerical methods: numerical continuation of periodic orbits, numerical integration of the vector field semi-analytically reduced to the centre manifold of the equilibrium point, and direct numerical continuation of invariant tori. We have stressed the range of validity of each methodology, and how they complement each other in order to produce global dynamical descriptions. This was one of the main motivations of the paper. We have recovered some of the known families of periodic orbits and invariant tori around the libration points of the RTBP and the classical Hill problem, but we have found that here they interact in different ways than in those problems due to the presence of the sail. In particular, we have found that several sail orientations give rise to a 1:1 resonance around the equilibrium point, which we have partially (numerically) unfolded. This is another main point of interest of the paper, that will be the basis of future work. The paper concludes with some comments on the usefulness of the trajectories obtained in possible asteroid missions.