Dynamical properties of the Molniya satellite constellation: long-term evolution of the semi-major axis

We describe the phase space structures related to the semi-major axis of Molniya-like satellites subject to tesseral and lunisolar resonances. In particular, the questions answered in this contribution are: (1) we study the indirect interplay of the critical inclination resonance on the semi-geosynchronous resonance using a hierarchy of more realistic dynamical systems, thus discussing the dynamics beyond the integrable approximation. By introducing ad hoc tractable models averaged over fast angles, (2) we numerically demarcate the hyperbolic structures organising the long-term dynamics via fast Lyapunov indicators cartography. Based on the publicly available two-line elements space orbital data, (3) we identify two satellites, namely Molniya 1-69 and Molniya 1-87, displaying fingerprints consistent with the dynamics associated to the hyperbolic set. Finally, (4) the computations of their associated dynamical maps highlight that the spacecraft are trapped within the hyperbolic tangle. This research therefore reports evidence of actual artificial satellites in the near-Earth environment whose dynamics are ruled by hyperbolic manifolds and resonant mechanisms. The tools, formalism and methodologies we present are exportable to other region of space subject to similar commensurabilities as the geosynchronous region.


Introduction
The present manuscript is part of a recent series of papers [1,6] dedicated to astrodynamical properties of Molniya spacecraft. It is well known that the Molniya orbit provides a valuable dynamical alternative to the geosynchronous orbit, suitable for communication satellites to deliver a service in high-latitude countries, as it is actually the case for Russia. In the present contribution, we focus on the decadal evolution of the semi-major axis. We approach the problem by studying the long-term and drag-free motion of a testparticle subject to the non-spherical geometry of the Earth and third-body perturbations due to the Sun and the Moon. The metrical Keplerian characteristic values of the Molniya-class, semi-major axis a (expressed in units of Earth radius, r E ), eccentricity e and inclination i, considered in this work are oe M = (a M , e M , i M ) ∼ (4.16 r E , 0.7, 63.4 • ). (1) To follow this goal, the zonal geopotential terms are first restricted to the second degree J 2 term. Molniya satellites have a mean motion close to 2 revolutions per day and thus are subject to a 2 : 1 resonant commensurability with the Earth's rotation rate (semi-synchronous orbits). Therefore, 12-hour resonant terms of the geopotential need to be taken into account to model the dynamics. The resonant terms are algebraically computed up to the 4th degree and order. Being interested in long-term dynamics, we deal with the various Hamiltonian contributions averaged over the fast variables, leading to the so-called secular dynamics. The fast timescales are connected to the mean anomaly of the test particle and the Moon and the Sun, denoted, respectively, by M, M M , M S . The averaged contributions are introduced as the direct computation of the integral with respect to the fast variables. For the zonal contribution, this averaging is performed in closed form with respect of the eccentricity. The quadrupolar lunisolar perturbations, depending, respectively, on M and M M or M and M S , are doubly averaged, also in closed form with respect of the eccentricity. For the resonant contribution of the geopotential, the averaging requires some extra care. First, the averaging is not performed in closed form over the eccentricity. Instead, we employ a truncated series expansion, which, considering the highly eccentric nature of the orbit, is given to 4th order in the eccentricity. Second, the averaging is not performed over the variable M directly, as it would not take into account accurately the resonant dynamics. Instead, this step calls for the introduction of new slow/fast variables taking into account the very resonant nature of the problem [2]. Once those variables are recognised and introduced explicitly, the averaged contribution is obtained in the usual way, i.e. by averaging over the (new) fast variables.
Molniya, Raduga, Gorizont and Ekran are Russian communication satellites inherited from the Soviet era. Molniya is the Russian word for lightning and, thus, given their interesting dynamical framework, is aptly named. The framework for the Russian communication satellite constellation was first presented by Bill Hilton in the British Interplanetary Society during the years 1959-60 [19] who suggested utilising highly inclined, highly eccentricity orbits for communication systems for high latitude regions. The Molniya constellation operates on a unique orbital configuration which is not exploited by any other type of satellite system. While typical communication satellites operate on a 24-h geosynchronous equatorial orbit, the high latitude of many Russian areas poses a problem for radio frequency transmissions from equatorially orbiting spacecraft. For example, the Russian republic spans a range of 40 degrees in latitude from North to South, with the northernmost point being located at 80 degrees North. The solution to the satellite communication problem for high-latitude regions is the Molniya orbital regime. Operating on a highly inclined and eccentric 12-hour orbit, the Molniya spacecraft reach geosynchronous altitude at apogee, providing access to Russian areas for over 8 hours per orbit, and reach approximately 600 km at perigee at much greater orbital velocities. Since their inception, over 160 Molniya spacecraft have been launched which have provided a platform for research on their unique dynamical framework for nearly 60 years. Molniya orbits gather two distinct resonant phenomena 1 , with quite distinct timescales, giving rise to interesting qualitative dynamical behaviour. Firstly, as we mentioned, they are affected by a 2 : 1 geopotential resonance. Secondly, their inclination close to the critical inclination value of 63.4 • places them near a socalled inclination dependent only lunisolar resonance [22]. While the first affects the semi-major axis of the orbit on a yearly timescale, the lunisolar effect manifests primarily on the eccentricity of the orbit, which exhibits large oscillations on a much longer timescale. These pulsations contribute to modulate the (no-longer constant) coefficients of the tesseral problem; henceforth, a coupling and indirect interplay between the two resonances might happen. The seminal contributions 1 The force model we employed is discussed in more details in "Appendix A". The resonant argument ω which appears in the expansion of the lunisolar Hamiltonian also appears in higher geopotential zonal terms. In this sense, Molniya orbits gather more than 2 resonant phenomena, being affected by zonal, tesseral and third-body resonances. Nevertheless, the effects of higher zonal terms on the semi-major dynamics are negligible for our study and timescale of interests as we will demonstrate later.
regarding the tesseral and lunisolar problems are gathered in [12,13], and in the PhD work of T. Ely [14], later extended to full papers [15,16]. F. Delhaise, J. Henrad and A. Morbidelli [12,13] have focused their study on the eccentricity, inclination and argument of perigee, without paying attention to the behaviour of the semimajor axis. T. Ely [14] connected the resonant problem with large-scale chaos affecting the semi-major axis, including the disturbing effects of the lunisolar perturbation, but for orbital parameters which differ quite significantly from Molniya orbit (in fact, he considered either moderate inclined orbits with i ∼ 20 • or inclinations in the vicinity of the 2g + h lunisolar resonance, i.e. at i ∼ 56 • ). Thus, the secular dynamics of Molniya semi-major axis remains partially unexplored.
The first contribution of this paper is to discuss the dynamics of the semi-major axis beyond the integrable picture. For this task, we rely on classical tools from nonlinear dynamics to portray the dynamical structures organising the long-term dynamics (Poincaré section, sections of finite-time variational indicators). The chaotic nature of eccentric and inclined orbits subject to tesseral resonances, often explained through an overlap of nearby resonances [9], has been known for some time in the context of tesseral resonances [8,10,14]. Nevertheless, as we will highlight, the extent of chaos affecting the semi-major axis phase space for Molniya satellites is much more limited in the range of i ∼ 63 • compared to the previously studied range of inclinations. In fact, large connected chaotic seas are absent from the dynamics. Yet, hyperbolic orbits still exist and surround the unperturbed separatrix as we will show.
The second contribution of this paper is to reveal the precise effects of this coupling on the dynamics of the semi-major axis. This is achieved via the introduction of several dynamical systems, aiming at isolating gradually the various effects and couplings. The driving principle is to introduce basic dynamical models, with the lowest number of degree of freedom (DoF) possible, which still encapsulate the physics and long-term qualitative features of the dynamics. Molniya orbits have also received attentions in [36,37], but predominantly oriented towards the description of the long-term evolution of the eccentricity. The authors have built simplified secular dynamical models, in the same spirit as "isolating" the building blocks of the dynamics and reconstructed the qualitative features of the eccentricity, inclination and argument of perigee observables. A few model generated orbits have been compared to the publicly available two-line element (TLE) datasets (see, for example, [7,33]) 2 . We underline that our contribution is paying particular attention to the orbit of Molniya 1-69 3 , left untouched in a previous study, as being "in the vicinity of the separatrix" [37]. TLEs remain mainly the sole reservoir of orbital data. The TLEs result from observational measures, coupled with an orbit determination process and numerical propagations performed with simplified theories of motion. In this respect, they form rather pseudo-observations instead of "pure" observational data. Approximately every 8 hours, the unclassified TLEs are released publicly. Molniya spacecraft have been tracked since the mid-70's, thus providing a sufficient long-time interval of TLEs to appreciate secular effects acting on the semi-major axis.
The third and final contribution of this paper is the clear connection of the dynamics of two satellites, Molniya 1-69 and Molniya 1-87 4 , with the fingerprints of the dynamics associated to the hyperbolic set. This last point sheds some light of the relevance of secular dynamical approaches and toolboxes for the field of space situational awareness and the continuing increasing space traffic. The patterns of the orbital semi-major axis time series (extracted from the corpus of TLEs, more details will be presented in the subsequent) of the two aforementioned satellites are convincingly approached under this umbrella.
The paper is organised as follows: -In Sect. 2, based on the Earth-only disturbing potential, a secular model is termed. A resonant integrable system is formulated from which analytical quantitative estimates (width of the resonance, characteristic timescales) are extracted. This integrable picture is altered by a multiplet of resonances producing a separatrix splitting phenomena, responsible for the apparition of a chaotic layer in the phase space. For Molniya parameters, the overlap of resonances is complete. The corresponding 2-DoF Hamiltonian and its phase space are described via Poincaré sections. -In Sect. 3, we introduce two models including lunisolar perturbations to overcome the limitations of the Earth-potential only based model. From these models, the effects of the lunisolar perturbations on the tesseral problem are studied. We use dynamical indicators to portray the phase space structures and reveal the hyperbolic set affecting the semi-major axis. The dynamics of the hyperbolic set is studied. -In Sect. 4, after providing more information about the TLEs datasets, we connect the dynamics of the data for satellites Molniya 1-69 and Molniya 1-87 with the dynamics of the hyperbolic set. Relying on our understanding of the underlying dynamics, we extract specific epochs and orbital parameters of the TLEs that spot the satellites within the hyperbolic tangle when computing their respective dynamical maps.
We close the paper by summarising our conclusions.

