A review of periodic orbits in the circular restricted three-body problem

ZHANG Renyong

Key Laboratory of Space Utilization,Technology and Engineering Center for Space Utilization,Chinese Academy of Sciences,Beijing 100094,China

Abstract: This review article aims to give a comprehensive review of periodic orbits in the circular restricted three-body problem(CRTBP),which is a standard ideal model for the Earth-Moon system and is closest to the practical mechanical model.It focuses the attention on periodic orbits in the Earth-Moon system.This work is primarily motivated by a series of missions and plans that take advantages of the three-body periodic orbits near the libration points or around two gravitational celestial bodies.Firstly,simple periodic orbits and their classification that is usually considered to be early work before 1970 are summarized,and periodic orbits around Lagrange points,either planar or three-dimensional,are intensively studied during past decades.Subsequently,stability index of a periodic orbit and bifurcation analysis are presented,which demonstrate a guideline to find more periodic orbits inspired by bifurcation signals.Then,the practical techniques for computing a wide range of periodic orbits and associated quasi-periodic orbits,as well as constructing database of periodic orbits by numerical searching techniques are also presented.For those unstable periodic orbits,the station keeping maneuvers are reviewed.Finally,the applications of periodic orbits are presented,including those in practical missions,under consideration,and still in conceptual design stage.This review article has the function of bridging between engineers and researchers,so as to make it more convenient and faster for engineers to understand the complex restricted three-body problem (RTBP).At the same time,it can also provide some technical thinking for general researchers.

Keywords: periodic orbit,restricted three-body problem,classification,orbit family,bifurcation.

1.Introduction

There is a growing need for space missions that utilize orbits in the cislunar space.In the past five years,several periodic orbits in the Earth-Moon system have been utilized,such as Artimes [1] inL1/L2Lagrange points,distant retrograde orbit (DRO) of asteroid redirect mission(ARM) [2,3],Lunar Orbital Platform-Gateway (LOP-G)[4–8] in near-rectilinear halo orbit (NRHO),transiting exoplanet survey satellite (TESS) [9–12] in a 2:1 resonant orbit in the cislunar space.Wind’s [13–18] original mission is to orbit the Sun at theL1-Lissajous orbit and also the first mission to utilize distant prograde orbit (DPO) concept.These missions brought new elements to aerospace engineering.

The limitations of many technical areas,including mission design capabilities,have hampered the thorough exploration of space.Therefore,there is a need for more innovative research in orbital dynamics design and its complex computational methods.The increased complexity of the dynamical model,including perturbation,provides numerous mission orbits for mission design.An in-depth study of these mission dynamics can contribute to scientific discovery and promote space exploration.

It has been historically proven that an effective method for acquiring new sciences in the solar system requires simultaneous consideration of the motion of the spacecraft associated with the Sun,Earth,and Moon.For example,the International Sun/Earth Explorer (ISEE-3) mission of National Aeronautics and Space Administration (NASA)provides valuable information on solar flares and Gammaray bursts [19].In addition,the Genesis spacecraft monitored solar wind particles within two years [20].The orbital shape used by the ISEE-3 and Genesis spacecraft facilitates the collection of scientific data and communication relays with the Earth.

The basic models of the above missions take into account the gravitational effect of the Sun,Earth,and Moon,which is based on the three-body problem.The threebody problem has always been the focus of mathematics and scientific research.On the other hand,the periodic orbits provide more insights to understand the complex dynamical system of the restricted circular three-body problem.The dynamics is chaotic and periodic orbits show order with different degrees of stability.Since the seminal work published by Poincare [21] in 1892,significant attention has been paid to finding periodic solutions.(Actually,Poincare believed that periodic orbits are the only way to understand complex three-body problems.)Early research quickly narrowed the study to periodic solutions for planar problems and identified many periodic orbit families.However,in 1920,Moulton [22] considered the three-dimensional problem,in which study he focused on the“oscillating satellite”near the collinear equilibrium point,and the ultimate goal was to calculate the periodic solution of the satellite.

From the 1930s to the 1960s,the three-body problem seemed to have been more or less forgotten.Since the 1960s,with the launch of the American Apollo program,research hotspots shifted from the initial celestial dynamics to the orbital dynamics of the Earth-Moon system.The study of restrictive three-body problems has become the focus of attention.

In 1966,Farquhar began to study the collinear translation point and explored the spacecraft oscillating at theL2point of the Earth’s Moon to provide services for the back of the Moon using communication satellites [23–25].In NASA’s research,the proposed solution involved a satellite that oscillates back and forth on the Earth and lunar motion planes.A solution of this type is unsatisfactory because the satellite periodically passes behind the Moon and cannot be seen from Earth.Therefore,an out-of-plane solution was sought to solve this problem.

Breakwell et al.[26] summarily calculated the energy consumption required fora continuous controller that forces equal frequencies,thereby ensuring an out-of-plane periodic solution that can always be seen from Earth.Therefore,the concept of“halo orbit”was advanced.In fact,this study of bounded motion near the collinear Lagrangian point opens up a new direction for spacecraft mission design.Breakwell et al.[27] gave higher-order periodic approximation solutions to the restricted threebody problem (RTBP),which makes it easy to calculate a halo orbital numerical solution.A naturally occurring solution was characterized by a specific relationship between in-plane and out-of-plane amplitudes.Unstable resident operations were required.Owing to the inherent instability of this type of orbit,orbit maintenance was required.

In the last century,most of the deep-space exploration missions were carried out in Sun-Earth system,mainly based on scientific observation.Since the U.S.President’s Vision for Space Exploration [28] announcement in January 2004,preparations have been under way for autonomous and manned activities on the lunar surface.For the development of the trajectory design options for a return to the Moon,incorporation of the gravity effects for not only the Moon,but also the Earth and Sun,offers very useful information for space craft motion in the vicinity of the Moon.The results allow flexibility for mission design and ultimately will facilitate exploration of Mars.

At present,the development of the near-Earth space mission has reached its end,and research on the Earth-Moon system is the most urgent task for humans to fly beyond near Earth space to deep space.The research on system dynamics of the Earth and Moon has also shifted from the initial theoretical research to engineering requirements.

Recently,NASA intended to build the Gateway as part of the since-cancelled Asteroid Redirect Mission [29,30].An informal joint statement on NASA-Roscosmos cooperation was announced on 27 September 2017 [31].Traveling to and from cislunar space (lunar orbit) would facilitate advancing the knowledge and experience necessary to journey beyond the Moon and into deep space.The project LOP-G would initially enter a NRHO around the Moon [32],which is a type of periodic orbit of the RTBP.

The periodic orbits are closely related to the type of above space missions,and the dynamics of the RTBP is the key to the study of the Earth-Moon system.It is necessary to have a systematic understanding and research of a periodic orbit.The implementation of the NASA LOP-G will push the dynamics of the Earth-Moon system studies to a climax.Meanwhile,the National Priority Project of the Chinese Academy of Sciences will support the Earth-Moon space exploration research in 2022.

This review is to present periodic orbits in the circular restricted three-body problem (CRTBP).A pair of celestial bodies include a primary and a minor body,and a spacecraft is considered a massless body.In the existing literature,the body pair include Earth-Moon,Sun-Earth,and Sunplanet,planet and its satellite.The examples of figures in this paper are all those in the Earth-Moon system.This work of the article has the function of bridging between engineers and researchers,so as to make it more convenient and faster for engineers to understand the complex RTBP.At the same time,it can also provide some technical thinking for general researchers.

2.CRTBP and periodic orbit

2.1 CRTBP

The CRTBP provides an autonomous approximation of the Earth’s lunar dynamics,enabling understanding of its underlying dynamic structure.The ephemeris model is a higher-fidelity simulation,and the CRTBP dynamic model is a reasonable approximation of the higher-fidelity dynamics model in the Earth-Moon system.In the application of the CRTBP,the movement of massless spacecraft under the influence of the gravitational pull of the Earth and the Moon is considered.Assume that the two main gravitational objects are mass points and make a circular motion around their common barycenter.The motion of the spacecraft is then described in a rotating coordinate system relative to the Earth and the Moon.By convention,the parameters in the CRTBP are dimensionless,and the dimensionless masses of the Earth and the Moon are equal to (1−µ) and µ,respectively,where the parameterµis equal to the ratio of the lunar mass to the total mass of the system.

In rotation coordinates,thexaxis is along the vector between the primaries,pointing to the small primary;thezaxis is parallel to the Keplerian primary orbital angular velocity vector;and theyaxis forms a Cartesian coordinate with thexandzaxes.Then,the motion equations of the third body are written as

wherer1andr2are the relative distances between the third body and two primaries,respectively.The equations are non-dimensional,and the characteristic quantities are the total mass,distance between the primaries,and magnitude of the system angular velocity.

There are five relative equilibrium points in the CRTBP system,including three collinear pointsL1,L2,andL3,and two equilateral pointsL4andL5.The collinear points located along the Earth-Moon line and the equilateral points form equilateral triangles with the two primaries.The CRTBP uses a Jacobi constant energy integral such thatJC=2U−v2,whereU(x,y,z)=(x2+y2)/2+(1−µ)/r1+µ/r2+µ(1−µ)/2 is the pseudo-potential function andvis the velocity magnitude relative to the rotating frame.

2.2 Mathematical modelling of periodic orbits in CRTBP

It is necessary to briefly introduce the concept of periodic orbits and quasi-periodic orbits in the CRTBP.In general,the periodic motion of a dynamic system repeats the same motion at equal intervals of time,including the motions that are repeated in a relative sense.Now have an Earth-Moon rotating coordinate system and let the system of differential equations [33] be written as

