## 1 Introduction  
On 23 April 1965, the first _Mohniya-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 Molniya s 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 _Mohniya 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\geq 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.