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

TL;DR

This paper develops an analytical model to analyze the equilibrium points and phase space of satellite dynamics influenced by planetary oblateness and solar radiation pressure, aiding understanding of orbital stability.

Contribution

It introduces a novel analytical framework combining planetary oblateness and solar radiation effects in a phase space model for satellite dynamics.

Findings

Derived equilibrium points and stability conditions in eccentricity-resonant angle space.

Provided Hamiltonian formulation for the coupled effects.

Applied model to Earth's case with visualization tools.

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.

Tables

Table 1. Table 1: Argument ψ j = n 1 ​ Ω + n 2 ​ ω + n 3 ​ λ S subscript 𝜓 𝑗 subscript 𝑛 1 Ω subscript 𝑛 2 𝜔 subscript 𝑛 3 subscript 𝜆 𝑆 \psi_{j}=n_{1}\Omega+n_{2}\omega+n_{3}\lambda_{S} of the periodic component in terms of n 1 , n 2 , n 3 subscript 𝑛 1 subscript 𝑛 2 subscript 𝑛 3 n_{1},n_{2},n_{3} . In the table, it is showed also the corresponding notation in the work by Cook gC62 .

$j$	$n_1$	$n_2$	$n_3$	Cook
1	1	1	-1	8
2	1	-1	-1	7
3	0	1	-1	15
4	0	1	1	14
5	1	1	1	6
6	1	-1	1	9

Table 2. Table 2: Coefficients of the quadratic equation c 1 ​ cos 2 ⁡ i + c 2 ​ cos ⁡ i + c 3 = 0 subscript 𝑐 1 superscript 2 𝑖 subscript 𝑐 2 𝑖 subscript 𝑐 3 0 c_{1}\cos^{2}{i}+c_{2}\cos{i}+c_{3}=0 associated with the equilibrium point at ψ j = 0 subscript 𝜓 𝑗 0 \psi_{j}=0 for j = 1 , 2 , 5 , 6 𝑗 1 2 5 6 j=1,2,5,6 . In the table, C ⊕ = μ ​ J 2 ​ r ⊕ 2 subscript 𝐶 direct-sum 𝜇 subscript 𝐽 2 subscript superscript 𝑟 2 direct-sum C_{\oplus}=\mu J_{2}r^{2}_{\oplus} , C S ​ R ​ P = 3 2 ​ P ​ c R ​ A m subscript 𝐶 𝑆 𝑅 𝑃 3 2 𝑃 subscript 𝑐 𝑅 𝐴 𝑚 C_{SRP}=\frac{3}{2}Pc_{R}\frac{A}{m} , β = 1 − e 2 𝛽 1 superscript 𝑒 2 \beta=\sqrt{1-e^{2}} , γ = cos 2 ⁡ ( ϵ 2 ) 𝛾 superscript 2 italic-ϵ 2 \gamma=\cos^{2}\left(\frac{\epsilon}{2}\right) , ρ = sin 2 ⁡ ( ϵ 2 ) 𝜌 superscript 2 italic-ϵ 2 \rho=\sin^{2}\left(\frac{\epsilon}{2}\right) .

$j$	$c_1$	$c_2$	$c_3$
1	$15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta +2C_{SRP}{a}^4{\beta }^4\gamma$	$2C_{SRP}{a}^4{\beta }^4\gamma \left(1-2e^2\right)-e\beta \left(3C_{\oplus }+4a^5{\beta }^{4}n{n}_S\right)$
2	$-15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta +2C_{SRP}{a}^4{\beta }^4\gamma$	$-2C_{SRP}{a}^4{\beta }^4\gamma \left(1-2e^2\right)+e\beta \left(3C_{\oplus }-4a^5{\beta }^{4}n{n}_S\right)$
5	$15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta +2C_{SRP}{a}^4{\beta }^4\rho$	$2C_{SRP}{a}^4{\beta }^4\rho \left(1-2e^2\right)-e\beta \left(3C_{\oplus }-4a^5{\beta }^{4}n{n}_S\right)$
6	$-15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta +2C_{SRP}{a}^4{\beta }^4\rho$	$-2C_{SRP}{a}^4{\beta }^4\rho \left(1-2e^2\right)+e\beta \left(3C_{\oplus }+4a^5{\beta }^{4}n{n}_S\right)$

