Research Article  Open Access
Yongjie Liu, Yu Jiang, Hengnian Li, "Analytical Propagation of Space Debris Density for Collisions near SunSynchronous Orbits", Space: Science & Technology, vol. 2022, Article ID 9825763, 17 pages, 2022. https://doi.org/10.34133/2022/9825763
Analytical Propagation of Space Debris Density for Collisions near SunSynchronous Orbits
Abstract
The increasing frequency of human launches has led to a dramatic increase in the amount of space debris, especially near sunsynchronous orbits. Most of the fragments are small in size, which may make tracking difficult. Therefore, characterizing the distribution, evolution, and collision risk of small debris has long been a difficult issue. This paper is aimed at investigating the orbital evolution and global dispersion behavior of debris clouds near sunsynchronous orbits. Firstly, the NASA breakup model is used to provide an initial distribution of small fragments after collision events. Secondly, the continuity equation is adopted to propagate the density variation analytically. Furthermore, we introduce some statistical quantities and the entropy of debris clouds to model the randomness and band formation. A theorem concerning the equivalence of the band formation and maximal entropy is presented. The accuracy of the band formation time estimation is also discussed. For noncatastrophic collisions at an altitude of 800â€‰km due to a projectile with a mass of 100â€‰g and a collision velocity of 1â€‰km/s, we compare the analytical and numerical results of space debris density. The results show that the maximal peak error is within 0.17, and the mean square error is about 0.25 at 400 days. Additionally, the entropy of right ascension of the ascending node is 8.5% less than that for debris clouds near an orbit with the same altitude and an inclination of 30â€‰deg. This indicates the concentrating behavior for debris clouds near sunsynchronous orbits.
1. Introduction
Over the past few decades, space debris has been growing rapidly due to the increasing launches and experiments. According to the debris information updated by ESA in April 2022, there exist about 30920 debris objects regularly tracked by Space Surveillance Networks. It is also estimated that there exist 36500 space debris objects greater than 10â€‰cm, 1 million debris objects from greater than 1â€‰cm to 10â€‰cm, and 130 million debris objects from greater than 1â€‰mm to 1â€‰cm. Some typical events have greatly contributed to the population of debris in the low Earth orbit (LEO) region. For example, the Cosmos2251 satellite collided with Iridium 33 satellite in 2009. This is the first hypervelocity collision of two intact satellites at an altitude of 790â€‰km. In recent years, India conducted an antisatellite experiment on March 27, 2019. Most recently, a Russian antisatellite test was finished on Nov 15, 2021, and the Cosmos1408 satellite broke up. These events have greatly increased the number of debris and the collision probability for spacecraft in LEO region. In particular, the debris density in sunsynchronous orbit (SSO) region is the highest [1]. Moreover, most of the fragments have characteristic lengths less than 10â€‰cm. This brings tremendous difficulties in tracking and cataloguing of objects. For this reason, investigating the evolution of debris clouds and probabilistic orbital collision hazards near SSO is a fascinating topic.
Many studies about the evolution of space debris clouds have been conducted. Barrows et al. [2] reviewed various methods applicable to the modeling of the dispersion process. For a single fragment, the orbital evolution can be analyzed in a standard way, such as the SGP4 method and numerical method. Semianalytical theory has also been used for orbit prediction and life time estimation of space debris [3]. However, a breakup event can generate thousands of fragments smaller than 10â€‰cm. It is difficult to monitor these smaller fragments. Additionally, predicting the longterm evolution of debris clouds fragment by fragment is really expensive in computation time. A phased approach was applied to analyze the dispersion of debris clouds by McKnight [4] and Jehn [5]. Based on the geometrical shapes of debris clouds, they divided the dispersion process into three typical phases: ellipsoid, toroid, and band. The first phase occurs when the fragments orbit several rounds, and the second phase occurs if the mean anomaly is randomly distributed in the interval . The third phase corresponds to the random distributions of the right ascension of the ascending node (RAAN) , argument of perigee , and mean anomaly . Letizia et al. [6, 7] adopted the continuity equation of hydrodynamics to characterize the variations of the debris density and the collision probability. Their approach not only modelled the orbital evolution of small fragments but also reduced the computational time efficiently. A multiple dimensional extension for the evolution of debris clouds was also created by these researchers [8]. Later, Letizia [9] considered the effects of solar radiation pressure and atmospheric drag on orbital evolution. Frey et al. [10] analyzed the evolution of debris clouds near highly eccentric Earth orbits (HEO). Wittig et al. [11] investigated the longterm density evolution of debris clouds based on semianalytical and differential algebra techniques. The body of these work in [6â€“11] made an extensive contribution to the research on orbital evolution of debris clouds. Recently, a complete description of the satellite breakup events based on the probabilistic analysis method was provided by Frey and Colombo [12]. The reachable domain theory was applied to model early mediumterm evolution and collision probability of debris clouds by Wen and Gurfil [13] and Wen and Qiao [14].
For collisions near SSO, Jiang et al. [15] analyzed the motion of debris under the perturbations including the nonspherical perturbation of the Earth, the third body perturbation, and the atmospheric drag perturbation. They applied the material point method to describe the impact process of a debris object and spacecraft. Efimov et al. [16] investigated the attitude evolution of space debris in SSO combining perturbation methods and a numerical approach. It was shown that the orbital precession and influence of orbital motion on induced eddy currents might make a significant impact on the attitude evolution. The flyby problem for large debris in SSO was discussed by Baranov et al. [17]. Olympio and Frouvelle [18] studied the mission design of a generic active debris removal spacecraft in SSO region based on lowthrust optimization. Zhao et al. [19] proposed a target sequence optimization strategy to remove dozens of debris from thousands of debris running on SSO.
In this work, we use the continuity equation to characterize the orbital evolution of debris clouds near SSO analytically. The harmonic term and atmospheric drag are considered the main perturbations for the evolution of space debris. The applicability of this method in many cases has been discussed in previous profound work in [6â€“9]. The validation of this approach for debris clouds near SSO is performed here. The main contribution of the current paper can be summarized into the following two aspects. (1) The equivalence between the band formation and maximal entropy is shown. (2) The accuracy of the previous band formation time estimation is analyzed for debris clouds near SSO. For debris clouds caused by noncatastrophic collisions near SSO, the entropy of the RAAN varies slowly, and the band formation is also much slower compared with those near lower inclination orbits. In order to characterize the global randomness of debris clouds, we use entropy and some other statistical quantities. Although entropy is an old concept that has been extensively used in many subjects, its application in modeling debris clouds is rare, as far as we can see. An estimate involving entropy and the mean square error is also presented. To estimate the band formation time, previous analytical estimations usually assume that the uniform distribution of apsidal and nodal lines is achieved when the fastest drift fragment reaches the slowest one. Concerning the evolution of debris clouds near SSO, we show that the previous estimation may not be applicable for this case. Significant concentrating behavior for the RAAN reflects the key influence of the SSO precession property.
This paper is organized as follows. Section 2 provides a brief introduction of the breakup model and orbital distribution characteristics of fragments after a noncatastrophic collision. Section 3 introduces some basic statistical quantities and entropy. The analytical propagation method is also briefly reviewed. Band formation is emphasized in this section. In each section, the numerical examples are presented to illustrate the related theory and methods.
2. Modeling the Breakup Event and the Characteristics of Space Debris
The initial condition of the analytical propagation for the debris clouds can be described with the breakup model. More precisely, the distributions of the physical properties (such as the characteristic length, area to mass ratio, ejection velocity) of the space debris are provided by the breakup model. Furthermore, the orbital information (state vectors or orbital elements) can be obtained using the variation of the velocity due to the collision event for the parent orbit. Although there are many different breakup models describing the impact and fragmentations such as SDM [20] and DAMAGE [21], the NASA breakup model [22] is still most extensively used. This model characterizes the fragmentation of spacecraft and rocket bodies after explosions and collisions. In the following discussion, we concentrate on the description of collisions.
2.1. Revisiting the NASA Breakup Model
In 1998, the NASA breakup model was built by analyzing historical collision events and impact tests. It proposed a framework to estimate the characteristic length, areatomass ratio, and ejection velocity of fragments. The areatomass ratio of a fragment is usually closely related to its life period due to the atmospheric drag. In this study, only fragments with sizes between 1â€‰mm and 8â€‰cm are matters of concern. Collisions with these fragments may disable an operational satellite and cause the explosion of an abandoned satellite or rocket body.
It is widely accepted that collision events can be divided into two types: catastrophic and noncatastrophic collisions. For catastrophic collisions, both the larger spacecraft and smaller projectiles are completely fragmented. Noncatastrophic collisions usually lead to the fragmentation of smaller projectiles and cratering of the larger spacecraft. According to the NASA breakup model, the threshold of catastrophic collisions is determined by the impact energytomass ratio. The critical impact energy per target mass is usually considered to be 40 . According to the reference [23], most of the recorded collision events are noncatastrophic. Our main focus is on noncatastrophic collisions.
The energytomass ratio is defined as
The number of space debris objects with characteristic lengths larger than obeys the power law, i.e., where the units of are and is the reference projectile mass with units of kg. An updated expression of by Krisko [24] is
Once the characteristic length is given, the areatomass ratio obeys a lognormal distribution. Specifically, satisfies a normal distribution with expectation and covariance as where . Similarly, the ejection velocity magnitude obeys a lognormal distribution with
As explained before, impacts by small fragments can lead to explosions of spacecraft. The ejection velocity of fragments due to explosions also obeys a lognormal distribution with
Here, we adopt expression (6) as in references [6, 8]. Moreover, the directions of velocity are assumed to be isotropically distributed in the spatial space. To avoid the generation of too high ejection velocities, an upper bound is usually set [6, 25]. For breakup events caused by explosions, a detailed description of the parameters can be seen in reference [22].
2.2. The Distribution of Orbital Elements after Fragmentation near SSO
Sunsynchronous orbits are usually nearly circular low Earth orbits, with altitudes lower than 1500â€‰km. Due to the nonspherical effects of the Earth, the RAAN () and argument of perigee () of an orbit are not fixed, but always change slowly. When only the effect is considered, the average variation rate can be expressed as [26]:
SSO are characterized by the same precession rate of the orbital plane as the Earthâ€™s revolution around the sun, i.e.,
In this work, we always assume that the parent orbit is sunsynchronous. The breakup event can be considered as an impulse at a fixed position on the parent orbit. Since only noncatastrophic collisions are concerned here, it is assumed that the magnitude of the impulse velocity is much lower than that of the parent spacecraft. Based on the Gauss planetary equation, the effect of on the variation of the orbital elements can be expressed as [26]: where is expressed in the satellite frame , , and . This frame is centered at the moving body with the axis aligned with the velocity direction. The axis is aligned with the normal direction of the orbital plane, while the axis is determined by the right hand law. From Equations (7) and (8), one can imagine that the impulse not only leads to direct variations of the orbit elements but also makes effects on the different drift rates of the angle variables.
2.3. Results
We consider a noncatastrophic collision between a projectile and a satellite with a relative collision velocity of 1â€‰km/s. The mass of the projectile is 100â€‰g and the parent orbital elements of SSO can be seen in Table 1.

