Higher Order Approximation to the Hill Problem Dynamics about the Libration Points.^{†}^{†}thanks: Commun Nonlinear Sci Numer Simulat (2018) https://doi.org/10.1016/j.cnsns.2017.12.007 (Preprint version)
Abstract
An analytical solution to the Hill problem Hamiltonian expanded about the libration points has been obtained by means of perturbation techniques. In order to compute the higher orders of the perturbation solution that are needed to capture all the relevant periodic orbits originated from the libration points within a reasonable accuracy, the normalization is approached in complex variables. The validity of the solution extends to energy values considerably far away from that of the libration points and, therefore, can be used in the computation of Halo orbits as an alternative to the classical LindstedtPoincaré approach. Furthermore, the theory correctly predicts the existence of the twolane bridge of periodic orbits linking the families of planar and vertical Lyapunov orbits.
1 Introduction
The restricted threebody problem, in which a body of negligible mass evolves under the gravitational action of two massive bodies whose relative motion is Keplerian (see [1], for instance), is a useful approximation to the real dynamics of planetary satellites and minor bodies of the solar system, and is customarily used in the first steps of mission designing for some artificial satellites [2, 3, 4]. Of specific interest in astrodynamics is the use of trajectories related to the libration point dynamics [5, 6, 7, 8]. These points are relative equilibria in a frame rotating with the Keplerian rotation of the primaries.
Rather than pursuing the general solution of the dynamics about these points, the focus has been put on the computation of particular solutions, either numerically, as is the case of the computation of periodic orbits [9, 10, 11, 12, 13] as well as other invariant manifolds [14, 15], or analytically, a case in which the LindstedtPoincaré method is the usual approach [16, 17, 18, 19].
In the case of the collinear points, due to the unstable dynamics, a partial reduction to the center manifold is customarily done to remove the hyperbolic directions, which are a consequence of the saddle center center character of the linear dynamics. After reduction, the center manifold Hamiltonian is of 2degrees of freedom (DOF) and, therefore, can be approached with the usual tools of nonlinear dynamics [18, 15].
A common case in the solar system is that the mass ratio between different pairs of celestial bodies (planetsun, satelliteplanet, etc.) is very small. Then, the Hill problem formulation —a limiting case of the restricted threebody problem in which the two massive bodies evolve in circular orbits about the mutual center of mass, and the body of negligible mass evolves close to the less massive body when compared to the distance between the massive bodies [1]— releases the dependency of the dynamics on physical parameters. This fact provides a wide generality to this particular formulation, introduced by Hill in his seminal investigations of the Moon’s orbit [20], whose application encompasses a variety of problems, ranging from the interaction between particles in planetary rings [21] to the study of the dynamics of planetary satellites [22] and about them [23, 24], or the modeling of tethered satellites applications [25].
The Hill problem dynamics is studied with the same techniques as the RTBP [26, 27, 28, 29, 30, 31, 32]. Alternatively, it has been recently proposed to study the dynamics about the libration points from the reduction of the Hill problem Hamiltonian to an integrable one [33]. Indeed, because the libration points are equilibria of the saddle center center type, the Hamiltonian can be transformed into one such that, after truncation of higher order terms, the hyperbolic part of the quadratic Hamiltonian has been converted to an integral, and the reduced Hamiltonian, which contains the resonant terms of the elliptic part of the quadratic Hamiltonian, is of 1DOF. The reduced phase space is the sphere [34, 35] and is customarily described in the Hopf variables [36]. However, the solutions of the reduced phase space are more insightfully described in Lissajous variables [37], which allow the reconstruction of the periodic solutions by simply evaluating the equilibria solutions for each value of the elliptic anomaly between and , on the one hand, and ease the computation of the period of the periodic orbit, on the other.
We focus on the Hill problem and extend the normalization of [33] to higher orders so that it can provide acceptable approximations of the solution for orbits far away from the libration points. Celestial mechanicians, as well as physicists, traditionally normalize a system by removing cyclic variables. Thus, the normalization is routinely split into the preliminary reduction to the center manifold and the following removal of shortperiod effects [38, 39]. However, the second reduction requires to handle long Fourier series which may become unwieldy at relatively low orders of the perturbation theory. On the other hand, mathematicians have a long experience dealing with perturbed harmonic oscillators in complex variables [34, 40]. With this alternative, the Hamiltonian reduction can be done with a single normalization. Furthermore, normalization in complex variables becomes a simple exercise of polynomial algebra and provides very simple expressions which only require the arithmetic operations that can be carried out directly by the computer’s hardware, in this way speeding notably evaluation of the perturbation solution.
The computation of higher orders of the normalization in complex variables allowed us to obtain a single analytical solution which is valid for all the periodic orbits of the main families of periodic orbits of the Hill problem originated from the libration points. Namely, the families of vertical and planar Lyapunov orbits, the family of Halo orbits, which bifurcates from the family of planar Lyapunov orbits, and, notably, the two lane bridge linking the families of planar and vertical Lyapunov orbits. The later, which exists only for energy values much higher than those of the libration points, is only achieved when the perturbation solution reaches the 6th order, yet corresponding solutions are just rough approximations of the partner periodic orbits. Acceptable approximations to orbits of this family come out only from the 14th order of the perturbation solution.
Computing higher orders in a perturbation approach may be questioned in two ways. On the one hand, the model approached by perturbations is always a simplification of the real dynamics, so the order of the perturbation solution in which the neglected effects of the dynamics would be apparent must be discussed in each particular application. On the other hand, by reasons of simplicity and efficiency, higher orders are customarily approached with the Lie transforms method [41, 42] using floating point arithmetic [43]. Because of that, the propagation of the truncation errors in the successive orders of the Taylor series expansion may increase nonnegligibly the numerical errors due to the number representation in the computer, a fact that would make nonsense trying to increase the accuracy of the solution by extending the computations beyond a certain order. To mitigate this last issue the perturbation solution is alternatively approached in integer arithmetic. However, while this last approach avoids the accumulation of truncation errors and, therefore, allows to progress exactly, the drawback of using integer arithmetic is the increasing size of the integers to be handled, which grows from order to order of the perturbation solution and may become enormous at relatively moderate orders. As a consequence, the time and memory requirements of the computation of successive orders grow high, thus making the computation of very high orders unpractical. Comparison of the solutions obtained using both techniques helps in estimating the growth of truncation errors of the floating point arithmetic perturbation solution, and provides a way of extrapolating practical limits for the applicability of a such kinds of solutions.
When the normalization is carried out in floating point arithmetic, the use of complex variables also helps in estimating the accumulation of truncation errors at each order of the perturbation theory. Indeed, when coming back from complex to real variables some residual complex terms will remain in the normalized Hamiltonian, and the size of the higher of the coefficients affecting this residual terms can be taken as an indicator of the truncation errors accumulated in the computations.
The paper is organized as follows. First, in Sect. 2, the Hill problem Hamiltonian is directly derived from the Newtonian dynamics [44], the origin is translated to a libration point, and the resulting Hamiltonian is expanded in Legendre polynomials in order to present a perturbative arrangement. Next, the linearized dynamics of the problem is discussed in Sect. 3; while this part is mostly borrowed from [33] additional details are given to show that the linear transformation that decouples the linearized dynamics is not unique; still, we adhere to tradition and use the usual transformation [18] in our computations. It follows the description of the normal form computation in Sect. 4, which is computed exactly using integer arithmetic up to the order 11, and approximately using floating point arithmetic up to the order 20; this section provides estimates of the accumulation of the truncation errors in the floating point arithmetic case due to the physical size of the computer’s registers. Also in this section, the equilibria of the reduced phase space are briefly discussed to show their correspondence with well known periodic solutions of the Hill problem. Finally, a variety of tests are presented in Sect. 5 to illustrate the performance, as well as the limits, of the analytical solution.
2 Hill problem Hamiltonian about the libration points
Let and be two massive bodies, of masses and , respectively, which, under their mutual gravitational attraction, are evolving in circular orbits about the system’s center of mass with constant angular velocity . Then, the distance between and remains constant. Let be a massless body evolving under the gravitational actions of and , and let define the position of with respect to the center of mass of the system. Then, from Newton’s gravitational law,
(1) 
where is the usual time, is the vector joining with , of modulus , is the vector joining , with , , and is the gravitational constant.
We study the motion of relative to in a rotating frame with rotation rate defining the axis direction , the axis direction is defined by the direction from to , and the axis direction completes a direct orthogonal frame. Then, , where, from the definition of the center of mass
and, from the derivative of a vector in a rotating frame,
(2) 
where dots over vectors mean differentiation in the rotating frame. Therefore, from Eqs. (1) and (2),
(3) 
where
(4) 
and, due to the circular motion of , the centripetal acceleration is
(5) 
Now, with Hill, we assume that and . More precisely, we assume that and . Then, and
Hence, after neglecting higher order terms, Eq. (3) is rewritten as the differential equation of the Hill problem
(6) 
The dot product of and Eq. (6) can be integrated, to yield
After scaling units of length by and time by , Eq. (6) reads
(7)  
(8)  
(9) 
revealing that the Hill problem does not depend on any parameter. Besides, it is simple to check that two equilibria, the socalled libration points, exist at the positions , , where .
The Hill problem accepts Hamiltonian formulation. Indeed, the flow in Eqs. (7)–(9) can be derived from the Hamiltonian
(10) 
where , , and are the conjugate momenta to , , and , respectively.
Because of the known symmetries of the Hill problem (see [32], for instance), it is enough to study the dynamics about one of the libration points, say . To do that, the origin is translated to . Since a translation is a canonical transformation, after neglecting constant terms the transformed Hamiltonian reads
(11) 
where
(12) 
Now,
where and is the distance from to the libration point. Then, for values , the term can be replaced by the usual expansion in Legendre polynomials , yielding
(13) 
in which
(14) 
whereas the other terms of the Hamiltonian comprise monomials of the form
(15) 
where are numeric coefficients and () are nonnegative integers.
3 Linear dynamics about the libration points
For small displacements about the libration point we can neglect terms of higher order than , and hence the zeroth order term (14) of the Hamiltonian (13) is representative of the dynamics. The last term in Eq. (14) has the form of a simple harmonic oscillator with frequency . Therefore, in the linear approximation, the motion in and decouples from the rest of the flow and is made of harmonic oscillations. That is, the equilibria of the Hill problem are of the center type relative to the direction.
On the other hand, the coupled  motion results from the integration of a linear differential system with constant coefficients. Indeed, from Hamilton equations,
(16) 
where means transposition, and
(17) 
The general solution of Eq. (16) is
(18) 
where is a 4 4 matrix of arbitrary coefficients, and is the vector of characteristic exponents , , which are the eigenvalues of . It is simple to check that
(19) 
are real numbers, thus giving place to an hyperbolic component of Eq. (18), a saddle direction, whereas
(20) 
with , are pure imaginary numbers, resulting in an elliptic or centertype component of Eq. (18). These exponents define the wellknown saddle center center type of the libration points of the Hill problem.
Note that the quadratic Hamiltonian
(21) 
is in separate variables and enjoys the same dynamical behavior as Eq. (14). The Hamiltonian flow stemming from Eq. (21) is
(22) 
with
(23) 
So it emerges the question if a canonical transformation can be found such that it transforms into . The answer is in the affirmative, and, because the equations of motion are linear, the linear transformation
(24) 
can be computed by solving from the underdetermined linear system
(25) 
which is obtained by equating the right sides of Eqs. (16) and (22), and replacing Eq. (24).
The solution of Eq. (25) is expressed as a function of 4 arbitrary coefficients, say
(26) 
where
(27) 
for an invertible transformation, and hence , , while and cannot vanish at the same time .
Additionally, the transformation in Eq. (24) must be canonical. Because of the linear character of the transformation, it happens that its Jacobian is also . Therefore, the canonicity is expressed as , where is the symplectic matrix of dimension 4, yielding the two additional constraints
(28) 
which make in Eq. (27), and left undetermined only two coefficients, say and .
Among the different possibilities, and in view of the similitudes of columns 1 and 3 of the matrix , it seems natural to choose , . In view of Eq. (28), it yields
with
which leads to the usual transformation matrix
(29) 
based on the eigenvector decomposition of in Eq. (17), cf. [18] (see, also, [45]).
The discussion of the particular merits of the different transformations derived from the biparametric family defined by Eqs. (26) and (28) is not made here. However to further stress that the use of Eq. (29) is not a requirement to achieve the reduction of the quadratic Hamiltonian (14) to the separable form in Eq. (21), we write explicitly an alternative choice. Namely,
that yields
(30) 
where .
4 Normalization
The instability of the libration points due to their saddle component can be skipped by choosing suitable initial conditions. Indeed, due to the center center part, when the hyperbolic part of the quadratic Hamiltonian (21) vanishes the resulting motion will be periodic or quasiperiodic. This dynamics is not constrained to the linear approximation and can be extended to the whole Hamiltonian
(31) 
which is obtained after applying the transformation (24) to Hamiltonian (13). The procedure for removing the hyperbolic components of Eq. (31) is called the reduction to the center manifold, and consist of a partial normalization that, after truncation to some order, converts into an integral, which is followed by constraining the motion to the manifold [46].
A following removal of the shortperiod terms normalizes the Hamiltonian [33]. However, the normalization of Eq. (31) does not need to be split into two different parts and, on the contrary, can be achieved with a single transformation. This is the approach we take here.
First of all, in view of the frequencies and are very close, we introduce a detuning parameter [47]
(32) 
Then, the zeroth order term (21) is split into , and we rearrange the Hamiltonian (31) as
(33) 
with
(34)  
(35)  
(36) 
where the comprise monomials in the subindex variables of the type given in Eq. (15). Note that, because of the detuning, is no longer an homogeneous polynomial of degree 3. With this artifact, the center components of Eq. (34) have the form of an elliptic oscillator of frequency .
The normalization of the Hamiltonian (33) is more easily achieved in complex variables [45, 34, 40]. Thus, the transformation
(37) 
with , is applied first. Hence,
(38) 
while the remainder terms of the Hamiltonian stay as monomials of the type of Eq. (15) in the new realcomplex variables. The form of the zeroth order Hamiltonian (38) converts the normalization process in an elementary problem of polynomial algebra (see [33] for details).
4.1 The reduced dynamics
The normalization is implemented using the Lie transforms procedure [41, 42], and yields a Hamiltonian in new variables in which, after truncation to some desired order, the hyperbolic and elliptic components of the quadratic Hamiltonian, viz.
(39) 
become formal integrals. Therefore, the normalized Hamiltonian is of one degree of freedom, and hence integrable. Furthermore, to skip the hyperbolic instability, we choose initial conditions in this way constraining the dynamics to the center manifold . In consequence, the final, normalized Hamiltonian is of the form
(40) 
where terms are polynomials in the complex variables, which are no longer homogeneous because of the detuning made in Eqs. (34)–(35). The first terms of the normalized Hamiltonian (40) are printed below:
(41)  
Note that the integral is not easily identified in the summands of the normalized Hamiltonian (40), which seems to remain as a 2 degrees of freedom Hamiltonian in the complex variables. The use of Hopf variables [36], given by the transformation
(45) 
with the constraint
(46) 
definitely helps in disclosing the formal integral , as well as in describing the reduced phase space, which is the sphere [34, 35]. Indeed, Eqs. (41)–(4.1) are trivially expressed in Hopf variables by using the relations
(47) 
Then, Eq. (40) takes the form
from which we derive the Hamiltonian flow , , where curly brackets represent the Poisson bracket operator. We find
(48)  
(49)  
(50) 
Hence, points
(51) 
on the sphere are always equilibria; the plus sign corresponds to vertical Lyapunov orbits and the minus sign to planar Lyapunov orbits. Besides, those points
(52) 
such that Eq. (49) vanish, that is
(53) 
are also equilibria. This new equilibria stem from the point in a pitchfork bifurcation at the value given by the root , and correspond to Halo orbits. Finally, points
(54) 
on the sphere such that Eq. (50) vanish, that is
(55) 
are equilibria as well. Computation of the roots and of the equation shows that they stem from in a pitchfork bifurcation, and collapse into . These equilibria on the sphere correspond to the twolane bridge of periodic orbits linking planar and vertical Lyapunov orbits. Interested readers are referred to [33] for full details on the discussion of the reduced phase space as well as basic references on the topic.
4.2 Computational issues
The normalization is computed exactly by avoiding decimal expansions of the involved rational and irrational numbers. The irrational numbers , , and are handled formally, an their respective powers are simplified as mach as possible. However, the size of the integer numbers involved in the rational coefficients of the monomials grows from order to order, as can be observed in Eqs. (41)–(4.1), soon causing memory allocation to become a serious issue, with the consequent rapid increase of computing time. In this way, we only succeeded in extending the computations to the eleventh order, in which the rational coefficients may involve integer numbers of more than 100 digits. We hasten to say that we relied on commercial, general purpose, symbolic algebra tools in our computations; development of specific manipulators by experts could, of course, ease considerably the task [48, 49, 50].
On the other hand, the use of floating point arithmetic expedites computations notably, but at the cost of introducing truncation errors due to the physical length of the computer’s registers. The computation time still grows exponentially with the order of the theory, but at a lower rate. This fact is illustrated in Fig. 1, where it is shown that when using floating point arithmetic the computation time grows with the order roughly as , , where is the time spent into the computation of the second order terms of the normalized Hamiltonian and generating function, whereas in the case of exact computations using integer arithmetic it grows as , . In fact, the time employed in the exact computation of the order 11 of the perturbation solution was almost times longer than , whereas was only in the floating point case; that is, approximately 16 times faster.
Memory handling issues are definitely less severe with the floating point approach. However, the number of terms to be evaluated by consecutive higher orders of the perturbation solution grows roughly with the quartic power of the order, and soon becomes enormous, as shown in Fig. 2. Hence, we did not progress in our computations further than the order 20.
The propagation of the truncation errors when using floating point arithmetic can be studied with the help of interval arithmetic [51, 50]. An alternative way of estimating these errors is as follows. On the one hand, the transformation from complex to real variables, which is exact when avoiding decimal expansions, will produce some residual complex terms due to the floating point arithmetic —which, of course, must be neglected. The absolute value of the greatest of the coefficients affecting these residual terms is an indicator of the truncation errors accumulated in the computations. On the other hand, the size of the coefficients of the monomials generally grows with the consecutive higher orders of the perturbation solution. Then, the ratio between the greatest coefficient of the spurious, complex monomials and the greatest coefficient of the true, real monomials can be taken as an estimator of the truncation errors introduced by the computer’s arithmetic. This is illustrated in Fig. 3, where we see that the growth rate of the complex residuals is higher than that of the coefficients of the real monomials, and Fig. 4, where we note the sharp growth of the truncation errors estimated with our criterion when passed the order 15th.
Besides, in view of the already mentioned rapid growth of the number of terms to be evaluated and because the dynamics about the libration points is generally highly unstable, thus making orbit propagation quite sensitive to the initial conditions, we propose values of this indicator as a practical limit for the validity of the analytical solution. From Fig. 4, this value would correspond to an order of the perturbation solution between, say, 10 and 16, yet these high orders are only needed for computing orbits far away from the libration points.
5 Performance of the analytical solution
The reduction carried out by the normalization is more insightfully appreciated when using Lissajous canonical variables [37]. In theses variables
(56)  
(57)  
(58)  
(59) 
with
(60) 
showing that the normalization removed the angle , the conjugate variable to the momentum . Then, the only variables of the reduced phase space are and , and the orbits are ellipses whose size, shape, and orientation evolve slowly.
The periodic orbits of the original space (the Hill problem) are computed analytically as follows. First, we compute the equilibrium in the Hopf variables representing the desired periodic orbit (planar or vertical Lyapunov orbits, Halo orbits, or periodic orbits of the bridge linking planar and vertical Lyapunov orbits). Then, and are trivially obtained from Eqs. (56) and (59), respectively, while is computed unambiguously from Eqs. (57) and (58). The choice of any particular value allows for the following computation of the complex variables
(61)  
(62)  
(63)  
(64) 
Next, the Lie transformation computed for achieving the reduction provides corresponding complex prime variables, from which the subindex 1 Cartesian variables are recovered using Eq. (37). Finally, Eq. (24) —with the matrix given by traditional choice in Eq. (29) or any other choice from the family represented by Eq. (26) with the constraints in Eq. (28) that could have been used alternatively— will provide the initial conditions relative to the libration point. Repetition of the procedure for different values of will give the desired orbit without need of integrating these initial conditions.
We explore the performance of the analytical solutions by comparing orbits predicted by different orders of the perturbation solution with their partner periodic orbits of the Hill problem computed numerically. We do the comparisons in three different scenarios. In the first one, we constrain the energy to values close to the energy of the libration point, a case in which only Lyapunov orbits exist. In the second case we explore higher energies, for which Halo orbits also exist but remain close to the libration point. Finally, we focus on the range of energy values in which, besides the Lyapunov and Halo orbits, the orbits of the twolane bridge linking planar and vertical Lyapunov orbits exist. The later is a quite challenging case because of the large size of the orbits, and, until our knowledge, has never computed before analytically.
5.1
For the vertical Lyapunov orbit, the order 4 of the perturbation solution is enough to mimic the true periodic orbit at the precision of the graphics, as shown in the left plot of Fig. 5. However, when initial conditions provided by the analytical solution are propagated in the original, Hill problem dynamics, the periodicity error , which is defined as
(65) 
where stands for any of the coordinates in the original phase space, is only of the order of . The accuracy increases in a continuous way up to the order 11, for which the propagation in the original model of initial conditions taken from the analytical solution result in a periodicity error better than after a period .
In the case of the planar Lyapunov orbit, the order 3 of the perturbation solution suffices for suplying initial conditions that close the orbit at the precision of the graphics (right plot of Fig. 5). However, is only of the order of . The periodicity error improves with higher orders of the solution and is of the order of for the order 11. Only very slight improvements are achieved when the perturbation solution is truncated to higher orders, which become negligible further than the order 13th.
5.2
Due to the fact that the orbits are much larger in this case, the 7th order of the perturbation theory is required to match a vertical Lyapunov at the precision of the graphics (left plot of Fig. 6), but the periodicity error is only . This value improves for increasing orders of the perturbation theory up to the order 15th, where , and does not improve with higher orders. For the planar Lyapunov orbit we needed to use the 9th order truncation of the theory to achieve initial conditions in the original problem leading to a periodicity error of the order of one thousandth, which is enough to close the orbit at the precision of the graphics, as shown in the right plot of Fig. 6. Like in the case of the vertical Lyapunov orbit, the order 15th of the perturbation theory improves the periodicity up to , but no further improvements are found with higher orders of the solution.
At this value of the energy, Halo orbits already bifurcated from the family of planar Lyapunov orbits. The analytical solution predicts them correctly starting from the ninth order truncation of the analytical solution, as shown in Fig. 7, where the periodicity error is . Successive higher orders of the perturbation solution succeed in gradually improving periodicity, but only up to , which happens with the order 17th of the analytical solution. Note that, while Lyapunov orbits always have the same initial conditions in the reduced phase space, as given by Eq. (51), the location of Halo orbits on the sphere is given by Eq. (52) after solving Eq. (53), and, therefore, depends on the order of the perturbation solution used in each case.
5.3
For this large value of , the periodic orbits are so large that in different parts of the orbits the radius to the libration point falls out the convergence region of the Legendre polynomials expansion, which defines the validity of the model. In particular, the Halo orbit surrounds the central body, and its distance to the libration point is always longer than the Hill radius. Therefore, it is not expected that initial conditions obtained from the analytical theory can be improved by differential corrections to converge to a true Halo periodic orbit.
Contrary to the Halo case, only parts of the other periodic orbits remain out the region where the Legendre polynomials expansion converges, and the analytical solutions succeeds in providing reasonable approximations to the true Lyapunov orbits. Thus, the order 13th of the perturbation solution is able to capture an approximation of the vertical Lyapunov orbit with periodicity , yet no further improvements are found for higher orders of the perturbation solution. Analogously, initial conditions of a planar Lyapunov periodic orbit with periodicity of the order of are obtained with the order 16th of the perturbation solution. In both cases the initial conditions are easily improved with differential corrections, leading to the true periodic orbits.
Furthermore, the second bifurcation of the family of Lyapunov planar orbits has already happened at this value of , and the perturbation solution is effective in capturing an orbit of the twolane bridge that links planar and vertical Lyapunov orbits. Indeed, as shown in Fig. 8, a truncation to the order 14th of the perturbation solution provides an orbit with periodicity . While the periodicity error is not improved with higher orders of the analytical solution, initial conditions provide by the order 14th are easily improved by differential corrections to get the true periodic solution.
6 Conclusions
A higher order normalization of the Hill problem Hamiltonian centered at a libration point eases the computation of a single perturbation solution that captures the four main families of periodic orbits of the Hill problem originated from the libration points. Namely, the families of planar and vertical Lyapunov orbits, the family of Halo orbits, and the twolane bridge of periodic orbits that connects both families of Lyapunov orbits. Planar and vertical Lyapunov orbits exist for all energies above the energy of the libration points, and, therefore, close to the libration points these kinds of orbits are accurately reproduced with the lower orders of the analytical solution. Still, higher orders of the solution are required if one wants to obtain vertical and planar Lyapunov orbits far away from the libration points within an acceptably accuracy. The initial orbits of the Halo family require, at least, the 5th order truncation of the perturbation solution, whereas the orbits of the twolane bridge of periodic orbits that connect the families of planar and vertical Lyapunov orbits require, at least, a 14th order truncation of the analytical solution to obtain a reasonable approximation of the corresponding real periodic orbits, whose initial conditions can be improved by means of differential corrections to get the true periodic solution.
The Hill problem has been chosen to study analytically the dynamics about the libration points because of its generality and simplicity. However, the procedures used in this research are general and can be analogously applied to the restricted threebody problem or variations of it including different perturbations.
Acknowledgemnts
Partial support by the Spanish State Research Agency and the European Regional Development Fund under Project ESP201676585R (AEI/ERDF, EU) is recognized. The first author (ML) was also supported by project ESP201341634P of the same agencies.
References

