Satellites’ orbital stability through normal forms

Abstract A powerful tool to investigate the stability of the orbits of natural and artificial bodies is represented by perturbation theory, which allows one to provide normal form estimates for nearly-integrable problems in Celestial Mechanics. In particular, we consider the orbital stability of point-mass satellites moving around the Earth. On the basis of the J 2 model, we investigate the stability of the semimajor axis. Using a secular Hamiltonian model including also lunisolar perturbations, the so-called geolunisolar model, we study the stability of the other orbital elements, namely the eccentricity and the inclination. We finally discuss the applicability of Nekhoroshev’s theorem on the exponential stability of the action variables. To this end, we investigate the non-degeneracy properties of the J 2 and geolunisolar models. We obtain that the J 2 model satisfies a “three-jet” non-degeneracy condition, while the geolunisolar model is quasi-convex non-degenerate.


Introduction
The study of the stability of celestial bodies is one of the major goals in Celestial Mechanics, especially in view of the growing importance taken, since the mid-20th century and in particular in the last decades, by the problem of the control of satellites and space debris orbiting around the Earth.
The dynamics of a small body moving around our planet is influenced by a number of factors, among which the most important are the attraction of the Earth, taking into account its non-spherical shape, and of the Sun and the Moon, which, for altitudes ranging from about 100 to 10 5 km, can be considered as third body perturbations to the dominant effect determined by the geopotential.
This work summarizes some stability results presented in De Blasi et al. 2021, where, to describe the small body's motion, two different models are considered: the J 2 model, which takes into account only the gravitational force of the Earth by including the associated geopotential, expanded in spherical harmonics and suitably truncated up to its dominant term, and the geolunisolar model, where the effects of Sun and Moon are added to the geopotential (see §2). The Hamiltonian formalism, as well as a particular set of actionangle variables (the so-called modified Delaunay variables), are used to describe both models.
Through a suitable sequence of canonical changes of coordinates (see §3), the Hamiltonians describing the J 2 and geolunisolar models, denoted respectively with H J2 and H gls , are transformed to assume the form of quasi-integrable Hamiltonian functions, where one or more actions can be considered as quasi-integrals of motion (namely integrals for a truncated normal form Hamiltonian). Their stability can be estimated by means of suitable perturbation theory procedures as in §4 (compare with Steichen & Giorgilli 1997). We also mention that other tools can be used to obtain stability estimates for quasi-integrable systems, such as Nekhoroshev's theorem (see Nekhoroshev 1962), which, under suitable hypotheses, provides exponentially long stability times. One of the main hypotheses required is a non-degeneracy of the integrable part of the Hamiltonian. In §5 we discuss the non-degeneracy condition of the normalized Hamiltonians describing the J 2 and geolunisolar models.

The models
Let us consider a point-mass particle, say a debris S, moving around the Earth along an elliptic orbit with parameters (a, e, i, M, ω, Ω), that is, semimajor axis, eccentricity, inclination, mean anomaly, argument of perigee and longitude of the ascending node. The motion of S is governed by the Earth's gravitational influence, whose associated potential can be expressed as an expansion in spherical harmonics (see Kaula 1966) which takes into account the non-spherical shape of our planet. Additionally to Earth's Keplerian attraction, the gravitational influences of Sun and Moon are taken into account.
Both models introduced in §1 are analysed by considering the associated Hamiltonian functions, denoted respectively with H J2 and H gls , first expressed in Cartesian coordinates and then in the modified action-angle Delaunay variables where μ E = GM E is the Earth's mass parameter.
The Hamiltonian H J2 can be expressed as the sum of a zero-order Keplerian term and the J 2 term, which indicates the deformation in the geopotential due to Earth's oblateness. In Delaunay variables, having defined δL = L − L * = √ μ E a − √ μ E a * with a * taken as a reference value for the semimajor axis, H J2 can be expanded in powers of √ δL, √ P and √ Q obtaining the Hamiltonian (δL, P, Q) cos (k 1 λ + k 2 p + k 3 q) , are polynomials of degree s in the actions and the sum is truncated up to order 2N † for computational reasons.
The geolunisolar Hamiltonian H gls is obtained by adding to H J2 the third-body gravitational potentials due the presence of Sun and Moon, assuming the latter to be strictly on the ecliptic plane. Since in this case the attention is focused on the secular stability of the parameters (e, i), an average over the satellite's, Sun's and Moon's fast angles is performed, implying the constancy of the semimajor axis a = a * . Moreover, the presence of two perturbing bodies on the ecliptic produces a shift in the equilibrium orbit of S (leading to the so-called forced elements), which modifies (e, i) = (0, 0) into (e, i) = (0, i * ), with i * = 0: a new set of action-angle variables (I 1 , I 2 , φ 1 , φ 2 ), centered in the corresponding forced values, is considered. Similarly to H J2 , the final geolunisolar Hamiltonian can be expressed as a trigonometric expansion in the new variables: where we notice that ν 1 ν 2 .

