Greenhouse gases reduce the satellite carrying capacity of low Earth orbit

Main

Most research into greenhouse gas (GHG)-driven climate change in Earth’s atmosphere has historically focused on the troposphere, because changes in this region are consequential for nearly every aspect of life on Earth. Above the troposphere, GHGs continue to have a pronounced effect, though the consequence of climate change in the upper atmosphere is typically perceived as relatively minor. Recently, however, a rapid expansion in the utilization of satellites in low Earth orbit (LEO) for communications, weather forecasting, navigation, defence and more has increased humanity’s reliance on the long-term sustainability of this region. Earth’s thermosphere extends far into LEO, and GHG-driven changes in the structure of the upper atmosphere need to be considered when planning satellite operations moving forward.

Atmospheric constituents are well mixed with nearly constant ratios up to the homopause at approximately 110 km (ref. ¹). Beyond this point, gravitational diffusion decreases the concentration of CO₂ with increasing altitude. GHGs such as CO₂ are particularly good at absorbing and re-emitting infrared (IR) radiation. In the troposphere, greenhouse gases trap and retain some of the IR radiation emitted from Earth’s surface, which leads to an overall warming effect in this region. In the stratosphere, mesosphere and thermosphere, infrared radiation from CO₂ plays an important role in the energy balance by counteracting the solar ultraviolet energy absorption from ozone and molecular oxygen. Conduction pulls heat from the upper thermosphere into the region between 90 and 135 km where it is radiated by CO₂ (refs. ^2,3). Thus, increasing concentrations of CO₂ inevitably leads to cooling in the upper atmosphere. A consequence of cooling is a contraction of the global thermosphere, leading to reductions in mass density at constant altitude over time^4,5.

The thermosphere regularly experiences expansion and contraction in response to solar activity⁶. Solar electromagnetic radiation cycles in intensity roughly every 11 years. The globally averaged thermospheric mass density in LEO can vary by an order of magnitude between solar minimum and maximum⁷. Density fluctuates more rapidly but by similar magnitudes in response to geomagnetic activity enhancements from Joule heating and energetic particle precipitation⁸. This space-weather-driven variation in thermospheric mass density is generally transient in nature but can severely degrade short-term satellite collision avoidance capabilities^9,10. GHG-induced thermosphere contraction is of particular interest in this work because it causes secular changes in density likely to last centuries.

Early efforts in the 1990s first modelled the impact of GHGs in the thermosphere^3,11. More recently, observational evidence derived from satellite drag decay studies have provided support for contraction predictions^4,5,6,12. The Sounding of the Atmosphere using Broadband Emission Radiometry instrument on NASA’s Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics spacecraft collected measurements of temperature and pressure in the mesosphere and lower thermosphere (up to ~105 km altitude) from 2002 to 2021. These data indicate cooling and contraction in the region that will contribute to reduced density at higher altitudes^2,13,14. Additional data from the Halogen Occultation Experiment on NASA’s Upper Atmosphere Research Satellite have provided further evidence of cooling and contraction^15,16. Declining mass density in LEO inevitably leads to a more rapid accumulation of orbital debris^17,18.

The objective of this work is to investigate the impact of thermosphere contraction on sustainable satellite populations in LEO at capacity to prevent runaway debris accumulation. The National Center for Atmospheric Research Whole-Atmosphere Community Climate Model with Thermosphere and Ionosphere Extension (WACCM-X)¹⁹ can incorporate changes in atmospheric composition, making it particularly useful for performing long-term studies simulating the response of the thermosphere to projected GHG emissions^20,21. The United States Naval Research Laboratory’s Mass Spectrometer and Incoherent Scatter Radar (MSIS 2.0)²² model is used to extend the WACCM-X density trends through the helium-dominated region of LEO up to 1,000 km. Beyond 1,000 km, thermospheric mass density is generally so low that it would take many decades for an object to deorbit. More details on density projections are provided in Methods.

