Phase space description of the dynamics due to the coupled effect of the planetary oblateness and the solar radiation pressure perturbations

Abstract

The aim of this work is to provide an analytical model to characterize the equilibrium points and the phase space associated with the singly averaged dynamics caused by the planetary oblateness coupled with the solar radiation pressure perturbations. A two-dimensional differential system is derived by considering the classical theory, supported by the existence of an integral of motion comprising semi-major axis, eccentricity and inclination. Under the single resonance hypothesis, the analytical expressions for the equilibrium points in the eccentricity-resonant angle space are provided, together with the corresponding linear stability. The Hamiltonian formulation is also given. The model is applied considering, as example, the Earth as major oblate body, and a simple tool to visualize the structure of the phase space is presented. Finally, some considerations on the possible use and development of the proposed model are drawn.

1 Introduction

The main objective of this work is to derive an analytical model in the three-dimensional case of the equilibrium points associated with the singly averaged dynamics induced by the coupled effect of the solar radiation pressure (SRP) and the planetary oblateness. As far as we know, such derivation does not exist at the moment and it can represent a fundamental tool in different fields. In particular, the stationary configurations, the invariant curves in the libration regions corresponding to the elliptic points and the hyperbolic invariant manifolds associated with the hyperbolic points can have several applications in planetary dynamics (e.g., Mignard 1982; Krivov et al. 1996) and in mission design (e.g., Colombo and McInnes 2012b; Wittig et al. 2017).

The literature shows that the pioneers in the identification of the role of the SRP effect coupled with the zonal harmonics $J_2$ of a primary body in the orbital evolution of a small body were Musen (1960) and Cook (1962). They developed the corresponding singly averaged equations of motion, for a satellite orbiting the Earth, in terms of Keplerian orbital elements, and remarked the existence of six orbital resonances involving the rate of precession of the longitude of the ascending node, of the argument of pericenter and the apparent mean motion of the Sun. Later on, Hughes (1977) expanded the disturbing function associated with SRP using the Kaula’s method up to high-order terms, and provided some examples on whether they can be relevant for satellites in Low Earth Orbit (LEO). Recall also that Breiter (1999, 2001) addressed the analytical treatment of lunisolar perturbations in canonical coordinates, considering not only the resonant condition, but more in general the critical points associated with the dynamical system, together with their stability.¹

In terms of applications, in the vast literature on the dynamics of dust in planetary systems, we can mention Krivov et al. (1996) and Hamilton and Krivov (1996), who studied the eccentricity oscillation associated with the singly averaged disturbing potential of SRP and $J_2$ for a dust particle orbiting a planet. Krivov and Getino (1997) applied the same model, but in the planar case and only considering one of the resonant terms, also for obtaining a phase space representation of the orbit evolution of low-inclination Earth orbits, in the perspective of high-altitude balloons.

Colombo et al. (2012) performed a parametric analysis of the SRP-$J_2$ phase space for different values of semi-major axis and area-to-mass ratio, identifying the equilibrium points corresponding to either heliotropic (apogee pointing toward the Sun) or anti-heliotropic (perigee pointing toward the Sun) orbits. The equilibrium conditions computed analytically for the planar case were numerically extended to nonzero inclinations. For such frozen orbits, different applications were proposed: the use of heliotropic orbits for Earth observation in the visible wavelength (Colombo et al. 2012; Colombo and McInnes 2012a), or the use of heliotropic orbits at an oblate asteroid (Lantukh et al. 2015), or also the use of anti-heliotropic orbits for geomagnetic tail exploration missions (McInnes et al. 2001; Oyama et al. 2008; Luo et al. 2018). Following (Krivov and Getino 1997) a particular application was proposed by Lücking et al. (2012, 2013) for passive deorbiting of spacecraft at the end of life. Given the mass parameter of the main body, the dynamics under consideration depends on three constants, the semi-major axis of the orbit, a given integral of motion and the area-to-mass ratio of the small body. In Lücking et al. (2012, 2013) and Colombo and de Bras de Fer (2016), the area-to-mass required to obtain a given eccentricity increase was computed by means of an optimization procedure, considering also the associated time.

With the same aim of looking for natural perturbations that can support a reentry from LEO at the end of life, Alessi et al. (2018) performed an extensive numerical mapping of the region, and identified deorbiting natural corridors, defined in terms of inclination and eccentricity for given semi-major axis. In Schettino et al. (2019), a numerical frequency portrait of the LEO region was presented, highlighting also the role that high-order resonances associated with SRP (Hughes 1977) might have.

The other typical example of mission that can be designed on the basis of the coupled SRP-$J_2$ dynamics is a mission to an asteroid or a comet (e.g., Russell 2012; Scheeres 2012).

In this work, the three-dimensional equations of motion which describe the variation in eccentricity, inclination, and a given angle accounting for the motion of the longitude of the ascending node and the argument of pericenter are analyzed, they are reduced to a two-dimensional case by means of a well-known integral of motion, and the analytical expressions for resonant conditions and equilibrium points are provided. The corresponding stability is given for the six main resonances, and the possible different phase space portraits for orbits at the Earth are described. Finally, the tools presented are analyzed in the perspective of a practical exploitation.

Throughout the work, it is always assumed that the major body is the Earth, but the analytical treatment developed is general and it can be applied to any oblate body. In our formulation, we adopted the Earth equatorial plane as reference plane, because our major interest is a possible exploitation for deorbiting. We notice, however, that for other applications, e.g., Scheeres (2012), the ecliptic reference system can be more suitable.

Moreover, the concept proposed will hold certainly in the cases described above, where the two perturbations—oblateness and solar radiation pressure—play a major role, but it will be of high interest to see how the solutions presented may persist under the effect of an additional effect and under which conditions they remain the skeleton of the long-term dynamics. For orbits at the Earth, the key case will be the dynamics due to lunisolar gravitational perturbations for Highly Elliptical Orbits (Colombo 2019) and in the geosynchronous region (e.g., Casanova et al. 2015; Gkolias and Colombo 2019).

2 Dynamical model

