Family of 2:1 resonant quasi-periodic distant retrograde orbits in cislunar space

Given the current enthusiasm for lunar exploration, the 2:1 resonant distant retrograde orbit (DRO) in Earth-Moon space is of interest. To gain an in-depth understanding of the complex dynamic environment in cislunar space and provide more options for parking orbits, this paper investigates the existence of quasi-periodic orbits near the 2:1 resonant DRO in the circular restricted three-body problem (CR3BP). Firstly, the numerical computation approach, continuation strategy, and stability analysis method of quasi-periodic orbits are introduced. Then, addressing the primary challenges in the continuation progress, we have developed an adaptive continuation algorithm with automatic adjustment of the step size and the number of discrete points that characterize the invariant torus. Finally, two types of 2D quasi-DROs and their linear stability properties are explored. Using Poincaré sections, we investigated the boundaries of the maximum extent attainable by both 2D quasi-DRO families in the CR3BP at a specific Jacobi energy, confirming that both types of quasi-periodic families have reached their respective boundaries. The algorithm described in this paper is beneficial for facilitating the computation of quasi-periodic families and aids in discovering additional potential dynamical structures.

1 Introduction

In recent years, there has been a resurgence of enthusiasm for lunar exploration, making cislunar space a focal point for human exploration and research. Several space missions, some successfully launched in the past years, and others scheduled for the near future, have been proposed for verification, scientific exploration, or space observation purposes. Notable missions include CAPSTONE (Gardner et al., 2023), KPLO (Song et al., 2021), ARTEMIS 1 (Williams et al., 2023), and VIPER (Colaprete et al., 2019). The exploration of cislunar space presents a prosperous prospect. As a stable orbit family suitable for long-term parking in cislunar space (Bezrouk and Parker, 2014), distant retrograde orbits (DROs) have garnered particular attention from institutions all over the world. The stability and strategic significance of DROs make them a compelling choice. According to existing literature, DROs have four primary applications: serving as a space harbor for scientific exploration, a candidate orbit for space domain awareness, a transit station for deep space exploration, and a relay orbit for inter-satellite links (Smitherman and Griffin, 2014; Stramacchia et al., 2016; Conte et al., 2018; Wang et al., 2019).

While scholars have extensively studied the dynamics and structural characteristics of DROs in the CR3BP, this knowledge falls short for practical engineering applications. For instance, the eclipse situation on DROs is suboptimal since DROs are nearly located on the Earth-Moon plane even when transitioned to the ephemeris model. This directly impacts the on-orbit stability of spacecraft. McCarthy B. P. and Howell (2021) have explored the potential for eclipse avoidance in quasi-DROs with $z$-axis components. Moreover, considering that DROs are quasi-periodic in the ephemeris model, incorporating quasi-DROs into the preliminary orbit design process for subsequent missions would be highly advantageous. This approach provides a more accurate initial guess, as demonstrated in McCarthy and Howell (2019). Additionally, from a dynamic perspective, the study of quasi-periodic motion near DROs contributes to a deeper understanding of the dynamical behaviors in cislunar space. This has sparked our interest in investigating the dynamical structures of quasi-periodic orbits. With numerous missions currently targeting the 2:1 distant retrograde orbit (2:1 DRO), this study primarily focuses on the analysis of quasi-periodic motion in the vicinity of the 2:1 DRO. The prevalence of quasi-periodic orbits and their distinct advantages over periodic orbits contribute to a broader range of orbit options, thereby expanding the design space (Olikara and Scheeres, 2012).

The CR3BP stands out as the simplest dynamic model for studying quasi-periodic motion in a multi-body system. In this model, the solution space is composed of four types of behaviors: equilibrium points, periodic orbits, quasi-periodic orbits, and chaotic motion (Folta et al., 2016). Olikara and Howell (2010) computed the family of Earth-Moon $L_1$ Lissajous tori, while Jorba et al. (2020) explored quasi-periodic motion in the vicinity of $L_1$. For quasi-periodic motion near $L_2$, one can refer to works in McCarthy B. and Howell (2021); Rosales et al. (2021); Lujan and Scheeres (2022). Furthermore, Jorba and Nicolás (2020) provide insights into dynamics around $L_3$. Quasi-periodic motion near equilibrium points has been well studied, as evidenced by the works of Hou and Liu (2010); Hou and Liu (2011). However, to the best of our knowledge, quasi-periodic structures near resonant DROs have not been thoroughly investigated.

The computation method of quasi-periodic orbits has undergone a transition from analytical or semi-analytical methods to full numerical methods, as discussed in Farquhar and Kamel (1973), Gómez et al. (1998), Jorba and Masdemont (1999), Hou et al. (2015). Analytical or semi-analytical methods like the Poincare-Lindstedt method and center manifold reduction often face challenges with small convergence regions. Therefore, we prefer the numerical method in this manuscript, specifically the GMOS (Gomez-Mondelo-Olikara-Scheeres) algorithm (Olikara and Scheeres, 2012), which iteratively computes the invariant curves represented by Fourier series on a stroboscopic map using Newton’s method. For additional details on other numerical methods, one can refer to Baresi et al. (2018) and the references therein. Numerous works in the literature have been dedicated to the calculation of quasi-periodic orbits or quasi-periodic families. Some focus on calculation skills, while others address computational efficiency (Schilder et al., 2005; Sánchez and Net, 2013; Olikara, 2016). Few studies address whether the continuation reaches the boundary, which demarcate the regions of bounded motion and chaos. As they do not clarify whether the failure of continuation is due to encountering resonance or reaching the continuation boundary, there is a high probability that the quasi-periodic families are incomplete. Nevertheless, it is crucial to recognize the boundaries of chaos and order.

Motivated by the above analyses, this paper delves into the quasi-periodic motion and its stability properties near the 2:1 resonant DRO in the CR3BP. However, The continuation process is not without challenges. Specifically, the quasi-periodic families encounter a series of resonant regions during the continuation process, requiring an increased number of sampling nodes to achieve a specified accuracy, as reported in Jorba and Olmedo (2009). If the step size is inappropriate, it can easily lead to singularities, resulting in program failure. Furthermore, the shape and size of the torus are factors determining the number of Fourier nodes, as discussed in Rosales et al. (2021). Considering these factors, we have developed an adaptive continuation framework that adjusts the continuation step size and the number of nodes to save computational time and enhance robustness. Pseudo-arclength continuation and natural parameter continuation schemes are employed for cross-verification. Additionally, with the aid of the Poincaré map and bifurcation theory, we have confirmed that the quasi-DRO families have reached their continuation boundaries.

The remaining structure of this paper is outlined as follows. Section 2 introduces the circular restricted three-body problem and presents a visual representation of the 2:1 resonant Distant Retrograde Orbit. Section 3 provides a detailed overview of the computational framework, discussing the challenges encountered in our calculations and introducing an adaptive algorithm. In Section 4, two types of 2D quasi-periodic orbit families are presented, and the boundaries of the orbit families are also verified. Stability analysis is conducted in Section 5. Lastly, conclusions are drawn in Section 6.

2 Preliminaries

2.1 Circular restricted three-body model

The circular restricted three-body problem (CR3BP) is a widely used model for studying spacecraft motion in multi-body systems. It describes the motion of an infinitesimal mass influenced by two massive bodies, such as the Earth and the Moon in cislunar space, within a rotating coordinate system centered at the Earth-Moon barycenter. The $x$-axis is aligned from Earth to the Moon, the $z$-axis is defined to coincide with the direction of the system’s angular momentum, and the $y$-axis is determined by the right-hand rule, as shown in Figure 1. The distance, mass, and time in the dynamical system are normalized such that

