$\displaystyle\gamma_{St}S_{t+1}$ $\displaystyle=0,$ (70)
$\displaystyle\gamma_{Dt}D_{t+1}$ $\displaystyle=0$ (71)
$\displaystyle\lim_{T\to\infty}\left(\frac{1}{1+r}\right)^{T}\lambda_{ST}S_{T+1}$
$\displaystyle=0$ (72)
$\displaystyle\lim_{T\to\infty}-\left(\frac{1}{1+r}\right)^{T}\lambda_{DT}D_{T+1}$
$\displaystyle=0.$ (73)

In what follows we drop time subscripts to reduce notational clutter. Period
$t$ values are shown with no subscript, period $t+1$ values are marked with a
′ after the variable, and period $t-1$ values are marked with a ′ before the
variable e.g. $S_{t-1}\equiv\leavevmode\nobreak\ ^{\prime}S$, $S_{t}\equiv S$,
$S_{t+1}\equiv S^{\prime}$. By (61),

$\displaystyle\lambda_{S}$
$\displaystyle=(1+r)(F+\frac{1}{1+r}m\lambda_{D}-\gamma_{X}+\gamma_{\bar{X}}).$
(74)

In the next period, this becomes

$\displaystyle\lambda_{S}^{\prime}$
$\displaystyle=(1+r)(F+\left(\frac{1}{1+r}\right)m\lambda_{D}^{\prime}-\gamma_{X}^{\prime}+\gamma_{\bar{X}}^{\prime}).$
(75)

By (62) and (63),

$\displaystyle\lambda_{S}$
$\displaystyle=\pi+(1+r)\gamma_{S}+\frac{1}{1+r}\\{\lambda_{S}^{\prime}(1-L(S^{\prime},D^{\prime})-S^{\prime}L_{S}(S^{\prime},D^{\prime}))-\lambda_{D}^{\prime}G_{S}(S^{\prime},D^{\prime})\\}$
(76) $\displaystyle\lambda_{D}$
$\displaystyle=\frac{1}{1+r}\\{\lambda_{D}^{\prime}(1+G_{D}(S^{\prime},D^{\prime})-\delta)+\lambda_{S}^{\prime}S^{\prime}L_{D}(S^{\prime},D^{\prime})\\}-(1+r)\gamma_{D}.$
(77)

Using (75),

$\displaystyle\lambda_{S}$
$\displaystyle=\pi+(1+r)\gamma_{S}-F(L(S^{\prime},D^{\prime})+S^{\prime}L_{S}(S^{\prime},D^{\prime})-1)-\frac{1}{1+r}\lambda_{D}G_{S}(S^{\prime},D^{\prime})$
(78)
$\displaystyle-\frac{1}{1+r}m\lambda_{D}(L(S^{\prime},D^{\prime})+S^{\prime}L_{S}(S^{\prime},D^{\prime})-1)+(L(S^{\prime},D^{\prime})+S^{\prime}L_{S}(S^{\prime},D^{\prime})-1)(\gamma_{X}^{\prime}-\gamma_{\bar{X}}^{\prime})$
(79) $\displaystyle\lambda_{D}$
$\displaystyle=FS^{\prime}L_{D}(S^{\prime},D^{\prime})+\frac{1}{1+r}\lambda_{D}^{\prime}(1+G_{D}(S^{\prime},D^{\prime})-\delta)+\frac{1}{1+r}m\lambda_{D}S^{\prime}L_{D}(S^{\prime},D^{\prime})$
(80)
$\displaystyle-\bigg{(}(1+r)\gamma_{D}+S^{\prime}L_{D}(S^{\prime},D^{\prime})(\gamma_{X}^{\prime}-\gamma_{\bar{X}}^{\prime}\bigg{)}$
(81)

Define

$\displaystyle\alpha_{1}^{\prime}$
$\displaystyle=\pi+(1-L(S^{\prime},D^{\prime})-S^{\prime}L_{S}(S^{\prime},D^{\prime}))F$
(82) $\displaystyle\alpha_{2}^{\prime}$
$\displaystyle=S^{\prime}L_{D}(S^{\prime},D^{\prime})F$ (83)
$\displaystyle\Gamma_{1}^{\prime}$
$\displaystyle=G_{S}(S^{\prime},D^{\prime})-m(1-L(S^{\prime},D^{\prime})-S^{\prime}L_{S}(S^{\prime},D^{\prime}))$
(84) $\displaystyle\Gamma_{2}^{\prime}$
$\displaystyle=1-\delta+G_{D}(S^{\prime},D^{\prime})+mS^{\prime}L_{D}(S^{\prime},D^{\prime})$
(85) $\displaystyle\kappa_{1}^{\prime}$
$\displaystyle=(1+r)\gamma_{S}-(\gamma_{X^{\prime}}-\gamma_{\bar{X^{\prime}}})(1-L(S^{\prime},D^{\prime})-S^{\prime}L_{S}(S^{\prime},D^{\prime}))$
(86) $\displaystyle\kappa_{2}^{\prime}$
$\displaystyle=(1+r)\gamma_{D}+S^{\prime}L_{D}(S^{\prime},D^{\prime})(\gamma_{X^{\prime}}-\gamma_{\bar{X^{\prime}}}),$
(87)

so that

$\displaystyle\lambda_{S}$
$\displaystyle=\alpha_{1}^{\prime}-\frac{1}{1+r}\lambda_{D}^{\prime}\Gamma_{1}^{\prime}+\kappa_{1}^{\prime}$
(88) $\displaystyle\lambda_{D}$
$\displaystyle=\alpha_{2}^{\prime}+\frac{1}{1+r}\lambda_{S}^{\prime}\Gamma_{2}^{\prime}-\kappa_{2}^{\prime}.$
(89)

Then,

$\lambda_{D}^{\prime}=\frac{\lambda_{D}-\alpha_{2}^{\prime}+\kappa_{2}^{\prime}}{\frac{1}{1+r}\Gamma_{2}^{\prime}}.$
(90)

Substitute (74) and (90) in (88) to get the following expression for
$W_{D}(S,D)$

$\frac{\frac{1}{1+r}\\{\Gamma_{1}^{\prime}(\alpha_{2}^{\prime}-\kappa_{2}^{\prime})+\Gamma_{2}^{\prime}(\alpha_{1}^{\prime}+\kappa_{1}^{\prime})\\}+\Gamma_{2}^{\prime}(\gamma_{X}-\gamma_{\bar{X}}-F)}{\frac{1}{1+r}(\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime})}.$
(91)

Iterate 91 to period $t+1$ and substitute into 90 to obtain

$\lambda_{D}^{\prime}=\frac{\frac{1}{1+r}\\{\Gamma_{1}^{\prime\prime}(\alpha_{2}^{\prime\prime}-\kappa_{2}^{\prime\prime})+\Gamma_{2}^{\prime\prime}(\alpha_{1}^{\prime\prime}+\kappa_{1}^{\prime\prime})\\}+\Gamma_{2}^{\prime\prime}(\gamma_{X}^{\prime}-\gamma_{\bar{X}}^{\prime}-F)}{\frac{1}{1+r}(\Gamma_{1}^{\prime\prime}+m\Gamma_{2}^{\prime\prime})}.$
(92)