Let us assume that a small body (e.g., a spacecraft) moves under the effect of the Earth’s gravitational monopole, the Earth’s oblateness and the solar radiation pressure. In particular, for the SRP it is assumed in the cannonball model that the orbit of the spacecraft is entirely in sunlight and that the effect of the Earth’s albedo is negligible (e.g., Krivov et al. 1996).

In the following, $(a,e,i,{\varOmega },\omega )$ are the Keplerian orbital elements of the small body measured with respect to the Earth equatorial plane, n the mean motion of the small body, $\mu$ the Earth gravitational parameter, $J_2$ the second zonal term of the geopotential, $r_{\oplus }$ the equatorial radius of the Earth, $C_{\oplus }=\mu J_2r^2_{\oplus }$, $\epsilon$ the obliquity of the ecliptic, P the solar radiation pressure at 1 AU, $c_R$ the reflectivity coefficient², A / m the area-to-mass ratio of the small body, $C_\mathrm{{SRP}}=\frac{3}{2}Pc_R\frac{A}{m}$, $\lambda _S$ the longitude of the Sun measured on the ecliptic plane, and the resonant argument is

$$\begin{aligned} \psi _j=n_1{\varOmega }+n_2\omega + n_3 \lambda _S, \quad \quad \quad \psi _j\in [0,2\pi ], \end{aligned}$$

(1)

with $n_1, n_2, n_3$ being as in Table 1, where we provide also the correspondence with respect to the notation used by Cook (1962).

Table 1 Argument $\psi _j=n_1{\varOmega }+n_2\omega + n_3 \lambda _S$ of the periodic component in terms of $n_1, n_2, n_3$

Under these hypotheses, following the description given in Krivov et al. (1996), Colombo et al. (2012), Colombo and McInnes (2012a), Alessi et al. (2018), the singly averaged equations of motion of the spacecraft can be written as

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}a}{\mathrm{d}t}&=0,\\ \frac{\mathrm{d}e}{\mathrm{d}t}&=-C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{na}\sum _{j=1}^6 {\mathcal {T}}_{j}{\partial {\cos {\psi }_j}\over \partial \omega },\\ \frac{\mathrm{d}i}{\mathrm{d}t}&=-C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}\sin i}\sum _{j=1}^6 {\mathcal {T}}_{j}\left( {\partial {\cos {\psi }_j}\over \partial {\varOmega }}-\cos i{\partial {\cos {\psi }_j}\over \partial \omega }\right) ,\\ \frac{\mathrm{d}{\varOmega }}{\mathrm{d}t}&={{\dot{{\varOmega }}}}_{J_2}+{{\dot{{\varOmega }}}}_\mathrm{{SRP}},\\ \frac{\mathrm{d}\omega }{\mathrm{d}t}&={{\dot{\omega }}}_{J_2}+{{\dot{\omega }}}_\mathrm{{SRP}}, \end{aligned} \end{aligned}$$

(2)

where

$$\begin{aligned} \begin{aligned} {\mathcal {T}}_1&=\cos ^2\left( \frac{\epsilon }{2}\right) \cos ^2\left( \frac{i}{2}\right) , \\ {\mathcal {T}}_2&=\cos ^2 \left( \frac{\epsilon }{2}\right) \sin ^2\left( \frac{i}{2}\right) ,\\ {\mathcal {T}}_3&=\frac{1}{2} \sin (\epsilon ) \sin (i), \\ {\mathcal {T}}_4&=-\frac{1}{2} \sin (\epsilon ) \sin (i),\\ {\mathcal {T}}_5&=\sin ^2\left( \frac{\epsilon }{2}\right) \cos ^2\left( \frac{i}{2}\right) ,\\ {\mathcal {T}}_6&=\sin ^2\left( \frac{\epsilon }{2}\right) \sin ^2\left( \frac{i}{2}\right) , \end{aligned} \end{aligned}$$

(3)

and

$$\begin{aligned} \begin{aligned} {{\dot{{\varOmega }}}}_{J_2}&= -\frac{3}{2}\frac{J_2r^2_{\oplus }n}{a^2(1-e^2)^2}\cos i,\\ {{\dot{{\varOmega }}}}_\mathrm{{SRP}}&=C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}\sin i}\sum _{j=1}^6{\partial {{\mathcal {T}}}_{j}\over \partial i}\cos {\psi }_j,\\ {{\dot{\omega }}}_{J_2}&= \frac{3}{4}\frac{J_2r^2_{\oplus }n}{a^2(1-e^2)^2}\left( 5\cos ^2i-1\right) ,\\ {{\dot{\omega }}}_\mathrm{{SRP}}&=C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{nae}\sum _{j=1}^6 {\mathcal {T}}_{j}\cos {\psi }_j-{{\dot{{\varOmega }}}}_\mathrm{{SRP}}\cos i. \end{aligned} \end{aligned}$$

(4)

If we assume that the dynamics is driven by a single term j at a time, that means that only one periodic component $\sin \psi _j$ (for e, i; $\cos \psi _j$ for ${\varOmega }, \omega$) affects the motion at a time, then we can simplify the description of the dynamics by analyzing the following system

$$\begin{aligned} \begin{aligned} \frac{\mathrm{d}e}{\mathrm{d}t} \Big |_{j}&=n_2C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{na}{\mathcal {T}}_{j}\sin {\psi }_j,\\ \frac{\mathrm{d}\psi _j}{\mathrm{d}t}&=n_1{{\dot{{\varOmega }}}}_{(J_2,j)}+n_2{{\dot{\omega }}}_{(J_2,j)}+n_3n_S, \end{aligned} \end{aligned}$$

(5)

where $n_S$ is the apparent mean motion of the Sun around the Earth and ${{\dot{{\varOmega }}}}_{(J_2,j)}$ and ${{\dot{\omega }}}_{(J_2,j)}$ are intended to include the contribution of $J_2$ and the jth-contribution of the SRP. In the following, the time derivative of the eccentricity due to the periodic term j, namely $\frac{\mathrm{d}e}{\mathrm{d}t} \Big |_{j}$, will appear as $\dot{e}_{|j}$.