Table 3. Table 3: Coefficients of the quadratic equation c 1 ​ cos 2 ⁡ i + c 2 ​ cos ⁡ i + c 3 = 0 subscript 𝑐 1 superscript 2 𝑖 subscript 𝑐 2 𝑖 subscript 𝑐 3 0 c_{1}\cos^{2}{i}+c_{2}\cos{i}+c_{3}=0 associated with the equilibrium point at ψ j = π subscript 𝜓 𝑗 𝜋 \psi_{j}=\pi for j = 1 , 2 , 5 , 6 𝑗 1 2 5 6 j=1,2,5,6 . In the table, C ⊕ = μ ​ J 2 ​ r ⊕ 2 subscript 𝐶 direct-sum 𝜇 subscript 𝐽 2 subscript superscript 𝑟 2 direct-sum C_{\oplus}=\mu J_{2}r^{2}_{\oplus} , C S ​ R ​ P = 3 2 ​ P ​ c R ​ A m subscript 𝐶 𝑆 𝑅 𝑃 3 2 𝑃 subscript 𝑐 𝑅 𝐴 𝑚 C_{SRP}=\frac{3}{2}Pc_{R}\frac{A}{m} , β = 1 − e 2 𝛽 1 superscript 𝑒 2 \beta=\sqrt{1-e^{2}} , γ = cos 2 ⁡ ( ϵ 2 ) 𝛾 superscript 2 italic-ϵ 2 \gamma=\cos^{2}\left(\frac{\epsilon}{2}\right) , ρ = sin 2 ⁡ ( ϵ 2 ) 𝜌 superscript 2 italic-ϵ 2 \rho=\sin^{2}\left(\frac{\epsilon}{2}\right) .

$j$	$c_1$	$c_2$	$c_3$
1	$15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta -2C_{SRP}{a}^4{\beta }^4\gamma$	$-2C_{SRP}{a}^4{\beta }^4\gamma \left(1-2e^2\right)-e\beta \left(3C_{\oplus }+4a^5{\beta }^{4}n{n}_S\right)$
2	$-15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta -2C_{SRP}{a}^4{\beta }^4\gamma$	$2C_{SRP}{a}^4{\beta }^4\gamma \left(1-2e^2\right)+e\beta \left(3C_{\oplus }-4a^5{\beta }^{4}n{n}_S\right)$
5	$15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta -2C_{SRP}{a}^4{\beta }^4\rho$	$-2C_{SRP}{a}^4{\beta }^4\rho \left(1-2e^2\right)-e\beta \left(3C_{\oplus }-4a^5{\beta }^{4}n{n}_S\right)$
6	$-15C_{\oplus }e\beta$	$-6C_{\oplus }e\beta -2C_{SRP}{a}^4{\beta }^4\rho$	$2C_{SRP}{a}^4{\beta }^4\rho \left(1-2e^2\right)+e\beta \left(3C_{\oplus }+4a^5{\beta }^{4}n{n}_S\right)$

