Semi-analytical estimates for the orbital stability of Earth's satellites

Normal form stability estimates are a basic tool of Celestial Mechanics for characterizing the long-term stability of the orbits of natural and artificial bodies. Using high-order normal form constructions, we provide three different estimates for the orbital stability of point-mass satellites orbiting around the Earth. i) We demonstrate the long term stability of the semimajor axis within the framework of the $J_2$ problem, by a normal form construction eliminating the fast angle in the corresponding Hamiltonian and obtaining $H_{J_2}$ . ii) We demonstrate the stability of the eccentricity and inclination in a secular Hamiltonian model including lunisolar perturbations (the 'geolunisolar' Hamiltonian $H_{gls}$), after a suitable reduction of the Hamiltonian to the Laplace plane. iii) We numerically examine the convexity and steepness properties of the integrable part of the secular Hamiltonian in both the $H_{J_2}$ and $H_{gls}$ models, which reflect necessary conditions for the holding of Nekhoroshev's theorem on the exponential stability of the orbits. We find that the $H_{J_2}$ model is non-convex, but satisfies a 'three-jet' condition, while the $H_{gls}$ model restores quasi-convexity by adding lunisolar terms in the Hamiltonian's integrable part.

1. Introduction One of the major goals of Celestial Mechanics is the analysis of the stability of the dynamics of celestial bodies. Knowing the behavior in time of the orbital elements of an object allows one to predict its future, in particular whether it will cross the orbit of other celestial bodies and eventually undergo collisions. When applied to artificial spacecraft and space debris orbiting around the Earth, the question of the stability becomes of crucial importance, especially in view of the problem of estimating the orbital survival times of operating satellites or space debris. It is therefore crucial to devise methods that allow to study the orbital stability of objects moving around our planet.
In this work we will not consider the complex dynamics of an artificial spacecraft, which should include the analysis of its shape, composition as well its rotational motion. We will rather consider a point-mass body around the Earth, that we can identify with one of the several millions of space debris orbiting our planet. In fact, the proliferation of artificial satellites in the last decades has led to the generation of an enormous amount of space debris with different sizes, from meters down to microns, and at different altitudes. Space debris are remnants of non operational satellites or the result of break-up events, either collisions or explosions. Since the altitude determines the contribution of the different forces acting on the object (the gravitational attraction of the Earth, its geopotential perturbation, the influence of Sun and Moon, the Solar radiation pressure, etc.), it is convenient to make a distinction in terms of the altitude. To this end, the space in the surrounding of the Earth is commonly split into three main regions: LEO ('Low Earth Orbit') denotes the region up to about 2 000 km of altitude in which the Earth's attraction, the geopotential as well as the atmospheric drag are the terms which greatly affect the dynamics of an Earth's satellite; MEO ('Medium Earth Orbit') refers to the region between 2 000 and 30 000 km in which the effects of Moon, Sun and Solar radiation pressure become important; GEO ('Geosynchronous Earth Orbit') refers to a thin (∼ 200 km) zone around the geostationary orbits (at 42 164 km from Earth's center), where the satellites are in synchronous resonance with the 24-hour rotation of the Earth around its spin-axis.
The huge amount of objects (up to millions) in LEO, MEO, GEO needs a careful control of their orbits and the analysis of their dynamical stability ( [6,3,5,4,15,31]), also in view of devising appropriate mitigation measures (see, e.g., [17,27,9]). For objects in LEO, it is of crucial importance to evaluate the orbital lifetime, which is strongly affected by the atmospheric drag which provokes a decay of the orbits ( [19,22,24,32,34]). In this work we focus on objects in MEO, GEO and beyond, thus not taking into account the dissipative effect of the atmosphere. Instead of using a propagation of the orbits to predict the stability time of the orbital elements, we propose a procedure based on analytical perturbative methods (see also [13]). More precisely, via a suitably defined sequence of canonical transformations, we construct a normal form of the Hamiltonian function, which enjoys the property that one or more of the Hamiltonian's Delaunay actions define quasi-integrals of motion (namely, integrals of the integrable part of the new Hamiltonian). Once the transformed Hamiltonian is obtained, the size of its remainder (which gives a control on the goodness of the approximation) can then be used to provide bounds on the time variations, and hence the stability time of the orbital elements (semimajor axis, eccentricity, inclination) as a function of the distance of the object from the Earth. We refer to this procedure as semi-analytical, which means that it uses an analytical method, precisely normal forms, whose coefficients are calculated numerically, namely with the aid of a computer.
We consider two different models to describe the motion of the debris around the Earth. The first model takes into account only the influence of the geopotential up to the term J 2 of its expansion in spherical harmonics; we refer to this problem as the J 2 model and denote the corresponding Hamiltonian as H J 2 , which results from truncating to a suitable power of the coordinates around reference values, and normalizing up to a suitable order, as described in Section 4.1. The second model, referred to as the secular 'geolunisolar' model (Hamiltonian H gls , truncated and normalized similarly to H J 2 , see Sections 2.2 and 5.1), includes also the effects of the Moon and the Sun, placing, for simplicity, the Moon strictly on the ecliptic; this last restriction means to omit from the Hamiltonian terms corresponding to lunisolar resonances other than the 'inclinationdependent' ones. The latter resonances, on the other hand, are those producing the most important effects as regards orbital stability (see [4,3] for a review). Furthermore, instead of formally eliminating the fast angle via canonical transformations (as we do in the pure J 2 problem), in H gls we just take the average of the Hamiltonian with respect to all fast angles, namely, the mean anomaly of the satellite as well as the fast angles of the Moon and Sun: this averaged model allows us to focus on the satellite's long-term dynamics (i.e. the secular one). The averaging is done in closed-form and leads to formulas equivalent to those described in [18]. Furthermore, we reduce this last Hamiltonian to actionangle variables around each forced equilibrium point, which corresponds to a non-zero inclination defining the so-called Laplace plane (see Section 2.2).
In summary, our stability estimates are obtained according to the procedure (i)-(iii) outlined below: i) Within the J 2 model, we make a formal elimination in the Hamiltonian of the fast angle (mean longitude); as a consequence, we get the preservation of the conjugate action variable corresponding to the semimajor axis. This allows us to compute the stability time for the semimajor axis at different altitudes, yielding stability times that increase with the altitude. ii) Using H gls , instead, the semimajor axis becomes a parameter (with a priori constant value), while we proceed to analyze the behaviour of eccentricity and inclination. The latter is obtained using a quasi-resonant normal form, which reflects the 1:1 near-resonance of the integrable part of H gls between the frequencies of the longitude of the ascending node and the sum of the argument of perigee and the longitude of the ascending node (see Section 2.2.1). This means that, close to the Laplace plane, the quasi-preserved secular quantities cannot be defined neither as the eccentricity e nor the inclination i alone, but rather by the Kozai-Lidov combination I = 1 − √ 1 − e 2 cos i ≈ (e 2 + i 2 )/2 (for e, i small). We then explore the dependence of the stability time of I on the altitude of the orbit. Our results show that the J 2 and lunisolar terms have an opposite effect on the time of stability as the distance from the Earth increases. As a by-product of this analysis, we also compute the so-called forced inclination (that is, the inclination of the Laplace plane), which corresponds to the shift of the secular equilibrium from a strictly equatorial orbit to an orbit with small positive initial inclination, an effect caused by the fact that the perturbing bodies (Moon and Sun) are in orbits inclined with respect to the Earth's equator. iii) Finally, as a first step towards obtaining exponential stability estimatesà la Nekhoroshev ( [26]), we check whether some so-called 'steepness' conditions are satisfied for the integrable part of both Hamiltonians H J 2 and H gls , namely whether the integrable parts are convex, quasi-convex, or satisfy the three-jet condition (see [8] and references therein). The results show that the J 2 model is three-jet non-convex, while the contribution of the lunisolar terms removes the intrinsic degeneracy of the J 2 part and allows us to conclude that the geolunisolar model is quasi-convex. A detailed application of the non-resonant form of Nekhoroshev's theorem in the Hamiltonian H gls is the subject of an independent paper (see [2]).
Summarizing, the previous strategy allows us to obtain three different stability results: one for the semimajor axis in the J 2 model, a second for the stability of the eccentricity and inclination in the geolunisolar model, and a third on the holding, altogether, of necessary conditions for Nekhoroshev-type stability of the satellite orbits. All three results point towards the same direction, i.e. that, at least far from exact resonances, orbital stability can be ensured at MEO, GEO and beyond for quite long times (10 4 − 10 6 orbital periods, 10 2 − 10 4 years). Besides these general numbers, one may remark that the calculation of the size of the remainder of the normal form actually provides an estimate of the rate of drift of the orbits in element space, an information required in orbital diffusion studies for defunct satellites and space debris.
This work is organized as follows. In Section 2 we present the J 2 and geolunisolar models. Section 3 briefly presents the method of the composition of Lie series, used in the computation of all our normal forms, along with some general estimates on the convergence of the normalizing transformation and the size of the normal form's remainder. Section 4 focuses on the stability estimates with the H J 2 model, while section 5 deals with the stability in the framework of the H gls model. Finally, the analysis of the steepness conditions for the Hamiltonians H J 2 and H gls is presented in Section 6.

