Stability of spatial orbits around Earth–Moon triangular libration points

ABSTRACT

The stability of spatial orbits around the triangular libration points in the ephemeris model of the Earth–Moon system is studied. Five contributions are made: (1) practical stable spatial orbits in the ephemeris Earth–Moon system lasting thousands of years or even longer are first reported, and spatial stable regions are identified. (2) The mechanism that shapes the boundaries of the spatial stable regions is investigated, and is found to be related to resonances among the precession rates of the lunar orbit, the precession rates of the small body, the mean orbital motion of the Sun, and the libration frequency of the co-orbital motion. (3) Influence on the spatial stable regions from the solar radiation pressure is studied. It is found that the spatial stable region generally shrinks with increasing solar radiation pressure strength. Dust grains with sizes of millimetres in magnitude or smaller generally escape in hundreds of years while objects with larger sizes can stay for thousands of years or even longer; (4) difference between the bi-circular problem model and the ephemeris model in describing the spatial stable regions is presented. (5) The observation of possible objects in the spatial stable orbits is discussed. With the public’s growing interest in the cislunar space, the current study is a good attempt to enhance the understanding of the practical orbital dynamics in the cislunar space.

minor planets, asteroids: general

1 INTRODUCTION

There is a long history of studies on the stability of motions around the triangular libration points (TLPs). In the circular restricted three-body problem (CRTBP), these points are stable in the linear sense if the mass parameter |$\mu$| is smaller than Routh’s critical value |$\mu _0$| (Szebehely 1967). In the planar case, the same stability condition holds in the non-linear sense except for two specific values of |$\mu \lt \mu _0$| (Alfriend 1970, 1971). To investigate the practical stable regions around these points, analytical solutions truncated at high orders are constructed in both Cartesian coordinates (Celletti et al. 1991; Lei et al. 2013; Liang et al. 2018) and polar coordinates (Giorgilli et al. 1997; Tan et al. 2021), from both the approach of directly solving the equations of motion (Deprit et al. 1965; Lei et al. 2014) and the approach of Hamiltonian normalization (Deprit, Henrard & Rom 1967; Gabern, Jorba & Locatelli 2005). Along with the analytical approach, families of periodic orbits around TLPs are also extensively investigated by researchers numerically (Henrard 2002; Hou et al. 2009).

The Sun–Jupiter system gained the most attention from researchers because natural bodies (the Trojans) are found moving around its TLPs. The existence of the Trojans holds important clues to the evolution history of the Solar system (Morbidelli et al. 2005; Nesvorný, Vokrouhlický & Morbidelli 2013). Compared with the idealized CRTBP model, the Sun–Jupiter system is perturbed by the orbit eccentricity of Jupiter and the gravities of the outer planets and Mars. Concerning orbital dynamics around the TLPs of the perturbed Sun–Jupiter system, a huge amount of work has been carried out. To list just a few, please see Di Sisto, Ramos & Beaugé (2014), Hou, Scheeres & Liu (2014a, 2016), Li et al. (2023a, b), and Vokrouhlický et al. (2024).

Much less work is carried out for the TLPs of the Earth–Moon system. Compared with the CRTBP model, perturbations from the Moon’s orbit eccentricity and the Sun’s gravity play important roles in shaping the dynamics of the TLPs in the real Earth–Moon system (Park et al. 2024). The lifetime of objects in the cislunar space is closely related to the obliquity of the Moon’s orbit as well (Cinelli et al. 2022). Assuming a coplanar configuration of the Sun–Earth–Moon system, Kamel et al. (1970) constructed analytical solutions up to fourth- order using the Hamiltonian approach and reported the existence of stable periodic orbits in resonance with the Sun’s orbital motion in the Earth–Moon synodic frame. Using the well-known bi-circular problem (BCP) model and numerical continuation, for each TLP, Gómez et al. (2001b) reported the existence of two large stable and one small unstable periodic orbits whose existence is closely related to the Sun’s perturbation. They are also called dynamical substitutes which replace the role of equilibrium points in the perturbed dynamics. These forced periodic orbits are actually quasi-periodic in the real Earth–Moon system due to the quasi-periodic motion of the Moon and the existence of other perturbations aside from the Sun’s gravity. Also assuming a coplanar configuration of the Sun–Earth–Moon system, Díez, Jorba & Simó (1991) gave analytical expressions of the small unstable dynamical substitute by directly solving the equations of motion. Hou et al. (2010) obtained semi-analytic approximations of the small unstable dynamical substitute in the ephemeris model of the Earth–Moon system using numerical iterations. Using the same method, the authors also obtained semi-analytic approximations of the two large stable dynamical substitutes in the ephemeris model of the Earth–Moon system and studied the observation of possible natural objects residing on these orbits (Hou et al. 2015). In a recent study (Park et al. 2024), the authors have assessed the influences from the Moon’s orbit and the solar gravity in different regions in cislunar space.

The above studies focus mainly on the planar motion. For spatial orbital motions, Simó, Jorba & Masdemont (1995) reported the existence of spatial stable quasi-periodic orbits in the BCP model which appear as tori surrounding the two large stable dynamical substitutes. Focusing on the vertical dimension, Jorba (2000) tested the stability of these tori in the ephemeris model using the ephemeris provided by the Jet Propulsion Laboratory (JPL) and found that the vertical family 3 (VF3) has the best stability in the ephemeris Earth–Moon system. However, his study was restricted to less than 1000 yr due to the length limitation of the JPL ephemeris used, and he did not study the resonance mechanisms that determine the boundaries of the stable regions. By using a low-order approximation of the tori in the ephemeris model and by numerical integrations, Hou et al. (2015) studied both the planar and vertical stable regions around the two large dynamical substitutes. Lissauer et al. (2008) reported that all except a trivially small fraction of orbits near the Earth–Moon triangular Lagrangian points are unstable on geologically short time-scales using a more realistic model incorporating the eight planets, the Sun, the Earth, and the Moon. However, they did not provide mathematical explanations to the origin of stability or instability. Except for the studies listed above, there are also some numerical analyses on the stability of the TLP in the Earth–Moon system perturbed by the Sun. For example, see Slíz-Balogh et al. (2022) and some recent ones (Gimeno et al. 2024).

One reason that researchers are interested in the practical stable orbits around the TLPs in the real Earth–Moon system is that small natural bodies or dust may reside temporarily or permanently on orbits around these points (Freitas et al. 1980), so an investigation of the inherent dynamics is also important in the sense of astronomical observations. One example is the ghost Kordylewski dust cloud. Polish astronomer Kordylewski (1961) made the first photography observation of two faint patches near the L5 point from the Polish mountain Kasprowy Wierch between 1961 March 6 and April 6. These patches were positioned near the L5 point with an angular diameter of 6°. The Kordylewski dust cloud remains a mystery due to a lack of observational evidence in the past decades. Although it has been reported in many works (Slíz-Balogh, Barta & Horváth 2019; Wang, Jiang & Hou 2021), researchers believe that these dust clouds, if they exist, are clusters of dust particles whose structures are only temporary and constantly changing with time. Simulations were carried out by Slíz, Kovács & Süli (2017), Salnikova, Stepanov & Shuvalova (2018), and correspondence has been found comparing the simulation result (Slíz-Balogh, Barta & Horváth 2018) and the imaging polarimetric observation (Slíz-Balogh et al. 2019). However, the simulations are restricted to low inclination angles, and some simulations only last thousands of years. In the current study, in the Sun-perturbed Earth–Moon system given by the long JPL DE431 ephemeris, we analyse the spatial stability of motions around the TLPs. Our studies find: (1) practical stable orbits lasting for more than 15 000 yr exist in the ephemeris model, and the spatial stable regions are separated into three parts; (2) boundaries of the spatial stable regions are closely related to the resonances among the orbital precession rates of the small body and the Moon, the orbital frequency of the Sun, and the libration frequency of the co-orbital motion. We have managed to find these resonances in this work: (3) objects such as dust particles may reside on or accumulate close to these stable regions. It is found that there is a higher probability of finding objects (if they exist) at the region close to the TLP points, or at two regions with about |$25^\circ$| difference in declination with respect to the TLP points; (4) the solar radiation pressure (SRP) influences the spatial stable region. Generally, the stable region shrinks with the increase of the strength of the SRP, and no longer exists for small bodies with sizes of millimeters or smaller in magnitude; and (5) the spatial stable region of the BCP model shows obvious differences from that of the ephemeris model, indicating that the precession of the Moon’s orbit plays an important role in shaping the boundaries of the spatial stable region.

The rest of the paper is structured as follows. Section 2 introduces the models and methods of our study. Section 3 shows the results. Section 4 makes some discussions on the resonances in the ephemeris model, astronomical observations, the SRP’s influence, and the difference in the stable spatial region between the BCP model and the ephemeris model. Section 5 concludes the study.

2 MODELS AND METHODS

2.1 Ephemeris model description

To study the influence of the Sun’s gravitational perturbation on orbital motions in the Earth–Moon system, many models have been proposed, including the well-known BCP (Simó et al. 1995), the quasi-BCP (QBCP, Andreu 1998), and the restricted Hill four-body problem (Scheeres 1998). In this study, using the DE431 ephemeris provided by JPL (Folkner et al. 2014), we calculate the motion of the small body under gravities of the Sun, the Earth, and the Moon. To study the stability of specific orbits near the TLPs of the Earth–Moon system, a choice of an Earth-centred frame is reasonable. The equation of motion is given as follows (Hou et al. 2015):

$$\begin{eqnarray} \ddot{\boldsymbol{R}} = -\mu _E\frac{\boldsymbol{R}}{R^3} -\mu _M(\frac{\boldsymbol{R} - \boldsymbol{R}_M}{|\boldsymbol{R}-\boldsymbol{R}_M|^3} + \frac{\boldsymbol{R}_M}{R_M^3}) -\mu _S(\frac{\boldsymbol{R} - \boldsymbol{R}_S}{|\boldsymbol{R}-\boldsymbol{R}_S|^3} + \frac{\boldsymbol{R}_S}{R_S^3}),\\ \end{eqnarray}$$

(1)

where |$\boldsymbol{R}$| and |$\ddot{\boldsymbol{R}}$| represent the position and acceleration of the small body with respect to the Earth respectively. |$\boldsymbol{R}_M$| and |$\boldsymbol{R}_S$| are position vectors of the Sun and the Moon from the Earth obtained from ephemeris.

To visualize the motion around TLPs, the synodic frame is usually utilized. Denoting the position vectors in the Earth-centred sidereal frame [in our work, the geocentric celestial reference system (GCRS) frame is used] and the Earth-centred synodic frame as |$\boldsymbol{R}$| and |$\boldsymbol{r}$|⁠, the transformation between the Earth-centred sidereal frame and the Earth-centred synodic frame writes:

$$\begin{eqnarray} \boldsymbol{R} = k\bf {C}\boldsymbol{r}, \ \ \ \boldsymbol{r} =k^{-1} \bf {C}^T\boldsymbol{R} , \end{eqnarray}$$

(2)

