Analysis on the Near-Polar and Near-Circular Periodic Orbits Around the Moon with J2, C22 and Third-Body Perturbations

1. Introduction

Nowadays, with the increasing demand for the lunar exploration and the development of the aerospace technology, the cislunar space is becoming a new frontier for human activities. In order to carry out long-term scientific research missions near the Moon and save fuel for placing a space station, it is necessary to use the stable and approximately stable orbits, such as distant retrograde orbits (DROs), near-rectlinear halo orbits (NRHOs) and near-polar frozen orbits [1]. It is a common idea to firstly study the periodic orbits of the restricted three-body problem (RTBP), then continue the solutions to a generalized time-periodic model, and lastly to a high-fidelity ephemeris model[2]. It is important to study periodic orbits as they help us to understand the real motions, and there is a long history on this study. Ever since the excellent work of H. Poincaré, periodic orbits have aroused the attentions of many mathematicians, astronomers, celestial mechanicians, and so on[3,4].

The family of the near-polar and near-circular periodic orbits are interesting. According to Poincaré’s classification, these periodic orbits belong to the third-type of the first kind. The existence of this family in the elliptic RTBP is shown by Xu and Fu [5]. Theoretically, these periodic orbits are very close to one primary. Numerically , Xu and Song [6] found that these orbits can be of high altitude, and some other interesting phenomena were found as well. The eccentricity and the inclination of such a periodic orbit vary periodically and satisfy the Lidov-Kozai effect. For some values of the eccentricity of the outer orbit and the ratios of the mean motion resonances between the inner and outer orbits, the periodic orbits can be approximately linearly stable. This assures the the judgement that the high-altitude and near-circular polar frozen orbits are suitable to place a lunar station[7]. However, Xu and Song [6] did not give enough details about the near-polar and near-circular lunar-type periodic orbits. For a further study, it makes sense to study the these long-period periodic orbits in a more accurate gravity field model.

The motion of a lunar orbiter with the altitude in the range $(2000,16000)$km is mainly affected by the gravitation of the Moon-Earth system. The Moon-Earth system can be simplified as a two-body system, while the Moon can be considered as a triaxial ellipsoid, as the main non-spherical perturbations come from the $J_2,C_{22}$ terms. The zonal coefficient $J_2$ represents the size of the oblateness, and the sectorial coefficient $C_{22}$ measures the equatorial ellipticity. The obliquity of the Ecliptic plane and Moon’s Path is neglected. There is an interesting phenomenon named tidal locking for the Moon, that is, the rotation period and the revolution period are nearly equal, and the longest axis always passes the Earth. For lower-altitude lunar orbits, Lara et al. [8] studied the long-term dynamics of the near-polar frozen orbits by reducing a 50-degree zonal model with the third-body effect of the Earth. Saedeleer [9] explained the Hamiltonian system of a lunar orbiter in the Moon-Earth circular RTBP with $J_2,C_{22}$ perturbations, and studied the averaged system by the Lie-Deprit method. With the same model, Nie and Gurfil [10] studied the lunar frozen orbits in the first-order double-averaged system in Delaunay elements by Zeipel’s method. For the frozen orbits, the slow mean orbital elements keep nearly fixed such that the costs of the orbital corrections are reduced.

With the $J_2,C_{22}$ perturbations and several orders of the Legendre expansions of the third-body perturbations, Carvalho et al. [11] studied frozen orbits and the critical inclinations of the lunar satellites by analyzing the double-averaged system. Considering the 3rd-degree gravity harmonics of the Moon, Tzirti et al. [12] investigated the Poincaré sections, the Fast Lyapunov Indicator Maps and some families of periodic orbits. Tzirti et al. [13] studied the secular dynamics of low-altitude lunar orbiters with high-degree gravity models by the frequency analysis and investigated the eccentricity-inclination space. Considering a lunar orbiter perturbed by the $J_2,J_3,J_4$ terms, El-Salam and El-Bar [14] investigated the families of frozen orbits. In Sirwah et al. [15], the perturbing function is considered up to the seventh zonal harmonic and the third-body perturbation of the Earth in an elliptic inclined orbit. Then they numerically studied the frozen orbits with the arguments of the pericenter at $\pi /2,3\pi /2$. An efficient approach based on the grid search, parallelization and the evolution strategy is introduced for computing periodic orbits in Dena et al [17]. Franz and Russell [18] introduced the database of the symmetric periodic orbits near Moon in the model of the circular RTBP via the grid search and unsupervised learning clustering algorithm. Legnaro and Efthymiopoulos [19] distinguished three types of lunar orbits by the range of altitudes, and studied the secular dynamics especially the secular resonances and the eccentricity growth of lunar satellites according to different models.

In Section 2, we provide a new approach to do the numerical continuation of the near-polar and near-circular periodic orbits of the elliptic RTBP, and give some numerical examples of the lunar-type periodic orbits. In Section 3, we study the existence and stability of the near-polar and near-circular periodic orbits in the elliptic RTBP with the $J_2,C_{22}$ perturbations. Some numerical examples are also given. Finally, Section 5 concludes this work.

2. Elliptic RTBP

2.1. Scaled Hamiltonian System

With the background of the motion of a lunar orbiter in the cislunar space, we study the periodic orbits around the smaller primary in the elliptic RTBP. Let ${\mathtt{P}}_1$ and ${\mathtt{P}}_2$ represent two mass points. The relative orbit from ${\mathtt{P}}_1$ to ${\mathtt{P}}_2$ is a Keplerian orbit, with the semi-major axis $a_p$, the eccentricity $e_p$, the mean motion $n_p$, and the time of periapsis passage ${\tau }_0$. Choose the units such that the gravitational constant $\mathtt{G}=1$, the distance unit $a_p=1$, the total masses $m_1+m_2=1$, and the mean motion $n_p=1$. Let $\mu =m_2/(m_1+m_2)<1/2$. Set the motion plane of primaries as the reference plane. The direction of the major axis from ${\mathtt{P}}_1$ to ${\mathtt{P}}_2$ is set as the $q_1$-axis. Set ${\mathtt{P}}_2$ as the origion and establish the right-handed Cartesian coordinate system ${\mathtt{P}}_2-q_1{q}_2{q}_3$. The initial time is set as ${\tau }_0=0$ or ${\tau }_0=\pi$. The sketch figure of the frame can be referred to [5].

In the coordinate system ${\mathtt{P}}_2-q_1{q}_2{q}_3$, the position of ${\mathtt{P}}_1$ is denoted as ${\mathit{X}}_p\in {\mathbb{R}}^3$, which is also a solution of the planar Kepler problem

$$\begin{array}{c}\hfill {\ddot{\mathit{X}}}_p=-{\mathit{X}}_p·{\parallel {\mathit{X}}_p\parallel }^{-3}.\end{array}$$

Here $\parallel ·\parallel$ represent the Euclidean distance norm, and the time is t. We have $\parallel {\mathit{X}}_p\parallel =1-e_{p}cos{E}_p$. Denote $E_p$ as the eccentric anomaly. The vector ${\mathit{X}}_p$ is

$$\begin{array}{c}\hfill {\mathit{X}}_p={\left(X_1,X_2,0\right)}^{\mathrm{T}}={\left(cos{E}_p-e_p,\sqrt{1-e_p^2}sin{E}_p,0\right)}^{\mathrm{T}},\end{array}$$

where the upper $\mathrm{T}$ represents transpose. The position of the infinitesimal body is $\mathit{u}\in {\mathbb{R}}^3$, and its conjugate momentum is $\mathit{v}=\dot{\mathit{u}}$. The Hamiltonian dynamical system of this problem is