The secular and geopotential based Hamiltonian
We present our steps and assumptions to recover a relevant secular Hamiltonian model for 12-hour orbits based on the geopotential only and we describe the associated dynamics. The disturbing potential of the Earth, in an Earthcentred and Earth-fixed frame, admits the following expansion [23] with the zonal and tesseral parts, respectively, given by where the vector (r, φ, λ) denotes the spherical coordinates (respectively, radius, latitude and longitude), r E denotes the mean Earth's radius, and μ is the gravitational parameter of the Earth. The P l,m are the Legendre polynomials of degree l and order m. The coefficients c l,m and s l,m are the harmonic coefficients describing Earth's gravity field where we denoted classically J l,0 = −c l,0 . Throughout this paper, we denote the Keplerian orbital elements in the usual way as (a, e, i, ω, Ω, M) where a denotes the semi-major axis, e the eccentricity, i the inclination, ω the argument of perigee, Ω the longitude of the ascending node and M the mean anomaly.

The secular zonal part
The zonal part is dominated by its quadrupole (l = 2) term and we therefore truncate V Z to l = 2. Being interested by secular properties, the M-average of the J 2 part, defining the secular J 2 contribution, is computed (in closed form over the eccentricity) using the differential relationship together with the formula r = a(1 − e 2 )/(1 + e cos f ) (see, for example, [23]): The classical final expression (5), expressed in terms of the orbital elements, reads Secular expressions of higher-order terms and their dynamical effects are discussed in Appendix 6. We note, for the remainder of the manuscript, we drop bars over averaged quantities, bearing in mind that we are dealing with secular functions in this study.

