Towards a sustainable exploitation of the geosynchronous orbital region

In this work the orbital dynamics of Earth satellites about the geosynchronous altitude are explored, with primary goal to assess current mitigation guidelines as well as to discuss the future exploitation of the region. A thorough dynamical mapping was conducted in a high-definition grid of orbital elements, enabled by a fast and accurate semi-analytical propagator, which considers all the relevant perturbations. The results are presented in appropriately selected stability maps to highlight the underlying mechanisms and their interplay, that can lead to stable graveyard orbits or fast re-entry pathways. The natural separation of the long-term evolution between equatorial and inclined satellites is discussed in terms of post-mission disposal strategies. Moreover, we confirm the existence of an effective cleansing mechanism for inclined geosynchronous satellites and discuss its implications in terms of current guidelines as well as alternative mission designs that could lead to a sustainable use of the geosynchronous orbital region.

contribution of satellites orbiting at geosynchronous altitude (GEO 1 ). GEO satellites have a semi-major axis of about R GEO = 42165 km and orbit about the Earth at the same rate that the Earth rotates around its axis. Due to this commensurability, the geosynchronous region provides us unique groundtracks that have been heavily exploited since the beginning of the Space Era. However, this came with a toll, the region about Earth's Clarke belt has been populated with artificial objects, most of them debris. GEO is nowadays the second most populated orbital region, after the low Earth orbits (LEO), with a couple of thousand catalogued objects and the list is continuously growing 2 .
The importance of the GEO region dictated to satellite operators to take measures at their missions' operational end-of-life. Indeed, even the very first GEO equatorial satellites applied some kind of re-orbiting manoeuvres to clear the longitude slots for future missions [24]. In the early 90s, when the first disposal studies were performed for equatorial GEO satellites [14], the use of super-synchronous graveyard orbits was proposed as an economical and effective solution to reduce the collision probability in the region. For the following years this was the norm followed by the operators, and is used up to now, modified to take also into account the perigee height variations due to solar radiation pressure perturbations [21,13]. The aforementioned ideas have sculpted the Inter Agency Debris Coordination Committee (IADC) set of mitigation guidelines for decommissioning GEO spacecraft [30] and the European Space Agency (ESA) instructions set for the ESA-operated GEO missions [23]. The guidelines suggest that a decommissioned GEO space system should: 1) increase the altitude of perigee by 235 km plus a factor accounting for the solar radiation pressure perturbations and 2) to circularise the graveyard orbit such that the eccentricity is of the order of 10 −3 .
However, from a sustainability point of a view, the use of graveyard orbits as suggested by the current mitigation guidelines, will keep increasing the collision probability in GEO [32]. Moreover, some further considerations support the need to investigate alternative ways of GEO exploitation. For instance, defunct satellites can act also as sources of secondary debris, even if there are no collisions. Satellites in graveyard orbits are subject to solar radiation torques that can constantly speed up their rotation [1]. The structural integrity of a space system rotating at several times per minute is not guaranteed. In fact, a population of high-area-over-mass (HAMR) object observed in GEO [52] was identified as multi-layer insulation (MLI) foils from satellites in graveyard orbits. The dynamics of these objects are quite different from the parent satellites [54] due to their higher Area-to-mass (A/m) ratio. Their long-term orbital evolution suggests that they cannot be contained in the nominal graveyard orbit and they could potentially cross the GEO protected region 3 .
Alternative solutions would require to either try to manoeuvre the satellite to a heliocentric orbit, but it is a cost-inefficient solution, or try to achieve atmospheric re-entry via natural perturbations. Unfortunately, the most effective perturbation that leads to the re-entry of close Earth satellites, namely the atmospheric drag, is not present at GEO altitude. Nevertheless, re-entry within reasonable time-scales could be achieved also by exploiting other type of perturbations. In their study of LEO satellites [3,4] suggested the possibility of re-entry due to solar radiation pressure resonances, but unfortunately this mechanism is also not effective for GEO. A more promising strategy would be to exploit the lunisolar gravitational perturbations [48,19], an idea that has been already been suggested for satellites in the Medium Earth Orbits (MEO) region [2,49,53,5] and Highly Elliptical Orbits (HEO) [31,15]. A living example of a mission that exploited natural lunisolar effects for its re-entry strategy is the INTEGRAL spacecraft in HEO. ESA operators manoeuvred INTEGRAL in 2015 such that it will perform a safe re-entry in 2029. The optimal manoeuvre design to enhance the effect of lunisolar perturbations [17,6] was used as starting point for the definition of the re-entry trajectory sequence [46].
In order to assess the possibility of re-entry due to lunisolar perturbations, the natural evolution of orbits at geosynchronous altitude has to be well understood. The dynamics of the geosynchronous equatorial orbits have been studied in the literature, both analytically and numerically [41,8,11,25]. Some aspects of the dynamics of the highly inclined GEO orbits have also been presented [20,55,56]. In this work, we will revisit the dynamics of geosynchronous region, in terms of searching for feasible re-entering highways or stable graveyard orbits when re-entry is not feasible. Having in mind a propagation span of at least 120 years, we employ an efficient semi-analytical propagator [16] that takes into account all the relevant forces at geosynchronous orbits in a single-averaged formulation (see Appendix A). The effect of Earth's precession, which is usually omitted, but turns out to be relevant for long-term propagations at high altitudes [29], has also been efficiently formulated and included (Appendix B).
Three different mapping methods are selected to highlight the different aspects of the dynamics. First, we study the contribution of the tesseral effects on the eccentricity growth for geosynchronous orbits. Then we employ dynamical maps over the angle-like variables, namely the argument of the perigee and the right ascension of the ascending node of the satellite, to explore the contribution of the different perturbations and their interactions. Finally, a global picture of the geosynchronous region in action-like variables (eccentricity and inclination) is presented to identify the orbital regimes where re-entry solutions are viable.
A natural separation is observed in terms of the natural evolution of initially low-and high-inclined orbits. Long-term stability and low eccentricity variations is the norm for low inclinations while an abundance of re-entering orbits exist in high inclinations. The eccentricity growth maps allow us to identify particular re-entry orbits that are interesting for future mission planning. We study the lifetimes of those orbits and try to identify the conditions that could lead to fast re-entry. Another interesting interaction that is revealed is the interplay of solar radiation and lunisolar perturbations for low-eccentric orbits. We present a case where, despite the usual behaviour where opening a solar-sail at the end-of-life enhances the de-orbiting process, this is not happening for GEO orbits. Finally, a general assessment of the current guidelines in GEO is made based on the current population and the underlying dynamics. The case of the Sirius constellation provides a strong case why a single equation guideline is not adequate to fully regulate the dynamically rich geosynchronous region.
The paper is organised in the following way, in Sec. 2 the single-averaged model is presented, in Sec. 3 we present the results of the dynamical mapping of the geosynchronous region, in Sec. 4 the findings in terms of post-mission disposal planning are discussed and in Sec. 5 we present our conclusions.