Since we can apply the same argument used in Daquin et al. (2016) for the coupled lunisolar gravitational perturbation and oblateness effect (see Alessi et al. 2018 and Sect. 2.2), the behavior in inclination can be recovered at any time by means of the integral of motion

$$\begin{aligned} {\varLambda }=(n_2\cos {i}-n_1)\sqrt{\mu a(1-e^2)}, \end{aligned}$$

(6)

which can be inverted to get

$$\begin{aligned} \cos {i}=\frac{{\varLambda }}{n_2\sqrt{\mu a(1-e^2)}}+\frac{n_1}{n_2}. \end{aligned}$$

(7)

Note that the evolution of the dynamics is determined by three parameters, the semi-major axis a, the area-to-mass ratio A / m and the integral of motion ${\varLambda }$ (which depends on the dominant $j-$resonant term).

2.1 Resonances and equilibrium points

A resonance takes place when the following condition is satisfied

$$\begin{aligned} {{\dot{\psi }}}_j= 0. \end{aligned}$$

(8)

Notice that this condition implies specific configurations, according to the periodic component j considered. In particular,

• for $j=1$, if Eq. (8) is verified at prograde orbits, then the longitude of the periapsis is sun-synchronous;

• for $j=2$, if Eq. (8) is verified at retrograde orbits, then the longitude of the periapsis is sun-synchronous;

• for $j=3$, the argument of periapsis is sun-synchronous;

• for $j=4$, the argument of periapsis is sun-antisynchronous;

• for $j=5$, if Eq. (8) is verified at prograde orbits, then the longitude of the periapsis is sun-antisynchronous;

• for $j=6$, if Eq. (8) is verified at retrograde orbits, then the longitude of the periapsis is sun-antisynchronous.

Table 2 Coefficients of the quadratic equation $c_1\cos ^2{i}+c_2\cos {i}+c_3=0$ associated with the equilibrium point at $\psi _j=0$ for $j=1,2,5,6$

Note that Eq. (8) depends on $(a,e,i,\psi _j)$, that is, also on $C_\mathrm{{SRP}}$, except if it is computed at $\psi _j=\pi /2$ or $\psi _j=3\pi /2$. In the latter case, the resonant condition can be computed only considering the effect due to $J_2$ on $({\varOmega }, \omega )$, that is, given $(a,e,n_1,n_2,n_3)$, by solving the following quadratic equation in $\cos i$

$$\begin{aligned} c_1\cos ^2{i}+c_2\cos {i}+c_3=0, \end{aligned}$$

(9)

where

$$\begin{aligned} \begin{aligned} c_1&=\frac{15n_2J_2r^2_{\oplus }n}{4a^2(1-e^2)^2},\\ c_2&=-\frac{3n_1J_2r^2_{\oplus }n}{2a^2(1-e^2)^2},\\ c_3&=n_3n_S-\frac{3n_2J_2r^2_{\oplus }n}{4a^2(1-e^2)^2}. \end{aligned} \end{aligned}$$

(10)

These are the resonant conditions depicted by Cook (1962) in his Fig. 5 and remarked in Alessi et al. (2018), Schettino et al. (2019) and Alessi et al. (2018).

On the other hand, an equilibrium point is computed if both equations in Eq. (5) cancel out. This must happen in particular at $\psi _j=0$ or $\psi _j=\pi$. For $j=1,2,5,6$, discarding the singularity at $e=0$ and $e=1$, the equation that gives the value of inclination corresponding to the location of the equilibrium points at $\psi _j=0$ and at $\psi _j=\pi$, given (a, e, A / m), is also a quadratic equation—Eq. (9)—in $\cos {i}$ that can be solved. The corresponding analytical coefficients are shown in Tables 2 and 3 for all the cases. For $j=3,4$, discarding also the singularity corresponding to $i=0$, the equation is instead a cubic equation of the type $\sin ^3{i}+s_1\sin ^2{i}+s_2\sin {i}+s_3=0$, that can also be solved analytically via classical techniques. In Tables 4 and 5, we show the coefficients of such cubic equations.

Table 3 Coefficients of the quadratic equation $c_1\cos ^2{i}+c_2\cos {i}+c_3=0$ associated with the equilibrium point at $\psi _j=\pi$ for $j=1,2,5,6$

Table 4 Coefficients of the cubic equation $\sin ^3{i}+s_1\sin ^2{i}+s_2\sin {i}+s_3=0$ associated with the equilibrium point at $\psi _j=0$ for $j=3$ and $j=4$

Table 5 Coefficients of the cubic equation $\sin ^3{i}+s_1\sin ^2{i}+s_2\sin {i}+s_3=0$ associated with the equilibrium point at $\psi _j=\pi$ for $j=3$ and $j=4$

The stability of the points can be evaluated by computing the eigenvalues of the matrix

$$\begin{aligned} J=\begin{pmatrix} {\partial {\dot{e}|j}\over \partial e} &{}\quad {\partial {\dot{e}|j}\over \partial {\psi _j}} \\ {\partial {{{\dot{\psi }}}_j}\over \partial e} &{}\quad {\partial {{{\dot{\psi }}}_j} \over \partial { \psi _j}} \end{pmatrix}, \end{aligned}$$

(11)

at the given equilibrium point, where

$$\begin{aligned} \begin{aligned} {\partial {\dot{e}_{|j}}\over \partial e}&=0,\\ {\partial {\dot{e}_{|j}}\over \partial {\psi _j}}&=n_2C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{na}{\mathcal {T}}_{j}\cos {\psi }_j,\\ {\partial {{{\dot{\psi }}}_j}\over \partial e}&=n_1{\partial {\dot{\varOmega }_{(J_2,j)}}\over \partial e}+n_2{\partial {\dot{\omega }_{(J_2,j)}}\over \partial e},\\ {\partial {{{\dot{\psi }}}_j} \over \partial { \psi _j}}&=0. \end{aligned} \end{aligned}$$

(12)

Note that the partial derivatives appearing in the third equation of Eq. (12) can be evaluated by taking into account Eq. (7), namely

