Drift and its mediation in terrestrial orbits

The slow deformation of terrestrial orbits in the medium range, subject to lunisolar resonances, is well approximated by a family of Hamiltonian flow with $2.5$ degree-of-freedom. The action variables of the system may experience chaotic variations and large drift that we may quantify. Using variational chaos indicators, we compute high-resolution portraits of the action space. Such refined meshes allow to reveal the existence of tori and structures filling chaotic regions. Our elaborate computations allow us to isolate precise initial conditions near specific zones of interest and study their asymptotic behaviour in time. Borrowing classical techniques of phase- space visualisation, we highlight how the drift is mediated by the complement of the numerically detected KAM tori.


Introduction
Various groups of scientists have become enchanted anew by the lunisolar resonances affecting the dynamics of terrestrial orbits. The study of them and the resurgence of their significance has not been visible since the notorious and colossal triptych of Breiter (1999Breiter ( , 2001a. Later rebranded by Rossi (2008) in the context of the medium-Earth orbits (MEOs), the study of their long-term dynamics, and in particular their eccentricity growths in the elliptic domain (Bonnard and Caillau, 2009), represent current deep motivations for the community. In our opinion, the most complete and up-to-date panorama of the literature is excellently presented by Armellin and San-Juan (2018). The existence of such a condensation allows us to adopt here a rather direct style in this present contribution. We are particularly interested by questions related to the stability of orbits. Based on the divergence of nearby trajectories, the existence of a mixed phase space where there is a cohabitation of stable and chaotic components has been recently pictured Gkolias et al., 2016;Celletti et al., 2017) and partially explained applying Chirikov's resonances overlap criterion (Chirikov, 1979). The Hamiltonian flow obtained under the simplest assumptions for the disturbing effects of the perturbers (i.e., a development restricted to its lowest order and averaged over fast variables), Moon and Sun, encapsulates all the details of the dynamics in which we are interested . In particular, the Hamiltonian possesses 2 degrees-of-freedom (DOF) and depends periodically on the time t (see Celletti et al. (2017) for omitted details). We recall in the next section how the Hamiltonian (1.1) H : DˆT 2ˆT Ñ R, px, y, tq Þ Ñ Hpx, y, tq " h 0 pxq`εh 1 px, y, tq with h 1 px, y, tq " ř mPAĂZ 3 ‹ h m pxq cospm¨py, tq`φ m q is obtained. The form of Eq. (1.1) is the standard form of a nearly-integrable problem written in action-angles variables. The non-linearity parameter ε belongs to a certain subinterval I of R`and is function of the semi-major axis, which is a first integral in the secular approximation. The functions th m u mPA are real valued functions of the sole action x P D Ă R 2 (and some constant physical parameters), the φ m are phase terms. When ε sweeps I, a transition from a globally ordered phase space to a mixed phase space is known to exist. It turns out that the presence of a chaotic regime for large values of ε, say for ε close to max I, corresponds to the range of semi-major axes where the navigation satellites are located. The occurrence of the two apparent antonyms, 'how awkward it is' 1 and 'how useful and fruitful it can be' 2 , crystallizes assuredly the challenges, implications and beauty of the dynamical and engineering problems we face. Gravitational problems are kaleidoscopes of pure and applied science. Our Solar System has been the source and receptacle of many theoretical and practical dynamical facets and aspects (KAM theory, hyperbolic dynamics, shadowing theory, numerical analysis, phase visualisation techniques). Spaceflight dynamics is not excluded and has gained from this rich heritage (Bollt and Meiss, 1995;Koon et al., 2000;Perozzi and Ferraz-Mello, 2010). We cannot take a definitive position on the space-debris mitigation via 'chaos targeting' and transfers in phase space, nevertheless, let us underline that the concept embraces the continuous necessary exchanges between (technological and scientific) communities.
In this paper, we depart from former goals where the main impetus was the explanation of the mechanisms supporting the apparition of chaos. Instead, we focus rather on i) the physical consequences in terms of transport in the phase space and ii) on the visualisation of these excursions via double sections in the action-like phase space. The techniques we used have been extensively employed in Dynamical Astronomy and overall in the context of the dynamics of quasi-integrable Hamiltonian systems and symplectic discrete maps (confer Lega et al., 2003;Todorović et al., 2008;Cincotta and Giordano, 2008;Páez and Efthymiopoulos, 2015;Lega et al., 2016, just to name a few). To achieve our tasks, we provide a cartographic view of the prograde and retrograde region in Section 3.2, based on a lighting-fast ad-hoc secular model that we recall in Section 3.1. The fine resolutions of the meshes used to discretize the domains D allow for highly detailed views of the phase space. We then focus on the computation of diameters-like quantities to relate the degree of hyperbolicity (a local property) with a more practical transport-like index (a global property). Thanks to our resolved grids, precise initial conditions (ICs) can be extracted, which lie near specific structures of interest, in particular where large diameters are expected. Once obtained, we proceed to their asymptotic analysis (in time) using ensemble orbit propagation (Section 4). We close with Section 5 where we summarise our contributions and discuss an open problem that inspires our future efforts.
1 The Lyapunov times τ L , which dynamically speaking constitute the barriers of predictability, are on the order of decades .
2 There is a birth of a new ideology to remedy the space-debris problem, based on a 'judicious' use of the instabilities to define re-entry orbits and navigate the phase space.