which possesses a particular solution,xi=φi(t),if the functions φihave the property that

In other words,when the initial value of φiis the same as its value at an epochTand consequently φi(t)=φi(t+T),then the functions φiare periodic intwith periodT,andxi=φi(t) is a periodic solution of the differential equations.

A continuous functionx(t) is almost periodic if,given ε>0,there exists anE(ε) such that for every real number α there is aTthat satisfies α≤T≤α+E(ε) and for which

Whenx(t) is a periodic function,ε=0 andTbecomes the period.An almost periodic function,therefore,is“periodic”with an“error”ε.A special class of almost periodic functions is called quasi-periodic,which can be seen in Farquhar (1973) [34] in the RTBP.In general,a quasi-periodic orbit is preferable to a periodic one,owing to the larger number of parameters that characterize quasiperiodic orbits [35].

2.3 Preliminary analysis by Poincare method

Poincare methods are commonly used to analyze periodic or quasi-periodic orbits.It is useful to employ Poincare methods to analyze more complete picture of the available libration point orbits at a particular energy level in the RTBP,and for potential design options.Through the use of a Poincare map,ann-dimensional continuous-time system is reduced to a discrete-time system of (n−1)-dimensions.By additionally constraining the Jacobi constantC,the problem is reduced to (n−2)-dimensions and,thus,the map for the CRTBP is represented in fourdimensional (for three-dimensional space).To generate a planar Poincare map,a surface-of-section,Σ,is defined such that Σ is transverse to the flow.A commonly usedΣ is one that represents a surface-of-section corresponding to crossings of thex-yplane.To compute the map,trajectories are integrated and crossings of Σ are recorded and displayed [36].

Consider the map in Fig.1 as projected into thex-yplane,produced to resemble the maps demonstrated by Folta [36].For the selected value ofC,several periodic orbits exist (see Fig.3),including a planar Lyapunov orbit (green,as shown in Fig.2),an axial orbit (see Fig.4)a vertical orbit (dark blue,as shown in Fig.6),and the northern and southern halo orbits.In thex-yprojection,the halo orbits share the same crossings of the map;the northern halo is featured in magenta in Fig.1.

Fig.1 Poincare map corresponding to crossings of the x-y plane,with sample orbits featured [36]

Fig.2 L1/L2 Lyapunov orbit families in the Earth-Moon rotating coordinate frame

Fig.3 Multiple families of periodic orbits in the Earth-Moon rotating coordinate frame

Fig.4 Axial periodic orbit family bifurcated form planer orbit in the Earth-Moon rotating coordinate frame

The first step of Poincare methods is to find a convenient Poincare sectionΣ.This can be either a section in time or space.It is assumed for visualization purposes that the section is taken in phase space.The main concern when choosing the plane of the section is ensuring that the velocity vector of the quasi-periodic orbit is as transverse to the plane of section as possible.This reduces the possibility that the integrated points will not return to the Poincare section.Thus,a good candidate for the Poincare section is the plane perpendicular to the velocity of the halo-orbit section.However,for the RTBP,a simple section on the ecliptic plane is also suitable.The Poincare section method can help guess the initial solutions of periodic or quasi-periodic orbits,and the precise solution can conveniently be computed by using the multiple-shooting method [37].

2.4 Methodologies and fundamental techniques

2.4.1 Differential corrections

Usually,when calculating a periodic orbit,one must first give an approximate first order linear or guessed approximation initial solution.Obviously,the evaluation of the approximate solution plays an important role in analyzing the properties of the periodic orbit around the equilibrium point within a certain accuracy requirement.However,the periodic orbit obtained by such an approximate solution is not the exact periodic solution of the actual restricted three-body model,so it cannot meet the accuracy requirements of the actual space mission orbit.Therefore,the analytical approximate solution method must be combined with the numerical method of differential correction to obtain an accurate initial value that satisfies the requirements of the periodic orbit.In this way,the problem of solving the orbital periodic solution is transformed into the problem of determining the initial value that satisfies the requirements of the periodic orbit;that is,the value of a certain moment determined by the linear approximation solution is used as the initial value of the periodic orbital guess.Through the method of differential correction,the initial value satisfying certain precision requirements is obtained through iteration as the initial value solution of the periodic orbit.

In general,differential correction can be implemented using the Newton method.The reference orbitx=φ(t,t0;x0),assuming that the expected state value at timet1isxd,deviates from the actual state valuex(t1) on the reference orbit at timet1as δx1,that is,

Suppose it is necessary to make a correction to δx0for the initial valuex0,so that the state value oft1isxd,that is,

Defining the functionF(x)=φ(t1,t0;x)−xd,related to the initial valuex0,the problem translates into using the Newton iteration method to find the appropriatexvalue that satisfiesF(x)=0.For a given reference valuex0,onlyF(x0)=φ(t1,t0;x0)−xd≠0 is needed to find δx0,which satisfiesF(x0+δx0)=φ(t1,t0;x0)−xd=0.Expanding with a Taylor series,

Ignoring the high-order small quantities,the differential corrections δx0=−DF−1(x0)F(x0) are obtained,that is,δx0=Φ(t1,t0)−1δx1.

2.4.2 Single and multiple shooting methods

Several techniques have been developed for the computation of periodic solutions.Single and multiple shooting methods are the most widely used techniques in existing software packages for the continuation of periodic solutions.A single-shooting technique solves the nonlinear system for the unknownx0andT.φ (x0,T,γ) is computed with an initial value solver.Since ordinary differential equation systems resulting from the discretization of parabolic partial differential equations are very stiff,implicit techniques should be used.The nonlinear system is usually solved by using the Newton’s method.This is an iterative method.At each iteration step,the dense linear system [38,39]

is solved,e.g.,using Gaussian elimination with pivoting,wheresis a phase condition and subscript ‘x’ and ‘T’denotes differentiation,and the variables are updated according to

This means,in single shooting,that the system over time τ represented by φ(x0) originated with initial statex0is uncovered by integrating the differential equations of the system with a standard explicit propagator.Generally,x0is adjusted with an update scheme such as the Newton’s method,to yield the desired final stateXf=φ(x0).Thus,for single shooting,all the problem sensitivities are associated on the initial statex0.This methodology for producing trajectories relies on an accurate initial vectorx0and is sometimes difficult for numerical convergence when the dynamics are sensitive to small adjustments in the initial state [40].

The single and multiple shooting methods could be implemented in solving an optimal control problem.Ultimately,single and multiple shooting achieve the same result.The multiple-shooting segments are introduced strictly for numerical reasons since the single shooting algorithm tends to fail when the time domain increases [41].However,a multiple shooting scheme offers better access to the design space,especially when the dynamics are sensitive,as often occurs in chaotic,nonlinear systems.Since sensitivities are distributed across the entire trajectory,a multiple shooting algorithm accommodates a larger convergence radius than an equivalent single shooting method [40].The extension to the multiple shooting method is obtained by dividing the entire time domain into sub-intervals and applying the simple-shooting procedure within each sub-interval.

2.4.3 Continuation method

A family of periodic orbits can be computed by continuation and differential corrections.The periodic orbit family of the RTBP is continuous.When the individual periodic orbits are computed,they can be used as the initial orbit to generate other new orbits by using the continuation method,by extending certain parameters.It is possible to quickly calculate all the periodic orbits,which are sought beyond the range of validity of the approximations near the initial orbit [42,43].Apply the continuation method to study the continuation of periodic orbits of the three-body problem [44–46].

3.Periodic orbits and their classification

3.1 Classifications of periodic orbit before the 1960s

The RTBP is one of the most famous problems in the modern history of astronomy and mathematics.The pioneers in the history of the problem include Newton,Lagrange,Jacobi,Poincare,Darwin,Moulton,and other famous mathematicians.The main contributions are as follows.

Most of the dynamical analysis of the three-body problem mechanics was developed by Newton (1643–1727) in 1687.Newton’s predecessor,Euler (1707–1783),was the first to simplify the problem by assuming the infinitesimal mass of the third body,thus introducing a“restricted”three-body problem.Euler is also the first person to model the dynamics in a rotating reference frame.Jacobi (1804–1851) first integrated the three-body problem model.Lagrange (1736–1813),at the same time,accompanied by Euler,identified five equilibrium solutions or liberation points.The Sun-Jupiter Trojan asteroid was correctly predicted by his solution.Almost 100 years later,Poincare (1854–1912) proved the non-integrability of the three-body problem.His work eventually became the basis of the modern dynamic system theory.The periodic orbit obtained by numerical analysis was first studied by Darwin (1845–1912) [47].

The aforementioned research was carried out during the 19th century.Contributions of the 20th century consist mainly of amplifications of Darwin’s work [48,49].In the early 1920s,various scholars began to study the method of solving the periodic orbit of the three-body problem and classifying the periodic orbits.Owing to the complexity of the problem,the research in this period mainly focused on the study of planar RTBPs.

Noteworthy,to correlate the Floquet theory and bifurcation analysis,see the details below,the concepts of families/classifications of periodic solutions,which are based on the definition by Goudas [50],are those that share a common hodograph.A hodograph is a continuous carve in phase space that consists of points belonging to different periodic solutions.In particular,the Floquet multipliers also change continuously.

3.1.1 Classifications of the Poincare problem

One of the many significant contributions of Poincare is the proof that there are periodic orbits in the three-body problem,and there are countless periodic solutions.Poincare divided the three-body problem’s periodic solutions into three categories.These three categories are the limiting problems of the two uncoupled Kepler problems [21],namely,the following:

(i) The first kind — two circular motions in the same plane.

In the first kind,the problem for µ=0 is extended by the method of analytic continuation to µ≠0.Research emphasis becomes periodic orbits for µ≠0.Direct and retrograde periodic orbits are demonstrated by this method.This problem is also used to analyze the characteristics of planets that have almost zero eccentricities and inclinations.In the limits at which the masses vanish,the orbits become circular,which is the singularity of the secular Hamiltonian.If the mean motionsn1andn2are the frequencies of the two circular Keplerian motions,the limit solution has periodT.The use of the method of analytic continuation to establish the existence of these orbits is synonymous with employing Poincare’s restricted problem,which treats“small”values of µ.Poincare’s problem is of primary interest in the dynamics of the solar system.

Hill’s inter mediate orbit used in his study of the Lunar problem is the most famous example of a period solution of the first kind.In his description of Hill’s result,Poincare corrected the incorrect guess by Hill regarding the global continuation of his solution.

(ii) The second kind — two elliptical motions in the same plane with resonant frequencies,and perihelia in conjunction or opposition and initial mean longitudes both equal to zero.

The second kind of problems are generated from elliptical two-body orbits in the plane of the primaries.The inclinations are still zero,but the eccentricities are finite(e≠0,i=0);in the limit,one obtains elliptical motions with the same direction of semi-major axes and conjunctions or oppositions at each half-period.

(iii) The third kind — two circular motions in different planes with resonant frequencies.

In the third kind,eccentricities are small,but inclinations are finite,and the limit motions are circular but inclined.Here,they face a bifurcation problem;indeed,they start from the degenerate and completely integrable situation in which families of periodic solutions of the reduced problem that exist for µ=0 are expected to break for µ≠0 and give rise to isolated periodic solutions.More precisely,they start from a resonant torus of dimension 2,corresponding to motions along two Keplerian ellipses.Poincare noticed that,provided one makes the ansatz that at some time the bodies are in symmetric conjunction,the dominant term of the expansion in eccentricities ofRis nothing but the reduced secular Hamiltonian [51].

3.1.2 Classifications of Moulton’s analysis

Moulton’s classification method originated from the“Periodic Orbit”published by Moulton et al.in 1920 [22,52].The main study was the analytical methods for periodic motion near the five equilibrium points.A power series method was proposed to make a third-order approximation of the solution near the triangular equilibrium point.His solution is similar to what is now called the vertical orbit [53].There are three types of finite periodic solutions near the collinear points that Moulton proved.The linearized motion near the collinear Lagrangian point and the differential equations relative to the equilibrium point are of the forms:

where ωxyis angular velocity of thex-yplane in rotating reference frame for the collinear lagrangian pointsL1,2,3and ωzis the angular velocity in thezdirection.Equation(8) will have periodic solution in three cases:

(i) ξ=η=0 and ζ is the form in (8c).

(ii) ξ=0 and η,ζ are the forms as (8b),(8c).

(iii) ξ,η,and ζ are the forms as (8a),(8b),and (8c),respectively,and ωxyand ωzare commensurable.

Moulton [22] proved that (8a),(8b) and (8c) can have finite periodic solutions in three-dimensional space,and gave corresponding numerical results.These numerical calculations depend on approximate analytical solutions near the equilibrium point.However,due to the limitations of the computing tools of the time,Moulton speculated that the solution of the third three-dimensional space could not be calculated numerically and proceeded with the analysis of planar orbits.In [33],14 different classes of periodic orbits in the three-body problem were identified;classesa,b,andcare defined as retrograde periodic orbits aroundL2,L1,andL3,respectively.

3.1.3 Classifications of Stromgren’s families

In the early 1930’s,the classical analyses by Sundman and Wintner concerning numerical experimentation in connection with the problem of three-body were highly praised,but the most outstanding example showing the importance of numerical undertakings is Stromgren’s [54,55] work performed between 1913 and 1939.The possible forms of motion of a third,indefinitely small mass,in the same plane as the two finite bodies,were to be investigated.A starting point for the research was found initially in certain known theorems of Lagrange that were adapted to the restricted problem.

Starting partly from the theorems of Lagrange and partly from other considerations,Thiekle,Burrau,and subsequently Darwin [48] (in the work“Periodic Orbits”)found a number of periodic solutions of the restricted problem.The main object of Stromgren’s research was to investigate the entire field of the restricted problem systematically for simple periodic solutions,and in practice this object falls into two parts:the discovery of new periodic solutions and the classification of old and new results in such a systematic way that the entire material is presented as a comprehensive survey of the possible forms of periodic motion.

The most complete work in the range0.1<µ≤1/2 was performed by Stromgren and his associates in Copenhagen.Darwin’s [48,49] and Moulton’s [22] contributions offer excellent commentaries and extensions of Stromgren’s work,but do not represent the completeness and precision achieved by the group at the Copenhagen Observatory.

The classification of Stromgren’s [55] families of orbits is based on the seven special points existing in the plane of the restricted problem:the five points of equilibrium and the two points of the primaries are located.Furthermore,while Stromgren discusses motions aroundL2,Szebehely [33] referred to motion aroundL3and compared the two systems further.In as much as in what follows the notation of the reference system,these aforementioned translation devices become important only when the reader wishes to study the original papers prepared by the Copenhagen Observatory.

Owing to the symmetry with respect to theyaxis,motions associated withm1andL3are the mirror images of the motions taking place aroundm2andL1,respectively.The expression“motion around”is widely and loosely used in the literature.When originally infinitesimal elliptical orbits around the collinear points become finite and,in fact,have large amplitudes,the word“around ”becomes meaningless.The proper use of“around”is guaranteed when it is restricted to the generating orbits.In this sense,they speak about seven special points and infinitesimal periodic generating orbits around them if these orbits exist.Because of the symmetry conditions in the Copenhagen problem,they haveL1,m2,L2(orL2,m1,andL3),andL4(orL5) as the only centers.The last point (L4) is not surrounded by infinitesimal periodic orbits for µ=1/2,the fact of which leaves three special points.There are no direct infinitesimal periodic orbits around the collinear points,only retrograde orbits,and there are direct and retrograde periodic orbits around the primaries.Consequently,the following families of periodic orbits were established by Stromgren,using his classification and the previously defined meaning of“around”[33].All of the families are shown in Table 1,and some families are shown in the following figures respectively.In Table 1,the acronym ‘RS’ represents reference system,‘FS’ represents fixed system.

Table 1 Classifications of Stromgren’s families

3.2 Classifications of periodic orbits after the 1960s

3.2.1 Classifications of Hill’s case

The general features of the solutions of Hill’s problem [56] can also be extended to study the plane-RTBP.The characteristic of Hill’s problem has the limiting case of the plane-RTBP when the mass of the second body tends toward zero.In a series of Henon’s work [57−60],the entirety of the solutions of the plane-RTBP were explored for a particular case that supposing that the masses of the first and second bodies are equal (µ=1/2).Henon [61] considered the range 0 ≤µ≤1/2 and that the case µ being small or zero covers most practical applications,µ of the Earth-Moon system being almost 0.012 of the RTBP.

It was shown that there were five families in Hill’s case,and Henon retained the nomenclature used in 1965a[57] and due to Stromgren [55]:familiesaandcoriginate in the Lagrangian orbits around the Lagrange equilibrium pointsL2andL1;familiesfandgbegin,respectively,as the retrograde and direct satellites of the second bodym2;and familyg′branches off familyg,at the critical orbitg1.Again,it cannot be completely excluded that other families of simple-periodic symmetrical orbits exist and have escaped their systematic search,but it is not very likely.

3.2.2 Classifications of three-dimensional periodic orbits

The three-dimensional problem gradually attracted scholars’ research interests after the 1960s.During the 20th century,researchers have been trying to find periodic solutions because they are the core of the study of understanding the non-integrable three-body problem.At the same time,many planar periodic orbital families are analyzed and calculated.With the dramatic increase in computing power,three-dimensional periodic families were emerging.Halo orbits are the closest to actual space mission design.Although it is proved that there are many three-dimensional periodic orbits,it is very difficult to calculate by grid search.

Owing to the development of computer technology,a larger range of three-dimensional periodic orbits can be calculated numerically.In the 1960s,Goudas [50] calculated 19 periodic orbital families based on the RTBP model,most of which had symmetry.Goudas extended Moulton’s study of relevant periodic orbits [22] but did not find the new familiar halo orbits.At the same time,the stability of periodic orbits is studied,and it is proved that the orbits of collinear Lagrange points are very unstable.Later,Bray [62] considered the large-scale threedimensional orbit problem again,because it could not continue to be linearized.They try to approximate analytical solutions as much as possible and use numerical computing to their advantage.

A comprehensive review of the restricted problem during this period was given by Szebehely [33] in 1967.For the related work that details more recent progress and potential future investigations into the problem,the reader is referred to Marchal [63].

(i) Axial and vertical orbits

Axial orbits are so named because the family appears to rotate about thex-axis from the bifurcating Lyapunov orbit.Vertical orbitals are called vertical because of their small amplitude motion in thezdirection [47].These families obtained numerical results near the collinear Lagrange points in 1966 [62].In this way,the study of three-dimensional orbits can be coplanar,and more attention is paid to out-of-plane perturbation and its influence on stability [64−67].The“vertical stability index”av was evaluated for planarorbit family.When the verticalstability index of the planar orbit is 1,it is a vertical critical orbit.It marks the bifurcation of planar and three-dimensional periodic orbits.Henon identified vertical critical orbits families a,b,and c which are aroundL3,L1,and ratio µ=0.00095 was explored by Kazantzis [68].L2,respectively,for mass ratio µ=1/2.The Sun-Jupiter The three-dimensional symmetric periodic orbits were classified by Goudas (1963) [50] as follows:

