On the Dominant Lunisolar Perturbations for Long-Term Eccentricity Variation: The Case of Molniya Satellite Orbits

Abstract

The aim of this work is to investigate the main dominant terms of lunisolar perturbations that affect the orbital eccentricity of a Molniya satellite in the long term. From a practical point of view, these variations are important in the context of space situational awareness—for instance, to model the long-term evolution of artificial debris in a highly elliptical orbit or to design a reentry end-of-life strategy for a satellite in a highly elliptical orbit. The study assumes a doubly averaged model including the Earth’s oblateness effect and the lunisolar perturbations up to the third-order expansion. The work presents three important novelties with respect to the literature. First, the perturbing terms are ranked according to their amplitudes and periods. Second, the perturbing bodies are not assumed to move on circular orbits. Third, the lunisolar effect on the precession of the argument of pericenter is analyzed and discussed. As an example of theoretical a application, we depict the phase space description associated with each dominant term, taken as isolated, and we show which terms can apply to the relevant dynamics in the same region.

1. Introduction

On 23 April 1965, the first Molniya-1 spacecraft was launched by the former Soviet Union [1]. After this, many others were set in orbit until 2004 and, currently, an total of 36 Molniyas are still in orbit as non-cooperant satellites [2]—that is, space debris. The orbits of these satellites were designed ad hoc for Russian communication needs and are now considered a class of special orbits around the Earth: the Molniya orbits.

As a matter of fact, the Molniya orbits have inspired particular interest among the scientific community for their dynamical features, namely the high eccentricity $e\ge 0.7$, the critical inclination $i\approx 63.4^{\circ }$ and the orbital period of approximately 12 h. The reason for this configuration can be ascribed to the need to cover the Russian territory, located at a high latitude. The chosen inclination not only ensures coverage of the range of latitude of interest, but also the removal of the effect of precession of the line of apsides induced by the oblateness of the Earth. In this way, the perigee and the apogee of the satellite remain almost frozen in time, according to the initial argument of pericenter $\omega ={270}^{\circ }$ chosen due to the Russian latitude [3,4]. Moreover, a Molniya satellite revolves two times around the Earth every day: its orbital period and the Earth’s rotation period are commensurable and this fact produces a tesseral resonance, sometimes called inappropriately mean motion resonance (e.g., [5]).

Because of its dynamical features, a Molniya satellite undergoes several orbital perturbations. The low value of the altitude of the perigee, approximately 500 km [6], gives a non-negligible effect due to the atmospheric drag, which deeply affects the evolution of the semi-major axis, besides the tesseral resonance. Moreover, the satellite spends most of the time at high altitudes; hence, the lunisolar effect plays a fundamental role in the eccentricity evolution on a timescale larger than the satellite orbital period.

From this brief overview, it is evident that Molniya orbits are placed in a particularly rich dynamical environment; therefore, caution is needed. Our intent is to provide a robust analytical insight that allows us to better understand the nature of the eccentric motion. Throughout the discussion, the semi-major axis will be assumed constant. Moreover, we stress that we are not interested in the role played by geopotential harmonics beyond the second zonal harmonics $J_2$. For further details on their role, the interested reader can refer to [7,8,9].

We will focus on a doubly averaged Hamiltonian model including the secular oblateness effect and the expansion of the lunar and solar disturbing functions up to the third order (Section 2.1), thus not considering the simplification of circular orbits for the perturbing bodies.

By exploiting an analytical approach based on the Hamiltonian theory (Section 2.2), it is possible to identify the lunisolar contributions that dominate the eccentricity (and inclination) dynamics in the long term. The analysis provided can be especially useful when perturbative models of different orders are built, particularly to estimate the possible emergence of chaos. In the literature, the choice of the perturbative terms is generally heuristic, based only on the amplitudes of the harmonic coefficients of the disturbing potential. Here, the amplitudes of the harmonic coefficients, the corresponding periods and the ratio between the amplitudes and the frequency are estimated for the specific case of the Molniya orbital region (Section 3.1), following [10]. These quantities will be our “observables”. The analysis is actually general and similar arguments can be applied to other initial conditions.

By identifying the main perturbing terms and the regime in eccentricity and inclination where they can play a role, we will be able to estimate a maximum overlapping region, comprising more than two resonant effects. This could be the basis for a future analysis of the possible chaotic evolution of the eccentricity.

Another novelty with respect to past works is that lunisolar effects are usually studied with a second-order doubly averaged model but the geocentric orbits of the third bodies are taken as circular. Under this assumption, the third-order disturbing function vanishes, thanks to the analytical expressions of the eccentricity functions appearing in it [11], and the argument of pericenter of the perturbing bodies is not well-defined; hence, they do not appear in the power series development of the disturbing functions. Here, we will provide the general treatment.

We also remark that most of the past literature on lunisolar perturbation effects approximated the slow frequencies of the problem with the precession rate caused by the Earth’s oblateness (see, e.g., [12]). Such approximation is generally both convenient and accurate enough but, as shown later in this paper, it seems not appropriate for the Molniya dynamics. Indeed, the lunisolar contribution on the slow frequencies is not negligible for the dynamics of the argument of the perigee of the satellite, because of the critical inclination.

Finally, this work is part of a wider project in which the goal is to analyze in depth the relevant astrodynamical properties of the Molniya satellites’ orbits. Specifically, the analysis presented here runs parallel to the one carried out in [2], whose numerical results agree with our general conclusions, as long as the assumption of a constant semi-major axis (i.e., outside the drag regime) holds. The interested reader can find there a more detailed numerical analysis; here, we will focus on the analytical part.

2. Theoretical Background

The theoretical background of our investigation is collected here for the sake of completeness, even if several concepts are well-known among the celestial mechanics community and the treatment may be long. In this way, the interested researcher has all the necessary information to follow our approach and results, and to reproduce them, if needed. In the following, the Keplerian orbital elements will be referred to as: the semi-major axis a, the eccentricity e, the inclination i, the argument of the perigee $\omega$, the longitude of the ascending node $\Omega$ and the mean anomaly M. Moreover, the subscripts ⊕, ☾ and ⊙ will denote the parameters of the Earth, the Moon and the Sun, respectively. The satellite’s elements will be denoted by no subscript.

2.1. The Disturbing Functions

Following [12], the solar disturbing function can be written as:

$$\begin{array}{c}{\mathcal{R}}_{\odot }=\sum _{l=2}^{\infty }\sum _{m,p,q=0}^l\sum _{j,r=-\infty }^{\infty }{\mu }_{\odot }\left(\frac{a^l}{a_{\odot }^{l+1}}\right){ϵ}_m\frac{(l-m)!}{(l+m)!}{F}_{lmp}\left(i\right)F_{lmq}\left(i_{\odot }\right)H_{lpj}\left(e\right)G_{lqr}\left(e_{\odot }\right)\times \hfill \\ cos\left[(l-2p+j)M-(l-2q+r)M_{\odot }+(l-2p)\omega -(l-2q){\omega }_{\odot }+m(\Omega -{\Omega }_{\odot })\right],\hfill \end{array}$$

(1)

where ${\mu }_{\odot }$ is the gravitational parameter of the Sun, the Keplerian elements of both the satellite and the Sun are written with respect to the equatorial reference plane, $F_{lmp}$ and $F_{lmq}$ are Kaula’s inclination functions, $H_{lpj}$ and $G_{lqr}$ are Hansen coefficients and ${ϵ}_m$ can be found in [12].

On the other hand, the motion of the Moon around the Earth is significantly perturbed by the gravitational attraction of the Sun; hence, the lunar inclination, node and argument of the perigee with respect to the celestial equator evolve as nonlinear functions of time. To overcome this issue, usually, a mixed reference plane is adopted, where the elements of the satellite are written with respect to the equatorial plane, while the elements of the Moon are referred to as the ecliptic, so that $i_{☾}$ is approximately constant and ${\omega }_{☾}$ and ${\Omega }_{☾}$ are approximately linear functions of time [12]. In this way, the lunar disturbing function reads

$$\begin{array}{c}{\mathcal{R}}_{☾}=\sum _{l=2}^{\infty }\sum _{m,p,s,q=0}^l\sum _{j,r=-\infty }^{+\infty }{(-1)}^{m+s}{(-1)}^{k_1}\frac{{\mu }_{☾}{ϵ}_m{ϵ}_s}{2a_{☾}}\frac{(l-s)!}{(l+m)!}{\left(\frac{a}{a_{☾}}\right)}^l{F}_{lmp}\left(i\right)F_{lsq}\left(i_{☾}\right)\times \hfill \\ H_{lpj}\left(e\right)G_{lqr}\left(e_{☾}\right)\{{(-1)}^{k_2}{U}_l^{m,-s}\left(ϵ\right)cos[(l-2p+j)M+(l-2q+r)M_{☾}+\hfill \\ (l-2p)\omega +(l-2q){\omega }_{☾}+m\Omega +s({\Omega }_{☾}-\frac{\pi }{2})-y_s\pi ]\hfill \\ +{(-1)}^{k_3}{U}_l^{m,s}\left(ϵ\right)cos[(l-2p+j)M-(l-2q+r)M_{☾}+\hfill \\ (l-2p)\omega -(l-2q){\omega }_{☾}+m\Omega -s({\Omega }_{☾}-\frac{\pi }{2})-y_s\pi \left]\right\},\hfill \end{array}$$