Use (91) and (92) in (89) to get

$\alpha_{1}^{\prime}=m(\alpha_{2}^{\prime}-\kappa_{2}^{\prime})-\kappa_{1}^{\prime}+\frac{1}{\frac{1}{1+r}}(\gamma_{\bar{X}}-\gamma_{X}+F)+\frac{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}{\Gamma_{1}^{\prime\prime}+m\Gamma_{2}^{\prime\prime}}\bigg{(}\Gamma_{1}^{\prime\prime}\frac{1}{1+r}(\alpha_{2}^{\prime\prime}-\kappa_{2}^{\prime\prime})+\Gamma_{2}^{\prime\prime}(\frac{1}{1+r}(\alpha_{1}^{\prime\prime}+\kappa_{1}^{\prime\prime})-F+\gamma_{X}^{\prime}-\gamma_{\bar{X}}^{\prime})\bigg{)}.$
(93)

Evaluate (93) in the previous time period as:

$\displaystyle\alpha_{1}=m(\alpha_{2}-\kappa_{2})-\kappa_{1}+\frac{1}{\frac{1}{1+r}}(^{\prime}\gamma_{\bar{X}}-^{\prime}\gamma_{X}+F)+\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\bigg{(}\Gamma_{1}^{\prime}\frac{1}{1+r}(\alpha_{2}^{\prime}-\kappa_{2}^{\prime})+\Gamma_{2}^{\prime}(\frac{1}{1+r}(\alpha_{1}^{\prime}+\kappa_{1}^{\prime})-F+\gamma_{X}-\gamma_{\bar{X}})\bigg{)}.$
(94)

Subtract $F(\frac{1}{\frac{1}{1+r}}+L^{\prime})$ from both sides and add
$F(L+SL_{S})$ to both sides to obtain

$\displaystyle\pi-rF-FL(S^{\prime},D^{\prime})$
$\displaystyle=F(L(S,D)+SL_{S}(S,D)-L(S^{\prime},D^{\prime}))+m(\alpha_{2}-\kappa_{2})-\kappa_{1}+\frac{1}{\frac{1}{1+r}}(^{\prime}\gamma_{\bar{X}}-^{\prime}\gamma_{X})+\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}$
$\displaystyle\bigg{(}\Gamma_{1}^{\prime}\frac{1}{1+r}(\alpha_{2}^{\prime}-\kappa_{2}^{\prime})+\Gamma_{2}^{\prime}(\frac{1}{1+r}(\alpha_{1}^{\prime}+\kappa_{1}^{\prime})-F+\gamma_{X}-\gamma_{\bar{X}})\bigg{)}$
(95) $\displaystyle\implies\xi(S^{\prime},D^{\prime})$
$\displaystyle=\underbrace{L_{S}(S,D)SF+\left(L(S,D)-L(S^{\prime},D^{\prime})\right)F}_{\text{\shortstack{Congestion
channel}}}+\underbrace{\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\Gamma_{2}^{\prime}(\frac{1}{1+r}\alpha_{1}^{\prime}-F)}_{\text{\shortstack{Pollution
persistence
channel}}}+\underbrace{\frac{1}{1+r}\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\Gamma_{1}^{\prime}\alpha_{2}^{\prime}}_{\text{\shortstack{Pollution
hazard channel}}}$
$\displaystyle\underbrace{+m\alpha_{2}}_{\text{\shortstack{Pollution\\\
hazard\\\
channel}}}+\underbrace{\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\bigg{(}\Gamma_{2}^{\prime}(\frac{1}{1+r}\kappa_{1}^{\prime}+\gamma_{X}-\gamma_{\bar{X}})-\frac{1}{1+r}\Gamma_{1}^{\prime}\kappa_{2}^{\prime}\bigg{)}-(m\kappa_{2}+\kappa_{1})+\frac{1}{\frac{1}{1+r}}(\gamma_{\bar{{}^{\prime}X}}-\gamma_{{}^{\prime}X})}_{\text{\shortstack{Adjustments
for prior or upcoming corner solutions}}}.$ (96)

Along an interior launch path, the MEC $\xi(S^{\prime},D^{\prime})$ reduces to

$\displaystyle\xi(S^{\prime},D^{\prime})$
$\displaystyle=L_{S}(S,D)SF+\left(L(S,D)-L(S^{\prime},D^{\prime})\right)F+\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\Gamma_{2}^{\prime}(\frac{1}{1+r}\alpha_{1}^{\prime}-F)+\frac{1}{1+r}\frac{\Gamma_{1}+m\Gamma_{2}}{\Gamma_{1}^{\prime}+m\Gamma_{2}^{\prime}}\Gamma_{1}^{\prime}\alpha_{2}^{\prime}+m\alpha_{2},$
(97)

and in an interior steady state the MEC further reduces to

$\displaystyle\xi(S,D)$
$\displaystyle=L_{S}(S,D)SF+\Gamma_{2}(\frac{1}{1+r}\alpha_{1}-F)+\frac{1}{1+r}(\Gamma_{1}+m)\alpha_{2}.$
(98)

###  B.3 The collision probability and new fragment formation functions

In this section we derive the functional forms of the collision probability
and new fragment functions, discuss the physical assumptions they encode, and
describe our process for calibrating the physical model in more detail.

For numerical simulations, we model the probability that objects of type $j$
are struck by objects of type $k$ as

$\displaystyle p_{jk}(k_{t})=1-e^{-\alpha_{jk}k_{t}},$ (99)

where $\alpha_{jk}>0$ is a physical parameter (“intrinsic collision
probability”) reflecting the relative mean sizes, speeds, and inclinations of
the object types (see Letizia (2016) for a derivation of the physical content
of $\alpha_{jk}$). The probability a satellite is destroyed is the sum of the
probabilities it is struck by debris and by other satellites, adjusted for the
probability it is struck by both. For satellite-satellite and satellite-debris
collisions, equation 99 gives us

$\displaystyle L(S,D)$ $\displaystyle=p_{SS}(S)+p_{SD}(D)-p_{SS}(S)p_{SD}(D)$
(100)
$\displaystyle=(1-e^{-\alpha_{SS}S})+(1-e^{-\alpha_{SD}D})-(1-e^{-\alpha_{SS}S})(1-e^{-\alpha_{SD}D})$
$\displaystyle\implies L(S,D)$
$\displaystyle=1-e^{-\alpha_{SS}S-\alpha_{SD}D}.$ (101)

We write the new fragment formation function as

$\displaystyle G(S,D)=F_{SD}p_{SD}(D)+F_{SS}p_{SS}(S)+F_{DD}p_{DD}(D),$ (102)