Hamiltonian normalization
Taking into account the models described in §2, we can obtain estimates on the stability times of the orbital parameters of the satellite, in particular of the semimajor axis in the J 2 model and of the quantity √ 1 − e 2 (1 − cos i) in the geolunisolar case. The approach taken to achieve this goal is heavily based on normal form methods (see Efthymiopoulos 2011, to which we also refer for the definition of Lie series transformations): in particular, using the Lie series technique, a formal elimination of the fast angle in the J 2 model and of a particular quasi-resonant combination of the angles within the geolunisolar framework is performed, reducing to new quasi-integrable Hamiltonian functions for which the stability estimates are obtained.
In the case of the model described by H J2 , expressed by the expansion 2.2, a composition of near-identity canonical transformations in the form of Lie series allows to remove the dependence on the fast angle λ up to a prefixed polynomial order in the actions' square roots, leading to the normalized Hamiltonian † where Z J2 , the so-called normal part, does not depend on λ and R J2 , the remainder, is of total order M = N − 3 in √ δL, √ P , √ Q. It is straightforward from the independence of Z J2 on λ that the variation of the first Delaunay action L (and, as a consequence, of the satellite's semimajor axis a) depends only on the remainder, as the former is a first integral for the normal part: in this sense, it can be considered a quasi-integral of the motion induced by the whole Hamiltonian function H norm J2 .
In the case of H gls in the form of 2.3, the presence of the 1 : 1 resonance determined by the relation ν 1 ν 2 translates in the unfeasibility of a normalization procedure which simply removes the dependence on the angles of selected terms, as it would imply the uncontrolled growth in the size of the remainder due to the presence of small divisors. We opt instead for a normalized Hamiltonian in quasi-resonant form, in the sense that the composition of Lie series transformations removes the dependence on all combinations of the angles φ 1 , φ 2 except for the resonant one φ 1 − φ 2 . The resulting normalized geolunisolar Hamiltonian takes then the form H norm gls (I 1 , I 2 , φ 1 , φ 2 ) = Z sec gls (I 1 , I 2 ) + Z res gls (I 1 , where the sum Z sec gls + Z res gls represents the normal part, divided into its secular and resonant terms, and the remainder R gls is again of total order M = N − 3 in √ I 1 , √ I 2 . The † From a rigorous point of view, the variables of H norm J 2 are not the original Delaunay variables (δL, P, Q, λ, p, q). However, since we are dealing with near-identity transformations, they differ from them only by short-term small oscillations which do not affect the secular stability of the orbital elements: for this reason, with an abuse of notation, we keep the original notation for the new Delaunay variables as well. dynamics induced by the normal part alone admits the integral which, similarly to the J 2 case, can be considered a quasi-integral for the whole Hamiltonian H norm gls and whose stability will be investigated in §4.

Stability estimates
Taking advantage of the peculiar structure of the normalized functions H norm J2 and H norm gls , stability estimates based on the size of R J2 and R gls can be produced. The size of the remainders can be quantified by means of the sup norm · ∞,D over a suitable bounded domain D in the variables, which depends on the model. In particular, if the remainders are small in the above sense with respect to the corresponding normal parts, the dynamics induced by H norm J2 (respectively H norm gls ) can be considered as a small perturbation of that defined by Z J2 (respectively, Z sec gls + Z res gls ). In the case of the J 2 model, the derivative of the first Delaunay action L depends only on the λ−derivative of R J2 : (4.1) Let us consider a bounded domain D ⊂ [0, 1) × [0, π/2] in eccentricity and inclination; recalling the dependence of the Delaunay actions on the orbital parameters, define the sup-norm i, λ, p, q) . (4.2) Suppose now that at time t = 0 the initial value of L is L 0 ; defining by L(T ) the evolution of L at time T , from (4.1) and the mean value theorem, one obtains Fixing a constant value ΔL, one has that a lower bound for the time T such that |L(t) − L 0 | ≤ ΔL is given by With the same reasoning, recalling the definition of Poisson brackets and their relation to the time evolution of functions under the dynamics induced by a Hamiltonian (see for example Giorgilli 2002), an analogous stability estimate can be obtained for the geolunisolar model: in particular, if Γ is the maximal variation allowed for the quantity where · ∞,D is the supremum of the absolute value of the argument for (e, i) ∈ D and (φ 1 , φ 2 ) ∈ T 2 . The theoretical considerations leading to the definition of the stability times T J2 and T gls are the basis of the numerical investigations on the stability for the J 2 and geolunisolar models for different reference values a * of the semimajor axis. Such estimates are obtained by implementing the following procedure:  • compute the initial Hamiltonians H N J2 and H N gls , obtained as finite trigonometric expansions in the action-angle variables; • normalize up to order M = N − 3 to obtain H norm J2 and H norm gls ; • compute the stability times T J2 and T gls on suitable domains and for chosen values of ΔL and Γ, provided that the remainders are small enough to represent small perturbations of the respective normal parts. Given a reference value a * for the semimajor axis, the maximal excursions for L(T ) and (I 0 + I 1 )(T ) are chosen to be where the first choice corresponds to a maximal variation of Δa = 0.1R E for the semimajor axis. Furthermore, for computational reasons, the supremum norm defined in 4.2 is replaced by an appropriate norm · ∞,D * over a finite grid D * ⊂ D, which turns out to be an upper bound of · ∞,D . To ensure the smallness of the remainders over the normal parts, D is set equal to [0, 0.15] × [0, π/2] for the J 2 model and to [0, 0.1] × [0, 0.1] in the geolunisolar case. Figure 1 shows the behaviour of T J2 for altitudes ranging from 10 3 to 10 5 km: as expected by the fact that the Earth's oblateness is more relevant for small altitudes, the stability time increases with a * , when the J 2 model approaches the classical Keplerian one. Table 1 lists the values of T gls for selected values from small to high altitudes; since the lunisolar influence is stronger far from the Earth, the stability times tend to decrease with the altitude.

