Attitude and orbit coupling of planar helio-stable solar sails

TL;DR

This paper investigates the coupled attitude and orbit dynamics of planar helio-stable solar sails, simplifying the complex system through averaging techniques to analyze their behavior and stability.

Contribution

It introduces a simplified model of solar sail dynamics by averaging out fast attitude phases, focusing on orbit perturbations caused by a fixed sail orientation.

Findings

The system can be modeled as a fast-slow dynamical system.

Averaging out the fast attitude phase simplifies the analysis.

Numerical simulations validate the simplified model.

Abstract

The coupled attitude and orbit dynamics of solar sails is studied. The shape of the sail is a simplified quasi-rhombic-pyramid that provides the structure helio-stablility properties. After adimensionalisation, the system is put in the form of a fast-slow dynamical system where the different time scales are explicitely related to the physical parameters of the system. The orientation of the body frame with respect to the inertial orbit frame is a fast phase that can be averaged out. This gives rise to a simplified formulation that only consists of the orbit dynamics perturbed by a flat sail with fixed attitude perpendicular to the direction of the sunlight. The results are exemplified using numerical simulations.

Tables

$a$  =	semi-major axis, km
$\alpha$  =	aperture angle of the sail, deg or rad
$A_s$  =	area of the panels, m2
$\beta$  =	frequency of the Sun-pointing direction, deg or rad
$C$  =	inertia moment around the $\zeta$ axis, kg m2
$d$  =	center of mass-center of pressure offset, m
$e$  =	eccentricity, [-]
$\epsilon$  =	ratio between characteristic time scales, [-]
${ℱ}_b$  =	body frame
${ℱ}_I$  =	Earth centered inertial frame
$𝐅$  =	force vector, N
$\gamma$  =	difference $\omega -\Omega$, deg or rad
$h$  =	height of the panels, m
$𝐢$  =	unit vectors of a basis
$J_2$  =	adimensional $J_2$ coefficient, [-]
$\lambda$  =	angle between Sun position and ${𝐢}_x$, rad or deg
$m_s/2$  =	mass of the panels, kg
$m_b$  =	mass of the bus, kg
$𝐌,M$  =	torque vector, component of torque vector, Nm
$\mu$  =	earth’s mass parameter, km3/s2
$𝐧$  =	normal vectors to panels
$n$  =	mean motion, deg or rad/s
$\Omega$  =	right ascension of the ascending node, deg or rad
$\omega$  =	argument of the perigee, deg or rad
$P$  =	sail panel
$p_{\mathrm{SR}}$  =	solar radiation pressure at 1 AU, N/m2
$r$  =	magnitude of the position vector of the spacecraft, km
$R$  =	Earth’s radius, km
$\Sigma$  =	surface of section
$\theta$  =	true anomaly, deg or rad
$𝐮$  =	unit vector
$v$  =	velocity components and its magnitude, km/s
$\phi$  =	Euler angle, rad or deg
$\Phi$  =	angular velocity of the attitude, rad or deg/s
$\varphi$  =	difference $\phi -\lambda$, rad or deg
$w$  =	width of the panels, m
$\xi ,\nu ,\zeta$  =	coordinates of ${ℱ}_b$
$x,y,z$  =	coordinates of ${ℱ}_I$
$X$  =	characteristic function of an interval
Subscripts
$\pm$  =	that refers to panels $+$ or $-$
$\odot$  =	that refers to Sun
$\mathrm{SRP}$  =	that refers to solar radiation pressure
$\mathrm{GG}$  =	that refers to gravity gradient
$\xi ,\nu ,\zeta$  =	in the direction of, in ${ℱ}_b$
$x,y,z$  =	in the direction of, in ${ℱ}_I$

Table 2. Table 2: Values of the physical parameters for α = 35 ∘ 𝛼 superscript 35 \alpha=35^{\circ} and α = 40 ∘ 𝛼 superscript 40 \alpha=40^{\circ} .

	$\alpha ={35}^{\circ },d=0$ m	$\alpha ={40}^{\circ },d=0$ m
$c_1$	$4.133317062536305\times {10}^2$	$5.747509656406245\times {10}^2$
$c_2$	$2.014647115843597\times {10}^0$	$1.923989341570575\times {10}^0$
$c_2$	$1.650597476175750\times {10}^{-4}$	$1.650597476175750\times {10}^{-4}$
$c_2$	$3.738547970136426\times {10}^{-6}$	$3.738547970136426\times {10}^{-6}$
$\epsilon$	$4.918703449585804\times {10}^{-2}$	$4.171191657263433\times {10}^{-2}$
Time unit	$2.203569462524180\times {10}^2$ s	$1.868685651104933\times {10}^2$ s

Table 3. Table 3: Values of the physical parameters for α = 45 ∘ 𝛼 superscript 45 \alpha=45^{\circ} and α = 60 ∘ 𝛼 superscript 60 \alpha=60^{\circ} .

	$\alpha ={45}^{\circ },d=0$ m	$\alpha ={60}^{\circ },d=0$ m
$c_1$	$7.624959636935995\times {10}^2$	$1.366246396170031\times {10}^3$
$c_2$	$1.811184377377631\times {10}^0$	$1.297157388066479\times {10}^0$
$c_2$	$1.650597476175750\times {10}^{-4}$	$1.650597476175750\times {10}^{-4}$
$c_2$	$3.738547970136426\times {10}^{-6}$	$3.738547970136426\times {10}^{-6}$
$\epsilon$	$3.621439426788271\times {10}^{-2}$	$2.705424915355282\times {10}^{-2}$
Time unit	$1.622397734086550\times {10}^2$ s	$1.212025036217823\times {10}^2$ s

Equations

$$R_{3}(\psi)=\left(\begin{array}[]{rr}\cos\psi&\sin\psi\\ -\sin\psi&\cos\psi\end{array}\right).$$

$$C\ddot{\phi }$$

$$C$$

$${\mathbf{F}}_{\mathrm{SRP}}^{\pm }$$

$$\ddot{\phi }$$

$$M_1(\psi )$$

$$M_0^{\pm }(\psi )$$

$$X_{J}:\mathbb{R}\to\{0,1\},\qquad X_{J}(\psi)=\left\{\begin{array}[]{rcl}1&\mbox{if}&\psi\in J\\ 0&\mbox{if}&\psi\notin J\\ \end{array}\right..$$

$$k_{1,1}$$

$$k_{2,0}$$

$$k_{0,2}$$

$$\displaystyle\left\{\begin{array}[]{rcl}\ddot{x}&=&\displaystyle-\frac{\mu x}{r^{3}}-\frac{3R^{2}\mu J_{2}}{2}\frac{x}{r^{5}}+\frac{A_{s}p_{\rm SR}}{m_{b}+m_{s}}a_{x},\\ \ddot{y}&=&\displaystyle-\frac{\mu y}{r^{3}}-\frac{3R^{2}\mu J_{2}}{2}\frac{y}{r^{5}}+\frac{A_{s}p_{\rm SR}}{m_{b}+m_{s}}a_{y},\end{array}\right.$$

$$(a_x,a_y{)}^{\top }$$

$$\displaystyle\begin{array}[]{rcl}a_{x}^{\pm}&=&\displaystyle\sin(\alpha\pm\lambda\mp\varphi)(\eta\cos(2\alpha\pm\lambda\mp 2\varphi)-\cos\lambda))\\ a_{y}^{\pm}&=&\displaystyle\sin(\alpha\pm\lambda\mp\varphi)(\mp\eta\sin(2\alpha\pm\lambda\mp 2\varphi)-\sin\lambda))\end{array}$$

$${\mathbf{n}}_{\pm }=R_3(\phi )(\sin \alpha ,\pm \cos \alpha {)}^{\top },{\mathbf{u}}_{\odot }=(\cos \lambda ,\sin \lambda {)}^{\top }.$$

$$-2\frac{A_s{p}_{\mathrm{SR}}}{m_b+m_s}(1+\eta )(\cos \lambda ,\sin \lambda {)}^{\top },$$