$\left\{\begin{array}{c}\left[M\right]=m_1+m_2\hfill \\ \left[L\right]=a_{12}\hfill \\ \left[T\right]=\sqrt{a_{12}^3/G\left(m_1+m_2\right)}\hfill \end{array}\right.(1)$

In Equation 1, $m_1$ and $m_2$ are, respectively, the mass of the Earth and the Moon. $a_{12}$ represents the Earth-Moon distance. Define $\mu =\frac{m_2}{m_1+m_2}$ as the mass ratio of the Earth-Moon system, then the Earth and the Moon are located at $(-\mu ,\mathrm{0,0})$ and $(1-\mu ,\mathrm{0,0})$ in the synodic frame, respectively. The dynamical equations that depict the motion of spacecraft in the Earth-Moon rotating frame are described as (Topputo, 2013).

$\begin{array}{cc}\hfill \ddot{x}-2\dot{y}& =\frac{\partial {\mathrm{\Omega }}_3}{\partial x}\hfill \\ \hfill \ddot{y}+2\dot{x}& =\frac{\partial {\mathrm{\Omega }}_3}{\partial y}\hfill \\ \hfill \ddot{z}& =\frac{\partial {\mathrm{\Omega }}_3}{\partial z}\hfill \end{array}(2)$

Figure 1

Figure 1. Diagram of the circular restricted three-body problem. The red curve represents the 2:1 DRO.

In Equation 2, the effective potential ${\mathrm{\Omega }}_3$ is defined as

${\mathrm{\Omega }}_3\left(x,y,z\right)=\frac{1}{2}\left(x^2+y^2\right)+\frac{1-\mu }{r_{ES}}+\frac{\mu }{r_{MS}}+\frac{1}{2}\mu \left(1-\mu \right)(3)$

In Equation 3, $r_{ES}{r}_{MS}$ respectively represent the distances from the spacecraft to the Earth and the Moon. Table 1 shows the relative value of physical constants in the Earth-Moon systems as well as their associated units and physical meanings.

Table 1

Table 1. Relevant physical constants in the Earth-Moon system.

The CR3BP is an autonomous system with a unique energy integral, i.e., the Jacobi constant as displayed in Equation 4,

$C=-\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right)+2{\mathrm{\Omega }}_3\left(x,y,z\right)(4)$

2.2 2:1 DRO in the CR3BP

The DROs are a special family of periodic orbits in the CR3BP, which moves clockwise (retrograde) around the Moon in cislunar space. Due to the symmetry of the DRO family about the $xz$ plane in the Earth-Moon rotating frame, the two intersections of one particular DRO satisfy the vertically crossing condition. By choosing the plane $\mathrm{\Sigma }:y=0$ as the Poincaré map, the particular DRO can be represented by a two-dimensional state, denoted as $\left(x_0,{\dot{y}}_0\right)$. Figure 2 shows the DRO family in the Earth-Moon system, obtained through differential correction and numerical continuation. The 2:1 resonant DRO is highlighted in orange, with an approximate period of 13.65 days. The differential correction equation can be formulated as

$\left[\begin{array}{c}\hfill \delta y_t\hfill \\ \hfill \delta {\dot{x}}_t\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill {\mathrm{\Phi }}_{\mathrm{2,1}}\hfill & \hfill {\mathrm{\Phi }}_{\mathrm{2,5}}\hfill & \hfill {\dot{y}}_t\hfill \\ \hfill {\mathrm{\Phi }}_{\mathrm{2,4}}\hfill & \hfill {\mathrm{\Phi }}_{\mathrm{2,5}}\hfill & \hfill {\ddot{x}}_t\hfill \end{array}\right]\left[\begin{array}{c}\hfill \delta x_0\hfill \\ \hfill \delta {\dot{y}}_0\hfill \\ \hfill \delta t\hfill \end{array}\right](5)$

where $\mathrm{\Phi }$ represents the 6 × 6 state transition matrix. The subscripts in Equation 5 specify the particular row and column elements of $\mathrm{\Phi }$ respectively.

Figure 2

Figure 2. DRO family in the Earth-Moon system. The orange curve is the 2:1 resonant DRO.

In general, the resonant ratio $p:q$ means that the spacecraft has made $p$ revolutions along the orbit in a total $q$ synodic periods, where the synodic period is about 27.3 days in the CR3BP.

3 Computation method of quasi-periodic orbits

A typical approach to exploring the phase space of a dynamical system involves the investigation of its invariant sets. The term “invariant sets” encompasses equilibrium points, periodic orbits, and invariant tori within the system (Rosales et al., 2021). These invariant sets, along with their associated invariant manifolds, provide valuable insights into the evolution of dynamical systems within phase space. Generally, quasi-periodic motion occurs on two- or higher-dimensional invariant torus surfaces. Equilibrium points and periodic orbits can be considered as special instances of invariant tori—specifically, zero- and one-dimensional invariant tori within phase space (Olikara and Scheeres, 2012).

3.1 Initial guess of invariant tori

For Hamiltonian systems, the linear approximation of the quasi-periodic invariant torus is often related to the center manifolds of the periodic orbits (Lujan and Scheeres, 2022). Considering periodic orbits as one-dimensional invariant tori, there exists a differential isomorphic mapping $\bar{u}\left({\theta }_1\right)=\left[\begin{array}{c}\hfill r\left({\theta }_1\right)\hfill \\ \hfill v\left({\theta }_1\right)\hfill \end{array}\right]:{\mathbb{T}}^1\to {\mathbb{R}}^n$ allowing the periodic orbits to be mapped from the state space with six degrees of freedom into a single-parameter phase angle space. In a Hamiltonian system, the monodromy matrix is symplectic, and its associated eigenvalues appear as reciprocal pairs. The eigenvalues of the monodromy matrix provide insights into the manifold structure near the periodic orbit. A conjugate complex eigenvalue pair with a modulus of 1 corresponds to the center manifold of the periodic orbit. Let ${\mathrm{\Phi }}_T\left(\bar{u}\left({\theta }_1\right)\right)$ denote the state transition matrix of a periodic orbit integrated from initial point $\bar{u}\left({\theta }_1\right)$. Representing the complex eigenvalue as ${\lambda }_{j+1}=e^{i{\alpha }_{j+1}}$ and the associated complex eigenvector as $y\left({\alpha }_{j+1}\right)$, the initial guess of the invariant torus can be expressed as (McCarthy B. P. and Howell, 2021):

$u\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)=\bar{u}\left({\theta }_1\right)+\kappa \mathrm{Re}\left[{\displaystyle \sum _{j=1}^m}{e}^{i{\theta }_{j+1}}y\left({\alpha }_{j+1}\right)\right](6)$

In Equation 6, $\kappa$ is a small quantity. The initial guess of rotation angle and frequency associated with each dimension in the phase space can be given as

$\begin{array}{cc}\hfill & \rho =\left[2\pi ,{\alpha }_2,\dots ,{\alpha }_{m+1}\right]\hfill \\ \hfill & \omega =\left[\frac{2\pi }{T},\frac{{\alpha }_2}{T},\dots ,\frac{{\alpha }_{m+1}}{T}\right]\hfill \end{array}(7)$

