Attitude and orbit coupling of planar helio-stable solar sails

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-stability properties. After adimensionalization, the system is put in the form of a fast–slow dynamical system where the different timescales are explicitly 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.


List of symbols a
Semimajor axis (km) α Aperture angle of the sail (deg or rad) A s Area of the panels (m 2 ) β Frequency of the Sun-pointing direction (deg or rad) C Inertia moment around the ζ axis (kg m 2 ) d Center-of-mass-center-of-pressure offset (m) e Eccentricity [-] ε Ratio between characteristic timescales [-] F b Body frame Difference ω − (deg or rad)

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 (Jorba-Cuscó et al. (2016)). 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, which 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 Colombo and Bras (2016), that consists of using the idea of deorbiting "outwards" on an elliptical orbit: The increase in the eccentricity of the orbit causes the perigee radius to progressively decrease. As justified in Lücking et al. (2012Lücking et al. ( , 2013, 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 Borja and Tun (2006), where the authors suggested maximizing the SRP acceleration when traveling toward the Sun and minimizing it when traveling 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 Colombo et al. (2019). This refinement reduces the number of attitude change maneuvers from twice per turn around the Earth to twice per year; see Colombo and Bras (2016).
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 Ceriotti et al. (2013). 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 Felicetti et al. (2016) 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 Felicetti et al. (2017).
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 Ceriotti et al. (2014), and later exploited in Heiligers and Ceriotti (2017) 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.
Although the contributions (Ceriotti et al. 2013(Ceriotti et al. , 2014Felicetti et al. 2016Felicetti et al. , 2017Heiligers and Ceriotti 2017) provide satisfactory results, there is, to the authors' knowledge, a lack of understanding of 1. The attitude dynamics, especially close to the Sun-pointing attitude, and 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 timescales. These two questions are addressed in this paper by considering a sail structure in between the single panel considered in Ceriotti et al. (2014) and Heiligers and Ceriotti (2017) and the full QRP (Ceriotti et al. 2013;Felicetti et al. 2016Felicetti et al. , 2017) 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 Miguel and Colombo (2018). 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., Lücking et al. 2012;Colombo et al. 2012), here the considered orbital dynamics are the J 2 -perturbed two-body problem 1 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 Miguel and Colombo (2018), 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 timescales 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, Sect. 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 timescales, 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 timescales are split, and related via a physical parameter that represents the ratio between timescales that only depends on the geometry of the spacecraft. In Sect. 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 Sect. 3 are tested in Sect. 4 with special emphasis on the physical interpretations just mentioned. The paper concludes in Sect. 5, where the main results of the contribution are summarized and different possible lines for future research are suggested.

Model
This section is devoted to providing the equations of motion of the planar dynamics of a heliostable 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 Sects. 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 Miguel and Colombo (2018). The section ends by putting the equations of motion in the context of dynamical systems with multiple timescales in Sect. 2.3.

Geometry of the sail structure
The spacecraft under consideration consists 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 (Ceriotti et al. 2013) that consists of two panels of equal size P ± and 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 F b attached to it. The coordinates in this frame are referred to as ξ, ν and ζ and the vectors of the basis i ξ , i ν , i ζ . The panels are attached to each other along an h-long side that lies on a line parallel to the ζ axis, and they form an angle α with respect to the plane ν = 0. The payload, whose mass is denoted as m b , is assumed to be on the ξ axis, at a distance d from the center of mass of the 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. Left: d < 0. Center: d = 0. Right: d > 0 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 F b . The main physical parameters of the system are: α, the aperture angle, and d, which accounts for the MPO. Sketches of top views of the spacecraft in F b are shown 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 α is measured in Fig. 2) did not produce any SRP acceleration. In Miguel and Colombo (2019), the attitude dynamic 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 Miguel and Colombo (2019).

Equations of motion
The considered planar orbit and attitude dynamics are a coupled system of differential equations in (S 1 × R) × R 4 , where S 1 := R/(2πZ): orientation and angular velocity for the attitude; and position and velocity of the spacecraft in an Earth-centered inertial reference frame 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 F I . The coordinates of F I are denoted x, y and z, and the vectors of the orthonormal basis are denoted i x,y,z . The vector i x points toward an arbitrarily chosen direction on the ecliptic (e.g., J2000), and since we are dealing with the planar problem, the vector i z is parallel to i ζ , and they are also perpendicular to the ecliptic plane. The triad is completed by choosing i y = i z × i x .
As the motion is planar the change of coordinates from F I to F b is done through R 3 (−ϕ), where ϕ ∈ S 1 is an Euler angle and R 3 is the rotation matrix around the z (and ζ ) axis. The rotation matrix R 3 reads, for any angle ψ ∈ [0, 2π), The Euler attitude equations in the present situation reduce to where M = (M ξ , M ν , M ζ ) is the torque due to the external forces considered and C is the inertia moment around the ζ axis in F b . Denote I ξ,b , I ν,b , I ζ,b as the inertia moments of the bus. Then one can see that (2)

SRP model
The force due to SRP exerted on each panel of the sail in F b is considered to be modeled as Markley and Crassidis (2014) where u is the unit vector in the Earth-Sun direction and n ± are the normal vectors to each panel; recall Fig. 1. Concerning the constants, η ∈ (0, 1) is the (dimensionless) reflectance of the sail and p SR = 4.56 × 10 −6 N/m 2 is the solar pressure at 1 AU which is considered to be constant.

Attitude dynamics
The effects taken into consideration are SRP and the non-symmetry of the spacecraft, so the total torque is M = M SRP + M GG , the sum of the SRP and gravity gradient torques. Let λ 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 ϕ relative to λ, so denote φ = ϕ − λ. These three angles are sketched in Fig. 3. If φ ∈ [−π + α, α], the panel P + produces torque (see Fig. 4a); and if φ ∈ [−α, π − α], it is the panel P − which produces torque (see Fig. 4b); in particular, if |φ| ≤ α, both panels face the Sun (see Fig. 4c). In all panels of Fig. 4, the sunlight direction is depicted as if it was in the direction (−1, 0) ; hence, λ = 180 • ; see Assume that the bus is symmetric in the sense that I ξ,b = I ν,b = I ζ,b . In this case, as derived in Miguel and Colombo (2018), the attitude equations of motion reduce to the following second-order ordinary differential equation (ODE) where D(α, d) is defined in Eq. 2; θ , ω and 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 Sect. 2.2.3, and where X J is the characteristic function of the interval J , The parameter μ = G M = 3.986 × 10 14 m 3 /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 η, and read 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, which takes into consideration that the panels have to face the sunlight to produce torque; this is why M 1 is a piecewise defined function.

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 Colombo et al. 2012;Lücking et al. 2013;Colombo and Bras 2016). The (dimensional) equations of motion of the orbit dynamics read 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 × 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 = 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. 7 are piecewise defined, as the SRP acceleration depends on ϕ − λ. Namely, a x and a y can be expressed as follows: where are the expressions of the adimensional factor of the accelerations due to panel P + (subscript +) and due to panel P − (subscript −). Equations 9 are obtained from Eq. 3, with n ± = R 3 (ϕ)(sin α, ± cos α) , u = (cos λ, sin λ) .
In the particular case that α = π/2 and ϕ = λ, 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 twice the acceleration of a flat panel of area A s = hw always oriented toward 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 (Jorba-Cuscó et al. (2016)).

Multiple timescale dynamics
The full set of 6 coupled differential equations, Eqs. 4 and 7, have two different timescales, 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 1" is devoted to the adimensionalization of the equations of motion, where the adimensionalization 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 (˙) denotes the derivative with respect to the adimensional time, τ . The full set of equations read where φ = ϕ − λ (depending on adimensional time), a x and a y are as in Eq. 8, and the constants are the adimensional quantities Written like this, the problem can be put in the form of a fast-slow system. To write Eq. 10 as an ODE of first order, one must introduce =φ =φ −λ, and v x =ẋ, v y =ẏ. Let us denote φ = (φ, ) and x = (x, y, v x , v y ). 2 Then the following holds.
Proposition 1 Assume thatλ = n = constant and that the torque coefficient for adequate f and g.
Proof If we define the proposition follows by considering the scalingφ The maps f and g correspond to the right-hand side of the resulting equations of motion forφ andx, respectively.
The parameter ε in Eq. 13 depends solely on physical quantities of the system and on the choice of the longitude scaling factor L in the adimensionalization procedure; see App. 1. If L is appropriately chosen, see Sect. 4.1, then ε is small, and Eq. 12 is a fast-slow system, written in the slow timescale τ . One can consider a time scaling t = τ/ε, the fast timescale, in which Eq. 12 reads ⎧ ⎨ Hence, ε is the ratio between the timescales in which the two characteristic motions, orbit and attitude, take place. The change to· variables can be extended to the argument of latitude of the Sunλ = λ which, in the fast timescale t, varies as dλ/dt = εn , similarly as for the orbital dynamics in Eq. 12. That is, from the point of view of the attitude, the position of the Sun varies in a slower constant rate, with factor ε.

Fast equations
The components of the vector field in Eq. 15 of the fast variables are those representing the evolution of the sail attitude where, if we denote we can write 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.

Slow equations
Concerning the components of the vector field in Eq. 15 of the slow variables, the scaling of the proof of Prop. 1 is such thatx = x, and hence, the form of the equations does not change. The vector field of the slow subsystem reads wherer = x 2 +ŷ 2 +ẑ 2 andr = x 2 +ŷ 2 and . Written as in Eqs. 12 and 14, 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 ε = 0, which can be interpreted as the spacecraft having an infinitely large area-to-mass ratio. Setting ε = 0 has two different meanings in each equivalent formulation, Eqs. 12 and 14: • In the slow timescale, Eq. 12, 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. 18 can be obtained as the derivatives of The time dependency comes from the SRP acceleration, as it depends on the position of the Sun. • In the fast timescale, Eq. 14, 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 ε = 0 both the slow and fast timescale dynamics are Hamiltonian: In the slow timescale, this holds for any fixed attitude, and in the fast timescale, it is not trivial and is justified in Sect. 3.1. However, the whole problem is not Hamiltonian due to the SRP and gravity gradient coupling.

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. 14, written in the fast timescale, the parameter ε is set to 0, 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 timescales, for small values of ε the variables that are frozen for ε = 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 ε > 0, via the so-called geometrical singular perturbation theory, based on Fenichel's theory on the preservation of normally hyperbolic invariant manifolds under perturbation. The expository text (Kuehn 2015) 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 Sunpointing attitude, which is expected to be preserved for long periods of time. Moreover, the fast dynamics can be averaged out.
First, Sect. 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 Sect. 3.2. This allows to obtain a physical interpretation of the results that are exposed in Sect. 3.3.

Dynamics of the fast subsystem
Consider Eq. 14, the full system written in the fast timescale t. By setting ε = 0, the orbit (slow) vector field vanishes, so the only non-trivial equations are those corresponding to the attitude (fast) dynamics, which read recall Eq. 15. This is not equivalent to the simplified model found and studied in Miguel and Colombo (2018), 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. 21 for ε = 0 is Hamiltonian with some Hamiltonian function K 0 and hence can be written as Proof If Eq. 21 was a Hamiltonian system with Hamiltonian K 0 , the equations of motion would be obtained as derivatives of K 0 as indicated in Eq. 22. As the first component of f isˆ , K 0 must haveˆ 2 /2 as summand, and as the second component of f is M 1 , another summand of K 0 would be an appropriate primitive of M 1 ; see Eq. 17. The function M 1 is 0 , so one has to make sure that this primitive is C 1 . If we denote then a solution is is the Hamiltonian function of a pendulum. 2. For ε > 0, the attitude equations of motion can also be written as the derivatives of a function K ε with respect toφ andˆ . Namely, if we extend Eq. 23 then for ε > 0 we recover the vector field Eq. 15 via The most relevant properties of the fast vector field f in Eq. 22 are: It is C 0 , the differentiability being lost at the switching manifoldsφ = ±π ∓ α, ±α; it is 2π-periodic inφ; and it is also symmetric with respect toφ = 0 as K ε ((φ,ˆ ),x) = K ε ((−φ,ˆ ),x). An example of the phase space can be seen in the left panel of Fig. 5, for α = 30 • and d = 0. In this figure, each curve represents an orbit of the fast subsystem Eq. 21 that is obtained as 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 (φ, 0) witĥ φ ∈ I − ∪ {0} ∪ I + , where I − = [−π, −π + α] and I + = [π − α, π) 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 ± = (±π ∓ α, 0), which 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 ± is sketched: E is the Sun-pointing attitude and H ± represent the angles of transition to from no reflective panel facing sunlight to one or vice versa. The rest of the equilibria, those whose abscissa isφ ∈ I − ∪ I + \ {−π + α, π − α}, have 0 as double eigenvalue.

Fast-slow dynamics and averaging
This section is devoted to the study of the dynamics of the full system for ε > 0. For small values of ε, one expects the structure found in Sect. 3.1 to be close to conserved in a sense that is made explicit here. More concretely, the function K ε is a first integral of the frozen-fast subsystem for ε = 0, so for ε > 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. 12 is not Hamiltonian): There is an equivalent formulation of the system in which one of the variables experiences oscillations of at most O(ε) for time intervals of length O(1/ε). 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φ = 0 are those for which there is numerical evidence of helio-stability properties of the spacecraft, see Miguel and Colombo (2018), 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 timescale in Eq. 14, let G be a bounded subset of the phase space such that, for all (φ,x) ∈ G, |φ| < α (that is, both panels face the sunlight) where all trajectories lie wholly in G. Note that, in G, for fixed values of x, the condition K ε = k defines one and only one trajectory of Eq. 14, and the restriction of G onto the (φ,ˆ ) variables consists of a continuum of closed and nested periodic orbits.
Lemma 1 For small enough ε > 0 and small enoughˆ (that is, close enough to the Sunpointing attitude), there exists a real analytic change of variables defined in the interior of G C : (φ,ˆ ,x,ŷ,vx ,vŷ) → (φ,˜ ,x,ỹ,ṽx ,ṽỹ) that transforms Eq. 14, with f and g as given in Eqs. 15 and 18, respectively, into a system of the form wheref ,g have zero average with respect toφ.
Proof This result fits within the scope of averaging theory, and the proof can be sketched as consisting of the two following steps.
1. Consider for the moment the dynamics of the fast subsystem Eq. 22. It depends onφ in a periodic way, but for each initial condition in G, the period and the range ofφ of points in this orbit are different. Restricting to the behavior close to the Sun-pointing direction, we can substitute the equations of motion by the Taylor series aroundφ = 0. Since the frozen subsystem is Hamiltonian, one can have K 0 in mind for the moment. The leading non-constant terms are This suggests to consider a first change of variables (· →·) ˆ ,x,ŷ,vx ,vŷ) → (φ,¯ ,x,ȳ,vx ,vȳ), where (¯ ,φ) are the usual Poincaré action-angle variableŝ the frequency being β = √ 2. The change is extended to the rest of the variables by simply choosingx =x,ȳ =ŷ,vx =vx andvȳ =vŷ. After this change of variables, the obtained equations of motion are 2π periodic inφ, 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 K ε . So, for orbits in G, we can obtain the equations of motion in the· variables as follows. Start by considering where ϑ := arctan(ŷ/x) −λ. Now, the expansion aroundφ = 0 reads where, if we introduce the· variables Eq. 26 with β = √ 2, it reads Now the attitude equations in· variables can be recovered via Concerning the orbital motion (or slow subsystem), it only depends onφ via the term due to SRP acceleration; see the rightmost summand in Eq. 18. In the region G, this acceleration (see Eq. 19) reads This can be treated as we did for K ε above, first separating the dependence onφ, substituting the sine and cosine functions by their Taylor expansions, and then introducing thē · 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. Let us focus for the moment on the differential equation of¯ , that is, after performing the first change of variables C 1 . This equation is obtained as −∂K ε /∂φ. Note that all summands except β¯ in the right-hand side of Eqs. 28 have a factor that depends onφ, only as an argument of a power of a sine function, so after deriving with respect toφ, the derivative of all terms appears in the right-hand side of the equation of motion of¯ . Moreover, in each summand there is a factor of the form of the right-hand side of d dφ sin lφ = l cosφ sin l−1φ , l ∈ Z, l > 0 that have zero average with respect toφ. Hence, one can get rid of the dependence onφ via some averaging steps as done, for instance, in the classical reference Neishtadt (1984). After ordering the terms (by orders in¯ , for instance), the jth 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φ) 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 β > 0 has to be bounded away from zero. This is satisfied as β = √ 2. On the other hand, the leading terms of the expansion are O(¯ 2 ) (see Eq. 28a); hence, we have to assume that¯ 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. 25, one has to perform a second change of variables that is a composition of N ≥ 0 averaging steps. The number of steps N depends on the size of |¯ | relative to ε. Note that there are terms in the expansions in Eqs. 28 that do not have ε as factor, namely those in the summatory in Eq. 28a, which start with¯ 2 . It may happen, if we are far enough from the equilibrium that O(¯ 2 ) > O(ε 2 ); that is, these terms dominate over those in Eqs. 28b and 28c. Let l ≥ 2, l ∈ Z be the first power such that O(¯ l ) < O(ε 2 ). To reach the desired form in the equations of motion in Eq. 25, we have to get rid of all terms in Eq. 28a that have¯ 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., Efthymiopoulos 2011). As the first term has order 2 in¯ , the number of averaging steps to perform is N = l − 2. In particular, if O(¯ 2 ) is already smaller than O(ε 2 ), all terms are small enough, and then, C 2 is the identity. For small enough¯ , C 2 is the composition of near-the-identity changes of variable such that, in the new variables, the lowest order term in¯ in the Taylor expansion that depends onφ is exchanged by a zero-average term with respect toφ of the same order in¯ . 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· variables, at the expense of having a change of variables that does not have zero average, plus higher-order terms that have ε as factor.
After performing C 1 and then C 2 , Eq. 12 has, in the· variables, the claimed form Eq. 25, wheref is the derivative with respect to¯ of the average with respect toφ of the terms in the sum Eq. 28a up to the first order comparable with ε in magnitude, with˜ in the place of¯ ;g consists of the slow equations of motion plus the average up to the same order of Eq. 29; and the functions R 1 , R 2 and the map R contain the non-averaged terms plus extra terms that appear as a by-product of the successive averaging changes of variables. Now we are in a position of using the averaging principle, which consists of replacing the right-hand side of Eq. 25 with its average with respect toφ and to neglect the terms O(ε 2 ). Let us refer to the averaged variables as·. The averaged equations of motion read where the angle brackets denote average with respect toφ. Recall thatφ indicates the average angle between the ξ axis in F b (i.e., the orientation of the sail) and the direction of sunlight. Consider the solutions of Eqs. 25 and 30 starting at the same initial conditions, and denote them by˜ (t),x(t) and˜ (t),x(t), respectively.
Theorem 1 Assume ε > 0 is small enough. Then, on the one hand,˜ is an adiabatic invariant, that is,