where $F_{jk}$ is the number of fragments produced in a collision between
objects of type $j$ and $k$. Letting $F_{SS}=\beta_{SS}S$,
$F_{SD}=\beta_{SD}S$, and $F_{DD}=\beta_{DD}D$ where $\beta_{jk}>0$ is a
physical parameter reflecting the physical compositions and masses of the
colliding objects, and using the forms in equation 99, we obtain

$\displaystyle
G(S,D)=\beta_{S}S(1-e^{-\alpha_{SS}S})S+\beta_{SD}(1-e^{-\alpha_{SD}D})S+\beta_{DD}(1-e^{-\alpha_{DD}D})D.$
(103)

The form in equation 101 is convenient as it allows us to solve explicitly for
the open access launch rate and is easy to manipulate. Similar forms have been
used in engineering studies of the orbital debris environment, and are
currently used by the European Space agency in developing indices to study the
long-term evolution of the orbital environment (Letizia, 2016; Letizia et al.,
2017; Letizia, Lemmens, and Krag, 2018).

To derive equation 101, we consider balls (satellites and debris) being placed
into bins (the set of all possible orbital paths within the shell of
interest). The probability of a specific satellite being struck by another
object is then equivalent to the probability that a randomly-placed ball ends
up in a bin containing the specific ball we are focusing on. This is a version
of the “pigeonhole principle”, used in Béal, Deschamps, and Moulin (2020) to
derive a similar form for satellite-satellite collisions.

Suppose we have $b$ equally-sized bins and $n+1$ balls in total, where $b\geq
n+1$. Without loss of generality, we label the ball we are interested in as
$i$. We will first place $i$ into an arbitrary bin, and then drop the
remaining $N$ balls into the $b$ bins with equal probability over bins. The
probability a ball is dropped into a given bin is $\frac{1}{b}$, and the
probability a ball is not dropped into a given bin is then
$\frac{b-1}{b}=1-\frac{1}{b}$. As we drop the remaining $n$ balls, the
probability that none of the balls is dropped in the same bin containing $j$
is

$Pr(\text{no collision with $i$})=\left(1-\frac{1}{b}\right)^{n}$ (104)

Consequently, the probability that any of the $n$ balls are dropped into $i$’s
bin is

$Pr(\text{collision with $i$})=1-\left(1-\frac{1}{b}\right)^{n}.$ (105)

Now suppose we are interested in the probability that members of a collection
of $j$ balls, $1\leq j<b$, end up in a bin with one of the remaining $n+1-j$
balls. The probability that any of the remaining balls end up in a bin with
any of the $j$ balls we are interested in is then

$Pr(\text{collision with $i$})=1-\left(1-\frac{j}{b}\right)^{n+1-j}.$ (106)

As the number of bins and balls grow large ($\lim_{b,n\to\infty}$), we obtain

$Pr(\text{collision with $i$})=1-e^{-j}.$ (107)

Though neither the number of objects in orbits nor the possible positions they
could occupy is infinite, the negative natural exponential form is likely a
reasonable approximation. If we suppose that we have two types of balls $j$
and $k$ of different sizes and bins the size of the smallest type of ball, we
get that the probability a ball of type $k$ is dropped into in a bin with a
ball of type $j$ as

$\displaystyle Pr(\text{$k$--$j$ collision})$
$\displaystyle=1-\left(1-\frac{\alpha_{jk}k}{b}\right)^{n+1-k}$ (108)
$\displaystyle\implies\lim_{b,n\to\infty}Pr(\text{$k$--$j$ collision})$
$\displaystyle=1-e^{-\alpha_{jk}k},$ (109)

which is the form in equation 99, where $\alpha_{jk}$ is a nonnegative
parameter indexing the relative sizes of objects $j$ and $k$. In the orbital
context, $\alpha_{jk}$ reflects not only the sizes of the objects but also
their relative speeds and inclinations. From here we obtain the form of $L$ by
applying standard rules of probability to satellite-satellite and satellite-
debris collisions. Equation 103 follows from the form of $L$.

This “kinetic gas-like” approximation is used extensively in the space debris
modeling literature as a tractable approximation of results from more complex
and computationally-intensive orbital mechanics simulators. It is most
suitable for long-term modeling studies with “large” (relative to the
timescale of orbital interactions) time steps. As described in Letizia (2016),
this approximation is equivalent to modeling collisions as a Poisson process.
The Poisson assumption that the number of events occurring in non-overlapping
time intervals are independent is equivalent to assuming that objects move
randomly throughout the shell volume. This assumption is clearly not true,
leading to our regularization approach described below. The assumption that
the probability of an event is proportional to the length of the interval
implies that fragment clouds are dispersed enough, and contain enough
fragments, to be considered a continuum. Since our model is solved at annual
timesteps while debris clouds evolve at much smaller timescales, this
assumption is reasonable for our purposes.

###  B.4 Modeling debris growth over the next century

As we note in the main text, truly “unbounded” growth is unphysical, as
collisional activity will reduce the fragments to smaller sizes and objects in
LEO will eventually decay due to drag, solar radiation pressure, and other
orbital perturbations. However, we follow the existing engineering literature
on source-sink evolutionary models of the debris environment in allowing
unbounded growth over the next century (Talent, 1992; Lewis et al., 2009;
Lifson et al., 2022). An example using empirical data from a collision and the
size-energy scaling law may help illustrate the underlying reasoning for this
modeling choice.

Consider a collision between two large intact bodies, e.g. an event like the
Iridium-Cosmos collision on February 10, 2009. Iridium 33 was an operational
US communications satellite ($S$) while Cosmos 2251 was defunct Russian
communications satellite ($D$). The table below from Kelso et al. (2009) shows
the relevant size and mass characteristics of the initial objects and
resulting fragments.

Table 1: “Table 1. Pre-Collision Satellite Characteristics.” from Kelso et al. (2009) Satellite | Number of Pieces | Total Volume (m3) | Dry Mass (kg) | Inclination (deg)
---|---|---|---|---
Iridium 33 | 386 | 3.388 | 556 | 86
Cosmos 2251 | 927 | 7.841 | 900 | 74

The average radii for fragments from Iridium and Cosmos were around 12.8 cm
and 12.6 cm, with average masses around 1.44 kg and 0.971 kg. These figures
imply that the tracked fragments larger than 10 cm radius account for most of
the initial body masses.31313110 cm is also the lower detection limit for
sensor systems, raising concerns about censoring. The mass accounting suggests
censoring may not be quantitatively large.

The relation between a uniform sphere’s kinetic energy and mass, given density
$\rho$ and velocity $v$, is

$KE(r)=\frac{1}{2}\underbrace{\rho\left(\frac{4}{3}\pi
r^{3}\right)}_{\text{mass = density$\times$volume}}v^{2}.$ (110)