The J 2 and geolunisolar models
Bodies orbiting around the Earth are primarily affected by the Keplerian attraction with our planet. However, for an accurate description of the dynamics it is mandatory to assume that the Earth is non-spherical. Beside the Earth, the satellite dynamics is subject to the gravitational influence of Sun and Moon. Section 2.1 describes the Hamiltonian model H J 2 kep , which includes the Earth's Keplerian term and the first nontrivial term in the expansion of the geopotential. Section 2.2 presents the Hamiltonian model H gls,sec , which includes J 2 and lunisolar terms, averaged over the fast angles.
2.1. The J 2 model. We consider a model describing the motion of a point-mass body, say a satellite S, under the effect of the Earth's gravitational attraction, including an approximation of the geopotential due to the oblateness of the Earth. Let r ≡ (x, y, z) be the position vector of S in a geocentric reference frame, with the plane (x, y) coinciding with the equatorial plane, and x pointing towards a fixed celestial point (e.g. the equinox). We consider the Hamiltonian describing the motion of S under the geopotential as the sum of two terms where is the Keplerian part (r = |r|), and is the J 2 potential term, arising from expanding the geopotential in spherical harmonics and retaining only the largest coefficient (see, e.g., [18]). The constants are the Earth's mass parameter µ E = GM E (G = Newton's constant, M E = Earth's mass), R E is the mean Earth's radius, and J 2 is a dimensionless coefficient describing the oblate shape of the Earth. The numerical values are: • µ E = 1.52984 × 10 9 R 3 E /yr 2 ; • R E = 6378.14 km; The Hamiltonian (1) is expressed in Cartesian coordinates. However, by a standard procedure, it can be transformed to an expression in the following set of modified Delaunay canonical action-angle variables where (a, e, i, M, ω, Ω) are the orbital elements of the satellite (semimajor axis, eccentricity, inclination, mean anomaly, argument of the perigee, longitude of the nodes). The passage is done by first expressing the Hamiltonian (1) in elements via the relations (see e.g. [25]) where f is the true anomaly and r, cos f and sin f are given by the series In the actual calculations, all series are truncated to order N = 15 in the eccentricity e. Finally, we pass from the elements (a, e, i, M, ω, Ω) to the canonical variables (L, P, Q, λ, p, q) by inverting Eqs. (4).
To perform the high order normal form computations described in Section 4, using computer algebra, it is convenient that the dependence of the Hamiltonian on the actionangle variables be expressed as a trigonometric polynomial. To this end, we first make a shift transformation L → δL around a reference value a * , with This means that the Hamiltonian found after expanding in powers of the quantity δL refers to the local dynamics of orbits with semimajor axis a ≈ a * . Every time when we change the reference value a * (i.e. the 'altitude' or 'distance' of the orbit from the Earth's center), we then perform the Hamiltonian expansion anew around L * and obtain the stability estimates corresponding to that reference value. One may also note that P = O(e 2 /2) and Q = O(i 2 /2), thus all three quantities δL, P and Q are small quantities for orbits not very far from the equator and not very far from circular. We then expand H J 2 kep (δL, P, Q, λ, p, q) in powers of √ δL, √ P , and √ Q up to the same order N = 15 as the original expansion in the eccentricity (this ensures missing no term in P, Q in the Hamiltonian up to the order N ). After this change, the truncated Hamiltonian takes the form (apart from a constant): k 1 ,k 2 ,k 3 ∈Z 0<|k 1 |+|k 2 |+|k 3 |≤s P s,k 1 ,k 2 ,k 3 (δL, P, Q) cos(k 1 λ + k 2 p + k 3 q).
The functions Z s and P s,k 1 ,k 2 ,k 3 (δL, P, Q) are polynomials of degree s/2 in the action variables (δL, P, Q). Finally, the frequencies n * , ω 1 * , ω 2 * are equal to: The Hamiltonian (8) is the point of departure for the stability estimates on the orbits' semimajor axes; one notices that ω 1 * = −ω 2 * = O (J 2 ), a fact implying that both these frequencies are way smaller than n * (µ E /a 3 * ) 1/2 (third Kepler's law). Accordingly, for all orbits the angle λ circulates at a rate which is O(1/J 2 ) faster than the rate of circulation of the angles p, q. Hence, λ constitutes the 'fast angle' of the Hamiltonian H ≤N J 2 . Its elimination through a suitable sequence of canonical transformations leads to the approximate constancy of the value of the semimajor axis, as detailed in Section 4.
2.2. The geolunisolar Hamiltonian. While stability estimates for the semimajor axis depend mostly on the Earth's J 2 term, the question of the long-term stability as regards secular variations in eccentricity and inclination requires considering the effects of the Lunar and Solar gravitational tides. Let us consider a celestial body B (either Moon or Sun) with mass M b moving around the Earth and whose orbit is exterior to that of the satellite. Let r = (x, y, z) and r b = (x b , y b , z b ) be the position vectors of S and B in a geocentric reference frame, with r = |r| and r b = |r b |. The tidal disturbance caused by B on S is described by the potential where µ b = GM b . The first term −µ b /r b in the multipolar expansion (10) does not depend on the coordinates of S, therefore it can be omitted from the Hamiltonian of motion of S. Thus, the tidal (or 'third body') perturbation terms in the Hamiltonian takes the form: (r · r s (t)) 2 r 5 s (t) where µ m , r m and µ s , r s are the mass and geocentric position vectors of the Moon and Sun respectively, and O 3 denotes octupolar or higher order terms in the expansion of the third body potentials. The exact form of the term H 3B depends now on the model adopted for the geocentric orbits of the Sun and Moon. In the framework of the present paper, we adopt the following models for Sun and Moon: (1) we suppose that the Sun's geocentric orbit is an ellipse lying in the Earth's ecliptic plane (i.e., with inclination i s0 = 23.43 • with respect to the equatorial plane), Ω s = 0 • , a s = 1.496 · 10 8 km and e s = 0.0167; (2) we assume the Lunar orbit as elliptic and also lying on the ecliptic plane, with a m = 384748 km, e m = 0.065 and i m0 = i s0 . Note that this assumption ignores the precession of the Lunar node (with period 18.6 yr) associated with the inclination of the Moon's orbit with respect to the ecliptic (by 5 • 15 ). While the precession of the Lunar node is important near secular lunisolar resonances 1 , it only has a minimal effect far from these resonances, as substantiated by numerical studies (e.g. [16], [30]). Thus, we ignore this effect in our present estimates (Section 6).
Under the above approximations, the satellite Hamiltonian H J 2 ls takes the form This is a Hamiltonian depending on three degrees of freedom (the coordinates and momenta of the satellite) as well as on time (through the vectors r m (t) and r s (t)). However, contrary to the case of the Hamiltonian H J 2 kep , in which we are interested in establishing the long-term stability of the semimajor axis over short-period oscillations, here we are interested in the question of the stability of the eccentricity and inclination of the satellite over secular timescales. Thus, as customary (see [18], [7]), we average H J 2 ls with respect to the mean anomalies of the satellite, Moon and Sun. The averaging can be done in closed form (see, for example, [18]), and leads to: Here, f, E are the satellite's true and eccentric anomaly, while f m , f s are the Moon's and Sun's true anomaly along their geocentric orbits. The averaged J 2 term takes the form (apart from constant terms): The terms H (e, i, ω, Ω), instead, turn out to be identical to those given in equations (3.6) and (3.7) of [3], setting i M = 0. Then, the Hamiltonian averaged over all short period terms, hereafter referred to as the secular geolunisolar Hamiltonian, takes the form which, in terms of the Delaunay modified variables, has two degrees of freedom.
2.2.1. Expansion around the forced inclination. As it was done in the case of the J 2 model (Eq. (8)), normal form computations for the Hamiltonian (14), expressed in Delaunay action-angle variables, require a polynomial expansion in the action variables around some preselected values. In the case of the secular geolunisolar Hamiltonian (14), a natural choice of the origin for such expansions is the forced element values: writing H gls,sec as a function of the Delaunay variables, say, H gls,sec (P, Q, p, q; a) (where the semimajor axis a is now a priori constant, hence, can be considered as a parameter), a forced equilibrium is defined as an equilibrium point of the secular Hamiltonian, i.e., a point (Q (eq) , P (eq) , q (eq) , p (eq) ) for which the following relations hold: where the subscript 'eq' denotes the condition Q = Q (eq) , P = P (eq) , q = q (eq) , p = p (eq) .
In the case of the Hamiltonian (14), a forced equilibrium solution can be computed by writing first H gls,sec in terms of the Poincaré variables as Expanding up to quadratic terms in the Poincaré variables, the truncated secular Hamiltonian has the form where the coefficients A 1 , B 1 are given by: The Hamiltonian (17) corresponds to two decoupled harmonic oscillators in the variables (X 1 , Y 1 ) and (X 2 , Y 2 ). The second harmonic oscillator (corresponding to the actionangle pair (P, p), hence, to the orbit's eccentricity vector) has an equilibrium point at X (eq) 2 , Y (eq) 2 = (0, 0), implying P (eq) = 0 and any value 0 ≤ p (eq) < 2π. This implies that the sub-manifold of circular orbits e = 0 (corresponding to P = 0) is invariant under the flow of the secular geolunisolar Hamiltonian. On the other hand, as regards the pair (X 1 , Y 1 ), Hamilton's equations for the Hamiltonian (17) yield: For i 0 = 0, the equilibrium point of (19) is given by More accurate expressions for the forced inclination i (eq) can be obtained by introducing (20) along with the remaining equilibrium values in the derivatives of the full secular Hamiltonian (14) and finding the roots of Hamilton's equations. One can readily verify that q (eq) = 0 at all orders, while i (eq) is subject to small corrections with respect to the expression (20). In physical terms, the forced inclination i (eq) defines the inclination of the Laplace plane: since the perturbing bodies (Moon and Sun) are in orbits inclined with respect to the equator, a satellite orbit can maintain its inclination constant when the latter has the value i (eq) . Inspecting the form of the coefficients (18), we find that i (eq) → 0 as a → 0, while it can be shown that i (eq) → i 0 for values of a greater than the GEO one (see for example [29]), reflecting the fact that the Laplace plane tends to coincide with the equator for satellite orbits close to the Earth (as imposed by the oblateness of the Earth), while it tends to coincide with the ecliptic at large distances from the Earth (where the Lunar and Solar tides dominate).
Returning to the expansion of the secular geolunisolar Hamiltonian, making the shift allows us to express the Hamiltonian as a polynomial in the variables (X 1 , δY 1 ) and (X 2 , Y 2 ). The Hamiltonian H gls,sec starts now with terms of second degree which we regroup in H 2 : where b 1 = 2B 1 and 1 , 2 , 3 , 4 are corrections of order O(µ b a 3/2 /(µ , with the index b referring to the Moon or Sun. All these corrections turn to be rather small, with relative size ∼ 10 −3 B 10 at semimajor axis a ∼ 10 4 km, where Thus, after a canonical rescaling , the secular lunisolar Hamiltonian resumes the form: This is the typical form of a secular Hamiltonian, consisting of linear oscillators (with frequencies ν 1 , ν 2 ) coupled with nonlinear terms. However, we have implying that the two frequencies are nearly equal This is a consequence of the axisymmetry of the J 2 model, implying that the secular frequenciesq = −Ω andṗ = −ω −Ω are equal for nearly equatorial orbits in this model. As we will see in Section 5, this near-equality implies that with the present normal form estimates one cannot establish independently the long term stability of the eccentricity and the inclination, but only the long-term stability of the Kozai-Lidov integral I =X 2 1 +Ỹ 2 1 +X 2 2 +Ỹ 2 2 , which couples oscillations between the eccentricity and the proper inclination of the satellite.
Finally, the Hamiltonian (22) can be written in action-angle variablesX The Hamiltonian (23) is the starting point for all normal form calculations in Section 5. For computational reasons, the expansion in (23) is truncated up to a maximal order N = 15, leading to the truncated form

Hamiltonian Normalization
In this Section we briefly recall some basic definitions related to normal form theory and its use in obtaining stability estimates based on the size of the normal form's remainder. In Sections 4 and 5 we will discuss the particular normalizations implemented on the Hamiltonians (8) and (24) respectively.

Normal form and remainder. Consider a Hamiltonian function of the form
where ω j are real constants, and (A, ϕ) ∈ R n × T n are action-angle variables. We assume that S σ is the complex strip for ρ, σ > 0. On D ρ,σ (U ) we define the norm of a function f = f (A, ϕ) as The aim of normalization theory is to introduce a near to identity canonical transformation Φ : (A, ϕ) → (A , ϕ ), so that in the new variables (A , ϕ ) the Hamiltonian (25) takes the form with the following properties: i) the transformation Φ is analytic in a domain D ρ ,σ (U ) with 0 < ρ < ρ, 0 < σ < σ, ii) the dynamics under Z(A , ϕ ), called the normal form, has some desired properties (see below), and iii) under the norm definition (28) one has R ρ ,σ Z ρ ,σ implying that the function R(A , ϕ ), called the remainder, introduces only a small correction with respect to the flow under the normal form term Z(A , ϕ ).
Regarding point ii) above, see, e.g., [10] for a definition of the properties of the normal form term in various contexts of perturbation theory (e.g. in the Kolmogorov-Arnold-Moser or Nekhoroshev theories). Here we mention three cases of particular interest, pertinent to our present work: Case 1: Birkhoff normal form. The function Z is independent of the angles ϕ . This kind of normalization allows us to prove the near-constancy of the action variables A.
Case 2: elimination of short-period terms. The real constants ω j in (25), called the unperturbed frequencies, are divided in two groups, 'fast' {ω 1 , . . . , ω K f }, 1 ≤ K f < n, and 'slow' {ω K f +1 , . . . , ω n }, such that min{|ω 1 |, . . . , |ω K f |} max{|ω K f +1 |, . . . , |ω n |}. In this case, it turns convenient to introduce a normalizing transformation Φ such that the normal form Z becomes independent of the 'fast angles' {ϕ 1 , . . . , ϕ K f }. Such is the case of the normal form encountered in Section 4, leading to estimates on the stability of the semimajor axis in the J 2 problem. The corresponding Hamiltonian is of the form (25), with n = 3, Case 3: resonant normal form. The frequencies ω j satisfy one or more quasi-commensurability relations of the form m · ω 0, with m ∈ Z n , |m| = 0. The maximum number of linearly independent and irreducible integer vectors m l , 1 ≤ l ≤ l max , yielding exact commensurabilities for a given set of frequencies ω j , satisfies 0 ≤ l max ≤ n. Since H 1 is analytic in D ρ,σ (U ) and periodic in ϕ, H 1 admits the Fourier decomposition where, according to Fourier theorem, the coefficients |h 1,k (A)| are bounded by exponentially decaying quantities O(e −|k|σ ). Then, it turns out that the appropriate normal form Z has the resonant form: for some Fourier coefficients ζ k (A ) and where M := {k ∈ Z : k · m l = 0 for all l = 1, . . . , l max } is the 'resonant module'. A normal form of the form (31) implies the existence of n − l max quasi-integrals of the form I i = K i · A, i = 1, . . . , n − l max , where the vectors K i satisfy the equations K i · m l = 0 for all l with 1 ≤ l ≤ l max . The quantities I i are called the resonant integrals of the Hamiltonian (31).