The NASA breakup model is applied at three different altitudes to investigate the initial distribution of the debris clouds. Due to the same mass of the projectile and collision velocity, there exist little difference in the physical characteristics of the fragments generated in these cases. As stated in reference [6], the total number of fragments with lengths lying between 1â€‰mm and 8â€‰cm is 2397. Figure 1 yields the distribution of the areatomass ratio and the characteristic length. It shows that about 70% of the fragments have lengths of less than 2â€‰mm (this corresponds to in Figure 1). About 81% of fragments own areatomassratios larger than 0.032â€‰m^{2}/kg (this corresponds to in Figure 1).
Figure 2 shows the initial distribution of the orbital elements as well as the perigee and apogee heights . From Figures 2(a), 2(c), and 2(e), it is known that presents a vshaped distribution with a vertex corresponding to the altitude of the parent orbit. This can be explained by Equation (10). Because the true anomaly of each parent orbit is set to be zero in Table 1, Equation (10) illustrates the fact that , where the subscript zero corresponds to the parent orbit. Thus, is approximately a constant, and the slope can be smaller for a higher orbit. It should be noted that takes on a rightanglelike shape. Moreover, some apogee heights may even be negative due to relatively large eccentricities. The corresponding debris may soon reenter the atmosphere and burn up entirely. From Figures 2(b), 2(d), and 2(f), we can conclude that for 95% of fragments generated above, inclinations lie within degree with respect to the inclination of the parent orbit.
For debris clouds near SSO, the precession rate of the RAAN is an important quantity. Figure 3 illustrates that the precession rates of most fragments are concentrated around 0.984â€‰deg/day. This almost equals â€‰deg/day. Due to the collision event, debris may lie in either quasisunsynchronous orbits or nonsunsynchronous orbits. After a simple analysis, we find that RAAN precession rates of more than 60% of the fragments fall within a 5% error with respect to the parent orbits in these cases.
3. Analytical Propagation of Debris Clouds
Before reviewing the analytical method, the density of the debris clouds and some important statistical quantities must be provided.
3.1. Initial Density Distribution and Statistical Properties of Debris Clouds
As discussed in previous section, the distribution of state vectors or orbital elements can be obtained. Therefore, it is necessary to carry out a transform from the set of orbits into density in a continuous form. Thanks to the work of Kessler [27], Izzo [28, 29], and Frey and Colombo [12], the density function can be expressed in terms of either the Cartesian coordinates or the orbital elements, which are equivalent in some sense. To be precise, one can perform a change of variables to transform from one form to the other. The details can be seen in reference [12]. In this work, we mainly use the density function expressed in terms of the orbital elements.
Based on the above simulations, it can be observed that the eccentricities and inclinations of the fragments are concentrated around those of parent orbits with proper assumptions. The perturbation effects on the eccentricities and inclinations of nearly circular low orbits are relatively small. To simplify considerations, we write the density as a function of the orbital elements and the time . Similar to the work in reference [8], we divide the variation interval of the semimajor axis, RAAN, and argument of perigee into 100 equivalent bins. The initial density at a specific element is approximately equal to the average density of the corresponding subinterval. The analytical method in Subsection 3.3 will present the routine to calculate efficiently. Thus, the density with a single element can be written as: where and correspond to the minimum and the maximum of the semimajor axis of the fragments, respectively.
Based on the density expressed in a single element, one can define the order moment of the debris clouds, as expressed in the probability theory. For example, the order moment of the semimajor axis can be formulated as:
Specifically, for the semimajor axis, we also extend this definition (14) to for its specific physical meaning. We denote . Because the Kepler energy for a massless particle in an elliptical orbit with the semimajor axis is proportional to , one may find that can represent the average energy for the whole debris cloud if the effect of the mass distribution is neglected. Another useful quantity is the mean value, which corresponds to the 1order moment. The quantity can be described as the mean radial distance in a statistical sense. With the basic Cauchy inequality, we can obtain an estimate as follows:
The equality is achieved if and only if is concentrated at a constant. This means that the semimajor axes of all fragments are equal. Therefore, the following proposition can be obtained:
Proposition 1. The inequality always holds: . The equality can only be reached when the semimajor axes of all fragments are equal.
Proposition 1 illustrates the fact that the product of the average scaled energy and semimajor axis has a lower bound 1. Furthermore, the lower bound can be reached if and only if the semimajor axes are equal. It is true that concentrating on a specific value for the semimajor axis of debris clouds is an ideal case. This may only rarely occur. However, the difference between and 1 can reflect the concentrating behavior in some sense. If we want to quantify the distribution of an element at a specific value, the 2order moment is also very convenient for describing to what extent the variable is concentrated. It is closely related to the covariance.
3.2. Entropy and Band Formation of Debris Clouds
In this subsection, we introduce the definition of entropy to provide a measure of the randomness and disorder of the whole debris cloud. It is a concept originating from thermodynamics and is now widely used in fields such as statistical mechanics and communication theory. To simplify our notations, we only define the RAAN entropy of the whole debris cloud as an example. The probability density function can be written as a normalization of the number density , i.e.,
Then, the RAAN entropy is defined as:
During the evolution of debris clouds, we may expect that the RAAN entropy is relatively small in a short time since the RAAN is concentrated around that of the parent orbit. As time goes on, the RAAN entropy will become larger. An ideal case occurs when the RAAN is uniformly distributed in the interval . When the RAAN, argument of perigee, and mean anomaly are all uniformly distributed, the band of debris cloud is considered to be formed. Based on the probability theory, we provide some novel results about the inner relationship between the entropy and randomness of a debris cloud.
Theorem 2. (1)The band of a debris cloud is formed if and only if both the RAAN, argument of perigee, and mean anomaly entropy have achieved their maximum .(2)The entropy can be bounded by the mean square error with respect to the uniform distribution from below. In fact, the following inequality holdswhere the function is expressed as
Proof. (1)By the definition of the band formation, we only need to show that for an orbital element RAAN or argument of perigee, the uniform distribution is equivalent to the maximal corresponding entropy. That is, the entropy achieves its maximum at exactly time when for . Note that satisfies , which implies that it is concave. By Jensenâ€™s inequality, the following inequality can be concluded:Thus, it directly follows that where the equality is achieved if and only if is a constant with respect to the RAAN. Because it is a probability density on , it can only be equal to . This means that the RAAN is uniformly distributed. (2)Denoting , the measure . Thus, we can see that . The inequality (18) is equivalent to the following inequality:with The inequality (22) can be found on page 93 of reference [30]. The whole proof is included in the appendix.
It is noteworthy that the function is simply the scaled mean square error between the distribution of RAAN and the uniform distribution. Furthermore, can be written in a cleaner form:
This theorem provides the entropy with bounds from both sides, i.e.,
Both the distance and the entropy can measure the difference in the distribution relative to the uniform distribution. These terms are related by this concise estimate. Of course, the above conclusions also hold for the argument of perigee. Although all of the above discussions are in a continuous form, the discrete form may be presented by writing the integration into a Riemann sum. We present related results in the remark below.
Remark 3. For discrete variable space , the corresponding probability is denoted as , and the entropy is denoted as , . Then, the following inequality holds where either equality is achieved if and only if the variables are uniformly distributed.
Finally, we comment on the estimate of the band formation time. The band formation not only implies the geometrical shape transition of debris clouds but also a distribution variation from high correlation to randomness. Due to the fast dispersion rate of the mean anomaly, only the distributions of the RAAN and the argument of perigee are examined. The conventional method uses the KolmogrovSmirnov test to compare distributions with the uniform distribution [31]. Ashenberg [32] raised an analytical estimate of the band formation time , as shown in formula (27) based on Equation (10). To derive this, the author assumes that the band formation is achieved once the fastest drift fragment reaches the slowest.
According to formula (27), the band formation time for nearly circular SSO can be calculated. The results are presented in Figure 4. In particular, it is known that the band formation time for the collision event with parent orbits in Table 1 is 145 days, 154 days, and 162 days. In the following sections, we discuss the accuracy of this estimate.
3.3. Analytical Propagation through the Continuity Equation
When the characteristic lengths of the space debris are relatively small, the continuity equation in fluid dynamics can be applied to provide an analytical propagation. As is shown in reference [8], during the first two phases, the nonspherical effects of the Earth play a dominant role. We consider the effect of both the oblateness of the Earth and atmospheric drag before the band formation. Once the band is formed, the angle variables of the debris clouds are randomly distributed, and the effect of the term can thus be neglected. Therefore, only the drag effect is considered after the band formation. One can expect that the atmospheric drag may make the band of debris clouds become â€śthinnerâ€ť and lower.
To describe the evolution of debris clouds, either the state vectors or orbital elements can be used as parameters. Their statistical behavior can be obtained by solving the continuity equation and a density transformation. In this subsection, we briefly introduce the method. The detailed theory can be seen in references [6, 8]. If are the parameters and is the timevarying density function, the continuity equation can be written as the following firstorder partial differential equation (PDE) [6, 8]: where is a vector value function, div is the divergence operator for the vector value function, and represents the variation rates of .
It should be noted that for the divergence structure in Equation (28), the Stokes theorem [33] can be used to obtain some conclusions. Denoting as the integral domain in the parameter space and S as the surface enclosing the domain, it can be observed that where is the infinitesimal area element of the surface and is the outer normal of the surface. The right hand side of Equation (29) is the flux of fragments across the surface . Therefore, we conclude the following proposition:
Proposition 4. The variation rates of the total fragment number in a fixed domain of parameter space equals the flux across the surface .
Following the above argument, the total number of fragments is conserved theoretically during the evolution process due to zero right hand side of Equation (28). However, the fragments are usually considered to reenter the atmosphere and be burned when the altitude is less than 50â€‰km. Therefore, the total number of fragments usually decreases during longterm evolution. These are some of the qualitative properties during the evolution of debris clouds.
In order to solve the PDE (28), the characteristic line method is usually applied [34]. Along the characteristic line, the PDE can usually be reduced to ordinary differential equations (ODE). Figure 5 illustrates that at time , the parameter can be connected to at the initial time through the line. The density at time can also be related to the initial density hypersurface.
By the characteristic line method, the solution can be obtained in a closed form [8]: where is the initial density and are determined by the characteristic line equation:
Therefore, the solving of Equation (31) plays a key role in deriving a fully analytical expression.
As in the previous explanation, we take , and the three parameters are the semimajor axis , RAAN , and argument of perigee . Based on Equations (30) and (31), the density can be formulated as [8] where , , , and . The expressions of are obtained by solving the characteristic equation after some simplified approximations [8], i.e., with , is set as 2.2 [35], and the function is defined by
The main simplified assumptions to obtain Equation (33) are as follows: (a) the parent orbits are nearly circular; (b) the eccentricities and inclinations are considered to be the same as the parent orbit and not affected by the perturbation forces during evolution. Therefore, the distribution of the eccentricities and inclinations can be approximately formulated as the Dirac distribution . It seems reasonable to make these assumptions due to our weakcollision consideration and previous simulation results. However, the eccentricities are surely to be varied by the atmospheric drag. This effect is neglected here. For extensions to modeling the evolution of the distribution of the eccentricity, one can refer to the work in [36]. In the following parts, we validate the analytical propagation method for evolution of debris clouds in the SSO region.
3.4. Results
In this subsection, we mainly focus on the collision event with the parent SSO at an altitude of 800â€‰km. The method can also be applied to other cases. There are other two reasons for choosing the altitude 800â€‰km as a numerical example. Firstly, according to the work in reference [23], the LEO spatial density of unclassified cataloged objects and cataloged intact objects (spacecraft and rocket bodies) is the highest around 800â€‰km. This can be seen in Figures 3 and 4 of reference [23]. Therefore, a parent orbit with an altitude of 800â€‰km is the most interesting case. Secondly, we want to compare the results of the SSO with those of a lower inclination orbit (30â€‰deg) in [6, 8]. In the numerical example of their work, the parent orbit is a LEO with an altitude of 800â€‰km.
To view the dispersion process of fragments in spatial space, we denote the parent orbit plane as the  plane with the axis pointing to the breakup point and the axis orthogonal to the axis forming a righthanded coordinate. The evolution process of debris clouds in both the  plane and the radialnormal plane is presented in Figure 6 with the numerical calculation. In the numerical propagation, the main perturbation forces that we include are  harmonic terms of the gravity field, the atmospheric drag, and the solar radiation pressure.
The red circle in the subfigures represents the parent orbit, and the large point marked in black is the breakup point. Other relatively small points represent the fragment projection on the plane. One can observe that during the first few days, the fragments are almost concentrated along the parent orbit. If the Kepler orbit is considered, the breakup point will be a pinch point. However, this is not the case here because both the term and the atmospheric drag are included. Due to the isotropic ejection velocity assumption, Figures 6(b), 6(d), 6(f), and 6(h) illustrate that the distribution of fragments in the normal direction is increasingly scattered. Figures 6(a), 6(c), 6(e), and 6(g) show that the debris clouds become more and more dispersed.
To validate the analytical propagation method, we compare the results with the numerical density propagation fragment by fragment. During analytical propagation, it is noteworthy that the parameter in formula (33) is not a constant but a function of the areatomass ratio. We adopt a strategy similar to that in references [6, 8]: dividing the sets of fragments into certain subsets according to the areatomass ratio. Here, we choose to divide them into 10 bins, and the areatomass ratio is set as the mean value in each bin. Firstly, the density is obtained from expressions (32) and (33). Then, formulas (11), (12), and (13) are applied to derive the densities . The evolution of the debris cloud density in with both analytical and numerical approaches is plotted in Figures 7â€“9.
From Figure 7, we can observe that the peak point always corresponds to 7173â€‰km approximately, namely, the semimajor axis of the parent orbit. Moreover, the heights of peaks are decreasing. Figure 8 shows the density evolution of the RAAN. We can verify that the locations of the peaks are moving almost consistently with formula (9). It is noteworthy that Figures 8(c) and 8(d) present some violent oscillations near the peaks for the analytical methods, which reflect the increased errors at 400 days and 600 days. In particular, there exist obvious differences between two methods at 600 days when the RAAN lies between 150â€‰deg and 300â€‰deg. This is mainly due to the simplified assumptions and approximations in deriving the analytical expression, such as the neglect of the solar radiation pressure. The discontinuity of the blue curve presented in Figure 8(d) also indicates the limitations of the currently employed analytical methods. How to improve the accuracy of analytical approach and avoid this situation is an interesting topic, which is left as our future work. Figure 9 yields that arguments of perigee are distributed more randomly. In general, the analytical results (the blue curve) are in good agreement with the numerical results (the red curve).
Based on the density, many statistical quantities can be calculated with formulas (11), (12), and (13). We here present the evolution of and the mean RAAN for the debris clouds shown in Figure 10. Figure 10(a) validates the idea that the are always larger than 1. This verifies the conclusion in Proposition 1. When the time is between 120 days and 360 days, this quantity is monotonically decreasing. At 360 days, it experiences a sudden increase. Figure 8(b) illustrates that the analytical approach can model the RAAN dispersion behavior very well, except for the time around 360 days. Note that 360 days is approximately equal to the nodal precession period of SSO, i.e., 365.24 days. The error around 360 days can be explained as a â€śresonanceâ€ť effect in some sense.
To analyze the accuracy of the analytical approach, we consider two kinds of indicators: where and correspond to the density with the analytical and numerical approaches, respectively. The indicator reflects the relative differences of the maximal value or the peaks, while actually measures the mean square error by integration. Table 2 shows the relative error for these two approaches. The relative error is usually less than 0.15 when the propagation time is less than 600 days (except for the case of the argument of perigee). This implies that the analytical approach can capture the characteristics of peaks well. However, the mean square root error can even reach 0.37 when the propagation time is 600 days. For other cases listed in Table 2, is usually less than 0.25. Therefore, the analytical approach is applicable when the propagation time is less than 400 days.