$$\begin{aligned} {\partial {{{\dot{{\varOmega }}}}_{(J_2,j)}}\over \partial e}= & {} -\,6\frac{J_2r^2_{\oplus }ne}{a^2(1-e^2)^3}\cos i+\frac{C_\mathrm{{SRP}}}{na(1-e^2)^{3/2}\sin i}{\partial {{\mathcal {T}}}_{j}\over \partial i}\cos {\psi }_j \\&-\,\frac{3}{2}\frac{J_2r^2_{\oplus }n}{a^2(1-e^2)^2} {\partial {\cos {i}}\over \partial e} + C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}}{\partial {{\mathcal {T}}}_{j}\over \partial i}\cos {\psi }_j{\partial \over \partial e}\frac{1}{\sin {i}}\\&+\,C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}\sin i}\cos {\psi }_j{\partial \over {\partial e}}{{\partial {{\mathcal {T}}}_{j}\over \partial i}},\\ {\partial {{{\dot{\omega }}}_{(J_2,j)}}\over \partial e}= & {} 3\frac{J_2r^2_{\oplus }ne}{a^2(1-e^2)^3}(5\cos ^2{i}-1) \\&-\, \frac{C_\mathrm{{SRP}}}{na\sqrt{1-e^2}}\left( \frac{{\mathcal {T}}_{j}}{e^2}- \frac{\cos {i}}{(1-e^2)^2\sin i}{\partial {{\mathcal {T}}}_{j}\over \partial i}\right) \cos {\psi }_j \\&+\,\frac{15}{2}\frac{J_2r^2_{\oplus }n\cos {i}}{a^2(1-e^2)^2} {\partial {\cos {i}}\over \partial e}+C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{nae}\cos {\psi }_j{\partial {{\mathcal {T}}_{j}} \over \partial {e}} \\&-\,C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}}{\partial {{\mathcal {T}}}_{j}\over \partial i}\cos {\psi }_j{\partial \over \partial e}\frac{1}{\tan {i}} \\&+\,C_\mathrm{{SRP}}\frac{e}{na\sqrt{1-e^2}\tan {i}} \cos {\psi }_j{\partial \over {\partial e}}{{\partial {{\mathcal {T}}}_{j}\over \partial i}}, \end{aligned}$$

where

$$\begin{aligned} {\partial {\cos {i}}\over \partial e}= & {} \frac{e{\varLambda }}{n_2\sqrt{\mu a}(1-e^2)^{3/2}},\\ {\partial {\sin {i}}\over \partial e}= & {} \frac{1}{\tan {i}}{\partial {\cos {i}}\over \partial e}, \end{aligned}$$

and ${\partial {{\mathcal {T}}}_{j}\over \partial i}$ and ${\partial \over {\partial e}}{{\partial {{\mathcal {T}}}_{j}\over \partial i}}$ can be derived from these expressions by applying the half-angle trigonometric formulae to Eq. (3).

2.2 Hamiltonian formulation

For completeness, we provide also the Hamiltonian formulation of the dynamics under study. To this end, we first recall that the equations of motion Eqs. (2) and (4) are obtained from the singly averaged disturbing potential associated with the oblateness and the solar radiation pressure effect, as described in Krivov et al. (1996), Colombo et al. (2012), Colombo and McInnes (2012a), Alessi et al. (2018). Notice that the disturbing potential associated with the SRP also corresponds to Eq. (17) in Casanova et al. (2015) and to Eq. (69) in Scheeres (2012). Then, we can apply the argument developed in Daquin et al. (2016) for lunisolar gravitational perturbations, adapted to the solar radiation pressure effect, under the assumption that only one periodic term j is relevant at a time for the motion of the small body. Considering, in particular, one of the two canonical transformations developed in Daquin et al. (2016) written in terms of classical Delaunay variables $\left( L, G, H, l, g, h \right)$, by introducing the canonical variables $\left( {\varGamma }, \tau \right)$ such that ${{\dot{\tau }}}\equiv {\partial {{\mathcal {H}}}\over \partial {\varGamma }}=\lambda _S$, we can formulate the Hamiltonian of the problem in action angle variables $(\Sigma , \sigma )$. We have³

$$\begin{aligned} \begin{aligned} \Sigma _1&\equiv \Sigma =n_2^{-1}G, \quad \quad \quad \quad \quad ~ \sigma _1\equiv \sigma =n_2g+n_1h+n_3\tau ,\\ \Sigma _2&\equiv {\varLambda }=-n_1G+n_2H, \quad \quad ~\sigma _2=n_2^{-1}h, \\ \Sigma _3&= -n_2^{-1}n_3G+{\varGamma }, \quad \quad \quad ~~ \sigma _3=\tau . \end{aligned} \end{aligned}$$

(13)

Note that $\sigma _1\equiv \sigma \equiv \psi _j$ and that Eq. (6) is derived from $\Sigma _2$.

In this way, we obtain

$$\begin{aligned} {\mathcal {H}}={\mathcal {H}}_{J_2}+{\mathcal {H}}_\mathrm{{SRP}}, \end{aligned}$$

(14)

where the part associated with $J_2$ is

$$\begin{aligned} {\mathcal {H}}_{J_2}=\frac{J_2 \mu ^4 r^2_{\oplus } \left( -3 \Sigma _1^2 n_1^2 n_2^2+\Sigma _1^2 n_2^4-6 \Sigma _1 \Sigma _2 n_1 n_2-3 \Sigma _2^2\right) }{4 L^3 \Sigma _1^5 n_2^7}, \end{aligned}$$

(15)

and the one due to SRP is

$$\begin{aligned} {\mathcal {H}}_\mathrm{{SRP}}=-\frac{C_\mathrm{{SRP}}}{2\mu } L^2 {\mathcal {S}}_j\sqrt{1-\frac{n_2^2\Sigma _1^2}{L^2}} \cos \sigma _1+n_S\left( n_3\Sigma _1+\Sigma _3\right) , \end{aligned}$$

(16)

where ${\mathcal {S}}_j$ depends on the resonance, namely