$$\displaystyle\left\{\begin{array}[]{rcl}\ddot{\phi}&=&\displaystyle c_{1}M_{1}(\phi)+\frac{c_{2}}{r^{3}}\sin(2\arctan(y/x)-2(\phi+\lambda)),\\ \ddot{x}&=&\displaystyle-\frac{x}{r^{3}}-c_{3}\frac{x}{r^{5}}+c_{4}a_{x},\\ \ddot{y}&=&\displaystyle-\frac{y}{r^{3}}-c_{3}\frac{y}{r^{5}}+c_{4}a_{y},\end{array}\right.$$

$$c_1=\frac{A_s}{m_b+m_s}\frac{p_{\mathrm{SR}}{k}_{1,1}{L}^3}{2C\mu },c_2=3\frac{D(\alpha ,d)}{C},c_3=\frac{3R^2{J}_2}{2L^2},c_4=\frac{A_s}{m_b+m_s}\frac{p_{\mathrm{SR}}{L}^2}{\mu }.$$

$$\displaystyle\left\{\begin{array}[]{rcl}\varepsilon\frac{{\rm d}\hat{\bm{\phi}}}{{\rm d}\tau}&=&f(\hat{\bm{\phi}},\hat{\bm{x}},\varepsilon)\\ \frac{{\rm d}\hat{\bm{x}}}{{\rm d}\tau}&=&g(\hat{\bm{\phi}},\hat{\bm{x}})\end{array}\right.,$$

$${\epsilon }^2$$

$$\displaystyle\left\{\begin{array}[]{rcl}\frac{{\rm d}\hat{\bm{\phi}}}{{\rm d}\tau}&=&f(\hat{\bm{\phi}},\hat{\bm{x}},\varepsilon)\\ \frac{{\rm d}\hat{\bm{x}}}{{\rm d}\tau}&=&\varepsilon g(\hat{\bm{\phi}},\hat{\bm{x}})\end{array}\right..$$

$$f(\hat{\mathbf{\varphi }},\hat{\mathbf{x}},\epsilon )$$

$$M_0(\psi )$$

$$M_1(\hat{\varphi })$$

$$g(\hat{\mathbf{\varphi }},\hat{\mathbf{x}})$$

$$({\hat{a}}_x,{\hat{a}}_y{)}^{\top }$$

$${\hat{a}}_{\hat{x}}^{\pm }$$

$${\hat{a}}_{\hat{y}}^{\pm }$$

$$\mathcal{H}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\hat{\mathbf{\varphi }}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}(\hat{\varphi },\hat{\Phi }{)}^{\top }=f(\hat{\mathbf{\varphi }},\hat{\mathbf{x}},0)=(\hat{\Phi },M_1(\hat{\varphi }){)}^{\top }={\left(\frac{\partial {\mathcal{K}}_0}{\partial \hat{\Phi }},-\frac{\partial {\mathcal{K}}_0}{\partial \hat{\varphi }}\right)}^{\top }.$$

$$K_0(\psi )$$

Full text

1. 1 Introduction

2. 2 Model

3. 2.1 Geometry of the sail structure

4. 2.2 Equations of motion

5. 2.2.1 SRP model

6. 2.2.2 Attitude dynamics

7. 2.2.3 Orbit dynamics

8. 2.3 Multiple time scale dynamics

9. 2.3.1 Fast equations

10. 2.3.2 Slow equations

11. 3 Dynamical aspects of the coupled system

12. 3.1 Dynamics of the fast subsystem

13. 3.2 Fast-slow dynamics and averaging

14. 3.3 Explicit equations and physical interpretation

15. 3.3.1 Stability of the Sun-pointing direction

16. 3.3.2 Averaged equations of motion

17. 3.3.3 Interpretation as an equivalent flat sail

18. 4 Numerical test cases

19. 4.1 Physical parameters and adimensionalisation factors

20. 4.2 Main numerical experiment

21. 4.3 Numerical results

22. 4.3.1 Shape, size and orientation of the osculating orbit

23. 4.3.2 The definition of equivalent effective area

24. 4.3.3 Full system versus averaged system

25. 5 Conclusions

26. A Adimensional model

Nomenclature

1 Introduction

Solar sails are a low-thrust propulsion system that takes the advantage of Solar Radiation Pressure (SRP) to accelerate a probe using a highly reflective surface. Even though the acceleration due to SRP is smaller than that achieved by traditional thrusters, it is continuous and enables possibilities that range from trajectory design to attitude control.

The effects of the SRP acceleration have been widely studied in the literature. In works that deal with the effect of SRP acceleration on the dynamics of spacecraft with a solar sail, its attitude is usually assumed to be fixed with respect to the direction of sunlight. This assumption simplifies the equations of motion and, in some cases, allows to deal with SRP acceleration as a perturbative effect and even to write the system as if it had Hamiltonian structure [13]. But maintaining the attitude fixed with respect to any direction, in particular that of sunlight, requires attitude control.

For the specific case of works that deal with flat sails whose surface is (theoretically) assumed to be always perpendicular to the sunlight direction, it is expected that in practice the attitude oscillates close to this direction if an appropriate center of Mass - center of Pressure Offset (that in this paper is referred to as MPO) is assumed.

In this work we investigate the possibility of considering a different structure that consist of a number of flat panels oriented in a way that cancels out some torque components and makes the Sun-pointing attitude stable. The motivation is to foster the use of such sails as passive deorbiting devices, that can be employed as end-of-life disposals that would reduce the attitude control requirements. Here the term passive is understood in the sense of the so-called passive deorbiting strategy, as defined in [5], that consists of using the idea of deorbiting “outwards” on an elliptical orbit: the increase of the eccentricity of the orbit causes the perigee radius to progressively decrease. As justified in [15, 16], this can be attained by orienting the sail panel always perpendicular to the sunlight direction. Other deorbiting strategies are the so-called “active” approaches (as opposed to “passive”) that consist of changing between maximal and minimal SRP acceleration along the motion. This was first studied in [2], where the authors suggested maximizing the SRP acceleration when travelling towards the Sun and minimizing it when travelling away from the Sun. A refined version of this approach is to maximize (resp. minimize) the SRP acceleration when the first averaged variational equation of the eccentricity is positive (resp. negative), see [7]. This refinement reduces the number of attitude change maneuvers from twice per turn around the Earth to twice per year, see [5].

Apart from considering an adequate MPO, helio-stability can be enhanced by means of a Quasi-Rhombic Pyramid (QRP) shape. This idea was first introduced in [3]. The suggested structure consists of 4 reflective panels that resemble the shape of a pyramid. In case the center line of the pyramid is oriented close enough to the sunlight direction and has an adequate MPO, this structure cancels out, on average along the motion, the components of the acceleration in other directions. For example, in [10] the authors study the linear stability of the Sun-pointing attitude; and this stability can be further enhanced by assuming a moderate spin around the adequate axis of inertia as proposed in [11].

A simplified version of the QRP that consists of a single triangular flat panel and an appropriately positioned payload of the spacecraft was considered in [4], and later exploited in [12] to design new periodic orbits in the circular restricted three body problem. The suggested spacecraft was shown to have undamped conservative oscillatory dynamics around the Sun-pointing direction.

Despite the contributions [3, 4, 10, 11, 12] provide satisfactory results, there is, to the authors’ knowledge, a lack of understanding of

The attitude dynamics, especially close to the Sun-pointing attitude, and 2. 2.

The attitude and orbit coupling: especially whether they can be analytically separated taking into account the fact that these two components have two characteristic time scales.

These two questions are addressed in this paper by considering a sail structure in between the single-panel considered in [4, 12] and the full QRP [3, 10, 11] whose orbit dynamics (without any attitude control) evolves strictly on the ecliptic plane if adequate initial attitude conditions are chosen. The structure consists of two panels with variable aperture and variable position of the payload, as introduced in [20]. The acceleration due to SRP assumes that photons are partially specularly reflected and partially absorbed. Building on previous contributions on the usage of the SRP effect for the design of end-of-life disposals (see, e.g., [15, 6]), here the considered orbital dynamics are the $J_{2}$-perturbed two-body problem111This is the motion of an artificial satellite around an oblate planet keeping only the $J_{2}$ term of the expansion of the perturbing potential in spherical harmonics. always perturbed by the SRP acceleration, that is, the spacecraft is considered to be always illuminated so the effect of eclipses are not taken into account. The SRP acceleration depends on the attitude of the spacecraft. The attitude dynamics are assumed to solely happen around an axis perpendicular to the ecliptic plane, and to be perturbed by the SRP and gravity-gradient torques. The effect of atmospheric drag is not taken into account. Hence, this study is relevant to higher Low Earth Orbits (LEO) (i.e. with altitude 700/800 km and above).

As numerically demonstrated in [20], this structure has the advantage that under some hypotheses related to the geometry of the sail -that are discussed later in this contribution-, the dynamics close to the Sun-pointing direction is close to a mathematical pendulum and hence the motion has a quasi-integral (adiabatic invariant) of motion that is almost preserved over some time interval. Also, the length of this time interval depends on the ratio between time scales as usually described by theorems concerning the accuracy of the averaging method.

The organization and presentation of the main results of this paper are as follows. First of all, § 2 is devoted to the review of the geometry of the spacecraft under consideration and to the derivation of the equations of motion. Despite having two characteristic time scales, the equations are not written in the form of fast-slow systems. The equations are put as a fast-slow system of differential equations, where the variables that evolve in different time scales are splitted, and related via a physical parameter that represents the ratio between time scales that only depends on the geometry of the spacecraft. In § 3, the dynamics of the system are studied in the context of fast-slow systems, and this includes the discussion of the possibilities of the separation of the motions. The system obtained by direct averaging of the fast phase (after adequate changes of variables) is related to the results of the averaging theorems. The section finishes with an enumeration of the physical interpretations of the results. The theoretical results and formulas of § 3 are tested in § 4 with special emphasis on the physical interpretations just mentioned. The paper concludes in § 5, where the main results of the contribution are summarized and different possible lines for future research are suggested.

2 Model

This section is devoted to providing the equations of motion of the planar dynamics of a helio-stable solar sail. These are a set of differential equations that govern the coupled attitude and orbit dynamics of the spacecraft under consideration. The content of § 2.1 and § 2.2 is a summary of the derivation of the equations of motion that is added for completeness. For further details the reader is referred to [20]. The section ends by putting the equations of motion in the context of dynamical systems with multiple time scales in § 2.3.

2.1 Geometry of the sail structure

The spacecraft under consideration consist of a payload or bus attached to two panels forming an angle. To avoid out-of-plane motion, one is lead to consider a simplification of a QRP [3] that consists of two panels of equal size $P_{\pm}$; of height $h$, width $w$, and area $A_{s}=hw$. Assume that the mass of each panel is $m_{s}/2$, so the mass of the whole sail structure is $m_{s}$. In the left panel of Fig. 1 a sketch of the sail structure is depicted.

The attitude dynamics of the spacecraft occurs in a reference frame $\mathcal{F}_{b}$ attached to it. The coordinates in this frame are referred to as $\xi,\nu$ and $\zeta$ and the vectors of the basis ${\bf i}_{\xi},{\bf i}_{\nu},{\bf i}_{\zeta}$. The panels are attached to each other along an $h$-long side, that lies on a line parallel to the $\zeta$ axis, and they form an angle $\alpha$ with respect to the plane $\nu=0$. The payload, whose mass is denoted as $m_{b}$, is assumed to be on the $\xi$ axis, at a distance $d$ from the center of mass of the two panels, see Fig. 2. The parametrization of the panels is chosen so that the center of mass of the spacecraft is at the origin of $\mathcal{F}_{b}$. The main physical parameters of the system are: $\alpha$, the aperture angle; and $d$, that accounts for the MPO.

Sketches of top views of the spacecraft in $\mathcal{F}_{b}$ can be seen in Fig. 2, where the bus is depicted as a black square, and the center of mass of the sail structure is depicted as a blue solid dot, added to visualize the parameter $d$. The left, center and right panels are sketches of spacecraft with $d<0$, $d=0$ and $d>0$, respectively.

The aim of this contribution is to study the oscillatory attitude dynamics close to the Sun-pointing direction of the spacecraft described in this subsection. For the purposes of this article and to simplify the exposition we considered that the back part of the panels (the side where the angle $\alpha$ is measured in Fig. 2) did not produce any SRP acceleration. In [19] the attitude dynamics model is extended to take into account the effect of the back side neglected here, and the dynamics close to the Sun-pointing attitude is shown to be exactly the same as the one obtained in this paper. Moreover, a numerical study of the consequences of considering different reflectance properties is performed in [19].

2.2 Equations of motion

The considered planar orbit and attitude dynamics, are a coupled system of differential equations in $(\mathbb{S}^{1}\times\mathbb{R})\times\mathbb{R}^{4}$, where $\mathbb{S}^{1}:=\mathbb{R}/(2\pi\mathbb{Z})$: orientation and angular velocity for the attitude; and position and velocity of the spacecraft in an Earth centered inertial reference frame $\mathcal{F}_{I}$.

Here SRP is considered to be the coupling effect between the attitude and orbit dynamics. It is then necessary to study the attitude dynamics in relation to the orbit dynamics that are considered to evolve in $\mathcal{F}_{I}$. The coordinates of $\mathcal{F}_{I}$ are denoted $x,y$ and $z$, and the vectors of the orthonormal basis are denoted ${\bf i}_{x,y,z}$. The vector ${\bf i}_{x}$ points towards an arbitrarily chosen direction on the ecliptic (e.g. J2000), and since we are dealing with the planar problem, the vector ${\bf i}_{z}$ is parallel to ${\bf i}_{\zeta}$, and they are also perpendicular to the ecliptic plane. The triad is completed by choosing ${\bf i}_{y}={\bf i}_{z}\times{\bf i}_{x}$.

As the motion is planar the change of coordinates from $\mathcal{F}_{I}$ to $\mathcal{F}_{b}$ is done through $R_{3}(-\varphi)$, where $\varphi\in\mathbb{S}^{1}$ is an Euler angle and $R_{3}$ is the rotation matrix around the $z$ (and $\zeta$) axis. The rotation matrix $R_{3}$ reads, for any angle $\psi\in[0,2\pi)$,

[TABLE]

The Euler attitude equations in the present situation reduce to

[TABLE]

where ${\bf M}=(M_{\xi},M_{\nu},M_{\zeta})$ is the torque due to the external forces considered, and $C$ is the inertia moment around the $\zeta$ axis in $\mathcal{F}_{b}$. Denote $I_{\xi,b},I_{\nu,b},I_{\zeta,b}$ as the inertia moments of the bus. Then one can see that

[TABLE]

2.2.1 SRP model

The force due to SRP exerted on each panel of the sail in $\mathcal{F}_{b}$ is considered to be modelled as [17]

[TABLE]

where ${\bf u}_{\odot}$ is the unit vector in the Earth-Sun direction, and ${\bf n}_{\pm}$ are the normal vectors to each panel, recall Fig. 1. Concerning the constants, $\eta\in(0,1)$ is the (dimensionless) reflectance of the sail and $p_{\rm SR}=4.56\times 10^{-6}\;{\rm N}/{\rm m}^{2}$ is the solar pressure at 1 AU which is considered to be constant.

2.2.2 Attitude dynamics

The effects taken into consideration are SRP and the non-symmetry of the spacecraft, so the total torque is ${\bf M}={\bf M}_{\rm SRP}+{\bf M}_{\rm GG}$, the sum of the SRP and gravity gradient torques. Let $\lambda$ be the argument of latitude of the apparent position of the Sun. The SRP torque has a different representation depending on the orientation of the sail structure with respect to the Sun, that is, it depends on the value of $\varphi$ relative to $\lambda$, so denote $\phi=\varphi-\lambda$. These three angles are sketched in Fig. 3.

If $\phi\in[-\pi+\alpha,\alpha]$, the panel $P_{+}$ produces torque, see Fig. 4, (a); and if $\phi\in[-\alpha,\pi-\alpha]$, it is the panel $P_{-}$ who produces torque, see Fig. 4, (b); in particular, if $|\phi|\leq\alpha$, both panels face the Sun, see Fig. 4 (c). In all panels of Fig. 4, the sunlight direction is depicted as if it was in the direction $(-1,0)^{\top}$, hence $\lambda=180^{\circ}$, see Fig. 3.

Assume that the bus is symmetric in the sense that $I_{\xi,b}=I_{\nu,b}=I_{\zeta,b}$. In this case, as derived in [20], the attitude equations of motion reduce to the following second order Ordinary Differential Equation (ODE)

[TABLE]

where $D(\alpha,d)$ is defined in Eq. 2, $\theta$, $\omega$ and $\Omega$ are the true anomaly, argument of perigee and Right Ascension of the Ascending Node (RAAN) of the osculating orbit; the last is being considered as it precesses due to the $J_{2}$ effect considered in § 2.2.3, and

[TABLE]

where $X_{J}$ is the characteristic function of the interval $J$,

[TABLE]

The parameter $\mu=GM=3.986\times 10^{14}\;{\rm m}^{3}/{\rm s}^{2}$ is the gravitational parameter of the Earth. The rest of the coefficients are physical parameters that depend on the geometry of the spacecraft and on the reflectance parameter $\eta$, and read

[TABLE]

The function $M_{0}$ in Eq. 5b can be interpreted as the scaled SRP torque due to a single panel. Since both panels are equal but oriented in a different way, one expects the expressions for both panels to be similar. Namely these can be found to differ only by a sign. The function $M_{1}$ in Eq. 5a, on the other hand, represents the joint SRP torque, that takes in consideration that the panels have to face the sunlight to produce torque; this is why $M_{1}$ is a piece-wise defined function.

2.2.3 Orbit dynamics

Inspired by previous studies on the effect of SRP for the design of end-of-life disposals, the considered orbit dynamics has been the two-body problem, perturbed by the $J_{2}$ term and the SRP acceleration due to the sail structure [15, 6, 16, 5]. The (dimensional) equations of motion of the orbit dynamics read

[TABLE]

where the first summand is the Keplerian term the second is the $J_{2}$ effect and the third one is the SRP acceleration. The constant $J_{2}=1.082\times 10^{-3}$ is the adimensional coefficient of the second order term in the expansion of the perturbing potential in spherical harmonics, $R$ is the radius of the planet, and $r=\sqrt{x^{2}+y^{2}}$.

It is important to remark that the factors $a_{x}$ and $a_{y}$ in the third summands of the right hand side of Eq. 9 are piecewise defined, as the SRP acceleration depends on $\varphi-\lambda$. Namely, $a_{x}$ and $a_{y}$ can be expressed as follows:

[TABLE]

where

[TABLE]

are the expressions of the adimensional factor of the accelerations due to panel $P_{+}$ (subscript $+$) and due to panel $P_{-}$ (subscript $-$). Equations 13 are obtained from Eq. 3, with

[TABLE]

In the particular case that $\alpha=\pi/2$ and $\varphi=\lambda$, that is, when the sail is completely opened and hence it is a rectangular flat panel with width $2w$ and height $h$, and the direction of the normal to the surface of the sail is parallel to the Sun-spacecraft direction, the SRP acceleration reads

[TABLE]

twice the acceleration of a flat panel of area $A_{s}=hw$ always oriented towards the Sun, maximizing the SRP acceleration.

The assumption that the Sun-spacecraft distance is constant and equal to the Sun-Earth distance is equivalent to considering a linear approximation of the potential of the Sun [13].

2.3 Multiple time scale dynamics

The full set of 6 coupled differential equations, Eq. 4 and Eq. 9, have two different time scales, attitude being faster than orbit dynamics. Hence these equations fit within the context of fast-slow dynamical systems. In this section we choose adequate variables to put them in a standard form for their treatment as a fast-slow system.

Appendix A is devoted to the adimensionalisation of the equations of motion, where the adimensionalisation factors $L$ (longitude) and $T$ (time) are introduced. For the sake of lightening the notation, we use the same notation for the adimensional variables, as the dimensional analogues will not be further used in this contribution. Also $(\dot{\,})$ denotes the derivative with respect to the adimensional time, $\tau$. The full set of equations read

[TABLE]

where $\phi=\varphi-\lambda$ (depending on adimensional time), $a_{x}$ and $a_{y}$ are as in Eq. 10, and the constants are the adimensional quantities

[TABLE]

Written like this, the problem can be put in the form of a fast-slow system. To write Eq. 17 as an ODE of first order, one must introduce $\Phi=\dot{\phi}=\dot{\varphi}-\dot{\lambda}$, and $v_{x}=\dot{x},v_{y}=\dot{y}$. Let us denote $\bm{\phi}=(\phi,\Phi)$ and $\bm{x}=(x,y,v_{x},v_{y})$ 222This notation will also be used for other sets of variables that are going to be introduced later: $\hat{\bm{\phi}}=(\hat{\phi},\hat{\Phi})$, $\bar{\bm{\phi}}=(\bar{\phi},\bar{\Phi})$, $\tilde{\bm{\phi}}=(\tilde{\phi},\tilde{\Phi})$, $\tilde{\tilde{\bm{\phi}}}=(\tilde{\tilde{\phi}},\tilde{\tilde{\Phi}})$, $\hat{\bm{x}}=(\hat{x},\hat{y},\hat{v}_{\hat{x}},\hat{v}_{\hat{y}})$, $\bar{\bm{x}}=(\bar{x},\bar{y},\bar{v}_{\bar{x}},\bar{v}_{\bar{y}})$. $\tilde{\bm{x}}=(\tilde{x},\tilde{y},\tilde{v}_{\tilde{x}},\tilde{v}_{\tilde{y}})$, $\tilde{\tilde{\bm{x}}}=(\tilde{\tilde{x}},\tilde{\tilde{y}},\tilde{\tilde{v}}_{\tilde{\tilde{x}}},\tilde{\tilde{v}}_{\tilde{\tilde{y}}})$. . Then the following holds.

Proposition 1**.**

Assume that $\dot{\lambda}=n_{\odot}={\rm constant}$, and that the torque coefficient $k_{1,1}>0$ (Eq. 6a). There exists $\varepsilon=\varepsilon(c_{1})>0$ and a phase scaling $(\bm{\phi},\bm{x})\to(\hat{\bm{\phi}},\hat{\bm{x}})$ such that in the $(\hat{\bm{\phi}},\hat{\bm{x}})$ variables Eq. 17 has the following form

[TABLE]

for adequate $f$ and $g$.

Proof. If we define

[TABLE]

the proposition follows by considering the scaling $\hat{\phi}=\phi,\quad\hat{\Phi}=\varepsilon\Phi,\quad\hat{x}=x,\quad\hat{y}=y,\quad\hat{v}_{x}=v_{x},\quad\hat{v}_{y}=v_{y}.$ The maps $f$ and $g$ correspond to the right hand side of the resulting equations of motion for $\hat{\bm{\phi}}$ and $\hat{\bm{x}}$, respectively. $\Box$

The parameter $\varepsilon$ in Eq. 22 depends solely on physical quantities of the system and on the choice of the longitude scaling factor $L$ in the adimensionalisation procedure, see App. A. If $L$ is appropriately chosen, see § 4.1, then $\varepsilon$ is small, and Eq. 21 is a fast-slow system, written in the slow time scale $\tau$. One can consider a time scaling $t=\tau/\varepsilon$, the fast time scale, in which Eq. 21 reads

[TABLE]

Hence, $\varepsilon$ is the ratio between the time scales in which the two characteristic motions, orbit and attitude, take place.

The change to $\hat{\cdot}$ variables can be extended to the argument of latitude of the Sun $\hat{\lambda}=\lambda$ which, in the fast time scale $t$, varies as ${\rm d}\hat{\lambda}/{\rm d}t=\varepsilon n_{\odot}$, similarly as for the orbital dynamics in Eq. 21. That is, from the point of view of the attitude, the position of the Sun varies in a slower constant rate, with factor $\varepsilon$.

2.3.1 Fast equations

The components of the vector field in Eq. 26 of the fast variables are those representing the evolution of the sail attitude

[TABLE]

where, if we denote

[TABLE]

we can write

[TABLE]

To obtain these expressions one has to use the double angle formulas in Eq. 5, using the symmetries of the involved trigonometric functions, arranged as in Eq. 5a.

2.3.2 Slow equations

Concerning the components of the vector field in Eq. 26 of the slow variables, the scaling of the proof of Prop. 1 is such that $\hat{\bm{x}}=\bm{x}$, and hence the form of the equations does not change. The vector field of the slow subsystem reads

[TABLE]

where $\hat{r}=\sqrt{\hat{x}^{2}+\hat{y}^{2}}$ and

[TABLE]

being

[TABLE]

Written as in Eq. 21 and Eq. 25, the equations of motion can be studied by separately considering only the attitude or the orbit dynamics. This is possible by dealing with the limit case $\varepsilon=0$, that can be interpreted as the spacecraft having an infinitely large area-to-mass ratio. Setting $\varepsilon=0$ has two different meanings in each equivalent formulation, Eq. 21 and Eq. 25:

• •

In the slow time scale, Eq. 21, the dynamics are the slow (orbit) system, constrained to the zeros of the function $f$. In practice, in the context of this paper, this means that the dynamics is constrained to a specific attitude. In fact, in this case the dynamics is that of a time-dependent Hamiltonian, and hence the components of $g$ in Eq. 41 can be obtained as the derivatives of

[TABLE]

The time dependency comes from the SRP acceleration, as it depends on the position of the Sun.

• •

In the fast time scale, Eq. 25, the position of the spacecraft is assumed to be constant and only the attitude evolves. The position of the spacecraft is, hence, a parameter of the system.

Moreover, for $\varepsilon=0$ both the slow and fast time scale dynamics are Hamiltonian: in the slow time scale this holds for any fixed attitude, and in the fast time scale it is not trivial and is justified in § 3.1. However, the whole problem is not Hamiltonian due to the SRP and gravity gradient coupling.

3 Dynamical aspects of the coupled system

In this section we highlight the most dynamically relevant aspects of the coupled adimensional attitude and orbit model. Written as a fast-slow problem, we can deal with the description of the dynamics as it is customary in this field: after splitting the equations of motion into fast and slow components as in Eq. 25, written in the fast time scale, the parameter $\varepsilon$ is set to [math], leaving only as non-trivial the equations that correspond to the fast dynamics. These are usually referred to as the fast subsystem, or frozen system. The dimension of the resulting problem is smaller than the original, and hence easier to study. The core idea is that, as there are two time scales, for small values of $\varepsilon$ the variables that are frozen for $\varepsilon=0$ evolve slower than the fast variables, so the dynamics of the full fast-slow system is expected to be close to the dynamics of the fast subsystem. More concretely, one can explain the dynamics of the full system by proving that some invariant manifolds of the frozen system are preserved when considering $\varepsilon>0$, via the so-called Geometric Singular Perturbation Theory, based on Fenichel’s theory on the preservation of normally hyperbolic invariant manifolds under perturbation. The expository text [14] is strongly suggested for a global overview on the field.

In the context of this contribution, the stability of the Sun-pointing attitude in the fast subsystem can be translated into the oscillatory motion of the spacecraft close to the Sun-pointing attitude, that is expected to be preserved for long periods of time. Moreover, the fast dynamics can be averaged out.

First, § 3.1 is devoted to the study of the dynamics of the fast subsystem assuming that the slow subsystem is frozen. After this, the averaging of the fast small oscillations of the sail around the Sun-pointing direction is studied in § 3.2. This allows to obtain a physical interpretation of the results that are exposed in § 3.3.

3.1 Dynamics of the fast subsystem

Consider Eq. 25, the full system written in the fast time scale $t$. By setting $\varepsilon=0$ the orbit (slow) vector field vanishes so the only nontrivial equations are those corresponding to the attitude (fast) dynamics, that read

[TABLE]

recall Eq. 26. This is not equivalent to the simplified model found and studied in [20], where the orbit dynamics was assumed to happen on a fixed Keplerian orbit. Recall also that the SRP acceleration due to the back part of the panels is neglected, as this effect is not relevant for the purposes of this paper.

The first relevant property of the fast subsystem is that it has Hamiltonian structure.

Proposition 2**.**

The system given in Eq. 44 for $\varepsilon=0$ is Hamiltonian with some Hamiltonian function $\mathcal{K}_{0}$ and hence can be written as

[TABLE]

Proof. If Eq. 44 was a Hamiltonian system with Hamiltonian $\mathcal{K}_{0}$, the equations of motion would be obtained as derivatives of $\mathcal{K}_{0}$ as indicated in Eq. 45. As the first component of $f$ is $\hat{\Phi}$, $\mathcal{K}_{0}$ must have $\hat{\Phi}^{2}/2$ as summand; and as the second component of $f$ is $M_{1}$, another summand of $\mathcal{K}_{0}$ would be an appropriate primitive of $M_{1}$, see Eq. 28. The function $M_{1}$ is $\mathcal{C}^{0}$, so one has to make sure that this primitive is $\mathcal{C}^{1}$.

If we denote

[TABLE]

then a solution is

[TABLE]

$\Box$

Remark 1**.**

For $|\hat{\phi}|\leq\alpha$, $\mathcal{K}_{0}=\hat{\Phi}^{2}/2-\cos(2\hat{\phi})/2-2K_{0}(-\alpha)$ is the Hamiltonian function of a pendulum. 2. 2.

For $\varepsilon>0$, the attitude equations of motion can also be written as the derivatives of a function $\mathcal{K}_{\varepsilon}$ with respect to $\hat{\phi}$ and $\hat{\Phi}$. Namely, if we extend Eq. 48

[TABLE]

then for $\varepsilon>0$ we recover the vector field Eq. 26 via

[TABLE]

The most relevant properties of the fast vector field $f$ in Eq. 45 are: it is $\mathcal{C}^{0}$, the differentiability being lost at the switching manifolds $\hat{\phi}=\pm\pi\mp\alpha,\pm\alpha$; it is $2\pi$-periodic in $\hat{\phi}$; and it is also symmetric with respect to $\hat{\phi}=0$ as $\mathcal{K}_{\varepsilon}((\hat{\phi},\hat{\Phi}),\hat{\bm{x}})=\mathcal{K}_{\varepsilon}((-\hat{\phi},\hat{\Phi}),\hat{\bm{x}})$. An example of the phase space can be seen in the left panel of Fig. 5, for $\alpha=30^{\circ}$ and $d=0$. In this Figure, each curve represents an orbit of the fast subsystem Eq. 44 that is obtained as a level set of $\mathcal{K}_{0}$, Eq. 48 (that is, points for which $\mathcal{K}_{0}={\rm constant}$). The origin $E$ represents the Sun-pointing attitude and the orbits around it represent oscillatory motion. The vertical dashed lines represent switching manifolds, that in this case represent physically that a panel either starts or ceases to face sunlight and hence starts or ceases to produce torque. For $|\hat{\phi}|<30^{\circ}$ both panels face sunlight, for $30^{\circ}<|\hat{\phi}|<180^{\circ}-30^{\circ}$ only one of them do (as explained in § 2.2.2) and for $|\hat{\phi}|>180^{\circ}-30^{\circ}$ the motion is completely rotational, i.e. the spacecraft tumbles, as no panel is assumed to produce torque.

The set of equilibria are an isolated point $E$ and a continuum: at all points $(\hat{\phi},0)$ with $\hat{\phi}\in I_{-}\cup\{0\}\cup I_{+}$, where $I_{-}=[-\pi,-\pi+\alpha]$ and $I_{+}=[\pi-\alpha,\pi)$ the vector field vanishes. Among these the most dynamically relevant are $E=(0,0)$, which is stable (provided $k_{1,1}>0$ in the original coordinates, see Eq. 6), and $H_{\pm}=(\pm\pi \mp\alpha,0)$, that are saddles whose invariant manifolds coincide, $W^{u}(H_{+})=W^{s}(H_{-})$ and $W^{u}(H_{-})=W^{s}(H_{+})$. The equilibria are indicated in the left panel of Fig. 5, and in the right panel the physical meaning of $E,H_{\pm}$ is sketched: $E$ is the Sun-pointing attitude and $H_{\pm}$ represent the angles of transition to from no reflective panel facing sunlight to one, or vice-versa. The rest of equilibria, those whose abscissa is $\hat{\phi}\in I_{-}\cup I_{+}\setminus\{-\pi+\alpha,\pi-\alpha\}$ have 0 as double eigenvalue.

3.2 Fast-slow dynamics and averaging

This section is devoted to the study of the dynamics of the full system for $\varepsilon>0$. For small values of $\varepsilon$ one expects the structure found in § 3.1 to be close to conserved in a sense that is made explicit here. More concretely, the function $\mathcal{K}_{\varepsilon}$ is a first integral of the frozen fast subsystem for $\varepsilon=0$, so for $\varepsilon>0$ and for adequate initial conditions, it is expected to vary, but slowly, along orbits. Hence, one expects the system to admit an analogue of an adiabatic invariant (recall that the full system in Eq. 21 is not Hamiltonian): there is an equivalent formulation of the system in which one of the variables experiences oscillations of at most $\mathcal{O}(\varepsilon)$ for time intervals of length $\mathcal{O}(1/\varepsilon)$. The goal of this section is to find the adequate change of variables that translates our problem in this context.

Attitudes close to the Sun-pointing direction $\hat{\phi}=0$ are those for which there is numerical evidence of helio-stability properties of the spacecraft, see [20], and the dynamics in this regime is that of a mathematical pendulum, recall item $1$ in Rem. 1. The following study is restricted to this situation.

For the full set of equations of motion written in the fast time scale in Eq. 25, let $G$ be a bounded subset of the phase space such that, for all $(\hat{\bm{\phi}},\hat{\bm{x}})\in G$, $|\hat{\phi}|<\alpha$ (that is, both panels face the sunlight) where all trajectories lie wholly in $G$. Note that, in $G$, for fixed values of $\hat{\bm{x}}$, the condition $\mathcal{K}_{\varepsilon}=k$ defines one and only one trajectory of Eq. 25, and the restriction of $G$ onto the $(\hat{\phi},\hat{\Phi})$ variables consists of a continuum of closed and nested periodic orbits.

Lemma 1**.**

For small enough $\varepsilon>0$ and small enough $\hat{\Phi}$ (that is, close enough to the Sun-pointing attitude) there exists a real analytic change of variables defined in the interior of $G$

[TABLE]

that transforms Eq. 25, with $f$ and $g$ as given in Eq. 26 and Eq. 41, respectively, into a system of the form

[TABLE]

where $\tilde{f},\tilde{g}$ have zero average with respect to $\tilde{\phi}$.

Proof. This result fits within the scope of averaging theory, and the proof can be sketched as consisting of the two following steps.

Consider for the moment the dynamics of the fast subsystem Eq. 45. It depends on $\hat{\phi}$ in a periodic way, but for each initial condition in $G$, the period and the range of $\hat{\phi}$ of points in this orbit are different. Restricting to the behaviour close to the Sun-pointing direction, we can substitute the equations of motion by the Taylor series around $\hat{\phi}=0$. Since the frozen subsystem is Hamiltonian, one can have $\mathcal{K}_{0}$ in mind for the moment. The leading non-constant terms are

[TABLE]

This suggests to consider a first change of variables ($\hat{\cdot}\to\bar{\cdot}$)

[TABLE]

where $(\bar{\Phi},\bar{\phi})$ are the usual Poincaré action-angle variables

[TABLE]

the frequency being $\beta=\sqrt{2}$. The change is extended to the rest of variables by simply choosing $\bar{x}=\hat{x}$, $\bar{y}=\hat{y}$, $\bar{v}_{\bar{x}}=\hat{v}_{\hat{x}}$, and $\bar{v}_{\bar{y}}=\hat{v}_{\hat{y}}$.

After this change of variables, the obtained equations of motion are $2\pi$ periodic in $\bar{\phi}$, which is still fast with respect to the rest of the variables.

It is worth noting that as the fast subsystem has Hamiltonian structure and the change of variables $C_{1}$ restricted to the attitude motion is canonical, one can still deal with these two equations as if they were a Hamiltonian system with Hamiltonian function $\mathcal{K}_{\varepsilon}$. So, for orbits in $G$, we can obtain the equations of motion in the $\bar{\cdot}$ variables as follows. Start by considering

[TABLE]

where $\vartheta:=\arctan(\hat{y}/\hat{x})-\hat{\lambda}$. Now, the expansion around $\hat{\phi}=0$ reads

[TABLE]

where, if we introduce the $\bar{\cdot}$ variables Eq. 58 with $\beta=\sqrt{2}$, it reads

[TABLE]

Now the attitude equations in $\bar{\cdot}$ variables can be recovered via

[TABLE]

Concerning the orbital motion (or slow subsystem), it only depends on $\hat{\phi}$ via the term due to SRP acceleration, see the rightmost summand in Eq. 41. In the region $G$, this acceleration (see Eq. 42) reads

[TABLE]

This can be treated as we did for $\mathcal{K}_{\varepsilon}$ above, first separating the dependence on $\hat{\phi}$, substituting the sine and cosine functions by their Taylor expansions, and then introducing the $\bar{\cdot}$ coordinates. These expressions are not added as only a part of them are useful for the next step, and we can refer to the expansions above to justify the form they have. 2. 2.

Let us focus for the moment on the differential equation of $\bar{\Phi}$, that is, after performing the first change of variables $C_{1}$. This equation is obtained as $-\partial\mathcal{K}_{\varepsilon}/\partial\bar{\phi}$. Note that all summands except $\beta\bar{\Phi}$ in the right hand side of Eqs. 60 have a factor that depends on $\bar{\phi}$, only as an argument of a power of a sine function, so after deriving with respect to $\bar{\phi}$, the derivative of all terms appears in the right hand side of the equation of motion of $\bar{\Phi}$. Moreover, in each summand there is a factor of the form of the right hand side of

[TABLE]

that have zero average with respect to $\bar{\phi}$. Hence, one can get rid of the dependence on $\bar{\phi}$ via some averaging steps as done, for instance, in the classical reference [21]. After ordering the terms (by orders in $\bar{\Phi}$, for instance), the $j$th step would consist of two substeps. Firstly the terms of the right hand side we want to get rid of have to be separated as the sum of periodic (with respect to $\bar{\phi}$) plus average parts. Secondly, a change of variables has to be constructed in a way that the non-periodic part remains intact and the periodic targeted terms do not appear in the equations written in the new variables. This is done at the expense of more terms appearing in higher orders that have to be dealt with in subsequent averaging steps.

To be able to perform these changes we have to require that, on the one hand, the frequency $\beta>0$ has to be bounded away from zero. This is satisfied as $\beta=\sqrt{2}$. On the other hand, the leading terms of the expansion are $\mathcal{O}(\bar{\Phi}^{2})$ (see Eq. 60a), hence we have to assume that $\bar{\Phi}$ is small enough so that the successive changes of variables to be performed are, in fact, invertible.

To reach the claimed form of the equations, Eq. 55, one has to perform a second change of variables

[TABLE]

that is a composition of $N\geq 0$ averaging steps. The number of steps $N$ depends on the size of $|\bar{\Phi}|$ relative to $\varepsilon$. Note that there are terms in the expansions in Eqs. 60 that do not have $\varepsilon$ as factor, namely those in the summatory in Eq. 60a, which start with $\bar{\Phi}^{2}$. It may happen, if we are far enough from the equilibrium, that $\mathcal{O}(\bar{\Phi}^{2})>\mathcal{O}(\varepsilon^{2})$, that is, these terms dominate over those in Eqs. 60b and 60c. Let $l\geq 2,l\in\mathbb{Z}$ be the first power such that $\mathcal{O}(\bar{\Phi}^{l})<\mathcal{O}(\varepsilon^{2})$. To reach the desired form in the equations of motion in Eq. 55, we have to get rid of all terms in Eq. 60a that have $\bar{\Phi}^{j}$, $j<l$ as factor. This is done via a sequence of averaging steps, each dedicated to eliminate the lowest power appearing (see, e.g., [9]). As the first term has order $2$ in $\bar{\Phi}$, the number of averaging steps to perform is $N=l-2$. In particular, if $\mathcal{O}(\bar{\Phi}^{2})$ is already smaller than $\mathcal{O}(\varepsilon^{2})$ all terms are small enough and then $C_{2}$ is the identity.

For small enough $\bar{\Phi}$, $C_{2}$ is the composition of near-the-identity changes of variable such that, in the new variables, the lowest order term in $\bar{\Phi}$ in the Taylor expansion that depends on $\bar{\phi}$ is exchanged by a zero-average term with respect to $\bar{\phi}$ of the same order in $\bar{\Phi}$. This can be done in such a way that after the $N$ steps the differential equations consist of the average of the original equation written in the $\bar{\cdot}$ variables, at the expense of having a change of variables that does not have zero average, plus higher order terms that have $\varepsilon$ as factor.

After performing $C_{1}$ and then $C_{2}$, Eq. 21 has, in the $\tilde{\cdot}$ variables, the claimed form Eq. 55, where $\tilde{f}$ is the derivative with respect to $\bar{\Phi}$ of the average with respect to $\bar{\phi}$ of the terms in the sum Eq. 60a up to the first order comparable with $\varepsilon$ in magnitude, with $\tilde{\Phi}$ in the place of $\bar{\Phi}$; $\tilde{g}$ consists of the slow equations of motion plus the average up to the same order of Eq. 65; and the functions $R_{1}$, $R_{2}$ and the map $\bm{R}$ contain the non averaged terms plus extra terms that appear as a byproduct of the successive averaging changes of variables. $\Box$

Now we are in position of using the averaging principle, that consists of replacing the right hand side of Eq. 55 with its average with respect to $\tilde{\phi}$ and to neglect the terms $\mathcal{O}(\varepsilon^{2})$. Let us refer to the averaged variables as $\tilde{\tilde{\cdot}}$. The averaged equations of motion read

[TABLE]

where the angle brackets denote average with respect to $\tilde{\phi}$. Recall that $\tilde{\phi}$ indicates the average angle between the $\xi$ axis in $\mathcal{F}_{b}$ (i.e. the orientation of the sail) and the direction of sunlight. Consider the solutions of Eq. 55 and Eq. 66 starting at the same initial conditions, and denote them by $\tilde{\tilde{\Phi}}(t)$, $\tilde{\tilde{\bm{x}}}(t)$ and $\tilde{\Phi}(t)$, $\tilde{\bm{x}}(t)$, respectively.

Theorem 1**.**

Assume $\varepsilon>0$ is small enough. Then, on the one hand, $\tilde{\Phi}$ is an adiabatic invariant, that is,

[TABLE]

On the other hand, starting with the same initial conditions, the solution of the slow variables in Eq. 55 is described by the solution of Eq. 66 with an accuracy of $\mathcal{O}(\varepsilon)$ over time intervals of length $\mathcal{O}(1/\varepsilon)$, that is,

[TABLE]

Proof. This follows from the theorem about the accuracy of the averaging method, see [1], Chapter 10, § 52. This applies as the equations of motion are analytic, so they meet the minimum (finite) differentiability requirements, and $|\beta|=\sqrt{2}>0$ is bounded away from zero. $\Box$

Remark 2**.**

The averaging steps performed have been done in such a way that the resulting equations of motion after the successive changes of variables were the average with respect to the fast angle $\bar{\phi}$. Although this is the expected observed mean motion, it adds difficulties in the equations as it forces the successive changes to not have zero average and makes the remainders of the expression Eq. 55 more involved. If instead one performs these steps getting rid of all possible periodic terms, the changes and the remainders remain $2\pi$ periodic through the whole process. Moreover, in this situation, the classical theorem of Neishtadt [21] is applicable, so there exists a change of variables that separates the phase $\hat{\phi}$ from the rest of the variables except from a remainder that has exponentially small bounds.

3.3 Explicit equations and physical interpretation

The benefits of the analysis performed in this section is that one can extract physical interpretations of interest for prospective applications: spacecraft design guidelines, simplification of the equations of motion and the interpretation of the averaged problem as an equivalent already studied problem.

3.3.1 Stability of the Sun-pointing direction

The first practical benefit can be extracted from the study of the fast dynamics. The condition $k_{1,1}>0$, that is a hypothesis of Prop. 1 is actually equivalent to

[TABLE]

This is a necessary constraint between the physical parameters of the system, $\alpha$ and $d$. Note that since $\alpha\in[0,\pi/2]$ and $\eta\in(0,1)$, the function $K(\alpha,\eta)<0$ we have the following result.

Corollary 1**.**

[20]** For each aperture angle $\alpha\in(0,\pi/2)$ there exists $d_{\min}$ such that if $d>d_{\min}$ the origin of Eq. 45 is locally stable.

This local stability can be physically understood as when the motion starts close enough to the Sun-pointing attitude the attitude librates around this state, as seen in Fig. 5.

This result relates the aperture angle with the MPO and can be understood as a guideline for construction and can be used for future control-related studies. Note that, in particular, $d_{\min}<0$ for $\alpha\in(0,\pi/2)$ and $d_{\min}=0$ for $\alpha=\pi/2$, the flat plate case; recall Fig. 2. In Fig. 6 the quantity $d_{\min}$ is displayed as a function of $\alpha$, for different values of the reflectance coefficient $\eta$.

3.3.2 Averaged equations of motion

The previous analysis justifies the averaging procedure with respect to the oscillations around the Sun-pointing direction. The main hypotheses for the theorems were, on the one hand, the closeness to the Sun-pointing direction; and on the other hand, the smallness of $\varepsilon$, that actually measures the separation of time scales between the attitude and orbit components. This was already observed in the numerical study performed in [20].

In the proof of Lemma 1, the changes of variable were chosen to justify the applicability of the well-known averaging results, but from a practical perspective, provided $\varepsilon$ and $\bar{\Phi}$ are small, the averaging process can be carried out up to any required order. Namely, the larger $\bar{\Phi}$ is, the higher the order of the averaging procedure has to be.

Concerning the attitude equations, in practice one can obtain the differential equations by first averaging the Hamiltonian $\mathcal{K}_{\varepsilon}$, and then computing the derivatives with respect to the new averaged action. Namely, using that

[TABLE]

the average of Eq. 60 reads (using the same names for the variables)

[TABLE]

The derivative of $\left<\mathcal{K}_{\varepsilon}\right>$ with respect to $\bar{\Phi}$ gives the approximation of the average frequency of rotation of $\bar{\phi}$.

Concerning the orbital dynamics, after the averaging procedure they are approximately decoupled from the attitude. The averaged effect of the SRP acceleration is found by computing the average of Eq. 65.

Let us proceed as done from Eq. 59 to Eq. 60, that is, first separating the dependence on $\hat{\phi}$ and then expanding the sine and cosine functions around $\hat{\phi}=0$ in Eq. 65. After this procedure, the only terms that will contribute to the average are with those that have $\cos(l\hat{\phi})$, $l\in\mathbb{Z}$ as factor. These can be written in the following compact form,

[TABLE]

where, if we expand the cosine terms and introduce the $\bar{\cdot}$ variables, we obtain

[TABLE]

whose average reads, re-using again the $\bar{\cdot}$ variables for their averaged analogue,

[TABLE]

Note that the sums in both Eq. 68 and Eq. 73 are convergent.

3.3.3 Interpretation as an equivalent flat sail

After the averaging procedure, the attitude and orbit dynamics become decoupled since the acceleration due to SRP becomes constant in the direction of Eq. 73, that in fact depends only on the attitude’s initial condition of the integration through the mean amplitude of the oscillations, as will be seen in § 4.

The first benefit of the analysis is that the equations of motion become Hamiltonian again as the averaging acts as fixing the attitude (recall the end of § 2.3). In fact, one can interpret the averaging as follows:

For each initial condition (that fixes the initial value of $\bar{\Phi}$) the dynamics of the averaged equations is that of a spacecraft with the same mass $m_{b}+m_{s}$ with a flat panel of effective area $A_{s}\cdot A_{\rm eff}$, always perpendicular to the Sun-spacecraft direction.

The expression of this effective area is obtained from Eq. 73. Comparing the acceleration due to SRP in our problem and the acceleration due to SRP that a sail with the same mass with area $A_{s}A_{\rm eff}$ would have yields

[TABLE]

and from this we get that

[TABLE]

This quantity is referred to as the area factor.

4 Numerical test cases

This section is devoted to exemplify numerically the analysis performed in § 3. Since the results refer strictly to the case in which both panels face sunlight, that is, $|\hat{\phi}|<\alpha$ in Eq. 17, the simulations are restricted to this case. The equations of motion used in this section are those obtained after the change of variables of Prop. 1, that in the fast time scale $t$ (denoting $(\;^{\prime})={\rm d}/{\rm d}t$) read

[TABLE]

where $(\hat{a}_{\hat{x}},\hat{a}_{\hat{y}})^{\top}$ read are given in Eq. 65.

Note that the theoretical restriction gives a stopping condition for the numerical integrations: if $|\hat{\phi}|>\alpha$ simulations are stopped, as in this case the sail is expected to tumble.

4.1 Physical parameters and adimensionalisation

factors

The system of ODE in Eq. 81 has its own interest for arbitrary choices of the parameters $c_{1},c_{2},c_{3}$ and $c_{4}$. Despite this, to justify the usefulness of the analysis performed in prospective real applications we are lead to choose values of the parameters that correspond to a structure that is constructible according to current technological boundaries. So, as done in [20] and following the guidelines in [8], the physical parameters of the structure are chosen to be $\eta=0.8$, $m_{b}=100$ kg, $w=h=9.20$ m and $m_{s}=3.60$ kg, that corresponds to an area-to-mass ratio of $A_{s}/(m_{b}+m_{s})=0.75$ ${\rm m}^{2}$/kg.

The results are exemplified for spacecraft with $\alpha=35^{\circ},40^{\circ},45^{\circ},60^{\circ}$, all of them with $d=0$, as for these parameters the Sun-pointing attitude is helio-stable, see Fig. 6. It is worth noting that smaller aperture angles $\alpha$ allow smaller oscillation amplitude and hence smaller angular velocity, and any value of $d$ with $|d|>0$ would give rise to a larger size of the gravity gradient perturbation.

Now, the values of the parameters still depend on the adimensionalisation quantities, $L$ for length and $T$ for time, see App. A. There it is justified that to obtain Eq. 81 one has to choose $T=\sqrt{L^{3}/\mu}$ so only $L$ has to be chosen. It has to be done in a way that $\varepsilon=1/c_{1}^{-1/2}$ is small, the gravity gradient torque has to have smaller size than the SRP torque, and the $J_{2}$ and SRP accelerations have to be also smaller than the term of the Kepler problem that is $\mathcal{O}(1)$. In Fig. 7 the values of $\varepsilon=1/c_{1}^{1/2}$, $c_{2}\varepsilon^{2}$ - as it appears as prefactor in the equation of $\hat{\Phi}^{\prime}$ in Eq. 81 -, $c_{3}$ and $c_{4}$ are shown as a function of the adimensionalisation length factor $L$, recall Eq. 18, for $\alpha=45^{\circ},d=0$ m (left) and $\alpha=60^{\circ},d=0$ m (right), as examples. From the two panels in Fig. 7 one infers that $L=20\,000$ km (that is highlighted as a vertical dashed line) is a proper choice for the purposes of this contribution according to the requirements listed above. The figures corresponding to $\alpha=35^{\circ}$ and $\alpha=40^{\circ}$ are qualitatively the same and quantitatively very similar.

With all these choices, the physical constants of the system that are used in the following sections are summarized in Tabs. 2 and 3. It is worth noting that, in both cases, $\varepsilon=\mathcal{O}(10^{-2})$ and that one unit of time $t$ of Eq. 81 is equivalent to $T\varepsilon\approx 2-3$ min.

4.2 Main numerical experiment

The attitude and orbit coupling and averaging results of § 3 are exemplified using one single orbit initial condition on which all the considered disturbing effects play a strong role. Namely, all motions considered start at the perigee of a Keplerian orbit characterized by

[TABLE]

Note that, in particular, the perigee radius is above the surface of planet Earth (with minimal altitude $350$ km) and this remains true along all propagation performed. Despite the low altitude, atmospheric drag is not considered as it was not taken into account in the study of § 3. Such considerations are left for future contributions.

Concerning the attitude initial conditions, a sampling between the Sun-pointing attitude and close to the limit of the validity of the theoretical results, $|\hat{\phi}|<\alpha$, has been performed. More concretely, for each aperture angle $\alpha$, 480 initial conditions of the form

[TABLE]

have been considered. That is, $0<\hat{\phi}\leq 0.90\alpha$ is sampled. This is to make sure that most trajectories do not reach $|\hat{\phi}|>\alpha$ before the maximal integration time. But this may also happen before for these initial conditions, as will be observed later in the case $\alpha=35^{\circ}$.

Finally, the maximal integration time for all considered initial conditions has been 1 year (that corresponds to a different maximal $t$ for each value of $\alpha$).

The theoretical study of § 3 only justifies $\varepsilon$-closeness of the original and averaged system in intervals of length $\mathcal{O}(1/\varepsilon)$ in the adimensional scale $t$. In actual time units, using the data in Tabs. 2 and 3 this length is $\mathcal{O}(10^{2})$ times the time unit, that accounts for $10^{4}\;{\rm s}\approx 3\;{\rm h}$ in the studied 4 cases. Despite this, in the performed numerical experiments one observes that for the considered orbit initial condition this interval is larger. The study of the time intervals where such $\varepsilon$-closeness holds in a large family of orbits is out of the scope of the present paper and hence it is not addressed here.

As a final consideration, as the flow is always transversal to $\hat{x}=0$, instead of the full 6D ODE we have considered the Poincaré section defined by

[TABLE]

Hence, the following is a study of a 5D discrete map on this surface. The reason for this choice is twofold: on the one hand, this is a reduction of the dimension of the phase space by one and this eases the analysis; and on the other hand it allows for better comparison when studying the system in Eq. 81 and its averaged analogue: non-averaged and averaged equations have been integrated with different numerical schemes (an implicit Runge-Kutta Gauss of 2 stages and order 4 and a Runge-Kutta-Fehlberg 7(8), respectively) with automatic stepsize control. The Poincaré section sets a fixed position where to compare the integrated orbits of both equations. Note that the integrated equations are different, and hence the times where the orbits intersect $\Sigma$ do not necessarily coincide, even if the integration starts at the same initial condition.

4.3 Numerical results

A selection of numerical results are shown and described here. These are related to the shape, size and orientation of the osculating orbit along the integration, to be able to compare between different apertures and oscillation amplitudes; to the assessment of the applicability of the equivalent flat sail in the averaged equations; and finally to the assessment of the difference between the original and averaged equations.

4.3.1 Shape, size and orientation of the osculating

orbit

The shape, size and orientation of the orbit of the spacecraft are studied via the Keplerian elements semi-major axis $a$, eccentricity $e$ and $\gamma:=\omega+\Omega$, the sum of the argument of the perigee and the RAAN. Recall that the latter has to be considered as the $J_{2}$ effect makes the ascending node precess. On a scale where the whole evolution along one complete year of integration is displayed, the differences between these three observables is qualitatively the same for all the considered values of $\alpha$.

In Fig. 8 an example of such evolution is displayed. The shown evolution is obtained with the attitude initial condition with $j_{\alpha}=0$ in Eq. 82, that is, the initial attitude closest to the Sun-pointing direction considered. Top, middle and bottom panels show the evolution of $a$, $e$ and $\gamma$ respectively. The differences due to choosing different aperture angles are highlighted in the right column zooms.

This shows that despite the sail has the same shape and size in all cases, the aperture angle of the sail, even if it oscillates mildly around its helio-stable attitude, produces a non-negligible impact on the orbit. Namely in this case we observe $\mathcal{O}(0.1)$ variations of $a$ measured in km, that is, $\mathcal{O}(10^{-2})$ in the adimensional variables; $\mathcal{O}(0.001)$ variations in $e$; and $\mathcal{O}(1)$ variations in $\gamma$, measured in degrees. Other attitude initial conditions show differences of the same orders of magnitude.

4.3.2 The definition of equivalent effective area

From the results concerning the variation of osculating orbit one reads that the fact that the oscillations close to the Sun-pointing direction are fast with respect to the orbit dynamics implies that the sail produces an average uniform effect. This can be interpreted as the QRP sail behaving as a flat sail with area $A_{s}\cdot A_{\rm eff}$, recall Eq. 74.

In Fig. 9 the area factor $A_{\rm eff}$ is plotted, as a function of $\bar{\Phi}$ for different values of $\alpha$ that include the study cases.

Recall that $\bar{\Phi}=(2\hat{\phi}^{2}+\hat{\Phi}^{2})/(2\sqrt{2})$, so this action gives information about the amplitude of the oscillations around the Sun-pointing direction. Figure 9 shows that one can choose the aperture angle in such a way that for small oscillations the effect of the QRP on the SRP acceleration magnitude is, on average, that of more than 2 flat panels of area $A_{s}$. In other words, for adequate $\alpha$, the oscillations are equivalent as to having a larger flat panel always oriented towards the Sun.

The question now is how to measure if this effect can be recovered in the attitude dynamics and SRP perturbation in the performed simulations. Recall that the results of the simulations are data on $\Sigma$ of the propagation of Eq. 81. To test the formulas provided in § 3, one can proceed as follows:

Measure of the average of the action. The quantity $\bar{\Phi}=(2\hat{\phi}^{2}+\hat{\Phi}^{2})/(2\sqrt{2})$ evaluated along the data appears to oscillate in what seems a quasi-periodic manner. Since we are interested in the average, for each initial condition, we denote $\left<\bar{\Phi}\right>$ the time average of $\bar{\Phi}$ along the whole 1-year long integration. 2. 2.

Evaluation of the area factor. With each obtained value of $\left<\bar{\Phi}\right>$ we can evaluate $A_{\rm eff}(\left<\bar{\Phi}\right>,\alpha)$, the area factor, see Eq. 74. This is referred to as the “theoretical” value. 3. 3.

Measure of the average SRP acceleration. The values of the components of the SRP acceleration in Eq. 65 are also averaged along the 1-year long data for each initial condition. As given in this equation, the obtained average values correspond to the area factor divided by either $\cos(\hat{\lambda})$ (along $\hat{x}$) or $\sin(\hat{\lambda})$ (along $\hat{y}$). The mean is denoted as $A_{\rm eff}^{\star}(\left<\bar{\Phi}\right>,\alpha)$, and referred to as the “numerical” value.

The results of this study are summarized in Fig. 10. To allow comparison for different aperture angles, since each aperture angle has a different maximal oscillation, $j_{\alpha}$ is chosen as abscissa instead of $\hat{\phi}_{0}$ itself, recall Eq. 82. The top panel displays the time average of $\bar{\Phi}=(2\hat{\phi}^{2}+\hat{\Phi}^{2})/(2\sqrt{2})$. This shows that the tendency is that the larger $\alpha$ is the faster $\left<\bar{\Phi}\right>$ increases in the $|\hat{\phi}|<\alpha$ region. In this figure three phenomena are highlighted as they require further clarification. Both (d) and (e) correspond to the same phenomenon that occur for $\alpha=40^{\circ}$ and $45^{\circ}$, respectively, that is a plateau. This can be explained by the fact that the corresponding initial conditions do not evolve on a torus or close to a torus but on a large width chaotic region where the attitude can range freely. Even though the phase space of our problem is 5D and confining manifolds in this context have dimension 4, a geometrical analogy would be the dynamics of a chaotic orbit in a Birkhoff region of an area-preserving map (2D): a region that is bounded by invariant curves (confining manifolds of dimension 1) with no other invariant curves inside, where chaotic orbits are confined and their iterates eventually fill a positive measure set that range the whole width between the two confining curves. The variations in (f) that occur for $\alpha=35^{\circ}$ are spurious data, as the corresponding sail reached $|\hat{\phi}|>\alpha$ before 1 year.

The bottom left panel of Fig. 10 shows the theoretical and numerical values of the area factor. The bottom right panel shows the difference between the theoretical and numerical values of $A_{\rm eff}$ (as defined in the enumeration in the beginning of this section)

[TABLE]

On the left the theoretical and numerical values are displayed in dashed and solid lines, respectively. As expected from the hypotheses of the theory, the fit is better the closer to the Sun-pointing direction we are, but even for larger values of $\left<\bar{\Phi}\right>$ Eq. 74 still gives a good first approximation of the area factor.

4.3.3 Full system versus averaged system

The application of the averaging method converts Eq. 81 into the simplified version that only contains orbit dynamics whose equations read

[TABLE]

recall Eq. 66. Note that here the motion can be integrated in the slow time scale $\tau$ so the notation $(\dot{\;})={\rm d}/{\rm d}\tau$ has been used. To test numerically how close the orbits of Eq. 88 are to those of Eq. 81, we use the numerically evaluated area factor, that is, we use Eq. 88 with $A^{\star}_{\rm eff}$ instead of $A_{\rm eff}$.

The initial conditions for this orbit propagation are the same ones as explained in § 4.2. That is, for each value of $\alpha$ considered, we integrate 480 orbits (one per value of $j_{\alpha}$), for 1 year, keeping track of the intersections of each orbit in $\Sigma$, see Eq. 83.

As above, we study the shape, size and orientation of the osculating ellipse via $a$, $e$ and $\gamma$ and the results are compared with the orbits explained in § 4.3.1. On a scale where the whole evolution along one complete year of integration is seen, the difference between the evolution of $a,e$ and $\gamma$ in the averaged and original equations are not noticeable. Their behaviour is such as that shown in Fig. 8. To study the discrepancies between the numerical results of these two equations, it is convenient to compare the adimensional semi-major axis, eccentricity and orientation, that shall be denoted as $\tilde{\tilde{a}}$, $\tilde{\tilde{e}}$ and $\tilde{\tilde{\gamma}}$ for the averaged equation Eq. 88 and $\hat{a}$, $\hat{e}$ and $\hat{\gamma}$ for the original one Eq. 81. A sample of the numerical results obtained can be seen in Fig. 11. It is important to stress that the comparison is done at the intersection of the orbits with $\Sigma$. Hence the abscissas are the number of iterates in $\Sigma$ (a total of $1789440$ in all cases), not the integration time as the intersection times in $\Sigma$ are not necessarily the same in the two compared problems.

The shown results correspond to 19 different values of $\left<\bar{\Phi}\right>$: $j_{\alpha}=0:25:480$. Concerning the differences in semi-major axes, along this first integration year, they are $\mathcal{O}(10^{-2})$; concerning the differences in eccentricity, they are $\mathcal{O}(10^{-5})$; and those for the orientation, they are at most $\mathcal{O}(10^{-4})$. Taking into account that $\varepsilon=\mathcal{O}(10^{-2})$, these numerical results fit within the accuracy expectations but the time span where these seem to be valid, in the chosen example, exceed the predictions done by classical theorems that are summarized in Th. 1.

As a final remark, note that these numerical results depend strongly on the chosen orbit, and a specific numerical study to investigate this fact should be done to assess this dependence. For the studied orbit (as it is indicated in the numerical results for $A_{\rm eff}$) one can find a non-negligible region of practical stability around the Sun-pointing attitude so that it makes sense to consider such structures as potential candidates for passive deorbiting devices. The studied orbit initial condition example is highly eccentric and produces a large gravity gradient torque perturbation. Attitude stability is expected to be enhanced in high altitude orbits, even all the way down to low Medium Earth Orbits (MEO), i.e. 2000 km of altitude, and for less eccentric orbits.

5 Conclusions

The coupled attitude and orbit motion is known to be a problem with two characteristic time scales. Although this property is not necessarily obvious in the form of the equations of motion, in particular for the class of spacecraft studied in this paper, -a simplified Quasi-Rhombic-Pyramid (QRP) that consists of a symmetric payload attached to an already deployed sail that is composed of two reflective panels forming an angle that endows it with helio-stability properties- the ratio between time scales, $\varepsilon$, can be found explicitly. An adequate time and phase scaling that depend on $\varepsilon$ can be performed to highlight the two different characteristic times of the motion. This parameter, in turn, depends solely on physical quantities that describe the shape and mass distribution of the spacecraft under consideration.

The attitude dynamics can be understood via a formal treatment of the equations by studying the fast subsystem of the motion by artificially setting the parameter $\varepsilon$ to [math] with equations written in the fast time scale. In the specific case of the simplified QRP it is a Hamiltonian system of one degree of freedom with $\mathcal{C}^{0}$ equations of motion that resembles a pendulum. The Sun-pointing attitude is a stable equilibrium provided an adequate center of mass - center of pressure offset and aperture angle are chosen. The phase space inside the separatrices is foliated by periodic orbits with different periods but both the Hamiltonian and their equations of motion are only analytic in a neighborhood of the Sun-pointing attitude. A change of variables of Poincaré type has to be performed inside this neighbourhood to force all attitude periodic orbits to have the same period.

The separation of the attitude and orbit dynamics by means of direct averaging of the orientation of the sail with respect to the sunlight direction provides a reasonable first approximation of the orbit dynamics. The theoretical justification of the applicability of averaging results rely mainly on two independent requirements: smallness of $\varepsilon$ and initial closeness to the Sun-pointing attitude. Also, the smoothness hypotheses of classical theorems only allow to restrict to the case where both panels face sunlight.

Even though in the studied practical examples, where the chosen physical data is that of a constructible structure, the value of $\varepsilon$ is not necessarily small, the oscillations around the Sun-pointing direction are fast enough so that the effect along the whole integration is perceived in the orbit dynamics as uniform.

This uniform effect can be interpreted as the sail structure consisting of a flat panel with the same reflectance as the original panels; and whose area is the area of one of the panels times a factor that depends on the amplitude of the oscillations and the aperture angle. By tuning these two parameters, this factor can be made larger than 2, hence giving rise to an effective area-to-mass ratio that is at least double the amount if we only considered one of the panels of the structure, but with the advantage of the enhanced stability provided by the oscillating character.

The formulas for the area factor could be improved by performing the changes of variable in a way that they have zero average. This would introduce, in particular, terms due to the gravity gradient in the area factor. The great disadvantage is that the kind of expressions one would be forced to handle would be increasingly involved as one goes further in the sequence of averaging steps. So, one should take into account the trade off between these tedious computations and the possible improvements in the formula, as the way they are presented in this contribution may be enough as a first approximation for practical purposes.

This work is intended to be a first step towards the comprehension of the long term dynamics of QRP. Future lines of research include the extension of the results of this paper to the 3D structure. In this situation it makes sense to consider eclipses as a moderate spin along one of the axes of inertia can be considered to enhance attitude stability, as done in [11]. Such assumption would reduce the rotation dynamics from 3 to 2 degrees of freedom. Other interesting lines of research would include the consideration of damping effects to enhance the stability properties, and the possible transition scenarios from SRP dominated to atmospheric drag dominated regions where such as structure should be used as a drag-sail.

Concerning applications, the long term stability of these families of spacecraft make them suitable for cargo transportation missions and interplanetary transfers, see e.g. [18].

Acknowledgments

The research leading to these results has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme as part of project COMPASS (Grant agreement No 679086). The authors acknowledge the use of the Milkyway High Performance Computing Facility, and associated support services at the Politecnico di Milano, in the completion of this work. The datasets generated for this study can be found in the repository at the link www.compass.polimi.it/$\sim$ publications. Fruitful conversations with I. Gkolias, M. Lara and C. Simó are also acknowledged.

Appendix A Adimensional model

In this section the adimensional set of differential equations are derived. This is done for the sake of providing a set of equations written in the form of a fast-slow system where averaging theorems can be formally applied.

To fix notation, here the superscript $\tilde{\cdot}$ is used to refer to the original variables with dimensions, not to be confused with those introduced in the averaging Lemma 1. The variables without $\tilde{\cdot}$ refer to the adimensional analogues. Let $L$ and $T$ denote the length and time units, measured in km and sec, respectively. Finally,

[TABLE]

are used to denote derivatives with respect to dimensional and adimensional time, respectively.

Consider the relations

[TABLE]

The length $L$ and time $T$ units are chosen to normalize the Kepler problem in such a way that it simplifies to the vectorial equation

[TABLE]

The first component of the dimensional Kepler problem reads $\ddot{\tilde{x}}=-\mu\tilde{x}/\tilde{r}^{3}$, where $\tilde{r}=\sqrt{\tilde{x}^{2}+\tilde{y}^{2}}$. From Eq. 89, $\tilde{r}=Lr$, and also

[TABLE]

if we choose

[TABLE]

It is clear that for this choice the second component of the adimensional Kepler problem reads $y^{\prime\prime}=-y/r^{3}$. The choice of $L$ is critical so that the magnitudes of the constants that appear in the problem are adequate for the applicability of the known averaging results. This is done in § 4.

Before scaling the equations of motion, it is convenient to shift the angle $\tilde{\varphi}$: let $\tilde{\phi}=\tilde{\varphi}-\tilde{\lambda}$, where $\tilde{\lambda}$ is the argument of latitude the Sun (i.e. the angle between the $x$ axis and its position vector in Earth-Centered coordinates, measured counter-clockwise) in its apparent motion around the Earth.

In the adimensional variables, the attitude equations of motion read

[TABLE]

as the right hand side of Eq. 4 has to be multiplied by $T^{2}$. On the other hand, the translational equations of motion read

[TABLE]

where $a_{x}$ and $a_{y}$ are given in Eq. 10 putting $\varphi=\phi+\lambda$.

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.