Dynamical taxonomy of the coupled solar radiation pressure and oblateness problem and analytical deorbiting configurations

Recent works demonstrated that the dynamics caused by the planetary oblateness coupled with the solar radiation pressure can be described through a model based on singly-averaged equations of motion. The coupled perturbations affect the evolution of the eccentricity, inclination and orientation of the orbit with respect to the Sun--Earth line. Resonant interactions lead to non-trivial orbital evolution that can be exploited in mission design. Moreover, the dynamics in the vicinity of each resonance can be analytically described by a resonant model that provides the location of the central and hyperbolic invariant manifolds which drive the phase space evolution. The classical tools of the dynamical systems theory can be applied to perform a preliminary mission analysis for practical applications. On this basis, in this work we provide a detailed derivation of the resonant dynamics, also in non-singular variables, and discuss its properties, by studying the main bifurcation phenomena associated to each resonance. Last, the analytical model will provide a simple analytical expression to obtain the area-to-mass ratio required for a satellite to deorbit from a given altitude in a feasible timescale.


Introduction
The effect of the Solar Radiation Pressure (SRP) on Earth satellites was recognised since the first space flights. The orbital evolution of Vanguard I [16] and the Echo balloons [20] was found to be significantly influenced by SRP, and singly-averaged equations were used to study the motion [15]. Treated as a perturbation problem, the singly-averaged contribution of SRP is integrable and analytical solutions can be obtained [14,17,19]. Nevertheless, when coupled with the effect of Earth's oblateness J 2 the system becomes a 2.5 Degrees-of-Freedom (DoF). An analytical insight can be recovered when treating locally the semi-secular SRP resonances. Namely, the singly-averaged perturbing function can be decomposed in six distinct terms [8], each of them dominating in a particular range of orbital elements. The dynamics arising by combining each of the harmonics with the secular evolution due to J 2 can be reduced to a 1 DoF resonant model [10,12,2].
The derivation of the six different resonant models in the three-dimensional case and their effect on the long-term evolution of resident space objects has been recently discussed in the literature [2,1]. In this work, we re-derive the resonant models in the Hamiltonian framework, providing also a non-singular representation of the resonant dynamics. We exploit the information of the analytical model in [1] to obtain further insight in the resonant structures. We employ the equations for computing the equilibria for each resonance and compute the number and their stability for each set of the dynamical parameters of the system. This allows us to construct bifurcations diagrams which give the main transitions in the phase space. Particular focus is given to the structure of the phase space about the vicinity of each SRP resonance. The effect of the engineering parameter, the spacecraft area-to-mass ratio, is also discussed.
Using the phase space portraits, useful information for mission design concepts is retrieved via the tools of dynamical system theory. In particular, a rigorous procedure to obtain deorbiting conditions along the resonances is presented. As for the planar case [12], we show how the deorbiting can occur.
The paper is organised in the following way: in Sec. 2 the basic force model is recovered in its singly-averaged formulation, in Sec. 3 we discuss the analytical description in the vicinity of SRP resonances, in Sec. 4 we provide a phase space analysis based on the main bifurcations associated to each resonance, in Sec. 5 we demonstrate the use of the models to analytically obtain feasible deorbiting configurations, and in Sec. 6 we present our conclusions.