$$\begin{aligned} \begin{aligned} {\mathcal {S}}_1&=\cos ^2\left( \frac{\epsilon }{2}\right) \left( \text {ci}+1\right) ,\\ {\mathcal {S}}_2&=\cos ^2\left( \frac{\epsilon }{2}\right) \left( -\text {ci}+1\right) ,\\ {\mathcal {S}}_3&=\sin \epsilon \sqrt{-\text {ci}^2+1},\\ {\mathcal {S}}_4&=-\sin \epsilon \sqrt{-\text {ci}^2+1},\\ {\mathcal {S}}_5&=\sin ^2\left( \frac{\epsilon }{2}\right) \left( \text {ci}+1\right) ,\\ {\mathcal {S}}_6&=\sin ^2\left( \frac{\epsilon }{2}\right) \left( -\text {ci}+1\right) , \end{aligned} \end{aligned}$$

(17)

with $\text {ci}=\frac{n_1n_2\Sigma _1+\Sigma _2}{n_2^2\Sigma _1}$.

This is also to say that the equations of motion can be also written as

$$\begin{aligned} \begin{aligned} {{\dot{\Sigma }}}&=-{\partial {{\mathcal {H}}\left( \Sigma ,\sigma \right) }\over \partial {\sigma }},\\ {{\dot{\sigma }}}&={\partial {{\mathcal {H}}\left( \Sigma ,\sigma \right) }\over \partial {\Sigma }}. \end{aligned} \end{aligned}$$

(18)

Note that the Hamiltonian function is a constant of motion and so it is $\Sigma _3$, since ${\mathcal {H}}$ does not depend on $\sigma _3$. In other words, as already noticed in Breiter (1999) for lunisolar perturbations, in Eq. (16) $\Sigma _3$ is a dummy variable and its value is not important.

Fig. 1

As a function of semi-major axis and eccentricity, for $j=1$ we show the inclination solutions (colorbar in degrees) of the quadratic equation, which gives the equilibrium points associated with $\psi =0$ (left two columns), and with $\psi =\pi$ (right two columns). The coefficients of the quadratic equation are given in Tables 2 and 3, respectively. Top row: $A/m=0.012 \, {\hbox {m}}^2/\hbox {kg}$; central row: $A/m=1 \, {\hbox {m}}^2/\hbox {kg}$; bottom row: $A/m=20 \, {\hbox {m}}^2/\hbox {kg}$

3 Dynamical description

In Figs. 1, 2, 3 and 4, we show the location of the equilibrium points at $\psi _j=0$ and $\psi _j=\pi$, associated with each term $j=1,2,5,6$, as a function of (a, e, i) for three values of area-to-mass ratio. These are $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$, $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ and $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$, characteristic of a standard satellite, a satellite equipped with an area-augmentation device with the present technology (Colombo et al. 2018) and a high area-to-mass ratio debris fragment, respectively.

Fig. 2

The same as in Fig. 1, but for $j=2$

Fig. 3

The same as in Fig. 1, but for $j=5$

Fig. 4

The same as in Fig. 1, but for $j=6$

Fig. 5

As a function of semi-major axis and eccentricity, for $j=3$ we show the inclination solutions (colorbar in degrees) of the cubic equation, which gives the equilibrium points associated with $\psi =0$ (left two columns), and with $\psi =\pi$ (right two columns). The coefficients of the cubic equation are given in Tables 4 and  5, respectively. Top row: $A/m=0.012 \, {\hbox {m}}^2/\hbox {kg}$; central row: $A/m=1 \, {\hbox {m}}^2/\hbox {kg}$; bottom row: $A/m=20 \, {\hbox {m}}^2/\hbox {kg}$

Fig. 6

The same as in Fig. 5, but for $j=4$

The solutions corresponding to inclination equal to zero, for $j=1$ and $j=2$, are equivalent to the ones described in Colombo et al. (2012). For $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ and $j=1$ or $j=2$, the solutions presented in Colombo and McInnes (2012a), for the inclined case when considering only $J_2$, are equivalent to the ones shown here, because for low A / m values the effect of the solar radiation pressure on $({\varOmega }, \omega )$ is negligible.

Notice how the inclination changes as the semi-major axis increases, increasing or decreasing depending on j and on whether the equilibrium is at $\psi _j=0$ or $\psi _j=\pi$. Note also that, by increasing the area-to-mass ratio, the main effect is to enlarge the range of semi-major axis where the effect can be detected. In general, for quasi-circular orbits, the dynamics is coupled up to about $a=15000\, {\hbox {km}}$, except for the cases $\psi _1=0$ and $\psi _2=\pi$ of the debris fragment, for which the effect can be considerable up to Geosynchronous Earth Orbits (GEO) altitudes. In Figs. 5, 6, we show the inclination values corresponding to the equilibrium points for $j=3$ and $j=4$, as a function of (a, e) for the same three values of area-to-mass ratio. The inclination shown is computed by considering one of the solutions of the cubic equations displayed in Tables 4 and 5. The other solutions correspond to the opposite value shown and to the singularity $i=0.$

The dynamics under consideration is organized according to the equilibrium points which exist for a given $(a,{\varLambda },A/m,j)$ combination. In other words, to obtain a single portrait of the behavior of the phase space, the information displayed in the previous figures has to be accompanied by the information on ${\varLambda }$. In particular, given (a, A / m, j), it is possible to distinguish regions in terms of number and type of equilibrium points, by computing the minimum and maximum values attained by the curves, which represent the evolution of ${\varLambda }$, as a function of e, at $\psi _j=0$ and $\psi _j=\pi$. The type of equilibrium point is determined by its linear behavior: from the analysis performed for this work for orbits around the Earth (from LEO to GEO altitudes), the point can be either a center or a saddle. Note that A / m is constant as given by the problem, a is constant under the effects considered (see Eq. (2)), while the dominant resonant term j might change along a curve, as it will be shown in the following.

Fig. 7

${\tilde{{\varLambda }}}$ as a function of eccentricity, at the equilibrium points existing for $j=1$, $A/m=1 \, {\hbox {m}}^2/\hbox {kg}$ and $a=8078 \, \hbox {km}$. On the right a close-up view, showing also the number of equilibrium points in the different regions. Blue: stable point at $\psi _1=0$; red: unstable point at $\psi _1=0$. Green: stable point at $\psi _1=\pi$; ocher: unstable point at $\psi _1=\pi$