where |$\bf {C}$| is the transformation matrix from the synodic frame to the sidereal frame and k is the instantaneous distance between the Earth and the Moon. Calculation of the matrix |$\bf {C}$| is easy and readers can refer to references for more details (Gómez et al. 2001a).

The following dimensionless units for length, mass, and time are used throughout our computation:

$$\begin{eqnarray} [L]=a_{\text{EM}}, \ [M]=m_{\text{E}} + m_{\text{M}}, \ [T]=\sqrt{[L]^3 / G[M]} \end{eqnarray}$$

in which |$a_{\text{EM}}$| is the mean semimajor axis of the lunar orbit, |$m_{\text{E}}$| and |$m_{\text{M}}$| are the mass of the Earth and the mass of the Moon, respectively, and G is the gravitational constant. In our work, |$a_{\text{EM}} = 3.84747981 \times 10^8 \ \mathrm{m}$|⁠, while |$m_{\text{E}}$| and |$m_{\text{M}}$| take values from the DE431 ephemeris (Folkner et al. 2014).

2.2 Vertical periodic family

In our work, we start from the vertical periodic family of the CRTBP model. Previous studies tell us that this family starts from infinitesimal vertical librations around the TLPs and terminates onto a vertical Lyapunov family of the collinear L3 point (Doedel, Romanov & Paffenroth 2007; Hou et al. 2008). For the Earth–Moon system, Fig. 1(a) shows the characteristic curve of the vertical periodic family around one TLP, along with that of the vertical periodic family around the point L3. Fig. 1(b) displays some example orbits in this family. The thick line in Fig. 1(a) indicates stable members of the vertical family.

Figure 1.

Part (a) shows characteristic curves (orbit period versus Jacobi constant) of the vertical periodic family around the TLP and the vertical periodic family around the point L3. ‘nd’ in the abscissa indicates the dimensionless units. The thick line begins at the L4 point indicates stable members of the vertical periodic family around the TLP. Part (b) displays some example orbits in the vertical periodic family.

Open in new tab Download slide

2.3 Numerical refinement

The second step of our study is to numerically continue periodic orbits of the vertical family to the ephemeris model. In previous works, a homotopy method is implemented to design periodic orbits in the elliptic restricted three-body problem from periodic orbits in the CRTBP model (Zheng & Zhao 2024). However, due to perturbations of different frequencies, the orbits are no longer periodic in the ephemeris model. As a result, we employ the multiple shooting method (Gómez, Masdemont & Simó 1998) to do the numerical continuation. The procedure is as follows:

$$\begin{eqnarray} t_0 = 0, \ t_1 = \Delta T, \ t_2 = 2\Delta T, \ \ldots , \ t_{N-1} = (N - 1)\Delta T. \end{eqnarray}$$

Map the above epochs into one orbital period of the periodic orbit by taking modulus with respect to the orbit period |$T_p$|⁠, i.e.

$$\begin{eqnarray} t_i^{\prime }= \mathrm{mod}(t_i, T_p), \ i=0, \ \ldots , N-1. \end{eqnarray}$$

By integrating the periodic orbit from |$t_0^{\prime } = 0$| to |$t_i^{\prime }$| in one single orbit period, we get a series of node points in the synodic frame:

$$\begin{eqnarray} \boldsymbol {X}_{\mathrm{ B}}(t_i^{\prime })= \begin{pmatrix}\boldsymbol{r}_i \\ \dot{\boldsymbol{r}}_i \end{pmatrix} , \ i=0, \ldots , N-1. \end{eqnarray}$$

Due to the periodicity, these node points are also the node points of the orbit at the epochs |$t_i$|⁠, i.e.

$$\begin{eqnarray} \boldsymbol {X}_{\mathrm{ B}}(t_i) = \boldsymbol {X}_{\mathrm{ B}}(t_i^{\prime }), \ i=0,\ldots ,N-1. \end{eqnarray}$$

We further transfer the node points from the barycentric synodic frame to the Earth-centred synodic frame by the following equation:

$$\begin{eqnarray} \boldsymbol {X}_{\mathrm{ E}}(t_i) = \boldsymbol {X}_{\mathrm{ B}}(t_i) + (\mu , 0, 0, 0, 0, 0)^T. \end{eqnarray}$$

(2) At each epoch |$t_i$|⁠, we get the transformation matrix |$\bf {C}_i$| from the synodic frame to the sidereal frame using the Moon’s orbit ephemeris. Using the relation given by equation (2), we transform the series of node points from the Earth-centred synodic frame to the GCRS frame.

(3) Using the node points in the GCRS frame and the equations of motion by equation (1), we employ the multiple shooting algorithm described in detail by Gómez et al. (1998) to numerically refine the node points to get a continuous orbit in the ephemeris model within a prescribed precision.

(4) After the continuous orbit in the GCRS frame is obtained, using the inverse transformation of equation (2), we get the position and speed of the orbit in the barycentric synodic frame.

In our work, the time interval is set as |$4 \ \mathrm{d}$|⁠, and the total number of node points is |$N=457$|⁠, which results in a numerically refined orbit of 5 yr. Fig. 2 shows two example orbits. One is practically stable (Fig. 2a) and one is unstable (Fig. 2b). We will come back to these two example orbits in Section 3.1. For readers to check the results, initial conditions of some practically stable example orbits are given in Table 1. The initial epoch for these orbits is J2000.0.

Figure 2.

Part (a) shows a practically stable orbit with a vertical amplitude of 0.547148. Part (b) shows an unstable orbit with a vertical amplitude of 0.597987. Both orbits are integrated for 2000 yr and it can be seen that the stable one keeps its original track, while the unstable one has a tendency to get away from its initial trajectory.

Open in new tab Download slide

Table 1.

Open in new tab

Initial conditions of some practically stable orbits around the TLP L5 in the GCRS frame. The initial epoch is J2000.0.

No.	\|$x \ \mathrm{(km)}$\|	\|$y \ \mathrm{(km)}$\|	\|$z \ \mathrm{(km)}$\|	\|$\dot{x} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{y} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{z} \ \mathrm{(km\,s^{-1})}$\|
1	\|$-373636.68606$\|	113900.29464	\|$-43970.01196$\|	\|$-0.330516340$\|	\|$-0.940562107$\|	0.0092540216
2	\|$-373756.25070$\|	112860.36004	\|$-45772.10566$\|	\|$-0.329154942$\|	\|$-0.940951838$\|	0.0167955237
3	\|$-373889.80225$\|	111487.20017	\|$-48298.84612$\|	\|$-0.327528490$\|	\|$-0.941278717$\|	0.0277596829
4	\|$-373888.37926$\|	107292.96127	\|$-57878.41273$\|	\|$-0.324618266$\|	\|$-0.939856389$\|	0.0738932978
5	\|$-374281.07285$\|	103184.71851	\|$-64289.96423$\|	\|$-0.320654488$\|	\|$-0.936923736$\|	0.1116355423
6	\|$-374316.45613$\|	103003.08692	\|$-64485.76785$\|	\|$-0.320427432$\|	\|$-0.936805596$\|	0.1129496806
7	\|$-377960.93185$\|	76908.840559	\|$-76258.29539$\|	\|$-0.292751626$\|	\|$-0.913497015$\|	0.2809561778
8	\|$-378493.35582$\|	73599.131859	\|$-76281.13755$\|	\|$-0.288173041$\|	\|$-0.907551961$\|	0.3067527971
9	\|$-381841.02136$\|	68560.410238	\|$-76160.61576$\|	\|$-0.282940263$\|	\|$-0.870909104$\|	0.3916785607
10	\|$-386421.29393$\|	54534.001480	\|$-68787.69384$\|	\|$-0.261471728$\|	\|$-0.842851726$\|	0.4595124446

No.	\|$x \ \mathrm{(km)}$\|	\|$y \ \mathrm{(km)}$\|	\|$z \ \mathrm{(km)}$\|	\|$\dot{x} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{y} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{z} \ \mathrm{(km\,s^{-1})}$\|
1	\|$-373636.68606$\|	113900.29464	\|$-43970.01196$\|	\|$-0.330516340$\|	\|$-0.940562107$\|	0.0092540216
2	\|$-373756.25070$\|	112860.36004	\|$-45772.10566$\|	\|$-0.329154942$\|	\|$-0.940951838$\|	0.0167955237
3	\|$-373889.80225$\|	111487.20017	\|$-48298.84612$\|	\|$-0.327528490$\|	\|$-0.941278717$\|	0.0277596829
4	\|$-373888.37926$\|	107292.96127	\|$-57878.41273$\|	\|$-0.324618266$\|	\|$-0.939856389$\|	0.0738932978
5	\|$-374281.07285$\|	103184.71851	\|$-64289.96423$\|	\|$-0.320654488$\|	\|$-0.936923736$\|	0.1116355423
6	\|$-374316.45613$\|	103003.08692	\|$-64485.76785$\|	\|$-0.320427432$\|	\|$-0.936805596$\|	0.1129496806
7	\|$-377960.93185$\|	76908.840559	\|$-76258.29539$\|	\|$-0.292751626$\|	\|$-0.913497015$\|	0.2809561778
8	\|$-378493.35582$\|	73599.131859	\|$-76281.13755$\|	\|$-0.288173041$\|	\|$-0.907551961$\|	0.3067527971
9	\|$-381841.02136$\|	68560.410238	\|$-76160.61576$\|	\|$-0.282940263$\|	\|$-0.870909104$\|	0.3916785607
10	\|$-386421.29393$\|	54534.001480	\|$-68787.69384$\|	\|$-0.261471728$\|	\|$-0.842851726$\|	0.4595124446

Table 1.

Open in new tab

Initial conditions of some practically stable orbits around the TLP L5 in the GCRS frame. The initial epoch is J2000.0.

No.	\|$x \ \mathrm{(km)}$\|	\|$y \ \mathrm{(km)}$\|	\|$z \ \mathrm{(km)}$\|	\|$\dot{x} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{y} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{z} \ \mathrm{(km\,s^{-1})}$\|
1	\|$-373636.68606$\|	113900.29464	\|$-43970.01196$\|	\|$-0.330516340$\|	\|$-0.940562107$\|	0.0092540216
2	\|$-373756.25070$\|	112860.36004	\|$-45772.10566$\|	\|$-0.329154942$\|	\|$-0.940951838$\|	0.0167955237
3	\|$-373889.80225$\|	111487.20017	\|$-48298.84612$\|	\|$-0.327528490$\|	\|$-0.941278717$\|	0.0277596829
4	\|$-373888.37926$\|	107292.96127	\|$-57878.41273$\|	\|$-0.324618266$\|	\|$-0.939856389$\|	0.0738932978
5	\|$-374281.07285$\|	103184.71851	\|$-64289.96423$\|	\|$-0.320654488$\|	\|$-0.936923736$\|	0.1116355423
6	\|$-374316.45613$\|	103003.08692	\|$-64485.76785$\|	\|$-0.320427432$\|	\|$-0.936805596$\|	0.1129496806
7	\|$-377960.93185$\|	76908.840559	\|$-76258.29539$\|	\|$-0.292751626$\|	\|$-0.913497015$\|	0.2809561778
8	\|$-378493.35582$\|	73599.131859	\|$-76281.13755$\|	\|$-0.288173041$\|	\|$-0.907551961$\|	0.3067527971
9	\|$-381841.02136$\|	68560.410238	\|$-76160.61576$\|	\|$-0.282940263$\|	\|$-0.870909104$\|	0.3916785607
10	\|$-386421.29393$\|	54534.001480	\|$-68787.69384$\|	\|$-0.261471728$\|	\|$-0.842851726$\|	0.4595124446

