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

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.


Introduction
The main objective of this work is to derive an analytical model in the threedimensional 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 hyperbolic points can have several applications in planetary dynamics [1,2] and in mission design [3,4].
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 [5] and Cook [6]. 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 [7] 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 [8,9] 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 stability 1 .
In terms of applications, in the vast literature on the dynamics of dust in planetary systems, we can mention Krivov et al. [2] and Hamilton and Krivov [10], 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 [11] applied the same model, but in the planar case and only considering one of the resonant terms, 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. [12] 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 [12,13], or the use of heliotropic orbits at an oblate asteroid [14], or also the use of anti-heliotropic orbits for geomagnetic tail exploration missions [15,16,17]. Following [11] a particular application was proposed by Lücking at al. [18,19] 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 [18,19,20] 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 which can support a reentry from LEO at the end-of-life, Alessi et al. [21] 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 [22], a numerical fre-quency portrait of the LEO region was presented, highlighting also the role that high-order resonances associated with SRP [7] might have.
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.

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., [2]).
In the following, (a, e, i, Ω, ω) 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, µ the Earth gravitational parameter, J 2 the second zonal term of the geopotential, r ⊕ the equatorial radius of the Earth, C ⊕ = µJ 2 r 2 ⊕ , the obliquity of the ecliptic, P the solar radiation pressure at 1 AU, c R the reflectivity coefficient 2 , A/m the area-to-mass ratio of the small body, C SRP = 3 2 P c R A m , λ S the longitude of the Sun measured on the ecliptic plane, and 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 [6]. Under these hypotheses, following the description given in [2,12,13,23], the singly-averaged equations of motion of the spacecraft can be written as 2 We always assume c R = 1 in this work. Table 1 Argument ψ j = n 1 Ω + n 2 ω + n 3 λ S of the periodic component in terms of n 1 , n 2 , n 3 . In the table, it is showed also the corresponding notation in the work by Cook [6]. where T 1 = cos 2 2 cos 2 i 2 , andΩ If we assume that the dynamics is driven by a single term j at a time, that means that only one periodic component sin ψ j (for e, i; cos ψ j for Ω, ω) affects the motion at a time, then we can simplify the description of the dynamics by analyzing the following system where n S is the apparent mean motion of the Sun around the Earth, andΩ (J 2 ,j) andω (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, de dt j , will appear asė |j . According to [24], the behavior in inclination can be recovered at any time by means of the integral of motion which can be inverted to get 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 Λ (which depends on the dominant j−resonant term).

Resonances and Equilibrium Points
A resonance takes place when the following condition is satisfieḋ Note that Eq. (8) depends on (a, e, i, ψ j ), that is, also on C SRP , except if it is computed at ψ j = π/2 or ψ j = 3/2π. In the latter case, the resonant condition can be computed only considering the effect due to J 2 on (Ω, ω), that is, given (a, e, n 1 , n 2 , n 3 ), by solving the following quadratic equation in cos i where These are the resonant conditions depicted by Cook [6] in his Fig. 5 and remarked in [21,22,23].
On the other hand, an equilibrium point is computed if both equations in Eqs. (5) cancel out. This must happen in particular at ψ j = 0 or ψ j = π. For j = 1, 2, 5, 6, discarding the singularity at e = 0 and e = 1, the equation which gives the value of inclination corresponding to the location of the equilibrium points at ψ j = 0 and at ψ j = π, 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. 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 ψ j = 0 for j = 1, 2, 5, 6. In the table, 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 ψ j = π for j = 1, 2, 5, 6. In the table, 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 ψ j = 0 for j = 3 and j = 4. In the table, 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 ψ j = π for j = 3 and j = 4. In the table, The stability of the points can be evaluated by computing the eigenvalues of the matrix at the given equilibrium point, where Note that the partial derivatives appearing in the third equation of Eqs. (12) can be evaluated by taking into account Eq. (7), namely, and ∂T j ∂i and ∂ ∂e ∂T j ∂i can be derived from these expressions by applying the halfangle trigonometric formulae to Eqs. (3).

Hamiltonian formulation
For completeness, we provide also the Hamiltonian formulation of the dynamics under study. To this end, we can apply the argument developed in [24] 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 [24] written in terms of classical Delaunay variables, we can formulate the Hamiltonian of the problem in action angle variables (Σ, σ). We have 3 3 In our formulation, n 1 and n 2 are inverted with respect to the formulation in [24] and we call Σ 1,2,3 what they call Λ 1,2,3 .
In this way, we obtain where the Keplerian part is the part associated with J 2 is and the one due to SRP is where S j depends on the resonance, namely, being ci = n 1 n 2 Σ 1 +Σ 2 n 2 2 Σ 1 . This is also to say that the equations of motion can be also written aṡ Note that the Hamiltonian function is a constant of motion and so it is Σ 3 , since H does not depend on σ 3 .

Dynamical Description
In Figs. 1-4, we show the location of the equilibrium points at ψ j = 0 and ψ j = π, 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 areaaugmentation device with the present technology [25] and a high area-to-mass ratio debris fragment, respectively. The solutions presented for j = 1 are consistent with the ones found in [12] for the planar problem and in [13] for the inclined case  considering only J 2 . Notice how the inclination changes as the semi-major axis increases, increasing or decreasing depending on j and on whether the equilibrium is at ψ j = 0 or ψ j = π. 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 ψ 1 = 0 and ψ 2 = π of the debris fragment, for which the effect can be considerable up to Geosynchronous Earth Fig. 3 The same as in Fig. 1, but for j = 5. The dynamics under consideration is organized according to the equilibrium points which exist for a given (a, Λ, 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 Λ. 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 Λ, as a function of e, at ψ j = 0 and ψ j = π. 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Λ ≡ Λ/ √ µ computed at the equilibrium points ψ 1 = 0 and ψ 1 = π 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Λ instead of Λ for the sake of clarity to avoid large numbers. In blue, the stable equilibrium points at ψ 1 = 0, in red the unstable ones. In green the stable equilibrium points at ψ 1 = π, 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 ∈ [0 : 1], i ∈ [4 : 64] deg. Note the narrow range of inclination, i ∈ [39.8 : 40.8] deg, within which the equilibrium point at ψ 1 = 0 can be unstable. Moreover, up toΛ −20.55 km 1/2 there exists only one equilibrium point at ψ 1 = 0, which is stable (see Fig. 9 top left); from Λ −20.55 km 1/2 toΛ −20.48 km 1/2 there are 2 stable equilibrium points at ψ 1 = 0 and one unstable equilibrium point at ψ 1 = 0 (see Fig. 9 top right, note also the homoclinic connection); fromΛ −20.48 km 1/2 toΛ −20.44 km 1/2 there exists one stable equilibrium point at ψ 1 = 0 (see Fig. 9 bottom left); for values ofΛ higher thanΛ −20.44 km 1/2 , there exist one stable equilibrium point at ψ 1 = 0 and one stable and one unstable equilibrium points at ψ 1 = π (see Fig. 9 bottom right). At the transitions, namely, atΛ −20.55 km 1/2 ,Λ −20.48 km 1/2 andΛ −20.44 km 1/2 there is a bifurcation: on the horizontal lines labeled as 2, for increasingΛ, 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Λ. In other words, the inclination corresponding to a given phase space portrait, such as the ones in Fig. 9  to Fig. 8. A preliminary definition of this domain was considered in [26], 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,Λ) curves and the phase space portrait associated with 5 equilibrium points for a = 12078 km. Moreover, we note that (e,Λ) curves can be not continuous, because Eq. (7) sets specific constraints on the range ofΛ (cos i ∈ [−1, 1]).
Looking to the (e,Λ) 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Λ, the equilibrium point existing at the lowest eccentricity can disappear, and there persist one unstable equilibrium at ψ 1 = 0 and one stable equilibrium at ψ 1 = π corresponding to a same eccentricity, and one stable equilibrium at ψ 1 = 0 and one unstable equilibrium at ψ 1 = π 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 ψ 2 = 0 is higher than the curve corresponding to ψ 2 = π, 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 ψ 2 = π and one stable and one unstable at ψ 2 = 0, that is, as in Fig. 9 on the bottom right, but the curves are shifted by π in ψ. At a = 9878 km, it can appear also one unstable equilibrium at ψ 2 = π and the 3 equilibrium points can be as before or all associated to ψ 2 = π, as in Fig. 9 on the top right, but again shifted by π in ψ. 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 π in ψ. 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 ψ 5 = 0 and one unstable at ψ 5 = π, 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 ψ 3 = 0 at a low eccentricity, and one unstable at ψ 3 = 0 and one stable at ψ 3 = π. 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 ψ 4 = 0 at a low eccentricity, one unstable at ψ 4 = 0 and one stable at ψ 4 = π, and then one stable at ψ 4 = 0 and one unstable at ψ 4 = π 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.

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., [27,28]) 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 ±iν. The closer the libration curve to the elliptic point, the closer its period to T = 2π/ν. 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 ψ j = π/2 and ψ j = 3π/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 Λ (orΛ), 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 [29], by simply looking for the following conditioṅ which assumes that the maximum occurs at either ψ j = 0 or ψ j = π.
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 ψ = 0, we can compute the eccentricity corresponding to ψ = π, by inverting Eq. (14) as a function of Σ 1 , or √ 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) 4 . 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. On the other hand, along an invariant curve of the stable ψ=π unstable ψ=π Fig. 12 For a = 10078 km, the curves showed represent the location of the equilibrium points, as a function of inclination and eccentricity for the six resonant terms. Blue: stable point at ψ j = 0; red: unstable point at ψ j = 0. Green: stable point at ψ j = π; ochre: unstable point at ψ j = π. Left: A/m = 1 m 2 /kg; right: A/m = 20 m 2 /kg. phase space (recall Figs. 9-10) it is possible to estimate the excursion in eccentricity using Eqs. (20) (and an analogous one to compute the minimum eccentricity achieved), and then apply Eq. (7) to see the corresponding excursion in inclination [22]. 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. Also, specific practical applications for both bounded and escape trajectories for the six resonant terms will be investigated. For instance, as shown in [12], 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 [15,16]. The design of disposal trajectories with a solar sail, so far performed analytically only for planar orbit [18] and numerically for inclined orbits [20], can be fully solved analytically with the method presented here.