The first example is presented in Fig. 7. For $j=1$, $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ and $a=8078\, {\hbox {km}}$, the value of ${\tilde{{\varLambda }}}\equiv {\varLambda }/\sqrt{\mu }$ computed at the equilibrium points $\psi _1=0$ and $\psi _1=\pi$ is depicted as a function of the corresponding value of eccentricity, considering the solution for $\cos {i}$ on the first column of Fig. 1. We display ${\tilde{{\varLambda }}}$ instead of ${\varLambda }$ for the sake of clarity to avoid large numbers. In blue, the stable equilibrium points at $\psi _1=0$, and in red the unstable ones. In green the stable equilibrium points at $\psi _1=\pi$, and in ocher the unstable ones. On the right, a closer view is given, indicating also the number of equilibrium points existing in the different regions. In Fig. 8, the same information is shown, but displays also the corresponding value of inclination (dashed lines, secondary x-axis). Following the two figures, we can see that for $j=1$, $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ and $a=8078\, {\hbox {km}}$, the phase space can be structured according to a minimum of 1 equilibrium point and to a maximum of 3 equilibrium points, which can exist in the range $e\in [0:1]$, $i\in [4:64]$ deg. Note the narrow range of inclination, $i\in [39.8:40.8]$ deg, within which the equilibrium point at $\psi _1=0$ can be unstable. Moreover, up to ${\tilde{{\varLambda }}}\simeq -20.55 \, {\hbox {km}}^{1/2}$ there exists only one equilibrium point at $\psi _1=0$, which is stable (see Fig. 9 top left); from ${\tilde{{\varLambda }}}\simeq -20.55 \, {\hbox {km}}^{1/2}$ to ${\tilde{{\varLambda }}}\simeq -20.48 \, {\hbox {km}}^{1/2}$ there are 2 stable equilibrium points at $\psi _1=0$ and one unstable equilibrium point at $\psi _1=0$ (see Fig. 9 top right, note also the homoclinic connection); from ${\tilde{{\varLambda }}}\simeq -20.48 \, {\hbox {km}}^{1/2}$ to ${\tilde{{\varLambda }}}\simeq -20.44 \, {\hbox {km}}^{1/2}$, there exists one stable equilibrium point at $\psi _1=0$ (see Fig. 9 bottom left); for values of ${\tilde{{\varLambda }}}$ higher than ${\tilde{{\varLambda }}}\simeq -20.44 \, {\hbox {km}}^{1/2}$, there exist one stable equilibrium point at $\psi _1=0$ and one stable and one unstable equilibrium points at $\psi _1=\pi$ (see Fig. 9 bottom right). At the transitions, namely at ${\tilde{{\varLambda }}}\simeq -20.55 \, {\hbox {km}}^{1/2}$, ${\tilde{{\varLambda }}}\simeq -20.48 \, {\hbox {km}}^{1/2}$ and ${\tilde{{\varLambda }}}\simeq -20.44 \, {\hbox {km}}^{1/2}$ there is a bifurcation: on the horizontal lines, for increasing ${\tilde{{\varLambda }}}$, either there start appearing 2 additional equilibrium points, one stable and one unstable, or 2 existing equilibrium points with a different linear stability occur to coincide and then disappear.

Fig. 8

${\tilde{{\varLambda }}}$ as a function of eccentricity, at the equilibrium points existing for $j=1$, $A/m=1 \, \hbox {m}^2/\hbox {kg}$ and $a=8078 \, \hbox {km}$. On the right a close-up view. Blue: stable point at $\psi _1=0$; red: unstable point at $\psi _1=0$. Green: stable point at $\psi _1=\pi$; ocher: unstable point at $\psi _1=\pi$. The dashed curves show the value of inclination corresponding to the given ${\varLambda }$, the color code is the same

Fig. 9

The possible phase space portraits for $j=1$, $A/m=1 \, {\hbox {m}}^2/\hbox {kg}$ and $a=8078 \, \hbox {km}$, as indicated in Fig. 7 right. Top left: ${\tilde{{\varLambda }}}=-20.6 \, {\hbox {km}}^{1/2}$; top right: ${\tilde{{\varLambda }}}=-20.5 \, {\hbox {km}}^{1/2}$; bottom left: ${\tilde{{\varLambda }}}=-20.45 \, {\hbox {km}}^{1/2}$; bottom right: ${\tilde{{\varLambda }}}=-20.3 \, {\hbox {km}}^{1/2}$

Fig. 10

On the left a close-up view of the ${\tilde{{\varLambda }}}$ behavior as a function of eccentricity, at the equilibrium points existing for $j=1$, $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ and $a=12078\, {\hbox {km}}$. Blue: stable point at $\psi _1=0$; red: unstable point at $\psi _1=0$. Green: stable point at $\psi _1=\pi$; ocher: unstable point at $\psi _1=\pi$. On the right, the phase space portrait structured in 5 equilibrium points for the same values of a, assuming ${\tilde{{\varLambda }}}=-10 \, {\hbox {km}}^{1/2}$. The right plot is split into two to highlight the phase space at $e\approx 0.955$

Notice that Fig. 8 can give also the information on the range of inclination within which we can have a given number of equilibrium points with the corresponding stability, and how the inclination changes according to the eccentricity, for given semi-major axis and ${\tilde{{\varLambda }}}$. In other words, the inclination corresponding to a given phase space portrait, such as the ones in Fig. 9, can be inferred by looking into Fig. 8. A preliminary definition of this domain was considered in Schettino et al. (2017), taking into account only the strictly resonant behaviors for almost zero eccentricities.