On the other hand, starting with the same initial conditions, the solution of the slow variables in Eq. 25 is described by the solution of Eq. 30 with an accuracy of O(ε) over time intervals of length O(1/ε), that is,
Proof This follows from the theorem about the accuracy of the averaging method; see Arnold (1988), Chapter 10, Sect. 52. This applies as the equations of motion are analytic, so they meet the minimum (finite) differentiability requirements, and |β| = √ 2 > 0 is bounded away from zero.

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φ. 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. 25 more involved. If instead one performs these steps getting rid of all possible periodic terms, the changes and the remainders remain 2π periodic through the whole process. Moreover, in this situation, the classical theorem in Neishtadt (1984) is applicable, so there exists a change of variables that separates the phaseφ from the rest of the variables except from a remainder that has exponentially small bounds.

Explicit equations and physical interpretation
The benefits of the analysis performed in this section are 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.

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, which is a hypothesis of Prop. 1, is actually equivalent to This is a necessary constraint between the physical parameters of the system, α and d. Note that since α ∈ [0, π/2] and η ∈ (0, 1), the function K (α, η) < 0 we have the following result. 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 shown in Fig. 5.
This result relates the aperture angle to 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 α ∈ (0, π/2) and d min = 0 for α = π/2, the flat plate case; recall Fig. 2. In Fig. 6, the quantity d min is displayed as a function of α, for different values of the reflectance coefficient η.

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 ε, which actually measures the separation of timescales between the attitude and orbit components. This was already observed in the numerical study performed in Miguel and Colombo (2018).
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 ε and¯ are small, the averaging process can be carried out up to any required order. Namely, the larger 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 K ε and then computing the derivatives with respect to the new averaged action. Namely, using that 1 2π 2π 0 sin n ψ dψ = 1 2 n n n/2 if n is even 0 i fn is odd , the average of Eq. 28 reads (using the same names for the variables) The derivative of K ε with respect to¯ gives the approximation of the average frequency of rotation ofφ.
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. 29.
Let us proceed as done from Eqs. 27 to 28, that is, first separating the dependence onφ and then expanding the sine and cosine functions aroundφ = 0 in Eq. 29. After this procedure, the only terms that will contribute to the average are with those that have cos(lφ), l ∈ Z as factor. These can be written in the following compact form, where, if we expand the cosine terms and introduce the· variables, we obtain whose average reads, reusing again the· variables for their averaged analogue, Note that the sums in both Eqs. 32 and 33 are convergent.

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. 33, which 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 Sect. 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 Sect. 2.3). In fact, one can interpret the averaging as follows: For each initial condition (that fixes the initial value of¯ ) 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 · A eff , always perpendicular to the Sun-spacecraft direction.
The expression of this effective area is obtained from Eq. 33. 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 eff would have yields and from this, we get that This quantity is referred to as the area factor.