(2)

where ${\mu }_{☾}$ is the gravitational parameter of the Moon, $ϵ$ is the obliquity of the ecliptic, $y_s=0$ for s even while $y_s=\frac{1}{2}$ for s odd. The analytical expansions of the Kaula inclination functions and of $U_m^{l,\pm s}$, the Hansen coefficients and the coefficients ${ϵ}_m$, ${ϵ}_s$, $k_1$, $k_2$, $k_3$ can be found in [12,13,14].

We are interested in a model including the oblateness effect and the third-body perturbation up to the third order, i.e., including the harmonics in Equations (1) and (2) with $l=2,3$. In order to investigate the long-term evolution, a doubly averaged model is used, i.e., the disturbing potential is averaged over the mean motion of the satellite and over the mean motion of the third body. The disturbing potential $\mathcal{R}$ can be written as the sum of the contribution due to the second zonal harmonics, say ${\overline{\mathcal{R}}}_{J_2}$, and the one due to Moon and Sun, say ${\stackrel{=}{\mathcal{R}}}_{☾}$ and ${\stackrel{=}{\mathcal{R}}}_{\odot }$, respectively. We have

$$\mathcal{R}={\overline{\mathcal{R}}}_{J_2}+{\stackrel{=}{\mathcal{R}}}_{☾}+{\stackrel{=}{\mathcal{R}}}_{\odot },$$

(3)

where the first term in Equation (3) is the singly averaged secular oblateness effect [15]

$${\overline{\mathcal{R}}}_{J_2}=\frac{1}{4}{J}_2\frac{{\mu }_{\oplus }{R}_{\oplus }^2}{a^3{(1-e^2)}^{\frac{3}{2}}}\left(1-3{cos}^{2}i\right),$$

(4)

where $J_2$ is the second-order zonal coefficient of the Earth’s gravity field harmonic expansion, ${\mu }_{\oplus }$ is the Earth’s gravitational parameter and $R_{\oplus }$ represents the equatorial mean radius of the Earth. On the other hand, the terms ${\stackrel{=}{\mathcal{R}}}_{☾}$ and ${\stackrel{=}{\mathcal{R}}}_{\odot }$ in Equation (3) are, respectively, the lunar and solar doubly averaged disturbing functions.They can be obtained as the collections of all the terms in Equations (1) and (2), respectively, such that:

$$\left\{\begin{array}{c}l=2,3;\hfill \\ l-2p+j=0;\hfill \\ l-2q+r=0.\hfill \end{array}\right.$$

(5)

We recall that the averaging procedure is allowed whenever no mean motion resonance and no semi-secular lunisolar resonance occur, the latter arising from a commensurability between the slow frequencies of the satellite and the mean anomalies of Moon and Sun.

2.2. Hamiltonian Dynamics

In order to highlight the Hamiltonian structure of the problem, a coordinate change is required, by considering the Delaunay canonical variables $(L,G,H,\ell ,g,h)$ defined as [16]

$$\begin{array}{cccccc}\hfill & L=\sqrt{{\mu }_{\oplus }a},\hfill & \hfill & G=L\sqrt{1-e^2},\hfill & \hfill & H=Gcosi,\hfill \\ \hfill & \ell =M,\hfill & \hfill & g=\omega ,\hfill & \hfill & h=\Omega .\hfill \end{array}$$

(6)

In this way, the Hamiltonian describing the long-term lunisolar effect on a Molniya satellite can be written as

$$\mathcal{H}(L,G,H,\ell ,g,h)=H_{kep}\left(L\right)+{\mathcal{H}}_{J_2}(G,H;L)+{\mathcal{H}}_{☾}(G,H,g,h;L)+{\mathcal{H}}_{\oplus }(G,H,g,h;L),$$

(7)

where

$$H_{kep}=-\frac{{\mu }_{\oplus }^2}{2L^2},$$

(8)

and

$${\mathcal{H}}_{J_2}=-{\overline{\mathcal{R}}}_{J_2},{\mathcal{H}}_{☾}=-{\stackrel{=}{\mathcal{R}}}_{☾},{\mathcal{H}}_{\odot }=-{\stackrel{=}{\mathcal{R}}}_{\odot },$$

(9)

written in terms of Delaunay variables. It has to be pointed out that in Equation (2), the harmonic argument vanishes for $l=2$, $p=1$ and $m,s=0$. Since the corresponding term in ${\mathcal{H}}_{☾}$ only depends on the actions $(L,G,H)$, we will call this special harmonic the lunar mean term. As for the lunar case, for $l=2$, $p=1$ and $m=0$, the solar harmonic argument in Equation (1) disappears and thus we will refer to the corresponding harmonic term as the solar mean term.

Let us introduce the following notation:

$$\begin{array}{cc}\hfill & {\mathcal{H}}_{☾}(G,H,g,h;L)=C_0^{☾}{\mathcal{A}}_0^{☾}(G,H;L)+\sum _{\alpha }{C}_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}(G,H;L)cos\left({\phi }_{\alpha }^{☾}\right),\hfill \\ & {\mathcal{H}}_{\odot }(G,H,g,h;L)=C_0^{\odot }{\mathcal{A}}_0^{\odot }(G,H;L)+\sum _{\gamma }{C}_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }(G,H;L)cos\left({\varphi }_{\gamma }^{\odot }\right),\hfill \end{array}$$

(10)

where:

• $\alpha$ and $\gamma$ index the finite number of lunar and solar, respectively, harmonics retained in the model;
• $C_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}(G,H;L)$ is the $\alpha$-th lunar harmonic coefficient, and the constant term $C_{\alpha }^{☾}$ includes the lunar orbital parameters;
• $C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }(G,H;L)$ is the $\gamma$-th solar harmonic coefficient and $C_{\gamma }^{\odot }$ includes the solar orbital parameters;
• $\alpha =0$ and $\gamma =0$ denote the mean terms, i.e., $C_0^{☾}{\mathcal{A}}_0^{☾}(G,H;L)$ is the lunar mean term and $C_0^{\odot }{\mathcal{A}}_0^{\odot }(G,H;L)$ is the solar mean term;
• ${\phi }_{\alpha }^{☾}$ is the $\alpha$-th lunar argument and ${\varphi }_{\gamma }^{\odot }$ is the $\gamma$-th solar argument.

Generally, both sine and cosine trigonometric functions appear in the development of the lunar disturbing function. However, in our case, only cosine harmonics persist because of the values that s and $y_s$ in Equation (2) assume when $l=2,3$. The mean anomaly is cyclic in the doubly averaged Hamiltonian Equation (10); hence, the action L is a first integral, i.e., the semi-major axis is constant in the long term.

Given the above considerations, the dynamics of the $G,H,g$ and h are given by the following Hamilton equations:

$$\left\{\begin{array}{cc}\hfill \dot{G}& =\sum _{\alpha }\left[C_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}(G,H;L)\right]\frac{\partial {\phi }_{\alpha }^{☾}}{\partial g}sin\left({\phi }_{\alpha }^{☾}\right)+\sum _{\gamma }\left[C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }(G,H;L)\right]\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial g}sin\left({\varphi }_{\gamma }^{\odot }\right),\hfill \\ \hfill \dot{H}& =\sum _{\alpha }\left[C_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}(G,H;L)\right]\frac{\partial {\phi }_{\alpha }^{☾}}{\partial h}sin\left({\phi }_{\alpha }^{☾}\right)+\sum _{\gamma }\left[C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }(G,H;L)\right]\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial h}sin\left({\varphi }_{\gamma }^{\odot }\right),\hfill \\ \hfill \dot{g}& =\frac{\partial {\mathcal{H}}_{J_2}}{\partial G}(G,H;L)+C_0^{☾}\frac{\partial {\mathcal{A}}_0^{☾}(G,H;L)}{\partial G}+C_0^{\odot }\frac{\partial {\mathcal{A}}_0^{\odot }(G,H;L)}{\partial G}+\hfill \\ \hfill & \sum _{\alpha }\left[C_{\alpha }^{☾}\frac{\partial {\mathcal{A}}_{\alpha }^{☾}}{\partial G}(G,H;L)\right]cos\left({\phi }_{\alpha }^{☾}\right)+\sum _{\gamma }\left[C_{\gamma }^{\odot }\frac{\partial {\mathcal{A}}_{\gamma }^{\odot }}{\partial G}(G,H;L)\right]cos\left({\varphi }_{\gamma }^{\odot }\right),\hfill \\ \hfill \dot{h}& =\frac{\partial {\mathcal{H}}_{J_2}}{\partial H}(G,H;L)+C_0^{☾}\frac{\partial {\mathcal{A}}_0^{☾}(G,H;L)}{\partial H}+C_0^{\odot }\frac{\partial {\mathcal{A}}_0^{\odot }(G,H;L)}{\partial H}+\hfill \\ \hfill & \sum _{\alpha }\left[C_{\alpha }^{☾}\frac{\partial {\mathcal{A}}_{\alpha }^{☾}}{\partial H}(G,H;L)\right]cos\left({\phi }_{\alpha }^{☾}\right)+\sum _{\gamma }\left[C_{\gamma }^{\odot }\frac{\partial {\mathcal{A}}_{\gamma }^{\odot }}{\partial H}(G,H;L)\right]cos\left({\varphi }_{\gamma }^{\odot }\right).\hfill \end{array}\right.$$