Definition 1. A r-th step Hamiltonian normalization process is a composition of near identity transformations
mapping the initial action-angle variables to the r-th step normalized action-angle variables via the successive transformations ( for all s = 1, . . . , r are analytic and with inverse analytic within non-null domains D ρ (s) ,σ (s) = ∅, and the rth-step Hamiltonian takes the form: with R (r) ρ (r) ,σ (r) Z (r) ρ (r) ,σ (r) . The semi-analytical estimates of stability that we will develop in the next sections are based on defining a suitable r-step sequence of canonical transformations Φ 1 , Φ 2 , . . . , Φ r reducing the size of the remainder R (r) ρ (r) ,σ (r) as much as possible given the initial Hamiltonian model considered. The appropriate sequence is found using the method of Lie series (see Section 4). The obtained times of stability are of order R (ropt) −1 ρ (ropt) ,σ (ropt) , where r opt is the normalization order yielding the smallest possible remainder norm. The value of r opt can be obtained via theoretical estimates (see [11]), but in practice, it is also limited by the maximum order in which our computer-algebra normal form calculations can proceed. Theoretical estimates imply that the size of the remainder norm is exponentially small in the inverse of the size of the perturbation H 1 ρ,σ in Eq. (25). For example, in the simplest case of the Birkhoff normal form, we have the following theorem (see [12] for full details). Theorem 2. Consider the Hamiltonian expressed in action-angle variables H(A, ϕ) = ω · A + f (A, ϕ), where ω ∈ R n satisfies the following Diophantine condition: there exist τ, γ > 0 such that |k · ω| ≥ γ |k| τ ∀k ∈ Z n \{0} (34) and f is real analytic on D ρ,σ for some ρ, σ > 0. Consider two positive parameters δ < ρ/2 and ξ < σ/2, and for r ≥ 1, let * Then, for any there exists a real analytic canonical transformation Φ : D ρ−2δ,σ−2ξ → D ρ,σ such that the transformed Hamiltonian has the form where the remainder R (r+1) can be bounded as Casting together (35) and (38), one readily sees that the remainder grows more rapidly than any power of r, namely as (r τ +2 ) r−1 . Consequently, this procedure does not converge for r → ∞. In any case, we remark that, as the threshold value for the normalization order r is proportional to the inverse of f 1/(τ +2) ρ,σ , if we manage to reduce the size of the initial remainder function, then we can increase the maximum value of r for which Theorem 2 is satisfied.
Similar estimates hold in the case of the resonant normal form constructions (see [11]). The behavior of the size of the remainder as a function of the normalization order r will be examined in detail in our semi-analytical computations in Sections 4 and 5 below.
3.2. Book-keeping and construction of the normal form. Both Hamiltonians (8) and (24) are of the form (25), therefore the above results on Hamiltonian normalization apply. In order to compute the composition of canonical transformations required in Eq. (32), we implement the method of composition of Lie series, after introducing a suitable book-keeping (see [10]) to separate terms in the Hamiltonian according to estimates of their order of smallness.
The proof consists in implementing Proposition 1 of [12] with r = 1. See [14] for the proof. In simple words, the exchange theorem implies that the result of a Lie series canonical transformation onto a function depending on (A, ϕ) can be found by implementing the sequence of Poisson brackets of the exponential operator exp(L χ ) directly on the function f , and substituting, after this operation, the arguments (A, ϕ) with (A , ϕ ).
The above definitions allow us to establish an algorithm for the calculation of the sequence of canonical transformations (32) using Lie series. The algorithm is obtained recursively by defining the r-th step as follows. Assume the Hamiltonian after r − 1 normalization steps, denoted by H (r−1) , is in normal form up to the book-keeping order r − 1: Then, the r-th step Lie generating function χ r and Hamiltonian H (r) are computed as follows: where Z r = Z (r−1) r . Remark 7. In the above algorithm, the notation f (r) implies a function depending on the canonical variables (A (r) , ϕ (r) ), which are connected to the original variables (A, ϕ) via the composition of Lie series transformations For simplicity of notation, unless explicitly required in the sequel we do not write the superscripts in the canonical variables defined in every step, but only in the functions in which these variables are arguments of. r,k (A) k · ω exp(ik · ϕ).

Stability of the semimajor axis in the J 2 model
We will now implement the Hamiltonian normalization discussed in Section 3 to eliminate the short period terms (depending on the mean longitude λ) in the Hamiltonian (8), leading to estimates on the long-term stability of the orbits' semimajor axis.
With reference to the algorithm of Subsection 3.2, normal form terms are specified as those non-depending on the mean longitude λ. After M normalization steps, the Hamiltonian takes the form  (47) we can only compute a truncated remainder, we probe numerically that the finite sum of the leading terms in the remainder (up to order N ) yields a remainder norm close to the limiting one (which corresponds to the limit N → ∞). To this end, we take as maximum normalization order M = N − 3, ensuring that at least the three first leading terms are included in the remainder (see [33]). Also, in estimating the size of the remainder through a suitable definition of the norm, we compute the sup norm on a closed and bounded domain D ⊂ R 2 : In practical computations, we can replace in (48) the supremum with the maximum, since the functions that we consider contain only positive powers of √ 1 − e 2 and cos i, hence, they are continuous on the domain D.

4.2.
Numerical results: stability of the semimajor axis. Having fixed the procedure for the normal form and remainder computations, we proceed in deriving stability estimates based on the time variations of the value of the semimajor axis in the J 2 problem. Fixing a reference value a * of the semimajor axis, we assume that, at the time t = 0, we have L = L * = √ µa * , i.e. δL = 0. Our aim is to estimate the fluctuations of L as functions of the orbital parameters e and i. The first question to settle is that, for every value of the reference parameter a * we have to specify the range of values of the variables (e, i) for which the remainder R (M ) J 2 is small enough to represent only a perturbation with respect to the dynamics determined by the secular part. In applications, we compute the value of R With reference to the Hamiltonian (47), if we consider the dynamics induced only by the secular part, we obtain that which implies that δL (hence L) is a constant of motion. We remind that δL is not the original Delaunay variable, but rather the one obtained after M normalization steps. If we denote by δL (0) the original variable, then we have To obtain δL (0) as a function of the new variable δL, we need to invert the transformation (49); we observe that (exp(L χ )) −1 = exp(−L χ ) , implying Since we are dealing with near-identity canonical transformations, we realize that δL (0) is the sum of δL and short period (small) variations which do not affect its stability. If we consider the full Hamiltonian in (47), then L is not constant anymore because of the dependence of R (M ) J 2 on λ. Using again Hamilton's equations, we see that Then, for every set of values, say (e * , i * , λ * , p * , q * ) ∈ D × T 3 , we obtain d dt L(e * , i * , λ * , p * , q * ) ≤ sup Requiring that the right hand side of (50) is of order of unity, then the stability time T becomes order of O (1/ dL/dt ∞,D ). Let us fix a constant value ∆L and suppose that we want to estimate the minimal time T 1 up to which the variation of L(e, i, λ, p, q; T ) stays bounded by ∆L: L(e, i, λ, p, q; T ) − L * ∞,D ≤ ∆L .
Using (50) we obtain that T 1 is given by Equation (51) can be used to derive the stability time of the semimajor axis a: recalling that, in general, L = √ µa, one has that ∆L = ∆a/2 µ/a, which allows to obtain a lower bound for the stability time of a given by This estimate will be used in Section 5.3 to obtain results on the stability time at different altitudese; in particular, ∆a is set to be equal to 0.1 R E .
To check that the norm R , we compute its value by taking a set of samples for the reference value of the semimajor axis a * , that correspond to different distances from the Earth's center (the radius of the Earth is R E = 6378.14 km). Precisely, we consider the following semimajor axes:           4) * ; the left plots provide the graph of ||dL/dt|| ∞,D as a function of the eccentricity for fixed values of the inclination, while the right plots give the norm as a function of the inclination for fixed values of the eccentricity. We notice that the norms tend to decrease when the eccentricity and the inclination are smaller, although the effect is more evident in the GEO region than closer to the Earth.
We now examine how the stability time changes as a function of the semimajor axis a * : in this case, we consider 1000 values for a * uniformly distributed from a in = 1.15679 (corresponding to an altitude of 1000 km, which we take as the first reference value, although in this region weak dissipative effects are possibily affecting the dynamics) to a f = 16.6786 (corresponding to an altitude of 10 5 km), using Eq. (52) with ∆a = 0.1 R E . Figure 4 confirms that the stability time increases with the altitude also in the case of the semiajor axis. Indeed, while for a * = a in we can ensure the stability of the semimajor axis for a period of the order of years, in the case a * = a f we have a stability time of the order of 10 4 years. From an analytical point of view, this behaviour of the stability time can be explained by the fact that, for higher distances, our model can be approximated by Kepler's problem in which the semimajor axis is constant.

Secular stability in the geolunisolar model
Using the J 2 model, we have demostrated how the stability of the semimajor axis can be established against short-period perturbations (depending on the satellite's mean anomaly). In this section, we focus, instead, on the long-term variations in the eccentricity and inclination of the satellite's orbit, for orbits close to circular (e < 0.1 rad) and with small inclination (|i| < 0.1). One can easily verify that, within the geolunisolar problem (Hamiltonian H ≤N gls,sec , see Eq. (24)), the phase-space manifold e = 0, corresponding to I 2 = 0, constitutes an invariant manifold of the flow, implying that circular orbits remain so for infinitely long times independently of their variations in inclination and longitude of the node. On the other hand, for e small, but non-zero, long-term variations of both the eccentricity and inclination can occur on timescales given by the inverse of the frequencies ν 1 and ν 2 (Eq. (22)). Since ν 1 ν 2 T short , where T short is the characteristic time of the frequency associated to the fast angle. Since J 2 ≈ 10 −3 , the short and long periods are separated by three orders of magnitudes, a fact which justifies altogether the simple averaging over mean anomalies which leads to the model of departure H ≤N gls,sec for the analysis of the secular stability. On the other hand, the fact that ν 1 ν 2 implies that, near the equator (or, more precisely, for orbits near the Laplace plane, see Section 2.2.1), the eccentricity and inclination have coupled variations (the so-called ' Kozai-Lidov' mechanism). This fact implies that, close to the Laplace plane, the term 'secular stability' cannot mean the long-term preservation of the eccentricity and inclination one independently of the other, but only the approximate preservation of the combination I ≈ e 2 + i 2 (see below for exact expressions) known as the Kozai-Lidov integral. The normal form construction and remainder estimates in the present section reflect these basic properties of the dynamics.

Normal form. Starting with the model H ≤N
gls,sec given in Eq. (24), the construction of the normal form proceeds with the algorithm described in Section 3 and the following settings: i) The book-keeping rule (exponent s in Eq. (39)) is set as s = s 1 + s 2 − 2, where s 1 and s 2 are the exponents appearing in Eq. (24).
iii) The maximum truncation order is set to N = 15, while the maximum normalization order is set to M = 12.
With the following settings, the Hamiltonian after r normalization steps, where r can take the values r = 1, 2, ...M , resumes the form: gls,sec (I 1 , I 2 ), hereafter called the secular part, contains all terms independent of the angles (corresponding to the choice k 1 = k 2 = 0 in the resonant module). The dynamics of this term implies separate preservation of the eccentricity and inclination (the latter around the Laplace plane). Instead, Z We now look at the dynamics induced by the sum of secular and resonant parts: gls,res (I 1 , I 2 , φ 1 − φ 2 ), called, altogether, the resonant normal form H norm (for simplicity, we drop the dependence on the normalization order r from the notation). The quantity I 1 + I 2 is a first integral for the dynamics induced by H norm , which implies that the vertical component of the angular momentum, which coincides with Θ, is preserved 2 . Given that L is constant, say L = L * = √ µ E a * , the quantity is a first integral and, as a consequence, the quantity is constant for the dynamics induced by the normal form. This means that e and i can change only in such a way that the value of I(e, i) remains constant. The fact that the presence of resonant first integrals determines a locking in the values of e and i is at the basis of the so-called Lidov-Kozai effect ( [23,21]), which is common, in a wide range of resonant combinations, in many models of Celestial Mechanics.