i) ‘Type A’ symmetric with respect to thex-zplane;

ii) ‘Type B’ symmetric with respect to thexaxis;

iii) ‘Type C’ both previous symmetries.

It is notable that the vertical critical orbits themselves can be far from the equilibrium point.The out-of-plane stability of the planar orbits was evaluated by Henon and did not involve the three-dimensional families.Vertical critical orbits can exist far from the equilibrium point.Robin [69] then added the theory of planar tridimensional orbits and analyzed their symmetry.The three-dimensional and planar projections of the Earth-MoonL2axial and vertical orbit families are shown in Fig.5 and Fig.6,respectively.

Fig.5 Three-dimensional and planar projections of the axial orbit families in the Earth-Moon rotating coordinate frame

Fig.6 Three-dimensional and planar projections of the L2 vertical orbit families in the Earth-Moon rotating coordinate frame

(ii) Halo orbits

The name“halo”for these orbits were first used by Farquhar in 1968 [23].Some researchers have advocated deploying relay satellites from the Apollo missions in Earth-MoonL2halo orbit.A halo orbit is an ideal location for a lunar communications and control center [70].Examining the vertical stability index (|av|=1) for the familiesa,bandc,it is shown that some critical orbits of the families havex-zplane symmetry property.Generated from plane periodic orbits,the three-dimensional familiesAlv,Clv,and Blv from Zagouras [71] are,in fact,theL2,L1,andL3halo families,respectively,in the Sun-Jupiter system.The authors do detect some stable orbits in theL1family,but the families are only partially computed and the stable orbits,now known to exist in the other two families,were not observed.The three-dimensional and planar projections of the Earth-MoonL1andL2halo orbit families are shown in Fig.7 and Fig.8,respectively.

Fig.7 Three-dimensional and planar projections of the L1 halo orbit families in the Earth-Moon rotating coordinate frame

Fig.8 Three-dimensional and planar projections of the L2 halo orbit families in the Earth-Moon rotating coordinate frame

(iii) DROs

DROs are a family of periodic orbital solutions of the motion equation in CRTBP,and its motion direction is opposite to that of the second gravitational body,which is derived from the numerical exploration of the Hill limit case of CRTBP [61].In this paper,the results of extending Henon vertical critical orbits to a three-dimensional finite family were given in [72,73].

Bifurcation families of orbits exist,and it is proved that three-dimensional orbits can actually be connected to multiple planar families of orbits.The infinitesimal periodic oscillations around the collinear Lagrange points are continued along the familiesa,b,andcof plane DROs aroundL1,L2,andL3[71].DROs can be extended to the study of the stability of plane periodic satellite orbits with respect to perturbations perpendicular to the plane for0 ≤µ≤1,which allows a determination of the three-dimensional stability of the orbits [74].The distant retrograde orbit family of the Earth-Moon system is shown in Fig.9.

Fig.9 Distant retrograde orbit families in the Earth-Moon rotating coordinate frame

(iv) NRHO

The near-rectilinear orbits were first proposed by Breakwell [27] in 1979,towards which the halo family tends asxmaxdecreases are amenable to an analytical approximation used in the Sun-Earth system [75].In the Earth-Moon system,the position is measured from the Moon and starts by supposing thatx,y≪z≪1.In 1984,using a three-dimensional method of regularization originally developed by Howell [76],it was applied to the restricted problem and used to develop approximations for“almost rectilinear”halo orbits near the collinear points.A thorough numerical investigation of halo orbits was also completed [77].More study on the“almost rectilinear”halo orbits in the later years and these orbits were named NRHOs [47,78].The three-dimensional and planar projections of the Earth-MoonL1andL2NRHO families are shown in Fig.10 and Fig.11,respectively.

Fig.10 Three-dimensional and planar projections of the L1 NRHO families in the Earth-Moon rotating coordinate frame

Fig.11 Three-dimensional and planar projections of the L2 NRHO families in the Earth-Moon rotating coordinate frame

(v) Butterfly orbits

The family of butterfly orbits are bifurcated from a sixdayL2NRHO and exhibit characteristic similar to NRHO.These orbits present a butterfly shape,and the lunar south pole remains in view for significant intervals of time [78].Robin and Markellos [69] studied butterfly orbits near the Moon in the Earth-Moon system.The butterfly’s orbital movement looks like a figure eight;such orbits allow observation of the Moon’s south pole in most of the time,due to the spacecraft’s slow motion at the far Moon point.The three-dimensional and planar projections of the Earth-Moon butterfly orbit families are shown in Fig.12.

Fig.12 Three-dimensional and planar projections of the butterfly orbit families in the Earth-Moon rotating coordinate frame

(vi)N-periodic orbits

Henon named a family of periodic orbits H5 [60,61],in which the study of families of periodic orbits in Hill’s problem was initiated.Because of computer limitations at the time,essentially only simple-periodic symmetric orbits were considered in H5,that is,orbits that intersect the horizontal ξ axis in two points only and which are symmetric with respect to that axis.

Then theN-periodic orbit was defined by Henon as such an orbit if it intersected the ξ axis 2Ntimes.Research definition of the study was extended to doubleperiodic and triple-periodic symmetric orbits.The threedimensional and planar projections of the Earth-Moon of the 5-period family orbit families are shown in Fig.13.

Fig.13 Three-dimensional and planar projections of the 5-period family orbit families in the Earth-Moon rotating coordinate frame

(vii) Resonant orbits

The motion of objects in the solar system is extremely complex but with high regularity.According to Newtonian motion mechanics,the gravitational interaction between physical forces causes this law of motion.One of the interesting phenomenon is resonance orbitals.In the case of periodic motion,resonance occurs when there is a simple numerical relationship between orbital periods [79].Resonance orbits can have different phenomena such as mean motion,Laplace,secular and Kozai resonance.The phenomenon of orbital-orbital resonance,where the period in question represents the orbits of two or more objects,will be described below.

Consider two objects A and B of arbitrary mass and describe their possible resonant motion relations.In the two-body Kepler motion,the motion period of object B ispand the motion period of object A isq,so the orbitalorbital resonance is defined by the parameterp:q.In the case of the Earth-Moon system,object A represents the Moon and object B represents the spacecraft.When the spacecraft just completedporbits around the Earth,and the Moon needs to completeqorbits around the Earth,the spacecraft and the Moon formed orbital resonance.For example,a spacecraft with a 1:3 resonance with the Moon completed one orbit around the Earth,while the Moon completed three cycles around the Earth.

In a multi-body dynamic system,the ratio ofp:qresonance to the period of the resonant object is not exactly equal.When the dynamics model is affected by more than one gravitational force,the time of one full flight of the affected object is not always constant.However,for thep:qresonance in the CRTBP,for example,in the Earth-Moon CRTBP model,the time of spacecraft completing the Earth’spcircle is infinitely close to that of the Moon completing theqcircle.The ratio of orbital periods is not rational,but an approximate rational fraction.The orbit is still closed when viewed in a rotating coordinate system [79].

The planar of the 1:1,1:2,1:3,1:4,2:1,2:3 and 3:2 resonant orbit families in the Earth-Moon system are shown in Fig.14 and Fig.15,respectively.

Fig.14 Planar of the resonant orbit families in the Earth-Moon rotating coordinate frame (ratio of 1∶ n)

Fig.15 Planar of the resonant orbit families in the Earth-Moon rotating coordinate frame (ratio of n∶ n)

4.Stability and bifurcation of periodic orbits

4.1 Stability of periodic orbit

The study of the characteristics of periodic orbits focuses on their stability.A periodic orbitX(x0,t) is said to be stable when a small change,ζ,in the starting conditionsx0does not lead to an appreciably different orbit.If,however,the particle describing a periodic orbit drift away when slightly perturbed,then this orbit is said to be unstable.

The linear and nonlinear problems revolve around the behavior of the dynamical system at and in the neighborhood of the equilibrium solutions.The process of linearization presents no difficulty at the Lagrangian points,nor does the analysis of the linear system’s stability,in principle,lead to any complications.The well-known theory of the orbital stability of these linear differential equations is reviewed below to improve the precision of their application to nonlinear systems.The interpretation of Lyapunov’s ideas and the work by LaSalle and Lefschetz,Cesari,and Coddington and Levinson is followed in part and the anticipated nonlinear problems discussed by Szebehely [33].

Two standard methods to evaluate the stability of periodic orbits were discussed by Szebehely.The first,Poincare’s approach,is the evaluation of the characteristic exponents by the integration of the variational equations.The problem of characteristic exponents was discussed by Whittaker [80] and followed Poincare’s [21] (also see Ince [81] and Cesari [82]).The second approach is the use of Darwin’s [48] equation of normal displacement,the derivation of which was given by Birkhoff [83].Other forms were given by Message [84].A discussion of the relation between the second-order equation for the normal displacement and the fourth-order variational equations was given by Wintner [85].The variational equations using the Thiele-Burrau regularizing variables were given by Rosenthal [86].The physical picture is,of course,weakened by the fact that the method of characteristic exponents and of normal displacements is inherently a linearized process,which often furnishes only necessary conditions for stability.The mathematical problem is the solution of Hill’s differential equations.In actual problems,the nature of the characteristic roots is of interest only,and Floquet’s method is available [87].No attempt is made here to review the problem of stability in celestial mechanics or the solutions of the restricted problem.