Model derivation
Let us assume that a spacecraft moves under the effect of a planet's gravitational monopole, the planetary oblateness and the solar radiation pressure. Moreover, we assume that the sun-rays are always perpendicular to the surface of the satellite (cannonball model), that the effect of the planetary albedo is negligible, that the solar flux is constant at 1 AU and that the satellite is entirely in sunlight. Under these assumptions, SRP is modelled as a constant force in the direction of the Earth-Sun line, and it can be derived from a potential function. The dynamics of the satellite in a geocentric equatorial inertial frame can be modelled by the Hamiltonian where µ is the gravitational parameter of the Earth, and r, v are the geocentric distance and velocity of the satellite. The Earth's oblateness effect is modelled as where C J2 = µR 2 ⊕ J 2 with J 2 the oblateness parameter and φ the geographic latitude of the satellite. The sine of the latitude is expressed in terms of the orbital elements of the satellite via sin φ = z r , where i is the inclination, Ω the Right Ascension of the Ascending Node (RAAN), λ = ω + f the longitude, ω the argument of perigee and f the true anomaly of the satellite. The rotation matrices R 1 (u), R 3 (u) are The solar radiation pressure contribution is given by where C SRP = P c R A m , P being the SRP constant at 1 AU , c R the reflectivity coefficient and A m the area-to-mass ratio of the satellite. X is the coordinate of the satellite in an Earth-centred system with the X-axis pointed towards the Sun (see Fig. 1): in terms of the orbital elements of the satellite it reads with λ the ecliptic longitude of the Sun and ε the obliquity of the ecliptic. Both H J2 and H SRP are considerable smaller than H kep and can be treated as perturbations of the two-body problem where H per = H J2 + H SRP . Employing a Hori-Deprit approach [7,6], the homological equation up to the first order of the normalisation process reads: where n is the mean motion of the satellite, W 1 is the generating function and is the first-order normalized Hamiltonian. The integrals in Eq. (2) can be carried out in closed form for both the J 2 and SRP contributions. In the case of J 2 we use the differential relationship dM = r 2 a √ 1−e 2 df along with r = a 1−e 2 1+e cos f to obtain the classical result where c i = cos i. For SRP we express the integral with respect to the eccentric anomaly E using the relations r sin f = a √ 1 − e 2 sin E, r cos f = a(cos E − e), r = a(1 − e cos E) and dM = r a dE.
The averaged model is then [10,12,2] with j T j cos ψ j (n 1 , n 2 , n 3 ) where c i , s i , c ε , s ε are the cosine and sine of the inclination i and the obliquity of the ecliptic ε, respectively. The model consists of 6 distinct harmonics describing the semi-secular evolution of the system, all of them containing the ecliptic longitude of the Sun λ . A second averaging of over Sun's mean motion results in a null Hamiltonian, that is, the SRP perturbation does not give rise to secular effects. Moreover, we should mention that the Hamiltonian Eq. (4) is integrable, and an analytical solution is obtained using an ecliptic rotating frame [14].
The singly-averaged model of the coupled Earth's oblateness and solar radiation pressure effects reads The Hamiltonian can be expressed in terms of the canonical Delaunay elements (L, G, H, l, g, h)such that to obtainH where for the coefficients T j we use the relationships Due to the averaging process, Eq. (7) does not depend on the mean anomaly M = l and thus the Delaunay action L, as well as the semi-major axis, are constant. The system has two Degrees-of-Freedom (DoF) and an explicit time dependence through λ (t) = λ ,0 +n s t, where n s is the frequency corresponding to Earth's orbital period of 1 year. Considering an extended phase-space, a dummy action I s with frequency n s is added to the Hamiltonian which yields a three DoF autonomous system.