Figure 1 shows the anticipated reductions in mass density in the thermosphere spanning 2000–2100. The actual historical record of F_10.7, a standard proxy for solar activity recording the solar radio flux at 10.7 cm, is used to estimate the globally averaged thermospheric mass density for January 2000–October 2023 using MSIS 2.0. After that, F_10.7 and the baseline density are projected according to a least-squares sinusoid that was fit to log₁₀(F_10.7) over the last six solar cycles going back to 1957. This is an idealized projection to show the scale of variability with the solar cycle. Figure 1a–c shows the anticipated reduction in density for three scenarios compared to an atmospheric baseline, which assumes a constant ground-level CO₂ concentration at 369 ppm, the average measured concentration from the year 2000. The three scenarios are a subset of the Shared Socio-economic Pathways (SSP)²³: SSP1–2.6, SSP2–4.5 and SSP5–8.5. Their ground-level CO₂ concentration pathways are shown in Extended Data Fig. 1. SSP1–2.6 sees CO₂ emissions declining to net zero around 2070 followed by net negative emissions. SSP2–4.5 assumes emissions remain around current levels until 2050 before a decline is apparent. SSP5–8.5 is a very high-emissions scenario where CO₂ emissions increase to roughly double the current level. Two apparent trends are clear in Fig. 1a–c. First, there is a slightly larger relative density reduction over time at higher altitudes. Second, the relative density reductions are uniformly more severe during solar minimum than during solar maximum, which leads to oscillations in the density multiplier relative to the baseline case as a function of time. Figure 1d shows the projected mass density at 600 km for the three SSP emissions scenarios and the baseline.

Greenhouse gases reduce the satellite carrying capacity of low Earth orbit — **Fig. 1: Long-term density reductions from thermosphere contraction.**

Space debris accumulation

On-orbit satellite failures, explosions and collisions have contributed to a large population of non-manoeuvrable and often non-trackable debris objects. Meanwhile, decreasing launch costs and maturing satellite technology have created conditions favourable for rapid commercialization across orbital regimes, especially in LEO. Today, a small number of commercial entities operate the large majority of the world’s active satellites as part of proliferated LEO constellations. As the number of anthropogenic objects in orbit increases, the risk of collision grows with it. Kessler and Cour-Palais showed in 1978 that unstable debris growth is possible if debris populations get too large and cause frequent collisions in orbit²⁴. The potential problem of runaway debris growth has subsequently come to be known as Kessler Syndrome. Some altitude shells already have an object density high enough to start experiencing runaway instability, notably around 900 and 1,400 km (refs. ^25,26). Operators are taking note of the change in the safety of the LEO environment as debris accumulates. Most new satellites launched today have some collision avoidance manoeuvring capability and are reliant on this capability to reduce the risk of collision with tracked debris to an acceptable level.

Figure 2a shows the number of tracked objects in the North American Aerospace Defense Command (NORAD) satellite catalogue as a function of altitude and time. The orbit time-averaged altitude of each tracked object was computed from the catalogue of two-line element (TLE) tracks from 1970 to 2023. Altitude bins are 20 km wide. Clear discontinuities may be attributed to the 2008 Chinese Fengyun 1 C anti-satellite test near 900 km altitude (contributing 3,438 pieces of trackable debris)²⁷ and the 2009 collision between the derelict Russian Cosmos 2251 and the active Iridium 33 satellite near 800 km (1,800 pieces of trackable debris)²⁸. Much of this debris will remain in orbit for decades, posing a durable threat to the active satellites in LEO.

**Fig. 2: Orbital debris population growth and life cycle.**

Figure 2b also shows the time-resolved distribution of tracked objects but separated by inclination to help identify individual events. Figure 2c displays the tracked orbital decay path of NORAD object number 4006 using the time-averaged altitude over an orbit from TLEs recorded throughout the object’s orbital lifetime. NORAD 4006 is a debris object generated by the Thor-Ablestar upper-stage rocket body explosion at 900 km altitude in 1961. The background of the plot shows the globally averaged thermospheric mass density as a function of time and altitude, computed using MSIS 2.0 and verified using satellite drag decay studies²⁹. The cyclic variability in density is attributable to the 11-year solar cycle. The debris object decays much faster during solar maximum as compared to solar minimum, when the density is lower. The drag acceleration on a satellite at any time is proportional to the mass density of the atmosphere it travels through. As the thermosphere cools and contracts with increasing GHG concentrations, the lifetime of objects subject to passive atmospheric decay will increase.