Goudas’s [50] method proposed in 1963 and Henon’s work in 1965 are the significant studies of periodic orbit stability in the RTBP.The significance of characteristic exponents goes beyond the stability problem,and it is not inconceivable to expect that the organization of families of periodic orbits is greatly enhanced by values of the characteristic exponents as shown by Henon [57–60],who also settled the linear stability question of the Copenhagen problem.

Henon evaluated the linear stability characteristics of the periodic orbit families of the Copenhagen category.Henon defined the horizontal stability parameters’a,b,c,dto analyze the stability of the periodic orbits,noting that Darwin’s [49] stability indexcis connected with Henon’s parameter by the equationa=cos(πc).Some of Henon’s major conclusions are the following:

(i) The classes of the Copenhagen problem may be organized in two groups according to linear stability behavior.In the first group are class (f),class (g),class (I),and class (m),and so this group contains the direct and retrograde satellite orbits and superior planetary orbits.In the second group are all the other classes.Only the first group contains stable orbits in any appreciable number;the members of the second group are almost all unstable.

(ii) Preparing plots of the relationa=a(C),whereCis Jacobian constant,only the first group seems to show values of ∂C/∂aappreciably different from zero in the range of stability,−1

(iii) Collision orbits show no special behavior,and they are found on the unstable parts of the curvesa(C).

(iv) Henon determined the functiona(C) numerically for class (I) of the Copenhagen category.These synodically retrograde orbits are direct orbits in the fixed system.Henon also analyzed the stability of all other families of periodic.

4.2 State transition matrix and stability index

Let φ(x0,t) represent the solution that results from numerically integrating the system differential equation (1) starting with the initial statex0and using the time variatet.The state transition matrix Φ describes how the solution flow φ(x0,t) changes with respect to the initial statex0,and is defined [40] as

The state transition matrix can be computed by numerically integrating the matrix differential equationsincex0andtare independent.

It is worth noting that at the initial time,φ(x0,t0)=x0,and therefore the initial conditions of the state-transition matrix are

whereIis the identity matrix.In addition,the stability index was defined by Breakwell and Brown [27].The symmetric matrixUXXof the second partial derivatives ofUwith respect tox,y,andz,were evaluated along the orbit.The state transition matrix of Φ at a re-crossing of thex-zplane is used to assist the convergence toward an orbit with a second perpendicular crossing of thex-zplane,i.e.,a periodic orbit.Stability is determined by the eigenvalues of the full cycle transition matrix Φ(tF,0),also called the monodromy matrixM,which is the state transition matrix after one full orbital periodT,orM=Φ(T,0).Two of the eigenvalues are always 1.The other four are in reciprocal pairs (λi,1/λi).The 6×6 matrix Φ(tF,0) can be reduced to a 4×4 matrixMto produce the four critical eigenvalues [27].Two stability indices have been defined as the arithmetic mean of each pair,vi=1/2(λi,1/λi).The stability indices can be calculated directly [88] from

where tr(·) is the trace of the matrix.Stability requires real values ofvbetween −1 and +1.It means that,if the stability index is less than or equal to one,the orbit is stable,and if the stability index is greater than one,it indicates that the orbit is unstable.The larger the stability index is,the faster the dynamics deviates from the orbit,that is,the faster the orbit diverges.In general,the stability index is directly related to the orbit maintenance cost and inversely proportional to the transfer cost.In general,the stability index increases as the orbital distance from the Moon increases.

In this paper,the stability of the main periodic orbit families of the Earth-Moon system is analyzed.Fig.16 and Fig.17 show the stability index analysis diagram of each periodic orbit family.It can be concluded from Fig.16 that axial family,vertical family,butterfly family and halo family are all unstable periodic orbit families.NRHO family,DRO family,DPO family all have parts that are stable.It can be concluded from Fig.17 that the families of resonance ratiop/q=1:1,1:3,2:3 are unstable,and the family of resonance ratiop/q=1:2 is partially stable,while most of them are unstable.Families of resonance ratiop/q=2:1 and 3:2 are all stable.These results of the stability analysis can provide a good reference value for choosing the Earth-Moon space orbit of the mission.

Fig.16 Stability index of the axial, vertical, halo, DRO, DPO,NRHO and butterfly periodic orbit families in the Earth-Moon system

Fig.17 Stability index of the resonant orbit families in the Earth-Moon system

4.3 Bifurcation of periodic orbit

Many families of periodic orbits,especially three-dimensional ones,exist,but they are extremely difficult to locate and compute,and a random numerical search will never be successful.Thus,study of the bifurcations,in which several families come together,which is critical and used as the basis of the study,is required.

A bifurcation is a break in the structure,or behavior,along a set of solutions.The bifurcation points that coincide with a change in the order of instability of the original solution encompass the traditional definitions of bifurcations.The bifurcations can be classified as traditional bifurcations (see Fig.18),which occur with order of instability change in the original family and other type bifurcations.The other types,includingn-period bifurcation,quasi-periodic bifurcation and modified secondary Hopf bifurcation,occur without a change in the order of instability along the original.

Fig.18 Bifurcation occuring when two pairs of eigenvalues move relatively to the unit circle in the complex plane in the Earth-Moon rotating coordinate frame

A change in stability properties indicates the presence of a bifurcation and represents an intersection between one or more different families of solutions.The families of solutions considered in this analysis are periodic orbits and/or invariant tori.Each type of bifurcation corresponds to a distinct type of stability change,i.e.,the change in the eigen-structure is different for different types of bifurcations.For a given mass parameterμ,there are 15 different types of traditional bifurcations possible in the RTBP [40].The traditional bifurcations are classified as either common or rare.

The three common bifurcation types are fold (a kind of tangent),period-doubling,and torus bifurcation [40].Rare bifurcations are the trifurcations,or the intersections of three different families of periodic orbits,which are characterized as a combination of fold.

The characteristics of the three types of traditional bifurcations are as follows:

(i) Tangent bifurcations,which occurs when a pair of eigenvalues move to (or from) the unit circle from (or to)the positive real axis.There are three types of this bifurcations defined for periodic solutions:cyclic-fold,pitchfork,and transcortical.The cycle-fold bifurcation represents only a change in stability of the family and no new periodic solutions exist.The only qualitative change is the order of instability.The pitchfork bifurcation signals a change in stability along the original family,as well as the appearance of a new family of periodic solutions.Of course,this implies that a hodograph corresponding to a new family of periodic orbits necessarily crosses the original hodograph.The stability of the new family at the bifurcation point must be equivalent to the stability of the original family across the intersection of the hodographs.However,while there is a change in stability along the original family,there is no change along the new family as it passes through the intersection point.

(ii) Period-doubling bifurcations,which occurs when a pair of eigenvalues move to (or from) the unit circle from(or to) the negative real axis.The stability characteristics of the original family change as the eigenvalues split at−1.Periodic orbits that are members of the new family have periods that are twice the length of the period that corresponds to orbits from the old family,but the orbits along the new family do not undergo a stability change at the bifurcation point.The three-dimensional and planar projections of the Earth-MoonL1period-doubling bifurcations orbit is shown in Fig.19.

Fig.19 Three-dimensional and planar projections of the L1 period-doubling bifurcations families in the Earth-Moon rotating coordinate frame

(iii) Torus bifurcations,also called Hopf bifurcation,which occurs when two pairs of eigenvalues collide on the unit circle,but off the real axis,and thus branch off the unit circle into the complex plane.Consequently,the order of instability of the family changes by two at this type of bifurcation point.In general,these collisions yield a new family of invariant tori about the single original periodic solution [89],indicating the presence of new,higher-dimensioned solutions.Both fold and period-doubling bifurcations initially appear to be special cases of a torus bifurcation.For a more thorough discussion of bifurcations,see [90].

The bifurcation point is relative to the vertical stability indexv.Hennon [57] has satisfactorily studied the intersection and existence of bifurcation of families of periodic orbits,as well as the conditions for the appearance of various types of branching points.In the three-dimensional version of this problem,only a few numerical examples of such bifurcation points were given by [71,91]and [77].Markellos [92] focused on detecting threedimensional asymmetric periodic orbits and determined the typical three-dimensional asymmetric periodic orbits of the bifurcating families.The bifurcations of families of plane-asymmetric periodic orbits with families of spaceasymmetric periodic (three-dimensional) orbits occur when the stability index |v|=1.Whenever the stability properties change across a family,a bifurcation orbit is identified.The bifurcation signals the intersection of two or more families of periodic orbits [93].

Intersections of families of three-dimensional periodic orbits that define bifurcation points were studied by Papadakis and Zagouras [93].The conditions of existence for bifurcation points were discussed and an algorithm for the numerical continuation of such points was developed for different µ values.Extensive investigation of periodic motion and bifurcations in the Sun-Earth/-Moon system for the collinear Lagrangian points in the Earth-Moon vicinity were completed by Howell and Campbell in 1999[94].A number of bifurcations and intersections representing the existence of other three-dimensional families were identified and various orbits were numerically computed as members of these intersecting families.

