A frequency portrait of Low Earth Orbits

In this work we deepen and complement the analysis on the dynamics of Low Earth Orbits (LEO), carried out by the authors within the H2020 ReDSHIFT project, by characterising the evolution of the eccentricity of a large set of orbits in terms of the main frequency components. Decomposing the quasi-periodic time series of eccentricity of a given orbit by means of a numerical computation of Fourier transform, we link each frequency signature to the dynamical perturbation which originated it in order to build a frequency chart of the LEO region. We analyse and compare the effects on the eccentricity due to solar radiation pressure, lunisolar perturbations and high degree zonal harmonics of the geopotential both in the time and frequency domains. In particular, we identify the frequency signatures due to the dynamical resonances found in LEO and we discuss the opportunity to exploit the corresponding growth of eccentricity in order to outline decommissioning strategies.


Introduction
It is known that the proliferation of space debris in the Low Earth Orbit (LEO) region has already become a critical issue to handle. In this context, as part of the H2020 ReDSHIFT (Revolutionary Design of Spacecraft through Holistic Integration of Future Technologies) project [1,2], a deep analysis to search for passive deorbiting solutions in LEO was carried out, by performing an accurate mapping of the phase space in order to identify stable and unstable regions. A detailed description of the results of the LEO cartography was presented by the authors in [3], while in [4] a general analysis on the role that resonances induced by solar radiation pressure (SRP) can play in assisting the deorbiting was provided.
In general, the key idea investigated in those works was to identify the orbits, and the associated mechanisms, where dynamical perturbations can induce a significant growth of the orbital eccentricity, in order to facilitate passive disposal. Accordingly with our findings, we concluded that in the case of a typical intact object in LEO, with an area-to-mass ratio of the order of A/m = 10 −2 m 2 /kg, perturbations as SRP, lunisolar effects and high degree zonal harmonics cannot ensure the reentry on their own but only in combination with the atmospheric drag. If the spacecraft is, instead, equipped with an area augmentation device, which increases the effective A/m by, e.g., two orders of magnitude, then we concluded that SRP alone can drive the dynamics, if the initial inclination of the orbit is close enough to a resonant inclination, given semi-major axis and eccentricity.
In the present work, we make a deeper analysis of the role of resonances which act in the LEO dynamics, by characterising the eccentricity of a set of orbits in terms of periodic components. Starting from the quasi-periodic time series of eccentricity, computed for a dense grid of initial conditions, we decompose the series in the main spectral components by means of a numerical computation of the Fourier transform. Then, we link each frequency component with the dynamical perturbation responsible for that signature. In this way, we have an additional tool to explore the relative importance of each given gravitational or non-gravitational perturbation in LEO as a function of the initial orbital elements. The final goal of such analysis is to support the cartography in identifying the orbits where a significant growth of eccentricity, led by one or more perturbations, can assist the passive disposal of objects at their end-of-life. The same analysis can also serve to identify the periodic drifts that operational orbits could experience.
In the past, the study of the chaotic dynamics within the Solar System led [5] to devise a method for a numerical estimation of the size of the chaotic zones, based on the variation in time of the main frequencies of the system. Since then, the algorithm for the frequency analysis was developed to study the stability of the orbits in many multi-dimensional conservative systems, in order to provide a global representation of the dynamics [6,7]. The Frequency Map Analysis algorithm (Numerical Analysis of Fundamental Frequencies -NAFF) is based on a refined and iterative numerical search for a quasi-periodic Fourier approximation of the solution of the system over a finite time span [8]. Considering, in particular, the issue of optimal design of artificial satellite survey missions around a non-axisymmetric body, Noullez et al. [9] proposed an alternative method with respect to the standard Fourier transform approach to characterise satellite orbits by computing the periodic components in order to identify regular orbits, meant as orbits whose inclination and eccentricity do not vary significantly over a given time scale. Concerning in particular the LEO region, Celletti and Galeş [10] studied the dynamics of resonances in LEO with the aim of identifying the location of equilibrium position and their stability. Within the common scope of defining suitable post-mission disposal orbits, they studied analytically, by means of a toy-model, whether an object is located in a stable or chaotic region. In such a way, the identification of stable orbits in LEO suggests the detection of possible graveyard orbits. In this paper, we focus on the possibility of exploiting the eccentricity growth induced by one or more dynamical perturbations at given orbits to facilitate the end-of-life reentry and we deepen this analysis by characterising the eccentricity evolution in terms of its main frequency components. A comprehensive characterization of the dynamical evolution of the eccentricity is a key ingredient in order to identify, among other things, possible disposal strategies for operational and future spacecraft.
The paper is organised as follows: in Section 2 we briefly describe the dynamical model adopted for the numerical propagation and we introduce the method to identify the frequency signatures which characterise the eccentricity evolution of a set of LEO orbits. In Section 3 we outline the results of our analysis, comparing the results of numerical propagation in the time domain with the findings of the frequency characterisation. Finally, in Section 4 we draw some conclusions.