Table 4. Table 4: Coefficients of the cubic equation sin 3 ⁡ i + s 1 ​ sin 2 ⁡ i + s 2 ​ sin ⁡ i + s 3 = 0 superscript 3 𝑖 subscript 𝑠 1 superscript 2 𝑖 subscript 𝑠 2 𝑖 subscript 𝑠 3 0 \sin^{3}{i}+s_{1}\sin^{2}{i}+s_{2}\sin{i}+s_{3}=0 associated with the equilibrium point at ψ j = 0 subscript 𝜓 𝑗 0 \psi_{j}=0 for j = 3 𝑗 3 j=3 and j = 4 𝑗 4 j=4 . In the table, C ⊕ = μ ​ J 2 ​ r ⊕ 2 subscript 𝐶 direct-sum 𝜇 subscript 𝐽 2 subscript superscript 𝑟 2 direct-sum C_{\oplus}=\mu J_{2}r^{2}_{\oplus} , C S ​ R ​ P = 3 2 ​ P ​ c R ​ A m subscript 𝐶 𝑆 𝑅 𝑃 3 2 𝑃 subscript 𝑐 𝑅 𝐴 𝑚 C_{SRP}=\frac{3}{2}Pc_{R}\frac{A}{m} , β = 1 − e 2 𝛽 1 superscript 𝑒 2 \beta=\sqrt{1-e^{2}} .

$j$	$s_1$	$s_2$	$s_3$
3	$-\frac{2C_{SRP}}{15C_{\oplus }e}{a}^4{\beta }^3\sin ϵ$	$-\frac{4}{15C_{\oplus }}\left(3C_{\oplus }-a^5{\beta }^{4}n{n}_s\right)$	$\frac{2C_{SRP}}{15C_{\oplus }}e{a}^4{\beta }^3\sin ϵ$
4	$\frac{2C_{SRP}}{15C_{\oplus }e}{a}^4{\beta }^3\sin ϵ$	$-\frac{4}{15C_{\oplus }}\left(3C_{\oplus }+a^5{\beta }^{4}n{n}_s\right)$	$-\frac{2C_{SRP}}{15C_{\oplus }}e{a}^4{\beta }^3\sin ϵ$

Table 5. Table 5: Coefficients of the cubic equation sin 3 ⁡ i + s 1 ​ sin 2 ⁡ i + s 2 ​ sin ⁡ i + s 3 = 0 superscript 3 𝑖 subscript 𝑠 1 superscript 2 𝑖 subscript 𝑠 2 𝑖 subscript 𝑠 3 0 \sin^{3}{i}+s_{1}\sin^{2}{i}+s_{2}\sin{i}+s_{3}=0 associated with the equilibrium point at ψ j = π subscript 𝜓 𝑗 𝜋 \psi_{j}=\pi for j = 3 𝑗 3 j=3 and j = 4 𝑗 4 j=4 . In the table, C ⊕ = μ ​ J 2 ​ r ⊕ 2 subscript 𝐶 direct-sum 𝜇 subscript 𝐽 2 subscript superscript 𝑟 2 direct-sum C_{\oplus}=\mu J_{2}r^{2}_{\oplus} , C S ​ R ​ P = 3 2 ​ P ​ c R ​ A m subscript 𝐶 𝑆 𝑅 𝑃 3 2 𝑃 subscript 𝑐 𝑅 𝐴 𝑚 C_{SRP}=\frac{3}{2}Pc_{R}\frac{A}{m} , β = 1 − e 2 𝛽 1 superscript 𝑒 2 \beta=\sqrt{1-e^{2}} .

$j$	$s_1$	$s_2$	$s_3$
3	$\frac{2C_{SRP}}{15C_{\oplus }e}{a}^4{\beta }^3\sin ϵ$	$-\frac{4}{15C_{\oplus }}\left(3C_{\oplus }-a^5{\beta }^{4}n{n}_s\right)$	$-\frac{2C_{SRP}}{15C_{\oplus }}e{a}^4{\beta }^3\sin ϵ$
4	$-\frac{2C_{SRP}}{15C_{\oplus }e}{a}^4{\beta }^3\sin ϵ$	$-\frac{4}{15C_{\oplus }}\left(3C_{\oplus }+a^5{\beta }^{4}n{n}_s\right)$	$\frac{2C_{SRP}}{15C_{\oplus }}e{a}^4{\beta }^3\sin ϵ$

Equations

$${\psi }_j=n_1\Omega +n_2\omega +n_3{\lambda }_S,{\psi }_j\in [0,2\pi ],$$

$$\frac{da}{dt}$$