Typically,particular solutions were discovered by numerical grid searches through phase space,and a third dimension would clearly make the search unmanageably large.Once a significant number of periodic orbits were available,many authors described not only their periodic solutions,but also methods of determining the stability as well as a methodology for revealing an entire family of solutions generated from a single starting point.Additionally,families of such planar periodic solutions,as well as their stability,were often mapped out for the entire spectrum of possible primary mass ratios.In the early 1960s,in contrast to the two-dimensional trend,Goudas and Kazantzis numerically determined tens of seemingly unrelated doubly-symmetric three-dimensional periodic orbits and analyzed their stability [95].It was later discovered that the planar problem could be used to predict at which points three-dimensional periodic solutions would bifurcate out of the plane.Isolating these so called“vertical self-resonant orbits”became the method of choice for ultimately determining the more difficult threedimensional orbits [66,69,96].Similar to the planar case,families of these three-dimensional orbits were computed and mapped for different mass ratios of the primaries.In the 1970s,Markellos applied the mechanism for predicting and calculatingn-period bifurcating solutions from a family of particular one-period solutions for the planar RTBP [97].He used this information to compute many new two-and three-dimensional orbits.In the 1980s and most of the 1990s,the search for new particular solutions in the CR3BP was limited.

4.4 Out-of-plane bifurcation: from planar to three dimensions

However,there remain many unanswered questions.First,the full power of bifurcation theory has yet to be applied to the three-body problem.For instance,whilen-period bifurcations have been used to predict new particular solutions in the planar problem,they have not yet been applied to the three-dimensional problem (n-period bifurcations are scarcely mentioned in traditional texts on bifurcation theory) [98–101].Since many of the bifurcating three-dimensional families never degenerate into planar trajectories,they may not,in fact,have been analyzed to any extent,except in a few isolated cases.To date,all three-dimensional families without at least one known planar member can only be determined by a grid search in a three-dimensional phase space or by a systematic evaluation of three-dimensionaln-period bifurcating solutions.In this effort,the bifurcation method is used to determine a large number of solutions [89,94].

Papadakis et al.[93] studied bifurcations of families near the triangular Lagrangian points.In 1999,Howell et al.[94] completed an extensive investigation of periodic solutions and bifurcations in the Sun-Earth/-Moon system for the collinear Lagrangian points in the Earth/Moon vicinity.Dichmann et al.[102] obtained numerical results for connections between the Lyapunov,halo,axial,and vertical orbits in the vicinity of the collinear points in the Earth-Moon system using AUTO (a software for continuation and bifurcation problems).In 2003,Doedel et al.[103] further applied AUTO to successfully map periodic solutions near the collinear points to those in the vicinity of the triangular points [104].Their study also included numerical results for axial and vertical orbits in the vicinity of the triangular points.

Grebow [47] gave strategies for predicting initial guesses for solutions,based on first-order linear approximations,and locating bifurcations along a manifold were offered.Many different types of orbit families were generated,including planar,axial,and vertical families of orbits near all five Lagrangian points.A general assessment of the stability of these solutions is also summarized.The exact location of the bifurcation orbits can be computed using continuation and a method of bisections [91,105].

4.5 Multiplicity bifurcation: from simple to multi-period

Unlike the common bifurcation with order of instability change in original solution,which is discussed in Subsection 4.3,other types of bifurcations do not change the order of instability of the original family.This kind of the bifurcation is classified as rare which has three types.The first type of bifurcation occurs when the eigenvalue on the unit circle passes through annth unit root,λn=1(n>2)and is an integer [100,101].The three-dimensional and planar projections of then-periodic are plotted in Fig.20.These bifurcation values help to locate new,differentquality,n-periodic solutions [89].The second type is the generalization of then-period branch,which is represented here as a quasi-periodic branch.They produce new quasiperiodic solutions that appear on the annulus around the original periodic solution.This quasi-period type of bifurcation seems to appear more frequently than the second Hopf bifurcation;it is significantly different from the secondary Hopf because it does not trigger stability changes,which means eigenvalue collision does not occur on the unit circle.The three-dimensional and planar projections of the quasi-periodic are plotted in Fig.21.Finally,the third type of bifurcation is generated by collisions between two pairs of eigenvalues (excluding +1 and −1)on the real axis,which split into complex planes.In the last case,the order of instability does not change,but the transfer of eigenvalues from real to complex constitutes a qualitative change in the nature of the solution.The periodic orbit of the third type bifurcation did not plot here.The three-dimensional and planar projections of the Earth-MoonL2n-periodic orbit and quasi-periodic orbit are shown in Fig.20 and Fig.21,respectively.

Fig.20 Three-dimensional and planar projections of n-periodic orbit in the Earth-Moon rotating coordinate frame

Fig.21 Three-dimensional and planar projections of quasi-periodic orbit in the Earth-Moon rotating coordinate frame

5.Practical computation techniques of periodic orbits

In the 20 century,periodic orbit solutions have attracted interest because they are a key component in understanding behavior in the non-integrable three-body problem.It is also the most critical technical problem that human beings must solve in deep-space exploration.Several straightforward methods for computing periodic orbits are presented here.

In general,investigations concerning families of periodic orbits approach the problem in various ways [94]:generation of families of periodic solutions from a single arbitrary orbit (both two-and three-dimensional);computation of three-dimensional orbits from planar solutions;definition of a“family”based on mass ratio as the parameter of interest;and the construction of links between families by investigation of the relationships between families of planar orbits.All of the above approaches for periodic solutions are based on the following methodologies:analytical solution using the two-body approximation;analytical solutions near equilibrium points;bifurcation approach of computation;sweeping of computation areas;and Poincare method.

5.1 Organization of existing numerical results

Numerical methods for the computation of many different types of periodic solutions have evolved in the last 60 years owing to the increasing speed and accuracy of modern computers.Moreover,once a single periodic solution is determined,neighboring solutions in the same family are computed using a continuation method over an additional constraint parameter;this method is extrapolated to determine complete families of periodic orbits [78].

However,these calculations relied heavily on analytical approximations based on assumptions concerning the structure of generating solutions near equilibrium points.Thus,it was practically impossible to determine arbitrary three-dimensional solutions.Nevertheless,there has been an enormous amount of information published on periodic solutions since the early studies [94].

Definitions of the members of certain families are often questionable and quite arbitrary,designed only to satisfy the discoverer’s theorem about his orbits.In the following,the existing incomplete material of numerical results of periodic orbits in a systematic fashion was firstly attempted to organize.Then detail the existing methodologies for periodic solutions.Naturally,certain classifications or groups of orbits and/or computation approaches may have been missed because of limited capacity.

(i) Numerical studies of the case µ=0 are surprisingly rare;after the early work of Hill,Kelvin,and Jackson,a systematic search for periodic orbits was made only by Matukuma [106−108] and co-workers with the help of desk computers.

The first complete set of numerical results was published by members of Copenhagen Observatory under the direction of Stromgren during the period 1913–1939.The mass parameter for the work has the value of µ=1/2,indicating that the primaries are equal.Prior to this undertaking,Darwin,from 1897 to 1910 using µ=1/11,computed several systematic but incomplete families of periodic orbits.At around the same time (1900–1917),Moulton’s school performed analytical and numerical work,using in the latter mostly µ=0.5 and 0.2.

The above three undertakings all use values of the mass parameter between approximately 0.1 and 0.5,and the majority of these computations took place prior to the 1930s.All three aimed at establishing sets of periodic orbits,mostly by numerical,but also by analytical.For these reasons,the Stromgren-Darwin-Moulton results will form the first group of investigations to be presented.

(ii) The second category is characterized by µ≅0.012.The value of the mass parameter is approximate to the Earth-Moon system.The first systematic investigation of periodic orbits for the Earth-Moon system can be found in the work of Egorov [109].A more complete and,in certain aspects,more extensive work appeared later [110]by Broucke.Rather special types of periodic orbits were studied by Newton,Huang [111],and Arenstorf [112].

Therefore,the contributions in the second category are all concerned with the Earth-Moon model and correspondingly with a value of approximately 0.012.They establish families of periodic orbits and appear in the literature during the period 1957–1963.The orbits of this group are best represented by the description“periodic lunar orbits”.

(iii) Periodic orbits around the triangular Lagrangian points comprise the third category.Such orbits are discussed in considerable detail as continuations of the infinitesimal elliptic orbits in Szebehely [33] (Chapter 5).The two values of the mass parameter,for which such families are treated here,are µ≅0.012 14 (Earth-Moon system) and µ≅ 0.000 953 8 (Sun-Jupiter system).

(iv) The fourth category contains non-periodic orbits computed for the previously mentioned value ofµ≅0.012.The practical implications of these orbits are significant in space activities.Those investigations resulting in families of orbits will be our primary focus.Our scope is painfully limited by this restriction,but it is nevertheless believed that no other approach is feasible.The previously mentioned work by Egorov [109] is the first reference,followed by Thuring [113],Buchheim [114],and Szebehely [115] et al.,among others.These orbits will be referred to as lunar trajectories since either the entire family or at least some of its members can be used for Earth-to-Moon missions.

(v) The fifth category contains non-periodic orbits computed for purposes of stellar dynamics.The values of the mass parameter used are in the same range as for the Copenhagen group,namely 0.1 ≤µ≤0.5;however,the orbits are not periodic and are highly specialized.Proponents of these undertakings are Kuiper [116],Kopal[117],Abhyankar [118],and Gould [119],among others.This category is designated“applications to binary systems”.

(vi) The next category contains three special topics:asymmetric and symmetric periodic orbits in the sun-Jupiter system computed by Message,asymptotic periodic orbits at the collinear Lagrangian points established by Deprit and Henrard [120],and a three-fold analytic continuation of a second-kind periodic orbit.These are presented under the title“additional periodic orbits”.

