Long-Term Evolution of Highly-Elliptical Orbits: Luni-Solar Perturbation Effects for Stability and Re-entry

This paper investigates the long-term evolution of spacecraft in Highly Elliptical Orbits (HEOs). The single averaged disturbing potential due to luni-solar perturbations, zonal harmonics of the Earth gravity field is written in mean Keplerian elements. The double averaged potential is also derived in the Earth-centered equatorial system. Maps of long-term orbit evolution are constructed by measuring the maximum variation of the orbit eccentricity to identify conditions for quasi-frozen, long-lived libration orbits or initial orbit conditions that naturally evolve toward re-entry in the Earth’s atmosphere. The behavior of these long-term orbit maps is studied for increasing values of the initial orbit inclination and anomaly of the perigee with respect to the Moon’s orbital plane. In addition, to allow meeting specific mission constraints, quasi-frozen orbits can be selected as graveyard orbits for the end-of-life of HEO missions. On the opposite side, unstable conditions can be exploited to target Earth re-entry at the end-of-mission.


INTRODUCTION
Highly Elliptical Orbits (HEOs) about the Earth are often selected for astrophysics and astronomy missions, as well as for Earth missions, such as Molniya or Tundra orbits, as they offer vantage point for the observation of the Earth and the Universe (Draim et al., 2002). Moreover, elliptical Geostationary Transfer Orbits are commonly selected to inject spacecraft in geostationary orbit. HEOs guarantees spending most of the time at an altitude outside the Earth's radiation belt; therefore, long periods of uninterrupted scientific observation are possible. In addition, geosynchronicity is opted to meet coverage requirements, enhanced at the apogee, and optimize the ground station down-link. If the inclination is properly selected, HEO can minimize the duration of the motion inside the eclipses. This paper, whose preliminary version was presented at the 25th AAS/AIAA Space Flight Mechanics Meeting, in Williamsburg (VA) in January 2015 (Colombo, 2015), investigates the longterm evolution of spacecraft in HEOs through the exploitation and development of semi-analytical techniques. The dynamics of HEOs with high apogee altitude is mainly influenced by the effect of third body perturbations due to the Moon and the Sun, which induces long-term variations in the eccentricity and the inclination, corresponding to large fluctuations of the orbit perigee and the effect of the Earth's oblateness.
The semi-analytical technique based on averaging is an elegant approach to analyze the effect of orbit perturbations. It separates the constant, short periodic and long-periodic terms of the disturbing function. The short-term effect of perturbations is eliminated by averaging the variational equations, or the corresponding potential, over one orbit revolution of the small body. Indeed, averaging corresponds to filtering the higher frequencies of the motion (periodic over one orbit revolution), which typically have small amplitudes (Ely, 2014). The resulting system allows a deeper understanding of the dynamics (Shapiro, 1995;Krivov and Getino, 1997). Moreover, the use of the average dynamics reduces the computational time for numerical integration as the stiffness of the problem is reduced, while maintaining sufficient accuracy compatible with problem requirements also for long-term integrations.
The effect of third body is usually modeled as a series expansion of the potential with respect to the ratio between the orbit semi-major axis and the distance to the third body. In averaged development the potential is usually truncated to the second order. For example, Cook's formulation gives the secular and long-periodic perturbation due to luni-solar perturbation obtained through averaging over one orbit revolution of the satellite (Cook, 1962). It assumes circular orbit for the disturbing bodies and considers only the second term of a/a ′ , where a and a ′ are, respectively, the spacecraft and the disturbing body semi-major axis about the Earth (Blitzer, 1970). However, it does consider the obliquity of the Sun and the Moon over the equator and the precession of the Moon plane due to the Earth's oblateness (in a period of 18.6 years with respect to the ecliptic). For orbits above Low Earth Orbits, considering only the second order is not enough to accurately estimate the effect of lunisolar perturbations. For example, Lara et al. (2012) focused on orbital configurations, such as the Global Navigation Satellite Systems, for which the second order effects of J 2 can be of the same order of magnitude as perturbations due to the P 2 up to P 5 terms in the Legendre polynomials expansion of the third-body's disturbing function. For this reason, some recursive formulations of the third body potential were developed (Cefola and Broucke, 1975; Laskar and Bou, 2010) also for recovering the short periodic effects (El'yasberg, 1967).
In this paper, the single averaged disturbing potential due to luni-solar perturbations is developed in series of Taylor of the ratio between the orbit semi-major axis and the distance to the third body, following the approach by Kaufman and Dasenbrock (1972). The effect of other zonal harmonic of the Earth gravity field is also modeled up to order 6 considering also the J 2 2 term (Liu and Alford, 1980). As we want to focus here on the interaction between the terms of the Legendre polynomial of the luni-solar perturbation and the Earth's oblateness, the effect of solar radiation pressure and aerodynamics drag is neglected in this work. The perturbed dynamics is propagated in the Earth-centered equatorial frame by means of the single averaged variation (over the orbit revolution of the spacecraft) of the disturbing potential, implemented in the suite PlanODyn (Colombo, 2016). The long-term evolution of high elliptical orbits is characterized in the phase space of eccentricity, inclination and argument of the perigee with respect to the Earth-Moon plane. Maps of long-term evolution are constructed to assess the maximum and minimum eccentricity attained as function of the initial conditions. Through these maps, conditions for quasi-frozen, or long-lived libration orbits are identified. Recent projects funded by the European Space Agency on the design of disposal trajectories for Medium Earth Orbits (Rossi et al., 2015), Highly Elliptical Orbits and Libration Earth Orbits (Armellin et al., 2014;Colombo et al., 2014b) have demonstrated the possibility of exploiting orbit perturbations for designing of passive mitigation strategies for debris disposal. From the results of these projects, many work stemmed out proposing the use of maps for characterizing conditions for natural re-entry in low Earth orbit (Alessi et al., 2018), medium Earth orbit Daquin et al., 2016;Gkolias et al., 2016;Skoulidou et al., 2017), geostationary Transfer orbit (Srongprapa, 2015), geostationary orbit Gkolias and Colombo, 2017), and highly elliptical orbit (Colombo, 2015).
Third body perturbations show important effects on the long-term evolution of orbits within the Medium Earth orbit region and above. A relevant contribution to the understanding of high-altitude orbit dynamics was given by Kozai's (1962) analytical theory on secular perturbations of asteroid with high inclination and eccentricity. Assuming one perturbing body on a circular orbit, and considering the second term of the disturbing potential, an analytical solution was found, later named as Kozai-Lidov dynamics (Kozai, 1962;Lidov, 1962), which links the eccentricity of the orbit with the inclination, measured from the perturbing body plane, and the argument of the perigee. Their work was more recently extended to consider the slow variation due to the 4 th term of the potential and highlighting shirking characteristics in the phase space (Katz et al., 2011;Naoz et al., 2013). Always in the line of long term analysis of the third body effect, El'yasberg (1967) derived the double averaged equations of the second term of the disturbing potential and Costa and Prado (2000) continued on the effort by El'yasberg by expanding the derivation of the double averaged potential up to order 8th. Their interest concerned the critical value of the inclination between the perturbed and the perturbing body related to the stability of near-circular orbits. In other words for inclination higher than the critical values, circular orbits get very elliptic, while for lower values the orbit stays nearly circular. In El'yasberg (1967) and Costa and Prado (2000) the double averaging was performed using the orbital elements of the spacecraft as defined with respect to the third body plane. However, when considering both the third body effects of the Sun and the Moon, one has to assume that they orbit on the same plane; this reduce to consider that the 5.1 • inclination of the Moon's plane over the ecliptic is equal to zero.
In this work, to guarantee the consistency with the singleaveraged approach, we do not make this assumption. The double average variation is here obtained averaging on the fast variable describing the orbital motion of the perturbing body around the Earth as in El'yasberg (1967) and Costa and Prado (2000), but the different inclination of the perturbing bodies planes is retained. The double-averaged disturbing potential is derived in the Earth-centered equatorial reference system. The choice of this system has the advance of allowing the description of the perturbing effect of the Sun and the Moon, considering the actual ephemerides and their inclination with respect to the equator, and accommodating also the inclusion of the effect of the zonal harmonics of the Earth potential, which also affect the motion.
By using the single averaged and double averaged equations presented in this paper, the behavior of quasi-frozen solutions appearing for high inclinations orbits can be reproduced. The choice of representing the maps in terms of inclination and argument of the perigee with respect to the Moon plane is quite appropriate here as we want to study the effects on highly elliptical orbit for which the Moon third body perturbations is the most relevant effect. In addition, to allow meeting specific mission constraints, stable conditions for quasi-frozen orbits can be selected as graveyard orbits for the end-of-life of HEO missions, in the case re-entry option cannot be achieved due to propellant constraints, such as XMM-Newton, which is taken here as practical example. On the opposite side, unstable conditions can be exploited to target an Earth re-entry at the end-of-mission (Jenkin and McVey, 2008;Colombo et al., 2014a). This is the case of the end-of-life of INTEGRAL mission, requiring a small delta-v maneuvers for achieving a natural reentry assisted by perturbations. In this paper, maps of stable and unstable HEOs are built, to be used as preliminary design tool for graveyard or frozen orbit design or natural re-entry trajectories at the end-of-life. Given the available delta-v on-board, the reachable space of orbital elements can also be identified as in the case of XMM-Newton mission.