$$\begin{array}{c}\hfill {\mathcal{H}}^{\mathtt{P}\mathtt{2}}={\displaystyle \frac{1}{2}}{\parallel \mathit{v}\parallel }^2-{\displaystyle \frac{\mu }{\parallel \mathit{u}\parallel }}-{\displaystyle \frac{1-\mu }{\parallel \mathit{u}-{\mathit{X}}_p\parallel }}+(1-\mu ){\displaystyle \frac{{\mathit{u}}^{\mathrm{T}}{\mathit{X}}_p}{\parallel {\mathit{X}}_p{\parallel }^3}}.\end{array}$$

(1)

The canonical differential equation system is $\dot{\mathit{u}}=𝜕{\mathcal{H}}^{\mathtt{P}\mathtt{2}}/𝜕\mathit{v}$, $\dot{\mathit{v}}=-𝜕{\mathcal{H}}^{\mathtt{P}\mathtt{2}}/𝜕\mathit{u}$. The 2nd-order differential equation system can be written as

$$\begin{array}{c}\hfill \ddot{\mathit{u}}=-\mu {\displaystyle \frac{\mathit{u}}{{\parallel \mathit{u}\parallel }^3}}-(1-\mu ){\displaystyle \frac{\mathit{u}-{\mathit{X}}_p}{\parallel \mathit{u}-{\mathit{X}}_p{\parallel }^3}}-(1-\mu ){\displaystyle \frac{{\mathit{X}}_p}{\parallel {\mathit{X}}_p{\parallel }^3}}.\end{array}$$

(2)

When $\mu$ is very small, it is not convenient to get the numerical solution, as the orbit scale is too small.

The orbit of the infinitesimal body is called as the inner orbit, and the orbit of the relative orbit of ${\mathtt{P}}_1$ is the outer orbit. Denote the orbital elements of the inner orbit as $a_s$, $e_s$, $i_s$, ${\mathrm{\Omega }}_s$, ${\omega }_s$, ${\mathcal{l}}_s$, where ℓ represents the mean anomaly. According to the symplectic scaling method, the variables can be scaled as

$$\begin{array}{c}\hfill \mathit{u}={\epsilon }^2{\mu }^{1/3}\xi ,\mathit{v}={\epsilon }^{-1}{\mu }^{1/3}\eta ,t={\epsilon }^{3}s,\end{array}$$

(3)

where s is the new time, and the small parameter $\epsilon$ represents the closeness of the infinitesimal body to primary ${\mathtt{P}}_2$. The new Hamiltonian is

$$\begin{array}{cc}\hfill {\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}}(\xi ,\eta ,s)=& {\epsilon }^2{\mu }^{-2/3}{\mathcal{H}}^{\mathtt{P}\mathtt{2}}(\mathit{u},\mathit{v},t)={\displaystyle \frac{{\parallel \eta \parallel }^2}{2}}-{\displaystyle \frac{1}{\parallel \xi \parallel }}\hfill \\ \hfill & -{\epsilon }^2{\mu }^{-2/3}(1-\mu )\left({\displaystyle \frac{1}{\parallel {\epsilon }^2{\mu }^{1/3}\xi -{\mathit{X}}_p\left(s\right)\parallel }}-{\displaystyle \frac{{\epsilon }^2{\mu }^{1/3}{\xi }^{\mathrm{T}}{\mathit{X}}_p}{\parallel {\mathit{X}}_p{\parallel }^3}}\right).\hfill \end{array}$$

(4)

The new differential equation system becomes

$$\begin{array}{cc}\hfill {\xi }^{″}={\displaystyle \frac{{\mathrm{d}}^2\xi }{\mathrm{d}{s}^2}}=& -{\displaystyle \frac{\xi }{{\parallel \xi \parallel }^3}}-{\epsilon }^6(1-\mu ){\displaystyle \frac{\xi }{\parallel \mathit{u}-{\mathit{X}}_p{\parallel }^3}}\hfill \\ \hfill & +{\epsilon }^4(1-\mu ){\mu }^{-1/3}{\mathit{X}}_p\left({\displaystyle \frac{1}{\parallel \mathit{u}-{\mathit{X}}_p{\parallel }^3}}-{\displaystyle \frac{1}{\parallel {\mathit{X}}_p{\parallel }^3}}\right),\hfill \end{array}$$

(5)

where ${\xi }^{\prime }=\eta$, $\mathit{u}={\epsilon }^2{\mu }^{1/3}\xi$, and the prime represents the derivative about the scaled time s. Let ${\widehat{a}}_s$ be the scaled variable and $a_s={\epsilon }^2{\mu }^{1/3}{\widehat{a}}_s$. We have $n_s^2{a}_s^3=\mu$, and $n_s^2{\epsilon }^6\mu {\widehat{a}}_s=\mu$, so ${\widehat{n}}_s={\epsilon }^3{n}_s$. It is convenient to set ${\widehat{n}}_s=1$ if ${\epsilon }^3=n_p/n_s$. Then we get ${\widehat{n}}_p={\epsilon }^3$. The Kepler equation for the outer orbit satisfies

$$\begin{array}{c}\hfill E_p-e_{p}sin{E}_p=n_{p}t={\widehat{n}}_{p}s.\end{array}$$

The numerical solution of Eq. (5) can be calculated by the integration effectively.

2.2. Symmetry and Periodicity

The scaled Hamiltonian system $\widehat{\mathcal{H}}$ keeps the same symmetries as the original Hamiltonian system ${\mathcal{H}}^{\mathtt{P}\mathtt{2}}$. One time reversing symmetry ${\mathcal{R}}_1$ with respect to ${\xi }_1$-axis is recalled.

$$\begin{array}{c}\hfill {\mathcal{R}}_1:({\xi }_1,{\xi }_2,{\xi }_3,{\eta }_1,{\eta }_2,{\eta }_3,s)\to ({\xi }_1,-{\xi }_2,-{\xi }_3,-{\eta }_1,{\eta }_2,{\eta }_3,-s).\end{array}$$

There exist periodic third-body perturbations in the lunar-type orbits in the elliptic RTBP. The desired periodic orbits can be continued from the two uncoupled Kepler orbits. The ratio of the mean motions of the uncoupled inner and outer Kepler orbits is set to be ${\widehat{n}}_s/{\widehat{n}}_p=n_s/n_p={\epsilon }^3$, and ${\epsilon }^3$ should be a small rational number. ${\epsilon }^3$ is set to equal $k/j$, where $k,j\in \mathbb{N}$ and $k\ll j$. This means that the inner orbit revolves j circles while the outer orbit revolves k circles. In order to understand the ${\mathcal{R}}_1$-symmetric periodic solution, a proposition is summarized as follows.

Lemma 1 

([5]). For the Hamiltonian system (1) of the elliptic RTBP. The Lagrangian set ${\mathcal{L}}_1$ with respect to ${\mathcal{R}}_1$-symmetry is

$$\begin{array}{c}\hfill {\mathcal{L}}_1=\{(\mathbf{u},\mathbf{v},t):u_2=u_3=v_1=0,t=0mod\pi \}.\end{array}$$

If a solution $Z(Z_1,t)$ satisfies $Z_1\in {\mathcal{L}}_1$ and $Z(Z_1,k\pi )\in {\mathcal{L}}_1$ with $k\in \mathbb{N}$, then the solution is ${\mathcal{R}}_1$-symmetric and periodic with a period $T=2k\pi$. For the Hamiltonian system (4), the corresponding Lagrangian set ${\mathcal{L}}_1^{\left(1\right)}$ is