(vii) The final category is the modern three-dimensional method.Farquhar first used the name“halo”and advocated using spacecraft in these orbits.Then,“almost rectilinear”halo orbits,also called NRHOs,near the collinear by Howell [76],were investigated.Retrograde periodic orbits were studied by Zagouras [71] and Benest [121] for 0 ≤µ≤1.Meanwhile,the relevant database generation is based on a grid search that incorporates a new discretization scheme centered around fundamental periodic orbit families,and a robust differential corrector implemented with a full second-order trust region method,which was investigated by Restrepo [122,123].The previous contributions are summarized in Table 2.

Table 2 Previous contributions of the circular restricted three-body problem and periodic orbits

Continued

5.2 Two-body approximation (two/three dimensional)

The distinguishing characteristic of the problem with a small mass ratio is that it is in the neighborhood of 1 with a known analytical solution:the case µ=0,or the twobody problem referred to rotating axes.One of the approaches in the search for periodic orbits in Broucke’s [88]study was to use the initial conditions for the two-body problem to find good approximations of the initial conditions for the periodic orbits in the three-body problem.Theoretically,it should be possible to start from the known solution for µ=0 and make an analytical continuation study for positive values.This method has been effectively used by many authors,notably by Poincare[21] among the earlier workers,and by Wintner [124],Barrar [125],and Arenstorf [112] in more theoretical investigations.Then,the existence of periodic solutionsx(t),deriving Keplerian motion of equations of motion for the plane-RTBP for small µ>0 were shown,which are near the generating solutionsx(t) belonging to arbitrary integers.

The two-body approach was quite fruitful because the small mass ratio µ of the problem is close to zero.If the mass ratio is zero,the solution is known because this is simply the two-body problem represented in a system of coordinates that is rotating with a constant angular velocity.To obtain initial conditions,the initial conditions for the two-body problem were used and then a rotation was applied [88].

5.3 Analytic solution around Lagrange points

In 1920,Moulton [22,52] was solely devoted to analytical methods for approximating periodic motion near all five equilibrium points,near which is the available linearization.A power-series method to develop third-order approximations for solutions near the triangular equilibrium points was employed.

Owing to the increasing speed and accuracy of modern computers,many significant studies on stability and bifurcations for in-plane periodic solutions were undertaken in 1965 [58].The work established an important link between families of solutions in the RTBP.At the same time,numerical solutions for in-plane periodic motion in the vicinity of the triangular equilibrium points [126,127] and vertical orbits near the collinear Lagrangian points [62] in the Earth-Moon system were obtained.Furthermore,Szebehely [33] published“Theory of orbits:the restricted problem of three bodies”;included in his study is an extensive compilation of both analytical and numerical analyses in the RTBP.

In his approach in the search for periodic orbits,Brouck used the limiting cases of families that are known by analytical considerations.There are nine known situations in which the existence of periodic orbits can be predicted analytically.Three cases correspond to the retrograde infinitesimal elliptical orbits around the collinear Lagrangian pointsL1,L2,andL3.Brouck [88] described these solutions in a general approximate initial solution.The six other solutions are all of circular type and are justified by two-body considerations.In two solutions,the distance from the satellite to the main massesm1andm2is supposed to be large in comparison with the distance fromm1tom2,so that the motion can be considered around the systemm1+m2,which may,in a first approximation,be considered as one single body.In the other four solutions,circular motion around one main body,m1orm2,is considered,and it is neglected in the other one.The approximation is valid if the radius of the circular motion is small enough compared with the distance ofm1tom2.Taking the possible motions in both directions,direct and retrograde,one has six approximated solutions.

In 1973,Henon [66,67] expanded his analysis to include vertical stability of periodic solutions in the RTBP.Farquhar [34] obtained third-order approximations for quasi-periodic motion near the trans-lunar Lagrangian point in the Earth-Moon system using a Linstedt-Poincare method.A few years later,Richardson and Cary [128]developed,via a method of successive approximations truncated to the fourth order,a model for quasi-periodic motion in the vicinity of the interior Lagrangian point for the Sun-Earth-Moon system.

A corrections scheme in 1977 that quickly computed asymmetric periodic solutions for the Stormer problem was also created by Markellos [92].The scheme proves invaluable for obtaining solutions with no visible symmetries.Furthermore,the scheme was successfully applied to the CRTBP,and the stability of the solutions was assessed using the methods developed by Henon [129].

In 1984,a three-dimensional method of regularization originally developed by Howell [76],was applied in the restricted problem and approximations developed for“almost rectilinear”halo orbits near the collinear points.A thorough numerical investigation of halo orbits was also completed [130].The improved approximations that included fourth-order terms were developed by Zagouras [77] and a correction scheme also was applied to obtain numerical results for vertical and axial orbits near the triangular equilibrium points.

6.Grid search of planar periodic orbits/constructing database

A grid-search method was presented by Broucke [88],in which the initial conditions are a first approximation for periodic orbits in the RTBP with a small mass ratio.It was found that in some cases the periodic orbits corresponding to µ=0 do not exist when the Earth-Moon mass ratio is used,and that in other cases periodic orbits exist that have no corresponding orbits in the two-body problem.Broucke systematically explored some areas of initial conditionsx0and.The values ofx0from −4 to 4 with steps of 0.1 were explored in the earlier stages of the work.Upper and lower limits for the values ofwere determined in a different way.The so-called escape or parabolic velocities were taken as approximate limits.In other words,the initial velocitywas taken such that the total energy (relative to the inert ial axes) is negative.Using the synodical initial conditionsx0,0,0,,the inertial energy of the satellite may be written in the form

whereEis restricted to be negative,and then the upper and lower limits forwere obtained.

The grid search method was first introduced in 1974[131] as a systematic way to compute complete networks of families of periodic orbits of non-integrable systems of two degrees of freedom,contained in a given region of the relevant space of initial conditions.Since then,it has been used for the numerical treatment of a number of dynamical systems,among which are,notably,the RTBP and Hill’s problem including their variants (Markellocite[131–133],Henon [98,134],Kanavos [135],Russell [136,137]).

The grid-search method [138] for the three-dimensional case is straightforward but requires a twofold extension:one in the direction of scanning through two,instead of one,variables;and another in the direction of applying it to more types of symmetric motions.Arbitrary threedimensional periodic orbit calculations are almost impossible using a grid search method when considering computational complexity.The latter extension is due to the fact that two distinct types of symmetric motions are possible and have indeed been found.A hierarchical gridsearch method was developed for systematically searching three-dimensional periodic orbits around irregular bodies by Yu [37,139–142].Kazantzis [95] presented an approach dealing only with the case of doubly symmetric motions and to the grid search applied for their discovery.The grid-search method applied for locating the initial conditions of doubly-symmetric periodic orbits in three dimensions,namely an orbit of initial stater0=(x0,0,0,0,vy,vz),is periodic of double symmetry.

Recently,Restrepo [122] presented a broad database of planar,axially symmetric three-body periodic orbits for planets and main planetary satellites in the Solar,Sun-Earth and Earth-Moon system that has been generated and made available online.The database generation is based on a grid search that incorporates a new discretization scheme centered around fundamental periodic orbit families,and a robust differential corrector implemented with a full second-order trust region method.The solutions include periodic orbits in the vicinity of the secondary orbits that circle the primary,and more complex solutions that orbit both,allowing for transitions between them.Descriptive nomenclature was developed,and a detailed characterization of the solutions used.

With the development of research and the improvement of computing power,a large number of planar and threedimensional periodic solution families can be calculated by numerical methods.However,due to the existence of numerous periodic orbits in three-dimensional space,it is impossible to calculate all the orbits by the method of random numerical search,so it must be combined with approximate analytical method to greatly improve the calculation speed and find more periodic orbits.

7.Stationkeeping of periodic orbit

Since the collinear vibration points are unstable for all restricted three-body systems,orbital control is required.The primary goal of orbital control is to maintain orbital stability of a spacecraft near a nominal orbital,which may be a Lissajous reference trajectory or a periodic reference trajectory.For the mission of Lagrange point in the Earth-Moon system,stationkeeping operation must be carried out about once a week in order to keep the spacecraft in orbit for a long time due to the unstable orbit[143].Many researchers have explored the issue of liberation-point residents in the Sun-Earth and Earth-Moon systems [144−149].There are a variety of approaches to orbit maintenance in the Earth-Moon system,but all can be generally classified as“short-term”and“long-term”stationkeeping strategies,which braced on the design ideas of a periodic orbit maintenance strategy mainly including a Target-point approach and the Floquet-mode approach [150].Both use maneuvers executed at discrete time intervals.The analysis includes some investigations of a number of the problem parameters affecting the overall cost [151,152].

The target-tracking mode,as presented by Howell et al.[153–155],which is based in Breack-well’s ideas to track the nominal orbit,guides the detector to the desired position through a reasonable design control strategy.It computes correction maneuvers by minimizing a weighted cost function.The cost function is defined in terms of a corrective maneuver as well as position and velocity deviations from a nominal orbit at a number of specified future timesti.The non-final state vectors at each timetiare denoted as“target points”.The target points are selected along the trajectory at discrete time intervals that are downstream of the maneuver.

The floquet mode,as developed by Simo et al.[156,157],based on the invariant manifold structure of the periodic orbit,maintains the divergence tendency of the unstable orbit by maintaining the orbit [150],and then the spacecraft can finally be guided to the target orbit.The Floquet module associated with the monodromy matrix is used to determine the unstable components caused by local errors.Both of these methods have been well applied in many complex multi-body system models.