It should be noted that all frequencies are nonresonant in Equation 7, i.e., for any $\mathit{k}\in {\mathbb{Z}}^n\backslash \{0\}$, $\mathit{k}\cdot \mathit{\omega }\ne 0$. Thus, the family of quasi-periodic orbits is Cantorian. A quasi-periodic orbit densely covers the surface of the torus within a long evolution time. If resonance occurs between those frequencies, the torus will degenerate into periodic orbits or lower-dimensional tori. If the computation continues, the torus typically deviates from the original orbit family and becomes distorted.

3.2 Invariant constriant

Quasi-periodic motion occurs on an invariant torus in a specific manner. As mentioned earlier, for any point initially located on the surface of the torus, when it is integrated for a certain time, the final state is still on the torus but rotated by a fixed phase compared to the initial state. This property forms the basis for calculating invariant tori or quasi-periodic orbits. Since all phase components do not resonate with each other, the quasi-periodic orbit gradually fills the entire invariant torus as time evolves. Therefore, computing a quasi-periodic orbit is equivalent to computing an invariant torus. The computation algorithm used in this paper is based on the GMOS algorithm proposed by Olikara and Scheeres, as described in Olikara and Scheeres (2012). Pseudo-arclength continuation and natural parameter continuation are both explored for verification purposes.

For the dynamical system of interest, it can typically be expressed as an ordinary differential equation in state space form, as follows (Jorba and Olmedo, 2009):

$\dot{x}=f\left(t,x,\lambda \right)(8)$

where $x\in {\mathbb{R}}^n$, and $\lambda \in {\mathbb{R}}^q$ is the external disturbance parameter. For the CR3BP, the system dynamics is time independent, and the external disturbance parameter $\lambda$ is 0. Now, the paradigm of computing invariant torus can be described as follows. Suppose that the invariant torus is characterized by $m+1$ non-resonant periodic components, where the first periodic component is typically fixed to 0 so as to reduce the dimension of the invariant torus. There exists a map $\mathrm{P}:{\mathbb{T}}^{m+1}\to {\mathbb{R}}^6$ such that for any point on the torus, $u\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)$ is a differential isomorphic of the quasi-periodic motion. Suppose that ${\varphi }_T\left(x_0\right)$ is the flow function generated by Equation 8, then if the initial state $u\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)$ is propagated for a stroboscopic mapping time $T=\frac{2\pi }{{\omega }_1}$, it gives

$\begin{array}{cc}\hfill {\varphi }_T\left(u\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)\right)& =u\left(u\left({\theta }_1+2\pi ,{\theta }_2+{\omega }_{2}T,\dots ,{\theta }_{m+1}+{\omega }_{m+1}T\right)\right)\hfill \\ \hfill & =u\left({\theta }_1,{\widehat{\theta }}^T+{\widehat{\rho }}^T\right)\hfill \end{array}(9)$

where $\widehat{\theta }={\left[{\theta }_2,{\theta }_3,\dots ,{\theta }_{m+1}\right]}^T$, and $\widehat{\rho }=\frac{2\pi }{{\omega }_1}{\left[{\omega }_2,{\omega }_3,\dots ,{\omega }_{m+1}\right]}^T$ is the rotation vector. One can conclude from Equation 9 that there is an invariant map with respect to the rotation vector $\rho$. As such, if a rotation operator $R_{-\rho }$ is defined to eliminate the rotation along the directions of $\rho$, the invariant constraint can be formulated as

$R_{-\rho }\left[{\varphi }_T\left(u\left({\left[\cdot ,{\widehat{\theta }}^T\right]}^T\right)\right)\right]-u\left({\left[\cdot ,{\widehat{\theta }}^T\right]}^T\right)=0(10)$

therefore, from a numerical point of view, computation of quasi-periodic torus is equivalent to seeking the solution of Equation 10. A practical numerical method is to represent the high-dimensional reduced torus on the stroboscopic map in the form of (truncated) Fourier series,

$u\left(\theta \right)={\displaystyle \sum _{k_1=0}^{N_1-1}}\cdots {\displaystyle \sum _{k_m=0}^{N_m-1}}{C}_{\vec{k}}{e}^{i\left(k_1{\theta }_2+k_2{\theta }_3+\cdots k_m{\theta }_{m+1}\right)}(11)$

where $C_{\vec{k}}\in {\mathbb{C}}^6$ is the Fourier coefficient. $\vec{k}={\left[k_1,k_2,\dots ,k_m\right]}^{\text{T}}$ is the harmonic set of the high-dimensional Fourier series. $N_j$ is the number of discrete nodes for each phase component ${\theta }_{j+1}$, $(j=\mathrm{1,2},\dots ,m)$. To this end, ${\theta }_{j+1}$ is uniformly discretized into elements in the set shown as follows:

${\theta }_{j+1}\in \left[\begin{array}{cccc}\hfill 0\hfill & \hfill \frac{2\pi }{N_j}\hfill & \hfill \cdots \hfill & \hfill \frac{2\pi \left(N_j-1\right)}{N_j}\hfill \end{array}\right](12)$

The first phase angle component is typically fixed (usually set to 0), then for a ($m+$1)-dimensional invariant torus, it is parameterized into ${\prod }_{j=1}^m{N}_j$ nodes. The Fourier coefficient $C_{\vec{k}}\in {\mathbb{C}}^6$ can be obtained by the Fourier transform of all the discrete nodes, namely:

$C_{\vec{k}}\left(\vec{k}\right)={\displaystyle \sum _{n_1=0}^{N_1-1}}\cdots {\displaystyle \sum _{n_m=0}^{N_m-1}}u\left(0,{\theta }_{2,n_1},\dots ,{\theta }_{m+1,n_m}\right)e^{-i\left(k_1{\theta }_{2,n_1}+k_2{\theta }_{3,n_2}+\cdots k_m{\theta }_{m+1,n_m}\right)}(13)$

In Equation 13, for each phase angle element ${\mathit{\theta }}_{j+1,n_l}$, $(l\le m)$, the second subscript $n_l$ represents its index in the discrete set in Equation 12. Let $[D]$ and ${[D]}^{-1}$ be, respectively, the Fourier transform and the inverse Fourier transform operators in matrix form. According to the rotation characteristics of the Fourier transform, we have (McCarthy B. P. and Howell, 2021):

$\left[R_{-\rho }\right]=\left[D^{-1}\right]\left[Q_{-\rho }\right]\left[D\right](14)$

In Equation 14, $\left[Q_{-\rho }\right]$ is a diagonal matrix that rotates the Fourier coefficients $C_{\vec{k}}\in {\mathbb{C}}^6$ in the frequency domain by a fixed angle, i.e., $C_{\vec{k}}^{\prime }=C_{\vec{k}}{e}^{i⟨\mathit{k},-\mathit{\rho }⟩}$.

3.3 Continuation of invariant tori

The continuation process begins with the calculation of the first invariant torus. As is discussed at the beginning of this section, the initial guess is provided by linearizing the center manifold of the periodic orbit. Let $U_0$ be the column vector composed of all discrete states on the high-dimensional invariant torus in a certain order, and $U_1$ be the column vector composed of all the end states integrated from $U_0$ with a stroboscopic mapping time $T$. We then can incorporate the computation of quasi-periodic orbits into the standard predict-correction framework that can be solved by Newton’s method. Now, the free variable can be written as