The secular resonant tesseral part for 12-hour orbits
To compute the resonant secular contribution of the tesseral part with we first express it in terms of the orbital elements using a series of formal substitutions. The spherical coordinates are related to the orbital elements by: where α stands for the right ascension of the satellite (again, we refer to [23] for omitted details). The longitude λ is written as a function of α and the sidereal time θ as The sidereal time θ evolves linearly with time as θ = E t, with E = 2π/sidereal day. Writing the inverse of the radius as the quantities sin f and cos f are then written using their infinite series representation as a function of the mean anomaly M and the Bessel functions J s (see, for example, [28]) Applying the aforementioned substitutions into Eq. (7) transforms it into an expression dependent solely on the orbital elements (a, e, i, Ω, ω, M) and the sidereal time θ . The angles appear as linear combinations over the rationales of M, θ − Ω and ω [23]. Computing at this stage the brute-force M-average to derive the secular tesseral contribution would suppress the dynamical effects of the resonant terms for 12-hour orbits. In fact, in the vicinity of 12-hour orbits, the fast angle combines with the fast variable M as to form a slow varying quantity. Therefore, in the neighbourhood of 12-hour orbits, the variable u S needs to be considered as a slow and independent variable. Dealing therefore with the variables u F , u S , ω, there is one fast angle u F and two slow angles, ω and u S . The resonant tesseral contribution is therefore obtained by averaging over the fast angle u F as The final expression has the form where and λ lm is a constant phase term defined as c l,m = −J lm cos mλ lm , In Table 1, we provide the final formal expression of the secular resonant terms 5 for 12-hour orbits up to l = 4 appearing inV T for the uplets k ∈ K with Eq. (12) truncated to k max = 4. Now that we have at hand the secular disturbing functions, the dynamics is cast into a Hamiltonian framework accounting for the Keplerian central part, where T 2 is obtained from (7) by restricting the expansion to l = 2. The Hamiltonian must be a function of canonical variables which are presented hereafter.
We mention however that we might sometimes refer to quantities expressed in orbital elements (non-canonical elements) and it is understood that the elements are themselves function of canonical variables.

Dynamics
We start by introducing the following canonical resonant coordinates where (L , G, H, , g, h) denote the classical canonical Delaunay variables related to the Keplerian elements by Given that the Hamiltonian (19) is time-dependent, θ = E = 2π/sidereal day, we supplement the dynamics with 1-DoF given by the canonical conjugate variables denoted by (Γ, τ = θ), withθ =τ = E . Those variables enter into the definition of (I 4 , u 4 ). The autonomous dynamics (we still note H the new Hamiltonian) reads The Hamiltonian (22) written in terms of the resonant coordinates (20) reduces to a 2-DoF system as both u 3 and u 4 are ignorable. Consequently, their conjugate canonical actions, I 3 and I 4 , are constant over time (i.e. parameters). Note that whenu 2 = 0, the problem is a 1-DoF problem and is therefore trivially integrable.

