## 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 Molinya 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  
Figure 5: Phase space associated with the \(2\delta-h_{\xi}\) dynamics. Contour plot of the Hamiltonian obtained after performing the suitable coordinate change. The pictures show how the phase structure changes by using different initial conditions to evaluate the first integral; \(e=e_{mnl}\) is the initial eccentricity, while the initial inclination varies: \(i=64.3^{\circ}\)**(top left)**, \(i=64.7^{\circ}\)**(top right)**, \(i=64.9^{\circ}\)**(bottom left)**, \(i=66\)**(bottom right)**. Equilibria emerge by increasing the orbital inclination. _Purple_ lines identify the circulation orbits. The separatrices arising from the different unstable equilibria are drawn in _red_ and _black_, while the corresponding librating curves are denoted in _magenta_ and _orange_, respectively. _Cyan_ lines appear only on the bottom left, where there is no clear distinction between the libration region associated with different pairs of recently emerged equilibria. Figure 6: **Left**: in the \((e,i)\) plane, the second-order resonances that correspond to the dominant terms of Table 7 and that can occur in the maximum overlapping region identified (black box). Purple: \(\dot{h}-h_{\xi}\) ; green: \(2\delta\); cyan: \(2\delta-\dot{h}\); orange: \(2\delta+\dot{h}-h_{\xi}\) ; yellow: \(2\delta+\dot{h}_{\xi}\) ; blue: \(2\delta-\dot{h}_{\xi}\). **Right**: the third-order resonances with largest ratio (see Table 7) that can occur in the same neighborhood. In both plots, the red circle denotes the Molinya parameters. Purple: \(3\delta+\delta_{\xi}\) + \(\dot{h}+h_{\xi}\) ; green: \(\delta+\delta_{\xi}\) + \(\dot{h}+h_{\xi}\) ; cyan: \(\delta-\delta\xi\)  – \(h-\dot{h}_{\xi}\). terms associated with the satellite and the lunar ascending nodes. The amplitude of the second-order dominant terms (beyond \(2\beta\)) is at least \(1.89\times 10^{-6}\) for the Sun and at least \(1.16\times 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.