Resonances
We are now interested in studying the resonance phenomena occurring in the coupled system. A resonance is occurring when the following commensurability holds: whereġ = ∂H/∂G,ḣ = ∂H/∂H andλ = ∂H/∂I s . The integer coefficients n 1 , n 2 , n 3 can take the values of the 6 combinations associated to the harmonics appearing in Eq. (5). As a first approximation, the angular frequencies can be computed by taking into account only the H J2 part of the Hamiltonian, namely,ġ Substituting Eq. (9) in Eq. (8), we can derive an approximation of the loci in action-space of the center of each resonance [5,2]. In the vicinity of each resonance, the associate harmonic is slow and is dominating the dynamics.
Assuming an isolated resonance approximation, the system can be described locally from the Hamiltonian T j cos ψ j + n s I s .
For each resonance ψ j (with j = 1 . . . 6) we introduce a resonant set of variables (Ψ, ψ) via a unimodular transformation of the Delaunay elements which brings the Hamiltonian to a 1-DoF form of where the part associated with J 2 is and the one due to SRP is while the coefficients T j are expressed in terms of the new variables using the equations c i = n1n2Ψ +Π n 2 2 Ψ and s i = 1 − c 2 i . Due to the resonant transformation both π and κ are ignorable, therefore, Π and K are constants and the term n s K can be dropped from the Hamiltonian. The action variable Π is a resonant integral of the system; its value is dictated from the initial conditions and remains constant during the orbital evolution. Expressed in orbital elements represents coupled oscillations in the eccentricity and inclination of the orbit, in a similar fashion to the Lidov-Kozai constant in the third-body gravitational perturbations [11,9]. Note that Π corresponds to Λ in [1]. The formulation provided here is equivalent to the one given in the past paper, but n 1 and n 2 are exchanged. For a satellite close to a resonance ψ j with mean initial elements (a 0 , e 0 , i 0 , ω 0 , Ω 0 , λ ,0 ), the orbit evolution in the (Ψ, ψ) plane, or equivalently in the (e, ψ) plane if we substitute Ψ = µa(1 − e 2 )/n 1 , is given from the contour line of the implicit equation: with L = L(a 0 ), Π = Π(a 0 , e 0 , i 0 ). Notice that Π depends on the resonance j through n 1 and n 2 . The equations of the motion related the resonant models H ψj are and the associated equilibria are given from the equations Those stationary solutions of the resonant model represent periodic orbits of the full equations of motion. Their stability is determined from the Hessian eigenvalues of the Hessian of the Hamiltonian H ψj computed at each equilibrium point. In many applications the equations of motion with respect to the eccentricity are of interest, and thus, by applying the chain rule for the derivatives and using the relations Ψ = G/n 1 and de/dG = − √ 1−e 2 e √ µa we can derive the semi-canonical form: whereġ J2 ,ḣ J2 are given in Eq. (9) anḋ that is the model developed in [1]. In all the above formulation the angles g, h are not well defined when the eccentricity and/or the inclination are zero. To alleviate this burden we use the modified Delaunay variables in Eq. (10) while a unimodular transformation allows to derive a new set of resonant variables (Σ, σ), similar to the ones introduced in [4] for the case of lunisolar resonances, namely, where k 1 , k 2 , k 3 are integers appearing in combinations associated to each harmonic in Eq. (4) After a Taylor expansion in the resonant action Σ and dropping the constant terms, the resonant Hamiltonian up to a truncation order N takes the form where c n , d n are constant coefficients depending on the dynamical parameters (L, Φ), the particular resonance (k 1 , k 2 ) and the physical parameters (ε, µ). We notice that up to second order (N = 2) the resonant model is similar to the Second Fundamental Model for first-order in the eccentricity resonances (SF M I ). However, in practice a truncation order of N ≥ 4 is required to retrieve the qualitative features of the phase space, while higher orders are used for more accurate computations. Actually, the phase-space exploration can be done with the non-truncated resonant Hamiltonian and we retain the complete dynamical portrait of the resonance, a fact that highlights the power of the closed-form averaging process. Finally, it is common to express the resonant Hamiltonian in terms of Poincaré variables