Suppose a fragment of around 10 cm radius is a uniform aluminum sphere—a
common assumption in debris modeling given the prevalence of aluminum in
satellite construction, e.g. Letizia (2016). Aluminum has a mass of around 2.7
g/cm3, giving a volume of 4188 cm3 and mass of around 11 kg. Typical objects
in low-Earth orbit have velocities on the order of 10 km/s (Lifson et al.,
2022; D’Ambrosio et al., 2023).323232Velocity in orbit is linked with
altitude—accelerating or decelerating along its forward direction raises and
lowers altitude, respectively. Such a fragment will therefore have a kinetic
energy of roughly 550 megajoules , or approximately 131 kg of TNT (energy
equivalent of 1 kg of TNT is 4.184 megajoules). This is in the category of
“hypervelocity” impacts that can shatter the intact object
(esa_hypervelocity1). If the object is like Iridium or Cosmos—not-atypical LEO
satellites—it may produce hundreds of fragments.

Since mass scales cubically with object radius, a reduction in average
fragment size to 1 cm radius reduces the mass to 0.011 kg, producing an impact
energy of 0.55 megajoules—comparable to the force of a hand grenade (European
Space Agency, 2023). Even if it takes tens of collisions with fragments of
1-10 cm radius to overcome shielding on a large intact object, the resulting
tens or hundreds of fragments will ensure net growth. To the extent that these
objects move in debris “fields”—which may occur systematically due to orbital
mechanics factors, particularly when a larger body is struck by a smaller one,
e.g. Oltrogge et al. (2022); Oltrogge, Alfano, and Hall (2022); Pardini and
Anselmo (2023)— their lethal effects at these and even smaller sizes may be
amplified.

Suppose we take 1 cm to be a conservative “lethal size limit”. How long will
it take for collisional activity to reduce a fragment below this limit?
Suppose the average cumulative annual collision probability for an arbitrary
debris fragment is 25%—perhaps a high estimate, but again erring on the side
of caution. That fragment will go roughly 4 years between collisions. If
fragments are reduced to roughly 1 cm radius after only two collisions, it
would take about 8 years for that debris fragment and its children to be
rendered nonlethal. At 1% collision probability, the fragment’s lethal
lifetime is around 200 years.

At 575 km altitude, a large intact object has a residence time (i.e. time
before it falls back to Earth due to drag) on the order of 10 years, and a 10
cm fragment has a residence time on the order of a year, for an upper bound on
lethal lifetime of around 11 years. At 775 km altitude, the residence times
are around 190 years for an intact object and 10 years for a fragment, for an
upper bound on lethal lifetime of around 200 years. During their residence
times the objects slowly drift downwards, entering lower shells. Most
satellites are currently near or above 575 km altitude. Since plausible lethal
lifetimes are on the order of relevant residence times, debris are likely to
spend most of their lives at lethal sizes. Given a sufficiently large amount
of mass at currently-popular altitudes (e.g. 100,000 satellites at 250 kg each
spread over 550-800 km altitude), it seems reasonable to consider potential
growth to “unbounded” levels over the next century.

###  B.5 Open access with a finite horizon

We employ an infinite-horizon modeling approach in the general model. However,
one may reasonably wonder whether our conclusions regarding the open-access
equilibrium are sensitive to this point. In this section we show that a
finite-horizon problem with terminal period $T$ (where it either becomes
prohibitively costly to use the volume or Kessler Syndrome occurs or both)
produces the same equilibrium condition.

Suppose there exists a final period, $T$, such that the potential launchers
will all exit the market. We are agnostic as to why this may be the case,
except to note that if such a period exists, it must be that there are no
profits to be gained from launching after that period. In the final period,
the launcher’s value becomes

$V_{iT}(S_{T},D_{T},X_{T})=\max_{x_{iT}\in\\{0,1\\}}\\{(1-x_{iT})\frac{1}{1+r}V_{iT+1}(S_{T+1},D_{T+1},0)+x_{iT}\left[\frac{1}{1+r}Q(S_{T+1},D_{T+1})-F\right]\\}.$
(111)

There are two possible cases here for the value of launching in the final
period, $\frac{1}{1+r}Q(S_{T+1},D_{T+1})-F$:

  1. 1.

$\frac{1}{1+r}Q(S_{T+1},D_{T+1})-F=0$. In this case the potential launchers
are indifferent between launching in the final period or not launching. By
backwards induction the equilibrium path up to period $T$ will match the one
derived in the general model in the main text, with equation (15) being the
equilibrium condition.

  2. 2.

$\frac{1}{1+r}Q(S_{T+1},D_{T+1})-F<0$. In this case, firms would prefer not to
launch. Optimization by individual launchers therefore implies
$V_{iT}(S_{T},D_{T},X_{T})=0$. This matches equation (14), which yields (15)
after some algebra. So again by backwards induction the equilibrium path up to
period $T$ will match the one derived in the general model in the main text.

Indeed, it is possible to go one step further: the existence of such a
terminal period (where $X_{t}=0\leavevmode\nobreak\ \forall t\geq T$) is
possible if and only if Kessler Syndrome occurs along the equilibrium path.

######  Proposition 4.

A terminal period $T$ where $\hat{X}_{t}=0\leavevmode\nobreak\ \forall t\geq
T$ can exist for an open-access equilibrium path $\\{\hat{X}_{t}\\}_{t}$ if
and only if Kessler Syndrome occurs (i.e. $\lim_{t\to\infty}D_{t}=\infty$)
along the open-access equilibrium path.

###### Proof.

The proposition asserts that

$X_{t}=0\leavevmode\nobreak\ \forall t\geq T\iff\lim_{t\to\infty}D_{t}=\infty$
(112)

We first show the $\impliedby$ direction, then the $\implies$ direction.

The “only if” direction, $X_{t}=0\leavevmode\nobreak\ \forall t\geq
T\impliedby\lim_{t\to\infty}D_{t}=\infty$: If $\lim_{t\to\infty}D_{t}=\infty$,
then there is some period $\bar{t}$ such that $D_{t}>D_{\bar{t}}$ for all
$t>\bar{t}$. From the law of motion for $D$ and our assumption that
$\lim_{t\to\infty}D_{t}=\infty$, we can see that $D_{t}$ must be monotonically
increasing after $\bar{t}$. So there must exist a period $T\geq\bar{t}$ such
that $L(S_{t},D_{t})F>\pi-rF$ for all $t\geq T$, i.e. where it becomes
unprofitable to launch one more satellite at that or any future period. Thus,
$X_{t}=0\leavevmode\nobreak\ \forall t\geq T$. This completes the $\impliedby$
direction.

The “if” direction, $X_{t}=0\leavevmode\nobreak\ \forall t\geq
T\implies\lim_{t\to\infty}D_{t}=\infty$: If $X_{t}=0$ for all $t\geq T$ then
it must be the case that $\frac{1}{1+r}Q(S_{t+1},D_{t+1})-F<0$ for all $t\geq
T$, else some firm would find it profitable to launch. Note that it must be
unprofitable to launch at $t$ given that there are no launches occurring at
$t$.

