YORP Effect on Long-Term Rotational Dynamics of Debris in GEO

The Yarkovsky–O’Keefe–Radzievskii–Paddack (YORP) effect describes the torque induced on space objects produced by solar radiation and thermal re-emission. Previous analyses have demonstrated its influence on long-term rotational dynamics of space debris objects in Geostationary Orbit (GEO), where YORP becomes predominant with respect to other external perturbations (e.g., atmospheric drag, gravity gradient, eddy current torque), leading to a wide variety of possible behaviors. The capability of forecasting time windows of slow uniform rotation, if any, would bring significant advantages in operations of Active Debris Removal and on-orbit servicing, especially in the detumbling phase. Also, a non-negligible impact of the End-of-Life configuration, in terms of movable surfaces orientation and center of mass location, could lead to guidelines for future satellites to be easier targets in the disposal phase. In this work, a previously derived semi-analytical tumbling-averaged YORP rotational dynamics model is leveraged. Exploiting an averaged model, computational time is strongly reduced while maintaining sufficient accuracy compared to propagation of Euler’s equations of motion. First, a satellite of the Geostationary Operational Environmental Satellite (GOES) family is analyzed and compared to previous studies to verify the correct implementation of the model. A wider analysis is performed on simple geometric models, such as a box-wing satellite, a 3U CubeSat, and a rocket body. The impact of object size, surface optical properties, and center of mass position on long-term rotational behavior is investigated, providing a general insight into these phenomena with a possible future application to existing objects in GEO.


Introduction
As of April 2022 [1], more than 30.000 objects have been detected in geocentric orbit, with a steadily rising trend. The space traffic itself is also undergoing notable changes, particularly in Low Earth Orbits, fuelled by the miniaturisation of space systems and deployment of large constellations. Although less congested in terms of numbers, satellites in Geostationary Orbit (GEO) and services they provide us with (e.g., telecommunications, weather forecasting, television broadcasting) are also threatened by the problem of space debris, since the nearly absence of atmospheric drag makes impossible to rely on natural orbital decay of defunct objects.
Besides of debris mitigation guidelines, Active Debris Removal (ADR) and on-orbit servicing missions may represent a valid possibility in reducing orbital clutter. However, all these missions must counteract several technical challenges. Basically, a target with large angular momentum may lead to several problems, as an increase in the fuel necessary to circumavigate the object or in the danger of collision and further debris generation. Also, synchronizing with target's motion, chasing a time-varying direction and docking are challenging mission phases, especially when dealing with uncooperative targets which may be "tumbling" (i.e., in a non-principal-axis spin state). These considerations underline that forecasting the target spin state that the spacecraft might encounter is crucial for success.
Nowadays, ground and space-based observations of space objects are used to investigate the spin state of a Resident Space Object (RSO). Light curves reconstructed by photometric data obtained by telescopes [2] can be combined with active techniques such as Satellite Laser Ranging (SLR) [3], Inverse Synthetic Aperture Radar (ISAR) images [4], or Doppler radar measurements [5] to acquire more information about the RSO spin rates and spin axis direction evolution in time.
In the past years, several studies have been dedicated to the analysis of rotational dynamics of space debris objects. ESA's Environmental Satellite (Envisat) is probably one of the most studied in this sense, since it will likely be the target for an ADR mission in the next decades. In 2017, Lin et al. [6] found its spin axis to be stable in the orbital reference frame by optical observation data. Also, a secular increasing in the spin period was observed and justified by the combined effect of gravity-gradient and eddy current torques, dominant perturbations in lower orbits.
Decommissioned in 2006, TOPEX/Poseidon was observed by SLR in 2017 [7] and its inertial spin rate had increased from 0 to almost 6 rotations per minute over the period of 11 years. Among other perturbations, the Solar Radiation Pressure (SRP) is supposed to be the main factor responsible for this gain in rotational energy. When moving to GEO, SRP is often found to be the most significant external torque, especially in the case of objects equipped with wide solar panels. Earl et al. [8] provided a comparison among photometric observations of different box-wing satellites in GEO, namely Solidaridad 1, Telstar 401, EchoStar 2, and HGS-1. Although similar geometry and orbit, the RSOs were found to spin with unique spin periods with respect to one another. Secular spin period variations were found with different amplitude and timescale among the different satellites, sometimes in a quasi-cyclic trend. Such a strong variability suggested a large number of variables which may affect the rotational behavior of similar objects. In this sense, the family of retired Geostationary Operational Environmental Satellite (GOES) 8-12 represents an incredible case study. These five nearly identical defunct satellites mainly differ for the End-of-Life (EoL) configuration in terms of movable surfaces orientation. The objects show large diversity in spin state [9], with different spin periods and regimes, showing both principalaxis and tumbling spin states. This diversity was addressed to the large impact of movable surfaces orientation on the external torque generated by SRP and thermal re-emission, as modelled by the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect [10].
In a more general approach, some solutions have been proposed to serve as simulators to analyze different RSOs. In 2015, The In-Orbit Tumbling Analysis tool ( OTA) [11] was proposed by Hyperschall Technologie Göttingen GmbH and Astronomical Institute of the University of Bern. More recently, an open-source coupled orbit-attitude propagation software, the Debris SPin/Orbit Simulation Environment (D-SPOSE), was presented by McGill University [12]. Both these tools are based on numerical integration of sixdegrees-of-freedom coupled orbit and attitude motion of the object in Earth orbit. This approach is the most common but also the most expensive in terms of computational time and it may be not so accurate for long-term propagations because of accumulation of truncation errors over time. In this work, a different approach is presented, which relies on the tumbling-averaged model proposed in [13].
In Sect. 2, the theoretical framework is presented. Starting from the reference systems used in the work, the common Euler-based model and the above-mentioned averaged model are presented. In Sect. 3 the simulation set-up is described, with a presentation of geometries under analysis and estimates of their moments of inertia and Center of Mass location. In Sect. 4, results of long-term propagations on each of these objects are presented.