Dynamical model and methods
As mentioned before, within the scope of ReDSHIFT, we performed an extensive mapping of the LEO phase space by propagating more than 3 million orbits, as described in [11,12,3] 1 , spanning from 500 km to 3000 km of altitude over the Earth surface, considering a wide range of eccentricities, from 0 up to 0.28, and inclinations, from 0 • to 120 • , 16 different (Ω, ω) configurations and two initial epochs. In the following, we limit our analysis to the case of right ascension of the ascending node, Ω, and argument of perigee, ω, both equal to 0 • , with the initial epoch set to 21 June 2020. The orbital propagation was carried out over a time span of 120 years by means of the semi-analytical orbital propagator FOP (Fast Orbit Propagator, see [13,14] for details), which accounts for the effects of 5 × 5 geopotential, SRP (assuming the cannonball model), lunisolar perturbations and atmospheric drag (below 1500 km of altitude). Two possible values of the area-to-mass ratio were considered: A/m = 0.012 m 2 /kg, selected as a reference value for typical intact objects in LEO, and A/m = 1 m 2 /kg, a representative value for a small satellite equipped with an area augmentation device, as a solar sail [15]. More details on the adopted model can be found in [3]. The results of the cartography can be displayed in contour maps showing the lifetime or the maximum eccentricity over the propagation interval as a function of the initial inclination and eccentricity, for each initial semi-major axis. A large set of maps can be found on the ReDSHIFT website 2 . In the following Sections, some examples will be provided.