The integrable approximation
Whenu 2 = 0, we derive the resonant integrable approximation assuming that the resonances are isolated [27]. It amounts to take into account in (19), besides the action-only dependent part, the harmonic with the largest amplitude. For Molniya's orbital parameters, the numerical evaluations of the coefficients give The 1-DoF approximation is therefore built on the oneharmonic Hamiltoniañ The resonanceu 1 = 0 (let us recall u 1 = 2θ − − 2h) that we denote R u 1 occurs for Solving this equation in I 1 for I 2 , I 3 determined by oe M , we find a resonant action I 1 leading to the resonant semi-major axis The orbits of (23) coincide with the set of level curves. Yet, as it will be clear in the subsequent sections, the phase space is analogue to the classical pendulum dynamics. Analytical characteristics of the resonance might be derived from a low-order Taylor expansion of the Hamiltonian near I 1 = I 1 . By keeping only the Table 1 Formal coefficients and resonant angles of the 2 : 1 resonance up to l max = m max = 4 and up to the 4th order in eccentricity quadratic action term, it reduces the Hamiltonian tõ where J 1 = I 1 − I 1 and The equilibria are given by the solution of leading to two the equilibrium solutions The eigensystem of the Jacobian matrix associated to (28) evaluated at the equilibrium solutions (29) shows that x s is elliptic (stable fixed point) and x u is a saddle (unstable fixed point), from which emanates the separatrix (curve associated to the energy level of the unstable equilibria). It also provides further temporal characteristic timescales. The eigensystem of the Jacobian evaluated at x s provides the two complex conjugate eigenvalues (λ s ,λ s ) leading to the characteristic periods of libration in the harmonic regime The eigensystem of the Jacobian evaluated at x u is composed by two real eigenvalues (λ u , −λ u ) defining the e-folding time The resonance half-width ΔJ 1 associated to (26), i.e. the distance between J 1 = 0 and the apex of the separatrix satisfies Solving the last equality for ΔJ 1 , we find where the reader is referred to Fig. 1 for further qualitative details.