As a measure to control the rate of derelict (no longer functional, non-manoeuvrable) satellites introduced into the LEO environment, the United States Federal Communications Commission recently passed a rule that requires operators to deorbit a spacecraft ‘as soon as practicable and no more than five years after the end of their mission’³⁰. ESA has similarly recommended a five-year deorbit timeline as a step towards achieving their zero debris objective. Until recently, a 25-year deorbit guideline has been standard. For non-manoeuvrable satellites, the deorbit timeline is entirely dependent on passive drag decay. Figure 3 shows the altitude that a satellite would need to operate at to deorbit passively within five or 25 years. It assumes a ballistic coefficient of 0.02 m²kg⁻¹, typical for an average satellite. The ballistic coefficient is simply computed from (beta =frac{{C}_mathrm{d}A}{m}), where C_d is the drag coefficient, A is the exposed area in the direction of motion relative to the atmosphere and m is the satellite mass. Most spacecraft have a ballistic coefficient between 0.005 and 0.05 m²kg⁻¹, but the value can change substantially within this range as a function of spacecraft attitude relative to the velocity vector over time³¹. The operating altitude for passive deorbit within a specified time window is observed to oscillate with the solar cycle because the decay rate is much greater during solar maximum. The deorbit start altitude appears low before 2024 because solar cycle 24, which peaked in 2014, had a weak maximum compared to the history of recorded solar maxima (Extended Data Fig. 2). Satellite operators should consider this effect when planning robust disposal strategies, as it may play a role in selecting an operating altitude.

**Fig. 3: Reducing satellite altitudes for regulatory compliance.**

Orbital carrying capacity

Carrying capacity is a concept that has been used in conservation ecology for decades to compute the size of a population that an ecosystem may sustain indefinitely. Populations at capacity are long-term stable, subject to constraints from the environment³². Stable population sizes and optimal distributions are found by rigorously accounting for the sources and sinks of new members.

If no governance action is taken to manage the occupation of Earth’s orbit, the environment will very likely become over-utilized, diminishing the orbital resource and limiting future access. The satellite carrying capacity for LEO may be computed by considering the many potential constraints on space activity. Reasonable constraints include launch limitations³³, spectrum allocation³⁴, limited capabilities in tracking and conjunction assessment³⁵, operational requirements³⁶ and dynamical stability²⁴. Whereas each of these constraints may allow for a different number of satellites in LEO, only the limiting constraints actively play a role in determining the ultimate capacity. It may also be true that different constraints are active within different regions. This work focuses only on the dynamical stability constraint, which determines bounds on activity that prevent widespread runaway debris accumulation.

Any realistic projection of the orbital debris environment should consider that fragmentation events (collisions, explosions, anti-satellite weapon tests and so on) will occur. Because a secular decline in thermospheric mass density will increase the lifetime of debris in orbit, the consequence of a fragmentation event will grow with time. Assuming a constant risk tolerance, the increased consequence of collisions may be compensated for by decreasing their probability. The most reliable way to reduce the probability of collisions is to decrease the number of satellites in the population. Capacity trends should agree with this insight.

Some metrics for measuring the relative environmental impact, or criticality, of tracked objects in orbit have been proposed as a tool for measuring capacity. These include the number-time product³⁷, the criticality of spacecraft index³⁸ and the Environmental Consequences of Orbital Breakup index³⁹, among many others. Whereas these metrics are useful for comparing the relative environmental impact between proposed missions or identifying the best targets for active debris removal, they are poorly suited for defining capacity in a rigorous way. Most existing applications of these metrics to capacity are heuristic and make use of Monte Carlo modelling with trial-and-error initial conditions. In this work, we define an analytical metric that is well suited for understanding how orbital capacity varies over time with changes in the environment.