Dynamics in the time domain
We are particularly interested in studying the time evolution of the eccentricity. Indeed, within the search for passive disposal solutions in LEO, the identification of orbits which can experience a significant growth of eccentricity becomes crucial, since in this case the lowering of the orbital perigee helps drag in being effective. Moreover, a variation in eccentricity causes an altitude variation which could become an issue also at the operational stage, for instance in the case of a large constellation.
Lagrange planetary equations (e.g., [16]) show that SRP, lunisolar perturbations and high degree zonal harmonics 3 cause long term periodic variations in the evolution of eccentricity, which become quasi-secular in the vicinity of 1 All the papers related to the project are available on the ReDSHIFT website at http://redshift-h2020.eu/documents/.
2 http://redshift-h2020.eu/results/leo . 3 The oblateness of the Earth, J 2 , does not affect the evolution of the eccentricity over long term (e.g. [16]). Table 1: List of the main resonances expected to be found in LEO: argument ψ j , values of the coefficients α, β, γ and corresponding index j. Resonances from j = 1 to j = 6 are due to SRP; resonances 7 and 8 are singly averaged solar gravitational resonances; resonances from 9 to 11 are doubly averaged lunisolar resonances.
a resonance involving the rate of the right ascension of the ascending node, Ω, and the argument of perigee, ω. In particular, we can write the instantaneous variation of e due to a given perturbation in the general form: where T is a coefficient which depends on (a, e, i) according to the given perturbation and the argument ψ can be written in general terms as: where α, β, γ = 0, ±1, ±2 depending on the perturbation and λ S is the longitude of the Sun with respect to the ecliptic plane, set as λ S = 90.086 • at the starting epoch. A resonance occurs when the conditionψ 0 is satisfied.
The list of the resonances expected from the theory and found by means of the LEO cartography [3] are shown in Table 1, where the corresponding expression for ψ and the value of α, β, γ are highlighted, together with an index (j = 1, ..11) associated to each resonance. Resonances indexed from 1 to 6 correspond to the conditioṅ withᾱ = 0, 1, and are associated to the zero-order expansion of the SRP disturbing function (e.g., [17,18]). Resonances 7 and 8 are singly averaged solar gravitational resonances (e.g., [19,20]), while resonances from 9 to 11 are associated to doubly averaged lunisolar gravitational perturbations (e.g., [19]). The rate of Ω and ω can be found by applying the Lagrange planetary equations and accounting in principle for both the effects of J 2 and SRP, while the effect of lunisolar perturbations can be neglected (e.g., [21]). The explicit expressions have been given, for instance, in [4]. In practice, in [4] we have shown that, for an initial orbit with Ω = ω = 0 • , at the assumed initial epoch (which corresponds to λ S ≈ 90 • ), the rate of Ω and ω due to SRP vanishes. In Figure 1 we display the behaviour of |ψ j | for each perturbing term highlighted in Table 1 (j = 1, ..11) as a function of the inclination i ∈ [0 • : 120 • ] for the case of a quasi-circular orbit (e = 0.001) with a semi-major axis a = 7978 km. The figure shows that curves associated to different perturbations may intersect and overlap creating a dense network of resonances in the phase space. Thus, depending on the given inclination, it may be hard to distinguish between the concurrent effect of different perturbations and to link the dynamical effect to the perturbation which produces it. To overcome this problem, we can take advantage of the fact that the adopted orbital propagator is set up in such a way that each dynamical perturbation in the model can be individually turned on or off. Since we aim at identifying the specific effect of a given perturbation on the eccentricity evolution and at characterising it in the frequency domain, we consider two simplified models, which fit our purposes: • model I: SRP on; lunisolar perturbations and drag off; geopotential: only J 2 ; • model II: SRP off; lunisolar perturbations and drag on; geopotential: 5 × 5.
Model I is particularly suitable to study the SRP effects on the eccentricity in the case of high A/m objects, when only SRP and drag play a primary role in the evolution. In the case of A/m = 1 m 2 /kg, atmospheric drag is effective in driving a reentry within 25 years for pericenter altitudes up to 1050 km (see [3,22]). Since this is a relatively high value, in order to focus on the effect due to SRP, we have decided to switch off the perturbation due to the atmospheric drag.
Model II, instead, is appropriate to study the effects led by lunisolar perturbations and high degree zonal harmonics: removing from the model the presence of SRP, we avoid the chance of mismodelling, since the resonant inclinations corresponding to lunisolar perturbations and geopotential can be close to those associated with SRP, as appears from Figure 1. In this case, adopting the low or the high value of A/m does not affect the eccentricity evolution.

Frequency characterisation of the eccentricity
The starting point for the frequency characterisation is to process the discrete eccentricity time series of a given initial orbit to obtain the discrete Fourier transform through a standard Fast Fourier Transform (FFT) algorithm (e.g., [23]), based on the Cooley-Tukey algorithm [24]. The basic idea is to identify the frequency and the amplitude of the main spectral features in the frequency series. The criterion we adopt is to account for any signature whose amplitude is, at least, 10 times stronger than the mean value of the spectrum in the surrounding area.
A first issue to be considered concerns the time sampling ∆t of the input series to be transformed. Indeed, the sampling frequency is f s = 1/∆t and, from Nyquist theorem, it follows that f s /2 is the highest frequency we can capture from our analysis. Since the perturbations we are interested in have periodicity of the order of months to years 4 , the sampling ∆t = 1 day, adopted in [11,12,3], is fully reasonable. On the other side, a more critical issue involves the lowest detectable frequency by our analysis, which is limited by 2/T , where T is the duration of the time series. This means that with the adopted time span of 120 years, signatures with periodicity up to 60 years would be, in principle, identified. In practice, signatures due to perturbations with periodicity of more than some years are poorly sampled by definition. Thus, we propagate the set of orbits of interest for a longer time span, 600 years, in order to catch unambiguously signatures with periodicity of some tens of years, as expected in the vicinity of a resonance.