Numerical test cases
This section is devoted to the numerical exemplification of the analysis performed in Sect. 3. Since the results refer strictly to the case in which both panels face sunlight, that is, |φ| < α in Eq. 10, 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, which in the fast timescale t (denoting ( ) = d/dt) read where (âx ,âŷ) read are given in Eq. 29. Note that the theoretical restriction gives a stopping condition for the numerical integrations: If |φ| > α simulations are stopped, as in this case the sail is expected to tumble.

Physical parameters and adimensionalization factors
The system of ODE in Eq. 35 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 Miguel and Colombo (2018) and following the guidelines in Dalla Vedova et al. (2018), the physical parameters of the structure are chosen to be η = 0.8, m b = 100 kg, w = h = 9.20 m and m s = 3.60 kg, which corresponds to an area-to-mass ratio of A s /(m b + m s ) = 0.75 m 2 /kg.
The results are exemplified for spacecraft with α = 35 • , 40 • , 45 • , 60 • , 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 α 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 adimensionalization quantities, L for length and T for time; see App. 1. There it is justified that to obtain Eq. 35 one has to choose T = L 3 /μ, so only L has to be chosen. It has to be done in a way that ε = 1/c −1/2 1 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 O(1). In Fig. 7, the values of ε = 1/c 1/2 1 , c 2 ε 2 -as it appears as prefactor in the equation ofˆ in Eq. 35-c 3 and c 4 are shown as a function of the adimensionalization length factor L, recall Eq. 11, for α = 45 • , d = 0 m (left) and α = 60 • , 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 α = 35 • and α = 40 • 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 Tables 1 and 2. It is worth noting that, in both cases, ε = O(10 −2 ) and that one unit of time t of Eq. 35 is equivalent to T ε ≈ 2 − 3 min.