The 2-DoF picture
Whenu 2 = 0, the energy function (19) defines a 2-DoF problem with the multiplet of three resonances R u 1 , R u 1 +2g , R u 1 −2g . Each isolated resonant problem admits its own pendulum reduction, with the possibility to overlap [9]. Analytical insights might be gained by some simplifications. In fact, let us approximate (19) with the following 2-DoF problem (35) with φ = u 1 + 2λ 22 , where we have assumed the rate of variation of u 2 = g to be ruled by the H J 2 part, that isu Therefore, u 2 evolves linearly with time which we denote as τ . Using the canonical equations, we find the three resonance centres of R u 1 , R u 1 +2g , R u 1 −2g to be located, respectively, at The mutual distances of the centre of the resonances with respect to the centre of R u 1 , . The corresponding δa amounts to be less than 1 km. Treated as isolated, the As inferred from the numerical computation of h 2,−2 , the resonance R u 1 −2g is negligible for practical purposes. Due to the inequalities a complete resonance overlap takes place (i.e. the resonances are strongly overlapped), by which is meant that the widths of the resonances (treated as isolated) are much larger than their mutual separations. This paradigm is encapsulated into an analogue of the socalled modulated pendulum approximation (see, for example, [27]). From this analogy, we might infer the absence of large chaotic seas known to exist for similar eccentricity range but at lower inclination [10,14]. Instead, we expect chaotic motions to appear only in the vicinity of the unperturbed separatrix [27,29], with a librational region filled by stable orbits. This fact is indeed corroborated by computing the Poincaré map. Stroboscopic map. The Hamiltonian (35) is a 1-DoF system periodically perturbed. Its phase space can be described by computing the associated Poincaré map, which is, given the periodic nature of the forcing, a stroboscopic mapping [26,34]. Let us denote this mapping by P and by V(0) a neighbourhood of J 1 = 0. By defining the lift and projector operators, respectively, as and the stroboscopic map is defined as where Φ t is the flow at time t associated to (35) and T g = 2π/ g . Note that the lift is parameterised by the choice of τ (0) = g 0 . The "dummy" variable Γ does not enter into the equations of motion. For Molniyalike spacecraft, T g defines a period of about 100 years (i.e. the order of 10 4 orbital revolutions). The mapping P is constructed numerically based on the numerical propagation of the system (35). Given a value of g 0 , the coordinates of the fixed points of the mapping P (i.e. the periodic orbits of (35)) are determined using a Newton method. Due to the periodicity the domain of g can be restricted to [0, π]. Let us recall that a fixed point z of P, P(z ) = z , is hyperbolic when the linearisation has at least one eigenvalue with modulus greater than one. In case of complex eigenvalues, the fixed point is elliptic. For g 0 = 0, the two fixed points (semi-major axis given in km) read as Changing g 0 alters slightly those coordinates and the slope of the eigenvectors associated to the unstable periodic orbit, which may widen the aperture of the librational domain by a few kilometres. The stable and unstable manifolds associated to an hyperbolic point z , are grown by iterating points belonging to the fundamental domain I ⊂ E s,u , where E s,u are, respectively, the stable and unstable eigenspaces associated to z (and derived from the eigensystem analysis). Recall that W s,u are locally tangent to E s,u . In Fig.2, we show the Poincaré section containing a chaotic zone surrounding the "unperturbed separatrix". A smaller portion of the phase space shows the first lobes associated to the stable manifold. The analysis of the eigensystem associated to the linearisation of P at the saddle fixed point is enlightening in deriving the Lyapunov timescale analytically. Let us recall that a Floquet characteristic exponent μ is a complex number satisfying where λ is an eigenvalue associated to the linearisation DP about the fixed point. For the hyperbolic saddle, the two eigenvalues {λ 1 , λ 2 = 1/λ 1 } are real and so are the corresponding {μ 1 , μ 2 }, called in this case the Lyapunov exponents. From the timescale of 1/μ ∼ 17 years is derived for the largest eigenvalue. This timescale has been compared with a brute-force estimation of the maximal Lyapunov exponents χ based on the variational dynamics (and their associated Lyapunov time τ L = 1/χ ) in the vicinity of the hyperbolic saddle. For hyperbolic orbits, we found Lyapunov times in the range of 15−18 years, thus in very good agreement with the analytical timescale based on DP. Consequences of the chaotic layer. The presence of the thin chaotic layer surrounding the unperturbed separatrix brings important distinguishable qualitative features to the dynamics: the semi-major axis might display intermittency phenomena and the resonant angle alternates between librational and circulational regimes. We note that such features have been observed for simulated geosynchronous orbits [5,35]. More precisely, for initial conditions in the chaotic layer, the and the "outer-circulation" regime for which The alternation takes place when the orbit returns close enough to the hyperbolic saddle x u where the scattering takes place.
Remark 1 This mechanism is illustrated and summarised within the composite panel in Fig. 3 based on the Hamiltonian model For = 0, T is integrable and has a saddle structure at (Λ, ρ) = (0, π). The separatrix has a cat-eye topology with half-width Δ = 2 √ 2. When = 0, 1, the resonances R ρ and R ρ+ ρ 1 , -apart, produce a separatrix splitting. The resonant angle of an orbit with initial condition in the hyperbolic set alternates among libration, Λ ρ 0, and circulation, Λ ρ 0. The "projection" of one orbit with initial condition close to the saddle (trapped in the hyperbolic tangle) in the space (Λ, ρ) shows that the orbit remains mainly guided by the unperturbed separatrix. Under our selected initial condition, when the angle circulates, the action is trapped in the tangle, evolving here in the domain Λ − := {Λ, Λ < 0}. When the angle librates, the action experiences full homoclinic loops and evolve within Λ = Λ − ∪ Λ + This process continues and possibly alternates in the vicinity of the saddle, producing scattering and contributing to the growth of the tangent vector. (35) is based on geopotential perturbations only. To build a more realistic model, the lunisolar perturbations, Moon and Sun, need to be included. In its present form, model K is limited in two ways: 1. Under the lunisolar effects and due to the proximity to the critical inclination value, the hypothesis that the argument of the perigee (g = u 2 ) flows linearly with time at a (constant) rate given by the J 2 effect is violated (a fact also observed at the data level). 2. The assumption that both the eccentricity and inclination are parameters is no longer true under the influence of the lunisolar perturbation. Increasing the complexity of the model gradually, we overcome the first limitation by decoupling the equations of motions. We isolate a simplified energy function L that dictates the time evolution of the argument of perigee,ġ = ∂ G L, that we use to form a 6-dimensional dynamical system with constant eccentricity and inclination. The variables (J 1 , u 1 ) are then studied. The second limitation is raised by introducing a 3-DoF Hamiltonian system, where both the eccentricity and inclination vary according to the dynamics.