Non-degeneracy conditions
One of the most important results of the XX century on the stability of nearlyintegrable Hamiltonian systems is represented by Nekhoroshev's theorem (Nekhoroshev 1962;Poshel 1993), which is a fundamental tool to provide exponential stability estimates for quasi-integrable systems of the form where (I, φ) ∈ U × T n , U ⊂ R n , n being the number of degrees of freedom of the system and ε a small parameter. Typically, h is called the unperturbed function, and is straightforwardly integrable as it does not depend on the angles φ, while f is the perturbing function.
A fundamental assumption in the original version of Nekhoroshev's theorem is a non-degeneracy of the unperturbed function h(I) called steepness condition, which is a geometric assumption rather complex to verify in practice. However, the steepness condition is implied by stronger non-degeneracy conditions, much simpler to check, such as the convexity, the quasi-convexity and the three-jet non-degeneracy (see, e.g., Chierchia, et al. 2018).
As a preliminary investigation of the applicability of Nekhoroshev's theorem to our cases, we establish whether these conditions are verified or not in the J 2 and geolunisolar models: to this end, a suitable splitting of H norm J2 and H norm gls into unperturbed and perturbed parts is needed. In particular, as the non-degeneracy conditions involve only h(I), one has to define precisely the unperturbed functions, denoted respectively as h J2 and h gls . In the case of the J 2 model, we impose h J2 (P, Q) to be the sum of all the angle-independent terms of H norm J2 , where, given the practical stability of L shown in §4, δL is supposed to be constant and equal to 0. In the geolunisolar model, h gls (P, Q) is given by the sum of all the terms of H norm gls which are independent on the angles and at most quadratic in the actions.
The • convexity: check if the Hessian matrices associated to h J2 and h gls are positively (or negatively) defined; • quasi-convexity: it is equivalent to Arnold isoenergetic non-degeneracy condition, involving first and second order derivatives of the unperturbed Hamiltonian; • three-jet non-degeneracy: it is a condition involving up to third order derivatives and it was verified numerically on a grid of 10000 points in the actions. For the considered altitudes and in the regime of eccentricities and inclinations defined by D, we find that h J2 is three-jet non-degenerate, while h gls is quasi-convex; this result is of particular relevance, as it implies that the presence of the lunisolar part removes the degeneracy of the J 2 model.