Jiye Zhang, Huichang Yu, Honghua Dai*
School of Astronautics, Northwestern Polytechnical University, Xi’an, 710072, China
* Corresponding Author: Honghua Dai. Email:
Computer Modeling in Engineering & Sciences 2023, 135(1), 5-43. https://doi.org/10.32604/cmes.2022.022585
Received 16 March 2022; Accepted 30 June 2022; Issue published 29 September 2022
The Moon exploration has many scientific and engineering meanings, and trajectory design for the spacecraft transferring from the Earth to the Moon has been the subject of countless works. The conventional transfer trajectory design method began with the Hohmann transfer (Fig. 1) and was improved using the patched conic model (Fig. 2). With the development of the dynamical system, the low-energy ballistic capture trajectory (Fig. 3) was proposed to optimize the fuel requirement. Completely different low-thrust trajectories (Fig. 4) were further developed for other available propulsion systems, like an electric propellant or a solar sail, significantly enhancing the payload capability and enabling high ΔV missions.
The primary transfer approach to the Moon can be broadly classified into four categories: direct transfer, phasing loop transfer (also known as direct, staging), low-energy transfer, and spiral transfer . The direct Moon transfer trajectory is the most commonly used method in engineering tasks . In the framework of the two-body problem, the most straightforward direct transfer trajectory can be obtained by the Hohmann transfer and later a modified patched conic model. Despite the enormous energy requirement and difficulties in searching for proper patching points, the patched conic method was favored in many lunar missions, such as Apollo, Luna, and Chang’E missions, because of its short flight time and simple flight control system. The total ΔV acquired transferring from a 300 km circular Earth parking orbit to a 100 km circular cislunar orbit ranges between 3.5 and 4 km/s. And the typical flight time of a direct transfer varies from two to five days. In general, the most fuel-efficient direct transfer needs about 4.5 days. Any more time-consuming travel typically sends the spacecraft beyond the Moon’s orbit before it falls back and encounters the Moon .
In the phasing loop transfer, the trans-lunar injection or the lunar orbit insertion burns are divided into several minor burns to minimize gravity losses. The initial direct and phasing loop trajectories are designed through either the Hohmann transfer or patched conic technique based on an analytical solution to the two-body problem. This approach can provide a chance to verify the operating condition, address the status of the orbiter and correct any anomalies before arriving at the Moon which can be seen in the Chang’E-1 and Chandrayaan-1 missions .
However, the analytical two-body model inherently involves hyperbolic approaches upon arrival, which requires a high ΔV maneuver to inject the spacecraft into the final Moon orbit . With the development of the circular restricted three-body problem (CRTBP), some studies showed that hyperbolic velocity is not a necessary condition for a cislunar transfer . By exploiting the concept of temporary ballistic capture, or weak ballistic capture theory, low-energy transfer  was proposed to reduce the total energy cost and has been executed in the rescue of the Hiten mission [7,8], and, more recently, in the ARTEMIS  and GRAIL missions . For a transfer from a 167 km circular Earth parking orbit to a 100 km circular cislunar orbit, the low-energy transfer yields improvements over Hohmann by 18% but the flight time is quite long up to several months.
Currently, Earth-Moon low-energy transfer trajectories are mainly constructed based on the weak stability boundary (WSB) theory or Lagrangian points and manifold theory [8,11]. According to the trajectory geometry, the low-energy transfer can be classified into the interior (or cislunar) and exterior (or translunar) [5,12]. Trajectories are basically located within the Moon orbit in the interior transfer and are always designed through the Earth-Moon system L1 Lagrangian point or its periodic orbits. As the Lagrangian point and the hyperbolic invariant manifolds of periodic orbits are approached asymptotically, an interior transfer via invariant manifold usually has a pretty long flight duration [13–16]. Moreover, because of the very long flight time, the solar gravity perturbation would cause a significant impact, which can change the invariant manifolds topologically to fail in transiting the L1 point or the periodic orbits . For the exterior transfer, the spacecraft is injected into highly ellipse orbits whose apogee at approximately four Earth-Moon distances. In this region, a small maneuver, combined with the solar perturbation enables an approach to the Moon from the exterior. An exterior transfer could be constructed by using the WSB theory [8,17,18] or the patched manifold technique [6,11,19], which requires more fuel in the accelerating maneuvers but less in the decelerating maneuvers, making it more economical than the direct transfer. However, this kind of transfer heavily depends on the geometry of the Sun, Earth, and Moon. Compared to the interior transfer using the Earth-Moon L1 point, the exterior transfer has fewer opportunities to transit the L2 point. An effective way to increase the transfer opportunities is to transit a periodic orbit near the L2 point instead of the point itself .
The spiral approach, a type of low-thrust transfer, is the one subjected to the longest flight time of all methods. The SMART-1 spacecraft reached the Moon following a low-thrust 2-year trajectory, becoming the first successful case using this approach to the Moon . The trajectory, in this case, can be divided into three segments: escaping from the Earth, free flight segment, and capturing by the Moon. The probe can achieve Earth escape and Moon capture through the cumulative effect of low thrust. Defining an (M, N) -loops transit orbit of the ballistic capture, the Moon would capture the probe after it runs M circles from the Earth’s large ellipse with a large semi-major axis and N circles on the Moon’s large ellipse orbit. The whole process is free-flight movement. To this end, the essence of small thrust is to search an (M, N) -loops transit orbit to extend the arc of low-thrust action. Moreover, due to the high energy level of the patched orbit in an exterior transfer, a spacecraft cannot follow an elliptical orbit with a sizeable semi-major axis before escaping from the Earth, which makes it challenging to apply the low-thrust in the exterior transfer. In contrast, the interior transfer based on the L1 point has a stable (M, N) -loops transit orbit, which can make full use of the low-thrust to accumulate energy for a long time to achieve gradual orbit change . A comparison of these four types of transfer is listed in Table 1.
Several review papers are available that discuss the recently proposed low-energy and low-thrust transfer methods and the complex mechanisms of the trajectory design. However, there is a lack of a review enabling readers to grasp the overall development and design principles of the Earth-Moon transfer trajectory construction from a historical and scientific perspective. This paper begins by demonstrating the traditional transfer approach, containing the patched conic and pseudostate methods. Then the low-energy transfer approach is illustrated from two aspects: designing via manifold theory and weak stability boundary theory. Finally, the design of the Earth-Moon transfer trajectory by using the low-thrust technique is also discussed. With no more complex theories and techniques in this paper, essential mechanisms and characteristics of these design methods are discussed to help readers easily grasp the basic overview of the current Earth-Moon transfer orbit design methods.
The traditional direct transfer is the most commonly used method in practical engineering tasks. The mainstream design methods to construct this type of transfer, the patched conic method, and the pseudostate technique have been proposed since the era of Apollo and achieved significant development in the past few decades. The patched conic model features a good analytical characteristic due to the multiple-two-body assumption. However, the approach of neglecting the gravitational force of the central celestial body successively may cause non-negligible errors at the patched point. A modified method, i.e., the pseudostate technique was then proposed, whose accuracy is between the patched conic and numerical integration models, owning better practicability in practical engineering missions. In this part, the construction of a direct Earth-Moon transfer trajectory based on these two methods would be reviewed.
In 1925, Hohmann  discussed the double-impulse transfer method between two circular coplanar orbits and concluded that the minimum energy transfer conic would be elliptic, tangent to the inner and outer orbits at perigee and apogee, respectively. When designing an Earth-Moon transfer orbit, the Hohmann transfer can only be used to design a primary trajectory because of its too strict condition, maneuvering at the perigee and neglecting the Moon’s influence on the spacecraft. In 1956, Egorov  proposed the patched conic method in the frame of the two-body problem. By introducing the concept of sphere of influence (SOI), a multi-body system can be divided into several two-body systems. Each portion of the trajectory is treated as a two-body trajectory about the most influential body . Then at the boundary of SOI, a transition is made from one central gravity body to the other by matching the relative position and velocity in the appropriated conic trajectories. One of the most commonly used SOI models is the Laplace sphere . For the Earth-Moon system, the radius of the Moon’s SOI can be calculated by
where, m and M represent the mass of the Moon and the Earth, respectively, and D denotes the average distance between Earth and Moon which follows that RS = 66300 km.
A general interstellar transfer trajectory in the solar system can then be divided into three segments based on the above assumptions : (i) a hyperbola about the departure, (ii) an ellipse about the Sun, and (iii) a hyperbola about the target planet. In the case of the Earth-Moon system, the whole trajectory can be simplified into two segments, i.e., the geocentric segment and selenocentric segment since the Moon’s SOI lies within the Earth’s SOI .
Many studies have been carried out based on the patched conic method from the era of Apollo. Penzo  impelled the patched conic method further in both theory and computing techniques. He also proposed a free-return orbit design scheme, and determined the effect of the perilune distance and outward/return inclinations to the Moon’s plane on the flight time . It is concluded that the possible launch conditions would be constrained when the perilune distance, flight path, and landing location are fixed. Fig. 5 shows a free-return orbit studied in . It can be seen that the geocentric orbit is an ellipse orbit, and the selenocentric orbit is a hyperbolic one relative to the Moon. An interesting phenomenon is that the SOI entrance velocity is not equal to the SOI exit velocity, due to the velocity combination effect. The phenomenon is also known as the gravity-assist through Moon swing-by , which can be further used in an interstellar transfer mission by accelerating the spacecraft.
Moreover, Dallas  discussed the characteristics of Earth-Moon free-return orbits in more detail. By analyzing the flight time, fuel consumption, and the difficulty of orbit design, he concluded that the counterclockwise class of trajectories is superior to the clockwise ones. To verify the accuracy of the patched conic model, Gibson Jr.  compared several parameters (e.g., the velocity (V) and lead angle (ϕ) at injection) between the patched conic model and the 7-body model (considering the gravitational forces of Earth, Moon, Sun, Venus, Mars, Jupiter, and perturbations due to the non-spherical gravitational fields of the Earth and the Moon). The results showed that the patched conic elements can approximate precision elements adequately. He also found that an even better approximation may be obtained by adding constant correction terms to the conic element. In addition, he proposed an iteration scheme that combines the patched conic model in conjunction with the integrated n-body model to obtain any desired integrated trajectory, and a fast convergence can be achieved. Li et al.  proposed a multi-segment lunar free-return trajectory when studying the design of a crewed spacecraft's Earth-Moon orbit. They established an analytical model of the multi-segment free-return trajectory using the patched conic technique. The results demonstrated that multi-segment free returns allow a wider range of lunar approach inclinations and free launch windows than free-return profiles. They also found that retrograde Earth returns generally require a higher velocity impulse but lower flight time than posigrade Earth returns which can also offer a wider range of lunar approach inclination.
Salazar et al.  studied the problem of determining the optimal flight duration and energy consumption vector for a spacecraft transferring to the Moon. As shown in Fig. 6, λ1 denotes the Moon arrival angle, and the minimum time of flight (TOF) occurs when λ1 = 0° and the minimum total cost (ΔV) occurs when the λ1 = 44°. Qi et al.  investigated the two-impulse transfer from the 167 km LEO to the permanent lunar capture orbit with the minimum ΔV in the framework of the planar CRTBP and the two-body model. Different from the traditional method, the Moon ellipse target orbit rather than the circular orbit was used for the sake of less fuel consumption. As an important optimization variable, the height of the perilune of the hyperbolic orbit h determines the optimal transfer and the capture point. As illustrated in Fig. 7, ΔV generally decreases and is followed by a rise with the increase of h. The minimum of ΔV occurs at a height of around 1200 km. The TOF showed a different trend of monotonic increase as h increases. However, the basic assumptions of patched conic, i.e., (i) the Moon acting in a circular motion around the Earth, and (ii) neglecting the gravity of the central body successively, can cause considerable errors. Hence, a more accurate dynamical model containing ephemeris needs to be considered for the practical Earth-Moon mission [35,36]. It should also be noted that, in operations of practical engineering missions, midcourse impulses are always employed to eliminate the initial error in entering the Earth-Moon transfer orbit. For the Chang’E-1 mission, two essential midcourse impulses were planned to perform at the 17th and the 90th h, respectively. An alternative was arranged at the 41st h, which would be performed when the residual error exceeds a certain upper bound after the first midcourse correction. However, due to its ideal flight condition, only one alternative midcourse correction was given at the 41st h of flight duration, and the speed increment used was only 4.8 m/s . Moreover, an important method, the phasing loop method was utilized to link the Earth-Moon transfer orbits in the Chang’E-1 mission. Concretely, the lunar probe ran three circles in the super GTO orbits and operated a small maneuver in the second apogee to increase the height of the perigee to 600 km. The later three impulses were employed at the perigee, changing the orbital period to 24 and 28 h, respectively. A final impulse was given to direct the spacecraft into 114 h Earth-Moon transfer orbit . This method has been proposed since the 1990s which is suitable for limited engine thrust and can reduce midcourse correction increments and widen the launch window [39–41].
The patched conic method has a good analytical characteristic, giving a more clear and intuitive view in designing an Earth-Moon transfer trajectory . However, the basic assumptions adopted can cause large errors in both position and velocity at the patched point. Thus, an iterative determination of the intersection at the SOI is necessary, which results in a waste of computing time and instability of the algorithm . In 1970, Wilson  proposed the pseudostate model in the study of the Earth-Moon transfer missions, which reduced patched conic errors by 80% without employing a search enabled by numerical integration. As the basis of the pseudostate model, the overlapped conic approximation was firstly derived as a closed-form approximate solution for the restricted three-body problem. For the direct transfer missions, the rectification terms in overlapped conic equations can be omitted to achieve the desired long propagation step. The model then can be simplified into the pseudostate model. Similar to the SOI concept in the patched conic method, a concept of pseudostate transformation sphere (PTS) is introduced in the pseudostate model to define the vicinity of the Moon . Within the PTS, the superimposed effect of the Earth and the Moon is considered. Outside the PTS, only a simple geocentric conic is used .
For an Earth-Moon transfer trajectory, the state of the spacecraft’s position at its start and end can be approximated by two instantaneous conics arcs. And the pseudostate theory is then utilized to find a bridge connecting the instantaneous conics at both ends. Suppose the spacecraft’s geocentric states of motion at an initial time and an arbitrary time are and , respectively. The selenocentric state at is . Let be the geocentric state of the Moon at . The relationship between and can be represented by and .
1. Geocentric instantaneous conic: Ignoring the gravitational force of the Moon, then the geocentric motion can be viewed as a two-body problem about the Earth. The pseudostate is obtained by propagating from to .
2. Linear backward propagation section: Calculate the selenocentric state corresponding to . The can be obtained by linearly backward evolution to (the PTS spherical moment of arrival) through and .
3. Selenocentric instantaneous conic: The perilune state can be obtained by propagating the in the central gravitational field of the Moon from to . Then the geocentric state can also be acquired by its relationship with .
Based on the pseudostate theory, Luo et al.  proposed a preliminary design algorithm for the free-return lunar flyby trajectories. It is a simple iterative algorithm that does not require any prior information and can generate a better initial solution for the free-return lunar flyby trajectory. Zhang et al.  proposed a linear approximation-based fast prediction algorithm based on the improved pseudostate theory for the Earth-Moon transfer trajectory. Simulation results showed that the prediction error of this algorithm is 10% of that of the multi-conic technique, and 1% of that of the pseudostate theory. And the algorithm is about 10 times faster than the numerical integration.
Byrnes et al.  proposed a multi-conic method that is able to consider a number of complex perturbations in the whole travel. By using a relatively large time step, the method is fast enough and more accurate than the traditional patched conic or pseudostate technique. In this method, the transfer orbit is decomposed into a series of segments and each part could be solved with a similar pseudostate technique. They found that this method’s accuracy is an order of magnitude higher than that of the pseudostate technique, the computation time is two orders of magnitude faster than that of the numerical integration method, and the error is only 1%~5% of that of the patched conic method. For circumlunar flight, errors of the calculated pericynthion and return perigee altitude could be around 20 km; velocity errors be about 0.5 m/s . Byrnes  integrated the multi-conic technique into a one-step method, which can be used in a manner analogous to the Lambert method to solve the transfer initial value problem between two planets in an N-body environment. Unlike the multi-conic method, the one-step method associated the pseudostates with the departure and arrival planets through Lambert conic, where the secondary gravity force fields outside the central force field were assumed to be active for a certain flight time. Damario et al.  had also shown that such Lambert-like applications of the one-step method may be extended to a flyby problem by combining two such single-step transfers, placed on either side of the flyby celestial body.
To simplify the procedure of the one-step method and improve its computational speed. Sergeyevsky et al.  proposed a one-step impact pseudostate technique for direct interplanetary transfer by assuming that the planetary center departure and arrival phase trajectories are rectilinear hyperbolas. This method can avoid the iterations required to obtain approximate hyperbolic trajectory features and the computation of the state transition matrix. Ramanan  applied the one-step impact pseudostate technique to generate the one-way Earth-Moon transfer trajectory and developed a comprehensive algorithm. In his paper, the pseudostate was determined by assuming a rectilinear hyperbolic for the arrival phase while the translunar injection state was obtained by solving a Lambert problem. The algorithm can reduce the error of the spacecraft to the target position by more than 95%, which provided a better preliminary value relative to the Lambert conic method with minimal computational effort. However, it is difficult to locate the hyperbolic orbital plane with paralleled position and velocity vectors of the rectilinear hyperbolic employed in the aforementioned algorithm. An undesirable collision with the target could occur in an actual mission scenario. Ramanan et al.  developed the one-step pseudostate technique for designing a non-impact Earth-Moon direct transfer trajectory. A simple matching strategy was employed to synchronize the geocentric transfer trajectory with the selenocentric trajectory with no need for computation of the state transition matrix. Characteristics of the nonimpact approaching trajectory were computed from the excess velocity vector. Furthermore, a semi-analytical based targeting technique was proposed to obtain fast and near exact solutions in the frame of three-body problems.
Moreover, it should be noted that the accuracy of both the pseudostate technique and the one-step method depends on the sweep-back duration which is directly relative to the size of PTS. On the one hand, PTS with a small size corresponds to a short sweep-back duration, leading to premature neglect of the secondary gravity force field. On the other hand, too long a sweep-back duration can cause poor approximation accuracy. The selection of the sweep-back duration or the radius of the pseudosphere is vital for transfer trajectory design with the pseudostate technique . In the context of the interplanetary transfer problem, Sweetser  and Kiedron et al.  established an experienced relationship through multiple simulations and obtained the optimal value of this parameter. The duration of the lunar mission flight is only 3–6 days, so the optimal value can be easily obtained by grid search where a zero sweep-back duration corresponds to the solution of the Lambert conic method . Results show that the achieved periapsis time after propagation is very close to the expected periapsis time when choosing the sweep-back duration of about 75% of the total flight duration. And the optimal percentage value changes with the day of departure.
Classical approaches to design the direct transfer from the Earth to the Moon such as Hohmann transfer, patched conic, and pseudostate methods are introduced in the previous section. Although it is desirable to deal with the analytical solution, the fuel consumption is quite high which can render missions infeasible in some circumstances . Fortunately, a new class of low-energy trajectories has recently been discovered and employed, which makes it possible for missions that could not be realized by classical approaches. The low-energy transfer is preferred to direct transfer for its appealing features, such as lower cost, the extended launch window, and the property of allowing for a fixed arrival epoch for many launch dates within the launch window [56–58].
In this part, we will focus on Earth-Moon transfer trajectory design based on Lagrangian point and manifold theory. The trajectory is classified by the Lagrangian point being used. The transfer through the Earth-Moon L1 point is usually located within the Moon’s orbit, and hence it is called interior or cislunar transfer. The Earth-Moon L2 point can also be used to construct the interior transfer, but because the energy of the L2 point is greater than the L1 point, it is not optimal when designing a cislunar transfer trajectory. However, recent studies [6,11,19] showed that with the help of the solar perturbation, the L2 point can be used to design the optimal exterior transfer which is also known as the translunar transfer. Before systematically analyzing and summarizing these two low-energy transfer methods and their development, the fundamental of the problem, the circular restricted three-body problem (CRTBP) will be introduced briefly.
The three-body problem dynamic model is attracting ever-growing attention from both academia and industry as its excellent accuracy compared with the traditional two-body and the perturbed two-body models . A restricted three-body problem was then proposed, which can be applied well to practical space missions. A modified model is the circular restricted three-body problem (CRTBP), which provides a high fidelity approximation to the spacecraft motion and can underline the dynamic nature well . In the CRTBP, a spacecraft is considered to be massless. Assume that two primaries are mass points and make circular motions around their common barycenter. The spacecraft’s motion is then described in a rotating frame relative to the two primaries. The non-dimensional motion equations of the third body can be written as 
where the subscripts denote the partial derivatives of the function:
with the mass parameter
Five equilibrium points also known as the Lagrangian points L1, …, L5 exist in the CRTBP system, where the gravitational and centrifugal forces are balanced on the spacecraft. For the Earth-Moon system, the mass parameter μ = 0.01214. The Lagrangian points of the Earth-Moon system are summarized in Table 2.
In 1836, the Jacobi constant, C, was discovered under the rotating frame by Jacobi, which is the only constant of energy integral presented in the CRTBP so far . The value can be determined from the initial conditions of the spacecraft’s motion:
For a given value of C, the curves of zero velocity are determined by
By specifying the spacecraft energy equal to the value at each Lagrangian point, the curves of zero can be plotted as Fig. 9 shows.
It is shown that a spacecraft with C equal to that of the L1 point can cross over from Earth to Moon. The one with C of the L2 point can cross the boundary to the exterior, and with C of the L4/L5 point can reach any point in the Earth-Moon system . Moreover, there are abundant periodic and quasi-periodic orbits near the Lagrangian points in the Earth-Moon space, such as Lyapunov orbits, halo orbits, Lissajous orbits and distant retrograde orbits (DROs), etc. [64–66]. For more detail, one can refer to reference , where a comprehensive review of periodic orbits in CRTBP has been provided.
In 1968, Conley  investigated the phase space structure near Lagrangian points and categorized the motions as (1) periodic orbits, (2) stable/unstable manifolds of periodic orbits, (3) transit orbits, and (4) non-transit orbits, these orbits are shown in Fig.10. Additionally, he thought that the invariant manifolds of a periodic orbit separate transit orbits from non-transit orbits, and the former can be used to construct low-energy cislunar transfer trajectories. McGehee  investigated the dynamical behavior of CRTBP and indicated the role of the stable and unstable manifolds of the Lyapunov orbit in understanding the transiting trajectories. To acknowledge the contributions of these authors, Marsden et al.  denoted the invariant manifolds as Conley–McGehee tubes (also known as C-M tubes). Yamato  then studied the dynamical model considering perturbations and illustrated that most of the tubes would be distorted but few of them could be preserved by small perturbations from the perturbed gravitation of the third celestial body.
In the two-body-based classical approaches, direct transfer always needs high energy (C < 2). Theoretically, C < CL1 = 3.2003 is enough for an Earth-Moon transfer, with which a spacecraft could target the neck region around the L1 point. For the planar CRTBP, a planar periodic would separate the interior of the Moon’s region from the rest of the Earth-Moon region. However, the conception above only promoted the possibility of low-energy transfer in theory, and the transfer time tends to be infinite. Based on the author’s knowledge, the design method can be divided into three ways according to the used design tools, i.e., design with (i) chaotic control, (ii) manifold theory and (iii) numerical optimization.
As a systematic and classical method, the design of Earth-Moon transfer based on chaos has been studied since the era of Poincaré . However, very few comprehensive reviews on it can be found. A pioneering work done by Ott et al.  on the later developed chaotic-motion-based chaos control theory and illustrated that a small initial perturbation might cause a huge control effect in a chaotic system. The method was come up with to address the local issue of stabilizing the trajectory on unstable periodic orbits embedded into the chaotic system . A further development enables it to solve the global issue and construct a trajectory passing the desired target in a chaotic region known as ‘targeting’, consisting of finding a nearby short trajectory to satisfy the terminal conditions . However, achieving ‘targeting’ in chaotic Hamiltonian systems is always challenging because of their complicated structures .
Many efforts have been made in designing chaotic transfer orbits to the Moon via targeting. In 1995, Bollt et al.  first obtained a chaotic Earth-Moon transfer orbit which the Moon with this method can ballistically capture. They applied a sequence of small perturbations to the initial condition of the spacecraft motion. Setting the Jacobi constant of C = −3.17948, slightly higher than CL1, and employing eight tiny perturbations, a low energy transfer orbit can be constructed linking with two very high altitude circular orbits of hE = 53,291 km and hM = 12,232 km. Fig. 11 shows a transfer orbit resembling the figure demonstrated by Bollt . The total ΔV = 749.8 m/s whereas the flight time was up to 748 d. To reduce the flight time, Schroer et al.  modified the method. They obtained a low-energy transfer orbit that required at least four impulses by targeting a series of unstable periodic orbits. This approach can decrease the flight time to 377.5 d with a ΔV very close to the value in . Macau et al.  designed a transfer with a slightly higher ΔV of 767 m/s and considerably cut the transfer time to 287 d. Ross et al.  obtained a family of solutions with the best requiring of ΔV = 860 m/s, reducing the transfer time to only 65 d.
Recently, Salazar et al.  extended Bollt’s method by increasing the pre-assigned distance between two points in the Poincare section. Needless recurrent loops can be removed, and thus, total flight time can be considerably reduced. Moreover, the transfer to triangular equilibrium points was explored in the study. Some studies also combined the targeting method with an optimization approach, which treated the chaotic transfer orbits as a class of multi-constrain fuel problems with multi-dimension decision variables [72,77,78]. With this approach, periodic orbits or random mass search is unnecessary, significantly reducing the transfer time and effectively constructing the low-energy transfer. For example, Wang et al.  presented an elitist teaching-learning-based and optimization-based optimal targeting method that can rapidly direct the chaotic series toward the Moon. The optimal transfer orbit they obtained was ΔV = 748.6 m/s and Δt = 213.5 d. Additionally, the trajectory can also be constructed by a hybrid approach using chaos control and invariant manifold simultaneously .
The above studies demonstrate that low-energy Earth-Moon transfer is possible when the Jacobi energy of the spacecraft is just a bit higher than CL1, although a long flight time is required. It is worth noting that the low-energy orbit constructed is to link two extremely high looping orbits since the region of chaos could only exist in the vicinity of such two high orbits. This is easy to understand from the invariant manifold perspective. Suppose backward integrate the stable manifold of the L1 point to intersect the low Earth orbit (LEO). In that case, only one impulse from the LEO orbit is needed for sending the spacecraft into the transfer orbit wrapped by the stable manifold. Unfortunately, the stable manifold of the Earth-Moon L1 point cannot reach a position very close to the Earth (about 0.1 Earth-Moon distance). Hence, Conley  in 1968 proposed a natural idea of linking the L1 transit orbits with two “looping” orbits around the Earth and the Moon through an impulsive provided at the crossing point. The spacecraft was supposed to transfer from the Earth to the two-dimensional Lyapunov orbit and the Moon’s inner region. This work has laid a theoretical foundation for the low-energy transfer construction by utilizing the Lagrangian point and its invariant manifold. Pernicka et al. then  defined a low-energy Earth-Moon transfer from a 167 km Earth parking orbit to a 100 km Moon circular orbit into two parts. One was from the L1 point to lunar orbit, and another from LEO to the L1 point. He applied AVS at the perigee/perilune location to control the orbit period and found a transfer with ΔV = 3824 m/s, Δt = 292 d. Topputo et al.  put forward an approach of patching conic-manifold to design a low-energy transfer orbit linking the two same target parking orbits as Pernicka et al.  studied. The transfer trajectory was divided into three parts: Earth-L1 leg, transit orbit, and L1-Moon leg. The connection arcs, i.e., the Earth-L1 and L1-Moon leg, were modeled by the three-body Lamber problem and solved separately with the two-point boundary problem. In their study, the amplitude of the hyperbolic transit orbits near L1 was also studied. A patched transfer orbit from a hE = 167 km Earth orbit to a hM = 100 km Moon transfer was found requiring ΔV = 3894.9 m/s and the transfer time Δt = 255.5 d. Again assuming the same two target orbits, Mengali et al.  developed a two-impulse optimization method. First, a closed-form approximate expression for the total velocity variation was deduced under the assumption of minimum bi-impulse maneuver. For a given Jacobi constant, a transit point was selected whose abscissa coincided with the L1 point. Then the amplitude of hyperbolic transit orbits above L1 and the velocity angle was taken as optimization variables. A hybrid methodology of combining genetic algorithms with a deterministic simplex method was employed to refine the solution. The best results of ΔV = 3861 m/s and Δt = 85 d were achieved with the fuel consumption and TOF almost equivalent to those obtained based on WSB theory.
Pontani et al.  proposed a geometrical approach for low-energy Earth-Moon mission analysis. They illustrated how isomorphic mapping could be profitably employed to optimize specific cislunar transfers by geometrically determining the desired optimization parameters. Liang et al.  systematically discussed both direct and low-energy cislunar trajectories in the context of restricted three-body dynamics. Different transfer trajectories were assessed by considering Jacobi energy and time of flight. They also illustrated the relationship between Jacobi energy, perilune/perigee radius, and the total flight time. The two-body energy with respect to the Moon at perilune and global sketches were accomplished. Xu et al.  explored the low-energy transfer of the Earth-Moon equivalent L1 and L2 points under the spatial bi-circular model (SBCM). They divided the transfer into two types: the cislunar transfer transiting the L1 point or its halo orbit, and the translunar transfer transiting the L2 point or its halo orbit. Though the solar perturbation in this model plays a secondary role compared with the Hamiltonian value, it turns the transfer time targeting the L1 point from infinite into finite and makes the low-energy transfer based on a time-dependent manifold possible. They illustrated that the L1 periodic orbit could increase the transfer opportunity by introducing another variable, i.e., serial points of a periodic orbit, but this also means a higher energy consumption. Tan et al.  investigated low-energy bi-impulsive transfer using periodic orbits. Two strategies of the direct and indirect methods were proposed. For the direct method, the large amplitude Lyapunov orbits were used to provide the initial value and approximate flight time. The indirect way aimed to design a bi-impulsive low-energy lunar flyby trajectory since lunar gravity can reduce the velocity increments of the Earth-Moon transfer. Concretely, two reference orbits were employed in this strategy, one is the Earth-Moon transfer trajectory which has been designed by the direct method and the other one is a candidate periodic orbit. The basic idea was to patch these two orbits together through differential correction. Both methods were extended into the bi-circular model (BCM) and solved by the optimization algorithm. The results showed that solar perturbation has the potential to further reduce the total cost of the bi-impulsive transfer. Moreover, considering the solar perturbation, Sousa-Silva et al.  designed a cislunar transfer by patching the Sun-Earth quasi-periodic orbit with the stable manifold of Earth-Moon L1 periodic orbit which can achieve a faster Earth-Moon transfer with ballistic capture.
Recently, an emerged conception of Earth-Moon transit station in cislunar space such as the Deep Space Gateway brings a lot of attention to the stable periodic orbit in the cislunar space, such as near-rectilinear halo orbits (NRHOs) and distant retrograde orbits (DROs) . Therefore, a brief introduction is also given to the transfer problem based on these two kinds of orbits in the Earth-Moon space. Both NRHOs and DROs have perfect stability, but they also have distinct characteristics. The DROs are planar while the NRHOs are highly inclined. The DROs can cover the Moon entirely within the lunar orbital plane and the NRHO has a small minimum distance to the Moon . Therefore, the NRHOs and DROs have different roles in future space missions.
As for the transfer problem related to the NRHOs, not only the transfers between Earth and NRHOs but also the transfers from NRHOs to the lunar surface and lunar orbits were studied [89–91]. Because of the stability of DROs, the hyperbolic manifolds commonly used in the design of low-energy transfer orbits no longer exist for the DROs. Thus, it is necessary to investigate the transfer orbits based on DROs specifically. For the transfer orbits from the Earth to DROs, Ming et al.  designed an indirect transfer to a DRO around the Moon by exploiting the Lyapunov orbits tangential to the DRO. Capdevila et al.  presented more solutions for direct transfers from the Earth to the DROs and indirect transfers by Lyapunov orbits associated to L1 or quasi-periodic DROs as bridges. Moreover, Capdebila et al.  also obtained direct or indirect transfer orbits based on the optimization approach. For the transfer from DROs to the Moon. Zhang et al.  studied the transfer orbits from DROs to low lunar orbits with inclinations between 0° and 90° in both CRTBP and planar bicircular restricted four-body problems. Peng et al.  utilized the numerical method to acquire low-energy transfer trajectories to DROs with multiple powered lunar flybys (PLFs) and WSB-like ballistic transfers. The results showed that a single PLF and a WSB-like ballistic arc could further reduce the total cost at the expense of time of flight. The ΔV-Pareto optimal solution could be enhanced and more feasible solutions for low-energy transfer and rendezvous trajectories could be obtained by the double and triple PLFs.
In addition, some studies have been carried out to find multi-impulse (mainly bi-impulsive) transfers under either the three-or four-body model and focused much on dealing with the optimization problems [97–99]. Since our review paid more attention to the characteristic of the dynamical models, relevant results would not be detailly introduced here. Comprehensive reviews of the topic can be seen in [4,5,100].
The L2 point in the Earth-Moon system can also provide opportunities to construct low-energy trajectories. However, different from the only type of cislunar transfers based on the L1 point, the transfer trajectories design based on the L2 point contains both the interior cislunar ones and the exterior translunar ones. The former is essentially the cislunar transfer trajectories transiting the Earth-Moon L1 point through a heteroclinic connection which costs more fuel than the cislunar trajectories merely using the L1 point. The latter aims to take full use of the solar gravity and be captured by the Moon ballistically, which has the same geometrical shape in the inertial frame as Belbruno’s result [8,20]. Therefore, considering the former cislunar transfer by using the L2 point is not a fuel-efficient Earth-Moon transfer, the emphasis will be put on the exterior translunar transfer in this section.
The reason why the exterior transfer trajectories can be constructed is the potential of the Sun-Earth system’s stable manifold reaching very close to the Earth enabling a spacecraft to be sent directly to the manifold from the near-Earth space via a single impulse. Then, another impulse is conducted at the transition of the manifold to make the spacecraft enter the stable manifold around the L2 point and be captured by the Moon. This type of transfer trajectory was firstly given by Belbruno  based on the Weak Stable Boundary(WSB) theory whose detail will illustrate in the next part. Later in 2001, Koon et al.  explained the dynamic mechanism using the invariant manifold theory. Enlightened by the patched conic method to construct the direct transfer trajectories, they developed the patched three-body approach to take advantage of the dynamical structure of the CRTBP to obtain preliminary solutions to the Sun-Earth-Moon-Sc restricted four-body systems. And the foundational idea was to use two solution arcs to connect an initial geocentric orbit to a final selenocentric: one is a non-transit orbit associated to the Lyapunov orbit around L1,2 of the Sun-Earth system; another is a transit orbit associated to the Lyapunov orbit around L2 of the Earth-Moon system. By introducing the tool of the Poincaré section, potential patching points could be found, when integrating forward, the spacecraft will be guided by the L2 Earth-Moon manifold and get ballistically captured by the Moon; when integrating backward, the spacecraft will hug the Sun-Earth manifolds and return to the Earth. The transfer trajectory is shown in Fig.12.
The patched manifold methods introduced above are based on three-body models. However, an expected trajectory is under the four-body or the ephemeris model. Therefore, an effective method is recomputing the Lyapunov orbits and their manifolds in the four-body model to reproduce the low-energy transfer orbit . In addition, the patch point selection was always achieved by the trial and error method by which an optimal solution is not a guarantee. Peng et al.  presented an improved DE algorithm and found an optimal patch point. The solution obtained is 92 m/s less than the Hohmann transfer. Gong et al.  studied the low-energy Earth-Moon transfer considering the angle between the ecliptic-and lunar orbit planes. They derived the analytical expression describing the relation between the initial ΔV at the parking orbit, the parking orbit’s altitude h, and the Lyapunov orbit’s amplitude Ax. A curve of ΔV to h of the patched manifold was also provided. Compared with the Hohmann transfer, a 20% reduction in ΔV was achieved (3285 m/s vs. 3960 m/s) in a transfer departed from 200 km Earth parking orbit and captured by the lunar ballistically. The flight time in this case is 300 d. However, the utilized periodic orbits were still the planar Lyapunov orbits, and the corresponding invariant manifolds were two-dimensional. Fantino et al.  investigated the transfers between connections. Among the four combination types, only and can provide low or zero-cost (if there is a heteroclinic connection) transfers. Moreover, the points of the invariant manifold trajectories, of which the radius vector is orthogonal to the velocity vector of the smaller primary, can effectively confine the unstable points of the WSB region, revealing the relationship between the invariant manifold and the WSB. They also studied the second primary’s capability of temporary capture and illustrated the efficiency of the case when the Jacobi constant of the invariant manifold is large and the amplitude of the periodic orbit is small.
The above studies show that the exterior transfer still costs a long transfer time, usually more than 90 days. This can be attributed to the slow movement of the stable/unstable orbits on the hyperbolic invariant manifolds of Lyapunov orbits, by which the transfer trajectories followed closely, towards the periodic orbits as time goes to plus/minus infinity. To construct a faster Earth-Moon transfer orbit, Sousa-Silva et al.  proposed a new strategy that connected the quasi-periodic orbits with the hyperbolic manifold-guided solution. As shown in Fig.13, with the same energy level, the black lines represent the non-transit orbits that follow closely or loosely and before reaching the selected Poincaré section after about 80 to 300 days. The blue lines are quasi-periodic orbits on two-dimensional tori around the Earth that reach the Poincaré section after less than 40 days. This alternative solution can reduce the total transfer time to about 10 days, while still providing lunar ballistic capture at arrival . However, it is worth noting that the selected departure orbit was a GTO orbit whose energy is considerably higher than that of the traditional Earth parking orbit. Finally, they also explored the superiority of the new method to the conventional direct Earth-Moon transfer method in extending the launch window and decreasing the midcourse maneuver.
Since solar perturbation is essential in constructing a low-energy exterior transfer, applying a restricted four-body model rather than a three-body one is ideal for achieving higher accuracy of a designed transfer. However, the commonly used plane-restricted four-body model is a non-autonomous system, which challenges much on methods using periodic orbits and their invariant manifolds. In 2005, Shadden et al.  used the Lagrangian Coherent Structure (LCS) to find invariant manifolds of 2D non-autonomous systems. Short et al.  also calculated the LCS under four-body and ephemeris systems. Then they applied it to designing a transfer orbit from the LEO to the vicinity of the Earth-Moon L1 point. Based on that, Qi et al.  proved the existence of a tube-formed time-dependent invariant manifold in the phase space. They also constructed a low-energy transfer trajectory by patching the unstable invariant manifold under Sun-Earth restricted three-body system and the time-dependent stable manifold under the Sun-Earth-Moon system. Same based on the concept of LCS, Onozaki et al.  considered the planar four-body model as two coupled three-body with perturbations and constructed a low-energy Earth-Moon transfer trajectory by patching two time-depend manifolds.
Some studies also focused on the Earth-Moon three-dimensional (3D) low-energy transfer problem [57,110–112]. In the two-dimensional planar CRTBP, the interior and exterior regions of their three-dimensional energy surface are well defined. In the spatial case, however, the invariant manifolds of a halo orbit are two-dimensional while the energy space is five-dimensional. Such a manifold cannot define the interior or exterior topologically well, making three-dimensional transfer significantly complex [57,113]. In the Jovian moon system, Gomez et al.  designed a trajectory ending in a high inclination around the Europa. Howell et al.  obtained the transfers between the Earth-Moon system and Sun-Earth L2 point in three-dimension circumstances. They computed the three-dimensional transit orbit inside the stable Earth-Moon tube to link the exterior and interior regions. The primary solution was then reproduced under the full ephemeris model. In addition, the particular Lissajous trajectories in each system were sought to accomplish the transfer at a low cost. Ren et al. [114,115] analyzed the limitation of the planar manifold assumption and proposed the necessary and sufficient conditions for producing a transit orbit. The statistical character of the natural lunar low-energy transfer orbits and the optimal performance of the artificial lunar transfer orbits were numerically investigated. Similar to the conclusion drawn in  under the planar circumstance, in the spatial case, the L2 transit orbit has more opportunities to realize a lunar collision than the L1 transit orbit. They pointed out that the distribution of the elevation angle of the lunar collision velocity vector relative to the lunar surface can be regular with low energy. And individually increasing the energy level is insufficient in improving the transfer ratio. The success of this patched three-body approach depends greatly on the spatial configuration where the invariant manifolds of the three-body subsystems intersect within a reasonable time . Xu et al.  studied available transfer opportunities under the spatial bi-circular model. They found that in the case of transiting the L2 point, a whole Earth-Moon transfer trajectory can be achieved only within the solar phasic angle interval , while it can be achieved within the pairs of and in the case of transiting the halo orbit near the L2 point.
Although showing a relatively long flight duration, the low-energy transfer has many advantages, including the low cost, an extended launch window, and the capability of fixed epoch arrival . Taking the GRAIL mission as an example, the low-energy transfer provided a long launch period (at least 21 days), which is difficult with the direct lunar transfer, whose consecutive launch dates to target a fixed ascending node usually are lower than three. Moreover, in the low-energy case, any launch date within the launch window leads to a fixed lunar orbit insertion date with essentially one certain lunar orbit design, which is also impossible using the direct lunar transfer . Considering the practical trajectory design, Parker et al.  proposed a targeting scheme to build 21-day launch trajectories from a 185 km Earth parking orbit to the 100 km circular lunar polar orbit, using up to two midcourse maneuvers. The dynamical model used was JPL’s DE421 planetary Ephemerides. They found that the average cost of 2.5 m/s per day could add a launch period, and the cost of the 288 reference transfers studied was approximately 71.7 29.7 m/s (1σ).
Section 2 introduced the traditional direct transfer method based on the two-body Keplerian ellipse from the Earth to the Moon. Such transfers have a hyperbolic excess velocity V∞ relative to the Moon which determines another need to be given to render the spacecraft captured into an elliptic orbit about the Moon.
For the traditional direct transfer orbit, however, the spacecraft cannot reach the moon if the TLI is eliminated. The only available way is to eliminate the spacecraft velocity in the vicinity of lunar periapsis which means the TLI orbits energy should be similar to the Moon’s orbit. Considering the maneuver at the apogee of the transfer orbit is not provided by the main engine but by the perturbation of the Sun or Moon, theoretically, the total required can be further reduced. To achieve this, the spacecraft should be sent to the region where the gravitational interactions between Earth, Sun, and Moon tend to balance which is called weak stability boundary (WSB). In this region, a small maneuver can change the spacecraft motion significantly and return the spacecraft for ballistic captured by the Moon. This type of transfer was firstly constructed by Belbruno and utilized by the Japan Hiten probe . By using this method, only one impulse can send the spacecraft to the Moon directly, while two impulses are needed for Hohmann and three for bielliptic. And the energy-saving over Hohmann and bielliptic are 18% and 37% respectively but a longer transfer time is required .
Various definitions have been given to the WSB. The most basic definition according to Belbruno  is that WSB is a transition region in the phase space where gravitational interactions between Earth, Sun, and Moon tend to balance. If projected into a three-dimension space, WSB describes a continuous three-dimension space and can be seen as a generalized Lagrangian point or the concept of the sphere of influence. Another definition also given by Belbruno  regards WSB as a region in phase space supporting a special type of chaotic motion for a special choice of elliptic initial conditions with respect to m2. Or as “in the CRTBP, WSB is a boundary set in the phase space between stable and unstable motion relative to the second primary. Perturbed by the first primary, Keplerian orbits about the second primary could be stable only if they preserve the character of bound motion after a certain number of revolutions. Yagasaki  described WSB as “a transition region between the gravitational capture and escape from the Moon in the phase space”. In this part, we first give three classical definitions of the WSB, i.e., analytical definition, numerical algorithmic definition, and definition given by Garcia and Gomez.
Definition 1: Analytical definition
Reference  gives the mathematical definition of WSB in the frame of CRTBP: is captured at the time t = t1, if the solution of CPTBP, satisfies
Meanwhile, , where and are position and velocity with respect to in centered inertial coordinate, and the two-body Kepler energy of can be given by
We consider the CRTBP and determine the set where . In addition, those points are considered where , i.e., local periapsis or apoapsis point. Set
Then defines a special set where ballistic capture occurs in the restricted problem. W represents the weak stability boundary. The motion of P3 near W is sensitive.
Definition 2: Numerical algorithmic definition
Reference  expressed WSB by numerical method. Considering a ray from the Moon with perilune distance lying on the , the spacecraft performs a stable motion as it follows a full cycle around the Moon without going around the Earth; otherwise, the motion is unstable. When gradually increasing the initial perilune distance along the , a finite distance can be found. The eccentricity e is another vital parameter for . Define
where has two components. One corresponds to retrograde motion about the Moon and the other to direct motion about after propagation from .
Definition 3: Definition given by Garcia and Gomez
Reference  gave the initial conditions along the that must be integrated to determine the stable or unstable regions around .
Initial conditions with positive velocity (osculating retrograde motion about )
Initial conditions with negative velocity (osculating direct motion about )
In both cases, is the modulus of the initial velocity which is given by
When fixing the value of eccentricity and the angle , all the possible values of the distance along are searched for which there is a change in the stability property of the orbit. For a finite number of points (up to a certain precision) such that if then the motion is stable and otherwise unstable. The number of and their values depend on and computational precision. Then, WSB can be defined as
According to Dutt , there are two common ways to design a ballistic capture transfer trajectory based on WSB theory, i.e., the back-propagation method and forward targeting method . Belbruno et al.  first proposed a methodology for producing a WSB-based low-energy trajectory. The dynamical model used is a general restricted four-body problem as given by
where represents spacecraft, the Sun, the Earth, and the Moon, respectively. And denotes the inertial position of the mass point, . The initial value for can be obtained from an epoch . For convenience, consider a spacecraft’s motion relative to the Earth and introduce a transformation. Let denote the relative position of the mass point in the Earth-centered non-inertial coordinate. The motion of spacecraft is renewed as
Specifically, once the initial value in the WSB region of the Moon is determined, the equation of motion (17) can be integrated back to the trajectory arriving in the WSB3 at approximately 1.5 × 106 km from the Earth. Then for a fixed position relative to the Earth, varying the initial velocity to make the spacecraft leave the Earth and arrive at the target point obtained by previous backward integration. An energy augmentation can be achieved through a lunar gravity assist. The combination of the two trajectory segments is the aimed WSB transfer. Moreover, a differential correction is always needed for targeting, and the equations of motion are solved via a fourth-order Runge-Kutta numerical integration. Yamakawa et al. [118,119] provided a systematic method for trajectory construction and obtained a wide variety of ballistic lunar capture trajectories for different cases. These trajectories were classified into two types: (1) an Earth side approach to the Moon through a geocentric orbit with an initially small semi-major axis and (2) an anti-Earth side approach through a geocentric orbit with an initially large semi-major axis. An analysis of a geocentric orbit considering the solar effect was performed based on simplified Lagrangian’s planetary equation. The investigation of the solar gravity’s influence on the geocentric orbit from the angular momentum viewpoint shows that the perigee distance would increase when the spacecraft location in the second or fourth quadrant in the Sun-Earth rotating frame. A detailed illustration of the solar perturbation effect is shown in Fig.14. Besides, by introducing lunar swing-by, the total flight time can decrease due to the rise of the geocentric orbit’s initial perigee caused by the reduction of the local eccentricity.
Circi et al.  studied the configuration requirement of the WSB transfers. The possibility of the transfers is determined by two important angles: one is the angle between Earth-to-Sun lines and Earth-to-apogee and another is the Sun-Earth-Moon configuration at spacecraft arrival to the Moon. The lunar gravitation capture and solar gravitational effect were analyzed in terms of variations of the Jacobi constant in the Earth-Moon and Earth-Sun systems. They also studied the sensitivity of the lunar encounter region to capture through the parametric analysis. Specific angle values were revealed, under which the lunar capture can be stable. The condition of generating WSB transfers was estimated in the case of the “quasi-ballistically” Moon captured transfer orbits. And the classification substance of the WSB trajectories was explained.
Despite efforts on Earth-Moon transfer orbits design using the WSB theory, the nature of the WSB and associated dynamics has not been well understood until the work done by Garcia and Gomez in 2007 . They proposed a more precise definition of the nth WSB through numerical explorations. The generalized WSB was regarded as a union of the set for n = 1, 2, 3…. A comparison of the stable/unstable regions around the small primary to the stable/unstable invariant manifolds of the periodic orbit about the collinear Lagrangian points was conducted. In the framework of the planar CPTBP, they concluded that for some range of energies, the WSB is located in the closure of the stable manifolds’ union around two of the collinear Lagrangian points. Romagnoli et al.  further investigated the lunar WSB in the planar and spatial cases using both the three-and four-body models. The computation of WSB in the study can be divided into single and multiple stable/unstable transitions. The eccentricity of the lunar orbit and the presence of the Sun showed ignorance influence on the WSB geometry. Moreover, a comparison was made between the lunar polar capture and equatorial capture when a stable and safe orbit (no impact on the lunar surface after capture) is required. Results showed the superiority of the equatorial capture in the capture number, although a higher collision risk is found. In addition, the equatorial capture with low perilune altitude presented the ability to achieve the minimum energy transfer. Ceccaroni et al.  then proposed a new definition of the WSB consistent with the existing algorithmic definitions. Belbruno et al.  also illustrated that for a certain condition of mass ratio and energy level in planar CRTBP, the WSB near the heavier body coincided with a branch of stable manifolds associated to the Lyapunov orbit at one of the collinear Lagrangian points.
Topputo et al.  modified the numerical approach proposed in . They defined the motion equations under the polar coordinates to compute the stable sets, and calculate the boundaries through a bisection algorithm. The modified approach was also employed to produce the WSB of the Sun-Jupiter system. The geometry of the stable sets was compared to that in the Earth-Moon system [117,121] in terms of the mass parameter. The study also discussed the possible correlation between the weak stability boundary and the invariant manifolds associated to the Lyapunov orbits. Belbruno et al.  also excavated the relationship between weak stability boundary and invariant manifolds. A rigorous definition in the context of the planar CRTBP was given, and a geometric argument based on the separatrix property of the invariant manifolds associated to the Lyapunov orbits was provided. They summarized the following points from numerical experiments’ results. First, with the algorithmic definition of the WSB, a prior knowledge of the Lyapunov orbit is no longer needed in finding trajectories on the stable manifolds of the Lyapunov orbit. Second, they indicated that WSB contained the symmetric homoclinic points in the whole range of Jacobi constants where the WSB is defined, which is different from the results in . Moreover, WSB could also exist in models where the hyperbolic invariant manifolds are inapplicable, such as the elliptic restricted three-body problem and the bicircular restricted four-body problem. Hence, the WSB may be a good substitute for the hyperbolic invariant manifolds in such cases. Besides the expected transitions, Sousa-Silva et al.  also found spurious transitions in various types, which was attributed to the separatrix role of the invariant central manifold in the planar CRTBP. The study discussed the spatial distribution of the authentic and spurious transitions along the boundary for sets of initial conditions with high eccentricity. The results showed the regular appearance of spurious transitions and collisional trajectories. In addition, they proposed a new definition of the stability boundary through the investigation of the smooth and fractal-like portions of the boundary. In order to validate the stable solutions in Earth-Moon ballistic capture transfers, Sousa-Silva et al.  also investigated the applicability of the initial condition sets generated by the algorithmic definition of WSB within the lunar SOI in the Earth-Moon low-energy capture transfers’ construction. They found subsets of stable conditions which provide feasible Earth-Moon transfers. In the study, an external low-energy Earth-Moon transfer with zero midcourse correction at the patching section was constructed. Griesemer  proposed a numerical targeting scheme to design ballistic lunar capture transfers. The initial guesses obtained through a simpler dynamical system were accurate enough to converge to the solution in the numerical targeting scheme. In the algorithm, a particular member of a family of periodic orbits, documented by Makellos  as family f 16, was utilized as the initial guess for an Earth-Moon transfer. van der Weg et al.  analyzed the lunar orbits and some of the possible contingencies for the European Student Moon Orbiter (ESMO). The operation orbit selected was a highly eccentric one, and the 6 months WSB transfer was obtained. A parameter study was undertaken for the sensitivity investigation of the lunar orbit insertion. And a database of transfer solutions from 2014 to 2015 was used to evaluate the robustness of weak capture and the planetary geometry at lunar arrival. Dutt et al.  also constructed the WSB transfer via back-propagation in CRTBP. Concretely, two highly elliptical geocentric orbits about the Earth and the Moon were propagated forward and backward respectively and patched using the Fixed Time of Arrival targeting method. The genetic algorithm was used to find patching points to reduce the required number of patching points. Various selenocentric and geocentric orbits were considered, and eligible candidates for WSB transfers were plotted on the phase space. Moreover, they investigated the fly-by Moon trajectory to show the reduction of the total flight time using this approach.
Additionally, Belbruno  proposed a forward targeting algorithm, making the generation of WSB transfers more accessible. In his study, the software package STK/Astrogator was used, and an optimization technique was employed for improving launch opportunities and saving. Dutt et al.  developed a method to find Earth-Moon WSB transfers by forward propagation. The problem was considered as an optimization of the two loops. The outer loop varied the departure parameters from an Earth parking orbit while the inner loop obtained a capture trajectory to the Moon by varying the small impulse given at the WSB of the Earth. Once a lunar capture orbit was obtained, they were able to achieve the minimum eccentricity for circularization.
The low-thrust method typically uses low-thrust propulsion systems such as electric propellants or solar sails. Employing high specific impulse, the fuel efficiency of low-thrust transfers enables the enhancement of payload capability and the potential of achieving high fuel consumption missions. Different from ballistic trajectories mentioned before, a spacecraft that adopted low-thrust transfer would be propelled for long period or even persistent using the low-thrust engine. Therefore, it experiences the gravitational attraction of the primaries and other perturbations along with the change in energy due to thrusters. There are two key problems in the low-thrust method: the dynamic equation under the action of continuous low-thrust has no analytical solution, and the thrust direction, thrust magnitude, and timing of thrust application during a flight are difficult to determine; The low-thrust propulsion system is too powerless to operate instantaneous braking as the high-thrust propulsion system does, hence it is hard for the spacecraft to be captured when it reaches the vicinity of the Moon. The design of low-thrust transfer trajectories can be modeled as a variational problem of optimizing fuel consumption or transfer time, and the method to solve this problem can be divided into the indirect method, direct method, hybrid method, and evolutionary algorithms [134,135]. A graphical illustration of approaches used to solve the spacecraft low-thrust transfer trajectory optimization problems is displayed in Fig. 15.
The indirect method is commonly used in the early study of transfer trajectory optimization. Based on the constraint equations and optimization index of an original problem, the indirect method introduces Lagrangian multipliers to construct the co-state equation and new optimization index to transform the low-thrust transfer trajectory optimization design problem into a two-point boundary value problem (TPBVP) with nonlinear ordinary differential equations. Both the dynamic and boundary constraints include the first-order optimal necessary conditions derived from the calculus of variations and Pontryagin’s Maximum Principle (PMP) .
However, deriving the first-order optimal necessary conditions is challenging in this case . For some special low-thrust orbits, Lawden  described the first-order optimal necessary conditions based on the primer vector theory. Russell  simplified the low-thrust transfer trajectory problem to only a few parameters through the primer vector theory and calculus of variations. Grigoriev et al. [140,141] derived the first-order optimal necessity conditions for complex dynamic models and complex interior point constraints. They also introduced the complementary relaxation and rigidity conditions into inequality constraints.
Moreover, an initial guess is required for the co-state variable. However, in the indirect method, the convergence region is narrow and the results are sensitive to initial guesses [142,143]. Many auxiliary methods have been proposed to release this sensitivity, such as the Adjoint Control Transformation (ACT) [144–147], normalization techniques [148,149], and etc. [142,150,151]. The homotropy method is also emerging in the optimization design of low-thrust Earth-Moon transfer trajectories [152–154]. It was first introduced by Bertrand et al.  in 2002. This method can effectively overcome the discontinuous control problem in the fuel optimization problem and can increase the convergence area of the optimization algorithm.
Efforts are also made in the optimal design of the low-thrust transfer trajectory in the cislunar space based on the indirect method. The early studies always divided the whole trajectory into the Earth phase, the middle phase, and the Moon phase, respectively. Then they solved the optimal control problem of each phase separately and patched the trajectories together. Based on this approach, Golan et al.  studied the Earth to polar Moon orbit low-thrust transfer problem in both two-and three-dimensional cases. Ranieri et al.  used a homotopy process by multiplying a scale factor with the terms due to third-body perturbation and rotational effect and obtained the spatial trajectory in the Earth-Moon rotating frame. To reduce the computational burden, Lee et al.  proposed an initial co-state estimation method to construct spiral trajectories directly for an arbitrary initial radius, transfer time, and target-specific energy in the Earth-Moon rotating frame.
Recently, the study on designing the Earth-Moon low-thrust transfer trajectories based on the invariant manifold has been carried out to make full use of the characteristics of three-body dynamics. Senent et al.  adopted an indirect method combined with the adjoint control transformation (ACT) to transfer a power-limited spacecraft from an LEO to the stable manifold associated to an unstable periodic orbit in the CRTBP. The guidance algorithm developed was also proven to be capable of dealing with both periodic unmodeled gravitational perturbations and power fluctuations. Oshima et al.  studied the global search strategy for low-thrust transfer to the Moon in the planar CRTBP. By implementing the necessary condition of optimality based on the Pontryagin extremum principle, a wide range of Pareto solutions were obtained, many of which exploited high-altitude lunar fly-by to reduce fuel consumption. Pérez-Palau et al.  studied the design of minimum fuel low-thrust trajectories from an LEO to a lunar orbit with the BCM model. In their study, two methods were proposed to overcome the numerical difficulties that arose from the huge sensitivity of the state and costate equations. The trajectories found were classified into several families based on orbital shape, transfer duration, and fuel consumption, and a physical interpretation was also given to the different families of trajectories based on invariant manifold theory. Taheri et al.  used the indirect method along with Lawden primer vector theory to design low-to high-thrust trajectories from a GTO to a halo orbit of the Earth-Moon L1 point. Considering the challenge caused by the sensitivities of the dynamical models, they introduced the continuous State transition matrix (STM)  to calculate the constraint sensitivities. This approach enables them to solve fewer subproblems by changing the continuation parameter at a greater rate when using the homotopy method. Same targeting the extreme sensitivity of the co-state variable when using the indirect method, Hecht et al.  combined the heuristic optimization algorithm and the homotopy continuation technique to obtain the effective initial value of the co-state variable. After determining a guess of the co-state variables, a single shooting method was employed to converge to the fuel-optimal solution. The method was then successfully used in a GTO-L1 halo orbit transfer problem. By combining the halo orbit manifolds at the Earth-Moon L1 point, Singh et al.  studied the end-to-end low-thrust transfer problem from a GTO orbit to an LLO orbit by exploiting the manifolds associated to the Earth-Moon L1 halo orbit. The homotopy method was still adopted for the final boundary conditions to improve the convergence. Moreover, considering the effect of Sun’s perturbation, they made a comparison between BCM and CRTBP models and achieve some conclusions.
With the development of computer technology and the enhancement of numerical computing capabilities, the direct method has become a widely used method in the optimization design of low-thrust transfer trajectory [166,167]. The direct method  discretizes the variational problem and transforms it into a nonlinear optimization problem with multiple nodes and optimization variables. Then the low-thrust transfer trajectory can be obtained by minimizing the performance index by using the nonlinear programming (NLP) algorithm. Different from the indirect method, the direct method does not have high requirements on the initial value. And the decision variables in NLP are directly related to the physical meaning, so its initial value is easier to obtain . However, the computation burden can rapidly increase when dealing with long-term or many-revolution orbital transfers because a large number of nodes are required for numerical iteration .
Betts et al.  used the transcription method to compute an optimal low-thrust transfer trajectory from an Earth orbit to a specified lunar mission orbit. The SMART-1-like trajectory as shown in Fig. 4 was modeled with three distinct phases as Pierson et al.  did. The vehicle dynamics were defined using a set of modified equinoctial coordinates. The solution presented required the solution of a sparse optimization problem with 211031 variables and 146285 constraints. The trajectories presented required two very long thrust arcs which make the overall mission duration shorter than multi-burn trajectories. Herman et al.  used the direct method to solve the problem of the shortest flight time of the low-thrust Earth-Moon transfer orbit with the initial acceleration of the spacecraft greater than 1 × 10−3 m/s2. Because of the relatively large initial acceleration, the boundary constraints of the initial orbit were not considered in his study.
Combining low energy and low-thrust approach, Mingotti et al.  constructed a WSB-like transfer trajectory shown in Fig.16. They assumed that an initial impulse launched the spacecraft into the stable manifold of the Sun-Earth Lyapunov orbit. Then an attainable set could be defined in the Poincaré section by integrating the controlled dynamic model backward and the whole design can be simplified to the search for a single point on a suitable surface of a section. The first guess obtained from the coupled three-body approximation was then optimized in the controlled four-body dynamics using a direct multiple shooting strategy. The maximal value of the thrust was 0.5N and the backward integration time was set as 6.90 EM in this study. The study was then extended into a lunar fly-by type . Another study carried out by Mingotti et al.  also studied the low-thrust Earth-Moon transfer problem, however, the backward propagation time was set to free. The optimal result showed that different from reference. , for the global optimal solution, a thrust arc could also be employed when the spacecraft is far from the Earth about six Earth-Moon distances.
The hybrid method  is a mixture of direct and indirect methods. The method uses a co-state variable differential equation and control law in the first-order necessary conditions, but the transversal condition is not adopted. In this frame, the optimization variables are more than constraint conditions, making the optimization problem a non-complete boundary value problem, and generally solved through nonlinear programming. Many early studies are carried out based on the thrust-coast-thrust control strategy. Pierson et al.  and Kluever et al.  adopted the three-stage approach to decompose the optimal problem from a circular LEO to a circular LMO. The designing variables were parameterized from the indirect method, and sequential quadratic programming (SQP) was employed to solve the TPBVP problem. The discussion in  was focused on planar trajectories, while reference  also considered three-dimensional trajectories with a Moon orbit inclined 90° to the Earth-Moon plane. Besides the three-stage approach, Kluever et al.  also proposed a modified method to develop the minimum-fuel, three-dimensional trajectory from a circular LEO to an inclined circular LMO. Ozimek et al.  used both direct single-and multiple-shooting methods to carry out low-thrust transfer trajectory designs from the Earth parking orbit to the Lagrangian point periodic orbits of the Earth-Moon system (including L1 halo orbits, L1 and L2 vertical orbits, and L2 butterfly orbits). The manifolds are utilized in the coast phase of the transfer trajectory. Qi et al.  established a low-thrust trajectory optimization problem from both 500 km LEO and 3,629 km medium Earth orbit (MEO) to the L1 halo orbit. By combining the variations approach and direct parameter optimization, the direct and the indirect transfer trajectory considering lunar flybys are obtained, both of which include manifold coasting phases, as shown in Fig. 17. In general, the discretization scale of the hybrid method is much smaller than that of the direct method, and the hybrid method shows higher precision, wider convergence region, and stronger numerical stability than the indirect method. But like the direct method, the quality of the result depends on the size of the discrete and the strength of the optimization algorithm. Moreover, it is usually necessary to assume the control strategy in advance and derive complex first-order necessary conditions.
Evolutionary algorithms optimize discrete variables by simulating the evolutionary process of nature such as the genetic algorithm, ant colony algorithm, and particle swarm algorithm. This method is essentially different from the above three numerical methods: (i) a large number of initial populations are randomly generated as initial guesses during the solution process, with no need for specific initial guesses; (ii) the gradient information of the optimization index is not necessary, and the optimal global solution in the search space relies on the evolution process to iterate. Evolutionary algorithms are mainly used to solve trajectory optimization problems with fewer discrete variables and determine discrete parameters such as mission start/end times, gravitational assistance times, and sequences . It can also be combined with processes like the shape-based method for determining relevant parameters [181,182] and the indirect method for selecting a more reasonable initial guess [149,183].
Besides the optimization methods introduced before, a shape-based method was also presented for initial estimates of low-thrust, gravity-assist trajectories. For the shape-based method, the low-thrust trajectory is parameterized with the basis function. The constraint can be dynamics equations and boundary conditions. The low-thrust trajectory can be obtained by solving the unknown parameters in the transit orbit . Taheri et al.  used the finite Flourier series (FFS) to construct the low-thrust trajectories with thrust acceleration constraints in the CRTBP. In their work, the ballistic intermediate arc of a typical low-thrust problem was replaced with a Flourier series approximation. By applying various numbers of revolutions around the primaries, the method demonstrated capability in establishing both prograde and retrograde final orbits. In addition, the method can approximate various shapes of the intermediate phase that does not follow a spiral trajectory. Madhusudan  developed a robust Flourier series approach for the fast generation of Earth-Moon trajectories using low thrust. Each component of the position vector was approximated using FFS as a function of time. Approximations were then used to design a trajectory that satisfies the equation of motion and constraints at discrete points, and the problem boundary conditions. Taking advantage of CRTBP characteristics, the thrust-free phase trajectory near L1 was designed and optimized to achieve plane change. Three-dimensional transfers to high lunar orbits, low lunar orbits, and halo orbits can be generated.
This paper reviewed various Earth-Moon transfer trajectories modeling and design methods. The engineering favored direct transfer based on the patched conic and the pseudo-state model features high energy consumption but short flight time. Since the gravitational force of two celestial bodies is considered simultaneously in the pseudosphere, the pseudo-state presents higher design accuracy than the patched conic method but is still lower than the three-body model. In comparison, the low-energy transfer effectively reduces the fuel requirement but suffers from a significant drawback of relatively long flight duration. Two commonly used mechanisms of designing low-energy transfer, namely the manifold and the weak stability boundary theory, are presented. The low-energy transfer trajectories obtained through the manifold theory can be varied according to the utilized Lagrangian point. As demonstrated by many groups, transiting through the L1 point enables the energy-minimal interior or cislunar transfer to be acquired, while the energy-minimal exterior or translunar transfer can be obtained by transiting the L2 point, transiting a periodic orbit can increase the chance to realize the Earth-Moon transfer but suffer a high energy consumption. In addition, the chaotic control method of interior transfer basically relies on the exponential amplification of small disturbances in the chaotic motion to reduce the flight time in the chaotic region. Multiple impulses are required in the process, and a series of periodic orbits have to be passed through. A low-energy Earth-Moon interior transfer can be more easily constructed based on the transit orbits surrounded by invariant manifolds, but the method will limit the transfer mode, flight time, and energy consumption of the spacecraft. Since it’s impossible to reach the Earth’s vicinity for the invariant manifold associated to the periodic orbit of the Lagrangian point, neither of the above two methods can directly construct a transfer orbit from LEO to the Moon. Therefore, additional maneuvers are needed for the spacecraft in the LEO to enter the manifold tube. In comparison, with the employment of an exterior transfer, a spacecraft can reach the Moon from the translunar space via the Sun-Earth non-transit orbit as the Sun-Earth manifold can reach the vicinity of the Earth. This method does not require additional connecting arcs, but the results are highly reliant on the geometry configuration of the Sun-Earth-Moon-spacecraft system. Moreover, the solar gravity can cause the distortion of the manifold tube, and the concept of a time-dependent manifold should be introduced to refine the model. The weak stability boundary theory is another important method to construct the low-energy exterior transfer, which has been successfully applied in practical engineering tasks by combining with the planar restricted three-body or four-body problem. However, there is a lack of theoretical support in the spatial case. The direct transcription, which maps an optimal control problem into a nonlinear programming problem, illustrates robustness in solving the optimal trajectory in the multi-body system. More solutions can be obtained when using a global optimization method at the expense of computation time. A common approach is applying the local optimization method for the optimal solution search from a rough initial value. However, this method is not the main topic of this review. The low-thrust transfer method demonstrates efficiency in saving fuel consumption with the conjunction of the manifold theory. Although limited by the difficulties of parameters optimization and lunar capture, the method shows promise in future Earth-Moon exploration missions. It is worth noting that the low-thrust method cannot directly design low-energy external transfer, and an initial impulse is still needed.
Based on the detailed review of Earth-Moon transfer trajectory design methods, a few urgent research topics are identified.
Full lunar surface coverage, high efficiency, and reusability are the three main requirements for the future lunar exploration. Using the invariant manifold of the three-body system can achieve ballistic landing from the Lagrangian point periodic orbit, but the lunar surface coverage is limited using the natural manifold. Also, a contradiction between the transfer time and fuel consumption always exists in the present design methods. Therefore, more efforts are needed in designing the landing orbit to achieve global arrival with minimum fuel consumption. In addition, a flight-based transportation network and an emergency return track are required to realize the transportation of goods and personnel.
In terms of the low-thrust orbit designing method, there are still some unsolved astrodynamics problems for constructing an Earth-Moon transfer trajectory. Although considerable progress has been made in the optimization design of low-thrust transfer based on the two-body and two-body perturbation models under the inertial frame, the method has been barely employed in the three-body dynamic model due to the unsolved theoretical problems. And the optimization algorithm shows poor convergence and low computational efficiency in this case.
Funding Statement: This work was supported by the National Key Research and Development Program of China (No. 2021YFA0717100) and the National Natural Science Foundation of China (Nos. 12072270 and U2013206).
Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.