#  A Population Dynamics Approach to the Distribution of Space Debris in Low
Earth Orbit

John Jurkiewicz Department of Mathematical Sciences, University of Wisconsin -
Milwaukee, Milwaukee, WI 53201, USA; Corresponding author<EMAIL_ADDRESS>Peter Hinow∗

###### Abstract

The presence of debris in Earth’s orbit poses a significant risk to human
activity in outer space. This debris population continues to grow due to
ground launches, loss of external parts from space ships, and uncontrollable
collisions between objects. A computationally feasible continuum model for the
growth of the debris population and its spatial distribution is therefore
critical. Here we propose a diffusion-collision model for the evolution of
debris density in Low-Earth Orbit (LEO) and its dependence on ground-launch
policy. We parametrize this model and test it against data from publicly
available object catalogs to examine timescales for uncontrolled growth.
Finally, we consider sensible launch policies and cleanup strategies and how
they reduce the future risk of collisions with active satellites or space
ships.

> Dedicated to the memory of Professor Ching-Shan Chou.

Keywords: Population dynamics, space debris, Kessler Syndrome, diffusion, Low-
Earth Orbit

##  1 Introduction

Since the beginning of human space exploration in 1957, the regions of outer
space in proximity to Earth have been inhabited by artificial objects. While
the population of near-Earth space in the period immediately following this
time was comprised of active man-made bodies, the abandonment of satellites,
inter-object collisions or unitary explosions have led to an ever-increasing
presence of so-called space debris - inactive, non-controllable entities
evolving according to Keplerian laws. This filling population risks runaway
collision of existing debris objects and poses an increasing danger not only
to safe operation of active spacecraft, but to any object existing in near-
Earth space. Indeed, in 2022 a discarded rocket stage collided with the Moon
causing a noticeable crater (Space.com, 2022a, ). Investigation of the
activity and evolution of this debris environment is therefore an important
direction of research.

The theory of artificial satellites has a rich history beginning in the early
1950s. Traditionally, near-Earth space is separated into three canonical
altitude regions, the Low Earth Orbit (LEO, altitude $\leq$ 2,000 km), the
Medium Earth Orbit (MEO, 2,000 - 36,000 km), and the Geosynchronous Earth
Orbit (GEO, $\geq$ 36,000 km). In each of these three broad regions, the
dominant forces and characteristics of existing objects differ substantially,
and therefore each region requires different scientific treatment (Rossi,,
2005). Of these three the GEO region is the most resistant to analysis since
technological limitations mean that observations are limited to objects above
a size of about 30 cm. Further, debris in the GEO region is not as impacted by
the gravitational effects of the Earth and therefore does not exhibit decaying
orbital patterns. Instead the gravitational effects of sun and moon as well as
the solar radiation are the dominant influences on the behavior of objects in
this region. The orbits of objects in this region are by necessity rather
rather static and predictable, as active satellites must maintain a precise
geostationary position above a point on Earth’s surface. For this reason, and
the low density of cataloged objects in GEO (541 active satellites as of
February 2022 (Johnston,, 2022)) the GEO region density is generally
considered to be stationary.

Figure 1: (Left) Distribution of space debris of various sizes as a function
of altitude. The LEO (left), MEO (center) and GEO (right) peaks are clearly
visible. From (Rossi,, 2005); reproduced with permission from Cambridge
University Press. (Right) Simulated growth of the catalog size of known
objects under various mitigation strategies up to the year 2110. The “business
as usual” (BAU, top funnel) is the most pessimistic scenario. From (Nikolaev
et al.,, 2012); reproduced with permission of the Advanced Maui Optical and
Space Surveillance Technologies Conference (www.amostech.com). Note that this
is the “catalog size”, and it does not include spatial localization.