$$\frac{de}{dt}$$

$$\frac{di}{dt}$$

$$\frac{d\Omega }{dt}$$

$$\frac{d\omega }{dt}$$

$${\mathcal{T}}_1$$

$${\mathcal{T}}_2$$

$${\mathcal{T}}_3$$

$${\mathcal{T}}_4$$

$${\mathcal{T}}_5$$

$${\mathcal{T}}_6$$

$${\dot{\Omega }}_{J_2}$$

$${\dot{\Omega }}_{SRP}$$

$${\dot{\omega }}_{J_2}$$

$${\dot{\omega }}_{SRP}$$

$${\frac{de}{dt}|}_j$$

$$\frac{d{\psi }_j}{dt}$$

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

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

$${\dot{\psi }}_j=0.$$

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

$$c_1$$

$$c_2$$

$$c_3$$

$$J=\left(\begin{array}{cc}\frac{\partial \dot{e}|j}{\partial e}& \frac{\partial \dot{e}|j}{\partial {\psi }_j}\\ \frac{\partial {\dot{\psi }}_j}{\partial e}& \frac{\partial {\dot{\psi }}_j}{\partial {\psi }_j}\end{array}\right),$$

$$\frac{\partial {\dot{e}}_{|j}}{\partial e}$$

$$\frac{\partial {\dot{e}}_{|j}}{\partial {\psi }_j}$$

$$\frac{\partial {\dot{\psi }}_j}{\partial e}$$

$$\frac{\partial {\dot{\psi }}_j}{\partial {\psi }_j}$$

$$\frac{\partial {\dot{\Omega }}_{(J_2,j)}}{\partial e}$$

$$C_{SRP}\frac{e}{na\sqrt{1-e^2}}\frac{\partial {\mathcal{T}}_j}{\partial i}\cos {\psi }_j\frac{\partial }{\partial e}\frac{1}{\sin i}+C_{SRP}\frac{e}{na\sqrt{1-e^2}\sin i}\cos {\psi }_j\frac{\partial }{\partial e}\frac{\partial {\mathcal{T}}_j}{\partial i},$$

$$\frac{\partial {\dot{\omega }}_{(J_2,j)}}{\partial e}$$

$$\frac{15}{2}\frac{J_2{r}_{\oplus }^{2}n\cos i}{a^2(1-e^2{)}^2}\frac{\partial \cos i}{\partial e}+C_{SRP}\frac{\sqrt{1-e^2}}{nae}\cos {\psi }_j\frac{\partial {\mathcal{T}}_j}{\partial e}-$$

$$C_{SRP}\frac{e}{na\sqrt{1-e^2}}\frac{\partial {\mathcal{T}}_j}{\partial i}\cos {\psi }_j\frac{\partial }{\partial e}\frac{1}{\tan i}+C_{SRP}\frac{e}{na\sqrt{1-e^2}\tan i}\cos {\psi }_j\frac{\partial }{\partial e}\frac{\partial {\mathcal{T}}_j}{\partial i},$$

$$\frac{\partial \cos i}{\partial e}$$

$$\frac{\partial \sin i}{\partial e}$$

$${\Sigma }_1$$

$${\Sigma }_2$$

$${\Sigma }_3$$

$$\mathcal{H}={\mathcal{H}}_{J_2}+{\mathcal{H}}_{SRP},$$

$${\mathcal{H}}_{J_2}=\frac{J_2{\mu }^4{r}_{\oplus }^2\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)}{4L^3{\Sigma }_1^5{n}_2^7},$$

$${\mathcal{H}}_{SRP}=-\frac{C_{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),$$

Full text

††journal: Celestial Mechanics and Dynamical Astronomy

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., Mi ; REF9 ) and in mission design (e.g., REF14 ; REF13 ).

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 pM60 and Cook gC62 . 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 H77 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 sB99 ; sB01 addressed the analytical treatment of luni-solar 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 stability111Recall that the SRP effect can be written analogously to the solar gravitational perturbation, except for the amplitude of the perturbation and the first order of its expansion H77 ..