Orbit propagation with semi-analytical techniques
The methods used to analyse the long-term orbital evolution of geosynchronous orbits are briefly discussed here. A detailed description of the force model is given in Appendix A. Aiming for a 120-years integration time span, a singleaveraged semi-analytical propagation was opted for, which is a typical practise for Earth satellite orbits. The PlanODyn [16] orbital analysis suite was adopted and further developed to include the relevant perturbation effects. For the main force model, the contributions of Earth's geopotential, third body perturbation from the Sun and the Moon and the effect of solar radiation pressure were considered. Additionally, a revised version of Earth's general precession contribution to the secular evolution was developed and it is presented in the Appendix B.
For the geopotential force a fourth order and degree truncation is chosen, on the grounds that adding higher order harmonics did not produce any noticeable effects. For the zonal harmonics the first order averaged contributions were considered [36] and the second order contribution due to J 2 (J 2 2 ) [9] was also included. On the other hand, for the tesseral effects, only the resonant contributions relevant for geosynchronous orbits were considered [36].
The third body potential is expanded up to fourth order in the parallactic factor and is averaged in closed form over the mean anomaly of the satellite [34,15]. This is more efficient computationally [40] instead of using a series expansion representation (see [12] and references therein). The positions of the perturbing bodies, i.e. the Sun and Moon, are computed from analytical time-series truncated to a sufficient order for our work [45]. The solar radia- tion pressure is also included under a cannonball approximation [35,39]. The Earth's shadow is not taken into account since the constant Sun exposure is valid at geosynchronous altitude. Finally, Earth's general precession is also considered (Appendix B).
In Fig. 1 a validation of the full force model for the GEO region is presented, including all the above mentioned contributions. An initial condition for a GEO equatorial satellite (a = R GEO , e = 0.01, i = 0.1 • , Ω = 10 • , ω = 50 • , M = 0 • and initial epoch 21/06/2020 at 06:43:12.0) is propagated for 120 years with the semi-analytical method (red line) and is compared with a high-fidelity (grey line) integration. The high-fidelity model is in Cartesian coordinates, using Cunningham's algorithm [47] for the geopotential calculations and NASA's SPICE tool-kit 4 for the ephemerides of the Sun and the Moon and to retrieve the transformation matrices related to the motion of Earth's rotation axis. The semi-analytical propagation used in this work is in excellent agreement with the high-fidelity one, validating the force model selection and the use of a single-averaged formulation. Moreover, the singleaveraged formulation is a few orders of magnitudes faster, which allows us to proceed with a massive and accurate characterisation of the phase space.