Orbit Evolution With Single-Averaged Disturbing Potential
To analyze the long-term and secular effect of orbit perturbations, it is convenient to use an averaging approach. In the case the effect of perturbations is conservative, this can be described through a disturbing potential R. In this work we consider the perturbation to the two-body dynamics due to the zonal harmonics of the Earth gravity field, and the third-body perturbation of the Sun and the Moon as The variation of the orbital elements is described through the planetary equations in the Lagrange form (Battin, 1999): that can be written in condense form as where α is here used as the condensed form of the Keplerian elements α = a e i ω M T , where a is the semi-major axis, e the eccentricity, i the inclination, the right ascension of the ascending node, ω the argument of the perigee and M the mean anomaly. Through single averaging operation, the potential can be replaced by the orbit-averaged form of the disturbing function: obtained under the assumption that the orbital elements are constant over one orbit revolution of the spacecraft around the central planet. Therefore, the variation of the mean elements is described by: where nowᾱ is the vector of the averaged orbital elements.

Luni-Solar Averaged Potential
For describing the effects of luni-solar third body perturbations, we follow the approach proposed by Kaufman and Dasenbrock (1972). Their approach is summarized in this section as is fundamental to give an insight into the orbital dynamics that is exploited in section Luni-Solar and Zonal Harmonics Maps. The disturbing potential due to the third body perturbation is (Murray and Dermott, 1999): where µ ′ is the gravitational coefficient of the third body, r and r ′ are the position vectors of the satellite and the third body with respect to the central planet, respectively, as represented in Figure 1. Equation (5) can be expressed as function of the angle ψ between r and r ′ if the cosine rule is exploited in the first term of Equation (5) and the dot product in the second term is resolved (Murray and Dermott, 1999): where cos ψ = r · r ′ rr ′ Kaufman and Dasenbrock (1972) express the disturbing potential as function of the spacecraft's orbital elements, choosing as angular variable the eccentric anomaly E, the ratio between the orbit semi-major axis and the distance to the third body r ′ : δ = a r ′ and the orientation of the orbit eccentricity vector with respect to the third body (Kaufman and Dasenbrock, 1972): where the eccentricity unit vectorP, the semilatus rectum unit vectorQ, and the unit vector to the third bodyr ′ are expressed with respect to the equatorial inertial system, through the following composition of rotations (Kaufman and Dasenbrock, 1972;Lara et al., 2012): where R 1 represents the rotation matrix around the x axis, R 2 the rotation matrix around y axis and R 3 the rotation matrix around the z axis. The full expression ofP,Q, andr ′ in terms of Keplerian elements can be found in Chao-Chun (2005). The variables ′ , ω ′ , i ′ , and f ′ in Equation (7) are, respectively, the right ascension of the ascending node, the argument of the perigee, the inclination and the true anomaly of the perturbing body on its orbit (described with respect to the Earth-centered equatorial reference frame) and u ′ = ω ′ + f ′ . Under the assumption that the parameter δ is small (i.e., the spacecraft is far enough from the perturbing body), Equation (6) can be rewritten as a Taylor series in δ as Kaufman and Dasenbrock (1972): Note that the summation starts from k = 2 as the term zero of the series is influent as is a constant, while the term 1 simplifies with the second term of Equation (6). The average operation in eccentric anomaly, passing through the mean anomaly via dM = (1 − e cos E) dE can then be performed, assuming that the orbital elements of the spacecraft a, e, i, , and ω are constant over one orbit revolution: where the averaged termsF k (A, B, e) are reported in Kaufman and Dasenbrock (1972) and in Appendix in a more compact form. Note that the term 2 of Equation (8) is the one given by Chao-Chun (2005). Equation (8) can be now inserted into the Lagrange equations by computing the partial derivatives with respect to the orbital elements, considering that the dependences of the terms in Equation (8) The derivatives up to the 8 th order of the Taylor series are reported by Kaufman and Dasenbrock (1972); we report them in Appendix in a more concise form up to order 6 th .