In terms of applications, in the vast literature on the dynamics of dust in planetary systems, we can mention Krivov et al. REF9 and Hamilton and Krivov REF10 , 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 REF11 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. REF1 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 towards the Sun) or anti-heliotropic (perigee pointing towards the Sun) orbits. The equilibrium conditions computed analytically for the planar case were numerically extended to non-zero inclinations. For such frozen orbits, different applications were proposed: the use of heliotropic orbits for Earth observation in the visible wavelength REF1 ; REF2 , or the use of heliotropic orbits at an oblate asteroid REF12 , or also the use of anti-heliotropic orbits for geomagnetic tail exploration missions REF7 ; REF6 ; REF8 . Following REF11 a particular application was proposed by Lücking at al. REF4 ; REF5 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 REF4 ; REF5 ; REF15 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. Alessi_CMDA 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_freq , a numerical frequency portrait of the LEO region was presented, highlighting also the role that high-order resonances associated with SRP H77 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 ; Scheeres ).

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 , 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 Camilla2019 and in the geosynchronous region (e.g., CasanovaEtAl ; GkoliasColombo ).

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 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., REF9 ).

In the following, $(a,e,i,\Omega,\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_{2}r^{2}_{\oplus}$, $\epsilon$ the obliquity of the ecliptic, $P$ the solar radiation pressure at 1 AU, $c_{R}$ the reflectivity coefficient222We always assume $c_{R}=1$ in this work., $A/m$ the area-to-mass ratio of the small body, $C_{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

[TABLE]

with $n_{1},n_{2},n_{3}$ are as in Table 1, where we provide also the correspondence with respect to the notation used by Cook gC62 .

Under these hypotheses, following the description given in REF9 ; REF1 ; REF2 ; Alessi_MNRAS , the singly-averaged equations of motion of the spacecraft can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

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 $\Omega,\omega$) affects the motion at a time, then we can simplify the description of the dynamics by analyzing the following system

[TABLE]

where $n_{S}$ is the apparent mean motion of the Sun around the Earth, and $\dot{\Omega}_{(J_{2},j)}$ and $\dot{\omega}_{(J_{2},j)}$ are intended to include the contribution of $J_{2}$ and the $j$th-contribution of the SRP. In the following, the time derivative of the eccentricity due to the periodic term $j$, namely, $\left.\frac{de}{dt}\right\rvert_{j}$, will appear as $\dot{e}_{|j}$.

**Since we can apply the same argument used in D16 for the coupled lunisolar gravitational perturbation and oblateness effect (see Alessi_MNRAS and Sec. 2.2), ** the behavior in inclination can be recovered at any time by means of the integral of motion

[TABLE]

which can be inverted to get

[TABLE]

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 $\Lambda$ (which depends on the dominant $j-$resonant term).

2.1 Resonances and Equilibrium Points

A resonance takes place when the following condition is satisfied

[TABLE]

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.

Note that Eq. (8) depends on $(a,e,i,\psi_{j})$, that is, also on $C_{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 $(\Omega,\omega)$, that is, given $(a,e,n_{1},n_{2},n_{3})$, by solving the following quadratic equation in $\cos i$

[TABLE]

where

[TABLE]

These are the resonant conditions depicted by Cook gC62 in his Fig. 5 and remarked in Alessi_CMDA ; Schettino_freq ; Alessi_MNRAS .

On the other hand, an equilibrium point is computed if both equations in Eqs. (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 showed in Tabs. 2-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 Tabs. 4-5, we show the coefficients of such cubic equations.

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

[TABLE]

at the given equilibrium point, where

[TABLE]

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

[TABLE]

where

[TABLE]

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 Eqs. (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 REF9 ; REF1 ; REF2 ; Alessi_MNRAS . Notice that the disturbing potential associated with the SRP also corresponds to Eq. (17) in CasanovaEtAl and to Eq. (69) in Scheeres . Then, we can apply the argument developed in D16 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 D16 written in terms of classical Delaunay variables $\left(L,G,H,l,g,h\right)$, by introducing the canonical variables $\left(\Gamma,\tau\right)$ such that $\dot{\tau}\equiv{\partial{\mathcal{H}}\over\partial\Gamma}=\lambda_{S}$, we can formulate the Hamiltonian of the problem in action angle variables $(\Sigma,\sigma)$. We have333In our formulation, $n_{1}$ and $n_{2}$ are inverted with respect to the formulation in D16 and we call $\Sigma_{1,2,3}$ what they call $\Lambda_{1,2,3}$.

[TABLE]

Note that $\sigma_{1}\equiv\sigma\equiv\psi_{j}$ and that Eq. (6) is derived from $\Sigma_{2}$.

In this way, we obtain

[TABLE]

where the part associated with $J_{2}$ is

[TABLE]

and the one due to SRP is

[TABLE]

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

[TABLE]

being $\text{ci}=\frac{n_{1}n_{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

[TABLE]

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 sB99 for lunisolar perturbations, in Eq. (16) $\Sigma_{3}$ is a dummy variable and its value is not important.

3 Dynamical Description

In Figs. 1-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$ m${}^{2}/$kg, $A/m=1$ m${}^{2}/$kg and $A/m=20$ m${}^{2}/$kg, characteristic of a standard satellite, a satellite equipped with an area-augmentation device with the present technology REF3 and a high area-to-mass ratio debris fragment, respectively.

The solutions corresponding to inclination equal to zero, for $j=1$ and $j=2$, are equivalent to the ones described in REF1 . For $A/m=0.012$ m${}^{2}/$kg and $j=1$ or $j=2$, the solutions presented in REF2 , 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 $(\Omega,\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$ 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 showed is computed by considering one of the solutions of the cubic equations displayed in Tabs. 4-5. The other solutions correspond to the opposite value showed and to the singularity $i=0.$

The dynamics under consideration is organized according to the equilibrium points which exist for a given $(a,\Lambda,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 $\Lambda$. 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 $\Lambda$, 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 Eqs. (2)), while the dominant resonant term $j$ might change along a curve, as it will be showed in the following.

The first example is presented in Fig. 7. For $j=1$, $A/m=1$ m${}^{2}/$kg and $a=8078$ km, the value of $\tilde{\Lambda}\equiv\Lambda/\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{\Lambda}$ instead of $\Lambda$ for the sake of clarity to avoid large numbers. In blue, the stable equilibrium points at $\psi_{1}=0$, in red the unstable ones. In green the stable equilibrium points at $\psi_{1}=\pi$, in ochre 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 showed, but displaying 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$ m${}^{2}/$kg and $a=8078$ 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{\Lambda}\simeq-20.55$ km1/2 there exists only one equilibrium point at $\psi_{1}=0$, which is stable (see Fig. 9 top left); from $\tilde{\Lambda}\simeq-20.55$ km1/2 to $\tilde{\Lambda}\simeq-20.48$ km1/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{\Lambda}\simeq-20.48$ km1/2 to $\tilde{\Lambda}\simeq-20.44$ km1/2 there exists one stable equilibrium point at $\psi_{1}=0$ (see Fig. 9 bottom left); for values of $\tilde{\Lambda}$ higher than $\tilde{\Lambda}\simeq-20.44$ km1/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{\Lambda}\simeq-20.55$ km1/2, $\tilde{\Lambda}\simeq-20.48$ km1/2 and $\tilde{\Lambda}\simeq-20.44$ km1/2 there is a bifurcation: on the horizontal lines, for increasing $\tilde{\Lambda}$, 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.

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{\Lambda}$. In other words, the inclination corresponding to a given phase space portrait, such as the ones in Fig. 9, can be inferred by looking to Fig. 8. A preliminary definition of this domain was considered in SchettinoIAC , taking into account only the strictly resonant behaviors for almost zero eccentricities.

For $j=1$, $A/m=1$ m${}^{2}/$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$ km. In Fig. 10, we show the corresponding behavior of the $(e,\tilde{\Lambda})$ curves and the phase space portrait associated with 5 equilibrium points for $a=12078$ km. Moreover, we note that $(e,\tilde{\Lambda})$ curves can be not continuous, because Eq. (7) sets specific constraints on the range of $\tilde{\Lambda}$ ($\cos{i}\in[-1,1]$).

Looking to the $(e,\tilde{\Lambda})$ curves for $j=1$, always considering the solution for $\cos{i}$ on the first column of Fig. 1, for $A/m=20$ m${}^{2}/$kg it is noticed a similar behavior, while for $A/m=0.012$ m${}^{2}/$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$ m${}^{2}/$kg and given $\tilde{\Lambda}$, 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, it is stressed the following. 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$ m${}^{2}/$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$ km, it can appear also one unstable equilibrium at $\psi_{2}=\pi$ and the 3 equilibrium points can be as before or all associated to $\psi_{2}=\pi$, as in Fig. 9 on the top right, but again shifted by $\pi$ in $\psi$. From about $a=10078$ 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$ m${}^{2}/$kg, the same behavior is found, while for $A/m=0.012$ m${}^{2}/$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 towards $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$ m${}^{2}/$kg and $A/m=1$ m${}^{2}/$kg, this is always the case from about $a=19000$ 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$ km.

Resonant term #3. For $A/m=1$ m${}^{2}/$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$ km (about $a=15000$ km for $A/m=0.012$ m${}^{2}/$kg and about about $a=17000$ km for $A/m=20$ m${}^{2}/$kg), where the equilibrium at the lowest eccentricity disappears. From that value of semi-major axis on, the 2 equilibrium points drift towards 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$ m${}^{2}/$kg until about $a=10000$ km, for $A/m=1$ m${}^{2}/$kg until about $a=11600$ km, for $A/m=20$ m${}^{2}/$kg until about $a=19000$ km. From that value on, there persist 4 equilibrium points, which are associated to 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., PC ; W ) 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 showed in Figs. 8-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 $\Lambda$ (or $\tilde{\Lambda}$), to achieve a given objective, in particular to obtain a bounded eccentricity variation within a given timespan or an escape trajectory.

Moreover, from the analysis provided in Sec. 2, the maximum eccentricity corresponding to a given invariant curve can be computed, as suggested in AC_IAC , by simply looking for the following condition

[TABLE]

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-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-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)$444as done for the resonant, in the strict sense, curves, e.g., in gC62 ; Alessi_MNRAS .. In Fig. 12, we show an example for $a=10078$ km and $A/m=1$ m${}^{2}/$kg and $A/m=20$ m${}^{2}/$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 towards 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 Eqs. (19) (and an analogous one to compute the minimum eccentricity achieved), and then apply Eq. (7) to see the corresponding excursion in inclination Schettino_freq . 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.

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 REF1 , 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 REF7 ; REF6 . The design of disposal trajectories with a solar sail, so far performed analytically only for planar orbit REF4 and numerically for inclined orbits REF15 , can be fully solved analytically with the method presented here. In the field of space debris, similarly to what was done in Volker 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 .

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., CasanovaEtAl ; GkoliasColombo ; L1 ) and the role of the atmospheric drag in the final deorbiting phase (see, e.g., R1 ).

Acknowledgements.

The authors would like to acknowledge the funding received by the European Commission Horizon 2020, Framework Programme for Research and Innovation (2014-2020), under the ReDSHIFT project (grant agreement n. 687500) and the funding received by the European Research Council (ERC) through the European Commission Horizon 2020, Framework Programme for Research and Innovation (2014-2020), under the project COMPASS (grant agreement n. 679086). The authors are grateful to Josep Masdemont, Ioannis Gkolias and Kleomenis Tsiganis for the useful discussions.

Compliance with ethical standards

Conflict of interest: The authors declare that they have no conflict of interest.