Dynamical Mapping
In this Section, a dynamical study of the GEO region is performed, with the goal of presenting a detailed and complete study of post-mission disposal opportunities. Our work focuses on two main elements: the study of the GEO area in its entirety with a complete force model and the search not only for good graveyard solutions, but also for potentially exploitable re-entry trajectories.
In this work we considered only dynamical indicators associated with the eccentricity evolution of the orbit, which is mainly driving the perigee varia-tions and is the key in studying and planning post-mission disposal strategies. Moreover, other type of dynamical indicators based on the exponential divergence of nearby orbits, could fail to provide a clear picture for our purposes. Namely, in the region of the separatrix of the geosynchronous resonance, it is known that there exist chaotic behaviour [8,54,11]. However, in disposal design applications, we are not interested in the chaotic evolution with respect to any of the action variables, but rather in the chaotic behaviour manifesting as eccentricity growth. The fact that the use of chaotic indicators in disposal design can be missleading was also mentioned in a recent work [18]. Therefore, we focus our studies on the eccentricity variations and more specifically, two types of indicators are used. First the typical diameter of the eccentricity defined as the absolute difference between the minimum and the maximum values of the eccentricity obtained along the propagation time span. The second indicator, mainly used in Sec. 3.3, is a normalised eccentricity diameter and is designed to study the eccentricity variations along a wide range of initial eccentricities and/or semi-major axis. It is defined as [26]: where e 0 is the initial eccentricity, e max its maximum value along a propagation and e re−entry is defined as the value of eccentricity for which the perigee value becomes equal to a re-entry condition. Assuming the re-entry condition of 120 km above the surface of the Earth (a re−entry = R Earth + 120 km) for a satellite at the GEO region the re-entry value for the eccentricity is e re−entry ≈ 0.846. The behaviour of ∆e then is the following: when the eccentricity stays bounded about its initial value ∆e → 0, while if the eccentricity grows enough to reach the re-entry value ∆e → 1. In other words, ∆e gives the fraction of the physically available phase-space that the eccentricity evolution has covered. For the numerical exploration the single-averaged formulation was propagated numerically, using a BulirschStoer numerical integrator with automatic time-step control, imposing a local relative tolerance of 10 −13 . A propagation is terminated if the re-entry condition is met. Studying the full 6-dimensional space of orbital elements for all possible configurations could be a really involved task, so instead we have selected to study appropriate slices of it. Various 2-dimensional grids of initial conditions were chosen with a resolution of 201 × 201 orbits, yielding 40400 orbits per plot presented here. For each orbit, the computational time is a few seconds for a 120-years time span, giving less than a day computational time on average for each map 5 . A total amount of more than 50 million orbits were propagated in the GEO region, for the needs of the ReDSHIFT project [51].
For the selected grids of initial conditions, two values of the effective Areato-mass ratio was used, a low one A/m = 0.012 m 2 /kg, which is the typical value for decommissioned satellites, and a high one A/m = 1.0 m 2 /kg, which represents a satellite that carries an area-augmenting device (i.e. a solar sail kit) that is deployed at its end-of-life. The reflectivity coefficient was constant and equal to c R = 1. Finally, the initial epoch was selected on the 21/06/2020 at 06:43:12.0.

Tesseral maps
The first effect to be explored is the interaction between the GEO tesseral resonance and the other perturbations. For this reason, a series of eccentricity diameter maps on a two dimensional grid of initial semi-major a and the 24hour resonant angle λ = ω + Ω + M − θ g were produced, with ω, Ω, M the argument of perigee, the right ascension of the ascending node and the mean anomaly of the satellite respectively, and θ g is the Greenwich hour angle. There exist several ways to vary the resonant angle, but it was decided to do so only by varying the mean anomaly at the initial epoch. This selection allows to focus on the different dynamics induced purely from the tesseral contribution. In fact, varying the initial ω or Ω, would not complement the study as changing those two angles would also affect the initial configuration with respect to the lunisolar and solar radiation pressure perturbations and would contaminate the results.
In Fig. 2 the mapping for a typical satellite with A/m = 0.012 m 2 /kg is presented. Both a low-inclined (left column) and high-inclined (right column) configurations are presented. The first thing to mention is the difference in the order of magnitude of the eccentricity diameters with respect to the inclination. The low-inclined orbits exhibit a variation of the order of 10 −3 , while for the highly inclined ones it is two orders of magnitude larger. A low-eccentricity (top row) and a high-eccentricity (bottom row) initial condition is also presented for each case. In all cases we observe a clear distinction between the dynamics within the geosynchronous resonance and outside of it. The low-eccentricity and low-inclination is the only case where the eccentricity variation is larger than the surrounding altitudes. The separatrix of the tesseral resonance is also obvious in all the cases, with the two stable equilibria 6 at λ = 74.94 • and λ = 254.91 • and the unstable at λ = 161.91 • and λ = 348.48 • (red points in Figs. 2 and 3). Moreover, the width of the resonance is reduced from 80 km in the low-inclination case to about 50 km in the high-inclination case. The splitting of the separatrix is also apparent, particularly in the low-inclination high-eccentricity case, due to the inclusion of odd tesseral-harmonics in the equations of motion. Finally, for the high-inclination case, we observe smaller eccentricity variations in the low initial eccentricity case than in the high initial eccentricity case, which is an effect due to the lunisolar perturbations.   Fig. 2 Dynamical maps on the semi-major axis -geosynchronous resonant angle (a − λ) plane for the low area-to-mass case A/m = 0.012 m 2 /kg. The left column corresponds to low inclination orbits (i = 10 • ) and the right column to the high-inclination ones (i = 60 • ). The top row shows the evolution of low-eccentricity (e = 0.01), while the bottom row that of high-eccentricity orbits (e = 0.2). The colormap corresponds to the value of the eccentricity diameter over 120 years. Notice the narrow range between the minimum and maximum eccentricity diameter among all the maps. Fig. 3 shows a similar mapping for the case of a GEO satellite with an area-augmenting device that allows an A/m = 1.0 m 2 /kg at the end-of-life. In this case, the eccentricity variations in low-inclination case are an order of magnitude higher than in previous case, while in the high-inclination case the variations are similar. The other features are similar, again with the separatrix apparently dividing the phase-space, the width of the resonance decreasing with the inclination and the low-eccentricity high-inclination case exhibiting higher eccentricity variations than the high-eccentricity high-inclination case. Finally, the additional structures that appear as curvy lines above and below   Fig. 3 Dynamical maps on the semi-major axis -geosynchronous resonant angle (a − λ) plane for the high area-to-mass case A/m = 1.0 m 2 /kg. The left column corresponds to lowinclination orbits (i = 10 • ) and the right column to high-inclination ones (i = 60 • ). The top row shows the evolution of low-eccentricity (e = 0.01), while the bottom row that of high-eccentricity orbits (e = 0.2). The colormap corresponds to the value of the eccentricity diameter over 120 years. Notice the narrow range between the minimum and maximum eccentricity diameter among all the maps.
the separatrix are associated with secondary resonances related to the solarradiation pressure.
The phase-space exploration of the geosynchronous resonance reveals some interesting features of the phase-space, the most prevalent being the separatrix that links the two unstable equilibria. However, one should also notice that the overall differences in the eccentricity variations along the same map are very small. From our investigation it is apparent that the chaos detected in previous works [8,11,54] both for low and high area-to-mass objects is not resulting in any exploitable eccentricity growth with respect to nearby orbits with the same phasing with respect to lunisolar and solar radiation pressure perturbations. Therefore, the conclusion of the resonant angle investigation is that, placing a satellite on the unstable equilibria of the tesseral resonance does not present any significant re-entry opportunities.