Analysis of the numerical results
The general results of the LEO dynamical mapping was already extensively described in [3]. In the following, we present the results obtained by assuming the two simplified dynamical models, described in Section 2.1. First, we consider the case of model I, i.e., we focus on the effect of SRP in the case of the augmented A/m ratio: we briefly recall the main findings in terms of time evolution of the eccentricity, then we discuss the results of the characterisation in terms of frequency components. Next, we present the same analysis in the case of model II, focusing on the effects of lunisolar perturbations and high degree zonal harmonics.

Model I 3.1.1. Analysis in the time domain
We recall that the model accounts, in this case, only for the effect of SRP and J 2 , while drag and lunisolar perturbations are turned off. We propagate the orbits assuming A/m = 1 m 2 /kg and we look for the inclinations where a growth of eccentricity due to SRP occurs. Some illustrative results are shown in Figure 2: on the left we show the maximum eccentricity achieved over 600 years of propagation as a function of the initial inclination and eccentricity, for initial a = 7978 km (top) and a = 8578 km (bottom), respectively. On the right panels we display the corresponding lifetime, in years. We recall that the atmospheric drag is effective up to 1050 km of altitude for the adopted A/m ratio. Thus, we selected on purpose two reference values for the initial semi-major axis which are significantly above the region where drag plays a role. If the effect of SRP is able to lower the perigee below 1050 km, then the removal of the drag from the model allows to check if the chance to reenter or not can be ascribed solely to SRP.
The lifetime panels show that, in the case an area augmentation device is available on-board, even for high altitude quasi-circular orbits, a reentry driven by SRP alone is feasible for inclinations in the vicinity of 40 • , which corresponds to the resonant condition: In the case of initial a = 7978 km, reentry can be achieved in about 7 years for initial e ranging from 0.0001 to 0.009 thanks to SRP alone, for an initial orbit at i = 39.5 • . In the case of initial a = 8578 km, SRP allows to reenter within 10 years at initial i = 37.5 • and in about 16 years for i = 38 • . The other resonances due to SRP, although not able to drive a reentry, cause, anyway, a remarkable growth in eccentricity, as can be seen from the left panels of Figure 2, which can be exploited to lower the perigee of the orbit. Referring also to Figure 1 in [4], which shows the location of the 6 main SRP resonances as a function of i and a, we can identify the following resonances corresponding to the bright inclination "corridors": •ψ 3 0 andψ 5 0 around i = 58 • and i = 54 • , respectively, in the top panel (a = 7978 km), while they intersect around i = 56 • at a = 8578 km; •ψ 4 0 andψ 6 0, both occurring in the vicinity i = 70 • .
Moreover, we can recognise other features at specific inclinations, appearing as fainter, but still visible, signatures. They can be associated to higher-order terms in the expansion of the SRP disturbing function (e.g., [17]): in Section 3.1.2 their identification will be assisted by the analysis in terms of frequencies.
For completeness, turning on the contribution due to the atmospheric drag in the model, we find that the synergic effect of SRP and drag can support reentry also at different values of inclinations (resonances) but, typically, only over long time scales. This is shown in Figure 3, in the case of an initial orbit at a = 7978 km and e = 0.001, assuming now model I with the further contribution of the drag. For the same initial orbit, Table 2 shows the lifetime associated to the initial inclination corresponding to the six SRP resonances. The table points out that the addition of the drag in the model can assist the reentry at inclinations close to the resonant ones, but only in the case of resonance 2 (in addition to resonance 1) the reentry can take place in less than 25 years.