No.	\|$x \ \mathrm{(km)}$\|	\|$y \ \mathrm{(km)}$\|	\|$z \ \mathrm{(km)}$\|	\|$\dot{x} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{y} \ \mathrm{(km\,s^{-1})}$\|	\|$\dot{z} \ \mathrm{(km\,s^{-1})}$\|
1	\|$-373636.68606$\|	113900.29464	\|$-43970.01196$\|	\|$-0.330516340$\|	\|$-0.940562107$\|	0.0092540216
2	\|$-373756.25070$\|	112860.36004	\|$-45772.10566$\|	\|$-0.329154942$\|	\|$-0.940951838$\|	0.0167955237
3	\|$-373889.80225$\|	111487.20017	\|$-48298.84612$\|	\|$-0.327528490$\|	\|$-0.941278717$\|	0.0277596829
4	\|$-373888.37926$\|	107292.96127	\|$-57878.41273$\|	\|$-0.324618266$\|	\|$-0.939856389$\|	0.0738932978
5	\|$-374281.07285$\|	103184.71851	\|$-64289.96423$\|	\|$-0.320654488$\|	\|$-0.936923736$\|	0.1116355423
6	\|$-374316.45613$\|	103003.08692	\|$-64485.76785$\|	\|$-0.320427432$\|	\|$-0.936805596$\|	0.1129496806
7	\|$-377960.93185$\|	76908.840559	\|$-76258.29539$\|	\|$-0.292751626$\|	\|$-0.913497015$\|	0.2809561778
8	\|$-378493.35582$\|	73599.131859	\|$-76281.13755$\|	\|$-0.288173041$\|	\|$-0.907551961$\|	0.3067527971
9	\|$-381841.02136$\|	68560.410238	\|$-76160.61576$\|	\|$-0.282940263$\|	\|$-0.870909104$\|	0.3916785607
10	\|$-386421.29393$\|	54534.001480	\|$-68787.69384$\|	\|$-0.261471728$\|	\|$-0.842851726$\|	0.4595124446

2.4 Frequency analysis method

As we have mentioned, one goal of the current study is to find the resonances that determine the boundaries of the spatial stable region. Since we lack analytical solutions in the current study, the basic frequencies are obtained by numerical methods. To analyse the orbit data in the frequency domain, a standard fast Fourier transform (FFT) is implemented, i.e.

$$\begin{eqnarray} x(t) \xrightarrow {\text{FFT}} X(\omega). \end{eqnarray}$$

However, the FFT method has some drawbacks if implemented directly. Due to aliasing and the leakage effect, some frequencies with small amplitudes may be submerged in the sidelobes of the strong ones. As a result, it becomes necessary to extract each frequency individually and to perform refinements on them. Some methods are proposed using multiple iteration processes to refine the frequencies obtained (Laskar 1993; Gómez, Mondelo & Simó 2010). The method we use is more direct and easier to implement. Denoting the original signal as |$s^{(0)}(t)$| and acquiring |$N_s$| sample points with a sampling interval |$\Delta t$| at a series of points |$s^{(0)}(t_i)$| where |$i=1,\ldots , N_s$|⁠, we are able to extract all frequencies individually from the original signal following the procedures below:

(1) As a start, we obtain the frequency series from the time signal series using FFT. Denote the frequencies obtained as:

$$\begin{eqnarray} \omega _j = j\omega _r, \ \ \ \omega _r = \frac{2\pi }{N_s\Delta t}, \end{eqnarray}$$

where |$\omega _r$| is the resolution of the FFT frequency map which depends on the total number of signal points and the sampling time interval |$\Delta t$|⁠. These frequencies are discrete and only approximations of accurate frequencies in the signal. In our work, |$N_s=2^{22}$| and |$\Delta t=0.5 \ \mathrm{d}$|⁠.

(2) Starting from the frequency with the largest amplitude, we obtain accurate values of the selected frequency and its associated coefficients by continuous Fourier transformation (CFT). For example, we select one frequency |$\omega _j$| produced by FFT and we perform CFT using the following formulae for an arbitrary frequency |$\omega$| in the range |$[\omega _j - \omega _r, \omega _j + \omega _r]$|⁠:

$$\begin{eqnarray} C^{(0)}(\omega) &=\frac{2}{T}\int _0^{T} x^{(0)}(t)\cos {\omega t} \mathrm{d}t, \\ S^{(0)}(\omega) &=\frac{2}{T}\int _0^{T} x^{(0)}(t)\sin {\omega t} \mathrm{d}t. \end{eqnarray}$$

The amplitude of the coefficient is denoted as:

$$\begin{eqnarray} A(\omega) = \sqrt{[C^{(0)}(\omega)]^2+[S^{(0)}(\omega)]^2}. \end{eqnarray}$$

It is an obvious fact that the accurate value of the frequency |$\omega _j$| corresponds to the value that makes the amplitude maximum, i.e.

$$\begin{eqnarray} (\omega _j)_{\mathrm{accurate}} = \lbrace \omega | A(\omega) = \max A(\omega ^{\prime }), \omega ^{\prime } \in [\omega _j - \omega _r, \omega _j + \omega r] \rbrace . \end{eqnarray}$$

Once we get the accurate value of the frequency, we also obtain its accurate coefficients.

(3) Suppose we have obtained |$N_\omega$| frequencies and their associated coefficients. We repeat the above step (2) for each frequency for refinement, but the time signal series is updated as:

$$\begin{eqnarray} \begin{aligned} x_n^{(1)} (t_i) = &x^{(0)}(t_i) \\ &- \sum _{k=0, k\ne n}^{N_\omega - 1} [C^{(0)}(\omega _k) \cos {(\omega _k t_i)} + C^{(0)}(\omega _k) \sin {(\omega _k t_i)}] \end{aligned} \end{eqnarray}$$

in which n is the index of the selected frequency to be refined. The reason for us doing so is to remove as much the influence from other frequencies as possible when we look for refined estimates of the selected frequency and its associated coefficients.

(4) Repeat steps (2) and (3) until the selected frequencies no longer change within a prescribed accuracy. In our work, the value of |$N_\omega$| changes for different orbits. The criterion for this number is decided by a minimum amplitude threshold in our study, which is chosen to be |$A_{\mathrm{min}} = 5 \times 10^{-5}$|⁠, and the prescribed accuracy is |$\Delta (\omega _j)_{\mathrm{accurate}} \lt 10^{-5} \omega _r$|⁠.

Assuming that the Moon’s orbit is quasi-periodic in the Earth–Moon–Sun system, we can extract five basic frequencies (Gómez et al. 2001b) by analysing the Moon’s motion in the GCRS frame. These five frequencies are recovered by applying the above algorithm to the lunar orbit obtained from the JPL DE431 ephemeris. Using the dimensionless units described in Section 2.1, these frequencies are:

(1) the frequency of the mean longitude of the Moon |$\lambda _\mathrm{ M}$|⁠, |$n_1 = 1.0$|⁠;

(2) the frequency of the mean longitude of the Moon from the Sun |$\lambda _\mathrm{ M} - \lambda _\mathrm{ S}$|⁠, |$n_2 = 0.925198240139$|⁠;

(3) the frequency of the mean longitude of the lunar perigee |$\tilde{\omega }_\mathrm{ M}$|⁠, |$n_3 = 0.0084517202$|⁠;

(4) the frequency of the mean longitude of the ascending node of the lunar orbit on the ecliptic |$\Omega _\mathrm{ M}$|⁠, |$n_4 = -0.004021929$|⁠; and

(5) the frequency of the mean longitude of the solar perigee |$\tilde{\omega }_\mathrm{ S}$|⁠, |$n_5 = 0.00000097604$|⁠.

Because of the azimuthal invariance in the third-body perturbation expansion (Murray et al. 1999), there exists only four instead of five independent frequencies (Li B. et al. 2023). In the lunar theory (Deprit et al. 1967), four combinations are often used and they are denoted as |$\boldsymbol{W} = (\omega _1, \omega _2, \omega _3, \omega _4)$|⁠, whose values are recovered as:

$$\begin{eqnarray} \omega _1 &=& n_1 - n_3 = 0.99154827976 \\ \omega _2 &=& n_1 - n_2 - n_5 = 0.07480135273 \\ \omega _3 &=& n_2 = 0.925198240139 \\ \omega _4 &=& n_1 - n_4 = 1.00402079132 \end{eqnarray}$$

in our study.

2.5 The BCP model

In the discussion section, we will address the difference between the ephemeris model and the BCP model. We believe readers can find details of this well-known model in the literatures, so we only briefly introduce the BCP model here. In the BCP model, the Moon’s orbit with respect to the Earth is assumed as a circular orbit following the Keplerian law, and the barycentre of the Earth and the Moon is moving on a Keplerian circular orbit around the Sun. The two orbits are assumed coplanar. An intuitive figure of this model is shown in Fig. 3.

Figure 3.

The BCP model of the Earth–Moon system perturbed by the Sun.

Open in new tab Download slide

Using the following dimensionless units commonly used in the CRTBP model (Szebehely 1967)

$$\begin{eqnarray} [L]=a, \ \ [M] = m_\mathrm{ E} + m_\mathrm{ M}, \ \ [T] = \sqrt{[L]^3 / G[M]} \end{eqnarray}$$

(3)

in which a is the semimajor axis of the lunar orbit, |$m_\mathrm{ E}$| and |$m_\mathrm{ M}$| are masses of the Earth and the Moon, and G is the gravitational constant, equations of motion in the Earth–Moon barycentric synodic frame are (Simó et al. 1995):

$$\begin{eqnarray} &\ddot{\boldsymbol{r}} + 2 \begin{pmatrix}-\dot{y} \\ \dot{x} \\ 0 \end{pmatrix} = \frac{\partial \Omega }{\partial \boldsymbol{r}}, \\ &\Omega = \frac{\mu _\mathrm{ E}}{r_1} + \frac{\mu _\mathrm{ M}}{r_2} + \mu _\mathrm{ S} (\frac{1}{|\boldsymbol{r}-\boldsymbol{r}_\mathrm{ S}|} - \frac{\boldsymbol{r}\cdot \boldsymbol{r}_\mathrm{ S}}{r_\mathrm{ S}^3}) \end{eqnarray}$$

in which:

$$\begin{eqnarray} &\mu _\mathrm{ E}=\frac{m_\mathrm{ E}}{m_\mathrm{ E} + m_\mathrm{ M}}, \ \ \mu _\mathrm{ M} = \frac{m_\mathrm{ M}}{m_\mathrm{ E} + m_\mathrm{ M}}, \ \ \mu _\mathrm{ S} = \frac{m_\mathrm{ S}}{m_\mathrm{ E} + m_\mathrm{ M}}, \\ &\boldsymbol{r}_\mathrm{ S} = (a_\mathrm{ S} \cos {\theta _\mathrm{ S}}, a_\mathrm{ S}\sin {\theta _\mathrm{ S}}, 0)^T. \end{eqnarray}$$