The running time advantage is the most evident, which has already been pointed out by many researchers. The programs are built by Matlab2018b in a laptop with a CPU with 4 kernels at 1.80â€‰GHZ. The computation time for the numerical approach is longer than 14 hours if the propagation continues for 400 days. However, the analytical approach takes about 10 minutes. The time taken for the analytical approach is also closely related to the number of bins we choose for the division.
Here, we may compare the results with those in a previous work by Letizia et al. [8]. In that work, the parent orbit was also at 800â€‰km, but the inclination was 30â€‰deg. They showed that the band of debris clouds was approximately formed after 150 days. For debris clouds near SSO at 800â€‰km, Ashenbergâ€™s estimation concluded that the band formation time should be about 154 days, as shown in Figure 4. However, Figures 8 and 9 illustrate the concept that for debris clouds near SSO at 800â€‰km, the band is not formed even after 400 days. The debris cloud is still far from band formation after 154 days. In fact, it is questionable to regard the band formation as the search between the slowest and the fastest fragments. The band formation is a global property. Except for testing the distribution, the entropy may be a promising quantity due to its physical meaning.
In the following, we use entropy to explain the differences in band formation time for these two cases. Figure 11 plots the variation process of the RAAN and the argument of perigee entropy for the debris clouds near SSO at 800â€‰km. To be compatible with the definition of density, we choose to use the discrete form of entropy, as discussed in Remark 3. In formula (26), is set to 100, the number of intervals we specify. Thus, the maximal entropy of the uniform distribution is , corresponding to the red lines in Figure 11. One can observe that the RAAN entropy is increasing when the time is less than 400 days. For the debris clouds in reference [8], the RAAN entropy is larger than that of SSO and reaches a relatively stable value after 150 days. The convergence behavior is evident. In addition, it can be observed that the entropy always varies slowly, but its value is larger than in both cases. In terms of formula (26), the LHS provides a lower bound estimate for the entropy. For the RAAN, the LHS can be verified to be nonpositive for both cases. Therefore, it is not plotted in Figure 11(a). We can obtain that the entropy of RAAN at 400 days in these two cases is 4.366 and 3.995, respectively. In other words, the entropy of RAAN for debris clouds near SSO is about 8.5% less than that of the other case. This illustrates that RAAN for debris clouds near SSO is more concentrated. However, for the argument of perigee, the lower bounds are plotted as dashed curves in Figure 11(b). One can observe similar variation trends for the lower bound and entropy.
4. Conclusion
In this work, we utilize the continuity equation to give a fully analytical description of the evolution of space debris clouds near SSO, including the density in terms of the semimajor axis, RAAN, and argument of perigee. For a projectile with a mass of 100â€‰g and a collision velocity of 1â€‰km/s, we analyze the initial distribution of fragments after noncatastrophic collision events with satellites in SSO at 500â€‰km, 800â€‰km, and 1200â€‰km based on the NASA breakup model. Then, we concentrate on orbital evolution of debris clouds near SSO at 800â€‰km. After comparing with the numerical propagation method, we find that the analytical approach is applicable to model debris clouds within 400 days. The maximal peak error is within 0.17, and the mean square error is about 0.25. Furthermore, the entropy is applied to describe the global randomness of debris clouds. We show that band formation can be defined as the maximal entropy of debris clouds. It can be concluded that the band is not formed within 400 days. Compared with debris clouds near an orbit with an inclination of 30â€‰deg, the entropy of the RAAN is 8.5% less, and the concentrating behavior is more evident here. This reflects some special precession properties of debris clouds near SSO.
The possible aspects of future work include improving the analytical propagation method, estimating the band formation time based on the entropy, and calculating the collision risk between debris clouds and spacecraft before the band formation. Analytical approach considering more perturbation effects can model the evolution of debris clouds more accurately. Based on the debris density, Letizia et al. [7] evaluated the collision probability after the band formation through the Poisson distribution. However, the orbital elements of fragments can be highly correlated during the early phases of debris clouds. The Poisson distribution is not applicable in this case. For debris clouds near SSO, the band formation usually takes more time. Therefore, analyzing the collision probability before the band formation also deserves more attention.
Appendix
Proof of Theorem. (2) for completeness, we here present a clear discussion. When is less than 1 for , the right hand side of (22) is nonpositive. Without loss of generality, we assume that may be larger than 1 and denote . Therefore, we obtain that the left hand side of (22) can be estimated by the function as Moreover, For , we have For , we have Therefore, for , we can conclude that This implies the desired inequality.
Notations
:  Semimajor axis, km 
:  Eccentricity 
:  Inclination, deg 
:  Argument of perigee, deg 
:  Right ascension of the ascending node, deg 
:  Mean anomaly, deg 
:  True anomaly, deg 
h:  Angular moment, kgâ€§m^{2}/s 
:  Semiminor axis, km 
:  Distance from the center of the Earth, km 
:  The mass of target spacecraft, kg 
:  The mass of projectile, kg 
:  Mean value of a normal distribution 
:  Standard deviation of a normal distribution 
:  Relative collision velocity, km/s 
:  Area to mass ratio of debris, m^{2}/kg 
:  The radius of the Earth, 6378.136â€‰km 
:  Gravitational constant, 
:  Mean angular velocity of the Earth around the sun, rad/s 
:  Characteristic lengths of space debris, cm 
:  Argument of latitude, deg 
:  Number density of a variable at time 
:  Entropy of a variable at time 
:  The drag coefficient of fragments 
:  The reference atmospheric density, 
:  The altitude of the collision event, km 
:  The scale altitude in the exponential atmosphere model, 124.64â€‰km [35]. 
Data Availability
The data used to support the findings of this study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this article.
Authorsâ€™ Contributions
Yongjie Liu is responsible for overall coordination, manuscript writing, and data collection. Yu Jiang provided supervisory support and data acquisition and contributed to the writing of the paper. Hengnian Li is the supervisor of the investigation efforts and offered many valuable helps and suggestions. Yu Jiang is the corresponding author.
Acknowledgments
The authors would like to acknowledge the National Natural Science Foundation of China (grant no. U21B2050) for the funding of this work.
References
 L. Anselmo and C. Pardini, â€śRanking upper stages in low earth orbit for active removal,â€ť Acta Astronautica, vol. 122, pp. 19â€“27, 2016. View at: Publisher Site  Google Scholar
 S. P. Barrows, G. G. Swinerd, and R. Crowther, â€śReview of debriscloud modeling techniques,â€ť Journal of Spacecraft and Rockets, vol. 33, no. 4, pp. 550â€“555, 1996. View at: Publisher Site  Google Scholar
 P. Dutt and A. K. Anilkumar, â€śOrbit propagation using semianalytical theory and its applications in space debris field,â€ť Astrophysics & Space Science, vol. 362, no. 2, p. 35, 2017. View at: Publisher Site  Google Scholar
 D. Mcknight, â€śA phased approach to collision hazard analysis,â€ť Advances in Space Research, vol. 10, no. 34, pp. 385â€“388, 1990. View at: Publisher Site  Google Scholar
 R. Jehn, â€śDispersion of debris clouds from inorbit fragmentation events,â€ť in Dresden International Astronautical Federation Congress, Dresden, Federal Republic of Germany, 1991. View at: Google Scholar
 F. Letizia, C. Colombo, and H. G. Lewis, â€śAnalytical model for the propagation of smalldebrisobject clouds after fragmentations,â€ť Journal of Guidance Control, and Dynamics, vol. 38, no. 8, pp. 1478â€“1491, 2015. View at: Publisher Site  Google Scholar
 F. Letizia, C. Colombo, and H. G. Lewis, â€śCollision probability due to space debris clouds through a continuum approach,â€ť Journal of Guidance, Control, and Dynamics, vol. 39, no. 10, pp. 2240â€“2249, 2016. View at: Publisher Site  Google Scholar
 F. Letizia, C. Colombo, and H. G. Lewis, â€śMultidimensional extension of the continuity equation method for debris clouds evolution,â€ť Advances in Space Research, vol. 57, no. 8, pp. 1624â€“1640, 2016. View at: Publisher Site  Google Scholar
 F. Letizia, â€śExtension of the density approach for debris cloud propagation,â€ť Journal of Guidance, Control, and Dynamics, vol. 41, no. 12, pp. 2651â€“2657, 2018. View at: Publisher Site  Google Scholar
 S. Frey, C. Colombo, and S. Lemmens, â€śEvolution of fragmentation cloud in highly eccentric earth orbits through continuum modelling,â€ť in 69th International Astronautical Congress (IAC 2018), Bremen, Germany, 2018. View at: Google Scholar
 A. Wittig, C. Colombo, and R. Armellin, â€śLongterm density evolution through semianalytical and differential algebra techniques,â€ť Celestial Mechanics & Dynamical Astronomy, vol. 128, no. 4, pp. 435â€“452, 2017. View at: Publisher Site  Google Scholar
 S. Frey and C. Colombo, â€śTransformation of satellite breakup distribution for probabilistic orbital collision hazard analysis,â€ť Journal of Guidance, Control, and Dynamics, vol. 44, no. 1, pp. 88â€“105, 2021. View at: Publisher Site  Google Scholar
 C. Wen and P. Gurfil, â€śModeling early mediumterm evolution of debris clouds using the reachable domain method,â€ť Journal of Guidance Control Dynamics, vol. 39, no. 12, pp. 2649â€“2660, 2016. View at: Publisher Site  Google Scholar
 C. Wen and D. Qiao, â€śCalculating collision probability for longterm satellite encounters through the reachable domain method,â€ť Astrodynamics, vol. 6, no. 2, pp. 141â€“159, 2022. View at: Publisher Site  Google Scholar
 Y. Jiang, H. Li, and Y. Yang, â€śEvolution of space debris for spacecraft in the sunsynchronous orbit,â€ť Open Astronomy, vol. 29, no. 1, pp. 265â€“274, 2020. View at: Publisher Site  Google Scholar
 S. Efimov, D. Pritykin, and V. Sidorenko, â€śLongterm attitude dynamics of space debris in sunsynchronous orbits: Cassini cycles and chaotic stabilization,â€ť Celestial Mechanics & Dynamical Astronomy, vol. 130, no. 10, 2018. View at: Publisher Site  Google Scholar
 A. A. Baranov, D. A. Grishko, V. V. Medvedevskikh, and V. V. Lapshin, â€śSolution of the flyby problem for large space debris at sunsynchronous orbits,â€ť Cosmic Research, vol. 54, no. 3, pp. 229â€“236, 2016. View at: Publisher Site  Google Scholar
 J. T. Olympio and N. Frouvelle, â€śSpace debris selection and optimal guidance for removal in the SSO with low thrust propulsion,â€ť Acta Astronautica, vol. 99, pp. 263â€“275, 2014. View at: Publisher Site  Google Scholar
 S. Zhao, J. Zhang, K. Xiang, and R. Qi, â€śTarget sequence optimization for multiple debris rendezvous using low thrust based on characteristics of SSO,â€ť Astrodynamics, vol. 1, no. 1, pp. 85â€“99, 2017. View at: Publisher Site  Google Scholar
 A. Rossi, A. Cordelli, and C. Pardini, â€śModelling the space debris evolution: two new computer codes,â€ť Advances in the Astronautical Sciences, vol. 89, p. 1217, 1996. View at: Google Scholar
 H. G. Lewis, G. Swinerd, N. Williams, and G. Gittins, â€śDAMAGE: a dedicated GEO debris model framework,â€ť in 3rd European Conference on Space Debris, Darmstadt, German, 2001. View at: Google Scholar
 N. L. Johnson, P. H. Krisko, J. C. Liou, and P. D. AnzMeador, â€śNASAâ€™s new breakup model of evolve 4.0,â€ť Advances in Space Research, vol. 28, no. 9, pp. 1377â€“1384, 2001. View at: Publisher Site  Google Scholar
 C. Pardini and L. Anselmo, â€śReview of past onorbit collisions among cataloged objects and examination of the catastrophic fragmentation concept,â€ť Acta Astronautica, vol. 100, pp. 30â€“39, 2014. View at: Publisher Site  Google Scholar
 P. H. Krisko, â€śProper implementation of the 1998 NASA breakup model,â€ť Orbital Debris Quarterly News, vol. 15, no. 4, pp. 1â€“10, 2011. View at: Google Scholar
 A. Rossi, G. Koppenwallner, P. H. Krisko, M. Oswald, and M. Xu, â€śNASA Breakup Model Implementation,â€ť in 24th IADC Meeting, Tsukuba, 2006. View at: Google Scholar
 R. H. Battin, â€śAn introduction to the mathematics and methods of astrodynamics,â€ť in AIAA Education Series, New York, NY, USA, 1999. View at: Publisher Site  Google Scholar
 D. J. Kessler, â€śDerivation of the collision probability between orbiting objects: the lifetimes of Jupiterâ€™s outer moons,â€ť Icarus, vol. 48, no. 1, pp. 39â€“48, 1981. View at: Publisher Site  Google Scholar
 D. Izzo, â€śEffects of orbital parameter uncertainties,â€ť Journal of Guidance, Control, and Dynamics, vol. 28, no. 2, pp. 298â€“305, 2005. View at: Publisher Site  Google Scholar
 D. Izzo, â€śStatistical distribution of Keplerian velocities,â€ť Journal of Guidance, Control, and Dynamics, vol. 29, no. 1, pp. 212â€“215, 2006. View at: Publisher Site  Google Scholar
 L. ErdĹ‘s and H. T. Yau, Dynamical approach to random matrix theory, American Mathematical Soc., 2017.
 J. E. Gentle, Elements of Computational Statist, Springer, 2002.
 J. Ashenberg, â€śFormulas for the phase characteristics in the problem of lowearthorbital debris,â€ť Journal of Spacecraft and Rockets, vol. 31, no. 6, pp. 1044â€“1049, 1994. View at: Publisher Site  Google Scholar
 W. Rudin, Principles of Mathematical Analysis, McGrawHill, 2004.
 Q. Han, A Basic Course in Partial Differential Equations, vol. 120, American Mathematical Society, 2011.
 D. A. Vallado, Fundamentals of Astrodynamics and Applications, Microcosm Press, Hawthorne, USA, 2001.
 S. Frey, C. Colombo, and S. Lemmens, â€śExtension of the KingHele orbit contraction method for accurate, semianalytical propagation of noncircular orbits,â€ť Advances in Space Research, vol. 64, no. 1, pp. 1â€“17, 2019. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2022 Yongjie Liu et al. Exclusive Licensee Beijing Institute of Technology Press. Distributed under a Creative Commons Attribution License (CC BY 4.0).