The doubly averaged lunisolar Hamiltonian
We adopt a simplified sub-model of the quadrupolar doubly averaged formulation to model the external third-bodies perturbations. The quadrupolar approximation is commonly employed to study medium-Earth orbit dynamics and has already demonstrated its relevance (see, for example, [11,17] where r M , r S denote the geocentric vectors of the Moon and the Sun, r M , r S the corresponding geocentric distances and μ M , μ S their respective gravitational parameters, a Legendre-like expansion of H M and H S , truncated to l = 2 (quadrupolar hypothesis), and averaged over the mean anomalies (M, M M ) and (M, M S ), respectively, defines the so-called doubly averaged third-body model. This averaging is performed in closed form over the eccentricity. Contrarily to the inner-perturbative part (geopotential), it requires to use the differential relationship coming from Kepler's equation (E refers to the eccentric anomaly). The double-averaging reduces H M and H S to an expansion of the form where φ M i and φ S j are permitted linear combinations of (ω, Ω, Ω M ) and (ω, Ω), respectively. Given that the angle Ω M does not enter φ S j , the summations are homogenised by introducing φ q , where with the convention that q 3 = 0 for the permissible solar arguments. The quadrupolar doubly averaged lunisolar Hamiltonian reads therefore with (under the quadrupolar assumption) For the sake of concision, let us denote h MS j = h M j +h S j . The Hamiltonian (50) is in general non-autonomous as time enters through the ecliptic precession of the lunar node, well-approximated by the linear law (Moon's elements are referred to the ecliptic plane) where 2π/| Ω M | defines a period of about 18.6 years. The Moon's inclination to the ecliptic plane is set to i M = 5 • 15. However, as we will see hereafter, the simplified model allowed by Molniya's parameters leads to a model independent of the argument of the Moon and, therefore, to an autonomous model.

Model for the time evolution of ω
The so-called double resonance model employed in [36] based on the lunisolar harmonics cos 2g and cos 2g±h has shown to provide a realistic model to capture the time evolution of the argument of perigee. They compared orbits generated using this model against TLEs data on several cases and were able to reproduce qualitatively the time evolution of the argument of perigee on several decades. More recently, [31] advocated that the cos h term produces a significant contribution to the dynamics of ω, where the reader is referred to Appendix 6 for further details. We adopt the following 2-DoF Hamiltonian system to model the time evolution of ω. From the canonical equations derived from (53), we derive the dynamics of g. The formal coefficients appearing in (53), expressed using the Keplerian elements, are listed in Table 2. The terms h M 0 and h S 0 refer to the action dependent only terms of H M and H S and read

Effects of lunisolar perturbation on the tesseral dynamics
We investigate now how the lunisolar perturbation affects the hyperbolic structures of the tesseral problem for Molniya parameters. The basic model dictates the time evolution of g being established by (53); we focus now on two dynamical systems improving the caveats of (35): 1. First, we introduce the differential system in R 6 defined by the equations of motion (EoM): u 1 , g), u 1 , g), H, g, h), H, g, h), H, g, h), H, g, h).
where J is defined as, g(t)).
Here L is the basic lunisolar Hamiltonian function. For short, we refer to the EoM (54) as model J and we denote the right-hand side by v J . G, H, u 1 , g) + L(G, H, g, h),

Second, we consider the 3-DoF Hamiltonian
expressed within the resonant coordinates. With respect to the model J , model S allows the action terms to evolve under the correct dynamics.
In particular, the tesseral coefficients h 2,0 (a, e, i), h 2,±2 (a, e, i) appearing in the dynamics of (I 1 , u 1 ) are no longer frozen (instead, they vary according to the changes of the Delaunay action vector (L , G, H )). The right-hand side derived from S is denoted by v S .
Let us emphasise that both models are π -periodic in g.
For both models, in order to reveal the dynamical template on the (I 1 , u 1 )-plane, we compute the Fast Lyapunov Indicators [18,25] on a 500×500 Cartesian mesh of initial conditions. We use the following definition of the FLI. For an n-dimensional autonomous ordinary differential system defined on a open domain D ⊂ R n , x = f (x), the FLI at time t is derived from the linear map D x f at a point x: and the associated variational equations as The vector w ∈ R n denotes the tangent (or deviation) vector. The computation of the FLIs over resolved grid of initial conditions discriminates efficiently the structures of a given dynamical system, including the stable or unstable manifolds, ordered or chaotic seas. One advantage of the FLI over the characteristic Lyapunov exponent is to get rid of the time-average computation, thus speeding the stability determination process. For regular orbits, the deviation vector grows linearly with time and therefore the FLI on regular KAM tori is characterised by values close to log(τ f ). In hyperbolic regions, the norm of the tangent vector grows exponentially fast, and therefore the FLI displays a linear trend surpassing quickly the value taken on KAM objects (see [3], chapter 5, for perturbative estimates).

Remark 2
The parametric dependence on t in (59) is raised after a calibration procedure. In our case, integration of several single orbits showed that τ f = 20 years are sufficient to obtain a sharp distinction. With our choice of units, the FLI of regular orbits is characterised by the value log(τ f ) = 4.99.