$X=\left[U_0;T;\rho \right](15)$

and the invariant constraint can be given as

$F_{\text{Tori\hspace{0.17em}}}={\left[D\right]}^{-1}\left[Q_{-\rho }\right]\left[D\right]U_1-U_0(16)$

It should be noted that the dimension of free variable $X$ in Equation 15 is $m+1+{\prod }_{j=1}^m{N}_j$, while the dimension of invariant constraint in Equation 16 is ${\prod }_{j=1}^m{N}_j$. Thus, additional $m$ constraints $F_{\text{extra\hspace{0.17em}}}$ are required to ensure convergence of the first invariant torus to a member of the unique quasi-orbit family. These additional constraints $F_{\text{extra\hspace{0.17em}}}$ can be incorporated by referring to Lujan and Scheeres (2022).

Before proceeding the continuation process, some additional phase angle constraints should also be added to ensure the desired direction of continuation. This is because for an invariant torus, any point on the torus satisfies the above constraint equations, namely, $F_{\text{tori\hspace{0.17em}}}$ and $F_{\text{extra\hspace{0.17em}}}$. Therefore, if the phase angle of the invariant torus is not constrained, the new quasi-periodic solution may still lie on the same invariant torus, which is obviously trivial. Assuming that $\widehat{u}$ is the previously converged invariant torus, and $u$ is the current torus to be computed, then the $(m+1)$-dimensional phase angle constraints can be expressed as:

$F_{\text{phase\hspace{0.17em}}}={\left[⟨u,\frac{\partial \widehat{u}}{\partial {\theta }_1}⟩,⟨u,\frac{\partial \widehat{u}}{\partial {\theta }_2}⟩,\dots ,⟨u,\frac{\partial \widehat{u}}{\partial {\theta }_{m+1}}⟩\right]}^T(17)$

In Equation 17, the partial derivatives $\frac{\partial \widehat{u}}{\partial {\theta }_i}$, $(i=\mathrm{2,3},\dots ,m+1)$ can be directly obtained from the Fourier representation in Equation 11, while the partial derivative $\frac{\partial \widehat{u}}{\partial {\theta }_1}$ is given as (Baresi et al., 2018).

$\frac{\partial \widehat{u}}{\partial {\theta }_1}=\frac{1}{{\omega }_0}\left(f\left(t,x,\lambda \right)-{\displaystyle \sum _{i=2}^{m+1}}{\omega }_i\frac{\partial \widehat{u}}{\partial {\theta }_i}\right)(18)$

Equation 18 is derived through the total differential equation. By now, the augmented constraint matrix can be obtained as $F={\left[F_{\text{tori\hspace{0.17em}}}^{\text{T}},F_{\text{extra\hspace{0.17em}}}^{\text{T}},F_{\text{phase\hspace{0.17em}}}^{\text{T}}\right]}^{\text{T}}$, and the Jacobi matrix is written as $DF=\frac{\partial F}{\partial X}$. Since the phase angle constraint $F_{\text{phase\hspace{0.17em}}}$ only specifies the direction of continuation, the null space of the Jacobi matrix is still one-dimensional. A pseudo-arclength continuation or natural parameter continuation scheme can be applied to obtain a unique quasi-periodic orbit family. In the process of calculating the DF matrix, the FFT function of ${\mathrm{M}\mathrm{a}\mathrm{t}\mathrm{l}\mathrm{a}\mathrm{b}}^{Ⓡ}$ is used to speed up the calculation, and the related calculation of FFT and its derivation can be found in Johnson (2011). Also, unfolding parameters are leveraged to square the Jacobi matrix, which greatly improve the speed of matrix inversion, for more details, see Olikara (2016).

Two major challenges during the computation process are the determination of sampling nodes and overcoming the resonant singularity. As mentioned in Rosales et al. (2021), the convergence of Newton’s method does not necessarily imply that the solution accurately represents the torus. Assume $\mathit{u}\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)$ represents the exact invariant torus, which can be represented by a high-dimensional truncated Fourier series, denoted as $\widehat{\mathit{u}}\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)$. Thus, the convergence of the Newton’s method only imply that the function $\widehat{\mathit{u}}\left({\theta }_1,{\theta }_2,\dots ,{\theta }_{m+1}\right)$ satisfies the invariant constraints at the sampling points. It is essential to verify whether the accuracy of the Fourier approximation is maintained beyond the sampling points. Typically, the accuracy of the invariant torus constraints is tested at denser random sampling points, specifically:

$\underset{\widehat{\mathit{u}}}{\max }‖\left[{\mathit{D}}^{-1}\right]\left[{\mathit{Q}}_{-\mathit{\rho }}\right]\left[\mathit{D}\right]{\varphi }_T\left(\widehat{\mathit{u}}\right)-\widehat{\mathit{u}}‖<to{l}_2(19)$

In this manuscript, tol₂ is referred to as the torus verification accuracy, while the convergence accuracy for Newton’s method is denoted as tol₁. If Equation 19 is not satisfied, a higher-order Fourier series is added (with more samples accordingly), and then the differential correction is repeated. The number of sampling points is not fixed throughout the continuation process. Generally, the more complex, distorted, and larger the invariant torus, the more sampling points are required.

Since the quasi-periodic torus forms a Cantor set in the phase angle space, another challenge when resonance occurs is the singularity phenomenon as the continuation progresses. During singularity, the invariant torus becomes highly distorted, and nonlinearity sharply increases. A viable strategy is to dynamically adjust the step size within a specified range when singularity occurs until it crosses the resonant regions. Additionally, a hybrid pseudo-arclength continuation and natural parameter continuation are employed with various initial step sizes to enhance the reliability of results. The above analysis is summarized in Algorithm 1.

Algorithm 1.Adaptive continuation framework.

Require: The first converged invariant torus $X_0=\left[U_0;T;\rho \right]$

Ensure: A definite quasi-periodic family

1:  Let $to{l}_1$ and $to{l}_2$ be the convergence and accuracy tolerance of torus.

2:  Let $r_{l\text{,max\hspace{0.17em}}}$ be the maximum allowable ratio of the current torus accuracy to the average accuracy of the previous $l$ torus.

3:  $k\Leftarrow 1$; $k_1\Leftarrow 0$; $k_2\Leftarrow 0$; $r\Leftarrow []$.

4:  Obtain the initial quasi-periodic orbit, and check its accuracy $r_0$, if $r_0\le to{l}_2$ is satisfied, then $r\Leftarrow \left[r;r_0\right]$.

5:  while $k\le N$ do

6:   Update the initial guess of $k$th torus by moving a step size $ds$, and compute the new torus

7:  if $to{l}_1$ is not satisfied then

8:   scale $ds$ within a given range, $k_1\Leftarrow k_1+1$

9:   if $k_1\ge k_{1,\max }$ then

10:    stop

11:   end if

12:   continue

13:  else

14:   $k_1\Leftarrow 0$

15:  end if

16:  if $to{l}_2$ is not satisfied then

17:   $r_l\Leftarrow \frac{r(\text{\hspace{0.17em}end\hspace{0.17em}})}{mean(r(\text{\hspace{0.17em}end\hspace{0.17em}}-l+1:\text{\hspace{0.17em}end\hspace{0.17em}}))}$

18:   if $r_l\ge r_{l,\max }$ then

19:    turn to line 8

20:   else

21:    add $m$ sampling nodes for the current torus, $k_2\Leftarrow k_2+m$