Main numerical experiment
The attitude and orbit coupling and averaging results of Sect. 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 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 Sect. 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, |φ| < α, has been performed. More concretely, for each aperture angle α, 480 initial conditions of the form have been considered. That is, 0 <φ ≤ 0.90α is sampled. This is to make sure that most trajectories do not reach |φ| > α before the maximal integration time. But this may also happen before for these initial conditions, as will be observed later in the case α = 35 • . 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 α).
The theoretical study of Sect. 3 only justifies ε-closeness of the original and averaged system in intervals of length O(1/ε) in the adimensional scale t. In actual time units, using the data in Tables 1 and 2 this length is O(10 2 ) times the time unit, which accounts for 10 4 s ≈ 3 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 ε-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 tox = 0, instead of the full 6D ODE we have considered the Poincaré section defined by 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 in 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. 35 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 step size 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 do not necessarily coincide, even if the integration starts at the same initial condition.

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.  Fig. 8 Evolution of the semimajor axis a (top row), eccentricity e (middle row) and γ = ω + (bottom row) , the sum of the argument of the perigee and the RAAN of the osculating orbit for the initial condition j α = 0 in Eq. 36

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' semimajor axis a, eccentricity e and γ := ω + , 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 are qualitatively the same for all the considered values of α. In Fig. 8, an example of such evolution is displayed. The shown evolution is obtained with the attitude initial condition with j α = 0 in Eq. 36, that is, the initial attitude closest to the Sun-pointing direction considered. Top, middle and bottom panels show the evolution of a, e Φ [-] and γ , respectively. The differences due to choosing different aperture angles are highlighted in the right column zooms. This shows that although 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 O(0.1) variations of a measured in km, that is, O(10 −2 ) in the adimensional variables; O(0.001) variations in e; and O(1) variations in γ , measured in degrees. Other attitude initial conditions show differences of the same orders of magnitude.

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 · A eff ; recall Eq. 34.
In Fig. 9, the area factor A eff is plotted, as a function of¯ for different values of α that include the study cases.
Recall that¯ = (2φ 2 +ˆ 2 )/(2 √ 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 α, the oscillations are equivalent as to having a larger flat panel always oriented toward the Sun.
The question now is how to measure whether 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 of the propagation of Eq. 35. To test the formulas provided in Sect. 3, one can proceed as follows: 1. Measure of the average of the action. The quantity¯ = (2φ 2 +ˆ 2 )/(2 √ 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 ¯ the time average of¯ along the whole 1-year-long integration. 2. Evaluation of the area factor. With each obtained value of ¯ we can evaluate A eff ( ¯ , α), the area factor; see Eq. 34. This is referred to as the "theoretical" value. Average ofΦ vs j α Err [-] Theoretical A eff vs. numerical A eff The values of the components of the SRP acceleration in Eq. 29 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(λ) (alongx) or sin(λ) (alongŷ). The mean is denoted as A eff ( ¯ , α) 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 α is chosen as abscissa instead ofφ 0 itself; recall Eq. 36. The top panel displays the time average of = (2φ 2 +ˆ 2 )/(2 √ 2). This shows that the tendency is that the larger the α is, the faster the ¯ increases in the |φ| < α 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 α = 40 • and 45 • , 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 α = 35 • are spurious data, as the corresponding sail reached |φ| > α 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 eff (as defined in the enumeration in the beginning of this section).
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 ¯ Eq. 34 still gives a good first approximation of the area factor.

Full system versus averaged system
The application of the averaging method converts Eq. 35 into the simplified version that only contains orbit dynamics whose equations read recall Eq. 30. Note that here the motion can be integrated in the slow timescale τ , so the notation (˙) = d/dτ has been used. To test numerically how close the orbits of Eq. 38 are to those of Eq. 35, we use the numerically evaluated area factor; that is, we use Eq. 38 with A eff instead of A eff . The initial conditions for this orbit propagation are the same ones as explained in Sect. 4.2. That is, for each value of α considered, we integrate 480 orbits (one per value of j α ), for 1 year, keeping track of the intersections of each orbit in ; see Eq. 37.
As above, we study the shape, size and orientation of the osculating ellipse via a, e and γ and the results are compared with the orbits explained in Sect. 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 γ in the averaged and original equations is not noticeable. Their behavior 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 semimajor axis, eccentricity and orientation, which shall be denoted asã,ẽ andγ for the averaged equation Eq. 38 andâ,ê andγ for the original one Eq. 35. A sample of the numerical results obtained is shown in Fig. 11. It is important to stress that the comparison is made at the intersection of the orbits with . Hence, the abscissas are the number of iterates in (a total of 1789440 in all cases), not the integration time as the intersection times in are not necessarily the same in the two compared problems.
The shown results correspond to 19 different values of ¯ : j α = 0 : 25 : 480. Concerning the differences in semimajor axes, along this first integration year, they are O(10 −2 ); concerning the differences in eccentricity, they are O(10 −5 ); and those for the orientation, they are at most O(10 −4 ). Taking into account that ε = 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, exceeds the predictions done by classical theorems that are summarized in Theorem 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 eff ), one can find a nonnegligible 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.

Conclusions
The coupled attitude and orbit motion is known to be a problem with two characteristic timescales. 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 quasirhombic 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 heliostability properties-the ratio between timescales, ε, can be found explicitly. An adequate time and phase scaling that depend on ε 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 ε to 0 with equations written in the fast timescale. In the specific case of the simplified QRP, it is a Hamiltonian system of one degree of freedom with C 0 equations of motion that resembles a pendulum. The Sun-pointing attitude is a stable equilibrium provided an adequate centerof-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 neighborhood 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 relies mainly on two independent requirements: smallness of ε 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 are that of a constructible structure, the value of ε 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, 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 toward 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 Felicetti et al. (2017). 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 makes them suitable for cargo transportation missions and interplanetary transfers; see, e.g., Mengali and Quarta (2007). Council (ERC) under the European Unions Horizon 2020 research and innovation program as part of project COMPASS (Grant Agreement No 679086). The authors acknowledge the use of the Milkyway High Performance Computing Facility, and the 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/∼ publications. Fruitful conversations with I. Gkolias, M. Lara and C. Simó are also acknowledged.

Compliance with ethical standards
Conflicts of interest The authors declare that there are no conflicts of interest regarding the publication of this paper.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

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· 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· refer to the adimensional analogues. Let L and T denote the length and time units, measured in km and sec, respectively. Finally, 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 The first component of the dimensional Kepler problem readsẍ = −μx/r 3 , wherer = x 2 +ỹ 2 . From Eq. 39,r = Lr, and also, if we choose T 2 L 3 μ = 1, or, equivalently T = L 3 μ .
It is clear that for this choice the second component of the adimensional Kepler problem reads y = −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 Sect. 4. Before scaling the equations of motion, it is convenient to shift the angleφ: letφ =φ −λ, whereλ is the argument of latitude the Sun (i.e., the angle between the x axis and its position vector in Earth-centered coordinates, measured counterclockwise) in its apparent motion around the Earth.
In the adimensional variables, the attitude equations of motion read φ = A s p SR k 1,1 L 3 2Cμ(m b + m s ) M 1 (φ) + 3 D(α, d) C 1 r 3 sin(2 arctan(y/x) − 2(λ − φ)), 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 ⎧ ⎪ ⎪ ⎨ where a x and a y are given in Eq. 8 putting ϕ = φ + λ.