Theoretical Framework
This study makes use of different reference frames, as represented in Fig. 1 O completes the right-handed system. In case of circular orbits, as it is in this work, Ẑ O is in the radial direction, pointing towards the Sun. The angular velocity of O with respect to N is = nX O where n is the heliocentric mean motion. Therefore, the rotation matrix can be expressed as ON = R 2 (− ∕2)R 3 (nt) , denoting by R i a rotation about the i-th axis and assuming t = 0 when the radial direction coincides with X .
The angular momentum reference frame H:{X H ,Ŷ H ,Ẑ H } is centered in the satellite with Ẑ H along the satellite's rotational angular momentum vector H. Being and , respectively, azimuth and elevation of Ĥ in the O frame, rotation from O to H can be expressed via the rotation matrix HO = R 2 ( )R 3 ( ). The same convention is applied for the principal axes reference frame P:{X P ,Ŷ P ,Ẑ P } whose axes are aligned with the object's principal axes of inertia (i.e., axes which make the inertia tensor diagonal).

YORP Effect Model
The Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect [14] is one of the applications of the Yarkovsky effect [15]. It expresses the influence of thermal radiation force on asymmetric rotating bodies, provoking torques able to influence spin velocity and spin axis direction. Mainly analyzed for its capability of affecting spin period of small asteroids, its effect on defunct satellites in GEO has been investigated by Albuja et al. [16]. In this work, the YORP effect model derived in [17] is leveraged. It accounts for specular and Lambertian diffuse reflection, absorption, and instantaneous re-emission. The force acting on the ith facet of the satellite is expressed by where P SRP is the Solar Radiation Pressure, i is the facet reflectivity, s i is the fraction of reflectivity that is specular, n î n î n i is the unit vector normal to the facet and n î n î n inî n î n i is a matrix outer product, 3x3 is the 3-by-3 identity matrix, û û u is the satellite-Sun unit vector, A i is the area of the ith facet, and c d i is defined as for Lambertian reflection). The last portion of the expression, max(0,û û u ⋅n î n î n i ) is the illumination function, which ensures that only illuminated facets contribute to the total force.
Once the elementary forces have been estimated, the external torque due to YORP acting on the object is found by where n f is the number of facets used to discretize the object and r i r i r i is the position vector from the facet centroid to the object center of mass.