Disposal maps
After having excluded the position of the satellite within the geosynchronous resonance as a source of interesting re-entry possibilities, here we investigate the orbital configuration with respect to the Sun and the Moon. The study is performed over a grid of initial argument of the perigee and right ascension of the ascending node. A similar approach has been used for the disposal design in the MEO region [2,5]. The advantage of this kind of approach is that it allows, given a particular point in the action-like variables space (a, e, i), to study its re-entry properties based on the initial orientation of the orbit with respect to the perturbing bodies. In addition, it is a good tool to study the interactions between all the angle related effects caused by lunisolar perturbations and solar radiation pressure.
In Fig. 4 we present how the contributions of the different effects mesh up to create a final disposal map for the GEO region (a = R GEO ). The initial orbital eccentricity is e = 0.1 and the initial inclination is 40 • with respect to the equator. The driving force that shapes the eccentricity growth at this altitude is the gravitational lunisolar interaction. Indeed, the eccentricity variation due to the geopotential and the solar radiation pressure are 3 orders of magnitude smaller than the those from the third body perturbations. Basically, under the isolated third body dynamics, all orbits with initial node of 180 degrees are reaching eccentricity values close to e re−entry . The way, however, that the geopotential and solar radiation pressure affect the combined effect evolution, is by fixing the frequency of the perigee oscillations. The tuned frequency of perigee, could suppress the Lidov-Kozai type dynamics [42,37] induced by the combined solar and lunar attractions. Therefore, the result of the evolution under the full dynamical model is quite complex and produces interesting dynamical structures.
Another interesting feature of those maps is that the position of the instabilities is mainly associated with the orientation of the node of the satellite with respect to the node of the Moon at the starting epoch. Therefore, changing the starting epoch could horizontally shift the appearing structures [2] and this feature repeats itself with a period of about 18.6 years, which is the nodal precession period of the Moon (also known as the Saros cycle). This adds a third dimension in the post-disposal design scheme and opens up for interesting design opportunities [49]. Namely, at the end-of-life one could wait for the value of the node to take correct value to maximise the effect of the lunisolar contributions.
In Notice that there exist significant differences between the low and high A/m cases, showcasing the importance of its contribution for low-eccentricity orbits. More specifically, at low-inclination the increased A/m forces only two stable configurations, those for ω +Ω −λ sun = 0 or π at a geosynchronous altitude, where λ sun is the ecliptic longitude of the Sun. For moderate inclinations, we observe that the solar radiation pressure seems to enhance the instability domain induced by the lunisolar perturbations, but this is not always the case (see also Sec. 4.2). Finally, for the 75 • inclination, due to the solar-radiation pressure, the stable perigee configurations do not exist any more and have been replaced by two stable nodal configuration for each perigee.
It is interesting now to examine the case of a higher initial eccentricity e = 0.2. In the panels presented in Fig. 6, in the same fashion as in Fig. 5, the different columns correspond to increasing values of the initial inclination (left to right) and the different rows to increasing values of the A/m (top to bottom). The behaviour with respect to the inclination increase is similar to the low-eccentricity case, i.e. as soon as the inclination exceeds the 40 • the orbits exhibits large eccentricity growth due to third body perturbations. On the other hand, in this case the effect of the solar-radiation pressure seems to be less profound. Indeed, except for the slightly enhanced eccentricity variations in the low-inclination case, there do not seem to appear any other significant differences between the lower and the upper row maps. From the dynamical maps presented here, one can deduce that for equatorial GEO satellites graveyard disposal is the only option. In our disposal mapping it is easy to identify the lowest perigee variation corridors (dark blue lines in left panels in Figs 5 and 6). This set of orbits should be targeted with post mission disposal manoeuvres, although this is not enough for a long-term safe graveyard. An adequate spacing between the disposed satellites should be ensured, such that the collision probability becomes minimal. On the other hand, for inclinations higher than about 40 • there exist an abundance of reentering solutions. The angle dependence of the position of unstable structures is not trivial at all, and poses interesting problems in re-entry disposal design which we will discuss in Sec. 4.
Another interesting aspect that we would like to highlight, is the effect of the higher A/m ratio in the low-eccentricity region. This is connected with the existence of a stable equilibrium of the solar-radiation pressure force at low eccentricities. On the other hand, for high initial eccentricities and inclinations, the evolution in the two A/m cases is almost identical.