$$\begin{array}{c}\hfill {\mathcal{L}}_1^{\left(1\right)}=\{(\xi ,\eta ,s):{\xi }_2={\xi }_3={\eta }_1=0,s=0mod\pi \},\end{array}$$

and the scaled periodic solution has a period $\widehat{T}=2j\pi$ with $j\in \mathbb{N}$. The ratio of the mean motions of the outer and inner orbits is $k/j$. The periodicity conditions used in this paper can be written as

$$\begin{array}{c}\hfill {\widehat{Z}}_1\in {\mathcal{L}}_1^{\left(1\right)},\widehat{Z}({\widehat{Z}}_1,j\pi )\in {\mathcal{L}}_1^{\left(1\right)}.\end{array}$$

(6)

According to the description in Xu and Song [6], there exist both near-polar and planar ${\mathcal{R}}_1$-symmetric near-circular lunar-type periodic orbits in the elliptic RTBP. However, Xu and Song [6] did not give enough details of the numerical continuation of such periodic orbits. The difficulty lies at the fact that the accumulated integration errors are relatively large when the infinitesimal body is very close to primary ${\mathtt{P}}_2$ and the integration time is long. For the Hamiltonian system (1), the initial values to be continued can be written as

$$\begin{array}{c}\hfill Z_0={\left(u_1^0,u_2^0,u_3^0,v_1^0,v_2^0,v_3^0\right)}^{\mathrm{T}}={\left(\pm {\mu }^{1/3}{(k/j)}^{2/3},0,0,0,0,\pm {\mu }^{1/3}{(k/j)}^{-1/3}\right)}^{\mathrm{T}}.\end{array}$$

The scale of the lunar-type orbits are small but the velocities are relatively big. For the Hamiltonian system (4), the initial values to be continued can be written as

$$\begin{array}{c}\hfill {\widehat{Z}}_0={\left({\xi }_1^0,{\xi }_2^0,{\xi }_3^0,{\eta }_1^0,{\eta }_2^0,{\eta }_3^0\right)}^{\mathrm{T}}={\left(\pm 1,0,0,0,0,\pm 1\right)}^{\mathrm{T}}.\end{array}$$

In this paper, the scaled variables are used and more numerical continuation results can be achieved. The continuation scenario is still based on the combination of the periodicity conditions and Broyden’s method with a line search. This means that ${\widehat{Z}}_0$ is expected to be continued to

$$\begin{array}{c}\hfill {\widehat{Z}}_1={\left(\pm 1+{\delta }_1,0,0,0,{\delta }_2,\pm 1+{\delta }_3;s_0\right)}^{\mathrm{T}},s_0=0,j\pi ,\end{array}$$

such that

$$\begin{array}{c}\hfill {\xi }_2({\delta }_1,{\delta }_2,{\delta }_3,s_{T/2})=0,{\xi }_3({\delta }_1,{\delta }_2,{\delta }_3,s_{T/2})=0,{\eta }_1({\delta }_1,{\delta }_2,{\delta }_3,s_{T/2})=0,\end{array}$$

where ${\delta }_1,{\delta }_2,{\delta }_3\in \mathbb{R}$ are small quantities, and $s_{T/2}=s_0+j\pi$.

2.3. Some Periodic Orbits

The numerical solutions of Eq.(5) are calculated by the variable step-size Runge-Kutta 7-8 routine with the double precision. The routines of Broyden’s method are referred to Press et al. [16] and the precision is guaranteed at ${10}^{-8}$. So the periodic orbits can keep the precision of ${10}^{-8}$. Let ${\mathtt{P}}_1,{\mathtt{P}}_2$ represent the Earth and the Moon, respectively. Let $\mu =0.0121505843947$, $e_p=0.0549$. Set $j/k$ as the ratio of the mean motions between the uncoupled inner and outer orbits. The real length of the semi-major axis of the outer orbit is about $A=328900.5597$km. The real length of the semi-major axis of the inner orbit is about $A{a}_0$km with $a_0={\mu }^{1/3}{(k/j)}^{2/3}$. The altitude is defined as the difference of $A{a}_0$ and the real length of the longer equatorial semi-major axis $1738.1$km. The high altitude zone is defined in the range about between 5000km and 20000km. In this zone, the third-body perturbation of the Earth is dominant. In the theoretical proof of the existence of the near-polar and near-circular periodic orbits in the elliptic RTBP, the small parameter ${\epsilon }^3=k/j$ is supposed to be small enough. One question is that how small is ${\epsilon }^3$ is when the periodic orbits can be calculated. Such periodic orbits are numerically investigated with ${\epsilon }^3$ in the range $9\le j/k\le 150$. In Figure 1 and Figure 2, it is found that the argument of pericenter $\omega$ rotates and oscillates. Figure 3 is for the case $j/k=37$. All these orbital elements are very symmetric. The details of the figures are shown in the image captions. Some more initial values for the periodic orbits are shown in Table 1. The type of the initial values are defined by the signs of the values of ${\xi }_1,{\eta }_3$ and $cos{E}_p$.

3. Elliptic RTBP with $J_2,C_{22}$ Perturbations

3.1. The Model and Numerical Experiment

The usual orbital elements are the semi-major axis a, the eccentricity e, the inclination i, the longitude of the ascending node $\mathrm{\Omega }$, the argument of the periapsis $\omega$, and the mean anomaly M. The eccentric anomaly is notated as $E(e,M)$, and the true anomaly is notated as $f(e,M)$. If just the perturbations of $J_2,C_{22}$ are added to the RTBP, the application of this model is for the case of the altitude about in the range $(2000,5000)$km. The non-spherical perturbing function is usually expressed by orbital elements. It is necessary to transform the perturbing function into the form of Cartesian coordinates for the convenience of the computation of the near-circular periodic orbits. For a massless artificial satellite in the Moon-centered inertial coordinate frame, the position is notated as $\mathit{u}$, and the conjugate momentum is $\mathit{v}=\dot{\mathit{u}}$. The Hamiltonian system of a lunar orbiter can be written in a form of perturbations of the RTBP,

$$\begin{array}{c}\hfill {\mathcal{H}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}={\mathcal{H}}^{\mathtt{P}\mathtt{2}}(\mathit{u},\mathit{v},t)+{\mathcal{H}}_{J2}-{\mathcal{H}}_{C22},\end{array}$$

(7)

where

$$\begin{array}{cc}\hfill & {\mathcal{H}}_{J2}=J_2{\displaystyle \frac{\mu a_m^2}{r^3}}{\mathcal{P}}_2(sin\phi ),{\mathcal{H}}_{C22}=C_{22}{\displaystyle \frac{\mu a_m^2}{r^3}}{\mathcal{P}}_{22}(sin\phi )·cos2(\mathrm{\Omega }-f_p+\psi ),\hfill \\ \hfill & a_m\approx 1738.1/328900.5597,J_2\approx 2.0322356\times {10}^{-4},C_{22}\approx 2.2381388\times {10}^{-5},\hfill \\ \hfill & r=\parallel \mathit{u}\parallel ,n_p=1,sin\phi ={\displaystyle \frac{u_3}{r}},{\mathcal{P}}_2(sin\phi )={\displaystyle \frac{3u_3^2}{2r^2}}-{\displaystyle \frac{1}{2}},\hfill \\ \hfill & cos\phi cos\psi =cos(f+\omega ),cos\phi sin\psi =sin(f+\omega )cosi,\hfill \\ \hfill & {\mathcal{P}}_{22}(sin\phi )=3{cos}^2\phi ,cos2(\mathrm{\Omega }-f_p+\psi )=2{cos}^2(\mathrm{\Omega }-f_p+\psi )-1,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill & {\mathcal{P}}_{22}(sin\phi )cos2(\mathrm{\Omega }-f_p+\psi )=6{[cos\phi cos(\mathrm{\Omega }-f_p+\psi )]}^2-3{cos}^2\phi \hfill \\ \hfill & =6{\{cos\phi ·[cos(\mathrm{\Omega }-f_p)cos\psi -sin(\mathrm{\Omega }-f_p)sin\psi ]\}}^2-3{cos}^2\phi )\hfill \\ \hfill & =6{[cos(f+\omega )cos(\mathrm{\Omega }-f_p)-sin(f+\omega )cosisin(\mathrm{\Omega }-f_p)]}^2-3{cos}^2\phi \hfill \\ \hfill & ={\displaystyle \frac{6}{r^2}}{(u_{1}cos{f}_p+u_{2}sin{f}_p)}^2-3{cos}^2\phi .\hfill \end{array}$$