Earth Zonal Harmonic Potential
The disturbing function of the zonal harmonic potential is expressed also in terms of classical orbital elements. The zonal harmonics were modeled up to order 6 considering also the J 2 2 term (Blitzer, 1970;Liu and Alford, 1980). We report here only the term associated with J 2 that is the one most important one for the application considered in this work. However, all the terms were retained in the orbit propagation.
where W is the oblateness parameter where J 2 = 1.083 · 10 −3 denotes the second zonal harmonic coefficient and R Earth is the mean radius of the Earth. n = µ Earth /a 3 is the orbit angular velocity of the spacecraft on its orbit, with µ Earth the gravitational constant of the Earth.

Model Validation
The averaged dynamical model described in section Orbit Evolution With Single-Averaged Disturbing Potential was validated by comparison with the actual ephemerides of two artificial satellites in highly-elliptical orbit: INTEGRAL and XMM-Newton. The orbit of INTEGRAL was used in Colombo et al. (2014a) and validated against the ephemerides from the NASA-Horizon system. XMM-Newton orbit was propagated in the time span from 1999/12/15 to 2013/01/01 with the initial Keplerian elements on 1999/12/15 at 15:00 as: a = 67045 km, e = 0.7951, i = 0.67988 rad, Ω = 4.1192 rad, ω = 0.99259 rad, true anomaly f = 3.2299 rad, and radius of the perigee 13,737 km. Figure 2 shows the results of the validation with the actual ephemerides form ESA (2013) (blue line). Luni-solar and Earth zonal harmonics perturbations are included in this validation with PlanODyn (red line) (Colombo, 2016). The reference system used is inertial, centered at the Earth, on the Earth equator.

