Solar radiation pressure resonances in Low Earth Orbits

The aim of this work is to highlight the crucial role that orbital resonances associated with solar radiation pressure can have in Low Earth Orbit. We review the corresponding literature, and provide an analytical tool to estimate the maximum eccentricity which can be achieved for well-deﬁned initial conditions. We then compare the results obtained with the simpliﬁed model with the results obtained with a more comprehensive dynamical model. The analysis has important implications both from a theoretical point of view, because it shows that the role of some resonances was underestimated in the past, and also from a practical point of view in the perspective of passive deorbiting solutions for satellites at the end-of-life.


I N T RO D U C T I O N
It is known that the effect of the solar radiation pressure (SRP) on the orbital motion is a long-period variation in eccentricity e and inclination i, along with one in longitude of ascending node and argument of pericentre ω, and that the magnitude of this variation depends on the area-to-mass ratio of the body (e.g. Beutler 2005). If we assume a perturbing potential which is averaged over the orbital period of the body, the SRP effect can be subject to resonances, associated with a commensurability among the rate of precession of the ascending node, the one of the argument of pericentre and the mean apparent motion of the Sun, say n S . As well known, when a resonance occurs, the period of the variation can become so large to give rise to quasi-secular effects in the corresponding orbital element.
Motivated by the numerical results (Alessi et al. 2017a,b) obtained recently within the ReDSHIFT H2020 project , the aim of this work is to describe the role of SRP resonances in the Low Earth Orbit (LEO) region, showing that their contribution can provide relatively significant eccentricity variations.
In the past, the effect of the SRP for orbits around the Earth was considered mainly in the perspective of bodies characterized by a very high area-to-mass ratio. We can mention works focused on the behaviour of Geostationary Earth Orbits (GEO) (e.g. Valk, Lemaître & Anselmo 2007;Rosengren & Scheeres 2013;Casanova, Petit & Lemaître 2015), or on mission concepts aiming to exploit a solar sail to deorbiting non-operational objects from Medium Earth Orbits (MEO) (e.g. , or on the design of frozen orbits for small-size objects (e.g. Colombo Lantukh et al. 2015). , in particular, analysed SRP effects on high area-to-mass objects in order to design frozen orbits for a swarm of small spacecraft dedicated to Earth observation and telecommunication purposes. Starting from the work by Krivov et al. (1996), they studied the behaviour in eccentricity corresponding to˙ +ω − n S = 0, by computing the equilibrium points, and their stability, in a dynamical system including the oblateness of the Earth and SRP.
From a theoretical point of view, the first works on the resonances induced by SRP started from the analysis of the resonances induced by the solar gravitational attraction. As a matter of fact, the two disturbing potentials differ in the first order of the expansion and in the amplitude of the perturbation (Hughes 1977). If we focus on the resonances affecting the eccentricity evolution, Musen (1960) and Cook (1962) were the first to highlight the existence of six SRP resonances, and to show their location in the (i, a) plane. Cook (1962), in particular, noted that, contrary to lunisolar gravitational resonances, the SRP resonances are able to give a variation in eccentricity also in the case of circular orbits. Musen (1960) labelleḋ +ω − n S as the 'most interesting resonance', as done by all the following authors. In the case of the gravitational perturbation, this resonance is also known as evection resonance 1 (e.g. Brouwer & Clemence 1961;Touma & Widsom 1998;Frouard, Fouchard & Vienne 2010). Breiter (1999, basing his analysis on the work done by Cook (1962), noted that 'resonant solar perturbations can be much stronger than the lunar ones'. He also considered˙ +ω − n S as the dominant resonance for prograde orbits and˙ −ω + n S for retrograde orbits. In his work, the strength of a given resonance is determined by the magnitude of the libration region. As already mentioned, a fundamental contribution to the topic was given by Hughes (Hughes 1977(Hughes , 1980(Hughes , 1981, who tried to provide a more general treatment to the problem, by analysing also the effect due to high-order terms. 2 He stressed that SRP resonances give rise to variations not only in eccentricity, but also in inclination, longitude of the ascending node and argument of pericentre. He considered a resonance as dominant according to the magnitude of the corresponding amplitude, function of semimajor axis a and eccentricity. More recently, Celletti et al. (2017) recognized the occurrence of semisecular resonances, as the SRP ones, in the LEO region, but they stated that 'Since the semisecular resonances occur at specific altitudes, the width of such resonances is small, and since the air drag provokes a decay of the orbits on relatively short time scales, one expects that these resonances play a minor role in the long-term evolution of space debris'. We will show that this is partially true. Though the regions where SRP resonances are effective are indeed narrow, for high values of area-to-mass ratio, they represent the main mechanism to drive a satellite back to the Earth also when the LEO is sufficiently high to have a negligible interaction with the residual atmosphere, e.g. above an altitude h 1000 km. On the other hand, if we consider low values of area-to-mass ratio, the combined effect of the atmospheric drag and the SRP resonances can reduce significantly the lifetime of the object.
One of the fundamental issue to mitigate the space debris problem and control the growth of the debris population is to see if there exist natural perturbations which facilitate the reentry to the Earth or, when this is not possible, ensure a stable graveyard orbit in the long term. At altitudes when the atmospheric drag is not effective, the only way to lower the pericentre altitude is to obtain an eccentricity increase. For MEO, like those of the GNSS constellations, the gravitational lunisolar perturbations can help to this end too (Rosengren et al. 2015;Alessi et al. 2016). In LEO, apart from the atmospheric drag, the main natural effects that can be exploited are the lunisolar gravitational resonances, in particular the one corresponding to˙ + 2ω (Alessi et al. 2017a,b) and the SRP. In particular, as just mentioned, area-to-mass values representing small spacecraft equipped with a solar sail can experience a change as dramatic as to reenter to the Earth, even without the support of the atmospheric drag. These eccentricity variations due to SRP occur in correspondence of the˙ +ω − n S resonance, considered as dominant in the past, but not only. On the other hand, understanding how the eccentricity might change due to the SRP is crucial also to control the stability of a graveyard orbit.
As a final note, we acknowledge that looking for the values of area-to-mass ratio ensuring a reentry for a circular orbit, given initial semimajor axis and inclination, Colombo & de Bras de Fer (2016) found numerically the same reentry corridors represented by at least two of the six SRP resonances, that will be described here.

THE MODEL
Let us assume that the radiation coming from the Sun is directed normally to the surface of the spacecraft, that is, the so-called cannonball model. Moreover, let us assume the orbit of the spacecraft entirely in the sunlight, and the effect of Earth's albedo as negligible.
Neglecting the light aberration, the SRP is a conservative force. The corresponding disturbing potential, averaged over the orbital Table 1. Resonance ψ j = n 1 + n 2 ω + n 3 λ S in terms of n 1 , n 2 , n 3 .
motion of the spacecraft, can be written as (e.g. Krivov et al. 1996) where λ S is the longitude of the Sun measured on the ecliptic plane, = 23.439 • is the obliquity of the ecliptic, P the solar radiation pressure, C R the reflectivity coefficient and A/m the area-to-mass ratio, with n 1 = {0, 1}, n 2 = ±1, n 3 = ±1, according to j, following Table 1, and T 5 = sin 2 2 cos 2 i 2 , If we assume that the dynamics is driven by a single resonance, then we can write the variation in the orbital elements as due to only one, say j, of the terms in (1) at a time. By applying the Lagrange planetary equations, we get where n is the mean motion of the satellite and we have assumed that the only other effect exerting a perturbation on the spacecraft is the oblateness of the Earth. We havė where J 2 is the second zonal term of the geopotential and r ⊕ the equatorial radius of the Earth. The effect of lunisolar perturbations on˙ ,ω can be neglected, following, e.g. Milani, Nobili & Farinella (1987). Note that long-period and secular effects associated with the oblateness of the Earth affect only the behaviour of the longitude of the ascending node, of the argument of pericentre and of the mean anomaly. That is, the eccentricity does not change in the long term due to J 2 (see e.g. Roy 1982).

T H E E C C E N T R I C I T Y R E S O NA N C E S
From equation (4), we can recognize the following six resonances affecting the behaviour in eccentricitẏ In Fig. 1, we show their location as a function of the inclination and semimajor axis, for given values of e and A/m, assuming that only the oblateness of the Earth is responsible of a variation in and ω. This approximation can be considered valid for A/m = 0.012 m 2 kg −1 , which represents the average value of the intact objects orbiting in LEO. On the bottom panel, we show a close-up in the LEO region, defined here up to h = 3000 km, in order to account also for possible graveyard orbits. Note that we prefer to represent the resonances as a function of i, a instead of a, e, as done, e.g. in , because we focus on values of eccentricity which represent orbits of operational spacecraft in LEO. Realistic values cannot exceed the value e = 0.15, the majority of the population in LEO actually characterized by e < 0.02. In general, except in the cases where the spacecraft is orbiting in a region where the atmospheric drag is dominant, or at a given resonance affecting the eccentricity evolution, the eccentricity can be assumed to vary in a negligible way with respect to the initial value. Also, in case of resonance, the time-scale variation of the eccentricity is much slower than the typical period of the orbits considered.
In Fig. 2, we show the location of the same resonances for two different values of eccentricity, accounting also for the variation of and ω due to SRP. In this case, the location of the resonances depends also on , ω, λ S . In the figure, we compute the curves corresponding to ψ = 0 • , 90 • , 180 • for A/m = 1 m 2 kg −1 , which represents a feasible value, given the current level of technology, for a small spacecraft equipped with an SRP enhancing device (see Colombo, Rossi & Dalla Vedova 2017). For a given resonance, the middle curve corresponds to the case ψ = 90 • , which is equivalent to the case of Fig. 1. The displacement among the curves of a same resonance can be appreciated at higher LEO, for quasi-circular orbits for the first and the second resonance, following the convention in Table 1. Note, however, that the first equation in (4) states that ψ = 0 • and ψ = 180 • are stationary points for the eccentricity, and thus we expect to have a quasi-secular behaviour in eccentricity only at approaching the condition ψ = 90 • .
Increasing the eccentricity, the value of the semimajor axis at which the crossing appears also increases, e.g. for e = 0.1 the first crossing occurs at a ≈ r ⊕ + 2245 km, the second at a ≈ r ⊕ + 1220 km. From preliminary simulations, the dynamics at the overlapping does not manifest a chaotic behaviour, rather the two resonances concur −ω + n S = 0. For a given resonance, the middle curve corresponds to ψ = 90 • . If n 3 = 1 the right curve corresponds to ψ = 0 • , the left one to ψ = 180 • ; if n 3 = −1 vice versa. to a possible increase in eccentricity. In general, we can observe isolated points where the Lyapunov time is significantly lower than that in the neighbourhood. They do not affect the overall dynamics, but they can be, however, object of further dedicated studies.
To compute the width of the resonances considered, the same argument developed for lunisolar doubly averaged gravitational perturbations by Daquin et al. (2016) can be applied here. Let us assume that the Hamiltonian function of the system is where H K is the Keplerian part, H J 2 the one associated with the oblateness of the Earth, and H SRP = −R SRP . The curves defining the boundaries of a given resonance j at X * ≡ (e * , i * ) are defined by the condition and ∂ 2 H sec ∂X 2 as in Daquin et al. (2016). For the sake of clarity, 4 ∂ 2 H sec ∂X 2 = 3 2 J 2 r 2 ⊕ a 4 (1 − e 2 ) 5/2 n 2 2 2 − 15 cos 2 i + 10n 1 n 2 cos i − n 2 1 . The width computed in this way turns out to be very narrow. In Fig. 3, we provide an example for A/m = 1 m 2 kg −1 , with a closeup in the regions where the width is the largest. Note that they correspond to resonances #1 and #3.
With respect to a possible ranking of the six resonances in terms of the effect they might produce, the maximum eccentricity variation that can be achieved at a resonance j can be estimated by integrating the first equation in (4). This is, In Fig. 4, we show the corresponding values as a function of i, a for e = 0.01, = 0 • , ω = 0 • , λ S = 90.086 • and A/m = 1 m 2 kg −1 and A/m = 0.012 m 2 kg −1 , sampling the (i, a) plane at a step of i = 0.2 • and a = 20 km. In the former case, the colour bar representing the maximum eccentricity variation achievable is bounded to 0.3, which is the value required to reenter from the highest LEO considered. In the latter case, the limit is set to 0.1 just for the sake of clarity, that is, to be able to appreciate different behaviours, given the extreme narrowness of the resonances.
For A/m = 0.012 m 2 kg −1 , it should be noted that the variation that can be obtained is weak, and that, by definition of resonance, the time to get to the peak value is significantly long. Also, given that the range giving the required combination of inclination and semimajor axis is tapered, a possible exploitation of these corridors for an operational satellite would require an accurate manoeuvring capacity. A similar argument applies for resonance #6 for A/m = 1 m 2 kg −1 .
Note also that the same expression, equation (10), can be used to quantify the maximum possible change in eccentricity in the neighbourhood of a given resonance, eventually avoiding an unrealistic operational time. This application is faced in detail in . We remark here that, for the strongest resonances (#1, #2 and #3), the range of initial inclinations that can be targeted to achieve a considerable change in eccentricity for the high area-to-mass case is generally at most ±1 • with respect to the resonant value i * , for given initial semimajor axis and eccentricity.

N U M E R I C A L E V I D E N C E
The analysis described in the previous section was compared to the results of the numerical simulations performed under the ReD-SHIFT H2020 project (Alessi et al. 2017a,b;Rossi et al. 2017). The aim of the simulations was to map the whole LEO region in terms of initial orbital elements, considering two initial epochs and the same two values of area-to-mass ratio considered here, in order to look for natural deorbiting solutions that may facilitate the compliance with the international mitigation guidelines. A well-defined grid of initial conditions was propagated for 120 yr by means of the semi-analytical orbital propagator FOP (Anselmo et al. 1996;Rossi et al. 2009), including in the dynamical model the geopotential up  to degree and order 5, SRP, lunisolar perturbations and atmospheric drag below an altitude of 1500 km. 5 The grid in initial inclination, in particular, was set at steps of 2 • in the range [2 • : 120 • ].
In Fig. 5, we show the maximum eccentricity achieved in 120 yr, as a function of initial inclination and eccentricity, starting from three different values of semimajor axis (a = r ⊕ + 1200 km, a = r ⊕ + 1360 km, a = r ⊕ + 1580 km), for quasi-circular orbits, with a same initial configuration, and assuming A/m = 1 m 2 kg −1 . It is possible to appreciate that the eccentricity can experience a non-negligible variation, only in correspondence of well-defined inclination bands (the brighter ones). These inclinations are located in the neighbourhood of the resonant values, described in the previous section. In other words, within the limit of the grid adopted, only at resonant inclination values, the eccentricity can naturally increase.
It is also interesting to note that in the three examples the role of the dominant resonance, that is the one providing the greatest variation, changes. This is what we have called a resonance switching. In the top panel, this role is taken by resonance #2; in the middle by resonance #1; in the bottom by resonances #1 and #3.
In Fig. 6, we compare the maximum eccentricity computed with the numerical simulation with the maximum variation obtained using estimate (10), for the same three initial conditions as in Fig. 5, considering e = 0.01 as initial eccentricity. The analytical value depicted is the sum of the six variations forecast by equation (10) for each resonance. This choice turned out to approximate better the actual behaviour of the maximum eccentricity, especially in the inclination regions between two resonant values. Given that in the simulations carried out with FOP other perturbations played a role and that the atmospheric drag was included, we consider that there exists a reasonable consistency between the numerical model and the simplified estimation. The resonance switching can also find an explanation, looking to the relative height corresponding to the cyan points, in particular for the first two examples.
In general, for all the cases explored, it turned out that also resonance #4 can yield a significant eccentricity variation, for the high area-to-mass value. See Fig. 7 for two examples.
For A/m = 0.012 m 2 kg −1 , in Figs 8 and 9 we provide two examples of the eccentricity evolution and the consequent variation in the pericentre altitude by propagation with FOP, for altitudes where the atmospheric drag is not effective. They correspond to quasi-circular orbits, which are located at the resonance #2 and #1, respectively, at the initial epoch. We note not only the variation of hundreds of kilometres obtained, but also the long time required. In Fig. 10, we show an example when the SRP resonance aids the action of the atmospheric drag. Starting from i = 41.95 • , the reentry is achieved in about 55 yr; starting from i = 42 • , the lifetime is reduced by about 8 yr. This shows that if the satellite is located in the neighbourhood of an SRP resonance, 6 it could be worth to increase a little the eccentricity to clearing the region more rapidly. In this way, we can take advantage of both putting the spacecraft in a denser layer of the atmosphere and the push in eccentricity given by SRP, which can lower the pericentre altitude in the long term. An accurate modelling of the interaction between SRP and atmospheric drag will be faced in the future.

C O N C L U S I O N S
In this work, we have reviewed the main orbital resonances associated with SRP which affect the evolution in eccentricity, in order to see if they might be exploited to design advantageous deorbiting solutions for LEO. Starting from a well-known theory, we have    showed that at least four out of six resonances can be considered to increase the eccentricity of a small spacecraft equipped with a solar sail. The variation that can be achieved along a given resonance, that is, at well-defined values of semimajor axis, inclination and eccentricity, is such that a reentry to the Earth can be obtained, even from high altitudes. For typical values of operational satellites, SRP resonances can be considered, instead, in combination with the atmospheric drag, provided an accurate manoeuvring capability. The simplified model including only SRP and Earth oblateness has been validated against a more realistic dynamical model. The corresponding results turn out to be consistent. The novelty of the work regards the possible applications, but it is also theoretical, because we have shown the importance of resonances which were neglected in the past.
Future work will look into the role of the initial phasing provided by initial epoch, longitude of ascending node and argument of pericentre, which could ensure either an increase or a decrease in eccentricity. A detailed analysis on the characteristic frequencies for the eccentricity evolution in the resonant regions due to SRP is ongoing. That study will also include a comparison between the theoretical derivation explained here and the numerical results, without the contribution of the drag. Part of the results can be found in . Also, the inclination evolution was neglected in the analysis presented, and it will be evaluated in detail. From the point of view of a space operator, the model is still simplified in the sense that the orientation of a sail cannot be assumed as always directed normal to the solar radiation. This is another issue, which will deserve further research.