The region most important to human space activity is LEO, due to the abundance
of human activity. High-velocity collisions ($\approx 10$ km s-1) have
occasionally resulted in the destruction of active satellites such as the
Iridium 33-Kosmos 2251 collision in 2009. Unitary satellite breakup or
destruction is also a threat to safe operations. Examples are the aftermath of
the intentional destruction of the Fengyun 1C weather satellite by a Chinese
anti-satellite missile in 2007, or the destruction of Kosmos-1408 in an
Russian anti-satellite weapon test in November 2021 (Carnegie Endowment for
International Peace,, 2022). Each of these situations results in a net
addition to the debris population that, due to the density of extant object,
increases the probability of future debris interactions. As early as 1978,
Kessler and Cour-Palais, (1978) proposed a cascade scenario in which
collisions create new space debris that further endangers active objects,
ultimately resulting in runaway evolution of the orbital object population.
This is now known as the Kessler syndrome. A multitude of mathematical models
and simulation software exists to predict the behavior of catalogued objects
above a certain minimum size, e.g. the EVOLVE (Reynolds and Eichler,, 1995)
and LEGEND (Liou et al.,, 2004) model families. However, there are estimates
of millions of unobserved objects in the 1 - 10 mm size range that are still
capable of causing significant damage, and their undetected nature makes
prediction difficult.

At the level of tracking individual objects, the gold standard in analysis of
artificial satellites has been the Simplified General Perturbations model
(SGP, currently on SGP4) (Hoots,, 1981). This model, based on the work of
Brouwer, (1959); Lyddane, (1963) assumes that ideal Keplerian orbits of
satellites in an Earth-based reference system are perturbed by certain zonal
harmonics that cause significant deviations from the ideal. Analytical
solutions to the equations of motion for objects in such perturbed orbits
allow observed state vectors of orbital objects to be propagated in time. This
is an effective means of orbit determination for individual orbits over short
timescales, and hence ideal for short-term policy and individual tracking.
Nevertheless, this method suffers from certain shortcomings. This model
requires input data to be provided in the format of NORAD two-line element
(TLE) sets, which is a very particular and obfuscating data format to encode
and decode orbital state vectors. Secondly, this model is a single-particle
model, which requires independent application to each cataloged object in
orbit. Therefore, propagating the entire population of orbital debris
potentially requires tens of thousands of separate applications. Thirdly, over
long timescales the predicted orbit can vary significantly from reality
(Riesing and Kahoy,, 2015). Finally, the model itself does not directly
include the presence of inter-object collisions or single-object explosions,
which, though rare, cannot be neglected.

There have been many attempts to incorporate collisions and explosions into
applications of the SGP4 model, or similar. Nikolaev et al., (2012), for
example, utilizes a modeling approach, computing the trajectory of every known
object of size $\geq$10 cm. In the event of a detected collision, well-known
empirical distributions for number and size of particles created in collisions
are used to update the catalog with new objects, whose trajectories are
subsequently propagated along with the original objects as an enlarged cohort,
to account for population growth. NASA’s benchmark is the LEGEND model (Liou
et al.,, 2004; NASA,, 2022). The LEGEND approach is similar to that of
Nikolaev et al., (2012), but more all-encompassing as the entire space up to
GEO is simulated in this model. Additionally, proprietary propagators are used
in place of the standard SGP4. This has been the main tool used by NASA to
study the near-earth environment since its completion in 2004. While very
accurate, these methods still require propagation of thousands of individual
orbits and sophisticated collision detection. Hence, an analytical population-
dynamics model can possibly incorporate the key features of orbital evolution,
collisions and explosions combined with significant computational savings.

Analytical models at the population level were first introduced by McInnes,
(1993), see also (Zhang et al.,, 2019) for a more recent contribution.
Although space debris is by definition non-biological material, the
mathematical modeling approach to its behavior shares a surprising similarity
to that of spatially structured biological populations that are subject to
birth and death processes, diffusion, as well as inter-species interactions
(Okubo,, 1980). The model of McInnes, (1993) includes a convection term and a
growth term due to binary collisions. In addition, space debris is subject to
the effect of solar wind which is included as a diffusion term. In this work,
we propose a spherically symmetric diffusion model with interaction term for
space debris population evolution. We examine public databases of existing
TLEs to initialize the model and set parameters for the subsequent numerical
simulations by comparisons with TLEs, the literature, and the LEGEND
population predictions. The diffusion parameter varies with altitude, thereby
incorporating the effects of both atmospheric drag, gravitational
perturbations and solar wind. The use of a novel analytical model allows for
significant computational savings and avoidance of numerical instability while
capturing the important mechanisms of population growth and providing good
approximation to debris activity. Moreover, we allow for explicit control at
the level of launch policy and potential cleaning mechanisms. This provides
the framework for testing of different strategies to mitigate the growth of
the population density. Our work therefore has the potential to inform future
space law and policy (Gast,, 2022).

##  2 Introducing a spatial population model for LEO