In order to better understand the angles $\phi$ and $\psi$, Figure 4 is recomanded for reading. In the Lunar inertial coordinate frame $P_2-q_1{q}_2{q}_3$, the longitude of the ascending node of the infinitesimal satellite is $\mathrm{\Omega }$, and the latitude is $\phi$. The prime meridian is set at the direction of the longest semi-major axis of the ellipsoid. The longitude of the satellite is $\mathrm{\Omega }-f_p+\psi$ in the rotating frame $P_2-x_1{x}_2{x}_3$. The Earth is always on the $x_1$-axis. Besides, $f_p$ is notated as the true anomaly of the outer orbit, ${\mathcal{P}}_2(·)$ is the 2nd-order Legendre polynomial, and ${\mathcal{P}}_{22}(·)$ is the unnormalized associative Legendre polynomial.

With the help of the symplectic scaling,

$$\begin{array}{c}\hfill \mathit{u}\to {\epsilon }^2{\mu }^{1/3}\xi ,\mathit{v}\to {\epsilon }^{-1}{\mu }^{1/3}\eta ,t\to {\epsilon }^{3}s,J_2\to {\epsilon }^6{\tilde{J}}_2,C_{22}\to {\epsilon }^6{\tilde{C}}_{22},\end{array}$$

and the Legendre polynomial expansion[5], the Hamiltonian (7) becomes

$$\begin{array}{cc}\hfill {\tilde{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}& ={\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}}(\xi ,\eta ,s)+{\epsilon }^6({\widehat{\mathcal{H}}}_{J2}-{\widehat{\mathcal{H}}}_{C22})\hfill \\ \hfill & ={\widehat{\mathcal{H}}}_0^{\mathtt{P}\mathtt{2}}+{\epsilon }^6({\widehat{\mathcal{H}}}_1+{\widehat{\mathcal{H}}}_{J2}-{\widehat{\mathcal{H}}}_{C22})+\mathcal{O}\left({\epsilon }^8\right),\hfill \end{array}$$

(8)

where

$$\begin{array}{cc}\hfill & {\widehat{\mathcal{H}}}_1=-{\displaystyle \frac{(1-\mu )r^2}{\parallel {\mathit{X}}_p{\parallel }^3}}{\mathcal{P}}_2(cos\theta ),r=\parallel \xi \parallel ,\hfill \\ \hfill & {\widehat{\mathcal{H}}}_{J2}={\tilde{a}}_m^2{\tilde{J}}_2\left({\displaystyle \frac{3{\xi }_3^2}{2r^5}}-{\displaystyle \frac{1}{2r^3}}\right),\hfill \\ \hfill & {\widehat{\mathcal{H}}}_{C22}={\tilde{a}}_m^2{\tilde{C}}_{22}\left({\displaystyle \frac{6{\tilde{x}}_1^2}{r^5}}-{\displaystyle \frac{3}{r^3}}+{\displaystyle \frac{3{\xi }_3^2}{r^5}}\right),\hfill \\ \hfill & {\tilde{x}}_1={\xi }_{1}cos{f}_p+{\xi }_{2}sin{f}_p,{\tilde{a}}_m={\epsilon }^{-2}{\mu }^{-1/3}{a}_m,\hfill \\ \hfill & cos\theta ={\displaystyle \frac{\xi }{\parallel \xi \parallel }}·{\displaystyle \frac{{\mathit{X}}_p}{\parallel {\mathit{X}}_p\parallel }}={\displaystyle \frac{1}{r}}({\xi }_{1}cos{f}_p+{\xi }_{2}sin{f}_p)={\displaystyle \frac{{\tilde{x}}_1}{r}}.\hfill \end{array}$$

The relations between the scaled rectangular coordinates and the osculating orbital elements are