[1]
V. Szebehely,
Theory of
Orbits. The Restricted Problem of Three Bodies, Academic Press Inc., New
York and London, 1967.
doi:10.1016/B9780123957320.500192.
URL http://www.sciencedirect.com/science/book/9780123957320  [2] G. Gómez, A. Jorba, J. Masdemont, C. Simó, Study of the transfer from the Earth to a halo orbit around the equilibrium point L1, Celestial Mechanics and Dynamical Astronomy 56 (1993) 541–562. doi:10.1007/BF00696185.

[3]
G. Gómez, A. Jorba, J. Masdemont, C. Simó,
Study
of the transfer between halo orbits, Acta Astronautica 43 (9) (1998) 493 –
520.
doi:http://dx.doi.org/10.1016/S00945765(98)001775.
URL http://www.sciencedirect.com/science/article/pii/S0094576598001775  [4] W. S. Koon, M. W. Lo, J. E. Marsden, S. D. Ross, The Genesis Trajectory and Heteroclinic Connections, in: K. C. Howell, F. R. Hoots, B. Kaufman, K. T. Alfriend (Eds.), AAS/AIAA Astrodynamics Conference 1999, Vol. 103 of Advances in the Astronautical Sciences, American Astronautical Society, Univelt, Inc., USA, 2001, pp. 2327–2343.