Bifurcation analysis
From the equations written above it is possible to compute the equilibrium points associated with the dynamical system together with the corresponding stability. The isolated resonance system has two constants of motion (apart from H): the semi-major is constant due to the averaging and the second integral Π from Eq. (14). The values of Π can be labelled based on the inclination of the corresponding circular orbit Π ≡ Π circ = n 1 cos i circ √ µa − n 2 . Given a set of values of the dynamical parameters (a, i circ ) and the engineering parameter A/m, a phase space is uniquely identified. In this phase space the number   Fig. 2. The corresponding phase space is depicted in Fig. 3.
of equilibrium points and their stability can be defined. By varying the parameters and tracking the structural changes in the system, a bifurcation diagram can be computed. Elliptic and saddle points appear and disappear based on the classical bifurcation theory for 1-DOF systems.
For the two harmonics that dominate in the low Earth orbit region, namely, ψ 1 and ψ 2 [3], in Fig. 2 Fig. 4 the possible phase space structures are depicted in the range of semi-major axis a ∈ [7000 : 9400] km. The black curves define the boundaries of the regions characterized by a different number of equilibrium points and corresponding stability. In red, the curves that would be obtained by considering only the oblateness effect for the rate of precession of ω and Ω are shown, that is, only Eq. (9) instead of Eq. (17). For the resonant harmonic ψ 1 in Fig. 2 we present two bifurcation diagrams for A/m = 1 m 2 /kg (left) and A/m = 5 m 2 /kg (right). There are 5 distinct dynamical regions labelled with Latin numbers, defined in Tab. 1, the corresponding phase spaces being shown in Fig. 3. As just mentioned, the red line corresponds to the resonance location based on the J 2 rates. We observe that for the low A/m = 1 m 2 /kg value, this line almost coincides with the transition from region I to region III and from I to II. From Tab. 1 it is easy to see that the following bifurcations take place: ψ = 0:

and in
transition from region I to region III: saddle-node; -transition from region V to region II: saddle-node; ψ = π: transition from region III to region V: saddle-node; -transition from region I to region II: saddle-node. Fig. 4 The bifurcation diagram for resonance ψ 2 assuming A/m = 1 m 2 /kg. The total number of equilibria, reported in the colorbar, is presented in the dynamical parameter space (a, i circ ). The red line corresponds to the resonant condition computed from solely the J 2 rate. Inclination stands for inclination of the circular orbit. For the type of equilibria in each region see Tab. 2.

I II
Region total ψ 2 = 0 ψ 2 = π I 1 -1 s II 3 1 s & 1 u 1 s Table 2 Number of equilibria and their stability (s: stable, u: unstable) for the resonance ψ 2 , corresponding to the two regions of Fig. 4.
The effect of an increased A/m on the bifurcation diagram is shown in the right panel of Fig. 2. In this case, the exact position of the bifurcation is no longer Fig. 5 The bifurcation diagram for resonance ψ 1 (top), ψ 3 (middle) and ψ 4 (bottom). The red curves correspond to the bifurcation that is obtained considering only the oblateness effect inω andΩ.
predicted from the J 2 only rates as significant contributions are added to thė ω andΩ rates due to SRP. The number of regions and their structure remains the same, however, the boarders between the different regions considerably alter. Specifically, for a fixed i circ the bifurcation for region I to regions II and III seem to happen at higher altitudes. This is in agreement to previous findings on the planar case (i circ = 0) of ψ 1 [12]. A similar analysis is now performed for the resonance ψ 2 . Assuming the same range semi-major axis as in ψ 1 case, we expect this harmonic to dominate in a range of inclinations from i ∈ [76 : 84]. The bifurcation diagram is shown in Fig. 4 and the equilibria with their associated stability are given in Tab. 2. The situation in this case appears to be more straightforward. There are two distinct regions: in region I there exists only 1 stable equilibrium at ψ 2 = π while in region II two more equilibria, 1 stable and unstable, appear at ψ 2 = 0 after another saddle-node bifurcation. The position of the bifurcation for A/m = 1 m 2 /kg is approximately estimated from the J 2 resonant condition (red curve in Fig. 4). Figure 5 extends the analysis by displaying the bifurcation diagrams for resonances ψ 1 , ψ 3 and ψ 4 , assuming A/m = 1 m 2 /kg and a range of semimajor axis a ∈ [7000 : 20000]. As also found in [1], the maximum number of equilibrium points is 5 also for ψ 4 , and it is interesting to note that for ψ 1 it can occur another type of bifurcation apart from the saddle-node: the transcritical one beyond a ∼ 17000 km and i circ ∈ [130 • : 150 • ]. The bifurcation diagram for ψ 5 and ψ 2 are symmetrical with respect to i circ = 90 • , the same holds for ψ 6 and ψ 1 .