The instantaneous Kessler capacity (IKC) is defined here as the number of some characteristic satellites that may be placed within an altitude range of interest while avoiding unstable debris growth. This is an instantaneous capacity term because it is computed by taking a ‘snapshot’ of the atmospheric profile at a time of interest and then assuming that this profile is held constant for all time. The dynamics model used to determine this IKC accounts for two ‘species’ with characteristic properties important for computing the debris-generating potential of objects in orbit. The first species is active satellite, which has a mass of 223 kg, a hard-body radius of 0.745 m and is always presumed to be station keeping at constant altitude. The second species is debris object, which has a characteristic mass of 0.64 kg, a hard-body radius of 0.09 m and a ballistic coefficient of 0.0172 m²kg⁻¹. Characteristic properties for satellites and debris were derived from a survey of objects in the tracked catalogue. New satellites have started trending towards a lower hard-body radius and mass, but it is difficult to project this trend forward. An object’s hard-body radius relates to its probability of being involved in a collision. The mass helps determine how many pieces of debris would be generated in a collision, and the ballistic coefficient determines the rate of decay. Rather than tracking individual objects, the source–sink approach used here tracks the distribution of the populations for each species over time. Drag is considered here to be the only force that can cause an object to traverse shell boundaries. Other dynamical effects from Earth’s non-spherical gravity field, solar radiation pressure or luni–solar interactions would largely influence object motion within shells and thus do not influence the population-level source–sink modelling used in this work. Optimal distributions are computed by maximizing the number of active satellites within the region of interest while avoiding unbounded debris growth.

Formulating the dynamics according to this two-species model allows for a straightforward analysis of the stability of the system subject to variations in the input. However, simplified source–sink models have limitations. First, satellite and debris object orbits are modelled as completely random within altitude shells. Actual populations have some inherent structure that is not captured. Second, not all satellites or debris objects may be represented well by the characteristic species in the model. Varying properties such as mass and hard-body radius can influence the risk a satellite poses to the future debris environment. Finally, complex behaviours such as collision avoidance manoeuvres are not considered. For these reasons, the IKC computed in this work should not be considered an absolute limit on the true number of satellites that may occupy a region, but rather a common measuring tool used to observe how changes in the environment impact capacity. More detailed analysis is necessary to determine consistent capacity metrics that can be applied to policy and governance applications.

The total IKC is computed for a given time by considering a stationary atmospheric profile, computing population equilibrium points for each altitude shell within the tracked region, then optimizing weights over the set of altitude shells to maximize the total number of active satellites that may occupy the region of interest. The details for computing the IKC are discussed separately in equations (3)–(11) in Methods. Figure 4 shows the result from computing the IKC for the 200–1,000 km altitude range every three months from 2000 to 2100 for each of the emissions scenarios of interest. For the baseline atmosphere between 200 and 1,000 km altitude, IKC cycles between approximately 34,000,000 characteristic satellites at solar minimum and 72,000,000 satellites at solar maximum, assuming the optimal distribution. However, the region between 200–400 km holds 99.6% of the available capacity during solar minimum and 98.4% of overall capacity during solar maximum. If we instead consider the altitude regime between 400 and 1,000 km (where most satellites operate in practice), the available IKC for the baseline atmosphere drops precipitously to 148,000–1,140,000 characteristic satellites during solar minimum and maximum, respectively.

**Fig. 4: Orbital capacity reductions across potential emissions pathways.**

The inset figure compares the optimal altitude distribution of satellites at solar maximum and solar minimum for the baseline case. A large reduction in local capacity (>90%) happens in the higher altitude shells at solar minimum. It is important to remember that IKC considers only dynamical stability. The inset figure shows that millions of satellites may operate in the lowest altitude shells while avoiding Kessler syndrome instability. Debris from collisions is removed almost immediately at these low altitudes, so collisions can happen at a very high rate and the environment will remain stable. However, few operators would choose to place their satellites in such a region under those circumstances, so the operational requirements probably become the active constraint in the lower shells.

The SSP cases show that in addition to the cyclic variability in IKC that comes with the solar cycle, secular reductions in mass density from thermosphere contraction have a real impact on the overall capacity of this region. By 2100, considering SSP5–8.5 as the worst-case scenario, a 50% reduction in capacity is observed at solar maximum with a 66% reduction at solar minimum. In the more operationally feasible 400–1,000 km range, the worst-case emissions scenario sees a 60% capacity reduction during solar maximum and an 82% loss during solar minimum by 2100. The magnitude of the oscillation between maximum and minimum capacity is reduced over time with increasing GHG emissions as the thermosphere’s response to the solar cycle is dampened by the secular cooling and contraction.

Discussion