|$m_\mathrm{ S}$| is the mass of the Sun and |$a_\mathrm{ S}$| is the semimajor axis of the Sun’s orbit. Since we assume a coplanar configuration and a circular orbit of the Earth–Moon barycentre with respect to the Sun, the angle |$\theta _\mathrm{ S}$| actually retrogrades in the Earth–Moon synodic frame. Using the dimensionless units in equation (3), we have:

$$\begin{eqnarray} \dot{\theta }_\mathrm{ S}=\sqrt{\frac{1+\mu _\mathrm{ S}}{a_\mathrm{ S}^3}} - 1. \end{eqnarray}$$

We follow the same procedure as described in Section 2.3 to get quasi-periodic orbits in the BCP model. By nature, these quasi-periodic orbits can be viewed as numerical approximations of the ‘VF3’ family (Jorba 2000).

3 RESULTS

3.1 Two orbit types

Substitute the transformation equation (2) into equation (1), and we obtain the equations of motion in the Earth-centred synodic frame (⁠|$\bf {I}_3$| represents the identity matrix):

$$\begin{eqnarray} \ddot{\boldsymbol{r}} &=&-2(\frac{\dot{k}}{k} \bf {I}_3 + \bf {C}^T\dot{\bf {C}}) \dot{\boldsymbol{r}} -(2\frac{\dot{k}}{k}\bf {C}^T\dot{\bf {C}} + \bf {C}^T\ddot{\bf {C}} + \frac{\ddot{k}}{k}\bf {I}_3) \boldsymbol{r} \\ && -\frac{\mu _E}{k^3} \frac{\boldsymbol{r}}{r^3} - \frac{\mu _\mathrm{ M}}{k^3}(\frac{\boldsymbol{r} - \boldsymbol{r}_\mathrm{ M}}{|\boldsymbol{r} - \boldsymbol{r}_\mathrm{ M}|^3} + \frac{\boldsymbol{r}_\mathrm{ M}}{r_\mathrm{ M}^3}) - \frac{\mu _\mathrm{ S}}{k^3}(\frac{\boldsymbol{r} - \boldsymbol{r}_\mathrm{ S}}{|\boldsymbol{r} - \boldsymbol{r}_\mathrm{ S}|^3} + \frac{\boldsymbol{r}_\mathrm{ S}}{r_\mathrm{ S}^3}). \end{eqnarray}$$

(4)

Denoting the position vector of the TLP as |$\boldsymbol{r}_0$|⁠, the position vector from the TLP as |$\boldsymbol{\rho }$|⁠, we have |$\boldsymbol{r} = \boldsymbol{r}_0 + \boldsymbol{\rho }$|⁠. Substituting this relation into equation (4) by grouping different terms, we reach the following equations of motion (Hou et al. 2010, 2011):

$$\begin{eqnarray} \ddot{\boldsymbol{\rho }} = \boldsymbol{F}_1(\boldsymbol{\rho }, \dot{\boldsymbol{\rho }}, t) + \boldsymbol{F}_2(t) , \end{eqnarray}$$

(5)

where:

$$\begin{eqnarray} \boldsymbol{F}_1(\boldsymbol{\rho }, \dot{\boldsymbol{\rho }}, t) &=&-2(\frac{\dot{k}}{k} \bf {I}_3 + \bf {C}^T\dot{\bf {C}}) \dot{\boldsymbol{r}} - (2\frac{\dot{k}}{k}\bf {C}^T\dot{\bf {C}} + \bf {C}^T\ddot{\bf {C}} + \frac{\ddot{k}}{k}\bf {I}_3)(\boldsymbol{r}-\boldsymbol{r}_0) \\ && -\frac{\mu _\mathrm{ E}}{k^3} (\frac{\boldsymbol{r}}{r^3}-\frac{\boldsymbol{r}_0}{r_0^3}) \\ && -\frac{\mu _\mathrm{ M}}{k^3}(\frac{\boldsymbol{r} - \boldsymbol{r}_\mathrm{ M}}{|\boldsymbol{r} - \boldsymbol{r}_\mathrm{ M}|^3} - \frac{\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ M}}{|\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ M}|^3}) \\ && -\frac{\mu _\mathrm{ S}}{k^3}(\frac{\boldsymbol{r} - \boldsymbol{r}_\mathrm{ S}}{|\boldsymbol{r} - \boldsymbol{r}_\mathrm{ S}|^3} - \frac{\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ S}}{|\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ S}|^3}), \end{eqnarray}$$

$$\begin{eqnarray} \boldsymbol{F}_2(t)&=&-\frac{\mu _E}{k^3}\frac{\boldsymbol{r}_0}{r_0^3} \\ && -\frac{\mu _\mathrm{ M}}{k^3}(\frac{\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ M}}{|\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ M}|^3} + \frac{\boldsymbol{r}_\mathrm{ M}}{r_\mathrm{ M}^3}) \\ && -\frac{\mu _\mathrm{ S}}{k^3}(\frac{\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ S}}{|\boldsymbol{r}_0 - \boldsymbol{r}_\mathrm{ S}|^3} + \frac{\boldsymbol{r}_\mathrm{ S}}{r_\mathrm{ S}^3}). \end{eqnarray}$$

The second term at the right-hand of equation (5) is independent of the state of the small body, indicating that the TLPs are no longer equilibrium points. The forced motions caused by the second term are called dynamical substitutes. As already mentioned in the Introduction, for each TLP in the ephemeris Earth–Moon system, currently three such dynamical substitutes are reported. One is unstable and the other two are stable (Hou et al. 2010, 2015).

Denote the state vector of an orbit and the first-order derivative as:

$$\begin{eqnarray} \boldsymbol {X} = \begin{pmatrix}\boldsymbol{\rho } \\ \dot{\boldsymbol{\rho }} \end{pmatrix}, \ \ \boldsymbol {\mathcal {F}}_1= \begin{pmatrix}\dot{\boldsymbol{\rho }} \\ \boldsymbol{F}_1 \end{pmatrix}, \ \ \boldsymbol {\mathcal {F}}_2 = \begin{pmatrix}\boldsymbol {0} \\ \boldsymbol{F}_2 \end{pmatrix} \end{eqnarray}$$

and we can rewrite equation (5) as:

$$\begin{eqnarray} \dot{\boldsymbol {X}} = \boldsymbol {\mathcal {F}}_1(\boldsymbol {X}, t) + \boldsymbol {\mathcal {F}}_2(t). \end{eqnarray}$$

(6)

Denoting the dynamical substitute as |$\bar{\boldsymbol {X}}$| and an orbit around it as |$\bar{\boldsymbol {X}} + \Delta \boldsymbol {X}$|⁠, we can expand equation (6) around the dynamical substitute and retain the linear term:

$$\begin{eqnarray} \Delta \dot{\boldsymbol {X}} = \boldsymbol {A}\Delta \boldsymbol {X}, \ \ \boldsymbol {A} = \frac{\partial \boldsymbol {F}_1}{\partial \boldsymbol {X}} |_{\boldsymbol {X} = \bar{\boldsymbol {X}}}. \end{eqnarray}$$

(7)

Under the assumption that orbits of the Moon and the Sun are quasi-periodic, the Jacobian matrix |$\bf {A}$| in equation (7) is quasi-periodic, with four basic frequencies as reported in Section 2.4. Introducing a transformation |$\Delta \boldsymbol {X} = \boldsymbol {B}\boldsymbol {Y}$| in which |$\boldsymbol {B}$| is also quasi-periodic with the same basic frequencies as the matrix |$\bf {A}$|⁠, equation (7) is transformed to:

$$\begin{eqnarray} \dot{\boldsymbol {Y}} = \boldsymbol {E}\boldsymbol {Y}, \ \ \boldsymbol {E} = \boldsymbol {B}^{-1}(\boldsymbol {A}\boldsymbol {B} - \dot{\boldsymbol {B}}). \end{eqnarray}$$

(8)

According to the quasi-Floquet theorem (Jorba et al. 1996), a suitable choice of the matrix |$\bf {B}$| can make the matrix |$\bf {E}$| in equation (8) a constant matrix. Obviously, eigenvalues of the matrix |$\bf {E}$| determine the stability of the orbits around the dynamical substitute. Previous studies show that eigenvalues of the matrix |$\bf {E}$| for the small unstable dynamical substitute are of the following forms (Hou et al. 2010):

$$\begin{eqnarray} \pm \lambda , \ \ \pm i\omega _l, \ \ \pm i\omega _v \end{eqnarray}$$

and eigenvalues for the two large stable dynamical substitutes are of the following forms:

$$\begin{eqnarray} \pm i\omega _s, \ \ \pm i\omega _l, \ \ \pm i\omega _v. \end{eqnarray}$$

The subscripts l, s, and v indicate the long-period, the short-period, and the vertical components. Solutions to equation (8) are of the following forms:

$$\begin{eqnarray} \boldsymbol {Y} = (\alpha _1 e^{\lambda t}, \alpha _2 e^{-\lambda t}, \beta \cos {\theta _l}, \beta \sin {\theta _l}, \gamma \cos {\theta _v}, \gamma \sin {\theta _v}) \end{eqnarray}$$

for the unstable case, and of the following forms:

$$\begin{eqnarray} \boldsymbol {Y} = (\alpha \cos {\theta _l}, \alpha \sin {\theta _l}, \beta \cos {\theta _l}, \beta \sin {\theta _l}, \gamma \cos {\theta _v}, \gamma \sin {\theta _v}) \end{eqnarray}$$

for the stable case. In the above equations,

$$\begin{eqnarray} \theta _l = \omega _l t + \theta _{l0}, \ \ \theta _s = \omega _s t + \theta _{s0}, \ \ \theta _v = \omega _v t + \theta _{v0}. \end{eqnarray}$$

The angles |$\theta _{l0}$|⁠, |$\theta _{s0}$|⁠, and |$\theta _{v0}$| are integration constants which can be arbitrarily chosen. The amplitude parameters |$\alpha _1$|⁠, |$\alpha _2$|⁠, |$\beta$|⁠, and |$\gamma$| for the unstable case (⁠|$\alpha$|⁠, |$\beta$|⁠, and |$\gamma$| for the stable case) are also integration constants that can be arbitrarily chosen. Transforming back to the original coordinates, if we only focus on quasi-periodic motions, there are six frequencies (⁠|$\omega _1$|⁠, |$\omega _2$|⁠, |$\omega _3$|⁠, |$\omega _4$| plus two free frequencies |$\omega _l$|⁠, |$\omega _v$|⁠) in total for the unstable case and seven frequencies (⁠|$\omega _1$|⁠, |$\omega _2$|⁠, |$\omega _3$|⁠, |$\omega _4$| plus three free frequencies |$\omega _l$|⁠, |$\omega _s$|⁠, |$\omega _v$|⁠) for the stable case. As a result, the numerical tori constructed in Section 2.3 are of two types. One type which is unstable can be expressed as:

$$\begin{eqnarray} \left\lbrace \begin{aligned} x = \sum [C^{\lt x\gt }_{ijklmn} \cos {\theta _{ijklmn}} + S^{\lt x\gt }_{ijklmn} \sin {\theta _{ijklmn}}]\\ y = \sum [C^{\lt y\gt }_{ijklmn} \cos {\theta _{ijklmn}} + S^{\lt y\gt }_{ijklmn} \sin {\theta _{ijklmn}}]\\ z = \sum [C^{\lt z\gt }_{ijklmn} \cos {\theta _{ijklmn}} + S^{\lt z\gt }_{ijklmn} \sin {\theta _{ijklmn}}], \end{aligned} \right. \end{eqnarray}$$

(9)

where |$\theta _{ijklmn} = (i\omega _1 + j\omega _2 + k\omega _3 + l\omega _4)t + m\theta _l + n\theta _v$| and |$C^{\lt x\gt }_{ijklmn}$|⁠, |$S^{\lt x\gt }_{ijklmn}$|⁠, |$C^{\lt y\gt }_{ijklmn}$|⁠, |$S^{\lt y\gt }_{ijklmn}$|⁠, |$C^{\lt z\gt }_{ijklmn}$|⁠, and |$S^{\lt z\gt }_{ijklmn}$| are respective amplitude coefficients.

The other type which is stable can be expressed as:

$$\begin{eqnarray} \left\lbrace \begin{aligned} x = \sum [C^{\lt x\gt }_{ijklmnp} \cos {\theta _{ijklmnp}} + S^{\lt x\gt }_{ijklmnp} \sin {\theta _{ijklmnp}}]\\ y = \sum [C^{\lt y\gt }_{ijklmnp} \cos {\theta _{ijklmnp}} + S^{\lt y\gt }_{ijklmnp} \sin {\theta _{ijklmnp}}]\\ z = \sum [C^{\lt z\gt }_{ijklmnp} \cos {\theta _{ijklmnp}} + S^{\lt z\gt }_{ijklmnp} \sin {\theta _{ijklmnp}}] ,\end{aligned} \right. \end{eqnarray}$$

(10)

where |$\theta _{ijklmnp} = (i\omega _1 + j\omega _2 + k\omega _3 + l\omega _4)t + m\theta _l + n\theta _v + p\theta _s$|⁠. Considering motions far away from the dynamical substitutes, the free frequencies are no longer the eigenvalues of the matrix |$\bf {E}$| in equation (8) due to non-linear effects. According to previous studies (Hou et al. 2010, 2015), they should be expanded as:

$$\begin{eqnarray} \omega _l = \sum _{j(2)=0}\sum _{k(2)=0}\omega _{jk}^{\lt l\gt }\beta ^j\gamma ^k, \ \ \omega _v = \sum _{j(2)=0}\sum _{k(2)=0}\omega _{jk}^{\lt v\gt }\beta ^j\gamma ^k \end{eqnarray}$$

for the unstable case, and

$$\begin{eqnarray} \omega _s &=& \sum _{i(2)=0}\sum _{j(2)=0}\sum _{k(2)=0}\omega _{ijk}^{\lt l\gt }\alpha ^i\beta ^j\gamma ^k \\ \omega _l &=& \sum _{i(2)=0}\sum _{j(2)=0}\sum _{k(2)=0}\omega _{ijk}^{\lt l\gt }\alpha ^i\beta ^j\gamma ^k \\ \omega _v &=& \sum _{i(2)=0}\sum _{j(2)=0}\sum _{k(2)=0}\omega _{ijk}^{\lt v\gt }\alpha ^i\beta ^j\gamma ^k \end{eqnarray}$$

for the stable case. The ‘(2)’ symbol means that the numeration in the summation is even integer.

In our work, we mainly focus on the spatial stable regions around the TLP. That is to say, we investigate the mechanism through which the vertical component |$\gamma$| influences the stability of the orbit and neglect the contribution from the long- and the short-period components, i.e. |$\alpha =\beta =0$|⁠. In such a case, we have unified expressions for the three frequencies expressed as:

$$\begin{eqnarray} \omega _s &=& \sum _{k(2)=0}\omega _{00k}^{\lt s\gt } \gamma ^k \\ \omega _l &=& \sum _{k(2)=0}\omega _{00k}^{\lt l\gt } \gamma ^k \\ \omega _v &=& \sum _{k(2)=0}\omega _{00k}^{\lt v\gt } \gamma ^k \end{eqnarray}$$

(11)

In the frequency domain, the two different types of orbits show different frequency patterns. Parts of the frequency maps of the two orbits in Fig. 2 are displayed in Fig. 4. Fig. 4(a) is for the stable orbit and Fig. 4(b) is for the unstable orbit. Since we are using numerical approximations of the unstable tori, the hyperbolic component is not completely zero in our case. This leads to chaos in the frequency domain for unstable orbits, i.e. the frequencies are not isolated but appear as a bunch of mixed frequencies. In other words, the frequency map of the stable orbit is clean in the sense that the frequencies are isolated. The FFT results are calculated using a data number of |$2^{21}$| with a sampling interval of |$0.5 \ \mathrm{d}$|⁠. The vertical dashed lines in Fig. 4(a) are important frequencies that determine the boundaries of the spatial stable region. They will be addressed later.

Figure 4.

Part (a) displays part of the FFT map of the stable orbit shown in Figs 2(a) and (b) displays the same part of the frequency domain for the orbit shown in Fig. 2(b). The FFT map in part (a) indicates a stable orbit in the sense that the frequencies are isolated, and the FFT map in part (b) is considered to be a pattern of an unstable orbit in the sense that the frequencies appear as a bunch of mixed frequencies.

Open in new tab Download slide

3.2 Spatial stable region

Here, ‘spatial’ means that we mainly focus on the vertical component whose amplitude is given by the value of |$\gamma$|⁠. This value corresponds to the coefficient of the free frequency |$\omega _v$|⁠, which can be obtained from frequency analysis on the z coordinate of the quasi-periodic orbits. Theoretically, we can separate the stable orbits from the unstable ones by carefully analysing their frequency maps, because the two types of ‘quasi-periodic orbits’ have different patterns in the frequency domain. Nevertheless, in order to do so we need a long enough orbit because some orbits are only weakly unstable. This is usually difficult and costs a lot of time. As a result, we take a simple and intuitive approach. We directly integrate the quasi-periodic orbits obtained in Section 2.3 in the ephemeris model for a long time and observe whether it can survive after the integration. The initial epoch is set as J2000.0, and the ephemeris used is the DE431 ephemeris. Due to the length limitation of the ephemeris, the maximum integration time is 15 000 yr. We think the orbit is stable if it can survive after 15 000 yr. The criterion of ‘survive’ is that the x coordinate in the synodic frame is no smaller than |$-0.2$|⁠. Here is the procedure for us to search for the spatial stable regions:

(2) For each vertical periodic orbit |$\boldsymbol {P}_i$| in the CRTBP model, use the numerical algorithm described in Section 2.3 to obtain a 5-yr length quasi-periodic orbit |$\boldsymbol {QP}_i$| in the ephemeris model. For each |$\boldsymbol {QP}_i$| in the ephemeris model, get its value of |$\gamma$| by frequency analysis on its z-coordinate data. Generally, the value of |$\gamma$| is very close to the one of the periodic orbit |$\boldsymbol {P}_i$|⁠.

(3) Starting from the initial point of each quasi-periodic orbit |$\boldsymbol {QP}_i$|⁠, we integrate it in the ephemeris model within an integration time of |$15\,000 \ \mathrm{yr}$|⁠. We record the time |$\boldsymbol {T}_f$| when the orbit becomes unstable, i.e. |$x \lt -0.2$| in the synodic frame. We think the orbit is stable if |$\boldsymbol {T}_f = 15\,000 \ \mathrm{yr}$|⁠. In this way, we are able to describe the spatial stable region in the ephemeris model.

Fig. 5 shows the result. The lower abscissa is the |$\gamma$|-value, and the ordinate is the time |$\boldsymbol {T}_f$|⁠. There are three practical stable regions in the figure, indicated by different colours. From left to right, we label these stable regions as 1, 2, and 3, respectively. These stable orbits actually are elliptic orbits around the Earth which are perturbed by the Moon. In order to interpret the amplitude parameter |$\gamma$| with a physically more intuitive quantity, we compute the mean orbit inclination of these orbits with respect to the Moon’s orbital plane. The upper abscissa of Fig. 5 is the orbit inclination. An obvious fact is that the practical stable region can extend up to 50°.

$The spatial stable regions around the TLP are characterized by the amplitude of the vertical component $\gamma$. The three main stable regions are marked with colours. From left to right, they are stable region 1, 2 and 3 respectively. The upper abscissa shows corresponding inclination angles with respect to the $\gamma$-values.$

Figure 5.

The spatial stable regions around the TLP are characterized by the amplitude of the vertical component |$\gamma$|⁠. The three main stable regions are marked with colours. From left to right, they are stable region 1, 2 and 3 respectively. The upper abscissa shows corresponding inclination angles with respect to the |$\gamma$|-values.

Open in new tab Download slide

3.3 Resonances

3.3.1 Fitting of the free frequencies

As is described above, the motions of the stable orbits are combinations of proper motion with frequencies |$\omega _l$|⁠, |$\omega _s$|⁠, and |$\omega _v$| and forced motion with frequencies |$\omega _1$|⁠, |$\omega _2$|⁠, |$\omega _3$|⁠, and |$\omega _4$|⁠. Combinations of the two kinds of motions can be analysed in the frequency domain using stable orbits data with long life spans (⁠|$15\,000 \ \mathrm{yr}$| in our work). For the vertical component |$\omega _v$|⁠, it is easy to extract it from the z-coordinate data of the stable orbits by simply finding the frequency with the maximum amplitude in the frequency domain. However, when it comes to the planar-free frequencies |$\omega _l$| and |$\omega _s$|⁠, the procedure becomes a little more complicated.

Fig. 6 displays two slices of FFT results with respect to different |$\gamma$|-values. The abscissa is the |$\gamma$|-value, the ordinate is the frequency value |$\omega$| and the plotted colours represent amplitudes of the corresponding frequency. Orbits of some |$\gamma$|-values are unstable and the corresponding regions in Fig. 6 are left blank. The two slices are chosen based on the fact that the two planar free frequencies |$\omega _l$| and |$\omega _s$| vary within the range of |$(0.2, 0.3)$| and |$(0.9, 1.0)$|⁠, respectively, for the |$\gamma$|-values displayed in Fig. 6. In Fig. 6(a), obviously there are some structures in the contour map. Three important structures are labelled by red, orange, and green lines, respectively. For each frequency structure, there are actually more than one frequencies in it. Taking the one labelled as red lines for example, our analysis finds that the difference between each nearby frequency in the structure is a fixed value which equals |$(\omega _4 - \omega _v)$|⁠. We heuristically take the middle frequency (labelled as a solid line) as the free frequency |$\omega _l$|⁠. As a result, the two dashed red lines are |$\omega _l \pm 5(\omega _4-\omega _v)$|⁠, respectively. The same phenomenon happens for other frequency structures in Fig. 6(a). As for Fig. 6(b), we also label out three important frequency structures using red, orange, and green lines, respectively. Also, taking the structure labelled by red lines as an example, we heuristically take the middle frequency (labelled as a solid line) as the frequency |$\omega _s$|⁠. Judging from Fig. 6, it seems that some frequency structures intersect in the contour maps, i.e. resonances happen. These intersection regions are also the regions where motion becomes unstable. So it is reasonable to speculate that these resonances form the boundaries of spatial stable regions. We will come back to these resonances in the forthcoming subsection. In this subsection, we show how we get these lines in Fig. 6.