Orbit Evolution With Double-Averaged Disturbing Potential
Under the further assumption that the orbital elements do not change significantly during a full revolution of the perturbing body around the central body (i.e., Earth), the variation of the orbit over time can be approximately described through the disturbing potential double averaged over one orbit evolution of the s/c and over one orbital revolution of the perturbing body (either the Moon or the Sun) around the Earth. The terms of the disturbing potential due to the third body effect in Equation (3) can be substituted by the double-averaged onē R 3−Sun andR 3−Moon : In this work we decided to express the double-averaged potential with respect to the Keplerian elements described in the Earth's centered equatorial reference system. This will give a more complex expression for the potential in Equation (11), with respect to the one by El'yasberg (1967) and Costa and Prado (2000) but it has the advantage of avoiding the simplification that Moon and Sun orbit on the same plane and facilitating the introduction of the effect of the zonal harmonics. The terms of the third body potential are obtained starting from the single-averaged terms, performing the averaging operationF k (A, B, e). where = − ′ . The term up to order four of the potential were computed; only the second termF 2 is reported here asF 3 = 0: · cos 2i ′ 6 + 9e 2 + 90e 2 cos (2ω) sin 2 i + 2 + 3e 2 +15e 2 cos (2ω) 1 + 6 cos (2 ) sin 2 i ′ +3 cos (2i) 6 + 9e 2 cos 2i ′ + −2 − 3e 2 + 5e 2 cos (2ω) The evolution of the orbit in time can be then computed computing the partial derivatives of Equation (10)    Note that, with the same procedure, the double averaged disturbing potential can be also written in the form proposed by El'yasberg (1967) and Costa and Prado (2000). A different reference system needs to be used, which is still centered at the central body (i.e., Earth) but the x-y plane lay on the perturbing body orbital plane, with its x-axis in the direction of the perturbing body on its orbit, the z-axis in the direction of the perturbing body angular momentum, and the y-axis that completes the reference system. The corresponding orbital elements are α 3Bsys = a 3Bsys e 3Bsys i 3Bsys 3Bsys ω 3Bsys M 3Bsys T . It is interesting to note that this rotating reference system is equivalent to the synodic system used in the circular restricted threebody problem.
In this case the eccentricity unit vectorP 3Bsys , the semilatus rectum unit vectorQ 3Bsys and the unit vector to the third bodŷ r ′ 3Bsys have to be expressed with respect to the third body rotating system, through the following composition of rotationŝ where 3Bsys = 3Bsys, 0 − u ′ . By using Equations (13), the new expressions of A 3Bsys and B 3Bsys are found as that are function of A 3Bsys u ′ − 3Bsys , i 3Bsys , ω 3Bsys , B 3Bsys u ′ − 3Bsys , i 3Bsys , ω 3Bsys so that the doubly averaged potential loses the dependence on the right ascension of the ascending node, F 3Bsys, k e 3Bsys , i 3Bsys , ω 3Bsys = 1 2π 2π 0F 3Bsys, k A 3Bsys u ′ − 3Bsys , i 3Bsys , ω 3Bsys , B 3Bsys u ′ − 3Bsys , i 3Bsys , ω 3Bsys , e 3Bsys df ′ .
So equivalently to Equation (12), we can compute the terms 2 to 4 were we dropped the subscript 3BSys.
F 3Bsys, 2 = 1 32 2 + 3e 2 (1 + 3 cos (2i)) + 30e 2 cos (2ω) sin 2 i F 3Bsys, 3 = 0 F 3Bsys, 4 = 9 32768 8 + 40e 2 + 15e 4 (9 + 20 cos (2i) + 35 cos (4i)) +560e 2 2 + e 2 (5 + 7 cos (2i)) cos (2ω) sin 2 i + 5880e 4 cos (4ω) sin 4 i that are equivalent to the expression in El'yasberg (1967) and Costa and Prado (2000). If the partial derivatives of Equation (11) are inserted into the Lagrange form of planetary equations we get the expression of variation of elements double-averaged over one orbit evolution of the s/c and over one orbital revolution of the perturbing body (either the Moon or the Sun). Figure 3 compares the evolution of a high-altitude orbit using the single averaged dynamics (red line) and the double averaged dynamics (cyan line) in eccentricity and argument of the perigee with respect to the Moon plane. It is interesting to see the evolution in the phase space of eccentricity and argument of the perigee measured with respect to the Moon orbiting plane as done in Ely (2005) and Colombo et al. (2014a). For the initial conditions in Figure 3, the simplified model by El'yasberg (1967) and Kozai (1962) predicts a pure librational orbit (Ely, 2005) (magenta line) which, in reality, is corrupted by the coupling between Moon and Sun third-body effect and by the effect of J 2 (see the red and cyan lines). The double averaged propagation derived in Equation (12) is used for the Moon and the Sun in Equation (10) for obtaining the orbit propagation represented with the cyan line, while the single averaged propagation described in sections Luni-Solar Averaged Potential and Earth Zonal Harmonic Potential is used for obtaining the orbit propagation represented with the red line. The single and double average propagation are compared against the actual spacecraft ephemerides from the NASA Horizon system in blue. Both the single and the double averaged approach show very good accuracy against the real ephemerides. With respect to the pure librational loop predicted by the Lidov-Kozai dynamics, still a quasi-librational behavior can be noted and will be further studied in the next section through propagation via the single-averaged dynamics.