Environmental variability in carrying capacity makes communicating and regulating orbital occupation limits in LEO difficult. While the computed capacity available during solar maximum between 200 and 1,000 km is approximately double the capacity at solar minimum, major increases in the occupation of LEO during this period would be misguided. Many satellites have planned operational lifetimes of several decades, whereas the solar-cycle-driven fluctuations in capacity occur over a period of approximately 11 years. If LEO were filled to its solar maximum IKC during solar minimum, rapid runaway debris instability would very likely follow. For these reasons, it is generally wise to plan and regulate according to the worst-case capacity scenario. Considering only the solar cycle, this means regulating activity according to a solar minimum condition. When adding the contribution of additional GHG emissions, the worst-case scenario capacity carries many fewer satellites across broad swaths of LEO by the year 2100.

The IKC capacity metric computed in this work relied on several necessary assumptions, including simplified orbital dynamics, characteristic object species, a lack of coordination among spacecraft and optimal occupation of LEO. It is important to understand these assumptions and the limitations of the models used to compute capacity solutions. Whereas revisiting these assumptions inevitably impacts the magnitude of the instantaneous capacity for an epoch of interest, the secular trend of capacity loss over time with increasing GHG concentrations should be apparent regardless.

This effort has shown that, perhaps unintuitively, climate change and orbital debris accumulation are two pressing issues of inextricable global concern requiring unified action. To regulate and make equitable use of LEO, a prerequisite is understanding how much resource is available to be shared. The satellite carrying capacity of LEO is strongly influenced by variability in the environment, and continued emissions of GHGs will deplete our collective orbital resource. Whereas capacity is a measure of how much orbital resource is available, occupation considers how much of that resource is consumed. Most tools at our disposal for preventing the accumulation of debris work by reducing satellites’ occupation of the orbital environment. Collision avoidance manoeuvres, active debris removal, improved tracking and increased coordination between operators can each be impactful strategies for reducing the debris-generating potential of a satellite in LEO. Unlike these conventional interventions, GHG emissions reductions mitigate capacity loss at the source, because they have a direct impact on satellite drag across all of LEO. Understanding and respecting the influence that the natural environment has on our collective ability to operate in LEO is critical to preventing the exploitation of this regime and protecting it for future generations.

Methods

Projecting thermospheric mass density reductions

WACCM-X is used to model the thermosphere under increasing ground-level carbon dioxide (CO₂) concentrations and varying solar activity levels. The ground-level CO₂ concentrations used in this work are shown in Extended Data Fig. 1. WACCM-X solves for physics, chemistry and dynamics of the atmosphere from ground level through to the upper boundary, which varies in altitude between 400 and 700 km (dependent upon solar and geomagnetic forcing). From some initial condition, the model requires spin up to reach a steady state representative of the true state of the atmosphere. A step-change increase in ground-level CO₂ emissions would require a substantial number of model years (and therefore computing time) to reach the additional mass of CO₂ in the atmosphere required to achieve steady state, effectively becoming a transient run. Instead, to reduce the model time required to reach a spun-up state, a ‘snapshot’ approach⁴⁰ is taken where CO₂ and carbon monoxide (which exists in chemical equilibrium with CO₂) altitude profiles at each latitude–longitude grid point are scaled from the initial state by the ratio between the desired and initial ground-level CO₂ concentrations. WACCM-X is then run for four months to allow a steady state to be reached, then an additional model year is run to allow a global-mean annual-mean to be taken to remove seasonal variation. This is performed for the CO₂ concentrations relating to the year 2000 (369 ppm), and the concentrations from SSP5–8.5 for the years 2005 to 2095 in ten-year steps, chosen to give more snapshots at the lower concentrations for the lower emissions scenarios. These are run at both low and high solar activities assuming F_10.7 is 70 and 200 solar flux units (sfu), respectively. Additional runs at concentrations of 369 and 639 ppm (2065 SSP5–8.5) have been performed at F_10.7 values of 100, 135 and 170 sfu to further understand the impact of varying solar activity on the density reductions.

Helium is a major constituent of the upper thermosphere⁶. Whereas helium is not modelled by WACCM-X, due to its inert nature, it can be added as an uncoupled model⁴¹. MSIS 2.0 is used to achieve this by first calculating atmospheric constituent profiles at each grid point. An adapted version of MSIS 2.0 is used, where the WACCM-X upper thermospheric temperature directly replaces the MSIS exospheric temperature parameter and the Bates temperature profile parameters (including temperature and its gradient at a geopotential height of 122.5 km). This adaptation ensures that thermospheric cooling is accounted for within the calculated, empirical constituent number densities. Then, the MSIS 2.0 helium number density profile is scaled by the ratio of the sum of other constituents’ number density profiles from WACCM-X. This step is summarized by:

$${mathrm{He}}_{mathrm{WACCMX}}=frac{sum _{mathrm{constituents}}{n}_{mathrm{WACCMX}}}{sum _{mathrm{constituents}}{n}_{mathrm{MSIS}2.0}}{mathrm{He}}_{mathrm{MSIS}2.0}$$

(1)

The mass densities of the modelled constituents within WACCM-X are summed, and the same is done with the same constituents from MSIS 2.0. A ratio of these summed constituents is taken and used to scale the helium mass density profiles from MSIS 2.0, which is added to the WACCM-X mass density.

This study also covers altitudes up to 1,000 km, but WACCM-X’s upper boundary ends below this. Bates–Walker profiles⁴² (the backbone of MSIS) are used to calculate each atmospheric constituent’s number density at higher altitudes by assuming diffusive equilibrium, allowing thermospheric mass density to be calculated up to an assigned upper boundary of 1,000 km. This is done via:

$${n}_{i}left(zright)={n}_{i}left({z}_{mathrm{ref}}right)exp left[-frac{{m}_{i}{g}_{mathrm{ref}}}{k{T}_{{infty }}}frac{left(z-{z}_{mathrm{ref}}right)left(R-{z}_{mathrm{ref}}right)}{R+z}right]$$

(2)

Where n_i (z) is the number density of constituent i at an altitude of z, ref is a reference altitude (taken as the top level of WACCM-X), m_i is the atomic mass of the constituent, g_ref is the gravity at the reference altitude, k is the Boltzmann constant, ({T}_{infty }) is the exospheric temperature and R is Earth’s radius.

This maintains the variability in the numerically modelled atmosphere while providing self-consistent number densities for constituents through 1,000 km. Mass densities are then calculated from these. With helium added and mass densities extrapolated to obtain profiles up to 1,000 km at each latitude–longitude grid position, global-mean annual-mean mass densities are taken. These are then divided by the reference value, which is the similar densities at the same F_10.7 value, but a CO₂ concentration of 369 ppm, corresponding to the year 2000. This provides scaling factors (dependent upon ground-level CO₂ concentration and F_10.7) relative to the year 2000, with intermediate values obtained by linear interpolation within this parameter space. The scaling factors have been mapped onto the future CO₂ concentrations provided by the SSPs, such that they then depend upon year and F_10.7. These can then be further mapped onto a future solar activity projection. The historical and projected F_10.7 are shown in Extended Data Fig. 2. By assuming the empirical model (MSIS 2.0) output at some given F_10.7 is representative of the year 2000, applying these scaling factors to the empirical output allows the model to account for the CO₂-induced thermosphere contraction without the additional computation time of running a numerical model such as WACCM-X. More detail on these methods, including comparison of future density reductions against historic trends, is provided by Brown⁴³.

Computing stability-based capacity

A two-species source–sink model governing the populations of satellites (S) and debris (N) is used to investigate bounds on stable capacity. Populations of each species are tracked within a series of discretized altitude bins. For this analysis, satellites station-keep indefinitely and are immediately replaced if destroyed in a collision (that is, S is steady state by construction). Thus, only the population dynamics of the debris need to be tracked. The stability question has two components: (1) what is the maximum number of satellites that can be placed in an altitude bin that allows for long-term stability in the debris population and (2) what is the maximum steady-state debris population at this capacity solution? Question (1) allows for computing bounds on IKC, but question (2) helps to know whether potential operational (risk) constraints are satisfied. There may be a dynamically stable solution that is not operationally feasible and vice versa.

Defining stability based on a source–sink model means that population distributions are tracked instead of individual objects. Because the distribution of S across the altitude bins under study is the design variable of interest, the only population that needs to be tracked is N. In this effort, debris population dynamics within an altitude shell of interest are found using an approach similar to existing source–sink model implementations^44,45,46.