$Two slices of FFT contour maps with respect to the vertical amplitude parameter $\gamma$ and the frequency $\omega$. Different colours in the contour maps indicate the amplitude of each frequency. Part (a) shows the result for $\omega = 0.2 - 0.3$, while part (b) shows the result for $\omega = 0.9 - 1.0$. Some important frequency families are plotted in the contour maps, with the solid lines as the middle frequency and the dashed lines as the boundary frequencies for the structures. From the top to the bottom, in part (a), frequencies at γ ═ 0.4 are ω = 2ω1 – 4ω3 + 2ωv, ω = ωl, ω = ω1 – 4ω3 + ωs + 2ωv and their belonging structures. In part (b), they are ω = ω1, ω = 3ω1 – 2ωv, ω = ωs and their belonging structures. Vertical blank regions correspond to the unstable orbits for which the FFT analysis fails.$

Figure 6.

Two slices of FFT contour maps with respect to the vertical amplitude parameter |$\gamma$| and the frequency |$\omega$|⁠. Different colours in the contour maps indicate the amplitude of each frequency. Part (a) shows the result for |$\omega = 0.2 - 0.3$|⁠, while part (b) shows the result for |$\omega = 0.9 - 1.0$|⁠. Some important frequency families are plotted in the contour maps, with the solid lines as the middle frequency and the dashed lines as the boundary frequencies for the structures. From the top to the bottom, in part (a), frequencies at γ ═ 0.4 are ω = 2ω₁ – 4ω₃ + 2ω_v, ω = ω_l, ω = ω₁ – 4ω₃ + ω_s + 2ω_v and their belonging structures. In part (b), they are ω = ω₁, ω = 3ω₁ – 2ω_v, ω = ω_s and their belonging structures. Vertical blank regions correspond to the unstable orbits for which the FFT analysis fails.

Open in new tab Download slide

$$\begin{eqnarray} \omega _s(\gamma) &=& \omega _{000}^{\lt s\gt } + \omega _{002}^{\lt s\gt } \gamma ^2 + \omega _{004}^{\lt s\gt } \gamma ^4 + \omega _{006}^{\lt s\gt } \gamma ^6 \\ \omega _l(\gamma) &=& \omega _{000}^{\lt l\gt } + \omega _{002}^{\lt l\gt } \gamma ^2 + \omega _{004}^{\lt l\gt } \gamma ^4 + \omega _{006}^{\lt l\gt } \gamma ^6 \\ \omega _v(\gamma) &=& \omega _{000}^{\lt v\gt } + \omega _{002}^{\lt v\gt } \gamma ^2 + \omega _{004}^{\lt v\gt } \gamma ^4 + \omega _{006}^{\lt v\gt } \gamma ^6 \end{eqnarray}$$

(12)

By numerically fitting the frequencies, we are able to get the coefficients in equation (12). Fig. 7 shows the fitting results and they are:

$$\begin{eqnarray} \omega _l &=& 0.305445 - 0.140781 \gamma ^2 + 0.078525 \gamma ^4 - 0.052583 \gamma ^6 \\ \omega _s &=& 0.928087 + 0.133203 \gamma ^2 - 0.143235 \gamma ^4 + 0.094466 \gamma ^6 \\ \omega _v &=& 1.004013 - 0.003468 \gamma ^2 - 0.003019 \gamma ^4 + 0.001357 \gamma ^6 \end{eqnarray}$$

(13)

Since the value of forced frequency |$\omega _4$| is fixed, we can use equation (13) to draw the frequency lines in Fig. 6. Judging from Fig. 6, the fitted results agree with the numerical structures in a perfect way.

$Part (a) shows the fitting result of frequency $\omega _l$. Part (b) shows the fitting result of frequency $\omega _s$. Part (c) shows the fitting result of frequency $\omega _v$.$

Figure 7.

Open in new tab Download slide

3.3.2 Resonances

Fig. 8 is a simplified version of Fig. 6 where only the spatial stable regions and the frequency lines are displayed. Resonances may happen between the forced frequencies |$\omega _1, \omega _2, \omega _3, \omega _4$| and the free frequencies |$\omega _l, \omega _s, \omega _v$|⁠, in the form of:

$$\begin{eqnarray} i\omega _1 + j\omega _2 + k\omega _3 + l\omega _4 + m\omega _s + n\omega _l + p\omega _v = 0. \end{eqnarray}$$

(14)

Figure 8.

Stable regions characterized in Fig. 5 and important frequencies marked in Fig. 6 are replotted in this figure. Part (a) shows the frequency slice of 0.2–0.3, while part (b) shows the frequency slice of 0.9–1.0. Stable regions are marked with the same colours as in Fig. 5 in this figure. Important frequencies are marked in solid and dashed lines with the same meanings as introduced in Fig. 6.

Open in new tab Download slide

According to Fig. 8, the regions where the frequency structures intersect are the possible resonances that shape the boundaries of the stable regions. Since we have three spatial stable regions, we have four boundaries – the lower boundary, the gap between the stable islands 1 and 2, the gap between the stable islands 2 and 3, and the upper boundary, so we have four main resonances. In the following, we label them as resonance <1>, <2>, <3>, and <4>, respectively. We have to remark that there is a family of frequencies instead of only one in each frequency structure, so the resonances happen between two families of frequencies instead of two frequencies. However, in the following, in order to label each frequency structure in a unique way, we identify it using the middle frequency (labelled as solid lines).

All resonances are shown in Table 2. In the following section, we will try to give physical interpretations to these resonances.

Table 2.

Open in new tab

Possible resonances that cause the instability of the orbits near the TLPs.

Index	i	k	m	n	p
<1>	2	\|$-$\|4	0	\|$-$\|1	2
<2>	1	\|$-$\|4	1	\|$-$\|1	2
<3>	3	0	\|$-$\|1	0	\|$-$\|2
<4>	1	0	\|$-$\|1	0	0

Index	i	k	m	n	p
<1>	2	\|$-$\|4	0	\|$-$\|1	2
<2>	1	\|$-$\|4	1	\|$-$\|1	2
<3>	3	0	\|$-$\|1	0	\|$-$\|2
<4>	1	0	\|$-$\|1	0	0

Table 2.

Open in new tab

Possible resonances that cause the instability of the orbits near the TLPs.

Index	i	k	m	n	p
<1>	2	\|$-$\|4	0	\|$-$\|1	2
<2>	1	\|$-$\|4	1	\|$-$\|1	2
<3>	3	0	\|$-$\|1	0	\|$-$\|2
<4>	1	0	\|$-$\|1	0	0

Index	i	k	m	n	p
<1>	2	\|$-$\|4	0	\|$-$\|1	2
<2>	1	\|$-$\|4	1	\|$-$\|1	2
<3>	3	0	\|$-$\|1	0	\|$-$\|2
<4>	1	0	\|$-$\|1	0	0

4 DISCUSSIONS

4.1 Interpretation of the resonances

The resonances in Section 3.3.2 are presented in the synodic frame while the physical meanings of these resonances may be better revealed in the GCRS frame. In the following, we will show that these resonances are closely related to the precession rates of the small body’s orbit, the Moon’s orbit, and the mean orbital motion of the Sun.

First, we explain the physical meanings of the forced frequencies. As already mentioned in Section 2.4, the frequency |$\omega _1=n_1-n_3$| represents the frequency of the angle |$\lambda _\mathrm{ M} - \tilde{\omega }_\mathrm{ M}$|⁠, so the frequency |$n_1 - \omega _1 = 1-\omega _1$| represents the precession rate of the Moon’s orbit argument of perigee. The frequency |$\omega _3=n_2$| represents the frequency of the angle |$\lambda _\mathrm{ M}-\lambda _\mathrm{ S}$|⁠, so the frequency |$n_1-\omega _3=1-\omega _3$| represents the mean orbital frequency of the Sun. The frequency |$\omega _4=n_1-n_4$| represents the frequency of the angle |$\lambda _\mathrm{ M}-\Omega _\mathrm{ M}$|⁠, so the frequency |$n_1-\omega _4=1-\omega _4$| represents the precession rate of the Moon’s orbit longitude of the ascending node.

Next, we explain the meanings of the free frequencies. For the long-period frequency |$\omega _l$|⁠, it is natural to interpretate it as the libration frequency of the |$1:1$| orbital resonance, i.e. |$\omega _l = \dot{\lambda } - \dot{\lambda }_\mathrm{ M}$| (Tan et al. 2023). For the short-period frequency |$\omega _s$|⁠, its physical meaning can be explained in the following way. Suppose that the small body is moving in a short-period orbit. If it is at the closest point to the Earth in the short period orbit, it is also at its perigee in the inertial frame. After a period of |$2\pi / \omega _s$| the small body returns to the point closest to the Earth, i.e. the perigee. However, the perigee itself precesses an angle of |$2\pi / \omega _s - 2\pi$| in the inertial frame. As a result, the precession rate of the small body’s argument of perigee is |$1-\omega _s$| (Hou, Scheeres & Liu 2014b). A similar explanation goes for the vertical frequency |$\omega _v$|⁠. Suppose that the small body is moving in a vertical periodic orbit. If it is at the |$z = 0$|⁠, |$\dot{z} \gt 0$| point in the synodic frame, it is also at its ascending node in the inertial frame. After a period of |$2\pi / \omega _v$|⁠, the small body returns to the point |$z=0$|⁠, |$\dot{z}\gt 0$|⁠, i.e. the ascending node. However, the ascending node itself precesses an angle of |$2\pi / \omega _v - 2 \pi$| in the inertial frame. As a result, the precession rate of the small body’s ascending node is |$1 - \omega _v$|⁠.

With the physical meanings of the forced and the free frequencies at hand, we can explain the resonances found in Section 3.3.2:

(1) For resonance <1> that shapes the lower boundary of the stable region, it is:

$$\begin{eqnarray} 2\omega _1 - 4\omega _3 - \omega _l + 2\omega _v &=& - 2(1-\omega _1) + 4(1-\omega _3) - 2(1-\omega _v) - \omega _l \\ &=& (4\dot{\lambda }_\mathrm{ S} - 2\dot{\tilde{\omega }}_\mathrm{ M} - 2\dot{\Omega }) - (\dot{\lambda } -\dot{\lambda }_\mathrm{ M}) = 0. \end{eqnarray}$$

So the resonance happens between the |$1:1$| resonance angle and the angle |$4\lambda _\mathrm{ S} - 2\tilde{\omega }_\mathrm{ M} -2\Omega$|⁠.

(2) For resonance <2> that shapes the gap between the stable islands 1 and 2, it is:

$$\begin{eqnarray} && \omega _1 - 4\omega _3 +\omega _s - \omega _l + 2\omega _v \\ &=& -(1 - \omega _1) + 4(1-\omega _3) - (1-\omega _s) - 2(1-\omega _v) - \omega _l\\ &=& (4\dot{\lambda }_\mathrm{ S} - \dot{\tilde{\omega }}_\mathrm{ M} -\dot{\tilde{\omega }} - 2\dot{\Omega }) - (\dot{\lambda } -\dot{\lambda }_\mathrm{ M}) = 0. \end{eqnarray}$$

So the resonance happens between the |$1:1$| resonance angle and the angle |$4\lambda _\mathrm{ S} - \tilde{\omega }_\mathrm{ M} - \tilde{\omega } - 2\Omega$|⁠.

(3) For resonance <3> that shapes the gap between stable islands 2 and 3, it is:

$$\begin{eqnarray} 3\omega _1 - 2\omega _v - \omega _s &=& -3(1 - \omega _1) + 2 (1-\omega _v) + (1-\omega _s) \\ &=& 2\dot{\Omega } + \dot{\tilde{\omega }} - 3\dot{\tilde{\omega }}_\mathrm{ M} = 0. \end{eqnarray}$$

It is one secular resonance between the small body and the Moon.

(4) For resonance <4> that shapes the upper boundary of the stable region, it is:

$$\begin{eqnarray} \omega _1 - \omega _s = -(1-\omega _1) - (1-\omega _s) = \dot{\tilde{\omega }} - \dot{\tilde{\omega }}_\mathrm{ M}. \end{eqnarray}$$

It is also one secular resonance between the small body and the Moon, well known as the Kozai resonance (Tan et al. 2020).

As a result, the resonances in Section 3.3.2 happen among the precession rates of the Moon and the small body, the orbital mean motion frequency of the Sun, and the libration frequency of the |$1:1$| orbital resonance.

4.2 Astronomical observation

A practical concern is whether there are natural bodies (including dust particles) moving on these spatial stable orbits. Suppose that they do exist, we have to answer the question that in which regions they are most likely to be observed. To answer this question, the following analysis is carried out in two steps. First, we pick up one spatial stable orbit and integrate the orbit for 2000 yr with an output interval of 12 h. At each moment, we compute its right ascension |$\alpha _t$| and declination |$\delta _t$| in the equatorial coordinate system as follows:

$$\begin{eqnarray} \begin{aligned} \alpha _t &= \arctan {\frac{R_y}{R_x}} \\ \delta _t & = \arcsin {\frac{R_z}{\sqrt{R_x^2 + R_y^2 + R_z^2}}}. \end{aligned} \end{eqnarray}$$

(15)

|$\boldsymbol{R}=(R_x, R_y, R_z)$| is the position vector of the target in the GCRS frame. We discretize the celestial sphere into a grid of |$360 \times 180$|⁠, with a bin size of |$1^{\circ }$| for both right ascension and declination. For each bin, we attribute an integer index |$p_{ij}$|⁠. Every 12 h, we calculate the target’s |$\alpha _t$| and |$\delta _t$| values. If |$(\alpha _t, \delta _t)$| is located in one bin, then the bin index is updated by |$p_{ij}=p_{ij}+1$|⁠. After 2000 yr, the probability distribution then writes:

$$\begin{eqnarray} P_{ij} = \frac{p_{ij}}{\sum _{ij}p_{ij}}. \end{eqnarray}$$

(16)

Taking the deviation of right ascension |$\Delta \alpha _t$| and declination |$\Delta \delta _t$| between objects and L5 point as abscissa and ordinate, the probability distribution of some example orbits is shown in Fig. 9. Fig. 9(a) shows the results for an orbit with a vertical component |$\gamma = 0.42868926$| (⁠|$i = 25.0874146^\circ$|⁠), while Fig. 9(b) shows an orbit with a vertical component |$\gamma = 0.66828450$| (⁠|$i = 41.1797354^\circ$|⁠). They belong to stable regions 1 and 2, respectively. Orbit in stable region 3 is not displayed because the difference between stable regions 2 and 3 is small in observation. Judging from Fig. 9, it seems that the probability distribution is higher at |$\Delta \delta _t = 0$| and at the largest values of |$|\Delta \delta _t|$|⁠. This is understandable. Assuming that the distance r is fixed, the relation |$\mathrm{d} z=r \cos \Delta \delta _t \cdot \mathrm{d} \Delta \delta _t$| holds. This means the same |$\mathrm{d} z$| corresponds to a much smaller |$\mathrm{d} \Delta \delta _t$| when |$\Delta \delta _t = 0$|⁠. As for the high probability distribution at the largest values of |$|\Delta \delta _t|$|⁠, it is simply due to the fact that |$\dot{z} = 0$|⁠.

$Observation probability distribution of two example orbits. Part (a) corresponds to an orbit with a vertical component $\gamma = 0.42868926$. Part (b) corresponds to an orbit with a vertical component $\gamma = 0.66828450$.$

Figure 9.

Observation probability distribution of two example orbits. Part (a) corresponds to an orbit with a vertical component |$\gamma = 0.42868926$|⁠. Part (b) corresponds to an orbit with a vertical component |$\gamma = 0.66828450$|⁠.

Open in new tab Download slide

Moreover, assuming that dust particles have even chances to move on these stable spatial orbits, we plot the probability distribution of all the orbits in a single plot, and the result is displayed in Fig. 10. Fig. 10(a) shows the result of all orbits, while Fig. 10(b) shows the result only with orbits in stable region 1. The distribution of all orbits shows an hourglass-shaped pattern with a right ascension deviation of |$-40^{\circ }\,\mathrm{ to}\, 40^{\circ }$| and a declination deviation of |$-50^\circ \,\mathrm{ to}\, 50^\circ$|⁠. There is a high-probability region close to the L5 point, which lies on the lunar orbit plane. There are also two high-probability regions around |$\pm 25^\circ$| in |$\Delta \delta _t$| related to the spatial stable orbits with vertical amplitudes |$\gamma$| ranging from 0.4 to 0.6. That also corresponds to the spatial stable orbits in stable region 1. Another structure to notice is that there are two obvious gaps near |$\pm 35^{\circ }$| in |$\Delta \delta _t$|⁠. Considering the difference between Figs 10(a) and (b), it is clear that these two gaps are caused by the gap between stable regions 1 and 2 as is stated in Section 3.2.

Figure 10.

Observation probability distribution by assuming that objects have even chances to move in the spatial stable regions. Part (a) shows the result for all spatial stable orbits, while part (b) shows the result only with orbits in stable region 1.

Open in new tab Download slide

Noting that the results in Fig. 10 are based on the assumption that dust particles have even chances to move on different stable orbits, the actual distribution of dust particles may be different. However, our result gives a rough estimation of positions where captured objects are more likely to be found. According to Fig. 10, objects are more likely to be found at three places, |$\Delta \delta _t = 0^{\circ }, +25^{\circ },\,\mathrm{ and}\, -25^{\circ }$|⁠. Although the region |$\Delta \delta _t = 0$| has the highest probability distribution, objects close to the lunar orbit are harder to observe against the zodiacal light.

4.3 The influence from the SRP

If there are objects residing around the TLPs of the Earth–Moon system, they should be of small sizes because it is not likely that large objects are missed by current observation technologies. For small objects such as dust particles in space, the SRP is a vital perturbation, and an analysis of its influence on the spatial stable region is helpful. Some work has been done considering the SRP on quasi-periodic orbits near colinear libration points (Chujo 2024) and on the stability of quasi-periodic orbits near the TLPs in the BCP model (Gimeno et al. 2024). It is even found that the SRP has a stabilizing effect for some orbits in the BCP and QBCP models (Jorba-Cuscó, Farrés & Jorba 2021).

In our work, the cannonball model is chosen since we have little knowledge about shape and other physical properties of objects around the TLPs. The SRP exerted on an object writes:

$$\begin{eqnarray} P_{\mathrm{srp}} = (1 + \sigma) \frac{P_0}{c}(\frac{R_0}{R_\mathrm{ S}})^2= (1 + \sigma)4.57 \times 10^{-6} (\frac{R_0}{R_\mathrm{ S}})^2 \end{eqnarray}$$

(17)

in |$[\mathrm{N} \, \mathrm{m}^{-2}]$|⁠. |$P_0 = 1367 \ \mathrm{W}\,\mathrm{m}^{-2}$| is the solar radiation flux at |$1 \ \mathrm{au}$|⁠, |$c = 299792458 \ \mathrm{m}\,\mathrm{s}^{-1}$| is the speed of light, |$R_0$| is the Sun–Earth mean distance (which equals |$1 \ \mathrm{au} = 149597870700 \ \mathrm{m}$|⁠), and |$R_\mathrm{ S}$| is the distance between the object and the Sun. |$\sigma$| is the reflectivity of the particle. With the cannonball model, the SRP acceleration acting on a homogenous spherical object with radius |$r_{\mathrm{s}}$| and density|$\rho _{\mathrm{s}}$| is:

$$\begin{eqnarray} \boldsymbol{a}_{\mathrm{srp}} = P_{\mathrm{srp}} \frac{A}{m} \hat{n}, \end{eqnarray}$$

where |$\frac{A}{m} = \frac{3}{4}\frac{1}{\rho _{\mathrm{s}} r_{\mathrm{s}}}$| and |$\hat{n}$| is defined as the unitary vector from the Sun to the object, i.e. |$\hat{n} = \frac{\boldsymbol{R} - \boldsymbol{R}_\mathrm{ S}}{|\boldsymbol{R} - \boldsymbol{R}_\mathrm{ S}|}$|⁠.

Following the same procedures introduced in Section 2.3, we first get a group of quasi-periodic orbits near the TLPs with their vertical components varying from 0 to 0.9 and then numerically integrate them to observe their sustainability in a longer integration time. We take the total integration time as |$\boldsymbol {T}_f = 15\,000 \ \mathrm{yr}$| and the same stability criterion |$x \ge -0.2$| to define stable orbits is chosen as in Section 3.2. One thing to notice is that, the three parameters |$r_\mathrm{s}$|⁠, |$\rho _\mathrm{s}$|⁠, and |$\sigma$| has a combined effect on the results of our analysis for the equation of the SRP acceleration can be rewritten as:

$$\begin{eqnarray} \begin{aligned} \boldsymbol{a}_{\mathrm{srp}} &= (1 + \sigma) \frac{P_0}{c}(\frac{R_0}{R_S})^2 \frac{A}{m} \hat{n} \\ &= 4.57\times 10^{-6}(1+\sigma)(\frac{R_0}{R_\mathrm{ S}})^2 \frac{A}{m} \hat{n} \\ &= [\frac{3}{4} 4.57\times 10^{-6}(\frac{R_0}{R_\mathrm{ S}})^2] [\frac{1}{\rho _\mathrm{s} r_\mathrm{s}} (1+\sigma)] \hat{n}. \end{aligned} \end{eqnarray}$$