22:    if $k_2\ge k_{2,\max }$ then

23:     stop

24:    end if

25:    continue

26:   end if

27:   continue

28:  else

29:   $k_2\Leftarrow 0$

30:  end if

31:  save $k$th torus, $k\Leftarrow k+1$, $r\Leftarrow \left[r;r_l\right]$

32:  end while

It should be noted that we introduce $r_l$ as an indicator to distinguish whether to add nodes or adjust the step size. Specifically, $r_l$ measures the rate of change of the torus verification accuracy $to{l}_2$. In general, a larger invariant torus requires more sampling nodes for an accurate characterization. Therefore, under the condition that the number of discrete nodes remains constant, the accuracy of the invariant torus will gradually decrease until a critical scale is reached. At this point, $to{l}_2$ is no longer satisfied indicating the need for additional nodes. Moreover, if there is a sudden decrease in the accuracy of the torus, it is likely due to an increase in the nonlinearity of the torus caused by its proximity to the resonance region. In such cases, a larger step size is necessary.

3.4 Linear stability of invariant tori

The stability analysis of an invariant set can provide insight into the evolution of the manifold structure in phase space. The methods we adopt here are based on those described in Jorba, 2001. The first term of invariant constraint in Equation 10 has the differential form as follows:

$A\left({\left[\cdot ,{\widehat{\theta }}^T\right]}^T\right)=R_{-\rho }\left[\frac{\partial {\varphi }_T\left(u\left({\left[\cdot ,{\widehat{\theta }}^T\right]}^T\right)\right)}{\partial u}\right](20)$

Jorba shows that, in reducible case, the eigenstructure of matrix $A\left({\left[\cdot ,{\widehat{\theta }}^T\right]}^T\right)$ in Equation 20 is equivalent to the eigenstructure of the Floquent matrix. Jorba also indicates that if there exist $n$ unrelated eigenvalues such that $B=\mathrm{diag}\left({\lambda }_1,\dots ,{\lambda }_n\right)$, then for each eigenvalue ${\lambda }_i$, ${\lambda }_i{e}^{⟨k,\rho ⟩}$ are also eigenvalues of matrix $B$, leading to the full eigenstructure forming circles in the complex plane. It should be noted that if the stroboscopic map is autonomous, then 1 is always an eigenvalue, which corresponds to the tangent direction of the torus.

The magnitude of the eigenvalues reflects the stability information. A convenient metric to measure the stability of invariant objects is the stability exponent, which is defined as:

$v_i=\frac{1}{2}\left(‖{\lambda }_i‖+‖\frac{1}{{\lambda }_i}‖\right)(21)$

If all the stability exponents $v_i$ in Equation 21 are equal to 1, the quasi-periodic orbit is stable, otherwise, the orbit is unstable and there are hyperbolic invariant manifolds.

4 Quasi-periodic orbits in the CR3BP

4.1 The center manifolds of 2:1 DRO and resonance ratios near their phase angles

During the continuation process of quasi-periodic orbits, encountering resonance among the phase angles characterizing the torus is inevitable. Resonance phenomena often signify the presence of other periodic orbits in phase space. If resonance regions are not stepped over during the continuation process, quasi-periodic orbits will deviate from their initial continuation direction. Consequently, the torus will undergo distortion as the continuation progresses, leading to a significant increase in nonlinearity. In previous literatures, while resonance phenomena were acknowledged during the computation of quasi-periodic orbits, it failed to elucidate whether the failure in orbit continuation resulted from encountering strong resonance or if the orbit family reached the boundaries of continuation. Consequently, the computed orbit families in these studies are essentially incomplete. Considering these factors, in this section, we initially elaborate on the resonance ratios that quasi-periodic motion near the 2:1 DRO might encounter.

Figure 3 depicts the relationship between the phase angles of center manifolds and the Jacobi energy of the entire DRO family, where the phase angles of the center manifolds align with the phase angles of the eigenvalues ${\alpha }_j$ of the monodromy matrix, providing initial values for the continuation of quasi-periodic orbits. The blue and red curves in Figure 3 represent the in-plane and vertical center manifolds, which, upon continuation, give rise to the in-plane quasi-periodic DROs and vertical quasi-periodic DROs (collectively referred to as the 2D quasi-DROs family). With the increase in Jacobi energy, the phase angles of the center manifolds both exhibit a pattern of initial increase followed by a decrease. It is noteworthy that when the Jacobi energy of the DROs reach 2.36766 (corresponding to a period of 6.24192), the phase angle of the normal center manifold becomes 0. This indicates a transition to a pair of hyperbolic invariant manifolds, signifying the instability of the DRO family. The yellow dashed line in the figure represents the Jacobi energy of the 2:1 DRO, which is 2.93052. At this energy level, the phase angles of the in-plane and normal center manifolds of the 2:1 DRO are 2.35822 and 1.46995, respectively, as indicated in Table 2.

Figure 3

Figure 3. The relationship between the phase angles of center manifolds and the Jacobi energy of the entire DRO family.

Table 2

Table 2. The center manifolds associated with the 2:1DRO in Earth-Moon system.

Douskos et al. (2007) investigated the impact of resonances on the in-plane stability region in the vinctiny of DROs. Their findings revealed that, with the exception of the period-tripling DRO (P3DRO), all resonant orbits bifurcating from the DROs are situated within the stable region and gradually approach the boundary of the stable region, namely, P3DRO. Hence, for the in-plane 2D quasi-DROs family around the 2:1 DRO, as the continuation progresses, the quasi-periodic family will be impeded by the resonance ratio 1/3, and the phase angle gradually decreases and eventually approaches $2\pi /3$. This assertion will be substantiated in subsequent discussions. As for the normal 2D quasi-DROs family, resonances near the initial phase angles include 1/4, 1/5, and 1/6, among others. Therefore, the final resonance ratio that impedes this family depends on the variation in phase angle of the normal center manifold.

4.2 In-plane 2D quasi-DROs family

2D quasi-periodic orbits can be characterized by two phase angles. The first phase ${\theta }_1$ represents the phase of the periodic orbit, while the second phase, denoted as ${\theta }_2$, indicates the torus on the stroboscopic mapping map. Starting from any point on the torus, the final state will remain on the same torus after a stroboscopic mapping time $T$. This property holds for torus on the stroboscopic mapping map with different ${\theta }_1$ and serves as the foundation for computing quasi-periodic orbits. As outlined in Section 3, to obtain a family of 2D quasi-periodic orbits, an additional constraint is required. Given that verifying whether an orbit reaches the boundary often requires analysis using the Poincaré section, and the scattered points on the Poincaré section typically represent projections of trajectory flows at specific Jacobi energy levels in the CR3BP, we choose the Jacobi energy as the additional constraint in this paper.