We begin by stating the underlying assumptions on the structure of orbital
debris. It is a natural assumption supported by empirical observations, as in
Figures 2-4, that the orbital population is spherically symmetric. Thus let
$u(r,t)$ denote the density of orbital objects at time $t$ and a distance of
$r$ from the center of the Earth. Assuming the Earth to be perfectly
spherical, the total population size is obtained by integration,

$U(t)=4\pi\int_{r_{E}+200}^{r_{E}+2000}u(r,t)r^{2}\,\mathrm{d}r,$ (1)

where $r_{E}=6378$ km is the equatorial radius of the Earth.

Figure 2: (Left) A computer generated image of space debris as seen from high
Earth orbit. Clearly visible are the LEO cloud and the GEO ring. From NASA’s
Orbital Debris Program Office; available at orbitaldebris.jsc.nasa.gov/photo-
gallery. (Right)  Histogram snapshot of spatial density of orbital population
in LEO as of March 22nd, 2022. Histogram generated by obtaining TLEs of all
catalogued objects, and extracting altitudes by use of pyorbital. This is
distinct from the density functions used in Equation 2.

Figure 3: Cross-sectional views of the distribution of satellite orbits in LEO
as of March 2nd, 2022. Aside from minor conical areas near the poles, orbital
ellipses cover the volume of near-Earth space essentially uniformly.

Figure 4: Histograms of (Left) orbital inclination and (Right) right ascension
(longitude) of ascending node (the longitude at which the orbit intersects
Earth’s orbital plane) as of March 2nd, 2022. These together show that
spherical symmetry is a reasonable assumption to model the dynamics of orbital
debris.

LEO space is affected by several factors, namely gravitational effects from
the Earth-Moon-Sun system, solar wind, and, in lower regions, atmospheric
friction. To account for these factors on the behavior of the LEO population,
we postulate that the spatial density $u(r,t)$ evolves in time according to a
diffusion equation with source and collision term $\Phi$,

$\frac{\partial}{\partial t}u(r,t)=\frac{1}{r^{2}}\frac{\partial}{\partial
r}\left(D(r)r^{2}\frac{\partial}{\partial r}u(r,t)\right)+\Phi(u,r).$ (2)

The diffusion term in this equation represents the combined effects of
atmospheric drag, a multi-body gravitational system and solar wind. Solar wind
is assumed to cause a diffusive behavior on average, due to the periodic
nature of orbits. In order to capture the altitude dependence of these
perturbations, we select a piecewise radially dependent diffusion parameter