$$\begin{array}{c}\hfill \left\{\begin{array}{c}r=a(1-e^2)/(1+ecosf),\hfill \\ u_1=rcos(f+\omega )cos\mathrm{\Omega }-rsin(f+\omega )cosisin\mathrm{\Omega },\hfill \\ u_2=rcos(f+\omega )sin\mathrm{\Omega }+rsin(f+\omega )cosicos\mathrm{\Omega },\hfill \\ u_3=rsin(f+\omega )sini,\hfill \end{array}\right.\end{array}$$

Here, ${\tilde{x}}_1=rcos\theta$, and $cos\theta$ can be expressed as

$$\begin{array}{cc}\hfill cos\theta =& cos(f+\omega )cos(\mathrm{\Omega }-f_p)-cosisin(f+\omega )sin(\mathrm{\Omega }-f_p)\hfill \\ \hfill =& (I_1+I_2)cos(f+\omega )cos(\mathrm{\Omega }-f_p)-(I_1-I_2)sin(f+\omega )sin(\mathrm{\Omega }-f_p)\hfill \\ \hfill =& I_{1}cos(f+\omega +\mathrm{\Omega }-f_p)+I_{2}cos(f+\omega -\mathrm{\Omega }+f_p),\hfill \end{array}$$

where

$$\begin{array}{c}\hfill I_1={\displaystyle \frac{1}{2}}\left(1+cosi\right),I_2={\displaystyle \frac{1}{2}}\left(1-cosi\right).\end{array}$$

According to the knowledge of the complex variable functions, the real triangular functions can be replaced by the complex exponential functions, such that the powers of the real triangular functions are easy to be calculated. By this method, we get

$$\begin{array}{cc}\hfill {cos}^2\theta =& {\displaystyle \frac{I_1^2}{2}}cos2(f+\omega +\mathrm{\Omega }-f_p)+{\displaystyle \frac{I_2^2}{2}}cos2(f+\omega -\mathrm{\Omega }+f_p)+{\displaystyle \frac{I_1^2+I_2^2}{2}}\hfill \\ \hfill & +I_1{I}_2\left[cos2(f+\omega )+cos2(\mathrm{\Omega }-f_p)\right].\hfill \end{array}$$

Now, it is convenient to apply the Hamiltonian system (8) both for the numerical computation of the periodic orbits and for the analysis of the first-order perturbed system. By numerical experiment, It is found that $\omega$ is near 0 when $cos{E}_p=-1$ and $\omega$ is near $\pi$ when $cos{E}_p=1$. An example is given in Figure 5. Some more initial values of these frozen periodic orbits can be found in Table 2. It is interesting to explain the phenomenon in the first-order double-averaged system.

3.2. First-Order Averaged System

The scaled canonical Delaunay elements are introduced in order to do the Von Zeipel transform and get the first-order averaged system.

$$\begin{array}{c}\hfill \mathcal{l}=M,g=\omega ,h=\mathrm{\Omega },L=\sqrt{a},G=L\sqrt{1-e^2},H=Gcosi.\end{array}$$

The completely expansion of the perturbed system in the Delaunay elements is difficult, so the mean anomaly ℓ is usually contained in $E\left(e\right(L,G),M)$ and $f\left(e\right(L,G),M)$. Hansen coefficients are usually taken in use in the elliptic expansion[3],

$$\begin{array}{c}\hfill {\left({\displaystyle \frac{r}{a}}\right)}^{n}exp\mathrm{i}mf=\sum _{k=-\infty }^{\infty }{X}_k^{n,m}\left(e\right)exp\mathrm{i}kM.\end{array}$$

We remove the hats above the Hamiltonian notations in (8). Let $r_2=\parallel {\mathit{X}}_p\parallel$. After elimination of the fast variable M, we get

$$\begin{array}{cc}\hfill {\overline{\mathcal{H}}}_1& =-{\displaystyle \frac{(1-\mu )a^2}{r_2^3}}\{{\displaystyle \frac{15}{4}}{e}_1^2\left[I_1^{2}cos2(\omega +\mathrm{\Omega }-f_p)+I_2^{2}cos2(\omega -\mathrm{\Omega }+f_p)+I_1{I}_{2}cos2\omega \right]\hfill \\ \hfill & +\left(1+{\displaystyle \frac{3}{2}}{e}^2\right)\left[{\displaystyle \frac{3}{4}}(I_1^2+I_2^2)+{\displaystyle \frac{3}{2}}{I}_1{I}_{2}cos2(\mathrm{\Omega }-f_p)-{\displaystyle \frac{1}{2}}\right]\}.\hfill \\ \hfill {\overline{\mathcal{H}}}_{J2}& ={\displaystyle \frac{{\tilde{J}}_2{\tilde{a}}_m^2}{L^3{G}^3}}\left({\displaystyle \frac{3}{4}}{sin}^{2}i-{\displaystyle \frac{1}{2}}\right),{\overline{\mathcal{H}}}_{C22}=-{\displaystyle \frac{3{\tilde{C}}_{22}{\tilde{a}}_m^2}{2L^3{G}^3}}·{sin}^{2}i·cos2(\mathrm{\Omega }-f_p).\hfill \end{array}$$

Consider that the second fast variable is $f_p$, these terms above can be averaged again as

$$\begin{array}{cc}\hfill {\overline{\overline{\mathcal{H}}}}_1& =-{\displaystyle \frac{(1-\mu )L^4}{L_2^3{G}_2^3}}\left\{{\displaystyle \frac{15}{16}}\left(1-{\displaystyle \frac{G^2}{L^2}}\right)\left(1-{\displaystyle \frac{H^2}{G^2}}\right)cos2g+{\displaystyle \frac{1}{16}}\left(5-{\displaystyle \frac{3G^2}{L^2}}\right)\left({\displaystyle \frac{3H^2}{G^2}}-1\right)\right\},\hfill \\ \hfill {\overline{\overline{\mathcal{H}}}}_{J2}& ={\overline{\mathcal{H}}}_{J2}={\displaystyle \frac{{\tilde{J}}_2{\tilde{a}}_m^2}{L^3{G}^3}}\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{3H^2}{4G^2}}\right),{\overline{\overline{\mathcal{H}}}}_{C22}=0,L_2=1,G_2=\sqrt{1-e_p^2}.\hfill \end{array}$$

The canonical differential equations of the first-order averaged system are

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}\mathcal{l}}{\mathrm{d}s}}& ={\displaystyle \frac{1}{L^3}}+{\epsilon }^6{\displaystyle \frac{𝜕({\overline{\overline{\mathcal{H}}}}_1+{\overline{\overline{\mathcal{H}}}}_{J2})}{𝜕L}},{\displaystyle \frac{\mathrm{d}g}{\mathrm{d}s}}={\epsilon }^6{\displaystyle \frac{𝜕({\overline{\overline{\mathcal{H}}}}_1+{\overline{\overline{\mathcal{H}}}}_{J2})}{𝜕G}},{\displaystyle \frac{\mathrm{d}h}{\mathrm{d}s}}={\epsilon }^6{\displaystyle \frac{𝜕({\overline{\overline{\mathcal{H}}}}_1+{\overline{\overline{\mathcal{H}}}}_{J2})}{𝜕H}},\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}s}}& =0,{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}s}}=-{\displaystyle \frac{15(1-\mu )L^4}{8G_2^3}}{e}^2{sin}^{2}i·sin2g,{\displaystyle \frac{\mathrm{d}H}{\mathrm{d}s}}=0.\hfill \end{array}$$

(9)

For a further step, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}\mathcal{l}}{\mathrm{d}s}}=& {\displaystyle \frac{1}{L^3}}-{\epsilon }^6{\displaystyle \frac{(1-\mu )L^3}{G_2^3}}\left[{\displaystyle \frac{15}{8}}(1+e^2){sin}^{2}i·cos2g+{\displaystyle \frac{1}{8}}(7+3e^2)(3{cos}^{2}i-1)\right]\hfill \\ \hfill & -{\epsilon }^6{\displaystyle \frac{3{\tilde{J}}_2{\tilde{a}}_m^2}{4L^4{G}^3}}(1-3{cos}^{2}i),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}g}{\mathrm{d}s}}=& -{\epsilon }^6{\displaystyle \frac{3(1-\mu )L^4}{8G_2^{3}G}}\left[5(-{sin}^{2}i+e^2)cos2g+(1-e^2-5{cos}^{2}i)\right]\hfill \\ \hfill & -{\epsilon }^6{\displaystyle \frac{3{\tilde{J}}_2{\tilde{a}}_m^2}{4L^3{G}^4}}(1-3{cos}^{2}i)+{\epsilon }^6{\displaystyle \frac{3{\tilde{J}}_2{\tilde{a}}_m^2}{2L^3{G}^4}}{cos}^{2}i,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}h}{\mathrm{d}s}}=& -{\epsilon }^6{\displaystyle \frac{3(1-\mu )L^4}{8G_2^{3}G}}\left(-5e^{2}cos2g+2+3e^2\right)cosi-{\epsilon }^6{\displaystyle \frac{3{\tilde{J}}_2{\tilde{a}}_m^2}{2L^3{G}^4}}cosi.\hfill \end{array}$$

In the near-polar and near-circular frozen orbits, the orbital elements a, e, i, $\mathrm{\Omega }$ are almost constants with small periodic amplitudes. Substitute the initial unperturbed orbital elements into the differential equation about $\dot{g}$, we have

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}g}{\mathrm{d}s}}\approx 1.5{\epsilon }^{6}cos2g-0.75J_2{\tilde{a}}_m^2.\end{array}$$

If g is fixed at 0 or $\pi$, we have ${\epsilon }^6\approx 1.5\times {10}^{-5}$, but it can just happen for the low-altitude lunar orbits.

In order to understand the periodic behaviors of these frozen orbits, we resort to Poincaré-Delaunay elements, as they are effective for near-circular and near-polar orbits.