Eccentricity-inclination maps
Although the complex dependence on the initial angles has already been discussed, here we attempt a global characterisation of the geosynchronous orbital region. We study the action-like variable space (a, e, i) and we address the angles dependence in a statistical manner, like in [26]. The semi-major axis is considered fixed and equal to the geosynchronous value for this mapping, since the same investigation for even up to 1000 km above or below R GEO yields very similar results. First, we present a set of maps for a fixed set of angles and then proceed with an angles-averaged dynamical mapping, i.e. for each point in the action-like variable space (a, e, i) we randomly select a sufficient sample of angular configurations (Ω, ω) and average the normalised eccentricity diameters over all the angles dataset. In Fig. 7 the eccentricity-inclination study for a = R GEO and two different angular configurations is presented. From disposal design point of view, a general feature of those maps that we should pay attention to is the generalised instability appearing at higher inclinations. And not only that, embedded in the unstable domain there exist particular configurations for which the eccentricity variations are small. Those regions of the phase space present some intriguing scenarios from future GEO exploitation, since they provide a stable operational regime next to an unstable region which could be used for end-of-life disposal (see for example the blue curves at 50 − 60 • inclination in Fig. 7).
As a last step, we present in Fig. 8 the angles-averaged dynamics over the eccentricity-inclination maps. In those maps, each point of the action-like space is sampled with 50 randomly selected angular configurations and the value of the dynamical indicator ∆e is averaged over all the angles. The result is a global dynamical map of the geosynchronous region, where the stable and unstable regions are clearly separated, albeit the information for the initial angle dependencies is lost. The region of above 40 − 45 • inclination presents a richness of re-entry solutions. Of particular interest is also the region from 65 − 75 • inclination where almost every orbit is re-entering within 120 years. The effect of the higher A/m (right panel) is almost negligible at higher eccentricities in the angles-averaged map. However, it creates some angle-related differences at low eccentricities.
The results of the angles-averaged study, convincingly confirm the natural separation of the phase space between low-inclined (below ∼ 40 • ) and high-inclined (above ∼ 40 • ) orbits. Namely, for low-inclined orbits, there is a natural deficiency of eccentricity growth orbits, and the search for stable graveyard solution is the only possible post mission design plan. However, at higher inclination another opportunity presents itself, an abundance of re-entering orbits exists. The global map indicates that for Earth satellites at geosynchronous altitude, the third body perturbations are prevailing over the other perturbations leading to large eccentricity variations for inclined orbits. The characteristics of those orbital highways and possible exploitation scenarios are discussed in the following Section.

Post-mission disposal issues
The findings of Sec. 3 are further discussed here, considering also issues related to post-mission disposal planning. We have already mentioned the abundance of highly-inclined orbits that have a considerable eccentricity growth within the 120-year time-span of our propagation. However, there is another crucial dynamical information that is hidden in the discussion of the previous paragraph, that being the orbital lifetime of each orbit. Of special interest are short-lived solutions that naturally re-enter the atmosphere. Examples of short-lived orbits are presented and their properties are discussed. Moreover, we present an interesting case where, even though an increased A/m usually enhances the de-orbiting process, not only this does not happen, but in fact the solar-radiation pressure effect cancels the real eccentricity growth mechanism. Finally, in the light of the dynamical mapping results, the current guidelines are discussed, and alternative ways of GEO exploitation are proposed.