[5]
E. M. Alessi, G. Gómez, J. J. Masdemont,
Leaving
the moon by means of invariant manifolds of libration point orbits,
Communications in Nonlinear Science and Numerical Simulation 14 (12) (2009)
4153 – 4167.
doi:http://dx.doi.org/10.1016/j.cnsns.2008.09.016.
URL http://www.sciencedirect.com/science/article/pii/S1007570408003110 
[6]
G. Gómez, M. Lo, J. Masdemont,
Libration Point
Orbits and Applications, World Scientfic, 2003.
URL https://books.google.com.br/books?id=06zUCgAAQBAJ 
[7]
H. Lei, B. Xu,
Highorder
solutions around triangular libration points in the elliptic restricted
threebody problem and applications to low energy transfers, Communications
in Nonlinear Science and Numerical Simulation 19 (9) (2014) 3374 – 3398.
doi:https://doi.org/10.1016/j.cnsns.2014.01.019.
URL http://www.sciencedirect.com/science/article/pii/S1007570414000446 
[8]
F. Topputo,
Fast
numerical approximation of invariant manifolds in the circular restricted
threebody problem, Communications in Nonlinear Science and Numerical
Simulation 32 (2016) 89 – 98.
doi:https://doi.org/10.1016/j.cnsns.2015.08.004.
URL http://www.sciencedirect.com/science/article/pii/S1007570415002816  [9] M. Hénon, Exploration numérique du problème restreint. I. Masses égales ; orbites périodiques, Annales d’Astrophysique 28 (1965) 499.
 [10] M. Hénon, Exploration numérique du problème restreint. II. Masses égales, stabilité des orbites périodiques, Annales d’Astrophysique 28 (1965) 992.
 [11] J. V. Breakwell, J. V. Brown, The ’halo’ family of 3dimensional periodic orbits in the earthmoon restricted 3body problem, Celestial Mechanics 20 (1979) 389–404. doi:10.1007/BF01230405.
 [12] K. C. Howell, ThreeDimensional Periodic Halo Orbits, Celestial Mechanics 32 (1984) 53.
 [13] M. Lara, R. Russell, B. F. Villac, Classification of the Distant Stability Regions at Europa, Journal of Guidance Control Dynamics 30 (2007) 409–418. doi:10.2514/1.22372.
 [14] K. C. Howell, H. J. Pernicka, Numerical determination of Lissajous trajectories in the restricted threebody problem, Celestial Mechanics 41 (1988) 107–124.