Remark 3
We restricted the computations of the FLIs over Σ forward in time, i.e. on a time interval [0, τ f ], τ f > 0, to obtain "positive in time FLIs", FLIs + . It is therefore understood that, in the context of the existence of hyperbolic invariants, these computations on Σ would reveal the trace of the stable manifolds. Similarly, backwards in time FLIs computed over [−τ f , 0], FLIs − would reveal the trace of the unstable manifolds. Both manifolds can be displayed on Σ by plotting, for example, the standard average or any others weighted average (see, for example, [18,25]). We computed a few of those maps backwards in time, to display the averaged FLI. However, because we are not particularly interested of highlighting homoclinic connections, we present hereafter only the forward in time FLI maps (i.e. we display FLI + ).
The FLIs computation are performed for the vector fields v J and v S over a Cartesian discretisation of Σ ⊂ R 2 , where and an initial deviation vector w 0 chosen arbitrarily 6 . The neighbourhood of the resonant action V(I 1 ) represents typically a range of 70 km in the semi-major axis. Our numerical campaign is parametric through 4 allowed values of i 0 , namely "piercing" the critical inclination value. This choice enters Σ i 0 through Although aware that the precise geometry of the hyperbolic structures depend on the initial phasing (ω, Ω), our investigations focus on (ω, Ω) = (270 • , 0). All the resulting maps of this numerical survey for models J and S are reported in Appendix B to ease the readability. We show hereafter in composite panels only the relevant information for the analysis. From this survey, we observe that: 1. Both models display a saddle-like point in Σ. This suggests the existence of an unstable periodic orbit (although this invariant has not been computed) for both flows, similar to the unstable periodic orbit we computed for model K. Following this idea, the hyperbolic set (high values of FLIs with yellow colour) emerging from the saddle-type structure is very likely to represent the intersections of the stable manifold of the hyperbolic invariant with Σ. The fine mesh of initial conditions allows to recognise lobes distinctively for model S. 2. The model J is overall weakly perturbed, and the hyperbolic layer is very close to the unperturbed separatrix. 3. On the contrary, the hyperbolic layer of S is much more developed. This fact is imputable to the indirect modulation of the coefficients h 2,0 (a, e, i) and The dynamics associated to the hyperbolic layer of model S is similar to model T apart that the coefficients of the respective resonances are slowly modulated in time. This is exemplified for two orbits in the composite panel of Fig. 4, together with macro-and microviews of the phase space structures. The orbit immersed within the stable region displays oscillations, while the orbit trapped into the hyperbolic layer displays the characteristic intermittency. Let us underline that the hyperbolic orbit, on the 20 year timescale, displays U-turns (i.e. the alternation between libration and circulation regimes of the resonant angle u 1 ) always directed towards lower semi-major axis, with a timescale of about 1.5 year. The full homoclinic loop takes about 3 years. We integrated the same orbit on a time interval 10 times larger and we noticed the unevenly distribution between upper and lower U-turns, the latter being more frequent. This property is clearly inferred from the thorough inspection and detailed geometry of the hyperbolic foldings becomes larger along a solution x(t) = a(t), e(t), i(t)) for initial conditions near x u . The layer's dependence upon ω 0 , for model S, is shown for the fixed value of u 1 = 3.6 in the last map of the composite plot of Fig. 4 (bottom right). It reveals a much wider width (roughly speaking on the order of 10 km in the semi-major axis) for ω ∈ [π, 3π/2], with petals structures. For ω ∈ [3π/2, π], the hyperbolic structure is much simpler to apprehend. For Molniya's typical variation of ω ∼ 270 • ± 20 • , this corresponds to the rectangle materialised with white dashed lines.

Connections and links with the dynamics of Molniya 1-69 and Molniya 1-87
On inspecting the extracted semi-major axis using Molniya 1-69 and Molniya 1-87 TLE data, we notice that they display intermittency phenomena on their semimajor axis 7 as repeated in Figs. 5 and 6. The figures also incorporate the time evolution of the resonant angle u 1 . The relevant part of the data, in the light of the oscillating models previously derived, cover in both cases at least 2 decades. Both data contain a transitory period, possibly remnants of unknown manoeuvres. For Mol- Fig. 6 Time history of the semi-major axis and resonant angle u 1 extracted from the TLE data for the satellite Molniya 1-87 niya 1-87, after the epoch corresponding to mean Julian day (MJD) of 57, 400, the satellite experienced a significant semi-major axis reduction. We will not pay attention to this part of the data. In the exploitable window, the resonant signature, consisting of alternation between libration and circulation, is well-marked and in accordance with the U-turns intermittent semi-major axis variations. In both cases, the intermittency U-turns take place for a ∼ 26, 550 km, compatible with the locations of the hyperbolic foldings we located with model S close to the saddle. It is worth mentioning that Molniya 1-69 has been left untouched in [36], as judged to "locate in the vicinity of this separatrix". Below, we give more credit to this claim, and we show that it is also the case for its cousin Molniya 1-87. At the light of the dynamical mechanisms presented in Sect. 3 and the fingerprints just described, it is tempting to say that both satellites evolve within the hyperbolic layer. To give more weight to this claim, we performed the following steps: 1. At epochs t corresponding to the apex of the first U-turns, we extract from the TLEs the corresponding orbital parameters and we record the values of (a , u 1 ).
2. We compute the dynamical maps with the FLIs for the vector field v S , using as parameters and phasing for the section Σ those extracted from the respective TLE at epoch t . 3. On the obtained dynamical maps, the points of coordinates a(t ), u 1 (t ) are spotted.
The obtained dynamics maps shown in Fig. 7 convincingly demonstrate that the satellites reside within the hyperbolic tangle.

