A frequency portrait of Low Earth Orbits
Giulia Schettino, Elisa Maria Alessi, Alessandro Rossi, Giovanni B. Valsecchi
AA frequency portrait of Low Earth Orbits
Giulia Schettino a, ∗ , Elisa Maria Alessi a , Alessandro Rossi a , Giovanni B.Valsecchi a,b a IFAC-CNR, Via Madonna del Piano 10, 50019 Sesto Fiorentino (FI) - Italy b IAPS-INAF, Via Fosso dei Cavalieri 100, 00133 Rome - Italy
Abstract
In this work we deepen and complement the analysis on the dynamics ofLow Earth Orbits (LEO), carried out by the authors within the H2020 ReD-SHIFT project, by characterising the evolution of the eccentricity of a largeset of orbits in terms of the main frequency components. Decomposing thequasi-periodic time series of eccentricity of a given orbit by means of a nu-merical computation of Fourier transform, we link each frequency signatureto the dynamical perturbation which originated it in order to build a fre-quency chart of the LEO region. We analyse and compare the effects on theeccentricity due to solar radiation pressure, lunisolar perturbations and highdegree zonal harmonics of the geopotential both in the time and frequencydomains. In particular, we identify the frequency signatures due to the dy-namical resonances found in LEO and we discuss the opportunity to exploitthe corresponding growth of eccentricity in order to outline decommissioningstrategies.
Keywords:
LEO, frequency analysis, SRP, lunisolar perturbations
1. 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 con-text, as part of the H2020 ReDSHIFT (Revolutionary Design of Spacecraft ∗ Corresponding author
Email addresses: [email protected] (Giulia Schettino), [email protected] (Elisa Maria Alessi), [email protected] (AlessandroRossi), [email protected] (Giovanni B. Valsecchi) pre-print September 14, 2018 a r X i v : . [ phy s i c s . s p ace - ph ] S e p hrough Holistic Integration of Future Technologies) project [1, 2], a deepanalysis to search for passive deorbiting solutions in LEO was carried out,by performing an accurate mapping of the phase space in order to identifystable and unstable regions. A detailed description of the results of the LEOcartography was presented by the authors in [3], while in [4] a general anal-ysis 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 theorbits, and the associated mechanisms, where dynamical perturbations caninduce a significant growth of the orbital eccentricity, in order to facilitatepassive disposal. Accordingly with our findings, we concluded that in thecase of a typical intact object in LEO, with an area-to-mass ratio of theorder of A/m = 10 − m / kg, perturbations as SRP, lunisolar effects andhigh degree zonal harmonics cannot ensure the reentry on their own butonly 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 alonecan drive the dynamics, if the initial inclination of the orbit is close enoughto a resonant inclination, given semi-major axis and eccentricity.In the present work, we make a deeper analysis of the role of resonanceswhich act in the LEO dynamics, by characterising the eccentricity of a setof orbits in terms of periodic components. Starting from the quasi-periodictime series of eccentricity, computed for a dense grid of initial conditions,we decompose the series in the main spectral components by means of a nu-merical computation of the Fourier transform. Then, we link each frequencycomponent with the dynamical perturbation responsible for that signature.In this way, we have an additional tool to explore the relative importanceof each given gravitational or non-gravitational perturbation in LEO as afunction of the initial orbital elements. The final goal of such analysis is tosupport the cartography in identifying the orbits where a significant growthof eccentricity, led by one or more perturbations, can assist the passive dis-posal of objects at their end-of-life. The same analysis can also serve toidentify the periodic drifts that operational orbits could experience.In the past, the study of the chaotic dynamics within the Solar Sys-tem led [5] to devise a method for a numerical estimation of the size of thechaotic zones, based on the variation in time of the main frequencies of thesystem. Since then, the algorithm for the frequency analysis was developedto study the stability of the orbits in many multi-dimensional conservative2ystems, in order to provide a global representation of the dynamics [6, 7].The Frequency Map Analysis algorithm (Numerical Analysis of FundamentalFrequencies - NAFF) is based on a refined and iterative numerical search fora quasi-periodic Fourier approximation of the solution of the system over afinite time span [8]. Considering, in particular, the issue of optimal design ofartificial satellite survey missions around a non-axisymmetric body, Noullezet al. [9] proposed an alternative method with respect to the standard Fouriertransform approach to characterise satellite orbits by computing the periodiccomponents in order to identify regular orbits, meant as orbits whose incli-nation and eccentricity do not vary significantly over a given time scale.Concerning in particular the LEO region, Celletti and Gale¸s [10] studied thedynamics of resonances in LEO with the aim of identifying the location ofequilibrium position and their stability. Within the common scope of definingsuitable post-mission disposal orbits, they studied analytically, by means of atoy-model, whether an object is located in a stable or chaotic region. In sucha way, the identification of stable orbits in LEO suggests the detection of pos-sible graveyard orbits. In this paper, we focus on the possibility of exploitingthe eccentricity growth induced by one or more dynamical perturbations atgiven orbits to facilitate the end-of-life reentry and we deepen this analysisby characterising the eccentricity evolution in terms of its main frequencycomponents. A comprehensive characterization of the dynamical evolutionof 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 dy-namical model adopted for the numerical propagation and we introduce themethod to identify the frequency signatures which characterise the eccentric-ity evolution of a set of LEO orbits. In Section 3 we outline the results of ouranalysis, comparing the results of numerical propagation in the time domainwith the findings of the frequency characterisation. Finally, in Section 4 wedraw some conclusions.
2. Dynamical model and methods
As mentioned before, within the scope of ReDSHIFT, we performed anextensive mapping of the LEO phase space by propagating more than 33illion orbits, as described in [11, 12, 3] , spanning from 500 km to 3000 kmof 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 analysisto the case of right ascension of the ascending node, Ω, and argument ofperigee, ω , both equal to 0 ◦ , with the initial epoch set to 21 June 2020. Theorbital propagation was carried out over a time span of 120 years by meansof the semi-analytical orbital propagator FOP (Fast Orbit Propagator, see[13, 14] for details), which accounts for the effects of 5 × A/m = 0 .
012 m / kg, selected as a reference value fortypical intact objects in LEO, and A/m = 1 m / kg, a representative valuefor a small satellite equipped with an area augmentation device, as a solarsail [15]. More details on the adopted model can be found in [3]. The resultsof the cartography can be displayed in contour maps showing the lifetimeor the maximum eccentricity over the propagation interval as a function ofthe initial inclination and eccentricity, for each initial semi-major axis. Alarge set of maps can be found on the ReDSHIFT website . In the followingSections, some examples will be provided. We are particularly interested in studying the time evolution of the eccen-tricity. Indeed, within the search for passive disposal solutions in LEO, theidentification of orbits which can experience a significant growth of eccen-tricity becomes crucial, since in this case the lowering of the orbital perigeehelps drag in being effective. Moreover, a variation in eccentricity causes analtitude 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 pertur-bations and high degree zonal harmonics cause long term periodic variationsin the evolution of eccentricity, which become quasi-secular in the vicinity of All the papers related to the project are available on the ReDSHIFT website athttp://redshift-h2020.eu/documents/. http://redshift-h2020.eu/results/leo . The oblateness of the Earth, J , does not affect the evolution of the eccentricity overlong term (e.g. [16]). able 1: List of the main resonances expected to be found in LEO: argument ψ j , valuesof the coefficients α, β, γ and corresponding index j . Resonances from j = 1 to j = 6are due to SRP; resonances 7 and 8 are singly averaged solar gravitational resonances;resonances from 9 to 11 are doubly averaged lunisolar resonances. Argument ψ j α β γ index j Ω + ω − λ S − − ω − λ S − − ω − λ S − ω + λ S ω + λ S − ω + λ S − ω − λ S − ω − λ S − ω ω ω ω . In particular, we can write the instantaneousvariation of e due to a given perturbation in the general form: dedt = T ( a, e, i ) sin ψ (Ω , ω, λ S ) , (1)where T is a coefficient which depends on ( a, e, i ) according to the givenperturbation and the argument ψ can be written in general terms as: ψ = α Ω + βω + γλ S , (2)where α, β, γ = 0 , ± , ± λ S is thelongitude of the Sun with respect to the ecliptic plane, set as λ S = 90 . ◦ at the starting epoch. A resonance occurs when the condition ˙ ψ (cid:39) ψ and the value of α, β, γ are highlighted, together with anindex ( j = 1 , ..
11) associated to each resonance. Resonances indexed from 1to 6 correspond to the condition˙ ψ = ¯ α ˙Ω ± ˙ ω ± ˙ λ S (cid:39) , (3)5 igure 1: Behaviour of | ˙ ψ | for each perturbing term j = 1 , ..
11 as a function of theinclination, for e = 0 .
001 and a = 7978 km, in the case A/m = 1 m /kg . with ¯ α = 0 ,
1, and are associated to the zero-order expansion of the SRPdisturbing function (e.g., [17, 18]). Resonances 7 and 8 are singly averagedsolar gravitational resonances (e.g., [19, 20]), while resonances from 9 to 11are associated to doubly averaged lunisolar gravitational perturbations (e.g.,[19]). The rate of Ω and ω can be found by applying the Lagrange planetaryequations and accounting in principle for both the effects of J and SRP,while the effect of lunisolar perturbations can be neglected (e.g., [21]). Theexplicit expressions have been given, for instance, in [4]. In practice, in [4]we have shown that, for an initial orbit with Ω = ω = 0 ◦ , at the assumedinitial epoch (which corresponds to λ S ≈ ◦ ), the rate of Ω and ω due toSRP vanishes.In Figure 1 we display the behaviour of | ˙ ψ j | for each perturbing term high-lighted 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 . a = 7978 km. The figure shows that curves associated to different pertur-bations may intersect and overlap creating a dense network of resonancesin the phase space. Thus, depending on the given inclination, it may behard to distinguish between the concurrent effect of different perturbationsand to link the dynamical effect to the perturbation which produces it. Toovercome this problem, we can take advantage of the fact that the adoptedorbital propagator is set up in such a way that each dynamical perturbationin the model can be individually turned on or off. Since we aim at identifying6he specific effect of a given perturbation on the eccentricity evolution and atcharacterising 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 ; • model II: SRP off; lunisolar perturbations and drag on; geopotential:5 × A/m objects, when only SRP and drag play a primaryrole in the evolution. In the case of
A/m = 1 m / kg, atmospheric drag iseffective in driving a reentry within 25 years for pericenter altitudes up to1050 km (see [3, 22]). Since this is a relatively high value, in order to focuson the effect due to SRP, we have decided to switch off the perturbation dueto the atmospheric drag.Model II, instead, is appropriate to study the effects led by lunisolarperturbations and high degree zonal harmonics: removing from the modelthe presence of SRP, we avoid the chance of mismodelling, since the resonantinclinations corresponding to lunisolar perturbations and geopotential can beclose 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 eccentricityevolution.
The starting point for the frequency characterisation is to process thediscrete eccentricity time series of a given initial orbit to obtain the discreteFourier transform through a standard Fast Fourier Transform (FFT) algo-rithm (e.g., [23]), based on the Cooley-Tukey algorithm [24]. The basic ideais to identify the frequency and the amplitude of the main spectral features inthe frequency series. The criterion we adopt is to account for any signaturewhose amplitude is, at least, 10 times stronger than the mean value of thespectrum in the surrounding area.A first issue to be considered concerns the time sampling ∆ t of the inputseries to be transformed. Indeed, the sampling frequency is f s = 1 / ∆ t and,from Nyquist theorem, it follows that f s / , the sampling ∆ t = 1 day, adoptedin [11, 12, 3], is fully reasonable. On the other side, a more critical issueinvolves the lowest detectable frequency by our analysis, which is limited by2 /T , where T is the duration of the time series. This means that with theadopted time span of 120 years, signatures with periodicity up to 60 yearswould be, in principle, identified. In practice, signatures due to perturbationswith 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, 600years, in order to catch unambiguously signatures with periodicity of sometens of years, as expected in the vicinity of a resonance.
3. Analysis of the numerical results
The general results of the LEO dynamical mapping was already exten-sively described in [3]. In the following, we present the results obtainedby 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 SRPin the case of the augmented
A/m ratio: we briefly recall the main findingsin terms of time evolution of the eccentricity, then we discuss the resultsof the characterisation in terms of frequency components. Next, we presentthe same analysis in the case of model II, focusing on the effects of lunisolarperturbations and high degree zonal harmonics.
We recall that the model accounts, in this case, only for the effect of SRPand J , while drag and lunisolar perturbations are turned off. We propagatethe orbits assuming A/m = 1 m / kg and we look for the inclinations where agrowth of eccentricity due to SRP occurs. Some illustrative results are shownin Figure 2: on the left we show the maximum eccentricity achieved over 600years of propagation as a function of the initial inclination and eccentricity,for initial a = 7978 km (top) and a = 8578 km (bottom), respectively. On theright panels we display the corresponding lifetime, in years. We recall thatthe atmospheric drag is effective up to 1050 km of altitude for the adopted We recall that moving close to a resonant orbit, the period of the perturbation actingon the eccentricity becomes gradually longer, up to quasi-secular if the orbital inclinationcorresponds exactly to a resonant condition. igure 2: Maximum eccentricity (left column) and lifetime over 600 years (right column;in the color bar: in years) as a function of the initial inclination at steps of ∆ i = 0 . ◦ and e at steps of ∆ e = 0 . A/m = 1 m / kg, for the initial orbits at a = 7978 km (top) and a = 8578 km (bottom), with Ω = ω = 0 ◦ and initial epoch 21 June2020. A/m ratio. Thus, we selected on purpose two reference values for the initialsemi-major axis which are significantly above the region where drag plays arole. If the effect of SRP is able to lower the perigee below 1050 km, then theremoval of the drag from the model allows to check if the chance to reenteror not can be ascribed solely to SRP.The lifetime panels show that, in the case an area augmentation deviceis available on-board, even for high altitude quasi-circular orbits, a reentrydriven by SRP alone is feasible for inclinations in the vicinity of 40 ◦ , whichcorresponds to the resonant condition:˙ ψ = ˙Ω + ˙ ω − λ S (cid:39) . (4)In the case of initial a = 7978 km, reentry can be achieved in about 7 yearsfor initial e ranging from 0.0001 to 0.009 thanks to SRP alone, for an initial9rbit at i = 39 . ◦ . In the case of initial a = 8578 km, SRP allows to reenterwithin 10 years at initial i = 37 . ◦ and in about 16 years for i = 38 ◦ . Theother resonances due to SRP, although not able to drive a reentry, cause,anyway, a remarkable growth in eccentricity, as can be seen from the leftpanels 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 SRPresonances as a function of i and a , we can identify the following resonancescorresponding to the bright inclination “corridors”: • ˙ ψ (cid:39) i = 40 ◦ (and i = 113 ◦ ); • ˙ ψ (cid:39) i = 80 ◦ ; • ˙ ψ (cid:39) ψ (cid:39) i = 58 ◦ and i = 54 ◦ , respectively, in the toppanel ( a = 7978 km), while they intersect around i = 56 ◦ at a = 8578km; • ˙ ψ (cid:39) ψ (cid:39)
0, both occurring in the vicinity i = 70 ◦ .Moreover, we can recognise other features at specific inclinations, ap-pearing as fainter, but still visible, signatures. They can be associated tohigher-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 interms of frequencies.For completeness, turning on the contribution due to the atmosphericdrag in the model, we find that the synergic effect of SRP and drag cansupport reentry also at different values of inclinations (resonances) but, typ-ically, only over long time scales. This is shown in Figure 3, in the case of aninitial orbit at a = 7978 km and e = 0 . The analysis of the maximum eccentricity maps (Figure 2 - left panels)shows that, in addition to the six resonances due to the zero-order expansionof the SRP disturbing function, other fainter signatures can be observed10 able 2: Resonant inclination i res and lifetime (in years) for each of the six main SRPresonances, in the case of initial a = 7978 km and e = 0 . Resonance i res Lifetime (yr)1 39 . ◦ . ◦ . ◦ . ◦
995 53 . ◦ Figure 3: Lifetime (in the color bar: in years) as a function of the initial inclination atsteps of ∆ i = 0 . ◦ and e at steps of ∆ e = 0 . A/m = 1 m / kg, for the initial orbit at a = 7978 km, with Ω = ω = 0 ◦ and initialepoch 21 June 2020. able 3: List of the first-order terms, expanding the SRP disturbing function up to first-order (e.g., [17]): argument ψ j , values of the coefficients α, β, γ and corresponding index j . Argument ψ j α β γ index jω − λ S − ω + 2 λ S ω − λ S − ω − ω − λ S − − ω + 2 λ S a = 7978 km and a = 8578 km inthe case of initial e = 0 . A/m = 1 m / kg, are shown in Figure 4.Each square in the plot represents a detected frequency component; the colorbar refers to the relative amplitude of the frequency signature , intended asthe corresponding intensity peak in the computed Fourier spectrum. Eachcoloured curve represents the behaviour of the argument | ˙ ψ j | as a functionof the inclination, with a cusp at the resonant inclination. As it can be seen,the detected signatures match almost exactly the theoretical curves. We alsopoint out that the amplitude of the signatures gradually grows approachinga resonant inclination. In particular, the effects of SRP first-order terms atgiven inclinations, which could be only partially inferred from the maximumeccentricity maps, can be clearly identified in the frequency chart.The signatures detected by means of the frequency analysis match thebright corridors detected in the maximum eccentricity maps in the left ofFigure 2. In particular, resonances 3 and 5 (see Table 1), which intersectfor a = 8578 km, can be individually identified for a = 7978 km both in The amplitude of each signature is normalised to the maximum detected amplitude,found in this case at the resonance ˙ ψ (cid:39) igure 4: Frequency signatures (filled squares) detected at each inclination for the initialorbit at a = 7978 km (top) and a = 8578 km (bottom), for initial e = 0 . A/m = 1m / kg. The | ˙ ψ j | curves are those associated to SRP resonances, shown in Tables 1 and3. The color bar refers to the relative amplitude of the frequency signature normalised tothe maximum detected amplitude. a = 8578 km shows a signature around i = 86 ◦ corresponding to thefirst-order ˙ ψ term, which does not appear for a = 7978 km neither in thecontour map nor in the frequency chart. Finally, the a = 7978 km chartshows a signature with singularity at i = 90 ◦ which can be associated to therate of Ω, appearing in the second-order expansion of the SRP disturbingfunction (see, e.g., [17]).From the lifetime maps in the right panels of Figure 2, we know that onlyin the case of resonance 1 SRP alone can drive a reentry. Nevertheless, themaximum eccentricity maps show that in the vicinity of a resonance a certaingrowth of eccentricity occurs anyway. Thus, in the perspective of designingpassive disposals and when dealing with operational issues, it is crucial toconsider the timescale over which the eccentricity variation takes place. Withthis in mind, assessing the change in eccentricity led by a perturbation with-out performing the numerical propagation, i.e., by characterising the LEOphase space in terms of frequencies, represents a very powerful tool.In [4], starting from Eq. (1) we showed that the maximum eccentricityvariation achievable due to the zero-order SRP resonance j for a given initial( a, e, i ) can be estimated as: ∆ e j = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) T j ( a, e, i )˙ ψ j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (5)On the other side, the amplitude associated to each detected frequency signa-ture in the Fourier transform gives an estimate of the eccentricity increment,as well. Both these values can be compared with the numerically computedmaximum eccentricity over 600 years: the three estimates are expected tocomply with each other.A general comparison for the initial orbit at a = 7978 km and e = 0 . i , we show the theoretical amplitude | T j / ˙ ψ j | with j = 1 , .. e max , achieved over the numerical propagation (redcircles) and the amplitude of frequency signatures detected by our analysis(filled squares; the color bar refers to the corresponding periodicity, i.e. theinverse of the detected frequency). A similar example for an orbit at a = 8578km is shown in [25]. The match is very good; moreover, we can observethat, as expected, the brighter squares, associated to signatures with longerperiodicity, are found only in the vicinity of resonances.14 igure 5: Theoretical amplitude | T j / ˙ ψ j | with j = 1 , .. a = 7978 km and e = 0 . Looking at the maximum variation in eccentricity achieved during prop-agation with FOP (red circles), some fainter features can be noticed at in-clinations different from those corresponding to the six main resonances.Comparing the inclination of these signatures with the resonant inclinationscorresponding to the arguments shown in Table 3, these fainter features canbe associated to the first-order terms in the expansion of the SRP disturbingfunction.In Figure 6, we show a detailed ( i, e ) zoom around the two main reso-nances found at this altitude: ˙ ψ corresponding to i ∼ ◦ and ˙ ψ in thevicinity of i ∼ ◦ . The maximum eccentricity displayed on the y − axis cor-responds to the eccentricity needed to lower the perigee down to 120 km, e km = 0 . | T / ˙ ψ | and | T / ˙ ψ | . Thisfurther confirms that the three quantities (theoretical amplitude, numericalmaximum eccentricity and amplitude of the frequency signature) provide thesame 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, adisagreement between the theory and the numerical propagation: accordingto the theory, the maximum eccentricity variation for initial i = 79 ◦ should15 igure 6: Comparison between theoretical amplitude | T j / ˙ ψ j | ( j = 1 on the top, j = 2on the bottom), maximum variation in eccentricity over propagation computed with FOP(left panels) and the frequency amplitudes detected by our analysis (right panels) as afunction of the inclination, for the initial orbit at a = 7978 km and e = 0 .
001 in the caseof model I. be sufficient to lead to reenter, while the ∆ e max computed with FOP turnsout to be lower than e km. The explanation for such a behaviour is thatduring the propagation also the inclination experiences a variation whichmoves the object away from the resonance, making the SRP perturbationless effective. In Figure 7 we show the evolution of e and i over 100 yearsfor the initial condition a = 7978 km, e = 0 . i = 79 ◦ . Both eccentricityand inclination show a periodicity of about 28 years but they are out ofphase: the eccentricity starts to grow led by the SRP perturbation; at thesame time, the inclination starts to decrease so that when the eccentricityreaches the maximum value e max = 0 .
14, the inclination is at its minimum, i min = 78 . ◦ , where, as can be inferred from Figure 6, the perturbation dueto the resonant term ψ is no longer effective in driving the reentry.This fact shows that the rate of i should be taken into account to providea full description of this case based on the dynamics. It is beyond the scope ofthis work to provide a full description on this scenario based on the dynamicalsystems theory, but we can provide a basic tool to obtain an a priori indicationon whether the orbit will exit from the resonance domain before achieving areentry.In Figure 7, in the panel showing the evolution of the inclination, it isalso displayed the behaviour predicted by the theory developed in [26] for16 igure 7: Eccentricity (top) and inclination (bottom) evolution over 100 years for initialcondition a = 7978 km, e = 0 . i = 79 ◦ in the case of model I. The inclination computedby propagation (blue line) is compared with the theoretical inclination (red circles) derivedfrom Eq. (7). lunisolar gravitational resonances, which can be applied also in the case ofSRP, as shown in [4]. In particular, it is demonstrated that there exists anintegral of motion, corresponding to( β cos i − α ) (cid:112) µa (1 − e ) = constant , (6)where α, β are as defined in Eq. (2). In other words, assuming that the mo-tion of the spacecraft is governed only by the Earth’s monopole, the Earth’soblateness and the solar radiation pressure, at any instant we can recover theinclination value from i = ± arccos (cid:32) constant β (cid:112) µa (1 − e ) + αβ (cid:33) , (7)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 theeccentricity and inclination evolution computed by propagation assumingmodel I (blue curve) with the behaviour obtained by assuming the completedynamical model (red curve), which includes all the perturbations providedby FOP. The initial orbit is the same as in Figure 7. We can observe that17 igure 8: Comparison of the eccentricity (left) and inclination (right) evolution over 30years, computed by propagation assuming model I (blue curve) and including all theperturbations provided by FOP (red curve). The initial orbit is for both cases: a = 7978km, e = 0 . i = 79 ◦ .Figure 9: Predicted inclination variation as a function of the initial inclination, assumingmodel I, a = 7978 km, e = 0 . 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 byEq. (7), by assuming a maximum variation in eccentricity as in Eq. (5), forresonances 1 and 2. We can notice that in the first case, when we consideran initial inclination in the resonance domain, the variation is not relevantif compared with the curves in the top panel of Figure 6). In the secondcase, the variation is instead important, of about 1 ◦ and moves the dynamicstowards the edges of the interval where the resonance is effective (comparewith the curves in the bottom panel of Figure 6).The above discussion showed that the assumption that T j is a function18 igure 10: Evolution of ∆ e = T / | ˙ ψ | over 60 years, computed by means of Eq. (5) onthe e and i values obtained by propagation with FOP in case of model I, for an initial a = 7978 km. of the initial values of eccentricity and inclination may provide a misleadinginformation. Figure 10 shows the evolution of ∆ e = T ( a, e, i ) / | ˙ ψ | , accord-ing to Eq. (5), assuming the values of eccentricity and inclination computedat each given time by propagation with FOP, in case of model I, for initial a = 7978 km, e = 0 . i = 79 ◦ . The y − axis upper limit corresponds to aperigee altitude of 120 km. As it can be seen, for the initial value of e and i ,the growth of eccentricity ∆ e is such that the reentry driven by resonance2 is feasible (the curve is not visible in the figure because it is higher thanthe eccentricity required to reentry). On the contrary, after only 5 years, theinclination has moved from its initial value (compare with Figure 7) enoughthat the corresponding growth in eccentricity due to resonance 2 alone is nomore capable to assure the reentry.Finally, similarly to Figure 6, the comparison between theoretical ampli-tude, maximum variation in eccentricity computed with FOP and amplitudeof the frequency signatures for a = 7978 km and e = 0 .
001 in the cases ofresonances 3 , , , Model II is particularly suitable to study the perturbation on eccentricitydue to lunisolar effects and high-degree terms in geopotential, since SRP hasbeen removed in this case. The effective area-to-mass ratio of the object does19 igure 11: Comparison between theoretical amplitude | T j / ˙ ψ j | ( j = 3 , , , 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 previousmodel, thus we assume
A/m = 0 .
012 m / kg for simulations.In analogy to the left panels of Figure 2, Figure 12 shows the maximumeccentricity as a function of the initial inclination and eccentricity for anorbit at a = 7978 km (left) and a = 8578 km (right), respectively. In thiscase, we do not show the corresponding lifetime maps: at these altitudes andfor quasi-circular orbits the maps would result blank since neither lunisolarperturbations nor high-degree terms of geopotential are capable to inducea growth of eccentricity such that the perigee is lowered down to altitudeswhere drag becomes effective. The synergic effect of drag and other perturba-tions can be possibly exploited at these altitudes only for initial eccentricitieshigher than 0.1 . The most evident signatures in the maximum eccentricitymaps are those at i = 63 . ◦ , . ◦ , also known as critical inclinations (e.g.,[27]), which corresponds to the condition ˙ ω = 0 (resonance 9 in Table 1). Contour maps similar to Figure 12 including eccentricities up to 0.28 can be found onthe project website. igure 12: Maximum eccentricity as a function of the initial inclination at steps of ∆ i =0 . ◦ and e at steps of ∆ e = 0 . A/m = 0 .
012 m / kg, for theinitial orbits at a = 7978 km (left) and a = 8578 km (right), with Ω = ω = 0 ◦ and initialepoch 21 June 2020. Figure 13 depicts the time evolution of different orbits with initial a =7978 km, considering two different initial inclinations: i = 63 . ◦ (top), whichcorresponds exactly to the resonant inclination for the condition ˙ ω = 0, and i = 63 . ◦ (bottom), i.e., only 0 . ◦ degrees next to the resonant value. Theinitial eccentricity varies from 0 .
001 to 0 .
15: on the left, we show the evolutionof eccentricity over 200 years, in the middle, the pericenter altitude and onthe right, the apocenter altitude. As it can be seen, the behaviour is differentif the initial inclination corresponds exactly to the resonant value or not. Upto initial e = 0 .
1, for both inclinations the eccentricity does not experience asufficient growth to lower the perigee in order to reenter. Indeed, for the caseof an initial quasi-circular orbit ( e = 0 . i = 63 . ◦ 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 eccen-tricity evolution is clearly longer than in the neighborhood of the resonance.For example, for i = 63 . ◦ and e = 0 . i = 63 . ◦ it reduces to 76 years. For higher eccentricities,such as e = 0 .
13 and e = 0 .
14, at i = 63 . ◦ the initial growth of eccentricityinduced by the perturbation lowers the perigee down to an altitude whereatmospheric drag becomes effective. Conversely, for i = 63 . ◦ the apogeestarts to lower while the perigee is not low enough for drag to be effectivein less than 200 years. Finally, for e = 0 .
15 the perigee is low enough that21 igure 13: Time evolution of eccentricity (left), perigee altitude (middle) and apogeealtitude (right) over 200 years of propagation with FOP, for initial a = 7978 kmand i = 63 . ◦ (top), i = 63 . ◦ (bottom), for 7 different initial eccentricities: e =0 . , . , . , . , . , . , .
15 in the case of model II, assuming as initial epoch21 June 2020. reentry is feasible at both initial inclinations thanks to the atmospheric drag.Looking at Figure 12, other fainter signatures at given inclinations canbe recognised: • at i (cid:39) ◦ , ◦ , in the a = 7978 km panel, corresponding to the well-known evection resonance (e.g., [28]) ˙ ψ = 2( ˙Ω + ˙ ω − ˙ λ S ) (cid:39) • at i (cid:39) ◦ , visible in the a = 8578 km panel, corresponding to thecondition ˙ ψ = ˙Ω + 2 ˙ ω (cid:39) • at i (cid:39) ◦ , clearly recognisable at a = 8578 km while distinguishableonly for very low eccentricities at a = 7978 km, which corresponds tothe resonant condition ˙Ω − ω (cid:39)
0, as will be discussed in Sect. 3.2.2.
The frequency components detected at each inclination for the initialorbits at a = 7978 km and a = 8578 km, assuming initial e = 0 . A/m = 0 .
012 m / kg are shown in Figure 14, where each frequency signature22 able 4: List of the other detected resonances due to lunisolar perturbations [19]: argument ψ j , values of the coefficients α, β, γ and corresponding index j . Argument ψ j α β γ index j Ω + 2 ω + 2 λ S − ω − λ S − − − ω − λ S − − − ω − ψ 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 de-tectable, 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 ψ associated to doubly-averaged lunisolar perturba-tions [19].The main signature in both frequency charts is the one at i = 63 . ◦ associated to resonance 9, which corresponds also to the brightest corridor inthe eccentricity contour maps of Figure 12. Concerning resonance 8, around i = 40 ◦ , the contour maps showed that it is not expected to be relevantfor e = 0 . e = 0 .
001 frequency charts ofFigure 14, while its role becomes more evident in the frequency charts ofFigure 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 ofeccentricity, we can notice also that resonances 7 , ,
11 and the higher orderresonances shown in Table 4 are only partially detectable in the e = 0 . e = 0 .
01 ones. Inparticular, the signature due to the ˙ ψ term is clearly visible in the a = 8578km 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 reachedthanks to high degree zonal harmonics and/or lunisolar perturbations, forthe initial eccentricities considered, is, at most, one order of magnitude lessthan exploiting SRP in the case of an area augmentation device.The most favourable case is found in proximity of resonance 9 ( ˙ ω (cid:39) igure 14: Frequency signatures (filled squares) detected at each inclination for the initialorbits at a = 7978 km (top) and a = 8578 km (bottom), assuming initial e = 0 . A/m = 0 .
012 m / kg. The | ˙ ψ j | curves are those associated to lunisolar resonances, shownin Tables 1 and 4. The color bar refers to the relative amplitude of the frequency signaturenormalised to the maximum detected amplitude. igure 15: Frequency signatures (filled squares) detected at each inclination for the initialorbits at a = 7978 km (top) and a = 8578 km (bottom), assuming initial e = 0 .
01 and
A/m = 0 .
012 m / kg. The | ˙ ψ j | curves are those associated to lunisolar resonances, shownin Tables 1 and 4. The color bar refers to the relative amplitude of the frequency signaturenormalised to the maximum detected amplitude. igure 16: Comparison between the maximum variation in eccentricity over propagationcomputed with FOP (cyan squares) and amplitude of the frequency signatures detectedby our analysis (blue squares) in the case of resonance 9, for the initial orbit at a = 7978km and e = 0 . where ∆ e max (cid:39) .
02 can be achieved. As already noticed, the frequencyanalysis shown in Figure 14 confirms this finding for both altitudes: themain signature appears at i = 63 . ◦ , corresponding to the cusp of the | ˙ ψ | curve. Figure 16 compares the behaviour of the numerical maximum eccen-tricity over propagation (cyan squares) and the amplitude found through thefrequency analysis (blue squares) around i = 63 . ◦ for an initial orbit with a = 7978 km and e = 0 . i = 63 . ◦ is mainly due to the perturbing effect of J . Indeed,if we consider only a 3 × ×
5, the increment ofeccentricity decreases from ∆ e × = 0 .
017 to ∆ e × = 0 . × a = 7978 km and i = 63 . ◦ for three different models, allincluding drag and lunisolar perturbations: (i) 5 × × × e can be, anyway, not negligible. Indeed, the perigee andapogee of the orbit can experience an oscillation which should be taken intoaccount 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 × igure 17: Evolution of e for initial a = 7978 km and i = 63 . ◦ for three different models,all including drag and lunisolar perturbations: (I) 5 × × × undergoes a 76 years periodic evolution with a maximum oscillation of 130km; after 10 years it experiences a variation of 55 km, while as much as 115km after 25 years.
4. Conclusions
In this paper we studied the evolution of the eccentricity of a large set oforbits both in the time and frequency domains, deepening the work alreadypresented by the authors in [3, 4].First, we considered the role of SRP in driving the dynamics for an objectequipped with an area augmentation device. We found that, for quasi-circularorbits, SRP can be exploited, possibly in concurrence with the atmosphericdrag, to lead the disposal within 25 years, but only if the initial orbitalinclination is close enough to the resonant inclinations associated with thecondition ˙ ψ = ˙Ω ± ˙ ω − ˙ λ S (cid:39) e, i ) state could be coarse and that, for giveninitial orbits, also the role of the variation of inclination over time should beconsidered, to give a coherent picture of the dynamics.Then, we focused on the role of lunisolar perturbations and high degreezonal harmonics. In this case, the growth of eccentricity induced by theperturbations 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 ofthe well-known critical inclination, corresponding to the resonant condition˙ ω (cid:39)
0, for an initial quasi-circular orbit at a = 7978 km. We verified thatthe computed growth of eccentricity of about 2 orders of magnitude after 40years is mainly due to the J perturbation, confirming the results found in[3]. Acknowledgements.