Figure 4 depicts the relationship between the stroboscopic mapping period and the rotation angle ${\rho }_{\mathit{plane}}$ of the in-plane 2D quasi-DROs family with fixed Jacobi energy. As the continuation progresses, the stroboscopic mapping time $T$ of the 2D quasi-DRO orbit family gradually increases, ranging from $\pi$ to 3.15904. Simultaneously, the rotation phase angle ${\rho }_{\mathit{plane}}$ gradually decreases, reducing from 2.35822 to 2.21959. During the continuation process, the in-plane 2D quasi-DRO family encounters various resonance regions, as indicated by the red dashed lines in Figure 4. At these points, ${\rho }_{\mathit{plane}}$ resonates with $2\pi$ (the first phase angle of the 2D quasi-periodic orbit). As the quasi-DROs family approaches the continuation boundary, the density of resonance regions increases. This phenomenon may be related to resonant orbits bifurcating from the DROs. When encountering resonance regions like 7/19, 13/36, and 9/25, the continuation algorithm proposed in this paper easily navigates through these regions by dynamically adjusting the step size. These resonance regions have narrow “gaps” in the orbit family. Therefore, they are referred to as weak resonance regions. On the other hand, encountering resonance regions like 4/11, 5/14, and 4/17 is relatively more challenging. It requires more frequent dynamic adjustments, and these three resonance regions have wider “gaps” in the orbit family. Consequently, they are referred to as strong resonance regions. In this study, the in-plane 2D quasi-DROs family is eventually blocked by the resonance region of 6/17.

Figure 4

Figure 4. The relationship between the stroboscopic mapping period and the rotation angle ρ_in-plane of the plane 2D quasi-DRO family.

Figure 5 shows four in-plane 2D quasi-DROs with different rotation angles, where the red curves represent the invariant torus on the stroboscopic mapping defined at two vertically crossing points of 2:1 DRO, the yellow curve corresponds to the DRO orbit, and the blue region represents the densely covered area of the quasi-periodic orbits. It can be observed that the size of the invariant torus increases as the rotation angle ${\rho }_{\mathit{plane}}$ gradually decreases. The coverage of the in-plane quasi-periodic orbits in the $xy$-plane expands progressively, forming a “ring-shaped” region, and the invariant torus is tangentially aligned with the inner and outer boundaries of this “ring-shaped” region. When the size of the invariant torus is small, it appears nearly elliptical. As the torus size increases, the outer side of the invariant torus starts to flatten, gradually assuming a more elongated shape. Subsequently, the growth of the torus in the $x$-direction gradually slows down, but it continues to expand in the $y$-direction. The continuation of the torus comes to an end when encountering the resonance ratio 6/17. In this study, the number of sample points required to character the invariant torus increases from 11 to 251.

Figure 5

Figure 5. Four members of the in-plane 2D quasi-DROs family. (A) T = 3.14426, ρ_plane = 2.34052. (B) T = 3.14547, ρ_plane = 2.33223. (C) T = 3.15360, ρ_plane = 2.26970. (D) T = 3.14426, ρ_plane = 2.34052.

Record the intersections between the in-plane 2D quasi-DROs and the Poincaré section $\mathrm{\Sigma }:y=0$, and project the state variables $(x,\dot{x})$ onto the two-dimensional Poincaré section, as shown in Figure 6A, where the red pentagram in Figure 6B represents the P3DRO, and its trajectory is plotted in Figure 6C. It can be observed that the entire quasi-periodic orbits forms a “shell-like” family of curves on the Poincaré section. Some curves exhibit significant gaps between them, indicating their proximity to resonance regions. Additionally, a small gap exists between the in-plane 2D quasi-periodic family and the P3DRO. Further investigation reveals that within this gap, The quasi-periodic orbits are densely distributed along P3DRO, whereas the distribution is sparse in regions far from P3DRO. And the invariant torus becomes highly nonlinear and distorted, as depicted in Figure 6C, deviating from the original direction of the quasi-periodic orbit family. Due to its proximity to the 1/3 resonance, the continuation algorithm easily encounters singularities, leading to program failure. Therefore, in a loose sense, the planar 2D quasi-periodic orbit family does reach its maximum boundary, namely, P3DRO.

Figure 6

Figure 6. (A) The projection of the in-plane 2D quasi-DROs family on the Poincaré section, and the red pentagram represents the P3DROs. (B) the P3DROs with equal Jacobi energy as 2:1 DROs. (C) The twisted torus between the in-plane 2D quasi-DROs family and the P3DROs.

4.3 Vertical 2D quasi-DROs family

Figure 7 illustrates the relationship between the stroboscopic mapping time $T$ and the rotation angle ${\rho }_{\mathit{vertical}}$ of the vertical 2D quasi-DROs family when the Jacobi energy is 2.93052. As the continuation progresses, the stroboscopic mapping time gradually decreases, while the phase angle gradually increases. More specifically, the former decreases from $\pi$ to 3.12049, and the latter increases from 1.46995 to 1.56322. The phase space of the vertical 2D quasi-DROs orbit family also exhibits a series of resonance regions. It is important to note that although some resonance regions may appear closely spaced in the phase space, this does not necessarily imply a small variation in the size of the quasi-DROs before and after the resonance regions. In this case, the quasi-periodic orbits near the resonance orbits occupy the resonance region, hindering the continuation of the quasi-DROs and requiring dynamic adjustments in step size to overcome the resonance regions. As the continuation progresses, the resonance regions become more densely packed, and higher frequencies are encountered, which may be related to the resonant orbits bifurcated from the DROs. Additionally, based on simulation results, when ${\rho }_{\mathit{vertical}}<1.51$, the $z$-axis size of the vertical 2D quasi-DROs family increases rapidly, ranging from 0 to approximately 0.05. However, for ${\rho }_{\mathit{vertical}}>1.51$, the growth in the $z$-axis size of the quasi-DROs family becomes sluggish, only progressing from 0.05 to around 0.06. When the height of the quasi-DROs approaches 0.06, even though ${\rho }_{\mathit{vertical}}$ continues to rise, it nearly stagnates. The vertical 2D quasi-DROs is ultimately impeded by the resonance region 1/4, agreeing with Figure 3.

Figure 7

Figure 7. The relationship between the stroboscopic mapping period and the rotation angle ρ_vertical of the vertical 2D quasi-DROs family.

Figure 8 shows four quasi-periodic orbits at different rotation angles. As the continuation progresses, the size of the vertical quasi-periodic orbits in the $z$-direction gradually increases, forming a “striped” region. And the invariant torus is tangent to the upper and lower boundaries of the “striped” region. At the far earth side of the vertically crossing point of the 2:1 DRO, the invariant torus is figure-eight-shaped. With the increase in torus size, the upper and lower edges of the torus start to straighten, indicating a slowdown in the growth of the orbit family in the $z$-direction. Subsequently, as the continuation proceeds, the $z$-size remains relatively constant, but the torus becomes more nonlinear, resulting in a more twisted and complex three-dimensional spatial structure of the “striped” region. It is important to note that, with the increase in torus size and distortion, more sampling points are needed for accurate representation. In this study, the number of sampling points required for the red invariant torus in Figure 8 increases from 11 to 791. It can be noted that as the torus approaches resonance 1/4, the quasi-periodic orbits will no longer be uniformly distributed throughout the entire bounded region. Instead, they primarily concentrate around a four-loop periodic orbit, referred to as P4DRO in this paper.

Figure 8

Figure 8. Four members of the vertical 2D quasi-DROs family. (A) T = 3.13703, ρ_vertical = 1.48652. (B) T = 3.13549, ρ_vertical = 1.49223. (C) T = 3.12613, ρ_vertical = 1.53127. (D) T = 3.12049, ρ_vertical = 1.56322.