Analysis in the frequency domain
The analysis of the maximum eccentricity maps (Figure 2 -left panels) shows that, in addition to the six resonances due to the zero-order expansion of the SRP disturbing function, other fainter signatures can be observed   Table 3: List of the first-order terms, expanding the SRP disturbing function up to firstorder (e.g., [17]): argument ψ j , values of the coefficients α, β, γ and corresponding index j.
at given inclinations. Thus, to build a complete picture of the eccentricity evolution in the LEO phase space we need to include the first-order terms in the expansion of the SRP disturbing function (e.g., [17]), which are listed in Table 3.
Following the procedure depicted in Section 2.2, we identified the main frequency signatures associated to the eccentricity, at each initial condition available. The frequency components detected at each inclination for the two illustrative cases of an initial orbit at a = 7978 km and a = 8578 km in the case of initial e = 0.001, with A/m = 1 m 2 /kg, are shown in Figure 4. Each square in the plot represents a detected frequency component; the color bar refers to the relative amplitude of the frequency signature 5 , intended as the corresponding intensity peak in the computed Fourier spectrum. Each coloured curve represents the behaviour of the argument |ψ j | as a function of the inclination, with a cusp at the resonant inclination. As it can be seen, the detected signatures match almost exactly the theoretical curves. We also point out that the amplitude of the signatures gradually grows approaching a resonant inclination. In particular, the effects of SRP first-order terms at given inclinations, which could be only partially inferred from the maximum eccentricity maps, can be clearly identified in the frequency chart.
The signatures detected by means of the frequency analysis match the bright corridors detected in the maximum eccentricity maps in the left of Figure 2. In particular, resonances 3 and 5 (see Table 1), which intersect for a = 8578 km, can be individually identified for a = 7978 km both in the contour map and in the frequency chart. Moreover, the frequency chart for a = 8578 km shows a signature around i = 86 • corresponding to the first-orderψ 17 term, which does not appear for a = 7978 km neither in the contour map nor in the frequency chart. Finally, the a = 7978 km chart shows a signature with singularity at i = 90 • which can be associated to the rate of Ω, appearing in the second-order expansion of the SRP disturbing function (see, e.g., [17]).
From the lifetime maps in the right panels of Figure 2, we know that only in the case of resonance 1 SRP alone can drive a reentry. Nevertheless, the maximum eccentricity maps show that in the vicinity of a resonance a certain growth of eccentricity occurs anyway. Thus, in the perspective of designing passive disposals and when dealing with operational issues, it is crucial to consider the timescale over which the eccentricity variation takes place. With this in mind, assessing the change in eccentricity led by a perturbation without performing the numerical propagation, i.e., by characterising the LEO phase space in terms of frequencies, represents a very powerful tool.
In [4], starting from Eq. (1) we showed that the maximum eccentricity variation achievable due to the zero-order SRP resonance j for a given initial (a, e, i) can be estimated as: On the other side, the amplitude associated to each detected frequency signature in the Fourier transform gives an estimate of the eccentricity increment, as well. Both these values can be compared with the numerically computed maximum eccentricity over 600 years: the three estimates are expected to comply with each other. A general comparison for the initial orbit at a = 7978 km and e = 0.001 is shown in Figure 5: as a function of i, we show the theoretical amplitude |T j /ψ j | with j = 1, ..6 for the six zero-order SRP resonances, the maximum variation in eccentricity, ∆e max , achieved over the numerical propagation (red circles) and the amplitude of frequency signatures detected by our analysis (filled squares; the color bar refers to the corresponding periodicity, i.e. the inverse of the detected frequency). A similar example for an orbit at a = 8578 km is shown in [25]. The match is very good; moreover, we can observe that, as expected, the brighter squares, associated to signatures with longer periodicity, are found only in the vicinity of resonances. Looking at the maximum variation in eccentricity achieved during propagation with FOP (red circles), some fainter features can be noticed at inclinations different from those corresponding to the six main resonances. Comparing the inclination of these signatures with the resonant inclinations corresponding to the arguments shown in Table 3, these fainter features can be associated to the first-order terms in the expansion of the SRP disturbing function.
In Figure 6, we show a detailed (i, e) zoom around the two main resonances found at this altitude:ψ 1 corresponding to i ∼ 40 • andψ 2 in the vicinity of i ∼ 80 • . The maximum eccentricity displayed on the y−axis corresponds to the eccentricity needed to lower the perigee down to 120 km, e 120km = 0.185. Both the squares corresponding to the numerical maximum eccentricity (left panels) and the amplitude of the frequency signatures (right panels) lie on the theoretical curves for |T 1 /ψ 1 | and |T 2 /ψ 2 |. This further confirms that the three quantities (theoretical amplitude, numerical maximum eccentricity and amplitude of the frequency signature) provide the same information, thus one can be adopted in place of the other.
We can notice, however, in the bottom panel on the left of Figure 6, a disagreement between the theory and the numerical propagation: according to the theory, the maximum eccentricity variation for initial i = 79 • should be sufficient to lead to reenter, while the ∆e max computed with FOP turns out to be lower than e 120km . The explanation for such a behaviour is that during the propagation also the inclination experiences a variation which moves the object away from the resonance, making the SRP perturbation less effective. In Figure 7 we show the evolution of e and i over 100 years for the initial condition a = 7978 km, e = 0.001, i = 79 • . Both eccentricity and inclination show a periodicity of about 28 years but they are out of phase: the eccentricity starts to grow led by the SRP perturbation; at the same time, the inclination starts to decrease so that when the eccentricity reaches the maximum value e max = 0.14, the inclination is at its minimum, i min = 78.3 • , where, as can be inferred from Figure 6, the perturbation due to the resonant term ψ 2 is no longer effective in driving the reentry.
This fact shows that the rate of i should be taken into account to provide a full description of this case based on the dynamics. It is beyond the scope of this work to provide a full description on this scenario based on the dynamical systems theory, but we can provide a basic tool to obtain an a priori indication on whether the orbit will exit from the resonance domain before achieving a reentry.
In Figure 7, in the panel showing the evolution of the inclination, it is also displayed the behaviour predicted by the theory developed in [26] for lunisolar gravitational resonances, which can be applied also in the case of SRP, as shown in [4]. In particular, it is demonstrated that there exists an integral of motion, corresponding to (β cos i − α) µa(1 − e 2 ) = constant, (6) where α, β are as defined in Eq. (2). In other words, assuming that the motion of the spacecraft is governed only by the Earth's monopole, the Earth's oblateness and the solar radiation pressure, at any instant we can recover the inclination value from where the constant can be obtained by evaluating Eq. (6) at the initial epoch. For completeness, in Figure 8 we show a comparison over 30 years of the eccentricity and inclination evolution computed by propagation assuming model I (blue curve) with the behaviour obtained by assuming the complete dynamical model (red curve), which includes all the perturbations provided by FOP. The initial orbit is the same as in Figure 7. We can observe that  the two models predict the same behaviour, except that, in the second case, the reentry is ensured (in 13.6 years) by the atmospheric drag.
In Figure 9, we show the behaviour predicted for the inclination by Eq. (7), by assuming a maximum variation in eccentricity as in Eq. (5), for resonances 1 and 2. We can notice that in the first case, when we consider an initial inclination in the resonance domain, the variation is not relevant if compared with the curves in the top panel of Figure 6). In the second case, the variation is instead important, of about 1 • and moves the dynamics towards the edges of the interval where the resonance is effective (compare with the curves in the bottom panel of Figure 6).
The above discussion showed that the assumption that T j is a function of the initial values of eccentricity and inclination may provide a misleading information. Figure 10 shows the evolution of ∆e 2 = T 2 (a, e, i)/|ψ 2 |, according to Eq. (5), assuming the values of eccentricity and inclination computed at each given time by propagation with FOP, in case of model I, for initial a = 7978 km, e = 0.001, i = 79 • . The y−axis upper limit corresponds to a perigee altitude of 120 km. As it can be seen, for the initial value of e and i, the growth of eccentricity ∆e 2 is such that the reentry driven by resonance 2 is feasible (the curve is not visible in the figure because it is higher than the eccentricity required to reentry). On the contrary, after only 5 years, the inclination has moved from its initial value (compare with Figure 7) enough that the corresponding growth in eccentricity due to resonance 2 alone is no more capable to assure the reentry. Finally, similarly to Figure 6, the comparison between theoretical amplitude, maximum variation in eccentricity computed with FOP and amplitude of the frequency signatures for a = 7978 km and e = 0.001 in the cases of resonances 3, 4, 5, 6 due to SRP is shown in Figure 11. Also in these cases the agreement is noticeable.