For $j=1$, $A/m=1 \, {\hbox {m}}^2/\hbox {kg}$, by increasing the semi-major axis, the number of possible equilibrium points can increase. Spanning a range of a up to the GEO ring, the maximum number is 5, and it is found already at $a=8178\, {\hbox {km}}$. In Fig. 10, we show the corresponding behavior of the $(e,{\tilde{{\varLambda }}})$ curves and the phase space portrait associated with 5 equilibrium points for $a=12078\, {\hbox {km}}$. Moreover, we note that $(e,{\tilde{{\varLambda }}})$ curves cannot be continuous, because Eq. (7) sets specific constraints on the range of ${\tilde{{\varLambda }}}$ ($\cos {i}\in [-1,1]$).

Looking into the $(e,{\tilde{{\varLambda }}})$ curves for $j=1$, always considering the solution for $\cos {i}$ on the first column of Fig. 1, for $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$ it is noticed a similar behavior, while for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ the number of equilibrium points can be 4 also in non-degenerate configurations. Taking as reference the phase space portrait depicted in Fig. 10 on the right, for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ and given ${\tilde{{\varLambda }}}$, the equilibrium point existing at the lowest eccentricity can disappear, and there persist one unstable equilibrium at $\psi _1=0$ and one stable equilibrium at $\psi _1=\pi$ corresponding to a same eccentricity, and one stable equilibrium at $\psi _1=0$ and one unstable equilibrium at $\psi _1=\pi$ corresponding to a different eccentricity.

With respect to the other resonant behaviors, the following is stressed. Generally speaking, it appears that the dynamics associated with the resonant terms #1, #2 and #4 is richer than for the other three resonant terms, in terms of maximum possible number of equilibrium points. Also, the value of area-to-mass ratio can change the value of semi-major axis at which a bifurcation can occur.

Resonant term #2 For $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$, considering the solution for $\cos {i}$ on the second column of Fig. 2, the curve corresponding to $\psi _2=0$ is higher than the curve corresponding to $\psi _2=\pi$, contrary to what happens for $j=1$. For the lower values of semi-major axis, the maximum possible number of equilibrium points is 3: one stable at $\psi _2=\pi$ and one stable and one unstable at $\psi _2=0$, that is, as in Fig. 9 on the bottom right, but the curves are shifted by $\pi$ in $\psi$. At $a=9878\, {\hbox {km}}$, it can appear also one unstable equilibrium at $\psi _2=\pi$ and the 3 equilibrium points can be as before or all associated with $\psi _2=\pi$, as in Fig. 9 on the top right, but again shifted by $\pi$ in $\psi$. From about $a=10078\, {\hbox {km}}$, the phase space can be structured around 5 equilibrium points, as in Fig. 10 on the right, but also shifted by $\pi$ in $\psi$. For $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$, the same behavior is found, while for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$, the only difference is that, analogously to what happens for the resonant term #1, it is possible to have 4 equilibrium points, as the one associated with the lowest eccentricity drifts as much toward $e=0$ as to disappear.

Resonant term #5 The dynamics is analogous to what described for $j=1$ except that we do not ever notice 5 equilibrium points and that there exist large regions where the number is only 2: one stable at $\psi _5=0$ and one unstable at $\psi _5=\pi$, or vice versa. In particular, for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ and $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$, this is always the case from about $a=19000\, {\hbox {km}}$.

Resonant term #6 The dynamics is analogous to what described for $j=2$, but, analogously to the case $j=5$, the maximum number of equilibria is 3 and in many cases we can see only 2. In particular, for the three values of area-to-mass explored this is always the case from about $a=16500\, {\hbox {km}}$.

Resonant term #3 For $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$, there are 3 equilibrium points: one stable equilibrium point at $\psi _3=0$ at a low eccentricity, and one unstable at $\psi _3=0$ and one stable at $\psi _3=\pi$. This occurs until about $a=15500\, {\hbox {km}}$ (about $a=15000\, {\hbox {km}}$ for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ and about $a=17000\, {\hbox {km}}$ for $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$), where the equilibrium at the lowest eccentricity disappears. From that value of semi-major axis on, the 2 equilibrium points drift toward higher values of eccentricity, analogously to what happens for the other resonant terms.

Resonant term #4 In this case, the dynamics can be as rich as in the case of resonant terms #1 and #2, in the sense that there can exist up to 5 equilibrium points: one stable at $\psi _4=0$ at a low eccentricity, one unstable at $\psi _4=0$ and one stable at $\psi _4=\pi$, and then one stable at $\psi _4=0$ and one unstable at $\psi _4=\pi$ at increasing eccentricity. This occurs for $A/m=0.012 \, {\hbox {m}}^2/{\hbox {kg}}$ until about $a=10000\, {\hbox {km}}$, for $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ until about $a=11600\, {\hbox {km}}$, for $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$ until about $a=19000\, {\hbox {km}}$. From that value on, there persist 4 equilibrium points, which are associated with increasing values of eccentricity for increasing values of semi-major axis.

4 Discussion and future directions

The analysis proposed in this paper provides the fundamental ingredients to describe the dynamics of a body with a high area-to-mass ratio, which orbits a given major oblate body, and which is subject to solar radiation pressure effects.

In particular, it is possible to apply all the dynamical systems theory tools (e.g., Parker and Chua 1989; Wiggins 2003) which can characterize the elliptic (stable) and hyperbolic (unstable) behaviors associated with the equilibrium points of the system. The information given by the eigenvalues and eigenvectors of the Jacobian matrix, Eq. (11), computed at the corresponding equilibrium point, can be used to compute the invariant librating curves in the neighborhood of the elliptic points and the hyperbolic invariant manifolds in the neighborhood of the hyperbolic points. Concerning the time required to move along a curve, at an elliptic point the two eigenvalues are conjugate pure imaginary, say $\pm i\nu$. The closer the libration curve to the elliptic point, the closer its period to $T=2\pi /\nu$. In Fig. 11, we show the value of the period T in years for the cases shown in Figs. 8, 9 and 10. The asymptotic behavior that can be seen in the figure corresponds to the line of transition between different behaviors—the point of bifurcation described before. As the distance with respect to the elliptic point increases, the time needed to cover the curve also increases, in the limit to reach the hyperbolic manifold, which tends, by definition, asymptotically to the unstable point. Note that the resonant curves corresponding to $\psi _j=\pi /2$ and $\psi _j=3\pi /2$, Eqs. (9)–(10), always act as separatrices between libration and circulation motion, when an equilibrium point at low eccentricity exists.