An effective cleansing mechanism
In Sec. 3 we have encountered some orbits with exceptional short life-time. A typical example of this type of orbits is shown in Fig. 9, where the evolution of a trajectory with initial condition a = R GEO , e = 0.3, i = 63 • , Ω = 240.0 • , ω = 0.0 • , M = 0.0 • and A/m = 0.012 m 2 /kg is presented. The interesting feature is its orbital lifetime, which is less than 15 years. To exclude the chance that this is an outcome of the truncated force model or the single-averaged formulation, the same initial condition was also propagated under the highfidelity dynamics. The orbital evolution of two orbits coincides, suggesting that there is a quite effective cleansing mechanism at GEO altitude, that can make satellites re-enter even within the 25-year rule that is imposed for LEO orbits [30]. Moreover, the collision probability of an orbit like this is really minimal; for the solution shown in Fig. 9 the total time spend in the LEO protected region 7 as well as the dwell time in the GEO protected region is just a few days. Fast re-entering orbits were also reported in the literature [55,33] for high initial eccentricities and inclinations at geosynchronous altitude. Therefore, we would like to further explore under which conditions those orbits appear and study their properties. In Fig. 10 Fig. 11 The Ω, ω orbital lifetime maps for orbits at a = R GEO , eccentricity e = 0. Moreover, this set of orbits is not a local characteristic that happens only for the 63 • of inclination. Fig. 11 shows the lifetime disposal maps for initial eccentricity e = 0.2 and different values of the initial inclination. The fast re-entry patches exist in a range of initial inclinations from 50 − 90 • . However, their structure and their position with respect the initial value of the satellites node changes with varying inclination, due to the varying orientation of the perigee and node with respect to the perturbers' planes [37].
It is interesting now to understand the mechanism that leads those orbits to re-enter within 20-30 years. In order to further study this effect we select a set of orbits for i = 63 • , e = 0.2, ω = 60 • and values of the node Ω equally spaced every 10 • . The eccentricity evolution of this set of orbits is presented in This effect is even more clear if we look at the evolution of the eccentricity vector in the orbital plane through the set of variables e cos ω, e sin ω. It is now clear, that the evolution is following a Lidov-Kozai type of evolution, induced by the combined contribution of the Sun and the Moon. The interesting orbits, with fast re-entry times, are just those for which the eccentricity evolution allows them to reach the re-entry value within the first quarter of the dynamical evolution cycle. A suitable analytical method to check for the existence and a-priori locate their position is currently under development [27]. The insight developed from the study of the triple averaged Hamiltonian model suggests that, the in-plane dynamics for a range of inclinations with respect to the ecliptic become such, that the maximum eccentricity acquired during the Lidov-Kozai type of dynamics is equal to or larger the atmospheric re-entry value. Those initial conditions correspond to various sets of equatorial inclinations and node combinations, and could produce the complicated patterns that we see in the disposal maps in Fig 10 and 11.

Solar-radiation pressure implications
In Sec. 3 we concluded that the effect of the solar radiation pressure can be important for low-eccentricity orbits. Moreover, usually during a lunisolar driven re-entry, an enhanced solar radiation pressure would promote the deorbit process, as it has been observed for the transition region between LEO and MEO [50]. However, this is not always the case, especially for the inclined geosynchronous orbits where the the Lidov-Kozai type dynamics are driving the re-entry.
More specifically, by further inspecting both the disposal maps for loweccentricity and high-inclination and the eccentricity-inclination action maps,  Fig. 14 The 3-dimensional evolution of the orbits in Fig. 13 in the e, i , φ SRP space.
In the left panel, the orbit evolution for a standard spacecraft A/m = 0.012 m 2 /kg is presented, whereas in the right panel, the orbit evolution for a spacecraft equipped with an area-augmenting device A/m = 1.0 m 2 /kg is shown.
we encountered cases where opening a solar sail would considerably delay the de-orbit process. An example of this type of interaction is given in Fig.13. The low A/m orbit (blue curve) is re-entering within about 60 years of evolution, while the high A/m orbit (red curve) has almost double the lifetime. In an attempt to understand the delayed re-entry, we plot again the evolution in the orbital plane using the e cos ω, e sin ω variable for the two orbits. Immediately we recognise that the low A/m orbit directly follows a Lidov-Kozai type evolution. On the other hand, the high A/m orbit is trapped about the origin for a long time span, until it finally escapes and follows again a third-body dynamics dominated trajectory. This dynamical interaction is further explained in Fig.14 where we identify the main cause of the low-eccentricity trapping to be nothing else than the stable equilibrium of the solar radiation pressure resonance for the high A/m case. Namely, by defining the resonant angle for the solar radiation pressure as φ SRP = Ω + ω − λ sun , the evolution of the two orbits with respect to their eccentricity, inclination and φ SRP is presented. For the low A/m case, φ SRP is always rotating and the dynamics are following the eccentricity-inclination evolution dictated mainly by the third-body perturbations. On the contrary, in the evolution for the high A/m case, φ SRP is initially rotating but with the decrease in the inclination it is trapped into the resonance and is forced to librate about the stable equilibrium. The induced frequency in the argument of the perigee evolution, temporarily suppresses the Lidov-Kozai effect and the eccentricity stays bounded to low values. The further decrease of the inclination, finally drives the orbit out of the solar-radiation pressure resonance and it follows once again the Lidov-Kozai type of dynamics as we have also seen in Fig. 13.