Model II 3.2.1. Analysis in the time domain
Model II is particularly suitable to study the perturbation on eccentricity due to lunisolar effects and high-degree terms in geopotential, since SRP has been removed in this case. The effective area-to-mass ratio of the object does Figure 11: Comparison between theoretical amplitude |T j /ψ j | (j = 3, 4, 5, 6), maximum variation in eccentricity over propagation computed with FOP (left panels) and the frequency amplitudes detected by our analysis (right panels) as a function of the inclination, for the initial orbit at a = 7978 km and e = 0.001 in tha case of model I.
not play a role in driving the dynamics, contrary to the case of the previous model, thus we assume A/m = 0.012 m 2 /kg for simulations.
In analogy to the left panels of Figure 2, Figure 12 shows the maximum eccentricity as a function of the initial inclination and eccentricity for an orbit at a = 7978 km (left) and a = 8578 km (right), respectively. In this case, we do not show the corresponding lifetime maps: at these altitudes and for quasi-circular orbits the maps would result blank since neither lunisolar perturbations nor high-degree terms of geopotential are capable to induce a growth of eccentricity such that the perigee is lowered down to altitudes where drag becomes effective. The synergic effect of drag and other perturbations can be possibly exploited at these altitudes only for initial eccentricities higher than 0.1 6 . The most evident signatures in the maximum eccentricity maps are those at i = 63.4 • , 116.6 • , also known as critical inclinations (e.g., [27]), which corresponds to the conditionω = 0 (resonance 9 in Table 1).  middle, the pericenter altitude and on the right, the apocenter altitude. As it can be seen, the behaviour is different if the initial inclination corresponds exactly to the resonant value or not. Up to initial e = 0.1, for both inclinations the eccentricity does not experience a sufficient growth to lower the perigee in order to reenter. Indeed, for the case of an initial quasi-circular orbit (e = 0.001), at resonance the perigee lowers only by 70 km after 10 years and 177 km after 25 years, while for i = 63.5 • the decrease of the perigee is 58 km after 10 years and 115 km after 25 years.
At resonance we can observe that the characteristic period of the eccentricity evolution is clearly longer than in the neighborhood of the resonance. For example, for i = 63.4 • and e = 0.001, the eccentricity shows a period of 137 years, while for i = 63.5 • it reduces to 76 years. For higher eccentricities, such as e = 0.13 and e = 0.14, at i = 63.4 • the initial growth of eccentricity induced by the perturbation lowers the perigee down to an altitude where atmospheric drag becomes effective. Conversely, for i = 63.5 • the apogee starts to lower while the perigee is not low enough for drag to be effective in less than 200 years. Finally, for e = 0.15 the perigee is low enough that reentry is feasible at both initial inclinations thanks to the atmospheric drag.

