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 M = m 1 + m 2 L = a 12 T = a 12 3 / G m 1 + m 2 (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 μ = m 2 m 1 + m 2 as the mass ratio of the Earth-Moon system, then the Earth and the Moon are located at ( − μ , 0,0 ) and ( 1 − μ , 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). x ̈ − 2 y ̇ = ∂ Ω 3 ∂ x y ̈ + 2 x ̇ = ∂ Ω 3 ∂ y z ̈ = ∂ Ω 3 ∂ z (2)

In Equation 2, the effective potential Ω 3 is defined as Ω 3 x , y , z = 1 2 x 2 + y 2 + 1 − μ r E S + μ r M S + 1 2 μ 1 − μ (3)

In Equation 3, r E S r M S 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.

The CR3BP is an autonomous system with a unique energy integral, i.e., the Jacobi constant as displayed in Equation 4, C = − x ̇ 2 + y ̇ 2 + z ̇ 2 + 2 Ω 3 x , y , z (4)

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 x z plane in the Earth-Moon rotating frame, the two intersections of one particular DRO satisfy the vertically crossing condition. By choosing the plane Σ : y = 0 as the Poincaré map, the particular DRO can be represented by a two-dimensional state, denoted as x 0 , y ̇ 0 . 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 δ y t δ x ̇ t = Φ 2,1 Φ 2,5 y ̇ t Φ 2,4 Φ 2,5 x ̈ t δ x 0 δ y ̇ 0 δ t (5) where Φ represents the 6 × 6 state transition matrix. The subscripts in Equation 5 specify the particular row and column elements of Φ respectively.

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.

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 u ̄ θ 1 = r θ 1 v θ 1 : T 1 → 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 Φ T u ̄ θ 1 denote the state transition matrix of a periodic orbit integrated from initial point u ̄ θ 1 . Representing the complex eigenvalue as λ j + 1 = e i α j + 1 and the associated complex eigenvector as y α j + 1 , the initial guess of the invariant torus can be expressed as (McCarthy B. P. and Howell, 2021): u θ 1 , θ 2 , … , θ m + 1 = u ̄ θ 1 + κ Re ∑ j = 1 m e i θ j + 1 y α j + 1 (6) In Equation 6, κ is a small quantity. The initial guess of rotation angle and frequency associated with each dimension in the phase space can be given as ρ = 2 π , α 2 , … , α m + 1 ω = 2 π T , α 2 T , … , α m + 1 T (7)

It should be noted that all frequencies are nonresonant in Equation 7, i.e., for any k ∈ Z n \ { 0 } , k ⋅ ω ≠ 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.

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): x ̇ = f t , x , λ (8) where x ∈ R n , and λ ∈ R q is the external disturbance parameter. For the CR3BP, the system dynamics is time independent, and the external disturbance parameter λ 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 P : T m + 1 → R 6 such that for any point on the torus, u θ 1 , θ 2 , … , θ m + 1 is a differential isomorphic of the quasi-periodic motion. Suppose that ϕ T x 0 is the flow function generated by Equation 8, then if the initial state u θ 1 , θ 2 , … , θ m + 1 is propagated for a stroboscopic mapping time T = 2 π ω 1 , it gives ϕ T u θ 1 , θ 2 , … , θ m + 1 = u u θ 1 + 2 π , θ 2 + ω 2 T , … , θ m + 1 + ω m + 1 T = u θ 1 , θ ̂ T + ρ ̂ T (9) where θ ̂ = θ 2 , θ 3 , … , θ m + 1 T , and ρ ̂ = 2 π ω 1 ω 2 , ω 3 , … , ω m + 1 T is the rotation vector. One can conclude from Equation 9 that there is an invariant map with respect to the rotation vector ρ . As such, if a rotation operator R − ρ is defined to eliminate the rotation along the directions of ρ , the invariant constraint can be formulated as R − ρ ϕ T u ⋅ , θ ̂ T T − u ⋅ , θ ̂ T T = 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 θ = ∑ k 1 = 0 N 1 − 1 ⋯ ∑ k m = 0 N m − 1 C k ⃗ e i k 1 θ 2 + k 2 θ 3 + ⋯ k m θ m + 1 (11) where C k ⃗ ∈ C 6 is the Fourier coefficient. k ⃗ = k 1 , k 2 , … , k m T is the harmonic set of the high-dimensional Fourier series. N j is the number of discrete nodes for each phase component θ j + 1 , ( j = 1,2 , … , m ) . To this end, θ j + 1 is uniformly discretized into elements in the set shown as follows: θ j + 1 ∈ 0 2 π N j ⋯ 2 π N j − 1 N j (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 ∏ j = 1 m N j nodes. The Fourier coefficient C k ⃗ ∈ C 6 can be obtained by the Fourier transform of all the discrete nodes, namely: C k ⃗ k ⃗ = ∑ n 1 = 0 N 1 − 1 ⋯ ∑ n m = 0 N m − 1 u 0 , θ 2 , n 1 , … , θ m + 1 , n m e − i k 1 θ 2 , n 1 + k 2 θ 3 , n 2 + ⋯ k m θ m + 1 , n m (13) In Equation 13, for each phase angle element θ j + 1 , n l , ( l ≤ 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): R − ρ = D − 1 Q − ρ D (14) In Equation 14, Q − ρ is a diagonal matrix that rotates the Fourier coefficients C k ⃗ ∈ C 6 in the frequency domain by a fixed angle, i.e., C k ⃗ ′ = C k ⃗ e i ⟨ k , − ρ ⟩ .

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 = U 0 ; T ; ρ (15) and the invariant constraint can be given as F Tori  = D − 1 Q − ρ D U 1 − U 0 (16) It should be noted that the dimension of free variable X in Equation 15 is m + 1 + ∏ j = 1 m N j , while the dimension of invariant constraint in Equation 16 is ∏ j = 1 m N j . Thus, additional m constraints F extra  are required to ensure convergence of the first invariant torus to a member of the unique quasi-orbit family. These additional constraints F extra  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 tori  and F extra  . 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 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 phase  = u , ∂ u ̂ ∂ θ 1 , u , ∂ u ̂ ∂ θ 2 , … , u , ∂ u ̂ ∂ θ m + 1 T (17) In Equation 17, the partial derivatives ∂ u ̂ ∂ θ i , ( i = 2,3 , … , m + 1 ) can be directly obtained from the Fourier representation in Equation 11, while the partial derivative ∂ u ̂ ∂ θ 1 is given as (Baresi et al., 2018). ∂ u ̂ ∂ θ 1 = 1 ω 0 f t , x , λ − ∑ i = 2 m + 1 ω i ∂ u ̂ ∂ θ i (18)

Equation 18 is derived through the total differential equation. By now, the augmented constraint matrix can be obtained as F = F tori  T , F extra  T , F phase  T T , and the Jacobi matrix is written as D F = ∂ F ∂ X . Since the phase angle constraint F phase  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 M a t l a 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 u θ 1 , θ 2 , … , θ m + 1 represents the exact invariant torus, which can be represented by a high-dimensional truncated Fourier series, denoted as u ̂ θ 1 , θ 2 , … , θ m + 1 . Thus, the convergence of the Newton’s method only imply that the function u ̂ θ 1 , θ 2 , … , θ m + 1 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: max u ̂ D − 1 Q − ρ D ϕ T u ̂ − u ̂ < t o 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 = U 0 ; T ; ρ

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 t o 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, t o 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.

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 ⋅ , θ ̂ T T = R − ρ ∂ ϕ T u ⋅ , θ ̂ T T ∂ u (20)

Jorba shows that, in reducible case, the eigenstructure of matrix A ⋅ , θ ̂ T T 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 = diag λ 1 , … , λ n , then for each eigenvalue λ i , λ i e ⟨ k , ρ ⟩ 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 = 1 2 λ i + 1 λ i (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.

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 α 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.

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 π / 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.

2D quasi-periodic orbits can be characterized by two phase angles. The first phase θ 1 represents the phase of the periodic orbit, while the second phase, denoted as θ 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 θ 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 ρ 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 π to 3.15904. Simultaneously, the rotation phase angle ρ 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, ρ plane resonates with 2 π (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 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 ρ plane gradually decreases. The coverage of the in-plane quasi-periodic orbits in the x y -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.

Record the intersections between the in-plane 2D quasi-DROs and the Poincaré section Σ : y = 0 , and project the state variables ( x , 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 7 illustrates the relationship between the stroboscopic mapping time T and the rotation angle ρ 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 π 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 ρ 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 ρ 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 ρ 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 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.

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 , 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.

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, x ̇ and z ̇ are zero when z takes the maximum value. Therefore, the vertical 2D quasi-periodic orbit can be characterized by three characteristic parameters, z max and its corresponding x and y ̇ , labeled as X = z max , x , y ̇ . By the way, since z max = 0 for the 2D in-plane quasi-DROs, then X = x max , y ̇ 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 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.