$D(r)=\left\\{\begin{array}[]{cl}\alpha\exp(-\lambda(r-r_{E})),&r<r_{E}+1000\\\
\xi,&r\geq r_{E}+1000\\\ \end{array}\right..$ (3)

Here $\alpha$ and $\xi$ are setting the diffusion rate in the near and the far
field, respectively, and $\lambda$ sets the spatial scale. The choice of these
values is discussed in detail in Section 3. The functional form is inspired by
the presence of atmospheric drag which decreases exponentially with altitude.
Above a certain altitude only the solar wind effects the trajectories. In
accordance with the observations of (Lemaitre and Hubaux,, 2013), we choose a
critical altitude of $1000$ km as atmospheric drag works on a timescale of
103-104 years at this altitude.

The term $\Phi$ represents a “birth and death process” for the debris
population. Debris is created from two sources namely collisions between
present objects and deposition of new objects by ground launch. We set

$\Phi(u,r)=C(u)+\Delta(r,t)$ (4)

where $C(u)$ is the rate of object creation due to binary collisions between
orbital objects, and $\Delta(r,t)$ is the rate of deposition of objects into
orbit. We are at the moment ignoring explosions as a source but leave the
possibility that these be incorporated as stochastic impulses in the future.
In McInnes, (1993), a quadratic collision term of Smoluchowski type is
proposed for binary collisions of space debris subject to certain averaging
assumptions,

$C(u)=\frac{\partial u(r,t)}{\partial
t}|_{collision}=\frac{1}{\sqrt{2}}\beta\gamma\sqrt{\frac{GM}{r+r_{E}}}u^{2}(r,t)$
(5)

where $G$ is the gravitational constant, $M$ is the mass of Earth, $\beta$ is
the average number of new particles created in a binary collision, and
$\gamma$ is the mean cross-sectional area of space debris. We assume that the
deposition rate is separable, that is

$\Delta(r,t)=R(r)T(t).$ (6)

The idea is that there are “natural” locations for deposition of working
objects and that merely the amount varies over time. The functional form is
allowed to vary when we consider different strategies to mitigate debris
growth, but this separable form is the standard for business-as-usual
computations.

The model is completed by a homogeneous Dirichlet boundary condition at the
lowest altitude as objects that enter the denser layers of the atmosphere burn
up (although occasionally large pieces of space debris have reached the
surface of the Earth) and exit the population. At the upper altitude we impose
a homogeneous Neumann boundary condition since there is no debris migrating
from the MEO region. Thus we have

$u(r_{E}+200,t)=0,\quad\frac{\partial}{\partial r}u(r_{E}+2000,t)=0$ (7)

for all $t\geq 0$. Finally, we will discuss in Section 3 the selection of the
initial condition

$u(r,0)=u_{0}(r),\quad r\geq r_{E}+200.$ (8)

##  3 Parametrization of the model

To validate the model and aid in setting certain parameter values, we compare
our model’s predictions to those of NASA’s LEGEND model (Liou et al.,, 2004)
in the case of continued deposition with no mitigation (NASA,, 2022). As a
tool, LEGEND is presented to the public in a limited manner, and the available
predictions are only in terms of the total population. Hence, comparison to
our work can be done only by integrating the density at each time stamp and
comparing the evolution of this total population to the predictions of LEGEND.
We run our simulations with an artificially reduced launch rate to enable
accurate comparisons to LEGEND’s predictions obtained from Liou, (2011).

Free parameters and initial conditions are established using real-world
orbital data in the TLE format obtained from Space-Track (US Space Command,,
2010). As such, the data do not describe Keplerian orbits and are specialized
for use with the simplified perturbations models (Hoots,, 1981). This
necessitates the use of specialized software to decode the TLE data into
Keplerian elements that can be used to compute position vectors. We use the
open-source Python distribution pyorbital to extract the osculating orbital
data (pypi.org/project/pyorbital). The altitude at any given epoch is computed
from the most recent TLE preceding the chosen epoch for each cataloged object
to obtain the cumulative population function. This function is then
differentiated and divided by $4\pi r^{2}$, as $r$ ranges from $r_{e}+200$ to
$r_{e}+2000$, to obtain the initial condition $u_{0}(r)$ for Equation (8).

To select the values of diffusion parameters $\alpha$ and $\xi$ for Equation
(3) we use variability present in the available TLE databases, as this
variability is due to the effects of the same perturbing forces. In the same
manner as Riesing and Kahoy, (2015), we take the two most recent TLEs for
every tracked object currently in LEO; the older TLE for each satellite is
propagated to the epoch of the newer TLE, and the error in altitude
(propagated to observed) is recorded, assuming that the more recent TLE gives
the correct osculating value. We gather these errors into groups based on
altitude and time difference in epoch and compute the standard deviation of
each bin. For each group, we use the formula $\sigma=\sqrt{2Dt}$ to gather
diffusion rates. Finally at each altitude we take the average value of these
diffusion rates to model the “true” diffusivity. Choosing $\xi=10^{-4}$ km2
d-1 in keeping with the observations of Lemaitre and Hubaux, (2013) and the
cataloged data, we fit the above piecewise exponential function to this
scatter plot via nonlinear regression as implemented in Wolfram Mathematica,
giving $\alpha=0.5783$ km2 d-1 and thus forcing $\lambda=0.0086$ km-1. A plot
of the diffusion rate as a function of the altitude is shown in Figure 5.

Figure 5:  Diffusivity function for Equation 3 for individual orbits at
different altitudes.

The parameters $\beta$ and $\gamma$ in the collision term (5) represent the
number of new objects created in a binary collision and the average cross-
sectional area of a piece of space debris, respectively. We set $\beta=2000$,
taking the Iridium-Cosmos collision of 2009 as a representative
(ApolloSat.com,, 2022). In practice, we set $\gamma$ to ensure a good fit in
the $L^{1}$ and $L^{\infty}$-norms of our model to accepted standard
predictions, namely those presented in LEGEND. This agreement requires the
presence of a birth term as described in Section 4. The resultant value of
$\gamma$ is $9.98\cdot 10^{-8}$ km-2, corresponding to a cross-sectional
diameter of approximately $17$ cm. This is a sensible value, as the lower
bound for observable diameter is $5-10$ cm. The parameters used in Equations
(2)-(8) are collected in Table 1.

$\alpha$ | 0.5783 km2 d-1
---|---
$\xi$ |  $10^{-4}$ km2 d-1
$\lambda$ | 0.0086 km-1
$\beta$ | 2000
$\gamma$ |  $9.98\times 10^{-8}$ km2
Table 1: Collected parameters for our model Equations (2)-(8).

Figure 6: (Left) Net addition of objects to the LEO population since 1957,
exhibiting a period of approximately $T=5$ y. (Right) Number of satellites
launched into LEO each year from 2000 to 2021, maintained by the UN Office for
Outer Space Affairs (accessible at
https://www.unoosa.org/oosa/osoindex/search-ng.jspx). Note the increase in
launches in 2019, the year Starlink was first launched (Space.com, 2022b, ).

The space and time components of the deposition function from Equation (6) are
possibly the most difficult objects to estimate for our model. First we begin
with

$R(r)=\sum_{k=1}^{n}a_{k}\exp\left(-\frac{(r-r_{k})^{2}}{\sigma_{k}^{2}}\right),$

where the $r_{k}$ are the typical orbit locations and the $\sigma_{k}$ are the
widths of those regions. This is motivated by inspection of the densities from
2014 through 2022; there appear to be certain “natural” altitudes at which
most objects are deposited. Our choice of values was made to give $R(r)$ a
total mass of 1, and is given in Table 2.

To set the $T(t)$ factor, we use the catalog of total number of objects
launched for the last several years, maintained by the UN Office for Outer
Space Affairs (accessible at https://www.unoosa.org/oosa/osoindex/search-
ng.jspx). We interpolate this list of values via cubic spline interpolation,
giving the variable launch rate for the years from 2000 to 2022. For years
prior to 2022, this interpolation is the time-varying factor $T(t)$.
Additionally, we make the simplifying assumption of uniform deposition over
the year. Hence, if $n$ objects are deposited during a year, the effective
rate is $\frac{n}{365}$ per day. This choice results in an agreement between
the propagated and observed densities from 2014 through 2020, see Section 4.

$k$ | 1 | 2 | 3 | 4
---|---|---|---|---
$a_{k}$ (km-3) | $5.599\cdot 10^{-16}$ | $1.39\cdot 10^{-11}$ | $8.39\cdot 10^{-12}$ | $1.39\cdot 10^{-11}$
$r_{k}$ (km) | 200 | 500 | 700 | 850
$\sigma_{k}$ (km) | 7.07 | 20.09 | 100 | 9.98
Table 2: Peak locations, widths and strengths of the deposition function
$R(r)$.

For propagation of current launch policies into the future, we note the
essentially periodic behavior of net object addition into the LEO population,
see Figure 6 left panel. Therefore, for future prediction, we multiply the
current (2022) launch rate $I_{2022}$ by

$P(t)=k_{1}+k_{2}\sin(2\pi f(t-t_{0}))$

That is, we let the basic 2022 launch rate be perturbed periodically in future
epochs. We take $f=0.2$ y-1 based on Figure 6, left panel. $k_{1}$ and $k_{2}$
are dimensionless free parameters that can be used to emulated different
launch policies. In our analysis, we take $k_{1}=1$ and $k_{2}=5.02\cdot
10^{-4}$. Thus for future propagation our time-varying factor is

$T(t)=I_{2022}P(t).$ (9)

##  4 Simulation and comparison with data

We implement numerical simulations of the model (2)-(8) using the Crank-
Nicolson method as the numerical solver. The spherical Laplacian in Equation 2
is discretized by centered differences with an evenly-spaced sampling of
$\Delta r=2.4$ km. In the time direction $\Delta t$ ranges between $0.05$ d
and $3.65$ d. Due to the small magnitude of the diffusivity, the ratio
$\tfrac{D\Delta t}{(\Delta r)^{2}}\leq\tfrac{1}{2}$ for all of these values,
and thus the method is numerically stable. All code used in this paper can be
found at https://github.com/jurkiew4/Space-Debris-Analytic-Model.

The red curves in Figure 7 show the evolution over the past 8 years of radial
space debris density, as obtained from Space-Track.com (US Space Command,,
2010). Note the steady, slow growth prior to 2022. The presence of the
noticeable peak developing around 450 km at the 2022 epoch is due to the
Russian test of their anti-satellite weapon November 2021, creating upwards of
1500 traceable pieces of debris in this altitude range (Carnegie Endowment for
International Peace,, 2022). While such unitary and exceptional events are
difficult to incorporate into a deterministic model, the model simulations
show that these create differing initial conditions with drastically different
possibilities for catastrophic growth.

Figure 7:  Comparison of observed (red) population and propagation of 2014
(blue) populations under the action of ground launches in 2016, 2018 and 2020
respectively, along with the population growth (lower right panel). We observe
an acceptable agreement between model predictions and the observations in both
the $L^{1}$ and $L^{\infty}$-norms.

Figure 8: Initial densities (upper left) in 2020 (blue) and 2022 (red), along
with densities propagated to 2120 (upper right) in the event of no further
deposition The lower panels show the total populations for each scenario. Note
that the 2022 growth is substantially greater due to the Kosmos-1408 debris
spike at 500 km altitude. Note diffusion effects at the lower altitudes, and
the general growth of the population as diffusion and collision occur.  Figure
9: Comparison of total population forecasts for 2020 population by out model
(blue) to the predictions of LEGEND (red) from 2020 to 2120. The LEGEND data
are reproduced from figures presented in Liou, (2011). The launch rate for our
simulations was intentionally reduced from the 2020 value for comparison to
LEGEND.

Figure 9 compares the aggregate population projections of Equations 2-8 to
those of NASA’s LEGEND model (which only reports the total population in a
publicly) for the next 100 years. We observe a close agreement between the
predictions of the model versus those of LEGEND, deviating by a maximum of
approximately 9%. Note these simulations were prepared with the value of
$T(t)$ artificially reduced to match the lower assumed launch rate at the time
the LEGEND data were generated.

Figure 8 shows the evolution of the population starting from two different
epochs, 2020 and 2022, respectively, in the absence of ground launches. It is
clear that the presence of the debris left from the Russian weapons test in
2021 causes significant hazard in the form of a persistent debris band,
causing greater population growth through collisions.

To analyze the effect of space policy on the debris population, we first
assume that ground policy does not account for the removal of space debris.
That is, $\Delta(r,t)$ is assumed to be always positive and of the separable
form specified in Equation (6). Comparisons of the propagated and observed
populations for 2016, 2018 and 2020 are shown in Figure 7. Figure 10 show the
impact of continued deposition on the growth of the 2022 population.
Catastrophic growth is triggered within the next 50 years according to our
evolutionary model.

Figure 10: (Left) Evolution of the 2022 populations (initial in red,
propagated in blue) to 2062 under continued deposition at 2022 rate. The
growth in this scenario over 40 years is a full order of magnitude greater
than the no-launch scenario after 100 years, underscoring the danger of
unmitigated launching. (Right) Propagation of the 2022 density to 2072 with
$\alpha=0.878$ km2 d-1 and $\xi=10^{-3}$ km2 d-1, taking effect at $1500$ km.
Compare this to Figure 8, where the smaller diffusivity allows a greater
buildup of population at lower altitudes.

To investigate the effects of possible active removal strategies, we now allow
$\Delta(r,t)$ to potentially take on negative values. As a first attempt, we
assume the removal a certain fraction of the population, that is
$\Delta(r,t)=-\eta u(r,t)$. Figure 11 shows the effect of a 5% per year
overall removal on the population over the next 50 years. While a blanket
removal strategy is not a feasible policy, this simulation shows that removal
strategies are directly effective in the mitigation of catastrophic growth.

Figure 11: Propagation of the 2022 population to 2072 under action of active
removal of the form $\Delta(r,t)=-\eta u(r,t)$. The removal amounts to $5\%$
per year. Active removal strategies therefore curb the possibility for
catastrophic growth.

##  5 Discussion and conclusion

Model (2)-(8) is a first step towards an accurate continuum model and is
amenable to efficient numerical implementation. Further, the inclusion of the
$\Delta(r,t)$ term allows for in-depth analysis of the effect of ground launch
policy on future evolution. In particular, we see that continuing on the
current course of launching will cause cascading growth to occur at a sooner
epoch than if the debris population were left to evolve on its own. As only a
few examples have demonstrated, an important part of future research is a
careful determination of the diffusion parameter and, once the technology
becomes widely available, the debris removal rate.

Of particular interest to scientists and policymakers are strategies for
reduction of debris input and the active removal of space debris. Such
strategies are a subject of active research. Klima et al., (2018), for
example, approach the topic from a game theoretic perspective. As debris
removal reduces the risk to all active satellites, every stakeholder may delay
their action leading to a strategic dilemma. Especially due to the cost
incurred in their use, such removal strategies, while physically possible as
of writing, are still very much experimental. Still, despite these
shortcomings, these new technologies represent the best feasible means to
avoid the Kessler syndrome. A responsible stewardship of near-Earth space is
clearly a task for public and private institutions alike that benefit from
human space activity (Gast,, 2022).

The means of these removal strategies bears some mention. Bonnal et al.,
(2020) describe an orbiting laser station whose purpose is to improve
knowledge about the trajectories of debris objects. There is also the
possibility to impart velocity changes to objects using this orbital laser.
Takahashi et al., (2018) propose a concept for space debris removal by a bi-
directional plasma plume ejected from a satellite. One of the beams impacts
the debris object while the other keeps the cleaning satellite in place. This
has been demonstrated in a laboratory experiment. All of these removal
techniques act rather locally, as opposed to the global removal considered in
this paper, but could be incorporated into further refinements of the model as
time-varying velocity fields in addition to the passive removal of the air
drag diffusivity currently present.

There exist a number of directions for future refinement and expansion of the
model. The parametrization of the model may be fine-tuned as technology
increases for detecting orbital objects, or as further data comes to light.
The model simulations are sensitive to parameter values. Figure 10 (right)
shows the effect of changing the parameters $\alpha$ and $\xi$, as well as the
critical altitude in the diffusivity $D(r)$. Indeed, a greater diffusivity can
arrest the growth due to collision events.

Further, Figure 8 demonstrates that, despite the high accuracy of the
propagation model in accounting for collision events, certain events (such as
unitary breakup or the Russian anti-satellite test) can occur with no advance
warning and are not directly related to the population at any time. This
hampers the ability to account for such events in a purely deterministic
partial differential equation model. Therefore, the inclusion of stochastic
impulses into $\Delta(r,t)$ could further enhance the accuracy and predictive
power of the model. While the statistics of such events are at present not
widely understood, this presents a very interesting and fresh area for
investigation.

Figure 11 demonstrate the sensitivity of the debris population to the policy
encoded by $\Delta(r,t)$. Because of the importance of policy to the overall
evolution of the population, choosing an effective policy is a key strategic
decision. Mathematically, we can formalize this choice into choosing a policy
$\Delta$ that minimizes a functional

$I(\Delta)=\int_{t_{0}}^{t_{1}}\int_{r_{E}+200}^{r_{E}+2000}F(r,t,u,u_{r},u_{t},\Delta)r^{2}\,\mathrm{d}r\,\mathrm{d}t,$

where the debris populations $u(r,t)$ is constrained by the model (2) - (8).
By formulating a function $F$ that reflects some effective risk metric, the
optimization problem described above can be used to directly inform policy-
making.

## Author statement

We declare that we have no competing interests.

## Acknowledgments

JJ is supported by a graduate fellowship from the Department of Mathematical
Sciences at the University of Wisconsin - Milwaukee. We thank Drs. Kerri Cahoy
(Massachusetts Institute of Technology), Colin McInnes (University of Glasgow)
and Jer-Chyi Liou (NASA) for sharing literature and valuable comments. We also
thank two unknown readers for constructive remarks.

## References

  * ApolloSat.com, (2022) ApolloSat.com (2022).  U.S. Satellite Destroyed in Space Collision.  Access at https://apollosat.com/space-collision/.
  * Bonnal et al., (2020) Bonnal, C., McKnight, D., Phipps, C., Dupont, C., Missonnier, S., Lequette, L., Merle, M., and Rommelaere, S. (2020).  Just in time collision avoidance - A review.  Acta Astronaut., 170:637–651.
  * Brouwer, (1959) Brouwer, D. (1959).  Solution of the problem of artificial satellite theory without drag.  The Astronomical Journal, 64:378–396.
  * Carnegie Endowment for International Peace, (2022) Carnegie Endowment for International Peace (2022).  The dangerous fallout of russia’s anti-satellite missile test.  Access at https://carnegieendowment.org.
  * Gast, (2022) Gast, R. (2022).  Ein All für alle (Die Zeit, 07-23-2022, in German).
  * Hoots, (1981) Hoots, F. R. (1981).  Reformulation of the Brouwer geopotential theory for improved computational efficiency.  Celestial Mechanics, 24:367–375.
  * Johnston, (2022) Johnston, E. (2022).  List of satellites in geostationary orbit.  Access at https://www.satsig.net/sslist.htm.
  * Kessler and Cour-Palais, (1978) Kessler, D. J. and Cour-Palais, B. G. (1978).  Collision frequency of artificial satellites: The creation of a debris belt.  J. Geophys. Res., 83:2637–2646.
  * Klima et al., (2018) Klima, R., Bloembergen, D., Savani, R., Tuyls, K., Wittig, A., Sapera, A., and Izzo, D. (2018).  Space debris removal: Learning to cooperate and the price of anarchy.  Front. Robot. AI, 5:54.
  * Lemaitre and Hubaux, (2013) Lemaitre, A. and Hubaux, C. (2013).  Space debris long term dynamics.  In Celletti, A., Locatelli, U., Ruggeri, T., and Strickland, E., editors, Mathematical Models and Methods for Planet Earth, volume 6 of INdAM, pages 113–121. Springer, Cham, Switzerland.
  * Liou, (2011) Liou, J.-C. (2011).  An active debris removal parametric study for LEO environment remediation.  Adv. Space Res., 47:1865–1876.
  * Liou et al., (2004) Liou, J.-C., Hall, D. T., Krisko, P. H., and Opiela, J. N. (2004).  LEGEND – a three-dimensional LEO-to-GEO debris evolutionary model.  Adv. Space Res., 34:981–986.
  * Lyddane, (1963) Lyddane, R. H. (1963).  Small eccentricities or inclinations in the brouwer theory of the artificial satellite.  The Astronomical Journal, 68:555–558.
  * McInnes, (1993) McInnes, C. R. (1993).  An analytical model for the catastrophic production of orbital debris.  ESA J., 17:293–305.
  * NASA, (2022) NASA (2022).  LEGEND: 3D/OD Evolutionary Model.  orbitaldebris.jsc.nasa.gov, Accessed 11-20-2022.
  * Nikolaev et al., (2012) Nikolaev, S., Phillion, D., Springer, H. K., deVries, W., Jiang, M., Pertica, A., Henderson, J., Horsley, M., and Olivier, S. (2012).  Brute force modeling of the Kessler syndrome.  In Advanced Maui Optical and Space Surveillance Technologies Conference (AMOS), Maui, HI, September 2012.
  * Okubo, (1980) Okubo, A. (1980).  Diffusion and Ecological Problems: Mathematical Models.  Springer Verlag, Berlin, New York.
  * Reynolds and Eichler, (1995) Reynolds, R. C. and Eichler, P. (1995).  A comparison of debris environment projections using the EVOLVE and CHAIN models.  Adv. Space Res., 16:127–135.
  * Riesing and Kahoy, (2015) Riesing, K. and Kahoy, K. (2015).  Orbit determination from two-line element sets for ISS-deployed CubeSats.  In 29th Annual AIAA/USA Conference on Small Satellites.
  * Rossi, (2005) Rossi, A. (2005).  Population models of space debris.  In Knežević, Z. and Milani, A., editors, Dynamics of Populations of Planetary Systems, volume 197 of Proceedings IAU Colloquium, pages 427–438. International Astronomical Union.
  * (21) Space.com (2022a).  A rogue 3-ton piece rocket debris just collided with the moon.  Access at https://www.space.com/rogue-rocket-stage-hit-moon-today.
  * (22) Space.com (2022b).  Starlink satellites: Everything you need to know about the controversial internet megaconstellation.  Access at https://www.space.com/spacex-starlink-satellites.html.
  * Takahashi et al., (2018) Takahashi, K., Charles, C., Boswell, R. W., and Ando, A. (2018).  Demonstrating a new technology for space debris removal using a bidirectional plasma thruster.  Sci. Rep., 8:14417.
  * US Space Command, (2010) US Space Command (2010).  Space-track.org.  Space-Track.org, Accessed 08-08-2022.
  * Zhang et al., (2019) Zhang, B., Wang, Z., and Zhang, Y. (2019).  Discrete evolution model based on mean spatial density for space debris environment.  Astrophys. Space Sci., 364:70.