Population, dynamics and guidelines
Currently, space operations in GEO are heavily concentrated in the loweccentricity, low-inclination region. This is evident in the right panel of Fig. 15 where the real population is superposed to the angles-averaged eccentricityinclination map. In this regime, mission planning and operations are well established, however there is not an efficient re-entry mechanism for inclinations lower than 40 • inclination. Hence, it is inevitable that the population of the space debris in the region will keep increasing.
On the other hand, carefully selected highly-inclined geosynchronous orbits can re-enter in time-scales comparable to the 25-years, which is the current IADC guideline for the LEO region. One could argue that inclined geosynchronous orbits are not as useful as the equatorial ones, but it has been shown that small constellations of a few eccentric and inclined geosynchronous satellites could reproduce the same coverage as a geostationary satellite [10]. This kind of exploitation has been already implemented successfully with the Sirius constellation. Sirius-1, 2 and 3 were operating in eccentric and inclined geosynchronous orbits for several years providing satellite radio services in North America. Unfortunately, at the end of the operational lifetime of the constellation, the operators following the current guidelines removed the three satellites from the GEO region, by reducing the semi-major axis by almost 10000 km and circularising the orbits. In the right panel of Fig. 15 the time evolution of the publicly available two-line elements (TLE) of three constellation components are shown and the end-of-life manoeuvres are highlighted in yellow. Nonetheless, the constellation was operating in an orbital region were fast re-entering orbits existed. We believe that considering re-entry as an alternative disposal solution, should not only be included as option in the guidelines for inclined satellites but also should be promoted. In this sense, the Sirius constellation was a missed opportunity to showcase a long-term sustainable exploitation of the geosynchronous region.
Another interesting idea to be explored in future geosynchronous mission design concepts is certainly the interplay between the solar radiation pressure and lunisolar perturbations. As we have seen in Sec. 4.2, using an area augmenting device can suppress the Lidov-Kozai effect and provide some low eccentricity variation operational orbits. A mission that uses the solar radiation trapping to stabilise and proceed to retract the area augmenting device at the end of the operations, could also lead to an atmospheric re-entry of the defunct satellite in a reasonable time frame.

Conclusions
The GEO orbital region was historically and is foreseen to be, one of the most precious assets in space exploitation. As it should be the norm with all natural resources available to humankind, the geosynchronous orbital region should also be treated in a sustainable way. Unfortunately, current practises in the region do not definitely ensure this.
In this work, a detailed dynamical cartography of the geosynchronous orbital region was performed to identify interesting possibilities in post-mission disposal strategies. Some of the key findings from the eccentricity variation maps in GEO are: 1) the positioning of the satellite only with respect to the geosynchronous resonance does not create interesting re-entry scenarios, 2) solar radiation pressure is important in the evolution of low-eccentricity lowinclination orbits, 3) for highly eccentric and inclined orbits the third-body perturbations dominate the dynamics, 4) there is a clear separation in the long-term evolution of low-and high-inclined orbits, implying that at geosynchronous altitude the Lidov-Kozai type of dynamics are prevailing.
Moving towards a sustainable GEO environment, mission design and planning should focus on exploiting fast re-entering orbits. Those are associated with particular geometries with respect to the Sun and Moon and they could be analytically located and introduced in the trajectory design process. Of course, this would require a whole re-assessment of the operations in the region. However, satellites in eccentric and inclined orbits could provide similar services to equatorial ones, but with the benefit that could achieve atmospheric re-entry at the end-of-life. Designing autonomous-navigation and propulsion systems to reach and follow those pathways is also a technological challenge but is well within the capabilities of future astrodynamical applications.
One of the first steps that we need to take in this direction, is to redesign the guidelines for decommissioning GEO satellites. In fact, the opportunity to apply this kind of re-entering strategies was presented in the past with the Sirius constellation. Unfortunately, the operators decided to follow the current regulations and re-orbited the satellites. What could have been a pioneer example of a clean exploitation of GEO is not enforced by current guidelines.
Moreover, given the dynamical complexities in the region, even the graveyard selection for low inclinations cannot be efficiently reduced to a single equation rule. Specialised tools that exploit the dynamical mapping of the region and sophisticated optimization algorithms should be used to provide the best-case disposal for each individual of post-mission disposal scenario. Considering the surrounding population in order to minimize collision probabilities of a given graveyard orbit should be also part of the disposal design process.

Geopotential
The perturbing function can be obtained in orbital elements through the classical Kaula expansion [36]: In the above expressions (l, m, p, q) are integers, F lmp (i), G lpq )(e) are the Kaula F and G functions [36], C lm ,S lm the non-normalised spherical harmonics coefficients of Earth's gravitational field and θ g the Greenwich hour angle. In the cases where the function G lmn is not a closed form function of e, it is expressed as a series up to O(e 20 ) for the needs of this work.

Resonant Tesseral harmonics
The resonant tesseral harmonics for the GEO case are also obtained from Kaula's expansion. Terms associated with J 21 (i.e. C 21 and S 21 ) are omitted because their strength is considerably smaller. The resonant contribution due to tesseral harmonics yields: where for each harmonic only terms satisfying the condition are selected. The resonant condition for satellites at geosynchronous altitude yields: l − 2p + q = m, with m = 1, 2, 3, 4 . . . .
Introducing the astronomical longitude as λ = ω + Ω + M − θ g the resonant contributions read: Third body perturbations The third-body potential implemented in PlanODyn is expanded in powers of the parallactic factor (a/r b ) as in [34,15], where r b is the geocentric distance of the perturber. Terms up to the fourth order (P 2 , P 3 , P 4 ) in the expansion are retained for both the Sun and the Moon.
The perturbing functions are single-averaged in closed form over the satellite's mean anomaly. This operation yields the following expressions for the disturbing functions [34,15]: where A and B are given from [34]:
Earth's general precession The long-term contribution of Earth's general precession is described in detail in Appendix B. The perturbing function yields: R prec = a (1 − e 2 ) µ ⊕ (P x sin(i) sin(Ω) − P y sin(i) cos(Ω) + P z cos(i)) where P = (P x , P y , P z ) is the angular velocity of a precessing geocentric equatorial frame relative to an inertial (i.e. J2000).