The initial phase with respect to the Sun, that is, the longitude of the ascending node, the argument of pericenter and the epoch, is crucial, given initial (a, e, i, A / m) to understand the dynamics that the small body can follow. In the neighborhood of one or more equilibrium points, it could be possible to play with the initial phase to move along a given invariant curve in the libration region or along a given separatrix. According to this, the eccentricity evolution will change and so the corresponding time will. Notice (see, e.g., Fig. 9) that even in case of circulation the natural dynamics can provide an eccentricity increase that can be exploited conveniently. With the same idea, it can be envisaged to play with the constants which determine the dynamics, a, A / m or ${\varLambda }$ (or ${\tilde{{\varLambda }}}$), to achieve a given objective, in particular to obtain a bounded eccentricity variation within a given timespan or an escape trajectory.

Fig. 11

For the elliptic points existing for the resonant term #1 and $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$, we show the period T in years corresponding to the libration curves in their neighborhood. Left: $a=8078\, {\hbox {km}}$, ${\tilde{{\varLambda }}}\in [-21:-20] \, {\hbox {km}}^{1/2}$; right: $a=12078\, {\hbox {km}}$, ${\tilde{{\varLambda }}}\in [-20:10] \, {\hbox {km}}^{1/2}$

Moreover, from the analysis provided in Sect. 2, the maximum eccentricity corresponding to a given invariant curve can be computed, as suggested in Alessi and Colombo (2018), by simply looking for the following condition

$$\begin{aligned} \begin{aligned} \dot{e}_{|j}&=0\\ {\frac{\mathrm{d}\dot{e}_{|j}}{\mathrm{d}t}}&=n_2C_\mathrm{{SRP}}\frac{\sqrt{1-e^2}}{na}{\mathcal {T}}_{j}{{\dot{\psi }}}_j\cos {\psi }_j<0, \end{aligned} \end{aligned}$$

(19)

which assumes that the maximum occurs at either $\psi _j=0$ or $\psi _j=\pi$.

On the other hand, each invariant curve, as the ones depicted in Figs. 9 and 10, is characterized by a constant Hamiltonian—Eq. (14)—and this information can be considered to compute a priori the maximum and minimum eccentricity that can be achieved. For instance, departing from the unstable direction associated with a saddle point, corresponding say to $\psi =0$, we can see if an eccentricity corresponding to $\psi =\pi$ can be attained and, in case, its value, by inverting Eq. (14) as a function of $\Sigma$, or $\sqrt{1-e^2}$. A detailed analysis on how this can be done will be the focus of a future investigation.

As it is well known, the single resonance hypothesis assumed in this work can fail to hold, that is, it can occur that two resonant terms can play a comparable role at the same time. This can be inferred from Figs. 1, 2, 3, 4, 5 and 6, which show that, for a same (a, e, i, A / m) combination, there exist equilibrium points associated with different resonant terms. The same situation can be also depicted as an intersection, in the (e, i) plane, between curves corresponding to equilibrium points for given (a, A / m).⁴ In Fig. 12, we show an example for $a=10078\, {\hbox {km}}$ and $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$ and $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$. Note how the stability can change by changing the area-to-mass ratio and also the number of intersections. By increasing the value of semi-major axis, the curves tend to be parallel with respect to the x-axis toward higher values of eccentricity.

The behavior in the regions where two dominant terms can exist will be studied in the future by analyzing the value of the corresponding amplitude and period and also the possible occurrence of a chaotic behavior.

On the other hand, along an invariant curve of the phase space (recall Figs. 9, 10) it is possible to estimate the excursion in eccentricity using Eq. (19) (and an analogous one to compute the minimum eccentricity achieved), and then apply Eq. (7) to see the corresponding excursion in inclination (Schettino et al. 2019). From such estimate, it would be trivial to see if a transition between different resonant terms might take place, by comparison with plots like those in Fig. 12.

Fig. 12

For $a=10078\, {\hbox {km}}$, the curves shown represent the location of the equilibrium points, as a function of inclination and eccentricity for the six resonant terms. Blue: stable point at $\psi _j=0$; red: unstable point at $\psi _j=0$. Green: stable point at $\psi _j=\pi$; ocher: unstable point at $\psi _j=\pi$. Left: $A/m=1 \, {\hbox {m}}^2/{\hbox {kg}}$; right: $A/m=20 \, {\hbox {m}}^2/{\hbox {kg}}$

Future work will include the estimate of the size and location of the chaotic regions along with the corresponding Lyapunov time, not only in case of overlapping resonances but also in the neighborhood of the hyperbolic invariant manifolds associated with the hyperbolic equilibrium points. Also, specific practical applications for both bounded and escape trajectories for the six resonant terms will be investigated. For instance, as shown in Colombo et al. (2012), the equilibrium points corresponding to orbits with apogee pointing the Sun can be exploited for enhanced Earth observation in the visible wavelength, while the equilibrium point corresponding to orbits with perigee pointing the Sun can be exploited for exploration mission of the Earth magnetosphere (McInnes et al. 2001; Oyama et al. 2008). The design of disposal trajectories with a solar sail, so far performed analytically only for planar orbit (Lücking et al. 2012) and numerically for inclined orbits (Colombo and de Bras de Fer 2016), can be fully solved analytically with the method presented here. In the field of space debris, similarly to what was done in Schaus et al. (2019) for the case of the “resonant reentry corridors”, it will be of interest to analyze the space debris catalog to see whether there exists a specific case trapped in one of the resonances found. Other applications can be considered for different planetary systems, such as the frozen solutions found in Scheeres (2012).

Finally, as mentioned at the beginning, it will be our priority to see how the phase space may change when an additional perturbation comes into play. In the case of the Earth, lunisolar gravitational perturbations and the ellipticity of the Earth will be considered for high-altitude orbits (see, e.g., Casanova et al. 2015; Gkolias and Colombo 2019; Valk et al. 2009) and the role of the atmospheric drag in the final deorbiting phase (see, e.g., Vilhena de Moraes 1981).