7.1 Short-term stationkeeping

A loose fit reduces the total ∆Vcost of the stationkeeping corrections for ISEE-3.A precise reference orbit is not necessary,however,because the spacecraft is loosely controlled about the approximate halo path [151].The loose fit reduces the total ∆Vcost of the stationkeeping corrections.A numerical algorithm has been developed to compute a velocity correction at any point on the trajectory that will minimize the distance between the satellite orbit and the nominal halo path.

The“short-term”approach applies to the execution of short-term maintenance objectives around the Lagrange point.Farquhar [23,70] and Breakwell et al.[26,157,158]have conducted extensive studies on optimal short-term strategies,some of which are based on Floquet theory to eliminate unstable error parts.Janes [159] designed a global optimal control algorithm for the Sun-Earth systemL2.ARTEMIS mission is in operation[1],and L1 andL2quasi-periodic orbit are maintained only once for 1−2 revolutions through the direct optimization algorithm[160].It was initially known from orbit maintenance data that approximately 60 maintenance maneuvers were performed [161].Scheeres implemented the feedback law specified by the instantaneous eigenvalue and eigenvector structure of the trajectory to maintain the short-term dynamic control method [162].

A spacecraft is subject to the dynamical nonlinearities inherent in the problem when it is forced to follow a nominal path derived with linearized equations of motion.Such non-linearities were treated by Wie [163] as trajectory-dependent,persistent disturbance inputs because they are functions of spacecraft position,velocity,and acceleration.A disturbance accommodating control technique,which was successfully applied to a space-station control problem [163,164],as well as flexible spacecraft control problems [165–167],was also successfully utilized by Hoffman [168] to reduce fuel expenditure for the haloorbit control problem in an Earth-Moon-spacecraft system.It was also shown that the method can be used as an iterative method for generating a large,complex,quasiperiodic Lissajous trajectory starting with a first-order reference trajectory.The trajectory-dependent nature of the persistent disturbances,however,was not considered by Hoffman [168].

A variety of stationkeeping strategies were previously investigated for applications in the Sun-Earth system and near the Earth-Moon Lagrangian points [169].To be operationally useful for the ARTEMIS mission,a stationkeeping strategy must satisfy several conditions:utilization of high-fidelity ephemeris models,yielding of optimal solutions,and applicability to Earth-MoonL1orL2Lagrangian-point orbits and any transfer between them[170].Numerous references offer discussions of stability and control for vehicles at both collinear and triangular Lagrangian-point locations.Hoffman [168] and Farquhar[70] both generated analysis and discussions for stability and control of spacecraft in Earth-Moon collinearL1andL2locations,respectively,within the context of classical control theory and/or linear approximations.Renault and Scheeres offered a statistical analysis approach [171].Howell [172] addressed the use of Floquet theory to select stationkeeping maneuvers to eliminate unstable modes associated with a reference orbit.Gomez et al.[173] independently developed and applied the Floquetmode approach specifically to translunar Lagrangianpoint orbits.Marchand and Howell [174,175] discussed stability including the eigenstructures near the Sun-Earth locations.Folta and Vaughn [176,177] presented an analysis of stationkeeping options and transfers between Earth-Moon locations and described the use of numerical models that include discrete linear quadratic regulators and differential correctors.Based on a comparison of Lagrangian-point orbit stationkeeping methods by Folta et al.[160,178] and the imposed operational constraints,it was determined that the optimal continuation strategy offered a flexible method for maintaining the ARTEMIS Earth-Moon Lagrangian-point orbits.Fundamentally,the optimal continuation strategy is designed to maintain the spacecraft in the vicinity of a Lagrangian-point orbit for 1−2 revolutions downstream [36,179].

7.2 Long-term stationkeeping strategy

Pavlak and Howell [180] demonstrated a long-term orbit maintenance technique incorporating multiple shooting.Folta et al.[160,161,178,181],Woodard et al.[160],and Sibeck et al.[182] provided both a review of various Earth-Moon Lagrangian-point orbit stationkeeping methods,as well as detailed operational stationkeeping and transfer results for the ARTEMIS mission.

The long-term stationkeeping costs are also an important factor in determining the feasibility of such systems.In 1971,Farquhar’s [70] investigation of the use of halo orbits to maintain a continuous communications link between a lunar ground facility and the Earth also included an examination of the stationkeeping costs.Extensive work on optimal stationkeeping strategies using Floquet modes for halo orbits in the Earth-Moon system was later completed by Simo et al.[157].In the latter approach,the unstable subspace that is available from dynamic systems theory is used to develop a stationkeeping strategy.Stationkeeping analyses of Earth-Moon halo orbits were also completed by Howell et al.[153,155,183] and Gomez et al.[150,158,184].Scheeres et al.[185] and Renault [171] investigated the generalized optimal placement of statistical control maneuvers applied to orbits in the Earth-Moon RTBP.The orbits in these studies are typical of those that might be used for lunar coverage;they also provide an additional benchmark for stationkeeping costs[186].Grebow [78] adapted models created by Folta [176]in which the patch points that define a Lagrangian-point orbit from a generator are targeted in STK.

At the same time,there are some "long-term" maintenance methods,usually require the spacecraft to meet the terminal mission constraints,and at the same time,not the optimal solution.For example,the approach adopted by Grebow[78] and Folta[178] targets a strictly maintained baseline trajectory.The "long-term" maintenance method proposed by Pavlak[180] is flexible and can maintain mission orbit through multiple shooting.

8.Summaries of applications of periodic orbits

With the increasing demand for deep space exploration missions,it is very important to study the periodic orbit based on the model of three-body problem system,because the orbit is the core of rapid design and calculation of practical complex multi-body dynamic mission orbits.In the three-body problem,there are innumerable periodic orbits in the region of space.They include orbitals and all kinds of symmetry and variety.In order to make use of these orbital types in dynamical problems without analytic solutions,we must understand the spatial structure of the various orbital types as well as parameterize and classify them as much as possible.Designing practical tasks successfully and efficiently requires a new perspective and a more comprehensive understanding of the solution space,which is imperative and also a daunting task.There are seven Lagrangian points in the vicinity of the Earth:theL1andL2Lagrangian points of the Sun-Earth system and theL1−L5Lagrangian points of the Earth-Moon system[187−192].Owing to its unique spatial position and dynamic characteristics,the translation point can be used not only as an excellent position for space observation[193],but also as a transit point for interplanetary lowenergy transfer [194–196].

A summary of typical Sun-Earth and Earth-Moon system Lagrangian-point missions so far are presented in Table 3.In Table 3,the superscript ‘1’ represents Sun-Earth system,‘2’ represents Earth-Moon system,‘*’ represents cancelled mission,‘#’ represents unknown at this time.Several conclusions drawn from Table 3 are analyzed as follows.

Table 3 Applications of periodic orbit missions

First,NASA has an absolute leading role in the number of scientific missions,followed by ESA,then by countries such as China,Japan,and Russia.

Second,the scientific missions that have been launched or are currently running are mainly in the periodic or quasi-periodic orbits of the Sun-Earth system,or in theL1andL2Lagrangian periodic orbits.The number of missions in the Earth-Moon system is relatively small,but from the recent project plan it can be seen that the Earth-Moon system periodic orbit missions show an increasing trend,one that is getting closer to human life and science exploration and practice in the Earth-Moon space.

Third,each Lagrangian-point mission plan is a scientific exploration that is significant to human science and the search for the origin of life in the universe.It is particularly noteworthy that every scientific mission has had a long lifetime,from planning to final launch,some even for decades.Owing to the complexity of the mission design or the huge amount of funds to be invested,it is regrettable that many good scientific planning missions have to be cancelled after many years of development.

Finally,the mission design and control problems of the periodic orbits in the three-body system are very complicated technically.It is even necessary to consider the fourbody problem in a three-body system with strong nonlinearity [197–199].Earth-Moon space will be the first step for human beings entering deep space.In the near future,increasingly more manned space missions will be carried out in Earth-Moon space.The construction of the Earth-Moon space communications hub,navigation service system,space station,science laboratory,and habitation module will be fully considered.This will bring opportunities in the deep research and application of Earth-Moon system periodic orbits,also including near Earth technology research [200–205].At the same time,due to the extremely high reliability requirements of manned missions,it also brings higher technical challenges to the design,computation,and control of periodic orbit missions.Therefore,scholars need more systematic understanding of Earth-Moon system periodic orbit missions,which is the purpose of this article.

9.Conclusions

Recent research and studies in the field of periodic orbits of the RTBP are consolidated and deeply discussed in this article.The development of the discovery and classification of the periodic orbit,including planar and three-dimensional types,are first introduced and consolidated.This classification includes various classes of problems,such as the Poincare problem,Moulton’s analysis,Stromgren’s families,Hill’s case,and those of other three-dimensional periodic orbits.The computation methods of periodic orbits are also consolidated for all families.The stability and bifurcation of periodic orbits are discussed.Some families of periodic orbits can be discovered by their bifurcation characteristics.In addition,a modified station keeping method and its applications in periodic orbits of the RTBP are reviewed.In the end,the RTBP model is only an ideal theoretical model.Although it is very useful for human understanding of the cislunar space,it is still different from the actual mechanical model.Therefore,the orbit design research of cislunar space exploration engineering mission oriented to the actual ephemeris will be challenging research.

Acknowledgment

During the writing of this article,I would like to thank Professor Gao Yang,Master candidate Zhou Chenguang and others for their good suggestions and great help in the revision of this article.