We only takes one |$\sigma$| and |$\rho _\mathrm{s}$| setup (⁠|$\rho _\mathrm{s} = 1 \times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$| and |$\sigma = 0.44$|⁠) in our simulation. The density |$\rho _\mathrm{s}$| takes a value smaller than the averaged density of C-type asteroids that equals to |$1.7 \times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$| (Vernazza et al. 2021), and we choose |$\rho _\mathrm{s} = 1 \times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$| for convenience. The reflectivity value is chosen as |$\sigma = 0.44$| arbitrarily. The stable regions for objects of different sizes are displayed in Fig. 11. There are 12 stacks in the figure, and each stack displays the spatial stable region for a particle size. An obvious phenomenon is that the spatial stable regions shrink with decreasing object size. In our setup, for objects with sizes of 1 mm or smaller, practical stable regions lasting 15 000 yr no longer exist. Result under other circumstances can be easily derived from our setup. By keeping the value of |$[\frac{1}{\rho _\mathrm{s} r_\mathrm{s}} (1+\sigma)]$| unchanged, we are able to calculate the minimum |$r_\mathrm{s}$| for stable orbits to exist under other circumstances. Denoting the minimum |$r_{\mathrm{s}}$| value for stable orbits to exist as |$r_{\mathrm{s}}^{(\min)}$|⁠, and we have:

$$\begin{eqnarray} r_\mathrm{s}^{(\min)} = \frac{\rho _{\mathrm{s}0}}{(1 + \sigma _0)} \frac{(1+\sigma)}{\rho _\mathrm{s}} r_{\mathrm{s}0}, \end{eqnarray}$$

where |$\rho _{\mathrm{s}0} = 1\times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$|⁠, |$\sigma _0 = 0.44$|⁠, and |$\rho _{\mathrm{s}0} = 1 \ \mathrm{mm}$|⁠. As is shown in Fig. 12, different |$r_\mathrm{s}^{(\min)}$| values are shown with respect to different |$(\rho _{\mathrm{s}}, \sigma)$| setups. The value of |$\rho _{\mathrm{s}}$| varies from |$1\times 10^3 $| to |$8\times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$|⁠, which is between the density of water and that of some asteroids with the highest bulk densities, and the reflectivity |$\sigma$| varies from 0 to 1. Some ranges for C-, S-, and M-types, E-type of X-type asteroids by Tholen classification (Tholen 1984) are marked with dashed lines. Density data are from estimation by Carry (2012). Our studies indicate that objects of millimetre sizes in magnitude generally escape in less than hundreds of years due to the strong SRP influence.

$The spatial stable region around the TLP is characterized by the amplitude of the vertical component $\gamma$ for $\rho _\mathrm{s} = 1 \times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$ and $\sigma = 0.44$. The SRP is introduced in this part of our work. Homogenous spherical objects with different sizes in radius are analysed. The sizes of objects strongly affect the structure of stable regions and objects of millimetre size scale can hardly find any stable orbits to reside.$

Figure 11.

The spatial stable region around the TLP is characterized by the amplitude of the vertical component |$\gamma$| for |$\rho _\mathrm{s} = 1 \times 10^3 \ \mathrm{kg}\,\mathrm{m}^{-3}$| and |$\sigma = 0.44$|⁠. The SRP is introduced in this part of our work. Homogenous spherical objects with different sizes in radius are analysed. The sizes of objects strongly affect the structure of stable regions and objects of millimetre size scale can hardly find any stable orbits to reside.

Open in new tab Download slide

$The values of the minimum $r_{\mathrm{s}}$ for stable orbits to exist with respect to different $(\rho _{\mathrm{s}}, \sigma)$ setups. Values of $r_\mathrm{s}^{(\min)}$ are of millimetre sizes or less in magnitude. Ranges for C-, S-, and M-types, and E-type of the X-type asteroids by Tholen classification are marked with dashed lines.$

Figure 12.

The values of the minimum |$r_{\mathrm{s}}$| for stable orbits to exist with respect to different |$(\rho _{\mathrm{s}}, \sigma)$| setups. Values of |$r_\mathrm{s}^{(\min)}$| are of millimetre sizes or less in magnitude. Ranges for C-, S-, and M-types, and E-type of the X-type asteroids by Tholen classification are marked with dashed lines.

Open in new tab Download slide

Reconsidering the probability distribution with the SRP, we can still get similar structures as in Fig. 10. This is because the stable regions generally shrink with the decrease in the radius of particles, but similar boundaries and gaps between stable regions 1 and 2 still exist. That is to say, our proposed observational strategies hold considering the influence of the SRP.

4.4 Comparison with the BCP model

In this section, we display the differences in describing the spatial stable region between the BCP model and the ephemeris model. Fig. 13 shows the spatial stable region for the BCP model. The same procedure as the previous analysis is adopted to get this region. A number of spatial quasi-periodic orbits lasting 5 yr are obtained by numerically continuating from their counterparts vertical periodic orbits in the CRTBP model. Then, numerically integrate the quasi-periodic orbits for 15 000 yr to observe whether they can survive the integration.

Figure 13.

The spatial stable region around the TLP for the BCP model.

Open in new tab Download slide

4.4.1 The lower boundary

For the BCP, if we still use the criterion of |$x \gt -0.2$| as the stability criterion, we find that the lower boundary of the spatial stable region is zero. This is not right, because vertical orbits with small out-of-plane amplitudes are actually unstable. One example is displayed in Fig. 14(a). The out-of-plane amplitude of this orbit is 0.21084. With the time increasing, the orbit moves in a large bounded area of the x–y plane randomly, indicating that it is actually unstable. However, it does not escape, i.e. x is always larger than |$-0.2$| in 15 000 yr. To avoid this misjudgment, we use a more strict criterion of |$x \gt 0.2$| as the stability criterion for the lower boundary displayed in Fig. 13. The orbit at the lower boundary is displayed in Fig. 14(b). Obviously, the orbit does not expand within 15 000 yr. Comparing the lower boundary value of Fig. 13 with that of the ephemeris model (see Fig. 5), we find that the lower stable boundary of the BCP model is significantly smaller than that of the ephemeris model. This is because the lower stable boundary of the ephemeris model is determined by the resonance happens between the |$1:1$| resonance angle and the angle |$4\lambda _S - 2\tilde{\omega }_\mathrm{ M} -2\Omega$|⁠. However, the Moon’s orbit does not have a valid ascending node on the ecliptic due to zero inclination angle in the BCP model.

Figure 14.

(a) An unstable orbit with an out-of-plane amplitude of 0.21084. (b) The practical stable orbit at the lower boundary of the spatial stable region displayed in Fig. 13. The orbit in each frame is integrated for 1 yr.

Open in new tab Download slide

4.4.2 The gaps between stable islands

Judging from Fig. 13, it seems that there are no gaps in the spatial stable region for the BCP model, which is different from the result of the ephemeris model. Taking two example orbits (⁠|$\gamma = 0.600748$| and 0.712905), Fig. 15 displays their trajectories. For the ephemeris model, the two |$\gamma$|-values are located in the unstable gaps (see Fig. 5), but they are obviously stable in the BCP model. We believe this difference is due to the fact that the Moon’s orbit does not precess in the BCP model, while the resonances that determine the gaps between the stable islands in the ephemeris model are closely related with the Moon’s orbit precession.

Figure 15.

Two practical stable orbits with out-of-plane amplitudes of (a) 0.600748 and (b) 0.712905 for the BCP model. Both orbits are unstable in the ephemeris model.

Open in new tab Download slide

4.4.3 The upper boundary

The upper boundary is the same as that of the ephemeris model, which is determined by the Kozai resonance. No further details are presented.

5 CONCLUSION

In this study, for the first time, we described the distribution of spatial stable regions in the ephemeris model and investigated the resonance mechanisms that shape these stable regions.

Starting from vertical periodic orbits around the TLPs in the CRTBP model, we numerically continue them to quasi-periodic orbits in the ephemeris model. Numerical integration is implemented within 15 000 yr to characterize stable regions as shown in Fig. 5. A simple but efficient frequency analysis method is proposed in Section 2.4. Stable orbits with different vertical components are analysed using this method. Four basic frequencies of the system are recovered as |$(\omega _1, \omega _2, \omega _3, \omega _4)$| from the lunar orbit in DE431 and three free frequencies |$(\omega _l, \omega _s, \omega _v)$| are extracted for every stable orbit obtained. Relations between the vertical component value |$\gamma$| and the three free frequencies are calculated in Section 3.3.1. Resonances that shape the boundaries of the spatial stable region are found using the frequency analysis results in Section 3.3.2 and it turns out that these resonances happen among the precession rates of the Moon and the small body, the orbital mean motion frequency of the Sun, and the libration frequency of the |$1:1$| orbital resonance.

Objects such as dust particles may reside on or accumulate close to these stable regions. As a result, astronomical observation of potential objects in stable spatial regions is discussed. It is found that there is a higher probability of finding objects (if they exist) at the region close to the TLP points, or at two regions with about |$25^\circ$| difference in declination with respect to the TLP points. To avoid the influence from the zodiacal light, the two regions with about |$25^\circ$| difference in declination with respect to the TLP points are recommended. Further analysis indicates that the SRP has a negative influence on the size of the spatial stable regions. These regions gradually shrink with decreasing particle sizes, i.e. stronger SRP. The spatial stable region disappears for dust particles of millimeter sizes in magnitude or smaller.

We also compare the spatial stable regions in the ephemeris model with those of the BCP model. On the one hand, qualitative difference is found in the sense that there are no gaps in the spatial stable region in the BCP model. This difference is believed to be related to the non-existence of the Moon’s precession compared with the ephemeris model. On the other hand, a quantitive difference is found that the lower stable boundary of the BCP model is significantly smaller than that of the ephemeris model, which is believed to be related to the absence of the inclination of the lunar orbit on the ecliptic.

ACKNOWLEDGEMENTS

This work is supported by the National Natural Science Foundation of China (no. 12233003). We would like to express our sincere gratitude to the scientific editor, Dr Matija Cuk, for his valuable advice on the scope of the article and for providing us with the opportunity to resubmit. We also appreciate the reviewer’s insightful suggestion to expand the discussion on the parameter choices in the SRP section of our work.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

REFERENCES

Alfriend

K. T.

1970

Celest. Mech. Dyn. Astron.

351

Month:	Total Views:
November 2024	179
December 2024	74
January 2025	77
February 2025	64
March 2025	91
April 2025	58

Article Contents

Stability of spatial orbits around Earth–Moon triangular libration points

ABSTRACT

1 INTRODUCTION

2 MODELS AND METHODS

2.1 Ephemeris model description

2.2 Vertical periodic family

2.3 Numerical refinement

2.4 Frequency analysis method

2.5 The BCP model

3 RESULTS

3.1 Two orbit types

3.2 Spatial stable region

3.3 Resonances

3.3.1 Fitting of the free frequencies

3.3.2 Resonances

4 DISCUSSIONS

4.1 Interpretation of the resonances

4.2 Astronomical observation

4.3 The influence from the SRP

4.4 Comparison with the BCP model

4.4.1 The lower boundary

4.4.2 The gaps between stable islands

4.4.3 The upper boundary

5 CONCLUSION

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only