5.2.
Remainder and stability estimates. As already mentioned in Section 4.2, we need to guarantee that the remainder is small with respect to the normal part; we denote again by D the domain over which the norm R (M ) gls ∞,D is computed where, for a function f = f (e, i, φ 1 , φ 2 ), the norm of f is defined as There exists an optimal value of M that minimizes the estimate of the remainder's norm, as shown in Section 5.3 for GEO orbits.
Since I 1 + I 2 is a first integral for H norm , we have that To evaluate the stability of I(e, i), we use the relation: then, for every (e * , i * , φ * 1 , φ * 2 ) ∈ D × T 2 , we have the following estimate: Let us now consider an orbit with initial point (I 1,0 , I 2,0 ) such that the corresponding eccentricity and inclination belong to D; consider its evolution up to t = T . Using the mean value theorem, we have that Setting Γ to be the maximum value for the variation of I 1 + I 2 in time, let us denote by T the minimum time such that for every T ≤T gls } ∞,D , which gives an estimate for the stability time of I 1 + I 2 and, consequently, of I(e, i). The stability results for the quantity I can be translated in terms of the orbital elements (e, i) as follows: in view of (55), for small values of e and i we find hence, if we consider the variations of I, e and i, they are connected by the relation For the limit case of e or i fixed and small, one find then the relative variation of I (and, as a consequence, of I 1 + I 2 ) is proportional to the relative variations of the orbital elements by a factor 2. To make the stability results for the geolunisolar model consistent with the ones obtained in Section 4.2 for the J 2 model, in Section 5.3 we set namely, recalling that ∆L = ∆a/2 µ/a and ∆a = 0.1, the maximal variation of I 1 + I 2 in the geolunisolar model is equal to the maximal variation allowed for the action L in the J 2 model. Since the stability results strongly depend on the distance from the Earth, we select five different altitudes, that correspond to cases of interest for the satellite's problem: • h (1) = 3000 km, above the atmosphere; • h (2) = 20000 km, that is in MEO region; • h (3) = 35786 km, the altitude of GEO orbits; • h (4) = 50000 km, corresponding to far objects; • h (5) = 100000 km, that corresponds to objects which are very far from the Earth's surface. The value of the remainder's norm depends on the altitude of the orbit: in particular, we can state that the stability time decreases as the altitude increases. Table 3 provides the value of R   As mentioned in Section 5, the remainder's norm depends on the normalization order M . Although the norm does not converge to zero if M tends to infinity, there is a value of M , called the optimal normalization order, say M opt , for which the norm of the remainder is minimal. Typically, this optimal value is greater than the order of the Taylor expansions of the numerically computed functions, and the estimates for the remainder is so good that there is no reason to push further the order of the expansion; for example, this is the case for the normalized Hamiltonian function which describes the geolunisolar problem computed for the GEO altitude. As we can see from Figure 6, the optimal normalization order is greater than or equal to 11. Since the values of R gls ∞,D in D, we proceed to compute the stability time for the quantity I(e, i) = 1 − √ 1 − e 2 (1 − cos i). As we can see from Table 5, the stability times are extremely long: this fact depends on the model we considered, with the Lunar orbit in the ecliptic plane without precession effects. However, we can notice a relevant decrease in the stability time for distances greater than GEO. This behaviour is opposite to that of the J 2 model where the stability time was increasing with the altitude (see Figure 4). In fact, at low altitudes the J 2 model is strongly affected by the Keplerian part and the geopotential, while the geolunisolar model takes into account both the inner effect due to the Earth and the outer effect due to Moon and Sun. As a final remark, to show the importance of taking the right domain in eccentricity and inclination, let us assume h = h (5) Table 4). If we compute the stability time in the enlarged domain B, we obtain just the value T = 0.00164 years.