LUNI-SOLAR AND ZONAL HARMONICS MAPS
To the purpose of studying the long term evolution of many initial conditions, a grid was built in the domain of inclination, eccentricity and argument of the perigee. Equally spaced steps in initial eccentricity (19 steps between 0.05 and 0.9), initial inclination (20 steps between 0.5 • and 90 • ) and initial right ascension of the ascending node (36 steps between 0 • and 180 • ) were selected as starting points. Note that, inclination and argument of the perigee are here described with respect to the Moon plane reference system; in other words i is the inclination of the spacecraft orbit with respect to the Moon's orbit plane and ω is the argument of the perigee measured from the direction of the ascending node of the spacecraft's orbit with respect to the Moon's orbit plane. Each initial condition on the grid is propagated over 30 years with the tool PlanODyn (Colombo, 2016) using the single averaged dynamics. As mentioned before, only luni-solar and zonal harmonics perturbations are here taken into account as we want to analyze their interaction. For each initial condition we evaluate the change between the minimum and the maximum eccentricity that the spacecraft will attain during its motion (see Figure 4): with setting t disposal = 30 years. As introduced the idea will be then to select limited e for graveyard disposal orbits or, at the opposite, maximum e orbits can be exploited for disposal through re-entry or for passive orbit transfer by exploiting the effects of perturbations. Figure 5 represents the result of the forward (a) and backward (b) integration for 30 years using, as initial conditions, the semi-major axis of XMM-Newton orbit, i.e., 67045.39 km, and different values of initial eccentricity and argument of the perigee in the grid, measured in the Moon reference system. The initial condition in terms of inclination and eccentricity corresponding to the one of the XMM-Newton orbit is represented by a star symbol. The maps are colored according to the maximum change of eccentricity e = e max − e min during the 30-year forward or backward propagation, respectively. If, for some initial conditions, the maximum eccentricity reaches the value of the critical eccentricity e critical , corresponding to a perigee of h p, re-entry = 50 km e critical = 1 − R E + h p, re-entry a the integration is terminated and the corresponding initial condition is marked with a cross symbol in Figure 5. Note that, for those solutions, the actual orbit evolution should be computed considering the effect of aerodynamic drag. This was not done in the current work to limit the computational time, however, we expect that the effect of drag will act as a dumping of the dynamical system, decreasing the amplitude of the librational or rotational loops in the (e, 2ω) phase space and progressively decreasing the semi-major axis, as shown in Colombo and McInnes (2011) for the case of solar radiation pressure, Earth's oblateness and drag. This will be subject of future work. Now, if we compute the maps considering both a forward and a backward propagation, both for 30 years, the map becomes more symmetric, as no particular choice (in terms of symmetry in the position of the Sun and the Moon) was made for the initial epoch. This is due to the fact that, if in the propagation the eccentricity reaches the maximum value e impact corresponding to a perigee equal to the Earth radius, the integration is terminated. Performing, from each starting point in the grid both a forward and a backward integration in time, instead, allows characterizing the dynamics in the phase space (e, 2ω) in terms of reachable conditions. Figure 6A shows the forward-backward maps starting from a semi-major axis of 67045.39 km, where now e = e max − e min is computed as setting t disposal = 30 years. The map in Figure 6A present and island of low eccentricity variation close to the initial condition of XMM Newton's orbit. The condition of quasi-frozen orbit is located around 2ω 0 = 180 deg and e 0 ≃ 0.8. The solutions around this conditions librate around the limited-eccentricity orbit. Other two islands of small eccentricity are present in this map, one around 2ω 0 = 180 deg and e 0 ≃ 0.2 and another around 2ω 0 = 0 deg and high values of the initial eccentricities. Those solution corresponds to orbits which rotate is terms of ω but have a limited variation in terms of eccentricity. Solutions starting at eccentricity close to zero have the highest variation of eccentricity. This means that for this semi-major axis and inclination, spacecraft on circular orbit can naturally get on an elliptical orbit and eventually reach re-entry. In order to approximate the half-period of the oscillation in the (e, 2ω) phase space, the time interval between the points when the spacecraft attain the minimum and the maximum eccentricity is computed (as an averaged between the forward and the backward propagation) and shown in Figure 6B. From Figure 6B it is for example possible to see that the quasifrozen solution is stable for 30 years. Indeed, due to the oscillation in inclinations, which has longer period, the orbit may encounter some instability is a longer propagation time is chosen as the eccentricity may increase beyond the critical eccentricity [this can be noted for example in the case of the INTEGRAL spacecraft (Colombo et al., 2014a)]. The maximum and minimum inclination attained with respect to the Moon's plane are shown in Figures 6C, D, while the variation of inclination is represented in Figure 6E. As it can be seen from Figure 6E the limited eccentricity conditions do not have a zero variation of the inclination; therefore, the long-term stability over longer time period is not guaranteed and an analysis over longer timeframe is required.
Once the behavior on the (e, 2ω) phase space is understood by looking at the map showing the result of the backward and forward propagation as in Figure 6A, it is possible to study the sensitivity to the initial condition in inclination with respect to the Moon's plane. In Figure 7, the e maps are computed for different starting inclination with respect to the orbiting plane of the Moon to show how the phase space changes depending on the initial inclination. For low inclinations the solutions are all rotational in ω (measured with respect to the Moon's plane) and exhibits a higher variation of eccentricity the more the initial eccentricity is high (Figures 7A-E) and for a high initial eccentricity the solutions which exhibits a higher eccentricity variation are characterized by 2ω 0 = 180 deg. This behavior is clearer in Figures 7F-L: the solutions initiating from the yellow area of high initial eccentricities and 2ω 0 = 180 deg rotate in ω, while quasi-equilibrium solutions exists at high eccentricities and 2ω 0 ≃ 0 deg (red island around 2ω 0 ≃ 0 deg). This is also visible from Figure 8 that shows the orbit evolution in the (e, 2ω) phase space for many initial conditions all at starting inclination with respect to the Moon's plane of 45 • and semimajor axis equal to 67045.39 km but for different starting ω 0 (color code) and for different initial eccentricity (Figures 9A,L).
Going back at the e maps in Figure 7, a new island appear for inclinations above 45 • at low eccentricities and 2ω 0 ≃ 180 deg of quasi-equilibrium solutions and librational solutions (Figure 7L), while the island at high eccentricities and 2ω 0 ≃ 0 deg move up. This is clearly visible from Figure 9 that shows the orbit evolution in the (e, 2ω) phase space and the polar plots for many initial conditions all at starting inclination with respect to the Moon's plane of 64.28 • and semi-major axis equal to 67045.39 km but for different starting ω 0 (color code) and for different initial eccentricity (Figures 9A, L).
For the semi-major axis considered (67045.39 km) i 0 = 45 deg is approximately the critical inclination (Kozai, 1962;Costa and Prado, 2000). Indeed, for inclination higher than the critical values, circular orbits get very elliptic (see yellow path in Figures 7L-T). The more the initial inclination with respect to the Moon increases, the more the initial orbit reaches the critical eccentricity for re-entry (cross symbols). Finally, for initial inclinations above i 0 = 70 deg, the quasiequilibrium solutions in correspondence of 2ω 0 ≃ 180 deg are in correspondence of e 0 > e critic , therefore not feasible, however, the fast eccentricity solution starting from e 0 = 0 still exist. It is important to remember that, even if qualitatively similar, the behavior depends also on the orbit semi-major axis that here is kept constant.