The analysis process, utilizing the Poincaré section technique, for the vertical 2D quasi-DROs is the same as that for the in-plane 2D quasi-DROs. The difference lies in projecting the state variables $(x,\dot{z},z)$ onto the three-dimensional Poincaré section, as shown in Figure 9B, where the red pentagram in Figure 9B represents the P4DRO, and its trajectory is shown in Figure 9A. The same phenomenon, such as significant gaps between family members, has occurred. Additionally, we have observed that certain curves are discontinuous, suggesting the existence of other resonant DRO orbits. From Figure 9B, it can be seen that P4DRO is indeed the boundary of the 2D vertical quasi-DROs family, which demonstrates the effectiveness of our algorithm.

Figure 9

Figure 9. (A) The projection of the vertical 2D quasi-DROs family on the Poincaré section, and the red pentagram represents the P4DROs. (B) the P4DROs with equal Jacobi energy as 2:1 DROs.

Note that when the projections of the vertical 2D quasi-DROs and the P4DRO on the Poincaré section are positioned to the left of the Moon, $\dot{x}$ and $\dot{z}$ are zero when $z$ takes the maximum value. Therefore, the vertical 2D quasi-periodic orbit can be characterized by three characteristic parameters, $z_{\mathit{max}}$ and its corresponding $x$ and $\dot{y}$, labeled as $\mathrm{X}=\left[z_{\max },x,\dot{y}\right]$. By the way, since $z_{\mathit{max}}=0$ for the 2D in-plane quasi-DROs, then $\mathrm{X}=\left[x_{\max },\dot{y}\right]$ according to Figure 6A. To confirm that chaos emerges beyond P4DRO when the Jacobi energy is 2.93052, the Poincaré section technology is utilized again. We firstly discretize the parametersc $(x,z)$ directly, and $\dot{y}$ can be determined through the Jacobi energy defined in Equation 4. Then, for each initial state characterized by the aforementioned three variables, integrate it for 2,000 revolutions. If the orbit remains bounded or does not collide with celestial bodies, it is considered a quasi-periodic orbit.

Figure 10 validates our previous conclusion that, when the Jacobi energy is 2.93052, P3DRO and P4DRO are the boundaries of the in-plane 2D quasi-DROs and vertical 2D quasi-DROs, respectively. It should be noted that P4DRO only serves as the boundary of the vertical 2D quasi-DROs when the Jacobi energy is colse to 2.93052. The phase angle of the vertical center manifold of DROs undergoes significant changes when the Jacobi energy is in the range of 2.8–3.2. Therefore, if the phase angle of the vertical center manifold of DRO approaches 1/5 or 1/6, the boundary of vertical 2D quasi-periodic orbits might be associated with other resonant DROs. However, this aspect is beyond the scope of this study.

Figure 10

Figure 10. Characteristic parameter space of quasi-DROs at a Jacobi constant of 2.93052. Magenta and yellow scatters represent 2D vertical and in-plane quasi-periodic families, respectively. Red and green pentagram represent the P3DRO and P4DRO, respectively.

5 Stability analysis

The stability analysis of these two families is conducted in this section. For Hamiltonian systems with six-dimensional state space, the Floquent matrix of invariant objects is symplectic, leading to all eigenvalues appearing in the form of three recipocal pairs. Since the stroboscopic is autonomous in the CR3BP, all torus has a pair of eigenvalues equal to 1, and futher, equal to $e^{⟨k,\rho ⟩}$. Therefore, our focus will be on the types of the remaining two pairs of eigenvalues. Eigenvalues with modulo 1 represent either the torus itself or the direction of other center manifolds, while eigenvalues whose modulus is not 1 indicate the direction of hyperbolic manifolds. As stated in Olikara and Scheeres (2012), the maximum dimension of the torus in $n$-dimensional autonomous system is up to $n/2$, then we can induce that all the eigenvalues of the 3D torus (if exists) in the CR3BP are equal to 1. Therefore, there are always two pairs of eigenvalues equal to 1 for 2D torus in the CR3BP, and the remaining pair of eigenvalues, after our numerical computation, are conjugate complex whose modulus is 1. And the phase angles of the center manifolds of these two types of 2D quasi-DROs are displayed in Figure 11.

Figure 11

Figure 11. (A) Variation of the phase angles of the center manifolds for the in-plane 2D quasi-DROs with respect to ${\rho }_{\mathit{plane}}$. (B) Variation of the phase angles of the center manifolds for the vertical 2D quasi-DROs with respect to ${\rho }_{\mathit{vertical}}$.

Therefore, both types of quasi-periodic orbits are long-term stable, offering more orbit choices for mission design.

6 Conclusion

In this paper, we have explored the dynamical structure near 2:1 resonant distant retrograde orbit in the circular restricted three-body problem. An adaptive continuation framework with adjustable step size and sampling nodes is proposed to solve the challenges encountered in the calculation of quasi-periodic families. The 2D plane and vertical quasi-periodic families are computed, respectively. Since the family of quasi-periodic orbits is Cantorian, a series of resonant regions will be encountered as the family evolves, leading to some “gaps” in the continuation process and the phenomenon of singularity of the invariant torus. In addition, we have observed that denser regions will occur in the coverage of quasi-periodic orbits when the rotation angle is close to resonance.

As a preliminary attempt, the geometric boundary that can be reached by the 2D quasi-periodic families in the CR3BP has been analyzed using the Poincaré map method. It has been verified by numerical simulation that the P3DRO and P4DRO are geometric boundaries of the 2D in-plane and vertical quasi-DROs. These dynamical structures greatly enrich the knowledge of the phase space near 2:1 resonant DRO, which is a major destination for cislunar space exploration.

As the CR3BP represents the most simplified model of multi-body dynamics, it disregards the eccentricity of the Moon’s orbit and the perturbative effects of the Sun. Future research will extend to examine the existence of quasi-periodic orbit families when these factors are incorporated. Additionally, future work will focus on applying these dynamical structures to space missions, including formation flying, transfer trajectories, and eclipse mitigation.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

MW: Writing–original draft, Writing–review and editing, Methodology, Validation, Conceptualization. CY: Methodology, Writing–review and editing. YS: Methodology, Writing–review and editing. HZ: Conceptualization, Funding acquisition, Methodology, Supervision, Validation, Writing–review and editing.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA30010200).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Reference

Baresi, N., Olikara, Z. P., and Scheeres, D. J. (2018). Fully numerical methods for continuing families of quasi-periodic invariant tori in astrodynamics. J. astronautical Sci. 65, 157–182. doi:10.1007/s40295-017-0124-6

CrossRef Full Text | Google Scholar

Bezrouk, C. J., and Parker, J. (2014). “Long duration stability of distant retrograde orbits,” in AIAA/AAS Astrodynamics Specialist Conference, San Diego, CA, August 4–7, 2014.4424

CrossRef Full Text | Google Scholar

Colaprete, A., Andrews, D., Bluethmann, W., Elphic, R. C., Bussey, B., Trimble, J., et al. (2019). “An overview of the volatiles investigating polar exploration rover (viper) mission,” in AGU fall meeting abstracts, Bowie, MD, November 3–4, 2021.

Google Scholar

Conte, D., Di Carlo, M., Ho, K., Spencer, D. B., and Vasile, M. (2018). Earth-mars transfers through moon distant retrograde orbits. Acta Astronaut. 143, 372–379. doi:10.1016/j.actaastro.2017.12.007

CrossRef Full Text | Google Scholar

Douskos, C., Kalantonis, V., and Markellos, P. (2007). Effects of resonances on the stability of retrograde satellites. Astrophysics Space Sci. 310, 245–249. doi:10.1007/s10509-007-9508-6