(11)

Following, e.g., [15], the oblateness of the Earth does not produce any effect on the actions G and H, but it causes a precession, or regression, of g and h, which is usually used to approximate the evolution of the angles, as already mentioned before. From the last two equations of the system, Equation (11), we notice that the angles undergo secular drifts, caused by both the oblateness and the lunar and solar mean terms, and periodic effects, given by integrating the oscillating terms whose amplitude is proportional to the partial derivatives of the harmonic coefficients. Since the Laplace radius is around $7.7R_{\oplus }$ [17], higher than the Molniya one, and is the geocentric distance at which the order of magnitude of the precession caused by the lunisolar perturbation is equivalent to the one caused by the Earth’s oblateness, the following approximation

$$\left\{\begin{array}{c}\dot{g}\approx \frac{\partial {\mathcal{H}}_{J_2}}{\partial G},\hfill \\ \dot{h}\approx \frac{\partial {\mathcal{H}}_{J_2}}{\partial H},\hfill \end{array}\right.$$

(12)

is usually both convenient and accurate enough. However, in the particular case of the Molniya dynamics, the orbits are critically inclined; thus, the third-body perturbation might not be necessarily negligible at least for $\dot{g}$, as confirmed by numerical experiments that will be presented in Section 3.1.

From Equation (11), it is noteworthy that $\dot{G}$, and thus the mean eccentricity, only depends on the harmonics showing $\dot{g}$ in the argument. On the other hand, only harmonics depending on h affect $\dot{H}$; therefore, the orbital inclination depends on the harmonics having g or h in the argument.

Furthermore, the first two equations of the system (11) suggest that the larger the harmonic coefficient, the deeper the resulting fluctuations in $\dot{G}$ and $\dot{H}$. Following the approach presented in [10] for the dynamics in the main asteroid belt, we will analyze the role in the evolution of G and H also of the harmonic frequency.

Let us consider a first approximation of the system (11) such that:

• the actions are assumed as a constant, namely

  $${\mathcal{A}}_{\alpha }^{☾}(G,H;L)={\mathcal{A}}_{\alpha }^{☾},{\mathcal{A}}_{\gamma }^{\odot }(G,H;L)={\mathcal{A}}_{\gamma }^{\odot },$$
• the angles evolve linearly in time, namely

  $${\phi }_{\alpha }^{☾}\left(t\right)={\phi }_{\alpha ,0}^{☾}+{\dot{\phi }}_{\alpha }^{☾}t,{\varphi }_{\gamma }^{\odot }\left(t\right)={\varphi }_{\gamma ,0}^{\odot }+{\dot{\varphi }}_{\gamma }^{\odot }t,$$

  being ${\phi }_{\alpha ,0}^{☾}$ and ${\varphi }_{\gamma ,0}^{\odot }$ generic initial conditions and ${\dot{\phi }}_{\alpha }^{☾}$ and ${\dot{\varphi }}_{\gamma }^{\odot }$ constants.

Note that the last assumption holds especially in the neighborhood of a resonance.

In this way, from (11), we obtain

$$\left\{\begin{array}{cc}\hfill \dot{G}& =\sum _{\alpha }{C}_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}\frac{\partial {\phi }_{\alpha }^{☾}}{\partial g}sin({\phi }_{\alpha ,0}^{☾}+{\dot{\phi }}_{\alpha }^{☾}t)+\sum _{\gamma }{C}_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial g}sin({\varphi }_{\gamma ,0}^{\odot }+{\dot{\varphi }}_{\gamma }^{\odot }t),\hfill \\ \hfill \dot{H}& =\sum _{\alpha }{C}_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}\frac{\partial {\phi }_{\alpha }^{☾}}{\partial h}sin({\phi }_{\alpha ,0}^{☾}+{\dot{\phi }}_{\alpha }^{☾}t)+\sum _{\gamma }{C}_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial h}sin({\varphi }_{\gamma ,0}^{\odot }+{\dot{\varphi }}_{\gamma }^{\odot }t).\hfill \end{array}\right.$$

(13)

By integration, defining $\Delta G=G\left(T\right)-G\left(0\right)$ and $\Delta H=H\left(T\right)-H\left(0\right)$, we obtain

$$\begin{array}{cc}\hfill & \Delta G=\sum _{\alpha }\frac{C_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}}{{\dot{\phi }}_{\alpha }^{☾}}\frac{\partial {\phi }_{\alpha }^{☾}}{\partial g}\left[cos\left({\phi }_{\alpha ,0}^{☾}\right)-cos\left({\phi }_{\alpha }^{☾}\left(T\right)\right)\right]+\sum _{\gamma }\frac{C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }}{{\dot{\varphi }}_{\gamma }^{\odot }}\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial g}\left[cos\left({\varphi }_{\gamma ,0}^{\odot }\right)-cos\left({\varphi }_{\gamma }^{\odot }\left(T\right)\right)\right],\hfill \\ \hfill & \Delta H=\sum _{\alpha }\frac{C_{\alpha }^{☾}{\mathcal{A}}_{\alpha }^{☾}}{{\dot{\phi }}_{\alpha }^{☾}}\frac{\partial {\phi }_{\alpha }^{☾}}{\partial h}\left[cos\left({\phi }_{\alpha ,0}^{☾}\right)-cos\left({\phi }_{\alpha }^{☾}\left(T\right)\right)\right]+\sum _{\gamma }\frac{C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }}{{\dot{\varphi }}_{\gamma }^{\odot }}\frac{\partial {\varphi }_{\gamma }^{\odot }}{\partial h}\left[cos\left({\varphi }_{\gamma ,0}^{\odot }\right)-cos\left({\varphi }_{\gamma }^{\odot }\left(T\right)\right)\right].\hfill \end{array}$$

(14)

In other words, the long-term variation of G and H depends on the ratio between the amplitude of the harmonic coefficients and the corresponding frequency. In the following, this ratio will be denoted as R. The larger the ratio, the deeper the long-term effects. Roughly speaking, a large ratio R, that identifies a strong perturbing term, is due to a large amplitude or to a small frequency. In the first case, the term produces deep periodic oscillations of G and H; in the latter, we are dealing with a resonant or near-resonant term leading to significant variation over the time. Consideration of the ratio gives an indication of the perturbation when the distinction between the two types of dynamics is not so sharp. To identify the contributions that govern the dynamics of G and H, and thus of e and i, it is thereby interesting to evaluate the amplitudes, the periods and the ratio of all the terms involved in the model. An in-depth analysis of these terms will be presented in Section 3.1.

3. Numerical Evaluation of the Dominant Terms of the Doubly Averaged Lunisolar Perturbation

In the following, we will assume the values

$$a_{moln}=26,554.3\mathrm{km},e_{moln}=0.72,i_{moln}=63.{43}^{\circ }$$

and definitions

$$L_{moln}=\sqrt{{\mu }_{\oplus }{a}_{moln}},G_{moln}=L_{moln}\sqrt{1-e_{moln}^2},H_{moln}=G_{moln}cos{i}_{moln}.$$

We will refer to the above parameters as the Molniya parameters. For the sake of consistency, we also use the Delaunay angles for both the Moon and the Sun, namely