The model
We recall briefly, for the sake of completeness, under which hypotheses the 2.5-DOF Hamiltonian is obtained. After the presentation of the model, we present to the newcomers a few facets of the resonant aspects.
2.1. Derivation of the Hamiltonian model. Numerical evidence has shown that, for the range of the treated perturbation I (recall I » r2.2r C , 4.65r C s), refinements of the gravitational potentials beyond the quadrupolar level are not necessary to capture details of the global dynamics we are interested in, even on long timescales . It means that when the potentials of the Earth and those of the external bodies, Moon and Sun, are developed using Legendre expansions, terms with l ą 2 are disregarded. By recognising the timescales of the dynamics, further simplifications are even possible to get a more pertinent analytical model (and numerical as well; see also Section 3.1). Based on the Lagrangian averaging principle (Grebenikov, 1965;Mitropolsky, 1967;Ghys, 2007), the potentials are averaged over the mean anomaly of the test particle and those of the third bodies 3 , @ and K .
For an oblate Earth, we recall the classical averaged potential Here pG, Hq denotes the second and third variables of the Delaunay actions pL, G, Hq, µ C denotes the (Earth's) gravitational parameter. The canonically conjugated vector of angles is classically denoted p , g, hq. Omitting details that might be found in Celletti et al. ( , 2017, the disturbing function of the Sun's attraction, R @ , reads as The scalar α @ " µ @´a 2 a 3 @¯p 1´e 2 @ q´3 {2 has a constant magnitude of " 3.96ˆ10´1 4 in the international system of units. The coefficients s m are defined as s m " K m p2´mq!{p2`mq!. The functions F 2,m,p p‚q refers to Kaula's inclination function (Kaula, 1966) and H 2,p,2p´2 peq is related to the Hansen coefficients. For the disturbing function of the Moon, the following formula holds true The expressions U m,˘s 2 are function of the obliquity of the eccliptic, and are present due to a rotation of the spherical harmonics needed in this mixed-reference frame formalism. Note that the size of the coefficient α The coefficients m m,s are defined as m m,s " p´1q rm{2s K m K s p2´sq!{p2`mq!. It turns out that the time derivative of the angle h K is well approximated by a constant frequency defining a period of 18.6 3 In the following, we use the subscripts ‚ C , ‚ K , ‚ @ to denote parameters referring to the Earth, Moon and Sun, respectively.
years. In other words, we consider the explicit time dependence of the lunar potential as linear. At this stage, it is recognisable and transparent that the Hamiltonian formed on the perturbations, H " H J2 pG, Hq´R @ pG, H, g, hq´R K pG, H, g, h, tq, (2.7) possesses 2 DOF and is periodically-time dependent (i.e., a 2.5-DOF problem). The explicit time dependence due only to the node of the Moon 4 plays a fundamental role in shaping the dynamics. The well-known distinctive feature with the case of 2 DOF is that, a priori (in absence of additional known first-integrals apart the energy function itself), the tori cannot act as practical barriers preventing transport in the phase space (for an N -DOF autonomous problem with N ě 3, the codimension between the N -dimensional tori and the dimension of the phase space restricted to an energy surface (2N´1) is at least 2). The Delaunay variable being a cyclic variable, its conjugate variable L " ? µa is a constant of motion. Let us introduce normalised new actionsx " x{ ? µa. The reduced system is kept canonical as long as the new anglesỹ " ? µa¨y are introduced and the physical-time multiplied by the same factor. It is clear that the new Hamiltonian has the same form as in Eq. (2.7). The previous factor α J2 absorbs now a contribution from L and we get the new α J2 " J 2 r 2 C µ 4 {4L 6 . Factorizing the external perturbation by the greatest α K , the Hamiltonian can be rewritten as The hierarchy α K ! α J2 enables us to write α K " εα J2 , ε ! 1, and Eq. (2.8) becomes HpG,H,g,h, ? µatq " h 0 pG,Hq`εh 1 pG,H,g,h, tq. (2.9) Clearly H shares the form of the standard perturbed Hamiltonian system as announced in the introduction. The very useful information that we got from these manipulations is that the dimensionless perturbative parameter ε is proportional to the secular invariant semi-major axis, . (2.10) (The mean motions of the test particle and disturbing bodies are noted n and n K , respectively.) Note that this perturbing parameter is of the same nature as that introduced by Breiter (2001b), but we are treating herein the regime of the lunisolar secular (not semi-secular) resonances. The Hamiltonian model based on the quadrupolar level is physically relevant up to a semi-major axis close to a max " 6r C (beyond, octupolar refinements, l " 3, are needed) corresponding to εpa max q " 0.8. In the following, we will be interested in semi-major axes up to a " 29, 600 km, leading to " 0.22. From our numerical investigations, we noted that for a " 13, 600 km, the chaos is thin and confined to a few inclinationdependent-only resonances. These two constraints together define the subinterval I " r0.004, 0.22s Ă R`of interest for ε. Adding quite 'virtually' the point t0u to this set, ε P´t0u Ť I¯, we obtain when ε " 0 an integrable dynamics with a linear flow on a torus. The actions are constant and determine the invariant tori. On these tori the dynamics consist of a rotation at constant speed characterised by the vector of constant frequencies (the unperturbed frequency vector ) ΩpG, Hq " p g , h q given by κp5 cos 2 i´1qp1´e 2 q´2, h "´κ cos ip1´e 2 q´2, Kaula, 1966). Let us be more precise regarding the definition of the interval I and the energy function considered. An important factor leading to the 2.5 DOF model lies in the omission of the tesseral contributions in the Hamiltonian (2.7). Tesseral resonances occur when the commensurability 9 M { 9 θ " q{p takes place.
Given the upper and lower bounds of I, the main tesseral resonances affecting the motion are given by the set of commensurabilities T " t2 : 1, 3 : 1, 4 : 1u (see e.g., Ely, 1996;Celletti and Galeş, 2015). Near a tesseral resonance, the semi-major axis is no longer secularly invariant and, in fact, might experience (confined) chaotic variations near their corresponding resonant action value L p:q ‹ . Putting all of these together, a more precise definition of our interval of interest for the perturbing parameter ε is instead I Ť p:q P T rεpL p:q ‹´δp:q , εpL p:q ‹ q`δ p:q s where δ p:q characterise the strength (width) of the resonance. The eventual local coupling of the tesseral contributions and the lunisolar resonances on the evolutions of the actions pG, Hq, for perturbing values precisely within the tesseral resonant domain (i.e., for ε P rεpL p:q ‹´δp:q , εpL p:q ‹ q`δ p:q s), is not discussed in this contribution, but is the object of a further study (accordingly, tesseral contributions are disregarded in the energy function). Using estimations obtained in former works (Ely, 1996;Celletti and Galeş, 2015), let us stress that the ratio of the measures mpI 1 q{mpIq " 0.99 is very close to 1. Moreover, accurate low-altitude or very local studies could require additional refinements of the energy function (e.g., near the critical inclination, see Lara, 2018). Although our model is subject to these possible limitations, we should highlight that the considered dynamics are accurate enough to describe general MEO orbits, and from a mathematical point of view provide a simple testbed to investigate transport theories and capture the big picture.
We abused the vocabulary and treat the eccentricity e and inclination i as 'actions', instead of using the veritable actions variables pG, Hq. This is rather to stick to classical notations (Kaula, 1966). Nevertheless, these variables are functionally independent and the 'true' actions easily obtained as e 2 " 1´pG{Lq 2 and cos i " H{G. Let us precise that, when dealing with the autonomous Hamiltonian by introducing an extra conjugated variables pΓ, γq P RˆT for ε ‰ 0, one couple of actions x " pG, Hq characterise an invariant torus of T 3 since Γ does not enter into the equations of motion. In other words, we can consider the orbits in the reduced phase space defined by X " pG, H, g, h, γq, x " pG, Hq P D Ă R 2 , y " pg, h, γq P T 3 ( . In Section 3, we will offer views of the dynamics in actionaction sections, meaning that within the space X , we fix the angles to a specific vector to obtain the section S " pG, . Let us now discuss fundamental phenomenon for ε ‰ 0.

Secular Lunisolar Resonances.
A determinant feature in the long-term properties of nearlyintegrable systems of the form hpx, yq " h 0 pxq` h 1 px, yq is the presence of resonances 5 (Lochak and Meunier, 2012). They arrive when a vector k P Z n ‹ satisfy with the (unperturbed) frequency vector a commensurability condition over the rationals. The resonant condition reads k¨Ωpxq " 0. For a fixed vector k P Z n ‹ , the sets (potentially empty) of the actions x such that k¨Ωpxq " 0 form the resonant manifolds. The resonance under consideration is then characterised by an index, the resonance order, usually though the 1 -norm of k, k 1 " ř i |k i |. Under the quadrupolar assumption, the system (2.9) is prone to resonate with a maximal order of 6. Let us consider the frequency vector Ωpxq " p g pxq, h pxq, K q, then, as already recognised by Ely (1996), the resonant conditions read as These algebraic equations admit non-trivial solutions that define the lunisolar resonant manifolds. The resonant manifolds are mirrored with respect to the resonance p0, k, 0q¨Ωpxq " k h pxq " 0. (However, as we will clearly illustrate it, the symmetry of the resonant manifolds does not imply a mirroring of the geography of the KAM tori and hyperbolic structures.) In Daquin et al. (2016), the extent of the resonant zones have been estimated (in a subdomain of the prograde 0 ă i ď π{2 domain) by reducing the Hamiltonian to the first fundamental model of resonance, a pendulum. This procedure involved the introduction of resonant coordinates through canonical transformations T k P SLp3, Qq leading to an intuitive physical interpretation of Chirikov's overlap as a driver of chaos. However, since the work of , it has been observed that such a reduction does not always capture the features of the dynamics. In order to get a more refined and precise view of the extent of chaos, a superior way is instead to look at the destruction of KAM curves, e.g., using fast dynamical indicators.

Phase-space views
We revisit and complement the transition order/chaos in terrestrial orbits by scanning the dynamics under the rays of a first order variational chaos indicator. The motivation is twofold: (1) the resolutions used in former studies are generally sufficient to detect and isolate chaotic components; yet, they are too coarse to detect the eventual presence of structures immersed within them. In addition, the dissection of the dynamics with a fine resolution makes possible the extraction of the chaotic skeleton with surgical precision 6 . This property will be used to study transport properties (confer Section 4).
(2) Gkolias et al. (2016) claimed that 'the retrograde orbits are not intrinsically more stable then their prograde counterparts'. This diagnosis was established by scanning the region with an averaged FLI (over some angles) focusing on low eccentricity (up to e " 0.1). We feel necessary to investigate further this assertion beyond e " 0.1 without the averaged indicator (which naturally tends to smooth and absorb the details). To overcome and constrain these two symptoms, we first briefly recall how we efficiently deal with the equations of motions, after which we present and discuss our highly-resolved phase-space views on a macroscale.
3.1. Numerical treatment of the equations of motions. Ordinary and partial differential equations with disparate scales (spatial, temporal) are numerically challenging. The difficulty arrises from the fact that the inhomogeneities in the scale constrain parameters of the numerical methods employed (say, e.g., the size of the timestep, the discretisation of the mesh) to be small and highly resolved. In the present case, we deal with (highly) oscillatory ODEs. They are omnipresent in the context of Newton's equations 7 and are ubiquitous in the context of Celestial Mechanics. To circumvent the problem, effective models and model reductions techniques are often employed. The core idea is to substitute to the original dynamics a more amenable, numerical and/or analytical, dynamical system (Givon et al., 2004;Lesne, 2006;Hartmann, 2007). One method of choice to design effective dynamics relies on the Lagrangian averaging principle (Grebenikov, 1965;Mitropolsky, 1967;Ghys, 2007;Pavliotis and Stuart, 2008), which has a long-lasting tradition in Celestial Mechanics. This principle is usually used when the components of the equations themselves allow to recognise explicitly the time scales. When it is so, the fast dynamics is integrated into the slow variables to design an averaged approximation. In this setup, the new slow constituents somewhat incorporate the informations of the fast-dynamics and serve as a new input for the investigations. To deal efficiently with our problem at hand, we adopt here our own secular MILAN model, as presented by Gkolias et al. (2016). The MILAN model is based on the vectorial Milankovitch element and admits a minimal force model (consisting of the averaged J 2 contribution, to which is added the secular quadrupolar third-bodies perturbations). The MILAN formulation bears also net advantages compared to the numerical treatment of the Hamiltonian equations in forms of those given in Section 2.1. First, the formulation is free of singularity and, secondly, the averaging is done in a closed form in the eccentricity. The external third-bodies potentials are both averaged over the fast variables of the problem, i.e., over the mean-anomaly of the test particle and the mean anomalies of the third bodies. This 'doubly' averaged model allows the propagation of a test-particle over 10 6 to 10 7 orbital periods in a few seconds only. Such a performance is essential in 6 In a somehow different but connected context, Armellin and San-Juan (2018) have shown that fine discretisations are also needed for optimisers to operate properly.
7 They arrive also in Molecular Dynamics (see Allen (2004) and García-Archilla et al. (1998) for introductory papers and issues).
investigating properties of the phase space for range of parameters. When invoking effective models, we always face the question of the relevance of the reduced model (how sound are the qualitative or quantitative informations derived from it). Gkolias et al. (2016) established the testimony of this doubly-averaged model against a singly-averaged approach. By simulating the two dynamics on different domains of the action-action, action-angle and angle-angle spaces 8 , we showed that dynamical features of interests were reproduced and in perfect agreement. Even if the simulation of the full dynamics (i.e., the original, non-reduced and 'exact' dynamics) on such domains is still missing in the literature, recent reassuring numerical agreements have been presented by Armellin and San-Juan (2018). Namely, they presented nice agreements between their in-house doubly average model and the original non-averaged dynamics. All these together allow us to be confident enough on the numerical results presented hereafter.
3.2. Highly resolved phase-space views. We use the Fast Lyapunov Indicator (FLI), a first-order variational indicator initially introduced by Froeschlé et al. (1997), to discriminate orbit stability. This scalpel has been used extensively over the past decade across different dynamical problems, ranging from Symplectic Maps studies to Dynamical Astronomy, including Astrodynamical practical problems (Froeschlé et al., 2000;Todorović and Novaković, 2015;Guillery and Meiss, 2017;Rosengren et al., 2017). The work of C. Froeschlé, M. Guzzo and E. Lega over the last decade provides a good overview of its possibilities and range of applications. When the dynamics under consideration is written in first order and autonomous form as 9 x " f pxq, x P R n , the FLI is simply derived from the variational system in R 2n , # 9 x " f pxq, Contrarily to Lyapunov exponents, the FLIs (computed at some time τ run for a specific set of initial conditions x 0 , w 0 ) keep trace of the resonant nature of the orbits, while taking approximately the same value FLI " logpτ run q on KAM tori (Froeschlé et al., 2000;Guzzo et al., 2002). By computing the FLIs on a discretised specific 2d-section of ICs (e.g., related to the action-action, action-angle or angle-angle planes) on a domain D, we can reveal the geography of the survival KAM tori and their complement hyperbolic set. The information given by the FLIs ('intensity') is then color-coded to obtain a map of stability. Note that sometimes, in order to get a sharper visualisation, the FLIs that initially take variation in J Ă R`are restricted to a subinterval K of J (see e.g., Lega et al. (2003); Páez and Efthymiopoulos (2015)). This rescaling is achieved by fixing the two following thresholds. The notion of chaoticity is based on the exponential evolution of the norm between two nearby orbits. Therefore, to reveal anomalies with respect to the linear trend (log-scale of an exponential growth), the criteria FLIpτ run q ě logpτ α run q " α logpτ run q, α ą 1, (3.3) can be used to derive a lower threshold for chaotic orbits (i.e., all FLIs larger than α logpτ run q are assigned to α logpτ run q). Symmetrically, we obtain an upper threshold to judge regularity with the criteria FLIpτ run q ď logpτ run {βq, β ą 1.
(3.4) (And again, all the FLIs smaller than logpτ run {βq are assigned to logpτ run {βq.) The Figs. 1 and 2 resume in many ways the transition from order to chaos in the prograde and retrograde region of terrestrial orbits. This original unscrewed fence views of the action-action phase space are very illuminative (and pedagogical) to visualise, with respect to the non-linear parameter εpaq, the proliferation of chaos. Each map represents the result of 1, 000ˆ500 ICs propagated over a long timescale. Different simulation times τ run have been used according to the perturbing parameter (the stronger is the perturbation, the shorter is the time required to get a sharp contrast of the dynamics). For the smallest perturbing parameter, a " 18, 600 km, τ run represents 30 lunar nodes, while for a " 29, 600 km 16 lunar nodes are sufficient to get a sharp contrast (most probably those propagation times could be slightly shortened). It represents about 7ˆ10 5 and 1.8ˆ10 5 test particle revolutions. The 'actions' have been uniformly distributed along the rectangle r50˝, 130˝sˆr0, e max s, with e max determined by the apogee-altitude condition e max " 1´pr C`δ q{a, δ " 120 km. In all our maps, we have set the initial angles y 0 P T 3 to zero. Anticipating a bit the next section, we are here interested in the dynamical mechanisms leading to transport; in particular we have not discriminated collisional orbits as we did in earlier work. As it has already been discussed several times and pointed out in several contributions Gkolias et al., 2016;, the inclination dependent-only resonances widen and develop chaos when ε is increasing, letting less and less room for invariant KAM tori. Eventually for a " 29, 600, there is a macroscopic chaotic component. At this macroscale, we even have the feeling of a chaotic path-connected space (i.e., for every two points in the hyperbolic set, there exists a hyperbolic path connecting them). This property is not exactly true as isolated chaotic islands do exist. Nevertheless, the volume of such isolated chaotic sea is rather small. Let us precise that, given the fact that we used different τ run , the color palette has a symbolic meaning only. Also, in the same way, the z-scale which sets the different levels of the perturbing parameter εpaq is not a linear scale, and again, has only a schematic pictorial purpose.
For the two extreme perturbing parameters considered in this work, εpaq with a " 18, 600 or 29, 600 km, we have superimposed for the newcomers the resonant manifolds obtained under the quadrupolar assumption (confer Eq. (2.12)). It is interesting to notice that, despite the symmetry of the resonant manifolds along the p0, k, 0q resonance, the chaos is not mirrored at all in the retrograde region (see Fig. 3). The coefficients of each harmonic, excepting the critical inclination, are dependent on the cosine of the inclination and hence the resonant topologies for prograde and retrograde orbits are necessarily different. A further striking illustration of this fact, on a microscale, is exemplified in Fig. 4. Such fine resolutions allow to reveal incredible structures and details of the phase space. These two simulations clearly show us, at least for this realisation of angles, that the retrograde region is more stable than its prograde counterpart. Applying the criteria given by Eq. (3.3) with α " 1.1, we found on that domain that the volume of chaotic orbits is 4 times larger than in its retrograde counterpart. We further quantified this question by applying various criteria on our former simulations. Tab. 1 summarises our results by giving the volume of chaotic orbits in the prograde versus the retrograde region, for slightly different values of α on a macroscale. From our survey (which should be extended for completeness), the numerics tend to show that, for small to moderate values of the perturbing values of ε, the volume of chaotic orbits is roughly the same. However, for larger values of ε, the prograde region is more chaotic than its counterpart, the difference being now of several percent. (But again, we are aware of the dependence of our result against the choice of y 0 . Further numerical investigations could constrain even more the result.) At a smaller scale, as already recognisable in Fig. 4, the discrepancies may be largely more significant. It would be interesting to support or invalidate this phenomenology by characterising i) the widths of the resonances of the retrograde domain and ii) by exploring their numerical widths as a function of the angles. Such an enterprise is yet to be performed.

Drift and visualisation of transport
The computation of the FLIs provided a quantification of the degree of hyperbolicity and a discrimination of orbit stability. From a 'practical' perspective, one might be more interested in drift estimation and visualisation of transport to quantify changes of the unperturbed first integrals. This section is devoted to this task by investigating asymptotic properties of initial conditions close or immersed in hyperbolic structures. We base our approaches on individual propagations and on spatial ensemble averages. Figure 1. A highly-resolved fence view of the stability of the prograde and retrograde regions obtained under the FLIs. The three slides depict the stability for a particular value of the non-linearity parameter εpaq which depends on the secularly invariant semi-major axis. KAM tori correspond to white to light-red color, stable resonant orbits appear in blue while red colors correspond to chaotic orbits. The values εpaq correspond to the three semi-major axis t18.6, 22.6, 24.6uˆ10 3 km (the z-scale is only symbolic, in particular the scale is not linear). This unscrewed view presents in a global, original and concise way the transition from order to chaos. See text for comments.

Drift estimation.
There exist a tension between the local degree of hyperbolicity and the eventual large transport. In fact, the astronomical concept of stable chaos teaches us that positivity of a Lyapunov exponent does not necessarily implies large excursion in the phase space (Milani and Nobili, 1992). Large excursions in the phase space can be the signature of transport along the level curves of an integrable system. Nekoroshev's long-time stability theorem does not exclude the existence of chaotic variation. Finally, beyond a critical value, Chirikov's overlap criterion of resonances give rise to large connected chaotic domains, allowing possibly macroscopic transport (Chirikov, 1979).
The problem of chaotic transport (sometimes referred as chaotic diffusion) in nearly-integrable Hamiltonian systems and Dynamical Maps still occupy efforts of various dynamicists (see e.g., Lange et al. (2016); Guillery and Meiss (2017); Cincotta et al. (2018)). Given an orbit computed up to a final time τ run , γptq " `x ptq, yptq˘( 0ďtďτrun , we use the diameter along the action-variables to measure the drift of the unperturbed first-integrals. More precisely, given an action-like vector x P D Ă R n , the diameter D of the orbit is defined as Figure 2. The same as in Fig. 1 apart that the values εpaq correspond to the three semi-major axis t27.6, 28.6, 29.6uˆ10 3 km.
(4.1) For our computations we chose the 8 -norm and computed the drift along the normalised action variables, i.e., along x " pG,Hq withG " G{L andH " H{L (we recall that in the secular approximation, L is a constant parameter determined by the semi-major axis).
The results of the computation of the diameters, according to Eq. (4.1) for two-extreme non-linearity parameters are shown in Fig. 5. Comparing the results with the FLIs maps, we note that the relation between hyperbolicity and large transport is not that straightforward. For a " 18, 600 km, we remark that regions with the larger FLIs do not necessarily correspond to regions where the transport is maximal. Conversely, the almost vertical resonant manifold emanating near i " 56.1˝does not have the largest degree of hyperbolicity; yet it carries the largest transport index. Switching to a " 29, 600 km, we note that the lowest diameter is already one order of magnitude larger than in the former case. The largest diameter is also significantly larger which confirm the known fact of the instabilities in the MEOs. We emphasise that the diameters have been computed on the same predefined grids of ICs used to estimate the FLIs (i.e., a highly-resolved grid of ICs). The emanating feeling of a resolution deterioration in the maps is once again a nice testimony of the sensitivity of variational indicators.
For large perturbing parameters, globally speaking, large hyperbolicity corresponds to large diameters. This fact has to be nuanced slightly near e " 0.7 and i " 70˝. Using an empirical criterion, we extracted from the maps the actions that satisfy the condition FLIpx, τ run q ě 1.2 logpτ run q (i.e., chaotic orbits) as those satisfying Dpx, τ run q ě 0.35. The tracing orbits are shown in Fig. 6 and illustrate the  link between large hyperbolicity and large diameters, and the necessity of finely resolved meshes (thin stable structures stripe the chaotic domains and can be detected with the diameters also).
Let us now comment on the diameter indicator that we used. Very often diameters-like quantities in terrestrial dynamics have been estimated using a more restrictive definition, namely a one-dimensional Figure 5. Estimation of the diameters for the two perturbing parameters εpaq with a " 18, 600 km and a " 29, 600 km. Figure 6. Extraction of the ICs satisfying Dpx, τ run q ě 0.35 (left hand-side) and the chaotic ICs satisfying FLIpx, τ run q ě 1.2 logpτ run q. Figure 7. Ensemble integration of a cluster of k " 200 orbits in a neighborhood Vpx ‹ q of x ‹ " p0.616, 88˝q. The ensemble averaged orbit of the considered observables are shown in red. The eccentricities do not experience a net sensitive dependence to the ICs, contrarily to the inclinations. From this example, it can be easily inferred that a diameter-measure based solely on the eccentricity (or equivalently on G) would fail to capture properties of the dynamics. 0.14 0.08 Table 1. The estimation of the volume of chaotic orbits in the prograde and retrograde regions, for various perturbing parameters and on various domains. The domain D refers to the definition of the domain in the prograde region, in the eccentricityinclination action phase space. This domain is then mirrored in its retrograde counterpart to serve as a new domain to determine the volume of chaotic orbits in the retrograde region, Vά . The Eq. (3.3) is used as a discrimination criteria. All results have been established with a fine mesh (all domains have been uniformly discretised with a grid consisting of at least 500ˆ500 initial conditions). The prograde region appears to be slightly more chaotic than the retrograde counterpart on a macroscale the more we increase the perturbing parameter. Significant differences may also exist at smaller scales.
diameter of a specific observable f (see, e.g., Alessi et al., 2016;Rosengren et al., 2017). This strategy reduces to nothing else than the amplitude estimation, equivalent to the estimation of ∆f " max t f pxq´min t f pxq. For the MEO problem, the eccentricity diameter along the time is, rightly, tracked (as efforts are directed towards the perigee height and the need of re-entry solutions). However, when used as an empirical 'measure of chaos', this diameter may be too loose. In fact, having in mind the geography of the resonant manifolds derived from the resonant condition in Eq. (2.12) and the fact that two actions characterises an invariant torus of T 3 , it is easy to 'create' a quasi first-integral by choosing ICs near certain manifold. As an example, let us fix a " 29, 600 km and consider a cluster of ICs in a small neighbourhood of Vpx ‹ q, where x ‹ " pe ‹ , i ‹ q " p0.616, 88˝q. The time evolution (over 25 lunar periods) of the eccentricity and inclination for the whole cluster of orbits (k " 200 orbits) is displayed in Fig. 7. The spatial averaged orbit is displayed and superimposed with a bold red line 9 . Clearly, the eccentricities of the whole cluster evolve in an apparent regular fashion. All the orbits incorporate similar dynamical informations, both on the quantitative and qualitative point of view. On the contrary, the inclination time-histories experience significant variations and a net sensitive dependence upon the ICs. From this example, easily generalisable, we easily infer why a one-dimensional diameter (based on the eccentricity) would fail in capturing these particularities.
Pushing further the idea, we extended this approach on a grid of ICs near the point x ‹ by computing accordingly the diameters (and the FLIs). The obtained maps are presented in Fig. 8. They confirm the rationale behind the intuition developed through the former example. Whilst the diameter based 9 Let us consider a cluster of size k, that we propagate up to time τrun. We obtain k orbits γ k P`Cr0, τruns˘. Le us denote by x j pi, tq the instantaneous value of the j-th component of the orbit γ i at a specific epoch t. The spatial averaged orbit of the component j (1 ď j ď n), xx j y, is then defined through its components obtained at any time t by xx j ptqy " 1 k ř k i"1 x j pi, tq. on both actions is in agreement with the FLI map, the method based on the one-diameter approach give an irrelevant and uniform signal.
Having presented a general way to quantify the drift, let us focus now on how the drift is mediated in the phase space.

Visualisation of transport. In the previous sections, we computed FLIs and diameters in various sections
Spvq " px, yq P DˆT 3 | y " v, v P T 3 ( (4.2) with D Ă R 2 . By fixing y " 0, particular features in Sp0q have been depicted. In order to visualise transport properties, and to show how its mediation is related to the detected hyperbolic web, we use projection and visualisation techniques that have been extensively used over the past decade to study transport in nearly-integrable Hamiltonian system, symplectic Maps and in Dynamical Astronomy (Guzzo et al., 2002;Lega et al., 2003;Cincotta and Giordano, 2008;Páez and Efthymiopoulos, 2015;Guillery and Meiss, 2017). For a recent overview specifically around the FLIs and their applications, we advise the reader to consult Lega et al. (2016) for a pedagogical introductory note. The methodology consists in the following. First, we compute the FLIs over a section Spvq, say on Sp0q 10 . After this step, we are then able to recognise initial conditions close to hyperbolic borders or immersed within the chaotic sea. We then select one IC of interest in Sp0q. Let x ‹ P Sp0q denotes this IC. Next, we define a small neighbourhood Vpx ‹ q of k ICs of x ‹ . In theory, it would be sufficient to deal with the sole numerical propagation up to τ run of the orbit emanating from x ‹ . However, the procedure is computationally facilitated by considering a cluster of k orbits. From these computed orbits γ k ptq P C`r0, τ run s˘, we keep trace only of the points that return close enough to the section Spvq. For that purpose, we introduce the family of sections tS δ pvqu δ which are δ-close to Spvq. These sections are defined as When δ Ñ 0, we recover the 'exact' section Spvq. The introduction of this family of section is essentially to circumvent numerical limitations. Firstly, we deal with a finite time τ run (that we would like to keep 'as small as possible' but 'large enough' to extract dynamical mechanisms). Secondly, in practice we do not deal with an orbit γ k ptq P C`r0, τ run s˘, but with a discretised version of this orbit computed, say (to facilitate the exposition), at each multiple of the fixed step size ∆t, tγ k ptq, t " i∆tu n i"0 , n∆t " τ run . All points of the orbits γ k ptq P C`r0, τ run s˘that return during the simulation to a section of tS δ pvqu δ are identically projected into the exact section Spvq, on which the FLIs are used as a background. By doing that, we are able to relate transport with the web detected by the FLIs. In our computation we dealt with the 8 -norm, δ is problem dependent and best determined by a calibration procedure 11 . Finally we worked with a cluster of size k " 200 initial conditions. Fig. 9 presents results in the range of 'small' perturbation for two initial points of interest applying the methodology described previously. The clusters has been propagated up to a timescale of about 5.8ˆ10 5 orbits revolution (25 lunar nodes). The ICs serving a definition to the cluster are depicted in red. The points of the orbits that cross the double sections of the set tS δ u δ are depicted in green. (To facilitate the reading and interpretation of the figures, the FLIs background have been color coded with a light opacity-like filter. The points returning to the section are intentionally magnified.) The two clusters focus on thin manifold that still carry transport (see Fig. 5). As the transport index is rather small, excursions are modest and rather confined. The returning points are guided by the thin hyperbolic skeleton detected by the FLI computation.
In Fig. 10 we repeated the same experiment with a larger ε for 4 different scenarios. The approach enables us to visualise and quantify the spread of the actions in the regime of strong chaos. The orbits of the clusters have been propagated on about 1.4ˆ2.8ˆ10 5 orbits revolution. The spread of the orbits is well more appreciable and develops more drastically within the action-space. It covers a large portion of the connected chaotic domain. As it is observed for all scenarios, the change in inclination can be superior to 15˝, with extremely large variations for the eccentricity (namely, the mechanism allows nearly circular orbits to become very eccentric).
It would be extremely interesting to extend the approaches and visualisation of the diffusive properties by extending the dimension of the visualised space. Taking advantage of our model, we were able to extend the traditional stability maps in one more direction by stitching together ad-hoc others FLI sections. The results presented in Fig. 11 complement the global stability picture of the actions space by 'unrolling' the dynamics according to one angle, here Ω. The resonant manifolds computed using Eq. (2.12) are depicted in black in the "action-space". In the regime of small perturbation (left panel), a pendulum-like structure is clearly identifiable (minor structures can also be identified). By varying the size of the perturbation, a bifurcation-like phenomena occurred and the initially elliptic point becomes of a hyperbolic nature where collisional orbits develop. Such a systematic parametric 10 In this work we were interested in the action-action plane, but the approach can be extended to action-angle or angle-angle planes. For example, a angle-angle section can be defined as T " px, yq P DˆT 3 | py 1 , y 2 q P B Ă T 2 , x P D, y 3 " v 3 ( . 11 To give an idea of the size of δ, the results presented in this manuscript have been obtained with δ " 0.08 for the small range of ε, δ " 0.1 for larger range. Different admissible δ just change the number of points on the section collected, but leave invariant the transport properties (angles are expressed in radians). methodology would allow, besides the quantification of chaos and the determination of the resonant regime (cf. Froeschlé et al., 2000), the determination of precise perturbing parameters where such phenomena occur.

Discussion and conclusive remarks
Dynamical chaos indicators as the FLI are valuable and formidable allies to gain knowledge on the dynamical system under investigation. Their systematic use over nearly the past two decades in transverse fields has brought its share of results. Applications towards terrestrial dynamics are still at their early stage but the current situation seems to evolve positively. In this contribution, we complemented and refined our past studies related to the long-term dynamics of terrestrial orbits in the range 2.91 to 4.64 Earth radii (ε P r0.02 : 0.22s). We showed the complementarity and benefits of visualising the global dynamics via sections, corroborated with the computation of the FLIs and practical action-diameter quantities. From our numerical experiments, we have seen that when the detected hyperbolic manifolds are very thin (but still carry large diameters), the transport occurs precisely along them. For higher values of the non-linearity parameter, resonances do overlap significantly and the transport is across a large domain of the chaotic sea. This mechanism allows nearly-circular orbits to become highly eccentric on a few lunar nodes only. In the later case of strong chaos, preferred directions for the transport are hard to establish. The FLIs allow to follow and delineate the routes of transport where the spread in the phase space take place. The natural complementary step that deserves serious attention concerns the nature of the transport, the computation of diffusion-like coefficients and its scaling with ε. (Note that if in our actual set-up, we do have access only to a limited number of different order of magnitudes of ε. A theoretical possibility to extend its range is to artificially increase the semi-major axis -even if we know that physically the procedure is not that relevant as octupolar contributions should be incorporated.) Let us comment and relate recent difficulties that we encountered in investigating these last points. Transport properties are generally characterised through the computations of moments of different order q, M q pτ q " @ |xpτ q´xxpτ qy| q D . (5.1) Let us underline once more that when we deal with the dynamics numerically, we only have access to finite time moments. Usually, the second-order moment, i.e., the spread of the actions (the variance), is used to discriminate the case of diffusion we deal with. More precisely, under the explicit ansatz Figure 10. Diffusive scenarios illustrating transport scenarios within the hyperbolic web in the regime of large perturbative parameters and strong resonances overlap. The red spots indicate the initial condition where the ensemble of initial conditions are defined. Points of the orbits that return sufficiently close to the section (on which the FLIs appear as a background) are depicted with a green point. See text for comments. Figure 11. These two FLI cubes, computed for a " 18, 600 km (left) and a " 29, 600 km (right) highlight a bifurcation phenomena. The resonant manifolds appear in black in the 'action-action' space. that M 2 pτ q " @ |xpτ q´xxpτ qy| 2 D " D 2 τ ν , (5.2) the diffusion is called either subdiffusive (ν ă 1), diffusive (ν " 1) or superdiffusive (ν ą 1). (The particular case of superdiffusive behaviour with ν " 2 is referred to ballistic diffusion.) The real parameter D 2 is the estimated diffusion coefficient, and its sole determination can be sometimes tricky due to technical difficulties (see, e.g., Lega et al. (2003) and further references in Cincotta et al. (2018)). Anomalies to the strict diffusive case (ν " 1), i.e., aberrations with respect to Gaussianity, might be the results of the existence of a mixed phase space (cohabitation of regular and chaotic components in the phase space) and correlation effects (Varvoglis et al., 1997;Zaslavsky, 2002). Let us note that, to the best of our knowledge, the study of the correlation function Cpτ q (even at least for the specific observable of interest, the eccentricity) and its possible decay which give us the scale of the correlation time τ C (see discussions in Wiggins and Ottino, 2004;Varvoglis, 2005) has never been undertaken for the MEO problem. (The exception is found in Wytrzyszczak et al. (2007) where the autocorrelation function properties are used to discriminate regularity for geosynchronous objects.) Daquin et al. (2017) claimed the normal character of the diffusion for the eccentricity observable in the regime of strong chaos. We redid some experiments along those lines apart that we used the spatial averaging ideology (and no longer the temporal averaging) assuming that all ICs of the cluster are equivalent. We met difficulties to confirm our former conclusions and we stress here that they should be taken with a grain of salt. In fact, in our experiments, we noticed that such a conclusion depends strongly on the ansatz made on the evolution of the variance and the choice of the time-horizon investigated. Regarding the question related to the time-horizon, there might exist a transient time τ tr. that should be constrained first. Indeed, in order to derive meaningful statistical conclusions, we have to ensure that τ " τ tr. (as a transient time seems to exist) and τ " τ C . It is possible that, unfortunately, in our present setting, τ tr. " τ , making conclusions hard to reach.
Constraining those difficulties are the directions being taken by our current research.

Author Contributions
The paper has been written by the authors in equal parts.