RE-ENTRY OR GRAVEYARD DESIGN
The initial conditions characterized by limited e variation identified in Figures 7L-T could be selected as graveyard orbits (Ely, 2005), while the high eccentricity variation solutions as initial condition for passive eccentricity increase to target Earth re-entry. It must be stressed that, while a re-entry trajectory will remove completely the spacecraft from the space environment, a graveyard solution will leave the satellite on a long-term stable orbit. In this sense, however, the time within which the orbit stability has been verified becomes an important parameter as the orbit may de-stabilize afterwards . On the other side, it must be noted that re-entry need to be carefully designed according to casualty risk constraints on ground (Merz et al., 2015).
This section will exploit the findings from the previous Sections to design the end-of-life disposal for XMM-Newton mission by enhancing the effect of the natural dynamics or lunisolar and J 2 perturbation. The approach for the optimal v computation was detailed in Colombo et al. (2014a) where the reentry of the INTEGRAL spacecraft was designed. The approach is summarized here and extended to the graveyard disposal design. For the disposal a single maneuver is considered, performed during the natural orbit evolution of the spacecraft (computed under the effect of perturbations). The finite variation in orbital elements α achieved by an impulsive maneuver, is computed through Gauss' planetary equations written in finite-difference form and the new set of orbital elements after the maneuver is propagated with PlanODyn (Colombo, 2016) in single-averaged orbital elements considering luni-solar perturbation and the zonal terms of the Earth gravity potential for the maxim disposal time t disposal = 30 years. Re-entry transfer orbits are selected that achieve in the propagation time t disposal the critical eccentricity e critical as in Colombo et al. (2014a), while for selecting a suitable graveyard disposal, solutions that attain the minimum e over the propagation time are selected. A graveyard orbit is designed imposing that after the maneuver, the variation of the eccentricity in time stays limited, that is e in Equation (14) is minimized. In order to analyze a wide range of disposal dates, different starting dates for the disposal were selected, whereas, to determine the maneuver magnitude v and direction (α and β) and the point on the orbit where the maneuver is performed f, an optimization procedure with genetic algorithm was performed in order to find the optimal set of parameters x = [ v α β f ] that minimizes the cost function J = e+w· v where w is a weighting factor set equal to 1,000. Rather than optimizing the timing for the optimal disposal maneuver; it was chosen instead to leave the time for that maneuver as a parameter of a sensitivity analysis; in other words, the insertion to graveyard was optimized for many different starting dates to analyze how the natural evolution of the orbit can be exploited. A maximum magnitude for the v is considered based on the available on-board propellant.