Analysis in the frequency domain
The frequency components detected at each inclination for the initial orbits at a = 7978 km and a = 8578 km, assuming initial e = 0.001, with A/m = 0.012 m 2 /kg are shown in Figure 14, where each frequency signature Table 4: List of the other detected resonances due to lunisolar perturbations [19]: argument ψ j , values of the coefficients α, β, γ and corresponding index j.
corresponds to a filled square and the color bar refers to the relative amplitude found in the Fourier spectrum. The solid curves in the figure represent the resonant arguments ψ j , with j = 7, ..11, associated to solar gravitational and lunisolar perturbations, shown in Table 1; the dashed curves refer, instead, to fainter, but still detectable, signatures listed in Table 4. They correspond to the arguments ψ j with j = 18, ..20 associated to singly-averaged solar gravitational resonances, and to the argument ψ 21 associated to doubly-averaged lunisolar perturbations [19].
The main signature in both frequency charts is the one at i = 63.5 • associated to resonance 9, which corresponds also to the brightest corridor in the eccentricity contour maps of Figure 12. Concerning resonance 8, around i = 40 • , the contour maps showed that it is not expected to be relevant for e = 0.001, while it becomes more important for more eccentric orbits. Indeed, it is only partially detectable in the e = 0.001 frequency charts of Figure 14, while its role becomes more evident in the frequency charts of Figure 15, which correspond to the same initial orbits of Figure 14 but with e = 0.01. Comparing the frequency charts corresponding to the two values of eccentricity, we can notice also that resonances 7 , 8 , 11 and the higher order resonances shown in Table 4 are only partially detectable in the e = 0.001 frequency charts, while they are clearly recognizable for the e = 0.01 ones. In particular, the signature due to theψ 21 term is clearly visible in the a = 8578 km maximum eccentricity map of Figure 12 as the bright corridor at i = 69 • , and it appears also in the corresponding frequency chart. Figure 12 showed that the growth of eccentricity that can be reached thanks to high degree zonal harmonics and/or lunisolar perturbations, for the initial eccentricities considered, is, at most, one order of magnitude less than exploiting SRP in the case of an area augmentation device.
The most favourable case is found in proximity of resonance 9 (ω 0),   where ∆e max 0.02 can be achieved. As already noticed, the frequency analysis shown in Figure 14 confirms this finding for both altitudes: the main signature appears at i = 63.5 • , corresponding to the cusp of the |ψ 9 | curve. Figure 16 compares the behaviour of the numerical maximum eccentricity over propagation (cyan squares) and the amplitude found through the frequency analysis (blue squares) around i = 63.5 • for an initial orbit with a = 7978 km and e = 0.001. As for the case of model I, there is a very good match between the two quantities. We can notice that the growth of eccentricity at i = 63.5 • is mainly due to the perturbing effect of J 5 . Indeed, if we consider only a 3 × 3 geopotential instead of 5 × 5, the increment of eccentricity decreases from ∆e 5×5 = 0.017 to ∆e 3×3 = 0.002, while if only lunisolar perturbations and 2 × 2 geopotential are included in the dynamical model, the eccentricity does not experience any variation at this inclination. These results are shown in Figure 17, which displays the evolution of eccentricity for initial a = 7978 km and i = 63.5 • for three different models, all including drag and lunisolar perturbations: (i) 5 × 5 geopotential, (ii) 3 × 3 geopotential, (iii) 2 × 2 geopotential.
Although at high altitudes in LEO the growth of eccentricity induced by geopotential or lunisolar perturbations is not capable to drive the reentry, the variation in e can be, anyway, not negligible. Indeed, the perigee and apogee of the orbit can experience an oscillation which should be taken into account if we are dealing with issues as the stability of an operational orbit. This happens, for example, in the considered case of initial a = 7978 km and e = 0.001 and assuming a 5 × 5 geopotential as in model II-(i): the perigee undergoes a 76 years periodic evolution with a maximum oscillation of 130 km; after 10 years it experiences a variation of 55 km, while as much as 115 km after 25 years.