To be explicit in the next steps, we write the satellite and debris stocks
with the previous-period aggregate launch rate $X_{t}$ shown explicitly as an
argument, i.e. writing $S_{t+1}(X_{t})$ and $D_{t+1}(X_{t})$. Along a path
$\\{S_{t}(0)\\}_{t>T}^{\infty}$, clearly $S_{t+1}(0)\leq S_{t}(0)$. Now,
$\frac{1}{1+r}Q(S_{t+1}(0),D_{t+1}(0))-F<0$ for all $t\geq T$ implies that
$L(S_{t+1}(0),D_{t+1}(0))F>\pi-rF$ for all $t\geq T$. Monotonicity of $L$ and
$S_{t+1}(0)\leq S_{t}(0)$ then imply that $D_{t+1}(0)\geq D_{t}(0)$.

If there exists a threshold $D^{K}$ such that
$\lim_{t\to\infty}D_{t}(X_{t})=\infty$ when $D>D^{K}$ for any $X_{t}$, then
there are only two possible cases: either $\lim_{t\to\infty}D_{t}(0)<D^{K}$,
or $\lim_{t\to\infty}D_{t}(0)\geq D^{K}$. The first case is a contradiction
when $G$ is strictly convex increasing, as each increase in $D_{t+1}-D_{t}$
must be larger than $D_{t}-D_{t-1}$ so eventually $D_{t}$ must exceed $D^{K}$.
Only the second case is consistent with the general physical model. This
completes the $\implies$ direction. ∎

Finally, how large is the volume available to be filled? Recent analyses
estimate the maximum capacity consistent with stable orbital populations (i.e.
no Kessler Syndrome) between 200-900km altitude to be on the order of 1.8
million active satellites, assuming no debris (Lifson et al., 2022). Over the
next few decades, the total number of objects slated for launch is expected to
be on the order of 80,000 satellites (Patel, Samira and Koller, Josef S.,
2022). It is unclear whether there is sufficient demand to support hundreds of
thousands of satellites, let alone over a million. While we do not think the
maximum capacity described in the engineering literature will be realized due
to both the externalities described here and in the economic literature and
the aforementioned demand limitations, the large capacity available makes the
issue seem less one of filling the volume with satellites or debris than one
of operating in the volume becoming too costly due to risk.

##  Appendix C Calibration details

###  C.1 Data

We calibrate the economic parameters of our model using data collected by The
Space Report (Space Foundation, 2021) on the annual revenues accruing to each
sector of the space economy from 2006-2019. These data have been used in other
economic analyses of space and orbit use (Wienzierl, 2018; Rao, Burgess, and
Kaffine, 2020; Crane et al., 2020; Rao and Letizia, 2021). The data are not
ideal for our purpose as they are aggregates covering the entire space sector,
but more granular datasets describing specific LEO satellite operators’
revenues and costs are not available. To focus on revenues and costs relevant
to LEO satellite operators, we use only the variables which are plausibly
attributable to LEO satellite activities. We calculate total LEO satellite
operator revenues as the sum of the “Satellite communications” and “Earth
observation” variables, and total LEO satellite operator costs as the sum of
the “Ground stations and equipment”, “Space Situational Awareness” (SSA),
“Insurance premiums”, “Commercial satellite launch”, and “Commercial satellite
manufacturing” variables. We discard variables representing revenues to the
direct-to-home television, GNT (Geolocation, Navigation, and Timing), and
satellite radio sectors, as these are provided by satellites in higher orbits
beyond LEO. We also exclude suborbital commercial human spaceflight deposits
as they are by definition for transit to regions below orbital altitudes (e.g.
50-80 km above mean sea level). Since our data is recorded annually, we set
the period length to 1 year. We display the calculated variables in table 2.
Note that these are _not_ the revenues and costs accruing specifically to LEO
operators—a distinction not possible given our data. Rather, these variables
represent a superset of LEO operator revenues and costs, as they necessarily
include some geostationary satellites. We describe our strategy to account for
this issue during calibration in Appendix C.2.

Table 2: Economic data. Figures are in nominal billion USD. Data from Space Foundation (2021) and authors’ calculations. Year | Maximum total revenues | Maximum total costs
---|---|---
| attributable to all operators | attributable to all operators
| potentially using LEO | potentially using LEO
$2006$ | $13.800$ | $80.840$
$2007$ | $16.368$ | $92.956$
$2008$ | $18.104$ | $85.371$
$2009$ | $18.695$ | $69.270$
$2010$ | $19.570$ | $68.460$
$2011$ | $21.424$ | $83.853$
$2012$ | $22.747$ | $93.779$
$2013$ | $23.683$ | $108.199$
$2014$ | $24.002$ | $127.567$
$2015$ | $25.884$ | $87.222$
$2016$ | $26.087$ | $89.201$
$2017$ | $26.545$ | $95.857$
$2018$ | $28.420$ | $99.930$
$2019$ | $27.320$ | $119.160$

We calibrate physical parameters of our model using a kinetic gas
approximation of orbital mechanics and data from DISCOS (Letizia et al., 2017;
European Space Agency, 2021). These data describe the launch traffic, active
satellites, and tracked debris objects (i.e larger than 10 cm diameter) in the
600-650 km shell over the 2006-2020 period. These data aggregate over
different types of operators (e.g. commercial operators, civil government
operators, defense operators). We display these data in table 3, along with
the collision probability calculated from the kinetic gas approximation
assuming satellite operators avoid 99% of all collisions between satellites
and 95% of all collisions between satellites and tracked debris. Letting the
avoidance success rates be $\kappa_{SS}$ and $\kappa_{SD}$, the probability of
an unavoidable collision becomes

$L(S,D)=(1-\kappa_{SS})(1-e^{-\alpha_{SS}S})+(1-\kappa_{SD})(1-e^{-\alpha_{SD}D})-(1-\kappa_{SS})(1-\kappa_{SD})(1-e^{-\alpha_{SS}S})(1-e^{-\alpha_{SD}D}).$
(113)

Many ostensibly-non-commercial satellites are operated as joint ventures with
commercial enterprises and many commercial satellite operators serve primarily
civil government or defense customers, so we do not separate the satellite
data by operator type. Further, since all satellites contribute to debris and
collision probability regardless of their operator type, non-commercial
operators’ satellites ought to be included in the state vector. Non-commercial
operators may also contribute to the observed “occupancy elasticity”
(described precisely in the following section), further complicating efforts
to properly disentangle payoffs to different operator types from the available
data.

The DISCOS physical data also provide object characteristics such as mass and
cross-sectional area, which are necessary for the kinetic gas approximation.
We describe the details of the kinetic gas approximation of orbital mechanics
in Appendix C.3.