$$,dot{{N}_{k}}={phi }_{{k}_{{NN}}}{eta }_{{NN}}{N}_{k}^{2}+{phi }_{{k}_{{SN}}}{eta }_{{SN}}{S}_{k}{N}_{k}+{phi }_{{k}_{{SS}}}{eta }_{{SS}}{S}_{k}^{2}+frac{{N}_{k+1}}{{tau }_{k+1}}-frac{{N}_{k}}{{tau }_{k}},$$

(3)

where N_k and S_k are the number of debris fragments and satellites in shell k, ϕ is the collision frequency between species based on kinetic gas theory⁴⁷ (so ({phi }_{{k}_{{NN}}}) is the collision frequency for debris-on-debris, ({phi }_{{k}_{{SN}}}) is the collision frequency for satellites-on-debris and so on). ϕ is computed by

$${phi }_{k,{ij}}={rm{pi }}frac{{v}_{{r}_{k}}{left({r}_{i}+{r}_{j}right)}^{2}}{{V}_{k}},$$

(4)

where ({v}_{{r}_{k}}) is the average relative velocity in shell k, r_i and r_j are the hard-body radii of species i and j, respectively, and V_k is the volume of shell k. Similarly, η determines the number of fragments that are created by collisions between species, and is computed using the National Aeronautics and Space Administration (NASA) standard breakup model⁴⁸. τ_k is the average time it takes for debris objects to be removed from shell k by drag. Debris enters shell k at a rate of ({N}_{k+1}/{tau}_{k+1}) from above and leaves the shell at a rate of ({N}_{k}/{tau }_{k}). The time constant τ_k may be approximated for circular orbits by first computing the decay rate in the orbital radius from the middle of the shell:

$${left(frac{mathrm{d}{r}}{mathrm{d}{t}}right)}_{r=h+{r}_mathrm{e}}=frac{{rho }_{h}{v}^{2}{C}_mathrm{d}A}{m}sqrt{frac{{left(h+{r}_mathrm{e}right)}^{3}}{mu ,}}$$

(5)

Here h is the queried orbit altitude, r_e is Earth’s average radius (6,378.15 km), ρ_h is the mass density at altitude h, v is the object velocity, C_d is the drag coefficient, A is the exposed area, m is mass and (mu) is Earth’s gravitational parameter. Assuming constant dr / dt throughout the shell (only valid for small shell widths), the average time in shell k for a debris fragment should be approximately

$${tau }_{k}approx frac{1}{2}frac{mathrm{dh}}{{left(frac{mathrm{d}{r}}{mathrm{d}{t}}right)}_{r=h+{r}_mathrm{e}}}=frac{1}{2}mathrm{dh}{left(frac{{rho }_{h}{v}^{2}{C}_mathrm{d}A}{m}sqrt{frac{{left({h}_{k}+{r}_mathrm{e}right)}^{3}}{mu }}right)}^{-1}$$

(6)

where dh is the width of the shell and h_k is the altitude half-way through the shell.

An important assumption in equation (1) is that the debris flux rate into and out of shells is proportional to the debris populations at the time of interest and the decay rate of that debris. Because it takes some time for debris to move from shell to shell, this assumption breaks down when debris populations grow or shrink rapidly. The assumptions used to derive equation (1) are valid for equilibrium points because debris populations (and thus new debris production rates) will be stable long term. Equilibrium points are much easier to compute as they may be solved in the quadratic form

$$y=a{x}^{2}+{bx}+c$$

(7)

where in this case:

$$begin{array}{l}a={{rm{phi }}}_{{NN}}{{{eta }}}_{{NN}}\ {b}={{rm{phi }}}_{{SN}}{{{eta }}}_{{SN}}{S}_{k}+frac{{N}_{k+1}}{{{rm{tau }}}_{k+1}}-frac{1}{{{rm{tau }}}_{k}}\ c ={{rm{phi }}}_{{SS}}{{{eta }}}_{{SS}}{S}_{k}^{2}.end{array}$$

(8)

Now equilibrium points may be computed analytically based solely on the input constant satellite population size, S. Extended Data Fig. 3a shows the positive equilibrium points for the number of debris objects that are computed as a function of S for a single shell. For any input S, two potential equilibrium points exist. P₁ is a stable node, whereas P₂ is an unstable node. If the initial debris population starts out below the stable node, it will eventually converge on that node. If the debris population size starts above P₁ but below P₂, then it will again converge onto P₁. If, however, the initial population of debris exceeds P₂, then it will experience unbounded growth.