CrossRef Full Text | Google Scholar

Farquhar, R. W., and Kamel, A. A. (1973). Quasi-periodic orbits about the translunar libration point. Celest. Mech. 7, 458–473. doi:10.1007/bf01227511

CrossRef Full Text | Google Scholar

Folta, D., Bosanac, N., Cox, A., and Howell, K. (2016). “The lunar icecube mission design: construction of feasible transfer trajectories with a constrained departure,” in AIAA Space Flight Mechanics Meeting, Napa, CA, February 14–18, 2016

Google Scholar

Gardner, T., Cheetham, B., Parker, J., Forsman, A., Kayser, E., Thompson, M., et al. (2023). Capstone: a summary of a highly successful mission in the cislunar environment.

Google Scholar

Gómez, G., Masdemont, J., and Simó, C. (1998). Quasihalo orbits associated with libration points. J. Astronautical Sci. 46, 135–176. doi:10.1007/bf03546241

CrossRef Full Text | Google Scholar

Hou, X., and Liu, L. (2011). On quasi-periodic motions around the collinear libration points in the real earth–moon system. Celest. Mech. Dyn. Astronomy 110, 71–98. doi:10.1007/s10569-011-9340-8

CrossRef Full Text | Google Scholar

Hou, X. Y., and Liu, L. (2010). On quasi-periodic motions around the triangular libration points of the real earth–moon system. Celest. Mech. Dyn. Astronomy 108, 301–313. doi:10.1007/s10569-010-9305-3

CrossRef Full Text | Google Scholar

Hou, X. Y., Xin, X., Scheeres, D. J., and Wang, J. (2015). Stable motions around triangular libration points in the real earth–moon system. Mon. Not. R. Astron Soc. 454, 4172–4181. doi:10.1093/mnras/stv2216

CrossRef Full Text | Google Scholar

Johnson, S. G. (2011). Notes on fft-based differentiation. MIT Appl. Math. Tech. Rep.

Google Scholar

Jorba, A. (2001). Numerical computation of the normal behaviour of invariant curves of n-dimensional maps. Nonlinearity 14, 943–976. doi:10.1088/0951-7715/14/5/303

CrossRef Full Text | Google Scholar

Jorba, À., Jorba-Cuscó, M., and Rosales, J. J. (2020). The vicinity of the earth–moon l1 point in the bicircular problem. Celest. Mech. Dyn. Astronomy 132, 11. doi:10.1007/s10569-019-9940-2

CrossRef Full Text | Google Scholar

Jorba, A., and Masdemont, J. (1999). Dynamics in the center manifold of the collinear points of the restricted three body problem. Phys. D. Nonlinear Phenom. 132, 189–213. doi:10.1016/s0167-2789(99)00042-1

CrossRef Full Text | Google Scholar

Jorba, A., and Nicolás, B. (2020). Transport and invariant manifolds near l3 in the earth-moon bicircular model. Commun. Nonlinear Sci. Numer. Simul. 89, 105327. doi:10.1016/j.cnsns.2020.105327

CrossRef Full Text | Google Scholar

Jorba, A., and Olmedo, E. (2009). On the computation of reducible invariant tori on a parallel computer. SIAM J. Appl. Dyn. Syst. 8, 1382–1404. doi:10.1137/080724563

CrossRef Full Text | Google Scholar

Lujan, D., and Scheeres, D. J. (2022). The earth-moon l2 quasi-halo orbit family: characteristics and manifold applications. AIAA SCITECH 2022 Forum, 2459. doi:10.2514/1.G006681

CrossRef Full Text | Google Scholar

McCarthy, B., and Howell, K. (2021). “Quasi-periodic orbits in the sun-earth-moon bicircular restricted four-body problem,” in Proceedings of the 31st AAS/AIAA Space Flight Mechanics Meeting, AAS Paper.

Google Scholar

McCarthy, B. P., and Howell, K. C. (2019). “Trajectory design using quasi-periodic orbits in the multi-body problem,” in 29th AAS/AIAA Space Flight Mechanics Meeting, Hawaii, USA.

Google Scholar

McCarthy, B. P., and Howell, K. C. (2021). Leveraging quasi-periodic orbits for trajectory design in cislunar space. Astrodynamics 5, 139–165. doi:10.1007/s42064-020-0094-5

CrossRef Full Text | Google Scholar

Olikara, Z. P. (2016). Computation of quasi-periodic tori and heteroclinic connections in astrodynamics using collocation techniques. Boulder, USA: University of Colorado at Boulder.

Google Scholar

Olikara, Z. P., and Howell, K. C. (2010). Computation of quasi-periodic invariant tori in the restricted three-body problem. Pap. No. AAS, 1–15.

Google Scholar

Olikara, Z. P., and Scheeres, D. J. (2012). Numerical method for computing quasi-periodic orbits and their stability in the restricted three-body problem. Adv. Astronautical Sci. 145, 911–930.

Google Scholar

Rosales, J. J., Jorba, À., and Jorba Cuscó, M. (2021). Families of halo-like invariant tori around $$L_2$$ in the earth-moon bicircular problem. Celest. Mech. Dyn. Astronomy 133, 16. doi:10.1007/s10569-021-10012-0

CrossRef Full Text | Google Scholar

Sánchez, J., and Net, M. (2013). A parallel algorithm for the computation of invariant tori in large-scale dissipative systems. Phys. D. Nonlinear Phenom. 252, 22–33. doi:10.1016/j.physd.2013.02.008

CrossRef Full Text | Google Scholar

Schilder, F., Osinga, H. M., and Vogt, W. (2005). Continuation of quasi-periodic invariant tori. SIAM J. Appl. Dyn. Syst. 4, 459–488. doi:10.1137/040611240

CrossRef Full Text | Google Scholar

Smitherman, D. V., and Griffin, B. N. (2014). “Habitat concepts for deep space exploration,” in AIAA Space 2014 Conference and Exposition, San Diego, CA, August 4–7, 2014.

CrossRef Full Text | Google Scholar

Song, Y.-J., Kim, Y.-R., Bae, J., Park, J.-i., Hong, S., Lee, D., et al. (2021). Overview of the flight dynamics subsystem for korea pathfinder lunar orbiter mission. Aerospace 8, 222. doi:10.3390/aerospace8080222

CrossRef Full Text | Google Scholar

Stramacchia, M., Colombo, C., and Bernelli-Zazzera, F. (2016). Distant retrograde orbits for space-based near earth objects detection. Adv. Space Res. 58, 967–988. doi:10.1016/j.asr.2016.05.053

CrossRef Full Text | Google Scholar

Topputo, F. (2013). On optimal two-impulse earth–moon transfers in a four-body model. Celest. Mech. Dyn. Astronomy 117, 279–313. doi:10.1007/s10569-013-9513-8

CrossRef Full Text | Google Scholar

Wang, W., Shu, L., Liu, J., and Gao, Y. (2019). Joint navigation performance of distant retrograde orbits and cislunar orbits via liaison considering dynamic and clock model errors. Navigation 66, 781–802. doi:10.1002/navi.340

CrossRef Full Text | Google Scholar

Williams, J., Dawn, T. F., and Batcha, A. L. (2023). “A history of orion mission design, copernicus software development, and the artemis i trajectory,” in AAS/AIAA astrodynamics specialist conference, Big Sky, MT, August 13–17, 2023, 23–241.

Google Scholar