[15]
G. Gómez, J. Mondelo,
The
dynamics around the collinear equilibrium points of the {RTBP}, Physica D:
Nonlinear Phenomena 157 (4) (2001) 283 – 321.
doi:https://doi.org/10.1016/S01672789(01)003128.
URL http://www.sciencedirect.com/science/article/pii/S0167278901003128  [16] R. W. Farquhar, A. A. Kamel, QuasiPeriodic Orbits about the Translunar Libration Point, Celestial Mechanics 7 (1973) 458–473. doi:10.1007/BF01227511.
 [17] D. L. Richardson, Analytic construction of periodic orbits about the collinear points, Celestial Mechanics 22 (1980) 241–253. doi:10.1007/BF01229511.
 [18] À. Jorba, J. Masdemont, Dynamics in the center manifold of the collinear points of the restricted three body problem, Physica D Nonlinear Phenomena 132 (1999) 189–213. doi:10.1016/S01672789(99)000421.

[19]
J. J. Masdemont,
Highorder expansions
of invariant manifolds of libration point orbits with applications to mission
design, Dynamical Systems 20 (1) (2005) 59–113.
arXiv:http://dx.doi.org/10.1080/14689360412331304291, doi:10.1080/14689360412331304291.
URL http://dx.doi.org/10.1080/14689360412331304291  [20] G. W. Hill, Researches in the Lunar Theory, American Journal of Mathematics 1 (1878) 5–26.
 [21] J.M. Petit, M. Hénon, Satellite encounters, Icarus 66 (1986) 536–555. doi:10.1016/00191035(86)900898.
 [22] M. A. Vashkov’yak, Evolution of the orbits of distant satellites of Uranus, Astronomy Letters 25 (1999) 476–481.