Table 3: Orbital traffic in the 600-650 km shell. Collision probability is rounded. Data from European Space Agency (2021) and authors’ calculations. Year | Satellites launched | Active satellites | Tracked debris | Collision probability
---|---|---|---|---
|  | satellites |  |
$2006$ | $15$ | $25$ | $211$ | $2.95\times 10^{-6}$
$2007$ | $84$ | $31$ | $275$ | $3.84\times 10^{-6}$
$2008$ | $168$ | $47$ | $273$ | $3.85\times 10^{-6}$
$2009$ | $72$ | $43$ | $393$ | $5.48\times 10^{-6}$
$2010$ | $156$ | $53$ | $444$ | $6.20\times 10^{-6}$
$2011$ | $30$ | $56$ | $411$ | $5.76\times 10^{-6}$
$2012$ | $73$ | $53$ | $429$ | $6.00\times 10^{-6}$
$2013$ | $213$ | $64$ | $454$ | $6.37\times 10^{-6}$
$2014$ | $261$ | $97$ | $484$ | $6.87\times 10^{-6}$
$2015$ | $175$ | $122$ | $495$ | $7.09\times 10^{-6}$
$2016$ | $15$ | $114$ | $494$ | $7.05\times 10^{-6}$
$2017$ | $26$ | $122$ | $525$ | $7.49\times 10^{-6}$
$2018$ | $36$ | $139$ | $506$ | $7.28\times 10^{-6}$
$2019$ | $33$ | $155$ | $543$ | $7.83\times 10^{-6}$
$2020$ | $9$ | $158$ | $626$ | $8.97\times 10^{-6}$

###  C.2 Economic calibration

To calibrate our economic model, we make three modifications to the open-
access equilibrium condition in equation (17). First, we allow the per-period
satellite payoff and cost to vary over time, i.e $\pi\to\pi_{t}$ and $F\to
F_{t}$. This changes the equilibrium condition to

$\displaystyle\pi_{t+1}$
$\displaystyle=(1+r)F_{t}-(1-L(S_{t+1},D_{t+1}))F_{t+1}$ (114)
$\displaystyle\implies L(S_{t+1},D_{t+1})$
$\displaystyle=1+\frac{\pi_{t+1}}{F_{t+1}}-(1+r)\frac{F_{t}}{F_{t+1}}.$ (115)

This form is similar to the one described in equation (17) but for the time
subscripts and term $1-(1+r)\frac{F_{t}}{F_{t+1}}$. This term represents
capital gains accruing to a period $t$ launcher from increases in the cost of
building and launching a satellite in period $t+1$. We abstract from
operators’ expectations over economic variables and assume they perfectly
forecast all $t+1$ objects.

Second, we allow the per-period satellite payoff to depend on the current
stock of satellites in orbit, i.e $\pi_{t}\to p_{t}(S_{t})$. We use a constant
elasticity form with exponential factor productivity growth, $p_{t}(S_{t})=\pi
e^{at}(1+\eta)S_{t}^{\eta}$, where $\eta$ is the “orbital occupancy elasticity
of per-period satellite payoffs”. We assume that the downstream market for
satellite outputs is competitive such that operators do not internalize
$\frac{\partial p_{t}}{\partial S_{t}}$. The equilibrium condition becomes

$\displaystyle L(S_{t+1},D_{t+1})$
$\displaystyle=1+\frac{p_{t+1}(S_{t+1})}{F_{t+1}}-(1+r)\frac{F_{t}}{F_{t+1}}.$
(116)

Third, we incorporate exogenous limited satellite lifespans to allow for
natural depreciation and replacement of satellites. Specifically, we assume
each satellite is replaced with probability $\mu$ each period. We calibrate
this value explicitly to simulate object stocks (described in the Appendix
C.3); for now, we leave this to be adjusted in the regression-based
calibration approach described below. The final equilibrium condition for our
simulations is

$\displaystyle L(S_{t+1},D_{t+1})$
$\displaystyle=1+\frac{1}{1-\mu}\frac{p_{t+1}(S_{t+1})}{F_{t+1}}-\frac{1+r}{1-\mu}\frac{F_{t}}{F_{t+1}}.$
(117)

To simulate future periods under different returns growth rate and occupancy
elasticity assumptions, we estimate the growth rate of total LEO satellite
operator costs. We estimate

$\displaystyle\log(F_{t})$
$\displaystyle=\eta_{0}^{F}+\eta_{1}^{F}t+\nu^{F}_{t},$ (118)

where $\log(F_{t})$ is the natural log of total LEO satellite operator costs,
$t$ is the year, the growth rate (the object of interest) is
$g=\exp(\eta_{1}^{F})-1$, and the regression error is $\nu^{F}_{t}$. The
estimated growth rate is roughly 2.5%, which is consistent with Crane et al.
(2020).

There are two final steps to our procedure: ensuring consistency between the
occupancy elasticity and factor productivity parameters, and accounting for
unobserved variables. To ensure consistency between the assumed elasticity and
implied orbital slot factor productivity and match the final observed value of
LEO-using sector revenues ($27.32b in 2019, see table 2), we calibrate the
factor productivity term $\pi$ in equation 25. Specifically, letting $K$ be
the observed value to match for each assumed elasticity value $\eta_{j}$,
setting $t=0$ and $S_{0}$ to the shell-specific initial condition
($S_{0}=158$), the factor productivity term $\pi_{j}$ satisfies

$\pi_{j}=\exp(\log(K)-\log(1+\eta)+\eta\log(S_{0})).$ (119)

Finally, as mentioned in the previous section, using maximum total sector
revenues and costs directly from the data in table 2 as though the data
reflects only operators in the 600-650 km shell is challenging for two
reasons. First, the data in table 2 cover all satellite operators—our variable
selection step is the only thing restricting the set of operators included in
the data. Even if we were successful in removing all operators outside of LEO
through variable selection when calculating total LEO operator revenues and
costs, the revenue and cost variables will still include operators outside the
600-650 km shell. The data aggregation implies an unobservable “shell-share”
coefficient, $s\in[0,1)$, scaling observed aggregate revenues and costs to
reflect only the portion attributable to satellites in the 600-650 km shell.
Second, theory predicts that the discount rate used by operators is a critical
parameter in determining LEO use, but this parameter is unobserved.

Fortunately, equation 115 offers a way to address both challenges. Letting
$\pi_{t}$ be the total LEO satellite operator revenues and $F_{t}$ be the
total LEO satellite operator costs from table 2, and $L_{t}$ be the collision
probability shown in table 3, we estimate the following regression on data
from 2006-2019:

$\displaystyle L_{t}$
$\displaystyle=\gamma_{0}+\gamma_{1}\frac{\pi_{t}}{F_{t}}+\gamma_{2}\frac{F_{t-1}}{F_{t}}+e_{t}.$
(120)

The estimated “adjustment coefficients” $(\gamma_{0},\gamma_{1},\gamma_{2})$
reflect the shell-share $s$, the discount rate $r$, as well as the satellite
turnover $\mu$ (though they are not separately identified). We use the
adjustment coefficients to simulate the model in future periods given
projected growth in $\pi_{t}$ and $F_{t}$. If the shell-share coefficients
(labeled $s$ in the preceding discussion) are common to revenues and costs and
time-invariant (or “close” and “slowly-varying”), they will (almost) cancel
out of the ratios we use in equation (120) and our estimated adjustment
coefficients would only reflect satellite turnover and discounting.333333This
is not the only interpretation of our estimates—as described in Rao, Burgess,
and Kaffine (2020), the adjustment coefficients may also reflect unmodeled
frictions in satellite launching and operation.