Non-degeneracy conditions
Beside the orbital stability obtained using the Birkhoff normal form as in the previous Sections, analytical estimates of the stability time can be obtained by the outstanding theorem developed by Nekhoroshev ([26]). Under suitable assumptions, the theorem gives a confinement of the action variables for exponentially long times. In particular, the Hamiltonian must satisfy a non-degeneracy condition which, in the original formulation, is called steepness condition. The definition of the steepness condition is quite technical and typically not trivial to verify for a specific Hamiltonian system. However, there are some sufficient conditions which imply steepness, whose verification requires the resolution of algebraic equalities and inequalities. This motivates the introduction of the following definition (see [20,28]).
(3) h(J) is three-jet non degenerate in J ∈ B if ω(J) = 0 and We remark that the convexity condition is equivalent to require that the Hessian matrix Q(J) is positive (or negative) definite in J. We add also the following definition of isoenergetically non-degenerate which, for Hamiltonian systems with 2 degrees of freedom, implies quasi-convexity.
One can prove (see [1]) that, for every Hamiltonian system with n degrees of freedom, quasi-convexity implies isoenergetically non-degeneracy: as a consequence, for twodimensional Hamiltonian systems, the two conditions are equivalent.
6.1. Numerical verification of the non-degeneracy conditions. We now apply the above definitions to the Hamiltonian functions introduced in Section 2. We consider the following cases:  H gls (I 1 , I 2 , φ 1 , φ 2 ) = H gls,sec (I 1 , I 2 ) + H gls,res (I 1 , I 2 , φ 1 , φ 2 ) + R gls (I 1 , where H gls,res depends only on the quasi-resonant combination φ 1 − φ 2 . To analyze the non-degeneracy conditions, we write the Hamiltonian as the sum of two terms, namely an integrable Hamiltonian h and a perturbative function f . For the J 2 -Hamiltonian, we set h(P, Q) to contain all the terms of H J 2 that are independent on all angles, while the perturbing function f contains all other terms. For the geolunisolar case, we choose h(I 1 , I 2 ) to be the angle-independent part of the truncation up to order 2 of H gls : in this way, the Hessian matrix of h is independent of the actions, and the computations are easier 3 .
Since the Hamiltonian functions depend on the parameter L * = √ µa * , we select four reference values for the altitudes that correspond to distances of interest in satellite dynamics: • 3000 km, for near-Earth objects; • 20000 km, for distance of the order of MEO; • 35790 km, for GEO orbits; • 50000 km, for far objects. is non zero.
If the convexity and quasi-convexity tests fail, one can control the three-jet nondegeneracy condition, that we compute, again, numerically, checking that the system      ω(P, Q) · u = 0 (∂ 2 h(P, Q)u) · u = 0 ((∂ 3 h(P, Q)u)u) · u = 0 (62) evaluated on a grid of values (P, Q) ∈ D admits only the trivial solution u = (0, 0, 0). Since convexity implies quasi-convexity and quasi-convexity implies three-jet non-degeneracy, to identify which of the conditions is satisfied, we proceed in the following way:  • we begin with the convexity test on the product of the eigenvalues: if the product is positive for every value of (P, Q) ∈ D , then h(P, Q) is convex; • if the convexity test fails, we pass to the quasi-convexity condition, checking the criteria given in Definition 10 and Remark 12; • if the quasi-convexity test fails, we check the three-jet non-degeneracy through the numerical test based on Definition 10.
6.2. Non-degeneracy of the J 2 Hamiltonian. We start from the convexity test; we denote by λ 1 , λ 2 the eigenvalues of the Hessian matrix of h. Table 6 gives the numerical values of λ 1 λ 2 for different altitudes and (P, Q) in the domain D (we recall that, since the values in the Hessian matrix depend on P and Q, we have an interval for λ 1 λ 2 instead of a single value). As one can see, the product of the eigenvalues is always negative or zero within numerical precision level, leading to the conclusion that the Hamiltonian H J 2 is not convex in D for the considered altitudes.
We can then pass to the quasi-convexity test. We consider the determinant of the matrix A defined in (61) for (P, Q) ∈ D . As we can see from Table 7, for every considered altitude the values of det A are equal to zero within the numerical precision level, leading to the conclusion that the J 2 Hamiltonian is not quasi-convex in D .
The failure of the quasi-convexity for the J 2 problem is a relevant fact: as we will see in Section 6.3, the effects of the lunisolar attraction will eliminate such degeneracy, making the total Hamiltonian quasi-convex.
We conclude with the test on the three-jet non-degeneracy condition. To make the computations quantitative, we solved the system (62) for values (P i , Q j ) on a mesh of 10000 points in D . For every pair of values (P i , Q j ) the only solution of the system is the trivial one u = (0, 0), leading to conclude that the Hamiltonian of the J 2 model is three-jet non degenerate in D .
6.3. Quasi-convexity of the geolunisolar Hamiltonian. As for the J 2 model, we start from the convexity test. In this case, the unperturbed Hamiltonian is a polynomial of degree 2 in the actions; then, the Hessian matrix of h(I 1 , I 2 ) does not depend on the values of I 1 and I 2 , and the same holds for its eigenvalues. This makes the test on the convexity of the Hamiltonian easier. Table 8 shows the values of λ 1 and λ 2 for different altitudes. As we can see, in every case the eigenvalues of the Hessian have opposite sign, showing that the geolunisolar unperturbed Hamiltonian is not convex in R 2 , and hence in D .
As for the quasi-convexity, we check whether the matrix A defined in (61) is nondegenerate for every value (I 1 , I 2 ) ∈ D .
From Table 9 we can see that the determinant of A is strictly positive for every value of the selected altitudes and every (I 1 , I 2 ) ∈ D . Hence, we conclude that the Hamiltonian for the geolunisolar case is quasi-convex. As observed at the end of Section 6.2, this fact is highly nontrivial, since it means that the lunisolar perturbation of the J 2 model removes the degeneracy.