Conclusions
In this paper we studied the evolution of the eccentricity of a large set of orbits both in the time and frequency domains, deepening the work already presented by the authors in [3,4].
First, we considered the role of SRP in driving the dynamics for an object equipped with an area augmentation device. We found that, for quasi-circular orbits, SRP can be exploited, possibly in concurrence with the atmospheric drag, to lead the disposal within 25 years, but only if the initial orbital inclination is close enough to the resonant inclinations associated with the conditionψ =Ω±ω−λ S 0 (resonances 1 and 2). In the vicinity of the other zero-order resonances (indexed from 3 to 6), but also in correspondence of the first-order SRP resonances (indexed from 12 to 17), a growth of eccentricity due to SRP takes place in any case but over longer time scales, of the order of tens to hundreds of years. Although this variation of eccentricity cannot be exploited for disposal, it needs to be taken into account for operational purposes in the perspective of identifying long-term stable orbits within LEO.
Moreover, in [4] we presented a simplified theory to analytically evaluate the growth of eccentricity induced by the six main SRP resonances. Here, we showed that the assumption to consider the variation of eccentricity only as a function of the initial (e, i) state could be coarse and that, for given initial orbits, also the role of the variation of inclination over time should be considered, to give a coherent picture of the dynamics.
Then, we focused on the role of lunisolar perturbations and high degree zonal harmonics. In this case, the growth of eccentricity induced by the perturbations does not cause a lowering of the perigee leading to a reentry, in the case of quasi-circular orbits. In particular, we analysed the case of the well-known critical inclination, corresponding to the resonant conditioṅ ω 0, for an initial quasi-circular orbit at a = 7978 km. We verified that the computed growth of eccentricity of about 2 orders of magnitude after 40 years is mainly due to the J 5 perturbation, confirming the results found in [3].