###  C.3 Physical calibration

Here we describe key equations and the ridge regression approach to correcting
for non-random object paths. Readers interested in detailed explanations of
the physics-based elements of our calibration approach, including derivations
and validation, are referred to Letizia (2016).

We require physically-appropriate values for the following parameters:
$\delta$, $\mu$, $\alpha_{SS}$, $\alpha_{SD}$, $\alpha_{DD}$, $\beta_{SS}$,
$\beta_{SD}$, $\beta_{DD}$. Calibrating $\delta$ and $\mu$ (the mean debris
decay rate and mean satellite active lifetime time) are the most
straightforward. We take data from ESA regarding the residence time
$\delta_{r}$ of debris objects and lifetime of active satellites $\mu_{r}$ at
different altitudes (European Space Agency, 2021). We set the decay rate for
debris objects as $\delta=\min\\{1-\delta_{r}^{-1},1\\}$ and the natural
turnover rate for satellites as as $\mu=\min\\{1-\mu_{r}^{-1},1\\}$. For both
parameters we calculate share-weighted averages across object types within the
category to reflect the effects of heterogeneous object dimensions, e.g.
$\delta$ reflects the weighted average of decay times for rocket bodies,
fragments, and intact derelict objects.343434Cross-sectional area and mass are
key determinants of orbital residence times. Both can vary significantly
within object classes. We calculate the share-weighted decay rate in the
600-650 km shell is roughly 7% every year.353535At these altitudes, the decay
rates from higher shells rapidly approach zero. For the 650-700 km shell, the
share-weighted average decay rate is roughly 5%, and for the 700-750 km shell
the decay rate is 3%. We therefore neglect objects entering the 600-650 km
shell from higher altitudes as they are unlikely to significantly change our
results. The share-weighted average active LEO satellite lifetime is roughly
$6.71$ years. This implies roughly $15\%$ of active satellites in LEO turn
over every year on average, i.e the fraction remaining is $1-\mu=0.85$.

We calibrate the parameters of $L$ and $G$ in two steps. First, we compute the
collision probability and new fragment formation parameters using a kinetic
gas approximation similar to the one used in Letizia et al. (2017) and
Letizia, Lemmens, and Krag (2018) as well as analytical fragmentation formulas
from Krisko (2011) and Letizia (2016) calibrated to the NASA standard breakup
model. These formulas require data on object mass and cross-sectional area,
which we obtain from DISCOS. The DISCOS parameters describe average values
across different types of active satellites and debris objects, so we compute
share-weighted averages for active satellites and debris objects. The kinetic
gas approximation implies that objects within the shell are moving randomly,
leading to our next step. Second, to adjust for the non-random motion of
objects in the shell, we regularize the expected fragmentation components of
$G$ by estimating a ridge regression on the debris law of motion using data in
table 3 and the analytically-computed parameter values. We also use this
second step to jointly estimate the launch debris parameter $m$ from the ridge
regression. We describe our procedure for calibrating the physical model
parameters in more detail in Appendix section C.3. Table 4 summarizes the
calibrated parameter values.

Table 4: Summary of calibrated parameter values for the 600-650 km shell. Values are rounded to the nearest integer or second non-zero decimal place. Parameter | Value | Notes
---|---|---
$\eta_{1}^{F}$ | 0.025 | Total costs growth parameter. Standard error is 0.009.
$\gamma_{0}$ | 3.35e-06 | Equilibrium adjustment coefficient 1 (open-access capital gains).
$\gamma_{1}$ | 2.22e-05 | Equilibrium adjustment coefficient 2 (gross satellite rate of return).
$\gamma_{2}$ | -2.67e-06 | Equilibrium adjustment coefficient 3 (open-access capital gains).
$\delta$ | 0.074 | Annual fraction of debris decaying to lower shell.
$\mu$ | 0.15 | Annual active satellite turnover rate
$\alpha_{SS}$ | 2.73e-07 | Satellite-satellite collision rate parameter.
$\alpha_{SD}$ | 2.73e-07 | Satellite-debris collision rate parameter.
$\alpha_{DD}$ | 2.78e-07 | Debris-debris collision rate parameter.
$\kappa_{SS}$ | 0.99 | Fraction of satellite-satellite collisions successfully avoided.
$\kappa_{SD}$ | 0.95 | Fraction of satellite-debris collisions successfully avoided.
$\tilde{\beta}_{SS}$ | 1,800 | Expected number of fragments from satellite-satellite collision. (regularized).
$\tilde{\beta}_{SD}$ | 333 | Expected number of fragments from satellite-debris collision. (regularized).
$\tilde{\beta}_{DD}$ | 327 | Expected number of fragments from debris-debris collision (regularized).
$m$ | 0.013 | Expected number of launch debris remaining in shell after 1 year (regularized).

To calculate the intrinsic collision probabilities
$\alpha_{SS},\alpha_{SD},\alpha_{DD}$, we start with data regarding object
cross-sectional areas for active satellites (commercial, military, civil
government, and other) and intact debris objects. We assume debris fragments
are uniform aluminium spheres of diameter 10 cm, and treat all other objects
as uniform spheres as well. We compute the cross-sectional areas of active
satellites and debris within each shell as share-weighted averages over 2006
2019 across the types of objects within each class, e.g. if 20% of the debris
objects are intact and 80% are fragments we calculate the area as
$0.2*(intact\leavevmode\nobreak\ area)+0.8*(fragment\leavevmode\nobreak\
area)$. Under these assumptions the rate at which a reference object moving
randomly at speed $s$ in a closed space of volume $V$ is struck by an object
of cross-sectional area $a$ is

$\frac{sa}{V},$ (121)

where the volume is determined by the altitude and our assumption of that the
space is a spherical shell, and the speed is determined by the altitude, the
Earth’s gravitational constant, and our assumption that the objects are
uniform spheres.

To calculate the unadjusted fragmentation rates, we use data on average object
masses from ESA along with a formula found to fit the high-fidelity NASA
standard breakup model described in Krisko (2011). Letting the mass of the
object struck be $M$, and assuming the object is shattered into uniform 10 cm
spheres, the number of fragments from a catastrophic collision $n$ is

$n=0.1M^{0.75}0.1^{-1.71}.$ (122)

The only steps remaining are to adjust our estimate of the expected number of
fragments from collisions for the non-random motion of objects in the shell,
and to set the value of the launch debris parameter $m$. ODE-based engineering
models of the debris environment use such adjustment coefficients based
fitting the ODE model to results from many computationally-costly runs of
high-fidelity orbital environment models, e.g. as in Somma et al. (2017);
Somma (2019). This approach would be even costlier for our model, as the
launch rate is endogenous, and would not provide a useful estimate of the
launch debris parameter $m$. We instead perform the adjustment and estimate
$m$ jointly using historical data and ridge regression, a regularization
technique used to improve out-of-sample predictive performance at the expense
of in-sample fit. Ridge regression achieves this goal by exploiting the bias-
variance trade-off, shrinking parameter values toward zero in exchange for
reduced prediction variance (Hoerl, Kennard, and Hoerl, 1985; Zou and Hastie,
2005).