$$\begin{array}{c}\hfill \left\{\begin{array}{ccccccccc}{Q}_1& =& \mathcal{l}+g,& Q_2& =& -\sqrt{2(L-G)}sing,& Q_3& =& h,\\ P_1& =& L,& P_2& =& \sqrt{2(L-G)}sing,& P_3& =& H.\end{array}\right.\end{array}$$

The first-order double averaged system in the Poincaré-Delaunay elements can be written as

$$\begin{array}{c}\hfill {\mathcal{K}}_f=-{\displaystyle \frac{1}{2P_1^2}}+{\epsilon }^6\left(-c_1{\mathcal{K}}_1+c_2{\mathcal{K}}_2\right),c_1={\displaystyle \frac{(1-\mu )}{16G_2^3}},c_2={\displaystyle \frac{{\tilde{J}}_2{\tilde{a}}_m^2}{4}},\end{array}$$

(10)

where

$$\begin{array}{cc}\hfill {\mathcal{K}}_1& =P_1^4(15{\mathcal{A}}_1{\mathcal{A}}_2+{\mathcal{A}}_3),{\mathcal{K}}_2={\displaystyle \frac{8}{P_1^3{G}_{two}^3}}-{\displaystyle \frac{96P_3^2}{P_1^3{G}_{two}^5}},\hfill \\ \hfill {\mathcal{A}}_1& =1-{\displaystyle \frac{G_{two}^2}{4P_1^2}}-{\displaystyle \frac{4P_3^2}{G_{two}^2}}+{\displaystyle \frac{P_3^2}{P_1^2}},{\mathcal{A}}_2={\displaystyle \frac{2P_2^2}{P_2^2+Q_2^2}}-1,\hfill \\ \hfill {\mathcal{A}}_3& =-5+{\displaystyle \frac{3G_{two}^2}{4P_1^2}}+{\displaystyle \frac{60P_3^2}{G_{two}^2}}-{\displaystyle \frac{9P_3^2}{P_1^2}},G_{two}=2P_1^2-P_2^2-Q_2^2=2G.\hfill \end{array}$$

The canonical differential equations can be calculated with the help of a symbolic computation software like the wxMaxima.

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{Q}_1}{\mathrm{d}s}}& ={\displaystyle \frac{𝜕{\mathcal{K}}_f}{𝜕P_1}}={\displaystyle \frac{1}{P_1^3}}+{\epsilon }^6\left(-c_1{\displaystyle \frac{𝜕{\mathcal{K}}_1}{𝜕P_1}}+c_2{\displaystyle \frac{𝜕{\mathcal{K}}_2}{𝜕P_1}}\right),\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{Q}_i}{\mathrm{d}s}}& ={\displaystyle \frac{𝜕{\mathcal{K}}_f}{𝜕P_i}}={\epsilon }^6\left(-c_1{\displaystyle \frac{𝜕{\mathcal{K}}_1}{𝜕P_i}}+c_2{\displaystyle \frac{𝜕{\mathcal{K}}_2}{𝜕P_i}}\right),i=2,3,\hfill \\ \hfill {\displaystyle \frac{\mathrm{d}{P}_2}{\mathrm{d}s}}& =-{\displaystyle \frac{𝜕{\mathcal{K}}_f}{𝜕Q_2}}=-{\epsilon }^6\left(-c_1{\displaystyle \frac{𝜕{\mathcal{K}}_1}{𝜕Q_2}}+c_2{\displaystyle \frac{𝜕{\mathcal{K}}_2}{𝜕Q_2}}\right),{\displaystyle \frac{\mathrm{d}{P}_j}{\mathrm{d}s}}={\displaystyle \frac{𝜕{\mathcal{K}}_f}{𝜕Q_j}}=0,j=1,3.\hfill \end{array}$$

The key parts of the differential equations of the first-order double averaged Hamiltonian system ${\mathcal{K}}_f$ (10) are as follows.

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{Q}_1}{\mathrm{d}s}}=& {\displaystyle \frac{1}{P_1^3}}-{\epsilon }^6{\displaystyle \frac{3c_2}{P_1^4}}\left({\displaystyle \frac{8}{G_{two}^2}}-{\displaystyle \frac{96P_3^2}{G_{two}^5}}\right)+{\epsilon }^6{\displaystyle \frac{c_2}{P_1^3}}\left({\displaystyle \frac{960P_3^2}{G_{two}^6}}-{\displaystyle \frac{32}{G_{two}^3}}\right)\hfill \\ \hfill & -{\epsilon }^{6}4c_1{P}_1^3(15{\mathcal{A}}_1{\mathcal{A}}_2+{\mathcal{A}}_3)-{\epsilon }^6{c}_1{P}_1^4[15{\mathcal{A}}_2·({\displaystyle \frac{G_{two}^2}{2P_1^3}}+{\displaystyle \frac{16P_3^2}{G_{two}^3}}\hfill \\ \hfill & -{\displaystyle \frac{G_{two}}{P_1^2}}-{\displaystyle \frac{2P_3^2}{P_1^3}})-{\displaystyle \frac{3G_{two}^2}{2P_1^3}}-{\displaystyle \frac{240P_3^2}{G_{two}^3}}+{\displaystyle \frac{3G_{two}}{P_1^2}}+{\displaystyle \frac{18P_3^2}{P_1^3}}].\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{Q}_2}{\mathrm{d}s}}=& {\epsilon }^6{\displaystyle \frac{c_2}{P_1^3}}\left({\displaystyle \frac{32}{G_{two}^3}}-{\displaystyle \frac{960P_3^2}{G_{two}^6}}\right)P_2-{\epsilon }^6{c}_1{P}_1^4[15\left({\displaystyle \frac{G_{two}}{P_1^2}}-{\displaystyle \frac{16P_3^2}{G_{two}^3}}\right){\mathcal{A}}_2{P}_2\hfill \\ \hfill & +\left({\displaystyle \frac{240P_2{P}_3^2}{G_{two}^3}}-{\displaystyle \frac{60Q_2^2{\mathcal{A}}_1}{{(P_2^2+Q_2^2)}^2}}-{\displaystyle \frac{3G_{two}}{P_1^2}}\right)P_2].\hfill \end{array}$$

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}{Q}_3}{\mathrm{d}s}}=-{\epsilon }^6{\displaystyle \frac{192c_2{P}_3}{P_1^3{G}_{two}^5}}-{\epsilon }^6{c}_1{P}_1^4\left[15{\mathcal{A}}_2\left({\displaystyle \frac{2P_3}{P_1^2}}-{\displaystyle \frac{8P_3}{G_{two}^2}}\right)+{\displaystyle \frac{120P_3}{G_{two}^2}}-{\displaystyle \frac{18P_3}{P_1^2}}\right].\end{array}$$

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}{P}_2}{\mathrm{d}s}}=& -{\epsilon }^6{\displaystyle \frac{c_2}{P_1^3}}\left({\displaystyle \frac{32}{G_{two}^3}}-{\displaystyle \frac{960P_3^2}{G_{two}^6}}\right)Q_2+{\epsilon }^6{c}_1{P}_1^4[15\left({\displaystyle \frac{G_{two}}{P_1^2}}-{\displaystyle \frac{16P_3^2}{G_{two}^3}}\right){\mathcal{A}}_2{Q}_2\hfill \\ \hfill & +\left({\displaystyle \frac{240P_3^2}{G_{two}^3}}-{\displaystyle \frac{60P_2^2{\mathcal{A}}_1}{{(P_2^2+Q_2^2)}^2}}-{\displaystyle \frac{3G_{two}}{P_1^2}}\right)Q_2].\hfill \end{array}$$