Conclusions
The constructed dynamical models and their analysis allowed us to deepen the understanding of Molniya's semi-major axis dynamics. The hyperbolic structures organising the phase space have been portrayed via variational indicators through a series of compact, tractable and realistic secular models. The effect of lunisolar perturbations, on the 20-year timescale, needs to be taken into account to reconstruct the correct dynamical template. In fact, the induced modulations of the eccentricity and inclination contribute sufficiently to change the "parameters" of the tesseral problem; the coefficients we denoted by h 2,0 and h 2,±2 . We connected the 20-year-long fingerprints of two satellites, Molniya 1-69 and Molniya 1-87, with the hyperbolic layer surrounding the unperturbed cat-eye separatrix. This hyperbolic layer, in absence of lunisolar perturbations, would be too thin to sustain the dynamical signatures visible at the publicly available data level. By computing their associated dynamical maps, we provided evidence that the two satellites are trapped within the hyperbolic tangle. The secular dynamics umbrella provided a reliable and robust mould to approach and explain the semi-major axis patterns extracted from the TLE space datasets. As far as we are aware, this result is the first report of long time scale hyperbolicity corroborated by pseudo-observations in the near-Earth space environment. The mechanisms and tools depicted in this contribution have relevance for other dynamical regions, most notably for the geosynchronous altitude where similar patterns have been observed on simulated orbits [5,30,35].
J.D. and A.L. acknowledge several discussions with Alessandra Celletti and Cȃtȃlin Galeş on the resonant potential. J.D. acknowledges discussions with Ioannis Gkolias on the J 2 2 effect and useful references provided. J.D. acknowledges several discussions all along this research with Aaron Rosengren.
Data Availability Statement TLEs data are available at spacetrack.org.

Conflicts of interest
The authors declare that they have no conflict of interest.

Appendix A: Dynamical models employed & physical parameters
As we mentioned in the introduction, the aim of this study is not to study with the greatest accuracy possible Molniya dynamics, with a comprehensive force model including uncertainty modelling and Monte Carlo like approaches. Quite on the contrary, we leverage the understanding of the dynamics of the semi-major axis from the essential "building blocks" with tractable contributions. In that respect, we would like to provide more context to the dynamical model we have employed. We have approached the problem as a dragfree model, with no solar radiation pressure, based on a compact geopotential model including relevant terms of the disturbing lunisolar potentials. Higherorder zonal secular terms can be obtained in closed form over the eccentricity following the same formal procedure discussed in Sect. 2.1. In terms of the orbital elements, up to order l = 5, they read: It is worthwhile to note that the resonant argument of perigee also appears in the above secular contributions; hence the idea that Molniya orbits, besides tesseral and lunisolar resonances, gather also "zonal resonances".
To include the second-order part term with factor J 2 2 in the secular Hamiltonian, with the form V J 2 2 = J 2 2 A(a, e, i) cos 2ω + B(a, e, i) , we used the formula given in [4,24]. The relevance of our model S has been assessed by including those effects, and the lunisolar h MS h cos h to L. This model forms an "extended" Hamiltonian modelS. We computed the dynamical map for the Hamiltonian vector field vS with e = 0.7, i 0 = 64 • 3 and (ω, Ω) = (270 • , 0) and we did not noticed significant macroscopic changes in the obtained dynamical template, henceforth the relevance of the Hamiltonian model S. Let us mention that even if the macroscopic structures do not change drastically, hyperbolic orbits generated under model S andS will separate in time (sensitivity to the slight change of physics), and the hope to follow them beyond a few Lyapunov times is a useless effort.

B Dynamical maps
We computed dynamical maps for a fixed value of e = 0.7 and i 0 "piercing" the critical inclination. They are presented in Fig. 8 for model S. Given that model  Fig.9. The maps have been computed on a 500 × 500 grid of initial conditions, forward in time, and over a time interval of 20 years. We have considered 4 values of the initial inclinations, namely i 0 ∈ {62 • 5, 63 • 4, 64 • 3, 65 • 2}. The initial phasing is set as (ω, Ω) = (270 • , 0 • ). If a given initial condition in the map fall within the highest region of the FLIs (yellow tone), then the orbit is hyperbolic and exhibit sensitive dependence upon the initial con-dition (i.e. any orbit starting with an initial condition slightly different will have a long-term different future; the orbits will separate with time). We note that the i 0dependence of the J model is quasi-absent. The model J is very close to the integrable picture, in the sense that the splitting of the separatrix is weak. The latter is much more manifest for model S, where we recall, the eccentricity and inclinations variables are no longer frozen. For increasing values of i 0 , we underline the growing asymmetry of the foldings near the saddle- like structure for the model S. This particular structure transfers directly at the single orbit level: an orbit trapped within the hyperbolic layer is more likely to display U-turns intermittency phenomena towards the lower semi-major axis. This observation, based on the thin structures of the lobes detected with a variational dynamical indicator on our model, is also in agreement with the actual two-line elements datasets for objects M1-69 and M1-87.