Since satellites are specifically coordinated to reduce collisions, the
adjustment for non-random motion should involve shrinking the expected number
of fragments from a collision (with the expectation taken over the probability
of a collision) toward zero. Ridge regression achieves this goal.
Additionally, ridge regression is often used when the number of variables is
“large” relative to the number of observations or when parameter estimates are
known to be noisy due to (for example) high degrees of collinearity. Our model
and data satisfy the former condition (with 4 parameters to estimate from 14
observations), and our physical calibration approach (specifically the
assumption that all objects are uniform spheres) causes collinearity in our
collision probability values. Since our collision model prescribes the
functional form of the collision probability as $(1-\exp(-{\alpha}_{jk}k))$,
the effect of non-random motion on new debris growth cannot be separately
identified from $\alpha_{jk}$ and $\beta_{jk}$. This is convenient for our
regression-based adjustment, since it allows us to pose the ridge regression
as a linear model. Specifically, letting $\bar{x}$ denote a physically-
calibrated parameter value, we estimate the following regression:

$\displaystyle
D_{t+1}-(1-{\delta})D_{t}=\rho_{SS}{\beta}_{SS}(1-\exp(-{\alpha}_{SS}S_{t}))+\rho_{SD}{\beta}_{SD}(1-\exp(-{\alpha}_{SD}S_{t}))+\rho_{DD}{\beta}_{DD}(1-\exp(-{\alpha}_{DD}D_{t}))+m+\nu^{D}_{t},$
(123)

where $\rho_{SS},\rho_{SD},\rho_{DD},m$ are parameters to be estimated and
$\nu^{D}_{t}$ is the error term. The final regularized estimates of the
fragmentation and launch debris parameters are shown in table 4 as
$\tilde{\beta}_{SS}$, $\tilde{\beta}_{SD}$, $\tilde{\beta}_{DD}$, $m$.

##  Appendix D Algorithms for equilibrium and optimum

To describe how we generate initial guesses for the social planner’s problem,
it is useful to formally state a finite-horizon sequence version of the
planner’s problem. Letting $T$ be the final period, the planner’s finite-
horizon sequence problem is

$\displaystyle\max_{\\{X_{t},S_{t+1},D_{t+1}\\}_{t=0}^{T}}$ $\displaystyle
S_{t}Q(S_{t},D_{t},X_{t})+\frac{1}{1+r}\sum_{\tau=t}^{T-1}\frac{1}{1+r}^{\tau-t-1}X_{\tau}\left(\frac{1}{1+r}Q(S_{\tau+1},D_{\tau+1},X_{\tau+1})-F\right)$
(124) $\displaystyle\text{s.t. }Q(S_{t},D_{t},X_{t})$
$\displaystyle=\pi+\frac{1}{1+r}(1-L(S_{t},D_{t}))Q(S_{t+1},D_{t+1},X_{t+1})\text{
if }t<T$ $\displaystyle Q(S_{T},D_{T},X_{T})$ $\displaystyle=\pi$
$\displaystyle S_{t+1}$ $\displaystyle\leq S_{t}(1-L(S_{t},D_{t}))+X_{t}$
$\displaystyle D_{t+1}$ $\displaystyle\geq
D_{t}(1-\delta)+G(S_{t},D_{t})+mX_{t}$ $\displaystyle X_{t}$
$\displaystyle\in[0,\bar{X}]\leavevmode\nobreak\ \leavevmode\nobreak\ \forall
t$ $\displaystyle S_{t+1}$ $\displaystyle\geq 0,D_{t+1}\geq 0$ $\displaystyle
S_{0}$ $\displaystyle=s_{0},D_{0}=d_{0}.$

The guess generation in algorithm 1 uses a result from Easley and Spulber
(1981), that optimal plans generated from solving a finite horizon problem
with sufficiently-large $T$ closely approximate infinite-horizon optimal
plans. Algorithm 1 describes our solution procedure more precisely.

1Generate a sparse initial grid, $\mathcal{G}_{0}$, over
$(S,D)\in\mathbb{R}_{[0,a]}\times\mathbb{R}_{[0,b]},\leavevmode\nobreak\
\leavevmode\nobreak\ a,b>0$.

2At each point on $\mathcal{G}_{0}$, solve program 124 with $T$ equal to a
large number. Larger is better; we use $T=150$, which balances compute time
with guess quality. This produces an initial guess on a sparse grid,
$\tilde{v}_{0}$.

3Using linear interpolation, “infill” $\tilde{v}_{0}$ (defined on
$\mathcal{G}_{0}$) to $v_{0}$ (defined on $\mathcal{G}_{1}$).
$\mathcal{G}_{1}$ has the same boundaries as $\mathcal{G}_{0}$
($(S,D)\in\mathbb{R}_{[0,a]}\times\mathbb{R}_{[0,b]}$) but contains more
points. This gives an initial guess defined on a denser grid.

4Set $\delta$ to some large number (we use 10) and $\epsilon$ to some small
number (we use 1% of the mean value of $v_{0}$). Set $i=0$ and
$W_{0}(S,D)=v_{0}$.

5while  _$\delta >\epsilon$_ do

6 At each node in $\mathcal{G}_{1}$, solve program 22 with
$W(S_{t+1},D_{t+1})=W_{i}(S_{t+1},D_{t+1})$. Label the value function obtained
as $W_{i+1}(S,D)$, defined over $\mathcal{G}_{1}$. We use linear interpolation
to compute $W_{i}(S_{t+1},D_{t+1})$ when $(S_{t+1},D_{t+1})$ is between nodes
of $\mathcal{G}_{1}$.

7 $\delta\leftarrow||W_{i}(S,D)-W_{i+1}(S,D)||_{\infty}$.

8 i $\leftarrow$ i+1

9 end while

Algorithm 1 Solve the planner’s problem

Generating the open-access policy function is much simpler. At each node on a
grid over $S$ and $D$ values (e.g. $\mathcal{G}_{1}$ as in Algorithm 1), we
solve the open-access condition

$\displaystyle\pi-rF-L(S^{\prime}(X,S,D),D^{\prime}(X,S,D))F=0$ (125)

for the open-access launch rate $X$.

To generate the phase diagrams, we use solved policy functions $X$ to compute
the evolution of the satellite and debris stocks at each grid node. More
precisely, we compute $dS=(S^{\prime}(X,S,D)-S)/h$ and
$dD=(D^{\prime}(X,S,D)-S)/h$ for a fixed positive value $h$. The value of $h$
is chosen to make the plotting more stable; we use $h=10$, but other values
yield similar results. The nullclines are plotted as the zero-isoclines of
$dS$ and $dD$.