[23]
M. Lara, R. P. Russell, B. Villac,
Fast estimation of stable
regions in real models, Meccanica 42 (5) (2007) 511–515.
doi:10.1007/s110120079060z.
URL http://dx.doi.org/10.1007/s110120079060z 
[24]
M. Lara, J. Palacián, R. Russell,
Mission design through
averaging of perturbed Keplerian systems: the paradigm of an Enceladus
orbiter, Celestial Mechanics and Dynamical Astronomy 108 (1) (2010) 1–22.
doi:10.1007/s1056901092862.
URL http://dx.doi.org/10.1007/s1056901092862 
[25]
J. Peláez, M. Lara, C. Bombardelli, F. R. Lucas, M. SanjurjoRivo,
D. Curreli, E. C. Lorenzini, D. J. Scheeres,
Periodic Orbits of a
HillTether Problem Originated from Collinear Points, Journal of guidance,
control, and dynamics 35 (1) (2012) 222–233.
doi:10.2514/1.53097/1.53097.
URL http://arc.aiaa.org/doi/abs/10.2514/1.53097  [26] M. Hénon, Numerical Exploration of the Restricted Problem, V. Hill’s Case: Periodic Orbits and their Stability, Astronomy and Astrophysics 1 (1969) 223–238.
 [27] M. Hénon, Vertical Stability of Periodic Orbits in the Restricted Problem. II. Hill’s Case, Astronomy and Astrophysics 30 (1974) 317.
 [28] M. Hénon, New Families of Periodic Orbits in Hill’s Problem of Three Bodies, Celestial Mechanics and Dynamical Astronomy 85 (2003) 223–246.
 [29] M. Michalodimitrakis, Hill’s problem  Families of threedimensional periodic orbits. I, Astrophysics and Space Science 68 (1980) 253–268. doi:10.1007/BF00641660.
 [30] C. Zagouras, V. V. Markellos, Threedimensional periodic solutions around equilibrium points in Hill’s problem, Celestial Mechanics 35 (1985) 257–267. doi:10.1007/BF01227656.
 [31] C. Simó, T. J. Stuchi, Central stable/unstable manifolds and the destruction of KAM tori in the planar Hill problem, Physica D Nonlinear Phenomena 140 (2000) 1–32. doi:10.1016/S01672789(99)002110.