In the first-order double-averaged system, $P_1$ and $P_3$ are constants. Set $P_1=1$ and $P_3=0$, we have $Q_3\equiv Q_3\left(0\right)$. The differential equations about $Q_2,P_2$ can be simplified as

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}s}}\left(\begin{array}{c}{Q}_2\\ P_2\end{array}\right)={\epsilon }^6\left\{{\displaystyle \frac{32c_2}{G_{two}^3}}-c_1\left[15{\mathcal{A}}_2{G}_{two}-246+18(P_2^2+Q_2^2)\right]\right\}\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)\left(\begin{array}{c}{Q}_2\\ P_2\end{array}\right).\end{array}$$

Suppose $Q_2\left(0\right)=P_2\left(0\right)=0$, then $Q_2\left(t\right)=P_2\left(t\right)=0$. We have ${\displaystyle \frac{\mathrm{d}{Q}_1}{\mathrm{d}s}}=1+{\epsilon }^6(8c_1-10c_2)$. So there exist polar-type and circular periodic orbits in the first-order double-averaged system.

4. Linear Stability

The Hamiltonian system ${\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}$ is

$$\begin{array}{cc}\hfill {\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}=& {\displaystyle \frac{\parallel \eta \parallel }{2}}-{\displaystyle \frac{1}{\parallel \xi \parallel }}-{\epsilon }^2{\mu }^{-2/3}(1-\mu )\left({\displaystyle \frac{1}{r_3}}-{\displaystyle \frac{{\epsilon }^2{\mu }^{1/3}{\xi }^{\mathrm{T}}{\mathit{X}}_p}{r_2^3}}\right)\hfill \\ \hfill & +J_2{\tilde{a}}_m^2\left({\displaystyle \frac{3{\xi }_3^2}{2r^5}}-{\displaystyle \frac{1}{2r^3}}\right)-C_{22}{\tilde{a}}_m^2\left({\displaystyle \frac{6{\tilde{x}}_1^2}{r^5}}-{\displaystyle \frac{3}{r^3}}+{\displaystyle \frac{3{\xi }_3^2}{r^5}}\right),\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill & r=\parallel \xi \parallel ,r_2=\parallel {\mathit{X}}_p\left(s\right)\parallel ,r_3=\parallel {\epsilon }^2{\mu }^{1/3}\xi -{\mathit{X}}_p\parallel ,\hfill \\ \hfill & {\tilde{x}}_1={\xi }_{1}cos{f}_p+{\xi }_{2}sin{f}_p,f_p=f_p\left(E_p(n_{p}s,e_p),e_p\right).\hfill \end{array}$$

The 2nd-order differential equations are

$$\begin{array}{cc}\hfill {\displaystyle \frac{{\mathrm{d}}^2\xi }{\mathrm{d}{s}^2}}=& -{\displaystyle \frac{\xi }{r^3}}-{\epsilon }^6(1-\mu ){\displaystyle \frac{\xi }{r_3^3}}+{\epsilon }^4(1-\mu ){\mu }^{-1/3}{\mathit{X}}_p\left({\displaystyle \frac{1}{r_3^3}}-{\displaystyle \frac{1}{r_2^3}}\right)-J_2{\tilde{a}}_m^2{\mathit{V}}_1+C_{22}{\tilde{a}}_m^2{\mathit{V}}_2,\hfill \end{array}$$

where ${\tilde{a}}_m={\epsilon }^{-2}{\mu }^{-1/3}{a}_m$, and

$$\begin{array}{c}\hfill {\mathit{V}}_1={\displaystyle \frac{𝜕}{𝜕\xi }}\left({\displaystyle \frac{3{\xi }_3^2}{2r^5}}-{\displaystyle \frac{1}{2r^2}}\right),{\mathit{V}}_2={\displaystyle \frac{𝜕}{𝜕\xi }}\left({\displaystyle \frac{6{\tilde{x}}_1^2}{r^5}}-{\displaystyle \frac{3}{r^3}}+{\displaystyle \frac{3{\xi }_3^2}{r^5}}\right).\end{array}$$

The canonical differential equations of the Hamiltonian system ${\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}$ are

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}\xi }{\mathrm{d}s}}=\eta ,{\displaystyle \frac{\mathrm{d}\eta }{\mathrm{d}s}}=-{\displaystyle \frac{𝜕{\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}}{𝜕\xi }}={\displaystyle \frac{{\mathrm{d}}^2\xi }{\mathrm{d}{s}^2}}.\end{array}$$

(11)

The solution is notated as $\xi =\xi (s,{\xi }_0,{\eta }_0)$, $\eta =\eta (s,{\xi }_0,{\eta }_0)$. The fundamental solution matrix of the linear variational equations is notated as

$$\begin{array}{c}\hfill \mathrm{\Phi }=\mathrm{\Phi }(s,{\xi }_0,{\eta }_0)=\left(\begin{array}{cc}{\displaystyle \frac{𝜕\xi }{𝜕{\xi }_0}}& {\displaystyle \frac{𝜕\xi }{𝜕{\eta }_0}}\\ {\displaystyle \frac{𝜕\eta }{𝜕{\xi }_0}}& {\displaystyle \frac{𝜕\eta }{𝜕{\eta }_0}}\end{array}\right).\end{array}$$

The system of the linear variational equations is

$$\begin{array}{c}\hfill {\displaystyle \frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}s}}=\left(\begin{array}{cc}{\mathit{O}}_3& {\mathit{I}}_3\\ \mathit{B}& {\mathit{O}}_3\end{array}\right)\mathrm{\Phi },\mathit{B}=-{\displaystyle \frac{{𝜕}^2{\widehat{\mathcal{H}}}^{\mathtt{P}\mathtt{2}\mathtt{JC}}}{𝜕\xi 𝜕{\xi }_0}},\end{array}$$

(12)

where ${\mathit{I}}_3$ is a $3\times 3$ identical matrix, ${\mathit{O}}_3$ is a $3\times 3$ zero matrix, and

$$\begin{array}{c}\hfill \mathit{B}={\mathit{B}}_1+{\epsilon }^6(1-\mu ){\mathit{B}}_2-{\epsilon }^6(1-\mu ){\mathit{B}}_3-J_2{\tilde{a}}_m^2{\mathit{B}}_4+C_{22}{\tilde{a}}_m^2{\mathit{B}}_5,\end{array}$$

$$\begin{array}{c}\hfill {\mathit{B}}_1=\left(\begin{array}{ccc}{b}_{11}& b_{12}& b_{13}\\ b_{12}& b_{22}& b_{23}\\ b_{13}& b_{23}& b_{33}\end{array}\right),b_{ii}=-{\displaystyle \frac{1}{r^3}}+{\displaystyle \frac{3{\xi }_i^2}{r^5}},b_{ij}={\displaystyle \frac{3{\xi }_i{\xi }_j}{r^5}},\end{array}$$

$$\begin{array}{c}\hfill {\mathit{B}}_2=\left(\begin{array}{ccc}{\widehat{b}}_{11}& {\widehat{b}}_{12}& {\widehat{b}}_{13}\\ {\widehat{b}}_{12}& {\widehat{b}}_{22}& {\widehat{b}}_{23}\\ {\widehat{b}}_{13}& {\widehat{b}}_{23}& {\widehat{b}}_{33}\end{array}\right),{\widehat{b}}_{ii}=-{\displaystyle \frac{1}{r_3^3}}+{\displaystyle \frac{3u_i(u_i-X_i)}{r_3^5}},{\widehat{b}}_{ij}={\displaystyle \frac{3u_i(u_j-X_j)}{r^5}},\end{array}$$