Equations of motion
The complete long-term evolution at geosynchronous altitude is described bȳ The equations of motion in orbital elements are then derived via Lagrange's planetary equations [7]: Appendix B: Semi-analytical modelling of Earth's general precession One effect that is usually overlooked in analytical and semi-analytical propagations is the effect of Earth's precession in the orbital evolution. Namely, for the Earth's gravity field representation, the lunisolar contributions and solar radiation pressure, a mean equator and mean equinox of the day frame is usually considered (MOD). All the ephemerides are also computed in this frame. However, this frame is not inertial due to Earth's precession. This effect is negligible for low and medium Earth orbits, but becomes important for GEO and higher-altitude satellites, especially for long-term propagations [29]. Particularly, it produces a long-term contribution to the inclination and right high-delity semi-analytical with precession semi-analytical without precession Fig. 16 Time evolution of the inclination for an orbit at about GEO semi-major axis, showing the contribution of Earth's precession on the evolution. The force model included only the J 2 effect and the transformation from MOD to J2000 frames. The three lines corresponds to: a) a high-fidelity propagation (grey line), b) a semi-analytical with only J 2 (blue line) and c) a semi-analytical that takes into account the Earth's precession contribution on top of J 2 (red line).
ascension of the ascending node that we would like to include in our simulations. In the literature, there exist different methods to model this effect [38,22,29].
Here, an efficient and accurate way to model the Earth's general precession effect is presented. The disturbing function for the Coriolis contribution in the non-inertial MOD frame is given by [28]: where µ ⊕ is the Earth's gravitational parameter, a, e the satellite's semi-major axis and eccentricity,ŵ is the unit vector, normal to the orbital plane, and P = (P x , P y , P z ) is the angular velocity vector of the non-inertial frame relative to an inertial one. In orbital elements the normal vectorŵ becomes: where i and Ω are the satellite's inclination and the right ascension of the ascending node. The disturbing function then reads and its contribution to the classical orbital elements, is given by Lagrange planetary equations: di dt prec = −P x cos Ω − P z sin Ω dΩ dt prec = P x cot i sin Ω − P y cos Ω cot i − P z dω dt prec = cos i (P z + P y cos Ω cot i − P x cot i sin Ω) , while the rest of the orbital elements (a, e, M ) are not directly affected. Keep in mind that the elements appearing in this formulation are non-osculating or contact elements. Since short-period terms do not appear in Eqs. 5 the equations of motion can be coupled with the rest of the single-averaged effects. However, one should pay attention when recovering the short-periodic terms of the mixed formulation [22].
To estimate the effect of the Earth's precession, we need now to estimate the angular velocity vector P of the MOD frame with respect to an inertial frame, which we assume to be the one defined by the mean equator and mean equinox of the epoch J2000. The rotation matrix from the MOD to the J2000 inertial frame is given from where R y , R z denote the clockwise rotation matrices with respect to the y and z axis respectively. The angles ζ, θ, z represent the general Earth's precession, i.e. the combined effect of lunisolar attraction on Earth's equatorial bulge and the change of Earth's orbital plane due to planetary attractions, and are computed from Lieske's theory [43]: ζ = 2306.2181 · T + 0.30188 · T 2 + 0.017988 · T 3 θ = 2004.3109 · T − 0.42665 · T 2 − 0.041833 · T 3 z = 2306.2181 · T + 1.09468 · T 2 + 0.018203 · T 3 where T is the Julian time past the 01/01/2000 in centuries and the angles are given in arcseconds. Now it suffices to observe that the angular velocity tensor Π between the two frames is given by: where theṘ is the time derivative and R is the transpose of R. Thus An illustrative example of this effect is given in Fig. 16, where the time evolution of the inclination for a hypothetical orbit at geosynchronous altitude for a time span of 120 years is presented. In order to make the Earth's precession effect clearer, we further account only for the Earth's oblateness (J 2 ) contribution. Under the single-averaged J 2 formulation (blue line), the mean inclination is constant. The same initial conditions are propagated also using a high-fidelity modelling (grey line), i.e. using a J 2 force model in Cartesian coordinates and NASA's SPICE tool-kit for the Earth's rotation matrices. The evolution now significantly differs, and a non-negligible contribution of Earth's general precession can be observed, producing an oscillation with an amplitude of about 8 seconds of a degree in the inclination evolution. However, this effect can be accurately included in the semi-analytical propagation, by adding the Earth's precession contributions (Eq. 5) (red line).