XMM-Newton Re-entry Disposal
The time interval considered for the disposal design is from 2013/01/01 to 2035/01/01. The maximum v available for the maneuver sequences is estimated to be 40.5 m/s in 2013/01/01 (Colombo et al., 2014a). The re-entry can be considered satisfied when an altitude of 50 km or lower is reached from the Earth's surface. The natural orbit evolution of the spacecraft is shown in Figure 10 in the (e, 2ω) phase space. The black points represent the initial conditions (and corresponding starting time) considered for the maneuver of disposal. For each starting point an impulsive maneuver of maximum magnitude equal to two times the available v on board on 2013/01/01 (equal to 81.1 m/s) is optimized in direction in order to maximize the following variation of the orbit eccentricity in the available t disposal . The results of the v optimization for each initial starting condition along the natural orbit evolution are reported in Figure 10 in colored line. Moving toward the external part of the phase space correspond to target the yellow regions in Figure 7P (for that map the initial semi-major axis and inclination are the one of the XMM spacecraft). Corresponding to large eccentricity variation. However, as it can be seen from Figure 10B, re-entry is not achievable within the considered t disposal for the maximum v considered (equal to 81.1 m/s) as the minimum perigee reached (among all the possible starting dates) is equal or above 4,700 km (too high for re-entry). This demonstrates that reentry for XMM-Newton is not a feasible option during this time range as the propellant requirements for such a disposal would be over the actual propellant on-board the spacecraft. Figure 10A represents the maneuver in the eccentricity-perigee angle (measured with respect to the Earth-Moon plane) phase space. Future studies for XMM-Newton disposal through reentry could investigate the possibility to increase the propagation time to verify whether the interaction between Moon and Sun third body perturbation will cause a natural decrease in the perigee. However, it needs to be taken into account that the available v on-board the spacecraft decreases with time due to orbit correction.