Extended Data Figure 3b–d evolves populations within the single shell of interest forward over many scenarios with varied initial debris counts. They show the result with a steady-state population within the reference shell of 500, 1,500 and 2,500 satellites, respectively. The red and green lines correspond to slices of P₁ and P₂ in Extended Data Fig. 3a. It is important to note that as the steady-state satellite population size increases, P₁ and P₂ approach each other. At some point, P₁ and P₂ converge. In this scenario, any initial debris population will become unstable. Even initial populations below the equivalence point will grow towards P₁ and then immediately grow from that point because a stable node no longer exists. This equivalence point represents an upper bound on the number of satellites that may operate while preventing long-term instability within a shell of interest. Therefore, the capacity is the number of satellites in the shell that are required for no stable debris solutions to exist, and it occurs when P₁ = P₂. This condition is met when

$${b}^{2}-4{ac}=0$$

(9)

$${left({{rm{phi }}}_{{SN}}{{{eta }}}_{{SN}}{S}_{k}+frac{{N}_{k+1}}{{{rm{tau }}}_{k+1}}-frac{1}{{{rm{tau }}}_{k}}right)}^{2}-4{{rm{phi }}}_{{NN}}{{{eta }}}_{{NN}}{{rm{phi }}}_{{SS}}{{{eta }}}_{{SS}}{S}_{k}^{2}=0.$$

(10)

Solving the equation above for S_k determines the number of steady-state satellites for the shell of interest where a stable debris population is possible.

The capacity solution for each shell shows that stability is possible even when the number of debris objects substantially exceeds the number of satellites. The per-shell capacity is greater at lower altitudes because debris is removed very quickly and thus collisions are more acceptable while maintaining stability. At higher altitudes, however, the consequence of a single fragmentation event is much higher, and thus the capacity is reduced.

Optimizing population distributions

Accounting for debris flux from above complicates assessments of orbital capacity across many altitude bins. Filling upper regions to capacity diminishes the capacity of the lower shells as debris enters from above at a higher rate. If the upper shells are filled all the way to their stable capacity, the lower shells can sometimes even have zero capacity for new satellites as they are overwhelmed with the existing debris flux from above. An optimization approach is necessary for determining how much to fill each successive shell to maximize the total capacity across the altitude range of interest. This weighting problem is solved by a successive optimization scheme.

The design variable for the optimization, ({{mathbf{W}}}_{c}), is an array of length n where n is the number of shells and ({{W}}_{{{c}}_{{k}}}) ∈ [0, 1]. Each element of ({{mathbf{W}}}_{c}) represents a weight on the optimal stable capacity ({S}_{k}^{,!* }) for each shell k, which is computed using equation (8). Initial weights are uniformly set to zero. The full optimization is formulated as

$$frac{min }{{{boldsymbol{W}}}_{c}}mathop{sum }limits_{k=1}^{n}-{{W}}_{{{c}}_{{k}}}{S}_{k}^{,!ast }$$

(11)

$$s.t.,{W}_{{c}_{k}}in [0,1],,forall k.$$

Optimizing the entire weight vector is difficult because of its large dimensionality, especially when using many small bins. To reduce the dimensionality of the optimization, each element of ({{mathbf{W}}}_{c}) is optimized sequentially using gradient descent to produce the largest possible population of satellites across bins. Repeating the per-shell optimization process several times eventually leads to convergence for ({{mathbf{W}}}_{c}). Extended Data Fig. 4 shows the optimal weights (({{W}}_{{{c}}_{{k}}})) for each shell after five iterations. The first iteration computes a rough solution that is biased by the initial presumed ({{mathbf{W}}}_{c}). After several iterations, regardless of the initial condition, the weights converge on a single profile. The resulting stable distributions of satellite populations and the corresponding steady-state debris populations are also shown in Extended Data Fig. 4. When computing the instantaneous Kessler capacity (IKC) over time for an altitude region of interest, the optimal distribution of characteristic satellites is computed for each time step by the process outlined above. Then the distribution is integrated to find the overall number of satellites that the region is capable of sustaining when the distribution is optimized. This process is repeated over time to account for changes in the atmospheric state from the solar cycle and GHG emissions.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.