$$\begin{array}{c}\hfill {\mathit{B}}_3=\left(\begin{array}{ccc}{\tilde{b}}_{11}& {\tilde{b}}_{12}& {\tilde{b}}_{13}\\ {\tilde{b}}_{12}& {\tilde{b}}_{22}& {\tilde{b}}_{23}\\ {\tilde{b}}_{13}& {\tilde{b}}_{23}& {\tilde{b}}_{33}\end{array}\right),{\tilde{b}}_{ij}={\displaystyle \frac{3X_i(u_j-X_j)}{r_3^5}},\end{array}$$

$$\begin{array}{cc}\hfill {\mathit{B}}_4& =\left(\begin{array}{ccc}{\overline{b}}_{11}& {\overline{b}}_{12}& {\overline{b}}_{13}\\ {\overline{b}}_{12}& {\overline{b}}_{22}& {\overline{b}}_{23}\\ {\overline{b}}_{13}& {\overline{b}}_{23}& {\overline{b}}_{33}\end{array}\right),{\overline{b}}_{ii}={\displaystyle \frac{2}{r_3^5}}-{\displaystyle \frac{15{\xi }_3^2}{2r^7}}-{\displaystyle \frac{12{\xi }_i^2}{r^8}}+{\displaystyle \frac{105{\xi }_i^2{\xi }_3^2}{2r^9}},\hfill \\ \hfill {\overline{b}}_{12}& ={\displaystyle \frac{105{\xi }_1{\xi }_2{\xi }_3^2}{2r^9}}-{\displaystyle \frac{12{\xi }_1{\xi }_2}{r^8}},{\overline{b}}_{i3}=-{\displaystyle \frac{15{\xi }_i{\xi }_3}{2r^7}}-{\displaystyle \frac{12{\xi }_i{\xi }_3}{r^8}}+{\displaystyle \frac{105{\xi }_i^2{\xi }_3^2}{2r^9}},\hfill \\ \hfill {\overline{b}}_{33}& ={\displaystyle \frac{3}{r_3^5}}+{\displaystyle \frac{2}{r_3^6}}-{\displaystyle \frac{75{\xi }_3^2}{2r^7}}-{\displaystyle \frac{12{\xi }_3^2}{r^8}}+{\displaystyle \frac{105{\xi }_3^4}{2r^9}},i=1,2,\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\mathit{B}}_5& =\left(\begin{array}{ccc}{\check{b}}_{11}& {\check{b}}_{12}& {\check{b}}_{13}\\ {\check{b}}_{12}& {\check{b}}_{22}& {\check{b}}_{23}\\ {\check{b}}_{13}& {\check{b}}_{23}& {\check{b}}_{33}\end{array}\right),{\check{b}}_i={\displaystyle \frac{9}{r^5}}-{\displaystyle \frac{15{\xi }_3^2+30{\tilde{x}}_1^2+45{\xi }_i^2}{r^7}}+{\displaystyle \frac{105{\xi }_i^2{\xi }_3^2+210{\xi }_i^2{\tilde{x}}_1^2}{r^9}},\hfill \\ \hfill {\check{b}}_{11}& ={\check{b}}_1+{\displaystyle \frac{12{cos}^2{f}_p}{r^5}}-{\displaystyle \frac{120{\tilde{x}}_1{\xi }_{1}cos{f}_p}{r^7}},{\check{b}}_{22}={\check{b}}_2+{\displaystyle \frac{12{sin}^2{f}_p}{r^5}}-{\displaystyle \frac{120{\tilde{x}}_1{\xi }_{2}sin{f}_p}{r^7}},\hfill \\ \hfill {\check{b}}_{12}& ={\displaystyle \frac{6sin2f_p}{r^5}}-{\displaystyle \frac{60{\tilde{x}}_1({\xi }_{1}sin{f}_p+{\xi }_{2}cos{f}_p)+45{\xi }_1{\xi }_2}{r^7}}+{\displaystyle \frac{105{\xi }_1{\xi }_2{\xi }_3^2+210{\tilde{x}}_1^2{\xi }_1{\xi }_2}{r^9}},\hfill \end{array}$$

$$\begin{array}{cc}\hfill {\check{b}}_{i3}& =-{\displaystyle \frac{60{\tilde{x}}_1{\xi }_{3}cos{f}_p+75{\xi }_i{\xi }_3}{r^7}}+{\displaystyle \frac{105{\xi }_i{\xi }_3^2+210{\tilde{x}}_1^2{\xi }_i{\xi }_3}{r^9}},i=1,2,\hfill \\ \hfill {\check{b}}_{33}& ={\displaystyle \frac{15}{r^5}}-{\displaystyle \frac{120{\xi }_3^2+30{\tilde{x}}_1^2}{r^7}}+{\displaystyle \frac{105{\xi }_3^4+210{\tilde{x}}_1^2{\xi }_3^2}{r^9}}.\hfill \end{array}$$

The initial values for Eq. (12) is the $6\times 6$ identical matrix. The Equations (11) and (12) should be integrated together. The scaled time for the numerical integration is the period $T_{jk}=2j\pi /k$ of a periodic orbit. However, the whole information of a symmetric periodic orbit can be gotten by the integration of a half period.

the monodromy matrix can by calculated by the following formula[6]

$$\begin{array}{c}\hfill \mathrm{\Phi }\left(2j\pi \right)=\mathrm{\Phi }\left(j\pi \right){\mathrm{\Gamma }}_1{\mathrm{\Phi }}^{-1}\left(j\pi \right){\mathrm{\Gamma }}_1,{\mathrm{\Gamma }}_1=\mathbf{diag}\{1,-1,-1,-1,1,1\}.\end{array}$$

There are three pairs of conjugate eigenvalues for the matrix $\mathrm{\Phi }\left(2j\pi \right)$. If one periodic orbit is linearly stable, the eigenvalues are in the unit circle. If not, at least one eigenvalue will be far away from the unit circle. One index to describe the stability is the summation of the moduli of the eigenvalues. The eigenvalues of $\mathrm{\Phi }\left(2j\pi \right)$ are also called characteristic multipliers. For just a few examples, the summation of the moduli of the multipliers can be calculated directly from the eigenvalues of $\mathrm{\Phi }\left(2j\pi \right)$. With the $J_2,C_{22}$ and the third-body perturbations, the stability index for the periodic orbit of $(j/k=38,+++)$-type is about 24, and 114 for the case $(j/k=50,+++)$, 590 for the case $(j/k=60,+++)$ case, and 3778 for the case $(j/k=70,+++)$. In the elliptic RTBP, the linear stability index is about 7 for the case $(j/k=9,+++)$. This reveals that it saves fuel at the high-altitude orbits.

5. Discussion

The paper provides a lot information about the near-polar, near-circular, lunar-type periodic orbits in the Moon-Earth elliptic restricted three-body problem with $J_2,C_{22}$ and third-body perturbations. Some periodic orbits are calculated and the orbital dynamics are well explained by the first-order double-averaged system. The symplectic scaling technique is introduced, and the small parameter ${\epsilon }^3=k/j$ represents the small ratio of the mean resonances between the inner orbit of the infinitesimal body and the outer orbit of the relative motion of the Earth. The linear stability index is introduced and the linear variational equations are calculated. The scaled Hamiltonian system can also be applied to the study of the near-planar lunar orbits.

More work can be done based on this paper. For high-altitude orbits, the perurbation from the Sun can be added to the Hamiltonian system, and the $J2,C22$ perturbations can be neglected. The quasi-bicircular model of the Moon-Earth-Sun system can be applied. The asymmetric periodic orbits can be studied with the aim of finding stable orbits. It is interesting to study the analytical solutions and the evolution of the orbits. The formulas of this paper can also be adjusted to study the orbits around a satellite of one planet.