XMM-Newton Graveyard Disposal
Another disposal option that can be investigated for HEO in case the re-entry option is not feasible, is the option to transfer the spacecraft into a graveyard orbit. The existence of longterm stable orbits can be investigated, where the evolution of the orbital elements due to natural perturbation is limited. Such orbits can be chosen as graveyard orbits. Such orbits are visible in the map in Figure 7P (for that map the initial semi-major axis and inclination are the one of the XMM spacecraft) with a red color corresponding to a small variation of eccentricity. Importantly, note that the strategy would ideally aims at reaching the center of libration in the phase space; however, due to the limitation in the maximum available v, imposed as an upper bound for the global optimization, a more stable orbit cannot be reached, but only an orbit that is more stable than the nominal one. In other words, optimal solution is to move toward the center of the phase space loop. Also note that, due to the chaotic behavior of the orbit evolution under the effect of luni-solar perturbation, a driving factor is the time period t graveyard used for the propagation of the orbit after the maneuver to compute Equations (15) and (14). For this study t graveyard was set equal to 30 years, but this number can be easily increased for further work. Figure 11A shows the optimal maneuver for a transfer into a graveyard orbit for each starting time analyzed. The maneuver is represented in the phase space of eccentricity, inclination and anomaly of the percenter with respect to the Earth-Moon plane. The magnitude of the maneuver is always close to the upper bound of available v as is clear that a higher v would allow reaching a more stable orbit. However, the new graveyard orbit reduces at least the oscillations in eccentricity, preventing the spacecraft from an uncontrolled re-entry within the 30-year period. As an example, a disposal trajectory, whose maneuver is performed on 20/04/2016, is shown in Figure 11B.

CONCLUSION
This article analyzed the effect of luni-solar perturbations and the Earth's oblateness on the stability of highly elliptical orbits. The disturbing potential of the third body perturbation is written in Taylor expansion of the distances to the third body. The potential is firstly averaged over the revolution of the spacecraft around the Earth (mean anomaly) then is averaged again over the revolution of the perturbing body around the Earth. The existence of quasi-frozen, librational and rotational trajectories foreseen by the Kozai's analytical theory are found also when the Sun third body effect and the Earth's oblateness are included in the simulation. Maps are constructed over a wide domain of initial conditions in terms of eccentricity, inclination and argument of the perigee with respect to the Moon's plane. These maps represents the change in eccentricity over a long time space and in general can be used to study the orbit stability properties. These findings are finally used to design the end-of-life disposal for the XMM-Newton spacecraft as graveyard orbit injection or Earth re-entry. In this work, the effect of tesseral harmonics was not taken into account, while this could be important for creating maps in correspondence of lower semi-major axis, for example for Geostationary transfer orbits; the same approach however could be followed.

DATA AVAILABILITY
The datasets generated for this study can be found in the repository at the link www.compass.polimi.it/publications (Colombo, 2018).

AUTHOR CONTRIBUTIONS
CC contributed to the conception and design of the study by writing the PlanODyn code and simulating the maps contained in this paper. She wrote the first draft of the manuscript firstly presented at the 25th AAS/AIAA Space Flight Mechanics Meeting, in Williamsburg (VA) in January 2015.

Luni-Solar Perturbation Development
The average disturbing potential due to the third body effect can be written as Kaufman and Dasenbrock (1972): where µ ′ is the gravitational coefficient of the third body, r ′ is the spacecraft distance to the third body. The factors of the Taylor expansion can be written as function of the power of the ratio between the orbit semi-major axis and the distance to the third body δ = a/r ′ multiplied by the function R k (A, B, e), which is function of the orbit eccentricity and the function A and B (Chao-Chun, 2005) that are here written in terms of A Blizer and B Blizer given by Blitzer (1970) to show the link between them.
A = cos ωA Blizer + sin ωA Blizer B = − sin ωA Blizer + cos ωA Blizer A Blizer = cos cos u ′ + cos i ′ sin u ′ sin B Blizer = − cos i cos u ′ sin − cos i ′ cos sin u ′ + sin i sin i ′ sin u ′ The terms F k (A, B, e) in Equation (16) are given by Kaufman and Dasenbrock (1972) and are here reported in a more compact form up to order 6th as:  The derivatives of Equation (16) with respect to the orbital elements need to be computed to be inserted into the Lagrange planetary equations Equation (4) as in Equation (9) (Blitzer, 1970): C Blizer = sin i cos u ′ sin − cos i ′ cos sin u ′ + cos i sin i ′ sin u ′ .