Full Dynamics Model
The most common approach to propagate rotational dynamics of a rigid body relies on numerical integration of Euler's equation coupled with the kinematic differential equation.
In this work, numerical integration of Eqs. (3) and (4) will be referred to as propagation of "full dynamics" where [I] is the body inertia matrix, is the angular velocity of B with respect to N, M is the external torque found in Eq. (2), and [ 0 , x , y , z ] is the quaternion expressing body's inertial attitude (i.e., the relative attitude of B with respect to N). Some disadvantages can be found in this formulation. First, it is computationally expensive, especially with large spin rates or if a stringent tolerance is required. Moreover, being and quaternion components relatively "fast" variables to be propagated, the truncation error may affect the reliability of the results for the long-term simulations which are presented in this work.

Tumbling-Averaged Model
In this work, the tumbling-averaged model which was derived by Benson et al. [13] is leveraged and is here summarized for clarity. It is based on description of rotational dynamics by means of four state variables.
First, the pole (i.e., rotational angular momentum) direction is identified by its spherical coordinates in the O frame, namely and , as shown in Fig. 1.
Another quantity of interest is the effective spin rate e which scales linearly the components of and is defined as The last state variable is the dynamic moment of inertia I d , which is defined as so that H = I d e and T = 1 2 I d e . The instantaneous value of I d with respect to the object principal moments of inertia defines its rotation regime among the different torque-free rotation modes as described in [18]. In this work, the long-axis convention is assumed, so [I] = diag(I i , I s , I l ) with I l < I i < I s . In the regimes of principal-axis rotation, the body is rotating about the maximum inertia axis for I d = I s , about the intermediate one for I d = I i and about the minimum inertia axis for I d = I l . Instead, the body is tumbling with precessing about the maximum inertia axis (i.e., short-axis mode or SAM) for I i < I d < I s , while it is tumbling with precessing about the minimum inertia axis (i.e., long-axis mode or LAM) for The equations of the model are derived starting from an application of the transport theorem to evaluate the deriva-   3 have to be defined. Two methods are proposed to compute averaged torques, a numerical and a semi-analytical one. They both rely on the assumption that, being the external torque a small perturbation with respect to the satellite torque-free motion, average is computed over the torque-free motion (as expressed by the time evolution of the Euler's angles , , ). Euler analytical solution for torque-free rotation of a rigid body [18] is leveraged in this work and is reported in Appendix 1. According to this solution, both nutation and rotation angles are driven by the scaled time parameter , while an independent formulation is required for the precession angle . Therefore, the average for a generic quantity F can be expressed as being K the complete elliptic integral of the first kind [19] and 4K the period of the Jacobi elliptic functions sn and cn (see Appendix 1).

Numerical Tumbling-Averaged Model
In the numerical average approach, for each combination of the state variables ̄ , ̄ , H , Ī d , the external torque is averaged over a given number of body orientations obtained in a torquefree tumbling motion. It can be easily seen that ̄ and H do not impact the averaged torque, so the average is performed for ̄,Ī d combinations. In Fig. 2, a sequence of 10.000 attitudes is computed over 200 rotation periods in a possible LAM and reported both in ( , , ) and ( , ) phase spaces. As is evident, these orientations are more easily expressed as ( , ) combinations. Moreover, it seems that going further in time, the point will cover the whole ( , ) phase space (i.e., ergodicity); thus, averaging over the torque-free rotation is reduced to an area average over the ( , ) phase space.

Analytical Tumbling-Averaged Model
In the analytical average approach, a fully analytical expression is obtained starting from the YORP effect model of where d d d = r r r × n n n , c si = 2 i s i , c ai = 1 − i s i and all the averaged products can be expressed by elliptic function averages, which are well known in the literature [19]. A complete discussion can be found in [13].

Simulation Set-Up
This work deals with long-term YORP-driven rotational dynamics simulations of different Resident Space Objects (RSOs) in GEO. The objects are assumed to be placed in a circular heliocentric orbit at 1 AU, so P SRP = 4.56 × 10 −6 N∕m 2 in Eq. (1), while their geocentric orbit is neglected. Earth eclipses are also neglected, which is a reasonable assumption for satellites in GEO, since their maximum duration is estimated to be lower than 72 min per day. For all the objects under analysis, a simplified 3D model of their external shape was reconstructed by an external software and imported in MATLAB in the form of .stl file. It allowed to work on a triangular facet model with precise orientation of facet unit normal vectors. The obtained faceted models are shown in Fig. 3.
In Fig. 3a, a faceted model of a GOES 8-12 satellite is shown. Accurate dimensions were found in the literature [20], as well as estimates of its principal moments of inertia and center of mass. It represents an interesting test case under YORP effect especially for its asymmetric shape, composed by the solar sail, the almost cubic bus and two movable surfaces, a wide solar panel and a small trim tab.
A 3U CubeSat with deployed solar panels is displayed in Fig. 3b. No specific object is taken as a reference, due to the variety of CubeSat configurations and to lack of notable examples of this kind of objects in GEO so far.
In Fig. 3c, the rocket body under analysis is shown. In this case, the Titan III-C Transtage served as a reference. Recently, the interest in this RSO has increased because of some documented fragmentation events [21]. Characteristic dimensions and masses were found in the literature or assumed by technical drawings [22], allowing for more accurate estimates of moments of inertia and center of mass. In particular, the actual cylindrical structure is simplified as an octagon to limit the number of facets, as well as the two engines outcoming from the structure. The other visible part is the extremity of the oxidizer tank, while the smaller fuel tank is stored inside and not visible.
The box-wing satellite model is displayed in Fig. 3d. Even though specific dimensions are taken from Optus B [23], results of this analysis can be likely applied to the wide As it is well known, parameters like inertia and center of mass location are fundamental when dealing with rotational dynamics. In this work, a large attention was addressed to their estimation. Only in the case of GOES, principal moments of inertia and center of mass location were found in the literature [13]. For the other objects, once mass and dimensions were assessed, a reasonable asymmetry in the bus due to internal equipment positioning was assumed, while other components were assumed to be perfectly symmetric. At this point, the center of mass was obtained as a weighted sum, while the inertia matrix was assessed by applications of the parallel axis    Table 1. According to Eq. (1) other fundamental input parameters to evaluate YORP effect are surface optical properties, namely the facet total reflectivity and its specular fraction s. Reasonable values for these quantities for the different sections of the objects are taken from [9] and reported in Table 2.

Results
In this section, results of long-term propagations obtained by the model presented in Sect. 2.3 on geometries presented in Fig. 3 are reported.

GOES 8
In this work, simulations on GOES 8 have been performed to validate the correct implementation of the model.
Being extensively studied in [13], the GOES 8 satellite presents an EoL configuration with solar panel orientation sa = 17 • and trim tab orientation tt = −60 • , where sa > 0 for rotations about −Ẑ B and tt > 0 about Ŷ B . The impact of movable surfaces orientation on the inertia matrix was properly considered, leading to I = diag[3432, 3570, 980]kg ⋅ m 2 , which slightly differs from the one reported in Table 1.
In Fig. 4, the evolution of state variables for a 5-year propagation is shown. The assumed initial conditions are e 0 = 0.3 • ∕s , I d = 0.62I s , 0 = 95 • , 0 = 50 • . Results obtained with both the numerically averaged and the analytically averaged model are compared and they almost overlap: thus in the following figures, unless otherwise stated, results obtained by the analytical averaged model will be reported.
Results of this simulation are qualitatively similar to the ones obtained by Benson et al. [13]; thus, it is reasonable to assume that the models are correctly implemented. Particularly, this simulation shows a "tumbling cycle" performed by GOES 8. As the name suggests, this behavior is interesting for its unexpected periodicity. Starting from a tumbling condition, in the first 2 years, the satellite increases its angular velocity about the minimum inertia axis, so e increases, while I d decreases. In the same period, the angular momentum direction is performing a relatively fast precession about the Ẑ O direction (i.e., the Sun direction) while increasing its elevation as visible from ̇ (preferred to for the seek of clarity) and trend. As becomes larger than 90 • , the overall trend changes and the object starts to slow down while moving to a uniform rotation condition (i.e., I d = 1 ), which is reached in less than 4 years. At such a low spin rate, there is not enough pole stiffness, so YORP can drive the pole direction towards the Sun (i.e., decreasing) in a relatively short amount of time, before a new tumbling cycle starts.
The temporal evolution of the angular momentum direction in both O and N reference frames is displayed in Fig. 5. It can be noticed that an apparently chaotic behavior in the inertial frame is instead a quite ordinate precession about the Sun direction in the Orbiting frame. In Fig. 6, for a small part of the simulation shown in Fig. 4, the results of the averaged models were compared with the one obtained propagating the full dynamics as described in Sect. 2.2. It can be seen that the green curve shows qualitatively the same behavior, confirming the reliability of results obtained through the averaged models. Instead, a comparison among computational times underlines the large difference between the 75 min required to propagate the full dynamics for 6 months and the few seconds that are sufficient to propagate the averaged models.

3U CubeSat
The CubeSat model under analysis is representative of a typical End-of-Life configuration with fully deployed solar panels. Even though, nowadays, there are few examples of  such small objects in high orbits, the investigation was led to assess whether the combined effect of reduced exposed area and lower inertia led to a larger or lower impact of YORP on its long-term rotational dynamics.
A Monte Carlo analysis was performed to assess the probability of transition from uniform rotation about the maximum inertia axis to tumbling in less than 5 years. Transition to tumbling is supposed to be completed if the satellite goes below I d = I i (i.e., transition to SAM) at least once in the simulation. 100 random pole directions (i.e., 0 , 0 ) were extracted for different Center of Mass locations and e 0 and results are reported in Table 3.
It can be seen a large impact of the initial spin rate in the first line, where this possibility of having a transition to tumbling in 5 years decreases from 89 to 10%, whereas it becomes negligible if the center of mass is too displaced from the symmetry axis as in the third line, where probability is above 87% even for a large initial angular velocity. Thus, to keep the object in a stable uniform rotation spin state, keeping the center of mass close to the symmetry axis was found to be more important than providing it with a large initial spin rate.
More in det ail, for t he specif ic case of CoM = [3.66, 7.80, 46.6]mm , results of simulations for the only state variable of interest in this case, I d , are shown in Fig. 7, where the dashed line is representative of I d = I i .
A source of asymmetry can be introduced in the system by changing the optical properties of one of the bus facets, which was previously fixed to = 0.6 . This assumption seems reasonable, since one of these faces may be equipped with an antenna or a solar panel. In particular, the effect of a different reflectivity has been assessed in Fig. 8 on a 20-year propagation with initial conditions e 0 = 6 • ∕s , Results show that such an asymmetry is likely to produce a faster and more chaotic behavior in the state variables. In particular, a possible tumbling-to-uniform rotation transition regime may start if a bus facet is provided with a low value of reflectivity such as = 0.27.

Box-Wing Satellite
The model for box-wing satellite is to be intended as representative of the large number of objects composed by an almost cubic bus and wide solar panels. As found in [9], orientation of these movable surfaces ( sp 1 , sp 2 ) is supposed to have a large impact on long-term rotational dynamics. In this analysis, the satellite is assumed to be shut down with both solar panels pointing in the same direction ( sp 1 = sp 2 ) and the effect of is investigated. The impact of solar panel orientation on the object inertia was properly taken into  . 7 CubeSat I d evolution for 100 randomly sampled pole directions with: a e 0 = 6 • ∕ s, b e 0 = 8 • ∕ s, and c e 0 = 9 • ∕s account, leading to a slightly different inertia matrix with respect to the one presented in Table 1. In Fig. 9, a 20-year propagation with initial conditions Results show very different behaviors for different values of . The red line, representative of = 40 • , shows the satellite needing about 10 years to transition from the initial uniform rotation to tumbling. Then, the object spins up about the long axis as shown since the very first years in the case of = 240 • . Instead, for = 120 • , a tumbling cycles-driven dynamics similar to the one described in [13] for GOES 8 is found. Here, the object comes back periodically to the starting uniform rotation conditions for a quite long period of time, which can be further investigated in future works.
A deeper analysis can be performed on the drivers of the averaged state variables evolution, which are the averaged torque components as shown by Eqs. (13)- (16). After numerical average operations, the averaged torque components are available in form of lookup tables M * = f (̄,Ī d ) as it was described in Sect. 2.3.1. A larger dimensionality may be necessary to investigate the effect of another parameter, While the low variability of M Y H with in LAM will require a further investigation in future works, it can be noticed that torque approaches zero for unexpected symmetric configurations, such as ≈ ±40 • and ≈ 180 • ± 40 • . A similar result was found also in [13] with slightly different values. It confirms the idea that for objects of this kind, a preferable shut down epoch may be set to obtain a slower evolution in the post-disposal phase.

Rocket Body
Finally, the analysis on the rocket body model intended to explore the long-term behavior of such a large, quasi-oblate body. As for the CubeSat, the trade-off between area exposed to the Sun and inertia could lead to unexpected impact of YORP on this object. In Fig. 12, a possible 20-year evolution of this object is shown for initial conditions e 0 = 0.15 • ∕s , I d 0 = 0.8I s , 0 = 254 • , 0 = 6 • . A slow trend is observed, with e steadily increasing, while I d and seem to converge to a sort of equilibrium state.
To investigate this behavior, lookup tables of averaged parameters derivatives ̇̄ , ̇H , ̇Ī d as function of and I d were extracted and are reported in Fig. 13, where solid black lines represent zero-value iso-level curves. To remove dependency from , ̇̄ has been averaged over it. As it can be seen, the bottom-left corners of these contours show values close to 0, leading to the slow evolution which was observed in the simulation of Fig. 12. This underlines a possible low impact of YORP effect for this kind of objects, probably because of their intrinsic symmetry.
The evolution of angular momentum direction in the O frame is reported in Fig. 14, showing that pole evolution can be referred to a precession about a direction that is tilted with respect to the Sun direction (represented by the yellow dot) towards the out-of-plane direction. The green and the red dots are respectively the initial and final points of the evolution.

Conclusions
In this work, a previously derived tumbling-averaged YORP rotational dynamics model has been applied to different geometries representative of Resident Space Objects in GEO. Leveraging an averaged model, long-term propagations up to 20 years could be performed in a computational time in the order of seconds, in a much faster approach with respect to numerical integration of full dynamics.
Analyses performed on GOES allowed us to confirm the compliance between the averaged approach and the propagation of full dynamics. The asymmetric shape of the satellite makes it an interesting target under YORP torque, since it  Studies on CubeSat quantified the impact of YORP according to Center of Mass location and initial spin rate. A Monte Carlo analysis showed that a high probability for this object to remain in uniform rotation in the first 5 years of the post-disposal phase can be reached only if the Center of Mass is kept close to the symmetry axis. Also, a possible impact of asymmetry in terms of reflective properties was explored.
The box-wing satellite model showed the large impact of End-of-Life solar array orientation on the long-term dynamics. The analysis confirmed the existence of symmetric solar panels orientations for which torque approaches zero, guaranteeing a slower evolution of the object in the post-disposal phase. Prediction of these favourable conditions can be of help in the definition of a preferable shut down epoch.
Finally, the large inertia and the intrinsic symmetry of the rocket body were found to prevent it from a fast rotational evolution. It suggested a low impact of YORP effect on the rotational dynamics of these objects.
In possible future applications, the model can be extended including other sources of external perturbation, such as internal energy dissipation due to flexible structures or fuel sloshing, while including also gravity gradient and magnetic torques can be useful to apply the model to RSOs in lower orbits. The accuracy in predicting long-term rotational behavior of debris objects other than GOES can be tested if accurate information about inertia and initial conditions are available.

Rigid Body Torque-Free Analytical Solution
This formulation is extracted from [18] as one of the possible analytical descriptions of the torque-free motion of a rigid body. In this work, the long-axis convention is applied, so [I] = diag(I i , I s , I l ) with I l < I i < I s .
First, expressing the angular momentum H H H = [I] ⋅ in the P frame, being the rotation matrix PH = R 3 ( )R 1 ( )R 3 ( ) so and can be obtained easily from Eq. (20), while a more complex formulation is required to obtain

Short-Axis Modes' Solution
For SAM, body frame angular velocity = 1 ; 2 ; 3 can be expressed as where the ± is according to the sign of 2 and sn , cn , dn are the Jacobi elliptic functions with elliptic modulus k and the scaled time parameter is related to time t by Instead, solution for is given by (20)