[32]
G. Gómez, M. Marcote, J. M. Mondelo,
The invariant manifold
structure of the spatial Hill’s problem, Dynamical Systems 20 (1) (2005)
115–147.
doi:10.1080/14689360412331313039.
URL http://dx.doi.org/10.1080/14689360412331313039 
[33]
M. Lara, A Hopf variables
view on the libration points dynamics, Celestial Mechanics and Dynamical
Astronomy (on line).
doi:10.1007/s1056901797784.
URL https://doi.org/10.1007/s1056901797784 
[34]
M. Kummer, On resonant
non linearly coupled oscillators with two equal frequencies, Communications
in Mathematical Physics 48 (1976) 53–79.
doi:10.1007/BF01609411.
URL http://projecteuclid.org/euclid.cmp/1103899811  [35] R. Cushman, Geometry of the Bifurcations of the Normalized Reduced HenonHeiles Family, Proceedings of the Royal Society of London Series A 382 (1982) 361–371. doi:10.1098/rspa.1982.0106.
 [36] H. Hopf, Über die Abbildungen der dreidimensionalen Sphäre auf die Kugelfläche, Mathematische Annalen 104 (1931) 637–665.
 [37] A. Deprit, The Lissajous transformation. I  Basics, Celestial Mechanics and Dynamical Astronomy 51 (1991) 201–225. doi:10.1007/BF00051691.
 [38] A. Celletti, G. Pucacco, D. Stella, Lissajous and Halo Orbits in the Restricted ThreeBody Problem, Journal of NonLinear Science 25 (2015) 343–370. doi:10.1007/s0033201592322.
 [39] M. Lara, A Hopf variables view on the libration points dynamics, ArXiv eprintsarXiv:1703.02887.
 [40] A. Giorgilli, L. Galgani, Formal integrals for an autonomous Hamiltonian system near an equilibrium point, Celestial Mechanics 17 (1978) 267–280. doi:10.1007/BF01232832.
 [41] A. Deprit, Canonical transformations depending on a small parameter, Celestial Mechanics 1 (1) (1969) 12–30. doi:10.1007/BF01230629.
 [42] D. Boccaletti, G. Pucacco, Theory of orbits. Volume 2: Perturbative and geometrical methods, 1st Edition, Astronomy and Astrophysics Library, SpringerVerlag, Berlin Heidelberg New York, 2002.
 [43] A. Deprit, A. Rom, The Main Problem of Artificial Satellite Theory for Small and Moderate Eccentricities, Celestial Mechanics 2 (2) (1970) 166–206.
 [44] M. Lara, J. SanJuan, Dynamic Behavior of an Orbiter Around Europa, Journal of Guidance, Control and Dynamics 28 (2) (2005) 291–297. doi:10.2514/1.5686.
 [45] A. Deprit, A. Delie, Trojan orbits. I. d’Alembert Series at L, Icarus 4 (1965) 242–266. doi:10.1016/00191035(65)900023.
 [46] G. Gómez, A. Jorba, J. Masdemont, C. Simó, Study Refinement of SemiAnalytical Halo Orbit Theory, Technical Report Contract 8625/89/D/MD(SC), European Space Operations Center, RobertBoschStrasse 5, 64293 Darmstadt, Germany (1991).
 [47] J. Henrard, Periodic Orbits Emanating from a Resonant Equilibrium, Celestial Mechanics 1 (1970) 437–466. doi:10.1007/BF01231143.
 [48] J. Henrard, A survey of Poisson series processors, Celestial Mechanics 45 (1988) 245–253. doi:10.1007/BF01229007.
 [49] J. F. SanJuan, ATESAT: Automatization of theories and ephemeris in the artificial satellite problem, Technical Report CT/TI/MS/MN/94250, Centre National d’Études Spatiales, 18, avenue Edouard Belin  31401 Toulouse Cedex 9, France (May 1994).

[50]
A. Jorba, A methodology
for the numerical computation of normal forms, centre manifolds and first
integrals of Hamiltonian systems, Experimental Mathematics 8 (2) (1999)
155–195.
URL http://projecteuclid.org/euclid.em/1047477059  [51] A. Celletti, L. Chierchia, Construction of analytic KAM surfaces and effective stability bounds, Communications in Mathematical Physics 118 (1988) 119–161. doi:10.1007/BF01218480.