Deorbiting configuration
The phase space analysis can be exploited to compute the initial conditions, in terms of orbital elements and area-to-mass ratio, that can lead to an atmospheric reentry. Since the semi-major axis does not change under the dynamics considered, a natural deorbiting can take place only if the eccentricity increases as much as to attain the critical value e cr = 1 − r ⊕ /a.
We consider the case where this condition occurs at either ψ = 0 or ψ = π, meaning that the area-to-mass ratio required is the minimum. Moreover, mindful of the phase space behavior depicted in the previous section (see Fig. 3), the steepest eccentricity growth from a circular orbit takes place starting from ψ = π/2 ∨ 3π/2 and following the stable direction associated with a hyperbolic equilibrium point or a libration curve associated with an elliptic equilibrium point. Based on this and following the idea presented in [13] for the planar case, the initial conditions that can enable a natural deorbiting from a circular orbit in the spatial case can be obtained by solving the following equation where η cr = 1 − e 2 cr , and the coefficients T j are computed using c i =

Ln2ηcr+Π Ln1ηcr
and s i = 1 − c 2 i . The ± sign depends on whether the critical eccentricity is attained at ψ = 0 or ψ = π, respectively. Note that only the positive solutions is admissible, because it represents a physical area-to-mass ratio. Moreover, it should be noticed that the solution does not provide any information on the time required to cover the invariant curve up to e = e cr . In the above expression, the only unknown is Π that can be set considering that the deorbiting starts from e 0 = 0 and ψ 0 = π/2 ∨ 3π/2, following a resonant curve. Given the value of semi-major axis a, the value of the inclination can be found by computing the resonant conditionψ = 0 at that configuration. Following [1], the resonant conditionψ = 0 at ψ = π/2∨3π/2 does not depend on the solar radiation pressure, but only on the oblateness effect. In particular, i 0 can be obtained by solving the quadratic equation in cos i 0 [1] with

Results
We have applied the procedure just described to the range a ∈ [6978 : 15000] km, considering a discretization of ∆a = 10 km. The lower limit for the range in a reflects the fact that below about 600 km of altitude a circular orbit decays naturally in 25 years due to the effect of the atmospheric drag.
In Fig. 6, we show the area-to-mass ratio computed by means of Eq. (25), as a function of initial semi-major axis for the six resonant arguments, for prograde orbits. The color reports if the deorbiting occurs at ψ = 0 (red) or ψ = π (green). It should be noticed that, as already detected numerically [3,18], the resonant terms j = 5, 6 are the least effective ones, requiring very high values of A/m. Moreover, the conditions depicted always correspond to a libration curve.
In Fig. 7, we show, on the left, the value of initial inclination and semimajor axis corresponding to the conditions depicted in Fig. 6 (right). On the right, we show with a thicker point the configurations that actually attain e cr in less than 50 years. The cases where this does not occur are associated to two different behaviors: either the dynamics is too slow and this happens, in particular, after the cusps that can be noticed in the figure; -or there exist two curves, not connected, associated with the same value of H, for a same integral of motion Λ. This happens, in particular, for j = 1 and j = 4 before for values of semi-major axis lower than a = 7500 km.  Fig. 7 Left: semi-major axis and inclination for the circular orbits corresponding to the minimum area-to-mass depicted in Fig. 6 on the right. Right: the blue thicker points correspond to the conditions that can deorbit in less than 50 years.

Conclusions
In this work, we have extended the analytical development of the equations of motion associated with the coupled oblateness -solar radiation pressure effect. In particular, the Hamiltonian formulation, only partially explained in [2,1] has been here accurately presented. On this basis, the whole phase space in Earth orbit up to a = 20000 km has been described by means of a dynamical taxonomy, that shows where the main bifurcations occurs, in particular the saddle-node and the transcritical ones. Finally, the analytical value of the area-to-mass ratio required to de-orbit has been provided, for the first time in the three-dimensional case, and the corresponding applications for a feasible disposal strategy by means of a solar sail has been presented.