Deorbiting spacecraft with passively stabilised attitude using a simplified quasi-rhombic-pyramid sail

Abstract This paper studies deorbiting using an analogue to the quasi-rhombic-pyramid concept for planar motion. The focus is on maintaining a stable (meaning oscillatory) attitude close to the direction of the velocity of the spacecraft relative to the atmosphere. The study consists of a massive computation of deorbit times chosen in a region of the phase space where atmospheric drag plays a leading role. Here, no damping effects are considered. Thus, any passive stabilisation observed is either due to solar radiation pressure or atmospheric drag. The results show that such stable deorbiting is feasible up to a threshold that depends upon the physical parameters of the sail. This threshold is around 500 km of altitude. Stable deorbiting is also shown to reduce the unpredictability that appears due to tumbling.


Introduction
The growing population of space debris in the last decades has increased the interest and need for end-of-life disposal devices and strategies. The recent literature has focused on the idea of carrying an on-board deployable surface that increases the area-to-mass ratio of spacecraft to enhance either the effects of Solar Radiation Pressure (SRP) or the effects of atmospheric drag, see, e.g. Colombo et al. (2012,  and Lü cking et al. (2013). The enhancements of the SRP acceleration for deorbiting apply for high altitude orbits, whereas atmospheric drag is useful for lower altitudes, but in either case the maximisation of their effect on the orbit for any purpose requires attitude control.
The main approaches that exploit the SRP acceleration to deorbit spacecraft can be organized in two groups according to their attitude control requirements. On the one hand, there is the passive approach, that consists of increasing the eccentricity of the orbit, which results in a decrease of the perigee radius. Here the term passive is used in the sense that the attitude of the spacecraft is fixed with respect to the direction of sunlight, see . On the other hand, there is the active approach, that consists of reducing the semi-major axis of the orbit by maximising the cross area when the spacecraft travels towards the Sun and minimising it when the spacecraft travels away from the Sun, see Borja and Tun (2006). This approach requires to change the attitude of the spacecraft twice per revolution around the Earth. These attitude control requirements can be reduced by only changing from maximal to minimal exposition to sunlight or vice versa according to the sign of the secular and long-term evolution of the eccentricity. This control consists of minimizing the cross area when this term is negative and maximising the cross area when it is positive. This reduces the attitude control manoeuvres from twice per revolution to approximately twice per year. This approach is referred to as modulating, and falls inside the active category, see Colombo . These approaches were compared in Colombo and de Bras de Fer (2016), where the passive strategy was found to be the most effective technique for short deorbiting times, while active approaches were of interest for longer deorbiting times.
The main drawback of the approaches that exploit the SRP is that attitude control is required, either to periodically change the attitude or to maintain the attitude fixed with respect to a prescribed direction for some time or along the whole motion. In this paper we are interested in the shapes of the sail that reduce the attitude control requirements to maintain a specific attitude, on average. That is, shapes of the sail such that if the spacecraft is adequately oriented towards a specific direction, this is orientation is either sustained or the attitude oscillates around this direction, in a passive way. Note that here the word passive is used in the sense that attitude stability is solely because of the external torques that act on the spacecraft. This is the sense in which this word will be used from now on, unless we refer to the passive deorbiting approach. In  the idea of using a Quasi-Rhombic Pyramid (QRP) shape for the sail was introduced. Such a shape is meant to cancel, on average, the components of the acceleration that are not along the sunlight direction, and hence it endows the spacecraft with what we refer to as auto-stabilising attitude properties. This requires choosing appropriate initial conditions for the attitude. This structure and variations of it have been studied in recent literature, with works that range from the dynamics of spacecraft with this kind of structures, to applications, see, e.g. , Heiligers and Ceriotti (2017), Felicetti et al. (2016), Felicetti et al. (2017), Miguel and Colombo (2019a) and Miguel and Colombo (2018) or its extended version Miguel and Colombo (2019c).
Atmospheric drag is a non-conservative force and, as such, it leads to the decrease of the orbital altitude of spacecraft, regardless of their attitude. Its effect becomes stronger if the atmospheric density is larger, and it is the dominating effect at very low orbital altitudes, namely below 500 km of altitude. Atmospheric drag can be exploited either to control the spacecraft's orbit, see e.g. Virgili-Llop et al. (2014), Virgili et al. (2015) or its attitude, see Kumar et al. (1995), Kumar et al. (1996), Psiaki (2004) and Virgili-Llop et al. (2019). These cited works deal with passive aerostabilisation or passive attitude stabilisation using atmospheric drag forces: If the atmospheric drag is the dominant effect, and if spacecraft with a QRP shaped sail, or shaped like some of its variations, are oriented towards the velocity vector, one expects their attitude to be passively maintained close to the velocity vector, as the sail resembles a shuttlecock.
This paper is a contribution to the understanding of the dynamics and performance of QRP shaped (or QRP-like shaped) sails for deorbiting purposes, when the spacecraft is affected by the SRP effect and atmospheric drag.
A simplified setting where the orbital dynamics are assumed to happen in the ecliptic plane is considered. The shape of the sail taken into consideration is a modification of the QRP shape that consists of two rectangular panels attached to each other. If adequate initial conditions for the attitude are chosen, spacecraft with such a sail would ideally only rotate around an axis that would be perpendicular to the orbital (ecliptic) plane. Thus, the orbital dynamics would remain planar even if they were coupled with the attitude dynamics.

Brief literature review
The class of spacecraft considered here has been already studied in the recent literature. These works focus on one single effect, either the SRP or the atmospheric drag, and deal with the passive stability enhancements of the variations of the QRP shape with respect to either the sunlight direction (for the SRP) or the velocity vector (for the atmospheric drag).
Concerning the SRP, the attitude and orbital dynamics of the full QRP shape sail introduced in  were studied in Felicetti et al. (2016). Special emphasis was put on the preservation of the attitude of the spacecraft taking into account eclipses and other disturbance torques like gravity gradient. Eclipses were shown to affect substantially the auto-stabilising properties of the QRP shaped sail, and a moderate spin and ring dampers were justified to be feasible solutions to maintain heliostability.
The problem considered here is that already tackled in Miguel and Colombo (2019b), where the effects of eclipses were not taken into account. For a large set of initial conditions with initial a 0 ¼ R È þ 1000 km and different values of e 0 -0 the attitude dynamics were studied in the case where the initial conditions for the attitude were chosen to be u 0 ¼ k 0 . That is, the spacecraft was considered to be initially oriented exactly towards the direction of sunlight. The attitude stability around the direction of sunlight was shown to be eventually lost for all initial conditions while deorbiting. Taking into account that all initial conditions were chosen in the Low Earth Orbit (LEO) regime, and according to Felicetti et al. (2016), if eclipses were taken into account, the oscillations around the direction of sunlight without the aid of any damping or spin are also expected to be lost.
Concerning the atmospheric drag, the ''shuttercock" design proposed in Psiaki (2004) to exploit atmospheric drag as a passive stabilising effect was studied in Rawashdeh et al. (2009). The feasibility of attitude stabilisation by adding damping effects was tested by assuming that the orbit was fixed and always circular, and that the attitude dynamics only occured around an axis perpendicular to the orbital one. The provided results suggest that such stabilisation is possible for altitudes below 450 km and that a value of the aperture angle between panels around 100 was the most efficient choice to reach aerodynamic stabilisation.
The effectiveness of a shuttercock-like sail for passive aero-stable deorbiting was justified in Roberts and Harkness (2007). The considered sail structure in this contribution was considered to have a cone shape to ease the analysis as, in practice, constructible structures would be pyramidal, as the QRP shape. Damping torques were considered to stabilise the structure during deorbiting. For some test initial values of the altitude, the optimal value of the aperture angle between panels was discussed.
The effect of atmospheric drag as self-stabilising effect was also studied in Miguel and Colombo (2019c) in the context of deorbiting, considering the atmospheric drag and gravity gradient torques. The sensitivity analysis with respect to the aperture angle and the center of masscenter of pressure offset shows that larger aperture angles lead to shorter deorbiting times, while the center of masscenter of pressure offset does not affect the deorbiting times.
The joint effect of SRP and drag is expected to lead to complex attitude dynamics. In the orbital region where atmospheric drag plays a role, this becomes increasingly relevant as the semi-major axis of the orbit decreases, on average, and for lower altitudes the air density becomes larger. This makes the attitude dynamics and their effect on the orbit dynamics highly unpredictable. The main goal of this article is to study whether the considered analogue of the QRP shape can maintain a stable (meaning oscillatory around some direction) attitude along the orbit as it naturally deorbits, where is it possible, and how the deorbit time is affected by these stability properties. No damping torques are added; hence any stabilising effect is exclusively due to the considered forces.
The paper is organised as follows. First, in Section 2 the considered family of sails is described, together with the main physical parameters these depend upon. In Section 3 the model of the attitude and orbit evolution is provided, and a common framework that allows to study the torques due to SRP and drag is introduced. This common framework is used in Section 4 to perform an analysis of the attitude dynamics by describing the phase space of the attitude variables. This is used to define the concepts of stable attitude and stable attitude deorbiting. In Section 5 the numerical study is explained and justified, and from it the main results on auto-stabilisation are drawn. The contribution finishes in Section 6 with conclusions and future lines of research that emerge from this paper.

Geometry of the sail structure
The class of spacecraft considered here consists of a payload attached to an already deployed deorbiting device. This device is a sail that consists of two equal rectangular panels, P þ and P À , of width w and height h that are attached to each other on one of the h-long sides. Consider a body frame F b whose origin is the Centre of Mass (CoM) of the spacecraft and whose coordinates are n; m and f. Denote the vectors of the basis of F b by i n ; i m and i f . A sketch of the considered class of spacecraft in the frame F b can be seen in Fig. 1. Plot (a) is a 3D view of the sail structure; and plot (b) is the projection of the spacecraft onto the ff ¼ 0g plane where all its relevant parameters and elements can be visualised. The sail panel P þ is contained in fm > 0g and P À is contained in fm < 0g, and they form an angle 2a 2 ð0; pÞ. The CoM of the sail structure (and hence the Centre of Pressure, CoP) is depicted as a (solid) blue dot. The normal vector to the surface of the panel P AE is denoted by n AE . The vectors n þ and n À are perpendicular to i f , and they are oriented as depicted in Fig. 1(b). The sail structure is symmetric with respect to the axis n. The CoM of the payload is assumed to be on this n axis. This is depicted as a black square in Fig. 1(b). The basis vector i n is parallel to the vector n þ þ n À . The triad is completed by This is a simplification of the QRP shape such that the attitude dynamics only occur around i f , provided adequate initial conditions for the attitude are chosen and if the external forces act in a plane that is perpendicular to i f . In this paper the orbital dynamics are considered to happen in the ecliptic plane, which is assumed to contain the considered external forces. The initial conditions for the attitude are chosen so that the vector i f is always perpendicular to the ecliptic plane.
Apart from the size, mass and reflectivity parameters, that characterise the area-to-mass ratio and the size of the acceleration perturbation due to either the SRP or drag, the main physical descriptive parameters of the structure are the aperture angle a 2 ð0; p=2Þ, that is measured either in radians or degrees; and the CoM -CoP offset d, that is measured in meters. The parameter d is defined as the (signed) distance between the CoM of the payload and the CoM of the sail structure. In Fig. 1(b) an example with d < 0 is shown. If d ¼ 0, both the CoM of the payload and the CoP are at the origin of F b , see the sketch in Fig. 3. The CoM mass of the payload of a spacecraft with d > 0 would have a larger n component than the CoP.
Denote the mass of the whole sail structure by m s and the mass of the payload by m b . Assume that the payload is cubic with side length s b . If the moments of inertia around the axes n; m and f are denoted by A; B and C, respectively, in Miguel and Colombo (2018) it is justified that The moments of inertia A and B are not provided as their explicit expressions are not needed in this contribution. However, their difference D, that in this particular case is a part of C, is a factor in the gravity gradient torque. This is used in Section 3.3.

Model
The model used in this paper is a coupled system of orbit and attitude differential equations. The expressions of the accelerations and the torques due to the SRP and drag were derived in Miguel and Colombo (2018), and the reader is referred to this reference for further details. The orbit dynamics are that of the J 2 problem -the motion of a mass-less particle around an oblate planet keeping only the second-degree zonal harmonic of the expression of the perturbing potential-perturbed by the SRP and atmospheric drag forces. This is written in an Earth centred inertial frame, that is denoted by F I . The coordinates of F I are denoted by x; y and z, and the vectors of its basis are denoted by i x ; i y and i z . The vector i x points towards an arbitrarily chosen direction. As the problem is planar, the vector i z is parallel to i f and perpendicular to the orbital plane, that is assumed to be fz ¼ 0g. The triad is completed by choosing i y ¼ i z Â i x . In this frame r sc denotes the Earthspacecraft vector. The modulus and the direction of r sc are denoted by r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi x 2 þ y 2 p and u sc , respectively, so r sc ¼ ru sc . Similarly, r S ; r S ; u S denote the Earth-Sun vector, its modulus and direction, so r S ¼ r S u S ; and v rel ; v rel ; u rel denote the relative velocity vector, its modulus and direction, so v rel ¼ v rel u rel . Let k denote the angle defined by the position of the Sun with respect to i x . The apparent motion of the Sun around the Earth is considered to be circular, so k ¼ n S t þ k 0 , where n S is the mean motion of the apparent motion of the Sun, and k 0 is the initial phase. See Fig. 2.
The orbital dynamics are coupled with the attitude dynamics. The equations of motion of the attitude dynamics can be written with sole dependence on one single Euler angle u 2 ½0; 2pÞ rad that describes the orientation of the body frame F b with respect to F I as where M is the sum of the torques taken into consideration, and C is the moment of inertia along the f axis, whose expression is provided in Eqs.
(1). In this contribution the backside of the sail is also assumed to produce torque and acceleration due to both SRP and drag. Here the self-shadowing of the panels is taken into account, and the shadow that the payload may cast on any of the panels is neglected. Let d denote the flight path angle, see Fig. 2. The torques and accelerations due to the SRP and drag depend upon the orientation of the spacecraft with respect to the directions in which these forces act. In the explicit formulas, this dependence is expressed as the dependence on the value of u À k for the SRP and of u À d for drag. To simplify the notation let us use here the same symbol / to denote the relative orientations, either u À k for SRP or u À d for drag. A sketch of the frames F I and F b in the SRP case can be seen in Fig. 3. In this figure the angles u; k and / ¼ u À k are depicted. For the drag case, the sketch would be analogous, but replacing k for d and u S for u rel . Table 1 contains the values of the area A exp that is exposed to either the SRP or the air flow for each possible relative orientation.
As only rotations around the f axis are considered, the shadow that one panel casts on the other is a rectangle of height h and width w À w 0 , where w 0 is the width of the exposed area. Hence, the area that is partially exposed due to self-shadowing is hw 0 . This happens in the cases / 2 ðÀp þ a; Àp=2 and / 2 ½p=2; p À aÞ in Table 1. The width of the exposed area w 0 can be written in closed form and can be simplified to w 0 ð/Þ ¼ 2w cos / sin a sinða À /Þ : Fig. 4 is a sketch of the projection of a spacecraft with d ¼ 0 and / 2 ½p=2; p À aÞ onto ff ¼ 0g. For this value of /, the spacecraft shadows itself. In the projection, denote by q AE the endpoint of the panel P AE . Finding the distance w 0 reduces to finding the point q 0 , the intersection between the line that is parallel to u S or u rel that passes through the endpoint of the panel that is completely exposed (in Fig. 4, q À ) and the partially shadowed panel (in Fig. 4, P þ ).

Force models
The accelerations due to the SRP and atmospheric drag are considered. Following McInnes (1999), Montenbruck and Gill (2005) the force due to the SRP exerted on the panel P AE whose exposed area is A exp is modelled assuming that a fraction of the incident radiation g 2 ð0; 1Þ is specularly reflected and the rest (a fraction 1 À g) is absorbed. The parameter g is referred to as the reflectivity parameter 1 . The force due to the specularly reflected radiation Fig. 2. Sketch of the Earth-centred inertial frame F I . Fig. 3. Sketch of the frames F I and F b and the angles u, the orientation of F b with respect to F I ; k, position of the Sun in F I ; and / ¼ u À k, the orientation of the F b relative to the direction of sunlight. Table 1 Exposed area as a function of the relative orientations.
hw hw 0 ðÀ/Þ ðÀp=2; Àa hw 0 ½Àa; a hw hw ða; p=2Þ 0 hw ½p=2; p À aÞ hw 0 ð/Þ hw ½p À a; pÞ hw hw contributes in the direction of n AE while the absorbed radiation contributes in the direction u S . Hence the force due to the SRP exerted on panel P AE can be written as where p SR ¼ 4:56 Â 10 À6 N=m 2 is the solar pressure at 1 AU. The parameter p SR is assumed to be constant. Note that the case g ¼ 1 corresponds to a perfectly reflecting solar sail, and the case g ¼ 0 corresponds to panels for which all radiation is absorbed. The parameter g typically ranges between 0:2 and 0:9, see Montenbruck and Gill (2005). The value g % 0:88 is estimated for solar sails with a highly reflective aluminium-coated side, see Mengali et al. (2007), Heaton and Artusio-Glimpse (2015) and Heaton et al. (2017). For the numerical tests, the value g ¼ 0:8 is used.
The effect of eclipses is taken into account using an Earth cylindrical shadow model without penumbra. As the dynamics are planar, the spacecraft is considered to be in eclipse when it is inside the strip of shadow cast by the Earth.
Concerning the atmospheric drag, here it is assumed that the kinetic energy of all the air particles that hit the surface of P AE is totally transferred to the spacecraft, and this gives rise to a force in the direction of u rel . Under this assumption the force due to atmospheric drag exerted on the panel P AE can be written as where q is the atmospheric density (measured in kg/m 3 ), and C D 2 ð1:5; 2:5Þ is an empirically determined dimensionless drag coefficient, see Montenbruck and Gill (2005). For simulations, the atmospheric density q is retrieved using the Exponential Atmospheric Model, see Vallado and McClain (2013), page 510. The used value for the drag coefficient is C D ¼ 2:2, as the sail consists of 2 flat plates, see Vallado and McClain (2013), page 551.
The exposed area A exp in Table 1 is a factor in both force models in Eqs. (3) and (4). As A exp is piece-wise smooth as a function of /, so are the accelerations and torques due to the SRP and atmospheric drag. These torques and accelerations have / ¼ AEp Ç a; / ¼ AEp=2 and / ¼ AEa as switching manifolds, where the equations lose differentiability with respect to /, but not continuity.

Orbit dynamics
Building on previous contributions on the usage of the SRP to design end-of-life disposals, see Lü cking et al.  (2016), the J 2 problem perturbed by the accelerations due to the SRP and the atmospheric drag is considered. In Cartesian coordinates, the equations of motion read where l È ¼ 3:986 Â 10 14 m 3 /s 2 is the gravitational standard parameter of the Earth, R È ¼ 6378:1 km is the mean Earth radius, and a SRP and a drag refer to the accelerations due to the SRP and drag. These accelerations are obtained as the sums of the acceleration on each panel.

Attitude dynamics
The attitude dynamics are governed by the set of 2 firstorder differential equations given in Eq. (2), where The summands are the torques due to the SRP, drag and Gravity Gradient (GG). In Miguel and Colombo (2018) it is justified that the torque due to the SRP can be written as The summand M AE SRP is the torque due to the panel P AE , and it is obtained by computing The function M AE 0 ðw; gÞ reads M AE 0 ðw; gÞ ¼ À 1 2 k 1;1 ðgÞ sinð2wÞ Çk 2;0 ðgÞ cos 2 w Ç k 0;2 ðgÞ sin 2 w; where Fig. 4. Sketch of the top view of a spacecraft with d ¼ 0 that shadows itself. In this sketch k or d is 180 , and / 2 ½p=2; p À aÞ.
2 Recall that the CoM of the spacecraft is at the origin of F b k 2;0 ðgÞ ¼ sin 2 a½4dgm b cos a þwðm b þ m s Þð1 À g cosð2aÞÞ; and ð8bÞ The advantage of this formulation is twofold: on the one hand, the provided explicit expressions allow to see that for j/j < a, that is, when the spacecraft is oriented close to the sunlight direction, the equations of motion of the attitude dynamics can be arranged as a perturbed pendulum close to an equilibrium point that is at / ¼ 0. The stability of this equilibrium point depends solely on the physical parameters of the spacecraft and can be determined using a closed-form formula. This is discussed in Section 4. On the other hand the torque due to the panel P AE due to the atmospheric drag can be expressed using the functions M þ 0 and M À 0 in Eq. (7), k 1;1 ; k 2;0 and k 0;2 in Eqs. (8), but with d in the place of k and g ¼ 0. Namely the torque due to drag can be written as These expressions are an arrangement of the formulas obtained using the relations M AE drag ¼ r b;AE Â F AE drag . It is important to remark that the force due to drag F drag does not depend on g, and the fact that M drag can be written using the functions M þ 0 and M À 0 for g ¼ 0 is a convenient interpretation that is a consequence of the chosen force models in Eqs. (3) and (4): all kinetic energy of the air particles that impact the sail panels being transferred to the spacecraft is interpreted as the case g ¼ 0 where all incoming solar radiation is absorbed.
This interpretation is useful as the attitude dynamics, and in particular the stability of the equilibrium point at / ¼ 0, can be theoretically handled by only studying the dynamical system for g 2 ½0; 1Þ, where w is an adequately translated and scaled angle. This is the content of Section 4. Finally, concerning the torque due to the gravity gradient, as the payload is assumed to be symmetric, it reads recall Eqs. (1) for the definitions of C and D.

Description of the attitude phase space
For adequate values of the physical parameters, the considered spacecraft is expected to exhibit auto-stabilisation properties close to the sunlight direction if the SRP is the dominant effect, or close to the relative velocity vector if atmospheric drag dominates instead. In case both effects have similar magnitudes, one expects that, for most initial conditions, the motion would be complex and highly unpredictable. This is because both effects affect the dynamics in different directions that depend on time and on the position of the spacecraft in its orbit. The only theoretically approachable scenario is that in which either the SRP or drag dominates the other effects taken into consideration. In this case, the equations of motion of the attitude variables can be written as where I is either the SRP or drag, and the other two remaining effects (drag or SRP, depending of what I refers to, and gravity gradient) are expected to be smaller. As both the SRP and drag effects can be written as some constants times the function M 0 ðw; gÞ in Eq. (10b), this suggests to study the dynamical system w Prime ¼ M 0 ðw; gÞ as defined in Eq. 10a. In Eq. (10a), ð0Þ ¼ d=ds where s is a scaled time variable and w is an adequately scaled and translated angle whose origin is either the sun-pointing direction (in the SRP-dominated case) or the relative velocity vector (in the drag-dominated case). This reduction is fully justified in Miguel and Colombo (2019a) in the SRP-dominated case. In case the atmospheric drag dominates, the change of variables that would lead to the simplified form of Eq. (10a) would be dependent on the altitude, as M drag in Eq. (9) depends on the air density q.
The points E 0 : w ¼ 0; w 0 ¼ 0 and E 1 : w ¼ Àp; w 0 ¼ 0 are always equilibria of the system in Eq. (10a), regardless of the values of the physical parameters the system depends upon. The equilibrium E 0 represents the direction of the sunlight and the direction of the relative velocity vector. In Miguel and Colombo (2018) it was proven that a necessary condition for the stability of E 0 (and hence the instability of E 1 due to a change of sign) is that d > d min , where Kða; gÞ ¼ g cosð3aÞ À cos a 2g cosð2aÞ þ g þ 1 : Note that g 2 ½0; 1Þ; a 2 ð0; p=2Þ rad, and the rest of the quantities involved in Eq. (12) are positive. Thus, d min < 0. In particular, if in F b the CoM has a larger n component than the CoP (i.e., d > 0), E 0 is a stable equilibrium of Eq. (10a). The existence of other equilibria depends upon the values of all physical parameters. This problem has to be evaluated numerically. Eq. (10a) written as an order 1 ordinary differential equation reads Hence, the set of equilibria of Eq. 10a are the pairs ðw; WÞ ¼ ðw Ã ; 0Þ, where M 0 ðw Ã ; gÞ ¼ 0.
To illustrate the qualitative description of the phase space of Eq. 10a the following values for the physical parameters have been chosen. Concerning the spacecraft size and mass, h ¼ w ¼ 9:20 m, m s ¼ 3:6 kg and m b ¼ 100 kg. This corresponds to a constructible structure according to the guidelines of Dalla Vedova et al. (2018). Concerning the reflectivity parameter, the values g ¼ 0 (fully absorptive case, interpreted as drag case) and g ¼ 0:8 (SRP case) have been chosen. All parameters except for a and d are fixed, so after these choices the parameter space is 2-dimensional.
In what follows the condition d > d min is assumed. That is, E 0 is assumed to be a stable equilibrium point. The region of the ða; dÞ-parameter space that corresponds to pairs ða; dÞ for which d > d min can be divided into two sub-regions according to the number of equilibria of the system in Eq. (10a). As M 0 is anti-symmetric, i.e. M 0 ðÀwÞ ¼ ÀM 0 ðwÞ, if additional equilibria exist, they have to come in pairs. The numerical evidence provided here suggests that there are either 2, 4 or 6 equilibria. These three cases are exemplified in Fig. 5 for a ¼ 45 and for 3 adequately chosen values of d. A description of each of these 3 cases follows.
(a) Only E 0 and E 1 are equilibria. The point E 0 is stable and E 1 is unstable. This is true also in the other two cases. (b) The equilibria are E 0 ; E 1 , and two other points that appear at w ¼ AEw bif ða; d; gÞ due to a bifurcation. These occur in the points where M 0 is tangent to The equilibria are E 0 ; E 1 , and two pairs of stableunstable equilibria that bifurcate from the previous bifurcation points in (b).
The bifurcation is of the saddle-centre type, that is of codimension one so it can be locally explained by one single parameter. As all parameters except for a and d have been fixed, one expects to be able to find a curve of bifurcations in the ða; dÞ plane that separates the set of parameters for which there are 2 and 6 equilibria. The bifurcation curve would be the set of pairs ða; dÞ for which there are 4 equilibria.
This is the content of Fig. 6, where (d) is the bifurcation diagram for g ¼ 0 and (e) is the bifurcation diagram for  g ¼ 0:8. In (d) and (e) the plane is separated into 3 regions, which are limited by red (dashed) lines, that represent the condition d ¼ d min (recall Eq. 12), and blue (solid) lines, that are the curves of bifurcation points. For pairs ða; dÞ chosen below these red (dashed) lines, E 0 is an unstable equilibrium. This contribution focuses on pairs ða; dÞ above the curves d ¼ d min , as for these values the equilibrium E 0 is stable. For pairs ða; dÞ that lie on the blue (solid) lines, Eq. (10a) has 4 equilibria, as in the case (b) in the previous enumeration, and the attitude phase space is expected to be qualitatively as shown in Fig. 5 The value d 2 is a numerically approximated bifurcation point.
It is worth noting that this is the qualitative description for all values of a. This is true as the red (dashed) line lies below the blue (solid) line. Thus, for any value of a, one can always find 3 values of d for which the phase portrait of the attitude variables w; W is qualitatively the same as the examples shown in Fig. 5.
These remarks should be taken into account when designing any control strategy that includes the change of aperture angle a. In the examples shown in Fig. 5, the bifurcations occur inside the main region of stability. However, as this system is conservative, it may happen that this stability region separates in more than one connected component and hence the overall global dynamics would drastically change, see Miguel et al. (2013).
Note that this numerical study can also be performed by considering a variable parameter g. This can be of interest for control strategies based on electrochromic properties of the sail panels, see .

Area-to-mass ratio concept
In this paper the attitude dynamics are considered to evolve freely, subject to the torques of the SRP and drag forces, and gravity gradient. The description of the attitude dynamics performed above in this same section gives a first approximation of the expected behaviour in the case where either the SRP or drag dominates, excluding the other effects. Once these are added, it is not possible to provide a complete description of the long-term evolution of the attitude as the dynamics are strongly dependent on the initial conditions for the orbit and attitude. Also, one has to take into account that in the orbital region we are interested in, the effect of eclipses cancels the SRP acceleration and torque for a non-negligible fraction of each revolution around the Earth, and the magnitude of the atmospheric drag acceleration and torque depends on the orbital altitude.
Assume that a and d are chosen so that, as a first approximation, the direction of sunlight for the SRP torque and the relative velocity vector direction for the atmospheric drag are stable equilibria of the attitude dynamics. That is, d > d min for g ¼ 0 and g ¼ 0:8, recall Eq. (12). Even if the spacecraft was initially oriented towards either of these directions, the other effects can lead to a complex attitude behaviour that pulls the spacecraft out of this stable state.
The magnitude of the SRP and drag effects is measured using the area-to-mass ratio: the ratio between the area that is exposed to either sunlight or air flow and the overall mass of the spacecraft. This is a faithful measure of the magnitude of these two effects provided the attitude is fixed with respect to either u S or u rel . In the case where the attitude evolves freely, the concept of area-to-mass ratio does no longer make sense as in general, the orientation of the spacecraft with respect to the directions of the forces changes continuously. It is only in very specific cases where one can find a quantity that is analogous to the area-tomass ratio of the spacecraft that can explain the effect of the attitude on the orbital dynamics on the long term. From a mathematical standpoint, this can only be found by means of averaging techniques, that in the class of spacecraft under consideration can only be applied rigorously in case the attitude dynamics are oscillatory around a stable direction, as in the studies of , Heiligers and Ceriotti (2017) and Miguel and Colombo (2019a).
The class of spacecraft considered in Miguel and Colombo (2019a) is exactly the same as that considered here. In this reference, it was found that if the quotient between the area of the panels A ¼ hw and the overall mass of the spacecraft M ¼ m b þ m s ; A=M was large enough, and the spacecraft was initially oriented close enough to the attitude stable equilibrium, one could find an explicit formula for an effective area, that was denoted by A eff : the area of a flat sail that would produce the same force as the simplified QRP shape considered here, on average. The area A eff depends on A, the aperture angle a, the reflectivity g; and on the amplitude of the oscillations around the Earth-Sun vector. This means that each amplitude of oscillation around the stable equilibrium in ð0; 0Þ in Fig. 5 gives rise to a different value of A eff .
In the particular case in which the sail structure is oriented exactly in the direction of the sunlight, i.e. u ¼ k, the effective area reads It is important to remark that far from being theoretical minutia, the fact that the necessary conditions meet so that averaging results can be applied indicates that more faithful predictions can be performed, in the sense that close initial conditions give rise to qualitatively and quantitatively similar results. This will be exemplified in the numerical results in Section 5.

Attitude stability concept
In a realistic situation one cannot expect the attitude of spacecraft to be exactly fixed with respect to some direction. In the best possible scenario, the attitude dynamics would be oscillatory. The direction around which the oscillations would occur, in turn, depends upon the leading effect the spacecraft is subject to, and different effects may play leading roles in different regions of the phase space.
Therefore, the concept of attitude stability has to be understood in a practical way, and always linked to a specific effect. In the dynamical systems lingo, stability is referred to as being practical when the motion is oscillatory for, at least, a prescribed amount of time. The length of this time interval is chosen accordingly to the studied phenomenon. This concept is used when the asymptotic behaviour is not analytically predictable.
In the context of this paper, the attitude stability can be linked either to the SRP, in case the attitude dynamics are oscillatory around u S ; or linked to the atmospheric drag, in case the oscillations occur around u rel .
As in Section 3, denote by / the relative orientation with respect to the direction of interest, either u À k for SRP or u À d for drag. A spacecraft is considered to have stable attitude in some time interval ½t 0 ; t f if for all points of its trajectory during this time interval, the attitude dynamics are librational around the direction of interest. This is equivalent to the existence of an upper bound K; 0 < K < p rad, such that j/j < K in ½t 0 ; t f .
The concept of stable attitude deorbiting hence refers to having a stable attitude from the beginning of the tracked motion until the situation in which the spacecraft is considered to be deorbited.

Numerical study
In this section we study numerically under which conditions the modified QRP shape considered here maintains a stable attitude around u rel along the motion, until it is considered to be deorbited. To have a general perspective, the set of initial conditions for the orbit and the attitude have to be appropriately chosen.

Initial conditions
The performance of the modified QRP shape considered here as a deorbiting device is measured by studying how the deorbit time behaves among a large number of initial conditions, in relation to the evolution of the attitude along the motion.
The deorbit time is the amount of time that the spacecraft spends from an initial condition until it reaches an altitude of 120 km. When a probe reaches altitudes close to the Kármán line, the provided model is no longer useful as thermal effects have to be taken into account, see Roberts and Harkness (2007).Initial conditions for the orbit.
To ensure a reduction of the semi-major axis of the orbit, the spacecraft has to be affected by atmospheric drag from the beginning of the motion. The numerical study in Miguel and Colombo (2018) suggests that aerostability with the class of spacecraft considered here is only feasible below 850 km of altitude, as above this threshold the air density is too small. On the other hand, as the effect of eclipses is taken into account, the effect of the SRP depends strongly on the initial orientation of the spacecraft with respect to the Sun. To account for that, we considered the following orbit initial conditions: k 0 ¼ 90 , and a 0 ¼ 250 þ 10ikm; i 2 f1; 2; . . . ; 60g; ð14aÞ x 0 ¼ 6 j; j 2 f1; 2; . . . ; 60g; ð14cÞ This is a 2D discretization of ða 0 ; x 0 Þ 2 ½250; 850 Â ½0; 360 into 3600 points, in a 60 times 60 equispaced grid. The value e 0 ¼ 0 has been chosen as for larger initial values of the eccentricity, the altitude of the orbit exhibits larger oscillations. Hence, along the motion, the spacecraft would enter and exit regions with higher air density q and in this case the spacecraft would eventually tumble. Initial conditions for the attitude. The only feasible scenario in which there are self-stabilising attitude properties throughout the whole trajectory is that where the spacecraft is initially pointing towards u rel , see Miguel and Colombo (2019b). For all initial conditions for the orbit, we have chosen

Physical parameters
The size of the sail panels is a critical physical parameter for the prediction of deorbiting times, as for a fixed mass of the spacecraft, the larger the cross area is, the faster is the deorbiting. To illustrate how the deorbit time and the attitude stability of the orbits depend upon the size of the spacecraft panels, we have considered three different constructible structures according to the guidelines of Dalla Vedova et al. (2018). We denote them by sc 1 ; sc 2 and sc 3 . Assuming that the two panels P þ and P À are squared, i.e., h ¼ w, for a fixed mass of the payload m b ¼ 100 kg, we have chosen three different values of w imposing that That is, the side length of the panel is chosen in a way that if the bus was attached to it, the area-to-mass ratio of the spacecraft would be 1; 2 or 5 m 2 /kg. The corresponding numerical values are provided in Table 2. The values in Table 2 reflect that mass of the sail panel depends on its size. In Dalla Vedova et al. (2018) the authors provide a method to assess the constructability of a squared solar sail for given values of the payload mass and the area-to-mass ratio, with current technological restrictions. The method gives the mass of the sail as output. This mass is the sum of the masses of the booms and the membrane.
If at the beginning of the motion the QRP shaped sail points exactly towards u rel , the sail would act for some short period of time (before it starts oscillating) as if it had an effective area-to-mass ratio A eff =ðm b þ m s Þ, where A eff is given in Eq. (13). In the three cases of Table 2 these are r 1;eff j u¼d ¼ 1:3623 m 2 =kg; r 2;eff j u¼d ¼ 2:6889 m 2 =kg; r 3;eff j u¼d ¼ 6:5551 m 2 =kg: Note that this effective area-to-mass ratio is larger than the area-to-mass ratio of a single panel. This is counterintuitive fact is analytically justified in Miguel and Colombo (2019a).

Sail configurations
The parameter d does not play a role in the deorbiting times, see Miguel and Colombo (2019c). Moreover, as also justified in this reference, larger a lead to faster deorbiting. To allow for comparison, three different configurations have been chosen: If all physical parameters but d are fixed, d ¼ 0 minimizes the magnitude of the gravity gradient perturbation, see Eqs. (11) and (1b). The physical interpretation of this condition is that the CoM and the CoP are as close as possible.

Methods
Each of the 3600 initial conditions in Eq. (14) has been numerically integrated from t ¼ 0 s until r < R È þ 120 km was reached. This has been done for each of the three spacecraft sizes sc 1 ; sc 2 and sc 3 in Table 2 and for each configuration in Eqs. (17). This has been done considering that the problem has two timescales. To deal with the stiffness of the equations of motion a 2-stage Implicit Runge-Kutta-Gauss method with automatic step size control has been implemented. The stiffness of the equations can lead to very small integration step sizes. Instead of saving all points along the numerical integration it is advisable to study the problem by using passages through a surface of section R (i.e. the Poincaré map on R) instead of actual time. An adequate choice of surface of section for this problem is as orbits are always transversal to the Cartesian axes. See Fig. 2. One iterate of the Poincaré map on R corresponds to one revolution around the Earth. This contribution focuses on the performance of the modified QRP shape sail in terms of its stability towards u rel , not on the actual deorbit times. The results are therefore given in terms of the iterates on R instead of in terms of nominal deorbit times. This eases the understanding of the results. The provided results given in terms of time instead of the iterates are qualitatively the same. Therefore, in what follows, we refer to time and iterates on R indistinctly.
Observable 1: Drag-stable time. The spacecraft is initially oriented towards the velocity vector, that is a stable equilibrium if drag dominates. If the magnitude of the SRP force is large enough, one expects the spacecraft to be pulled out of any oscillatory motion around the velocity vector and to enter a tumbling state. For spacecraft that were tumbling with respect to both u S and u rel , in Miguel and Colombo (2019b) it was observed that in some cases the oscillatory motion around the velocity vector was recovered at low altitudes close to the threshold of 120 km. This is due to the air density becoming larger and hence the atmospheric drag becoming the dominant effect, together with the shuttercock-like shape of the sail. This suggests to measure the time interval in which the spacecraft oscillates around u rel not from the beginning of the motion forward, but backwards from the end of the trajectory, when the spacecraft reaches r < R È þ 120 km. The length of this interval is referred to as drag-stable time, and is denoted by t i;j dÀs . The indices i and j are those in the initial conditions in Eqs. (14). This is measured by counting iterates starting from the last point backwards, and looking for the first iterate on R such that ju À dj > 90 .
Orbits starting at higher altitude require more time to deorbit. To allow for comparison between different altitudes, one needs a relative measure of the drag-stable time. Let us denote the total deorbiting time by t i;j deor , where the indices i; j are again those in Eqs. (14). The chosen observable that is used here to inform about drag-stable time is Observable 2: Average deorbit time. The SRP effect is only added when the spacecraft is in the illuminated phase (i.e., not in eclipse). Hence, its effect is cancelled for some nonnegligible fraction of each revolution around the Earth. This makes the problem strongly dependent on the initial phase, and this is accounted for by choosing different initial values of x 0 þ X 0 . Minimal differences in the initial phase may lead to very different deorbit times. Thus, the study of the results per integrated orbit can only be seen as a particular case study that cannot be necessarily generalised.
A measure of the deorbit time understood solely as a function of the initial altitude of the orbit can be obtained by averaging with respect to all the considered initial phases. This suggests to introduce a second observable: for a fixed altitude with index i in Eq. (14), we consider the average deorbit time, that is denoted by t i deor . It is obtained by computing the average over the 60 considered initial values of x 0 þ X 0 : Observable 3: Dispersion of deorbit times. The dependence of the deorbit time on the initial phasing x 0 þ X 0 can be measured by its dispersion, i.e. the standard deviation. The third observable is, hence,

Results
The results on the performance of the modified QRP shape regarding attitude stability, measured as the fraction of time spent in drag-stable state are shown in Fig. 7. The observable O 1 is displayed as color maps: the first, second and third column correspond to sc 1 ; sc 2 and sc 3 ; and the first, second and third row correspond to each of the three considered aperture angles, a ¼ 30 ; 45 and 60 .
Each pixel is p i;j dÀs (Eq. (19)) where the indices i; j are those in Eq. (14), increasing from left to right and from bottom to top in each subplot of Fig. 7. White pixels indicate that the oscillations of u À d along the whole motion have had an amplitude of at most p=2 rad, so one can consider that drag-stable deorbiting without any control device has been achieved. Black pixels indicate that the dragstable time was zero. Vertical green lines at 475 km for sc 1 , at 485 km for sc 2 and at 500 km for sc 3 are added for reference. The reason why they are added is explained with the results of O 3 . For the studied cases, successful stable deorbiting happens roughly for all orbits up to an altitude threshold h Ã that seems to depend upon the physical parameters of the spacecraft. This is roughly close to 500 km of altitude and increases for larger sail sizes. The stable deorbiting seems to be less likely for altitudes between 500 and 550 km and drag-stable deorbiting without any external control seems to be completely unfeasible above 650 km of altitude.
The results suggest that successful drag-stable deorbiting is feasible from higher altitudes if the sail panels are larger -compare the rightmost column, panels (c), (f) and (i) with the other two-. Of special interest is the case sc 3 with a ¼ 30 , where all non-white pixels are completely above 550 km of altitude. For sc 3 , but with a ¼ 45 and 60 most non-white pixels appear to be above 550 km, which is still better than sc 1 and sc 2 .
The color-maps in Fig. 7 also provide information on the effect of a in the passive attitude stabilisation: smaller amplitude angles a also seem to lead to successful dragstable deorbiting from higher altitudes -compare the top row, (a), (b) and (c) with the two below. This suggests enhanced stability properties for sharper sails.
As these results are relative measures, they have to be interpreted in light of the behaviour of the mean and the dispersion of the deorbiting times, measured as the observables O 2 (Eq. (20)) and O 3 (Eq. (21)).
The results of O 2 are shown in Fig. 8. The top row shows the average deorbiting time as a function of the altitude t i deor in the full range of considered initial altitudes, and the bottom row are details of the pictures right above for the smallest initial altitudes.
The mean deorbiting time seems to depend exponentially on the altitude. This is expected, as the air density increases exponentially when the altitude decreases. For reference on the order of magnitude of the times we are dealing with, we highlight that the average deorbiting times for spacecraft initially at h 0 ¼ 850 km, with aperture angle a ¼ 30 are 3.2646 years for sc 1 , 1.6504 years for sc 2 and 0.67271 years for sc 3 . All the shown results in Fig. 8 give evidence of the fact that a larger aperture angle leads to a faster deorbiting, even for cases where most of the orbit is spent tumbling.
These results explain what was observed above for O 1 .
1. On average, for a fixed overall mass, having a larger panel size leads to a faster deorbiting, and hence to a smaller deorbiting time t i;j deor , see e.g Colombo et al. (2017). The effect of the SRP torque is a disturbance that moves the spacecraft out of the stable equilibria of the drag torque. The SRP effect is not applied uniformly due to the non-symmetric shape of the spacecraft and passages through eclipses. The chance these facts may eventually lead to a tumbling state are reduced if the integration time is smaller. Thus, for a larger area-tomass ratio one expects successful drag-stable deorbiting from higher altitudes to be more feasible.  2. For sharper sails, i.e. for sails with smaller values of a, the effective area A eff Eq. (13) is smaller. This is not only true in the case u ¼ d, but also for non-zero oscillation amplitudes, as shown numerically in Miguel and Colombo (2019a). As the effective area is smaller, so is the size of the SRP and drag torques and accelerations, and this endows the structure with enhanced stability properties, see Roberts and Harkness (2007). Thus, sharper sails allow for successful drag-stable deorbiting from higher altitudes.
The dependence of the deorbiting time on the initial phasing x 0 þ X 0 and on the fact that the trajectory is completely drag-stable or not is measured with O 3 ; r deor , see Eq.
(21). The corresponding results are shown in Fig. 9. For all the three spacecraft and all the chosen values of a, there seems to be a common qualitative behaviour of r deor : For the lowest initial altitudes, namely below 475 km for sc 1 , below 485 km for sc 2 and below 500 km for sc 3 , the standard deviation is strictly below 0.5 iterates. This threshold is indicated in each panel, and also in Fig. 7, in all plots with a vertical green line. For these cases, stable deorbiting leads to highly predictable motion. Note that there are still white pixels on the right hand side of the green lines in Fig. 7. For these cases, even though stable deorbiting is feasible in the sense that the amplitude of the oscillations is below p=2 rad, the increasing oscillatory motion, makes the orbit more sensitive to initial conditions (through the entrances and exits from not illuminated phases). For altitudes above the indicated thresholds, there is an abrupt increase of r deor , that is more prominent for a ¼ 45 and a ¼ 60 . For a ¼ 30 , the increase of r deor seems to be slower. This can be linked to the enhanced stability of sharper configurations. Despite this, all 3 aperture angles reach similar values for the highest altitudes considered.

Discussion and future work
In this paper, a modification of the QRP shaped sails of , that consists of only two panels, has been assessed as a feasible deorbiting device. The shuttercock-like shape endows the structures with attitude stability properties that depend solely on the physical parameters of the system. The study has been performed without considering any damping effect, so any evidence of stability is exclusively due to the external forces taken into consideration. The main physical parameters can be chosen in such a way that the attitude stability with respect to the sunlight direction for the SRP, and with respect to the velocity vector for the atmospheric drag, can be guaranteed at the same time.
Drag-stable deorbiting has been put to the test by initially aligning the spacecraft with the relative velocity vector. The results obtained suggest that successful drag-stable deorbiting is feasible for initially circular orbits for altitudes below a threshold h Ã that depends on the size of the spacecraft and on the aperture angle between the panels. Larger area-to mass ratios and larger aperture angles give rise to faster deorbiting. Smaller aperture angles provide enhanced stability properties that may be of interest to deorbit spacecraft from higher altitudes. The threshold h Ã is around 500 km of altitude in the studied cases and is larger for larger the area-to-mass ratios.
Moreover, attitude stability has been linked to more predictable deorbiting. This requires small amplitude oscillations. In the case the spacecraft tumbles, the shown numerical results suggest that the cannonball approximation can be used to obtain an approximation of the mean effect, yet accurate predictions of deorbit times cannot be expected when using such approximation, see e.g. Rosengren and Scheeres (2013).
The main future lines of research that emerge from this contribution are the following.
1. Sharper shapes give rise to enhanced stability properties, yet its effective cross-area is smaller. This results in slower deorbiting times. In Roberts and Harkness (2007) the authors suggest that one can find an optimal choice of a in the sense that the effective cross-area is maximized while maintaining the enhanced stability properties. A prospective line of research may include tackling this problem numerically and analytically. This is expected to depend both on the aperture angle a and the CoM-CoP displacement d. Fig. 9. Observable O 3 : Standard deviation of the deorbiting time r deor as a function of the initial altitude. Spacecraft sc 1 (first column), sc 2 (second column) and sc 3 (third column).
2. The study of the possibilities of the 3D structure, as proposed in the original reference , has already been initiated in the recent papers Felicetti et al. (2016), Felicetti et al. (2017). Another prospective study may include tackling the justification and measurement of non-negligible regions of stability around the sunlight and velocity directions, the dynamics of oscillatory states as those performed in Miguel and Colombo (2019a), and the effects of coupling SRP and drag and the possibilities for passively stabilised deorbiting.