$$\left\{\begin{array}{c}{\omega }_{☾}=g_{☾}\hfill \\ {\Omega }_{☾}=h_{☾}\hfill \end{array}\right.,\left\{\begin{array}{c}{\omega }_{\odot }=g_{\odot }\hfill \\ {\Omega }_{\odot }=h_{\odot }\hfill \end{array}\right.$$

(15)

Important results will be translated in terms of Keplerian elements in order to be more understandable.

3.1. The Dominant Terms in the Long-Term Dynamics

According to the theoretical considerations exposed in Section 2.2 and taking as a reference system (11), we evaluate the amplitudes of the harmonic coefficients (see

Table 1. Largest amplitudes (in km²/s²), in absolute value, of the solar harmonic coefficients, together with the corresponding argument (left column). The values are computed by evaluating the harmonic coefficients at the Molniya parameters, i.e., $|C_{\gamma }^{\odot }{\mathcal{A}}_{\gamma }^{\odot }(L_{moln},G_{moln},H_{moln})|$.

${\mathit{\varphi }}_{\mathit{\gamma }}^{\odot }$, $\mathit{l}=2$	$|{\mathit{C}}_{\mathit{\gamma }}^{\odot }{\mathcal{A}}_{\mathit{\gamma }}^{\odot }|$
$2g$	$8.29\times {10}^{-6}$
$2g+(h-h_{\odot })$	$6.42\times {10}^{-6}$
$h-h_{\odot }$	$5.44\times {10}^{-6}$
$2g-(h-h_{\odot })$	$2.45\times {10}^{-6}$
Mean Term	$1.89\times {10}^{-6}$
$2(h-h_{\odot })$	$1.18\times {10}^{-6}$
$2g+2(h-h_{\odot })$	$1.13\times {10}^{-6}$
$2g-2(h-h_{\odot })$	$1.64\times {10}^{-7}$

and

Table 2. Largest amplitudes (in km²/s²), in absolute value, of the lunar harmonic coefficients, together with the corresponding argument. The values are computed by evaluating the harmonic coefficients at the Molniya parameters.

${\mathit{\phi }}_{\mathit{\alpha }}^{☾}$, $\mathit{l}=2$	$|{\mathit{C}}_{\mathit{\alpha }}^{☾}{\mathcal{A}}_{\mathit{\alpha }}^{☾}|$
$2g$	$1.79\times {10}^{-5}$
$2g+h$	$1.39\times {10}^{-5}$
h	$1.18\times {10}^{-5}$
$2g-h$	$5.30\times {10}^{-6}$
Mean Term	$4.09\times {10}^{-6}$
$2g+h-h_{☾}$	$2.75\times {10}^{-6}$
$2h$	$2.55\times {10}^{-6}$
$2g+2h$	$2.43\times {10}^{-6}$
$h-h_{☾}$	$2.33\times {10}^{-6}$
$2g+h_{☾}$	$1.16\times {10}^{-6}$
$2g-h_{☾}$	$1.16\times {10}^{-6}$
$2h-h_{☾}$	$1.11\times {10}^{-6}$
$2g+2h-h_{☾}$	$1.06\times {10}^{-6}$
$2g-h+h_{☾}$	$1.05\times {10}^{-6}$
$h_{☾}$	$5.31\times {10}^{-7}$
$2g+h+h_{☾}$	$4.02\times {10}^{-7}$
$2g-2h$	$3.55\times {10}^{-7}$
$h+h_{☾}$	$3.41\times {10}^{-7}$
$2g-2h+h_{☾}$	$1.55\times {10}^{-7}$
$2g-h-h_{☾}$	$1.54\times {10}^{-7}$
$2h-2h_{☾}$	$1.20\times {10}^{-7}$
$2g+2h-2h_{☾}$	$1.15\times {10}^{-7}$
$2g+h-2h_{☾}$	$5.90\times {10}^{-8}$
$h-2h_{☾}$	$5.00\times {10}^{-8}$
$2h+h_{☾}$	$4.78\times {10}^{-8}$
$2g+2h+h_{☾}$	$4.56\times {10}^{-8}$
$2g-h+2h_{☾}$	$2.25\times {10}^{-8}$
$2g-2h+2h_{☾}$	$1.68\times {10}^{-8}$
$2g-2h_{☾}$	$1.13\times {10}^{-8}$
$2g+2h_{☾}$	$1.13\times {10}^{-8}$
$3g-g_{☾}+h-h_{☾}$	$5.92\times {10}^{-8}$
$g+g_{☾}-h+h_{☾}$	$5.60\times {10}^{-8}$
$g-g_{☾}+h-h_{☾}$	$5.60\times {10}^{-8}$
$3g-g_{☾}+2h-h_{☾}$	$5.45\times {10}^{-8}$
$g+g_{☾}-2h+h_{☾}$	$5.15\times {10}^{-8}$
$3g+g_{☾}+h_{☾}$	$3.97\times {10}^{-8}$
$3g-g_{☾}-h_{☾}$	$3.97\times {10}^{-8}$
$3g+g_{☾}-h+h_{☾}$	$2.26\times {10}^{-8}$
$3g+g_{☾}+h+h_{☾}$	$2.16\times {10}^{-8}$
$g-g_{☾}-h-h_{☾}$	$2.05\times {10}^{-8}$
$g+g_{☾}+h+h_{☾}$	$2.05\times {10}^{-8}$
$g-g_{☾}+2h-h_{☾}$	$1.97\times {10}^{-8}$
$g-g_{☾}+3h-h_{☾}$	$1.75\times {10}^{-8}$
$3g-g_{☾}+2h-2h_{☾}$	$1.03\times {10}^{-8}$
$3g-g_{☾}+3h-h_{☾}$	$1.00\times {10}^{-8}$

), their partial derivatives with respect to the actions (see

Table 3. Largest amplitudes (in rad/s), in absolute value, of the partial derivatives of the harmonic coefficients with respect to the actions, together with the corresponding argument. The values are computed by evaluating the terms at the Molniya parameters.

${\mathit{\phi }}_{\mathit{\alpha }}^{☾}$	$\left|{\mathit{C}}_{\mathit{\alpha }}^{☾}\frac{\mathit{\partial }{\mathcal{A}}_{\mathit{\alpha }}^{☾}}{\mathit{\partial }\mathit{G}}\right|$					${\mathit{\phi }}_{\mathit{\alpha }}^{☾}$	$\left|{\mathit{C}}_{\mathit{\alpha }}^{☾}\frac{\mathit{\partial }{\mathcal{A}}_{\mathit{\alpha }}^{☾}}{\mathit{\partial }\mathit{H}}\right|$
Mean Term	$1.25\times {10}^{-10}$					Mean Term	$3.85\times {10}^{-10}$
$2g+h$	$3.72\times {10}^{-10}$					$2g$	$2.80\times {10}^{-10}$
$2g$	$3.41\times {10}^{-10}$					h	$2.76\times {10}^{-10}$
h	$2.57\times {10}^{-10}$					$2g-h$	$1.76\times {10}^{-10}$
${\mathbf{\varphi }}_{\mathbf{\gamma }}^{\odot }$	$\left|{\mathbf{C}}_{\mathbf{\gamma }}^{\odot }\frac{\partial {\mathcal{A}}_{\mathbf{\gamma }}^{\odot }}{\partial \mathbf{G}}\right|$					${\mathbf{\varphi }}_{\mathbf{\gamma }}^{\odot }$	$\left|{\mathbf{C}}_{\mathbf{\gamma }}^{\odot }\frac{\partial {\mathcal{A}}_{\mathbf{\gamma }}^{\odot }}{\partial \mathbf{H}}\right|$
$2g+h-h_{\odot }$	$1.72\times {10}^{-10}$					Mean Term	$1.78\times {10}^{-10}$
$2g$	$1.58\times {10}^{-10}$					$2g$	$1.30\times {10}^{-10}$
$h-h_{\odot }$	$1.19\times {10}^{-10}$					$h-h_{\odot }$	$1.28\times {10}^{-10}$
Mean Term	$5.81\times {10}^{-11}$

) and the periods of the harmonic arguments (see Tables 5 and 6) at the Molniya parameters. Since the functions involved are sufficiently regular, the results provide an accurate estimate of the perturbing terms affecting a satellite in a Molniya regime.

Table 1 shows all the amplitudes of the second-order solar harmonics. The number of the third-order solar harmonics computed is 28, but the corresponding coefficients are too small to be displayed: the largest values are of the order of ${10}^{-11}$ km²/s². In the lunar case, the number of the second-order harmonics evaluated is 38, ranging from values of approximately ${10}^{-5}$ km²/s² to ${10}^{-11}$ km²/s². Instead, the third-order contribution consists of 196 harmonics, ranging from approximately ${10}^{-8}$ km²/s² to ${10}^{-17}$ km²/s². Of all these estimates, the largest amplitudes of both the second- and the third-order lunar potential are listed in Table 2.

We notice that the largest lunar amplitudes in Table 2 are only one order of magnitude lower than the $H_{J_2}$ contribution in Equation (7), which is of order ${10}^{-4}$ km²/s² when evaluated at the Molniya parameter. This is due to the fact that a Molniya satellite attains high altitudes.

Moreover, from the table, it is clear that, although considering the geocentric orbits of the Moon and of the Sun with their actual eccentricity ($e_{☾}=0.0549$, $e_{\odot }=0.0167$), the third-order contribution in the Hamiltonian given by a third body is quite small, even if Molniya are highly eccentric and high-altitude orbits.

In Table 3, we give the estimates of the largest values of the partial derivatives of the harmonic coefficients with respect to the actions, computed in the Molniya region. This analysis would indicate the dominant terms in the angular dynamics defined by the last two equations in Equation (11). Table 3 on the left shows the lunisolar contribution to $\dot{g}$, while, on the right, we report the terms determining the dynamics of h. Only the second-order lunisolar contribution was taken into account.

The precession caused by the lunar and solar mean terms (Table 3 on the left) is certainly larger than the secular drift of the argument of perigee due to the Earth’s oblateness, which cancels out, by definition, at the critical inclination.

Moreover, the amplitudes of the oscillations seem to be quite large. These facts imply that the third-body effects on the dynamics of the argument of perigee are small but not negligible if compared with the oblateness effect. If we change slightly the inclination—for example, $i\approx 63.43\pm 2^{\circ }$—the oblateness drift becomes only one order of magnitude larger than the ones produced by the third-body mean terms. For practical applications, because of the unavoidable errors in spacecraft maneuvering, the conclusion is that the gravitational attraction exerted by the Moon and the Sun is essential to study in depth the periods and the perigee resonance that characterize Molniya-like orbits (Section 4.1 and Section 4).

On the contrary, the partial derivatives of both the lunar and the solar mean terms in Table 3 on the right ensure that the oblateness effect is still the dominant perturbation for $\dot{h}$, as it usually happens in the case of a no frozen condition. Indeed,

$${\dot{h}}_{J_2}=\frac{\partial {\mathcal{H}}_{J_2}}{\partial H}{|}_{(L_{moln},G_{moln},H_{moln})}\approx -2.63\times {10}^{-8}.$$

(16)

Hence, to better capture how the lunisolar perturbation might affect the dynamics of the angular terms, the periods of the arguments involved in the doubly averaged model are computed assuming $\dot{h}\approx {\dot{h}}_{J_2}$ and different approximations of $\dot{g}$ taking into account also the the lunisolar terms listed in Table 3.

To handle the periodic effects, we need to set initial conditions also for the argument of pericenter and for the longitude of the ascending node of the satellite. An initial argument of perigee at ${270}^{\circ }$ is the best initial condition to allow the largest coverage of the Russian territory by a Molniya satellite [3,4]; thus, we focus only on how lunisolar effects on the angles vary with respect to the initial ascending node of the satellite. This choice is dictated by the fact highlighted, for instance, by Anselmo and Pardini in [18]: the initial ascending node is crucial for the satellite’s lifetime. Moreover, we can always fix $h_{\odot }=0$. The notation adopted to represent the different approximations of $\dot{g}$ and the details are arranged in

Table 4. Notation and details of the various approximations of $\dot{g}$ considered in the work. The first column is the notation; the second to the fourth columns show whether a perturbing effect is taken into account (✓) or not (−). The last column reports, where relevant, the fixed value (in degrees) of the initial longitude of the ascending node $h_{IC}$ chosen to compute the approximation of $\dot{g}$.

Notation	${\mathbf{J}}_2$-Term	Mean Terms	Periodic Terms in Table 3	${\mathit{h}}_{\mathbf{IC}}$
${\dot{g}}_{J_2}$	✓	−	−	Not relevant
${\dot{g}}_0$	✓	✓	−	Not relevant
${\dot{g}}_1$	✓	✓	✓	90
${\dot{g}}_2$	✓	✓	✓	180
${\dot{g}}_3$	✓	✓	✓	0
${\dot{g}}_4$	✓	✓	✓	270

.

In

Table 5. Largest periods (>7 years) for the lunar and solar arguments appearing in the second-order doubly averaged disturbing potential. The first column helps to identify which third body the arguments belong to. The periods (in years) are computed assuming $\dot{h}\approx {\dot{h}}_{J_2}$ and following the assumptions listed in Table 4 for $\dot{g}$, as specified at the top of each column.

	Argument	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_{{\mathit{J}}_2}$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_0$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_1$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_2$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_3$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_4$
$☾,\odot$	$2g$	9777.54	513.68	196.72	163.76	150.31	130.27
☾	$2h_{☾}-h-2g$	40.25	43.47	50.34	53.07	54.65	57.89
☾	$2h_{☾}-h$	40.08	40.08	40.08	40.08	40.08	40.08
☾	$2h_{☾}-h+2g$	39.92	37.18	33.30	32.20	31.64	30.65
☾	$h_{☾}+2g$	18.65	19.31	20.56	21.00	21.24	21.71
☾	$h_{☾}$	18.61	18.61	18.61	18.61	18.61	18.61
☾	$h_{☾}-2g$	18.58	17.96	17.00	16.71	16.56	16.28
☾	$h_{☾}-h-2g$	12.73	13.03	13.58	13.78	13.88	14.08
☾	$h_{☾}-h$	12.71	12.71	12.71	12.71	12.71	12.71
☾	$h_{☾}-h+2g$	12.69	12.40	11.94	11.79	11.72	11.58
☾	$2h_{☾}+2g$	9.31	9.47	9.77	9.87	9.92	10.02
☾	$2h_{☾}$	9.31	9.31	9.31	9.31	9.31	9.31
☾	$2h_{☾}-2g$	9.30	9.14	8.89	8.81	8.76	8.68
$☾,\odot$	$2g+h$	7.56	7.66	7.85	7.92	7.95	8.01
$☾,\odot$	h	7.55	7.55	7.55	7.55	7.55	7.55
$☾,\odot$	$2g-h$	7.55	7.55	7.27	7.21	7.19	7.13

, we list the largest second-order periods obtained, while

Table 6. Largest periods for the lunar and solar arguments appearing in the third-order doubly averaged disturbing potential. The first column helps to identify the third body. The periods (in years) are computed assuming the approximations of $\dot{g}$ listed in Table 4 and $\dot{h}={\dot{h}}_{J_2}$.

	Argument	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_{{\mathit{J}}_2}$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_0$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_1$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_2$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_3$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_4$
⊙	$g-g_{\odot }$	23,669.36	1036.82	394.83	328.48	301.43	261.16
⊙	$g-g_{\odot }$	16,659.31	1018.06	392.08	326.57	299.83	259.95
⊙	$3g-g_{\odot }$	6919.27	343.50	131.30	109.28	100.30	86.91
⊙	$3g+g_{\odot }$	6161.37	341.41	131.00	109.07	100.12	86.78
☾	$3g-g_{☾}-3h_{☾}$	184.42	367.55	488.04	279.03	227.11	168.40
☾	$g-g_{☾}-3h_{☾}$	181.00	217.27	329.57	396.41	444.55	575.42
☾	$g+g_{☾}+3h_{☾}$	177.71	152.69	123.19	115.89	112.33	106.23
☾	$3g+g_{☾}+3h_{☾}$	174.54	117.70	75.75	67.86	64.29	58.51
☾	$3g-g_{☾}-2h+2h_{☾}$	108.11	154.24	562.28	4104.67	1736.91	473.80
☾	$g-g_{☾}-2h+2h_{☾}$	106.93	118.62	145.73	157.48	164.55	179.68
☾	$g+g_{☾}+2h-2h_{☾}$	105.77	96.37	83.72	80.28	78.55	75.52
☾	$3g+g_{☾}+2h-2h_{☾}$	104.64	81.15	58.72	53.87	51.59	47.80
☾	$3g+g_{☾}+h+h_{☾}$	52.03	60.78	85.12	97.91	106.45	127.23
☾	$g+g_{☾}+h+h_{☾}$	51.75	54.35	59.41	61.27	62.32	64.37
☾	$g-g_{☾}-h-h_{☾}$	51.48	49.15	45.63	44.59	44.05	43.08
☾	$3g-g_{☾}-h-h_{☾}$	51.21	44.86	34.07	37.03	35.04	32.37

shows the third-order ones. We remark that the arguments that do not depend on $\dot{g}$ are also shown for the sake of completeness. As we will see in the second part of the work, to each resonant term including g, one can associate an integral of motion that comprises the value of the semi-major axis, eccentricity and inclination. Thus, the terms in the disturbing potential that explicitly cause the inclination to evolve implicitly change the long-term picture when combined with the eccentricity resonances.

In both tables, the arguments are grouped with respect to the associated periods to make the reading easier. In Table 5, macro periods of the order 40.08 years, 18.61 years, 12.71 years, 9.30 years and 7.55 years are highlighted, consistent with the frequency analysis [2]. The argument $2g$ represents the main resonant angle, because of the critical inclination, and thus we compute a very large period, for any assumptions we take for $\dot{g}$. Moreover, we find the well-known value 18.61 years in correspondence with the period of the lunar ascending node. Small periods are related to high frequencies that are not very sensitive to the value of $\dot{g}$. On the contrary, the largest periods strongly depend on the approximation chosen. The same feature also emerges from Table 6, where the groups are of four arguments in the lunar case and of two arguments in the solar case. The solar third-order critical arguments $g\pm g_{\odot }$ and $3g\pm g_{\odot }$ (top of Table 6) behave as the main resonant angle. Conversely, the third-order lunar arguments $3g-g_{☾}-2h+2h_{☾}$ and $3g+g_{☾}+h+h_{☾}$ behave in the opposite way: by increasing the lunisolar perturbation, the arguments may even become critical.

In

Table 7. Lunisolar harmonics with $R>1$ km²/s. The first column helps to identify the third body; the second one indicates whether the corresponding argument appears in the second- or in the third-order lunisolar potential. The ratioR (in km²/s) is computed by assuming the approximations of $\dot{g}$ listed in Table 4 and $\dot{h}={\dot{h}}_{J_2}$.

	l	Argument	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_{{\mathit{J}}_2}$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_0$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_1$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_2$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_3$	$\dot{\mathit{g}}\approx {\dot{\mathit{g}}}_4$
☾	2	$2g$	879,496.40	46,205.55	17,695.51	14,730.36	13,520.88	11,718.5
⊙	2	$2g$	407,137.87	21,389.55	8191.64	6819.00	6259.11	5424.75
☾	2	$2g+h$	526.48	533.92	547.08	551.51	553.90	558.45
☾	2	h	446.00	446.00	446.004	446.00	446.00	446.00
⊙	2	$2g+h$	243.72	247.16	253.25	255.30	256.41	258.52
⊙	2	h	206.46	206.46	206.46	206.46	206.46	206.46
☾	2	$2g-h$	200.75	197.99	193.47	192.05	191.29	189.89
☾	2	$2g+h-h_{☾}$	175.79	180.01	187.69	190.33	191.78	194.54
☾	2	$h-h_{☾}$	148.84	148.84	148.84	148.84	148.84	148.84
☾	2	$2g+h_{☾}$	108.85	112.73	120.00	122.58	124.00	126.75
☾	2	$2g-h_{☾}$	108.44	104.85	99.25	97.56	96.67	95.06
⊙	2	$2g-h$	92.93	91.65	89.56	88.90	88.55	87.90
☾	2	$2g-h+h_{☾}$	66.96	65.43	62.97	62.22	61.82	61.08
☾	2	$h_{☾}$	49.65	49.65	49.65	49.65	49.65	49.65
☾	2	$2h$	48.33	48.33	48.33	48.33	48.33	48.33
☾	2	$2(g+h)$	46.15	46.48	47.04	47.22	47.32	47.51
☾	2	$2h-h_{☾}$	26.42	26.42	26.42	26.42	26.42	26.42
☾	2	$2(g+h)-h_{☾}$	25.23	25.46	25.84	25.97	26.04	26.17
⊙	2	$2h$	22.37	22.3	22.37	22.37	22.37	22.37
⊙	2	$2(g+h)$	21.36	21.51	21.77	21.86	21.91	21.99
☾	2	$2h_{☾}-2g-h$	11.92	12.88	14.91	15.72	16.19	17.15
☾	2	$h_{☾}+2g+h$	10.85	10.96	11.15	11.21	11.24	11.31
☾	2	$2h_{☾}-h$	10.07	10.07	10.07	10.07	10.07	10.07
☾	2	$h_{☾}+h$	9.19	9.19	9.19	9.19	9.19	9.19
☾	2	$2g-2h$	6.72	6.68	6.60	6.58	6.56	6.54
☾	3	$3g+g_{☾}+h+h_{☾}$	5.65	6.60	9.24	10.63	11.56	13.82
☾	3	$g+g_{☾}+h+h_{☾}$	5.32	5.58	6.10	6.29	6.40	6.61
☾	3	$g-g_{☾}-h-h_{☾}$	5.29	5.05	4.69	4.58	4.53	4.43
☾	2	$2h_{☾}+2g-h$	4.51	4.21	3.77	3.64	3.58	3.47
☾	2	$h_{☾}-2g+h$	4.14	4.10	4.03	4.01	4.00	3.98
☾	2	$2h_{☾}-2h$	3.85	3.85	3.85	3.85	3.85	3.85
☾	2	$h_{☾}+2g-2h$	3.68	3.64	3.59	3.57	3.57	3.56
☾	2	$2(h_{☾}-g-h)$	3.67	3.72	3.79	3.82	3.83	3.86
⊙	2	$2(g-h)$	3.11	3.09	3.06	3.04	3.04	3.03
☾	3	$3g-g_{☾}-h-h_{☾}$	2.12	1.86	1.53	1.45	1.41	1.34
☾	3	$3g-g_{☾}-h_{☾}$	1.77	1.81	1.89	1.92	1.94	1.97
☾	3	$3g+g_{☾}+h_{☾}$	1.76	1.18	1.65	1.63	1.62	1.60
☾	3	$3g+g_{☾}+h$	1.27	1.18	1.05	1.01	0.99	0.96
☾	3	$g-g_{☾}-h$	1.21	1.25	1.31	1.33	1.34	1.36
☾	3	$3g-g_{☾}+h-h_{☾}$	1.21	1.23	1.25	1.26	1.26	1.27
☾	3	$g+g_{☾}+h$	1.21	1.18	1.13	1.11	1.10	1.09
☾	3	$g-g_{☾}+h-h_{☾}$	1.15	1.15	1.16	1.16	1.16	1.16
☾	3	$g+g_{☾}-h+h_{☾}$	1.15	1.14	1.13	1.13	1.13	1.13

, all the harmonics whose $R>1$ km²/s are shown, from the highest to the lowest, taking as reference the assumption of the fourth column. At the top of the table, we find clearly the main resonant argument $2g$, which is the lunisolar term dominating the dynamics of the eccentricity for Molniya orbits, as is well known.

Subsequently, the harmonics corresponding to $2g\pm h$ and h show $R>200$ km²/s; thus, they produce a significant contribution to the dynamics; precisely, $2g\pm h$ affects both e and i, while h only perturbs i. These arguments are far from being critical; the periods are around 7.55 years, but their amplitudes are quite large (see Table 1 and Table 2).

In [5], the role played by $2g\pm h$ was analyzed by means of a double resonance model, even though the resonance $2g+h$ does not produce any resonances overlapping. Moreover, in [2], the authors proved that the terms corresponding to $2g$, $2g\pm h$ and h capture the main features of the long-term eccentricity variation, by comparing the pseudo-observation data given by the Two-Line Element (TLE) with ad hoc numerical propagations. Following the results in [2], in this work, the dominant terms are chosen as the ones associated with $R>100$ km²/s. Indeed, as long as the hypothesis of a constant semi-major axis holds true, a model including the lunar harmonics with $l=2$, $m=0,1$ and $s=0,1$ in Equation (2) and the solar harmonics with $l=2$, $m=0,1$ in Equation (1) seems to be the “poorest” model capable of representing the third-body effect in the actual mean eccentricity of many Molniya spacecraft. Such terms of the quadrupolar expansion occupy the first 14 lines of Table 7, having only 2 of them $R<100$ km²/s significantly.

4. The Phase Space Structure of the Main Terms

As an example of a possible application of the analysis carried out so far, we show the phase space corresponding to the dominant terms that have been identified in the previous section. As mentioned above, we focus on the terms whose ratio R, as lunar or solar harmonics, is larger than 100 km${}^2/$s, except for the contribution corresponding to h and $h-h_{☾}$, and we aim to identify the neighborhood in eccentricity and inclination where they play a role for the long-term dynamics. As before, for the sake of completeness, we introduce the investigation by recalling the basic concepts.

4.1. The Resonant Dynamics

A non-autonomous dynamical system can be converted into an autonomous one by adding one dimension to the phase space. Therefore, without loss of generality, in what follows, we assume to have an autonomous 3-degree-of-freedom nearly integrable Hamiltonian system. Referring to the classical theory presented in [16], let us take into account as a concrete example, useful to our purpose, the resonant Hamiltonian:

$${\mathcal{H}}_{res}(\mathbf{I},\psi )={\mathcal{H}}_0\left(\mathbf{I}\right)+\epsilon f\left(\mathbf{I}\right)cos(\mathbf{k}·\psi ),$$

(17)

where $\epsilon$ is the small parameter, $\mathbf{I}\in {\mathbb{R}}^3$ and $\psi \in {\mathbb{T}}^3$ are the action-angle variables 1 for the unperturbed Hamiltonian ${\mathcal{H}}_0$, and $\mathbf{k}·\dot{\psi }=0$ in some region of the phase space. The Hamilton equations arising from ${\mathcal{H}}_{res}$ are:

$$\left\{\begin{array}{c}\dot{\mathbf{I}}=\epsilon \mathbf{k}·f\left(\mathbf{I}\right)sin(\mathbf{k}·\psi ),\hfill \\ \dot{\psi }=\frac{\partial {\mathcal{H}}_0}{\partial \mathbf{I}}\left(\mathbf{I}\right)+\epsilon \frac{\partial f}{\partial \mathbf{I}}\left(\mathbf{I}\right)cos(\mathbf{k}·\psi ),\hfill \end{array}\right.$$

(18)

where $\frac{\partial {\mathcal{H}}_0}{\partial \mathbf{I}}\left(\mathbf{I}\right)$ is the vector of the main frequencies. From the classical theory, by resonance, we mean a commensurability between the main frequencies for some value of $\mathbf{I}={\mathbf{I}}^{*}$. In this case:

$$\mathbf{k}·\frac{\partial {\mathcal{H}}_0}{\partial \mathbf{I}}\left({\mathbf{I}}^{*}\right)=0.$$

(19)

We refer to the relation in Equation (19) as an exact resonance, while we talk of real resonance, or simply resonance, by referring to the following relation:

$$\mathbf{k}·\dot{\psi }\left(\mathbf{I}\right)=\mathbf{k}·\frac{\partial {\mathcal{H}}_0}{\partial \mathbf{I}}\left(\mathbf{I}\right)+\epsilon \mathbf{k}·\frac{\partial f}{\partial \mathbf{I}}\left(\mathbf{I}\right)cos(\mathbf{k}·\psi )=0.$$

(20)

If the perturbation is sufficiently small with respect to the unperturbed dynamics, then $\mathbf{k}·\dot{\psi }\left({\mathbf{I}}^{*}\right)\approx 0$ and the exact resonance may well approximate the real resonance at least up to the first order in $\epsilon$. There always exists a canonical transformation $\Phi$ such that the critical argument $\mathbf{k}·\psi$ is a new angle, i.e.,

$$(\mathbf{I},\psi )\stackrel{\Phi }{\mapsto }(\mathbf{J},\theta ),{\theta }_1={\mathbf{e}}_{\mathbf{1}}^T·\theta =\mathbf{k}·\psi .$$

(21)

After performing a coordinate change $\Phi$, the new one-degree-of-freedom (1-DOF) Hamiltonian

$${\mathcal{H}}_{res}^{\prime }(J_1,{\theta }_1)={\mathcal{H}}_0^{\prime }\left(\mathbf{J}\right)+\epsilon f^{\prime }\left(\mathbf{J}\right)cos{\theta }_1,$$

(22)

describes a two-dimensional motion taking place along the level curves $J_2=J_2^{*}$ and $J_3=J_3^{*}$ in the $(J_1,{\theta }_1)$ plane, where $J_1$ is the action conjugate to the critical angle ${\theta }_1$ and ${\mathbf{J}}^{*}=\Phi \left({\mathbf{I}}^{*}\right)$.

According to the Standard Resonance Model (SRM), also known as integrable approximation in [16], the Hamiltonian ${\mathcal{H}}_{res}^{\prime }$ can be developed in Taylor series of $J_1$ around $J_1^{*}$, which, as a first approximation, represents a pendulum-like dynamics in the neighborhood of the exact resonance.

Following [16], the resonant region is the libration region around the elliptic equilibrium point and its maximum libration width measured at the apex of the separatrix is given by

$$|J_1-J_1^{*}|\le 2\sqrt{\left|\frac{c}{\beta }\right|}.$$

(23)

If there are two or more resonances, then we can separately study the dynamics corresponding to each one, making the assumption that they are isolated. The resulting motion is pendulum-like with an appropriate coordinate change for every single resonance and the pendulum-like model may give a good approximation as long as the libration regions remain isolated. If any resonance overlap occurs, then the pendulum-like model breaks down. The separatrices of different resonances are connected if two or more resonances overlap; therefore, an initial condition in this region may produce jumps from one libration region to one other, showing chaotic diffusion [16].

Another scenario in which the classical approach does not provide a reliable description of the real resonant dynamics occurs when the exact resonance is not a good approximation of the real resonance. The real equilibria arising from the suitable Hamiltonian ${\mathcal{H}}_{res}^{\prime }$ are solutions of:

$$\left\{\begin{array}{c}{\dot{J}}_1=\epsilon f^{\prime }\left(\mathbf{J}\right)sin{\theta }_1=0,\hfill \\ {\dot{\theta }}_1=\frac{\partial {\mathcal{H}}_0^{\prime }}{\partial J_1}\left(\mathbf{J}\right)+\epsilon \frac{\partial f^{\prime }}{\partial J_1}\left(\mathbf{J}\right)cos{\theta }_1=0.\hfill \end{array}\right.$$

(24)

As in the pendulum case, the condition ${\dot{J}}_1=0$ implies ${\theta }_1=n\pi$ with $n\in \mathbb{Z}$. By replacing this solution in the second equation of (24), then the latter splits into two different equations, namely

$$\left\{\begin{array}{c}\frac{\partial {\mathcal{H}}_0^{\prime }}{\partial J_1}\left(\mathbf{J}\right)+\epsilon \frac{\partial f^{\prime }}{\partial J_1}\left(\mathbf{J}\right)=0,\hfill \\ \frac{\partial {\mathcal{H}}_0^{\prime }}{\partial J_1}\left(\mathbf{J}\right)-\epsilon \frac{\partial f^{\prime }}{\partial J_1}\left(\mathbf{J}\right)=0.\hfill \end{array}\right.$$

(25)

Hence, the solutions of (25) are not necessarily the same: the stronger the perturbing effects, controlled by $\epsilon$, the more the solutions of the system (25) are separated in the phase space. In such a case, the Taylor approximation, which characterizes the classical approach, fails to capture a deep asymmetry.

In the case of deep asymmetry, the maximum libration width can be computed by taking advantage of the invariant condition associated with the Hamiltonian function, namely

$${\mathcal{H}}_{res}^{\prime }(J_1,{\theta }_{1s})={\mathcal{H}}_{res}^{\prime }(J_{1u},{\theta }_{1u}),$$

(26)

where $(J_{1s},{\theta }_{1s})$ is the stable equilibrium of the system and $(J_{1u},{\theta }_{1u})$ is the unstable one, such that $J_{1s}\ne J_{1u}$. In what follows, we will refer to this as the non-standard approach (NSA), only to distinguish it from the pendulum-like case.

4.2. Single Resonance Hypothesis Description

If we apply the above description to the Hamiltonian Equation (7), the actions are $\mathbf{I}=(L,G,H)$, the angles $\mathit{\psi }=(\ell ,g,h)$ and the unperturbed term ${\mathcal{H}}_0$ is given by the $J_2$-term and both the lunar and solar mean terms, namely

$${\mathcal{H}}_0={\mathcal{H}}_{J_2}(G,H;L)+C_0^{☾}{\mathcal{A}}_0^{☾}(G,H;L)+C_0^{\odot }{\mathcal{A}}_0^{\odot }(G,H;L).$$

(27)

The resonant perturbation is given by the Hamiltonian contribution of a lunisolar dominant harmonic related to the resonant angle, considered on a case-by-case basis. This is to say that, in order to start from a simplified analysis, we assume that only one lunisolar term affects the dynamics at a time. Under the hypothesis of isolated resonance, the dynamics in a small enough neighborhood are described by a resonant Hamiltonian of the form given in Equation (17).

For each case that we will show, after performing a coordinate change of the form shown in Equation (21), the motion evolves in the $(J_1,{\theta }_1)$ plane. The first integral of motion is Kozai-like, i.e., the evolution of e and i is coupled (recalling that the semi-major axis is constant in our assumptions). Thus, the phase portraits are computed assuming $L=L_{moln}$ and they correspond to the given integral of motion computed at the Molniya parameters.

The resonant dominant harmonics for which the SRM gives a reliable description of the phase plane structure are listed in

Table 8. Resonances whose dynamics are well described by the SRM. The first column identifies the resonances through the critical argument associated with it; the second column shows the first integral arising from the resonant Hamiltonian. $\Gamma$ is the dummy momentum introduced to add one dimension to the phase space in case of resonances involving the lunar ascending node. The center of libration is given in terms of Keplerian elements: $(e^{*},i^{*})$ identifies the libration center of the exact resonance while $(e_s,i_s)$ gives the center of libration related to the real resonance. $|J_{1s}-J_{1u}|$ (in km²/s) gives information about the asymmetry of the resonant region. The last two columns show the libration width in terms of e and i; the same values are obtained with the SRM Equation (23) and with NSA Equation (26). The values in inclination are reported in degrees.

Critical Argument	First Integral	${\mathit{e}}^{*}$	${\mathit{e}}_{\mathit{s}}$	${\mathit{i}}^{*}$	${\mathit{i}}_{\mathit{s}}$	$|{\mathit{J}}_{1\mathit{s}}-{\mathit{J}}_{1\mathit{u}}|$	$\mathbf{\Delta }\mathit{e}$	$\mathbf{\Delta }\mathit{i}$
$2g-h$	$\sqrt{a(1-e^2)}(cosi+\frac{1}{2})$	0.64	0.64	69.14	69.03	114.41	0.13	7.3
$2g+h$	$\sqrt{a(1-e^2)}(\frac{1}{2}-cosi)$	0.98	0.98	56.06	56.06	1.20	0.004	0.55
$2g+h-h_{☾}$	$\begin{array}{cc}& \Gamma -\frac{1}{2}\sqrt{a(1-e^2)}\hfill \\ & \sqrt{a(1-e^2)}cosi\hfill \end{array}$	0.52	0.52	62.79	62.80	337.35	0.15	0.29
$2g+h_{☾}$	$\begin{array}{cc}& \Gamma +\frac{1}{2}\sqrt{a(1-e^2)}\hfill \\ & \sqrt{a(1-e^2)}cosi\hfill \end{array}$	0.76	0.76	61.56	61.55	10.27	0.03	1.63

, while the resonances associated with the arguments in

Table 9. Resonances whose dynamics are not appropriately described by the SRM. The first column identifies the resonance through the associated argument; the second column shows the first integral arising from the resonant Hamiltonian. The equilibria are given in terms of eccentricity and inclination: $(e_s,i_s)$ stable ones, $(e_u,i_u)$ unstable ones. $|J_{1s}-J_{1u}|$ (in km²/s) gives information about the asymmetry of the resonant region. The last two columns show the maximum excursion, in terms of e and i, that may be attained in the libration region as computed with NSA through Equation (26). The values in inclination are reported in degrees. N.E. means “not evaluated”.

Critical Argument	First Integral	${\mathit{e}}_{\mathit{s}}$	${\mathit{e}}_{\mathit{u}}$	${\mathit{i}}_{\mathit{s}}$	${\mathit{i}}_{\mathit{u}}$	$|{\mathit{J}}_{1\mathit{s}}-{\mathit{J}}_{1\mathit{u}}|$	$[{\mathit{e}}_{\mathbf{min}},{\mathit{e}}_{\mathbf{max}}]$	$[{\mathit{i}}_{\mathbf{min}},{\mathit{i}}_{\mathbf{max}}]$
$2g$	$\sqrt{a(1-e^2)}cosi$	0.72	0.71	63.29	63.69	625.10	$[0.55,0.79]$	$[59.30,68.11]$
$2g-h_{☾}$	$\begin{array}{cc}& \sqrt{a(1-e^2)}cosi\hfill \\ & \Gamma +\frac{1}{2}\sqrt{a(1-e^2)}\hfill \end{array}$	$\begin{array}{cc}& 0.62\hfill \\ & 0.57\hfill \end{array}$	$\begin{array}{cc}& 0.60\hfill \\ & 0.54\hfill \end{array}$	$\begin{array}{cc}& 67.97\hfill \\ & 68.95\hfill \end{array}$	$\begin{array}{cc}& 68.48\hfill \\ & 69.57\hfill \end{array}$	N.E.	$[0.44,0.68]$	$[66.28,70.83]$

exhibit a non-standard behavior. In Table 8, we show the center of libration related to each resonant harmonic and the corresponding maximum width, as computed with the standard approach (SRM) through Equation (23). The maximum real excursion in eccentricity and in inclination, computed by using Equation (26), gives substantially the same width obtained with the SRM. These facts can be appreciated by looking at the phase portraits from Figure 1, Figure 2, Figure 3 and Figure 4. The dynamical structure arising from the pendulum-like Hamiltonian is depicted on the left, while, on the right, we show the results obtained from the resonant Hamiltonian not developed in Taylor series; the Y-axis is always converted in eccentricity or in inclination.

We note that the phase portraits for the cases where the angle does not include the lunar node have been computed considering the contribution of both Sun and Moon. Moreover, since the harmonic argument depending on the lunar ascending node leads to a non-autonomous resonant Hamiltonian, it is necessary to introduce a dummy momentum $\Gamma$ and a new conjugate angle depending on the lunar node in order to eliminate the explicit linear time dependency. For this reason, there are two first integrals in correspondence of $2g\pm h_{☾}$, $2g+h-h_{☾}$ in Table 8 and Table 9.

The lunisolar periodic component with argument $2g$ produces a non-negligible contribution to the dynamics of the argument of the pericenter, if compared with the oblateness effect and with the precession due to the lunisolar mean terms. Therefore, the ideal model SRM does not give a reliable estimate of the resonant region, in accordance with [5]. Figure 1 depicts the dynamics in the $(e,2g)$ plane and in the $(i,2g)$ plane around the main resonance. The maximum excursion in inclination, as computed with SRM, is $[58.03,67.{36}^{\circ }]$ and it is rather similar to the one obtained with NSA in Table 9, the difference being around one degree both for the minimum and the maximum inclination. The excursions in eccentricity given by the two models are instead quite different: the minimum value of the eccentricity reached in the libration region of the pendulum-like approximation is approximately $0.59$ and is quite different from $e_{min}=0.55$. The maximum values are both above the threshold of $e=0.76$2.

The second dominant term, namely $2\dot{g}+\dot{h}$, is not resonant by definition at the critical inclination. By using the Molniya parameters as initial conditions, the equilibrium point lies in the retrograde orbit region, at $i=110.{99}^{\circ }$ (see Figure 2).

This means that for a Molniya satellite with $(a_{moln},e_{moln},i_{moln})$, the argument $2g+h$ always circulates with a period of approximately $7.6$ years (see Table 5). In any case, the libration region of $2\dot{g}+\dot{h}$ is quite narrow and does not overlap with the other dynamics taken into account, especially with the main resonance and with $2\dot{g}-\dot{h}$, as already pointed out in [5].

In general, also for the other dominant terms identified (see Figures and Figure 4), the dynamics corresponding to $(a_{moln},e_{moln},i_{moln})$ are of circulation, except for the case $2\dot{g}+{\dot{h}}_{☾}$ (Figure ) if we consider a certain tolerance in the inclination, which is reasonable for practical applications [2].

Finally, this analysis shows that initial conditions corresponding to the Molniya parameters allow for a motion in which the width associated with the arguments $2g$, $2g-h$, $2g+h-h_{☾}$ and $2g+h_{☾}$ intersects. By putting together the maximum and the minimum i and e (see Table 8 and Table 9) that may be attained in the libration region of every single resonance, we obtain a maximum overlapping region defined as

$$e\in [0.44,0.79],i\in [59.{30}^{\circ },72.8^{\circ }]$$

(28)

that is widely extended both in eccentricity and in inclination. This is depicted in Figure 5, where we display also the resonance $\dot{h}-h_{☾}$, because it takes place in the same domain. In turn, it is not possible to encapsulate the complexity of the Molniya dynamics in a simple way—for instance, with a single-resonance model or a double-resonance model.

Only three third-order resonances show a ratio larger than a few second-order terms. From Figure 6, which depicts the location of such resonances, we expect that $3g+g_{☾}+h+h_{☾}$, $g+g_{☾}+h+h_{☾}$ and $g-g_{☾}-h-h_{☾}$ occur in the maximum overlapping region found in Equation (28).

5. Conclusions

In this work, we have considered a doubly averaged Hamiltonian model accounting for the Earth’s oblateness effect and the octupolar lunisolar perturbation to estimate the perturbing contribution caused by each term in the specific case of the Molniya satellites’ orbit.

Concerning the evolution of the mean eccentricity and inclination, the role played by the ratio between the harmonic amplitude and the harmonic frequency is crucial to set a “hierarchy” on the large amount of terms that we are dealing with. The analysis provided is general and can be exploited to approach the long-term evolution in other orbital regions. In the case considered here, in addition to the harmonics corresponding to $2g$ and $2g\pm h$, considered in the past to be the most important, our investigation shows that the long-term behavior in eccentricity is influenced, at the same level, also by perturbing terms associated with the satellite and the lunar ascending nodes. The amplitude of the second-order dominant terms (beyond $2\dot{g}$) is at least 1.89 × 10${☾}^{-6}$ for the Sun and at least 1.16 × 10${☾}^{-6}$ for the Moon. The periods of the most relevant terms are instead: 7.21, 7.55, 7.92, 11.79, 12.71, 16.71 and 21 years. In particular, the terms that exhibit a large ratio are also those identified as the dominant contributors through the numerical experiments given in [2], at least as long as the hypothesis of a constant mean semi-major axis holds. Moreover, we have shown that the third-body effect is crucial also for the evolution of the argument of the pericenter, because of the critical inclination.

Besides providing a robust basis for the choice of the terms to consider in an analytical modeling, a possible application of the analysis is to explain observational data or to estimate the maximum eccentricity growth in the long term due to lunisolar perturbations, as done in [19,20] in the case of solar radiation pressure.

Moreover, in the second part of the paper, we have described the phase space associated with the most important resonant terms causing explicitly an eccentricity variation. We have highlighted in particular when the pendulum-like approach does not give a reliable estimate of the resonance extension. In these cases, the critical inclination characterizing Molniya orbits deforms the lobe of the libration region.

The identification of a maximum overlapping region defined by the main resonant contributions, following the phase space investigation, could be a starting point for the further investigation of a possible chaotic behavior. Finally, the third-order lunisolar perturbation does not seem to be particularly significant as regards the dynamics, but it could play a more important role in relation to chaotic phenomena.