# Managing Space Debris in the Era of Large Constellations: A Dynamic Game Approach

Pierre Bernhard, Marc Deschamps, Georges Zaccour

<figure><img src="image_1.png" /></figure>

## To cite this version:

Pierre Bernhard, Marc Deschamps, Georges Zaccour. Managing Space Debris in the Era of Large Constellations: A Dynamic Game Approach. 2025. hal-05485645

HAL Id: hal-05485645

https://hal.science/hal-05485645v1

Preprint submitted on 9 Feb 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Managing Space Debris in the Era of Large Constellations: A Dynamic Game Approach*

Pierre Bernhard

Centre INRIA de l'Université de la Côte d'Azur, France

Marc Deschamps

CRESE UR3190, Université Marie et Louis Pasteur, Besançon, France

Georges Zaccour

GERAD, HEC Montréal, Canada

December 11, 2025

## Abstract

Space infrastructure and data have become indispensable to humanity at all scales, which is largely due to the advent of the New Space era, characterized by the creation of large satellite constellations. Unfortunately, these remarkable developments exacerbate an already pressing issue: the growing accumulation of space debris. In this article, we propose a dynamic game-theoretic model showing that this problem can be contained, if an international institution is created and empowered to levy taxes to finance active debris removal activities.

**Keywords:** Dynamic games; space debris, large constellations, active debris removal, international institution.

## 1 Introduction

Since its inception in the late 1950s, the space age has never been as advanced as it is today. Every day, humanity relies on space-based infrastructure and data at local, national, and international levels. Whether in communications, transportation, agriculture, security, defense, health, environmental monitoring, or other fields, space has become indispensable to science, commerce, and public policy.

This transformation, often referred to as New Space or Space 2.0, emerged in the early 21st century, driven by political decisions and public investment, scientific progress, advances in microelectronics and satellite miniaturization, reduced launch costs, digital innovation, and private sector engagement. Once the preserve of a handful of nations and a few large specialized corporations,

<footnote>*The authors would like to thank the Centre National d'Etudes Spatiales, France, for its financial support within project PERSEDS.</footnote>
space has become increasingly accessible: today, only a few countries lack a space policy, and thousands of private companies are active across the sector.

One of the most striking developments in the space sector, particularly since the 2020s, has been the emergence and expansion of large-scale projects, mainly in low Earth orbit (LEO), between 100 and 2,000 kilometers in altitude.<fnref n="1" /> These projects include several large constellations, i.e., groups of more than 100 satellites operating as an integrated system, to provide high-speed, low-latency Internet connectivity and global communications coverage.<fnref n="2" /> Prominent examples include *Starlink* (currently about 7,400 satellites, with authorization for up to 42,000), *Guowang* (expected to comprise 13,000 satellites), *Kuiper* (around 10,000 satellites), and *IRIS*<fnref n="2" /> (several hundred LEO and MEO satellites).<fnref n="3" /> To grasp the scale of this transformation, it is useful to recall that the major navigation constellations (GPS, GLONASS, Galileo, and BeiDou) each consist of roughly 30 satellites, and that as of May 2025, there were fewer than 12,000 active satellites in orbit. The magnitude of this change is also reflected in the rapid increase in annual launch activity: in 2024 alone, a record 250 launches placed more than 2,500 spacecraft into orbit (see Figure 1).

<figure><img src="image_3.png" /><figcaption>Figure 1: Evolution of number of launches (Liou (2025,p.3))</figcaption></figure>

Unfortunately, these new constellations are emerging in an already alarming context with regard to space debris.<fnref n="4" /> Since the dawn of the space age, humans have left debris in orbit, both during satellite launches and following the end of satellites’ operational lifetimes, particularly in LEO. This accumulation has led to severe congestion, increasing the likelihood of collisions between active satellites and existing debris, thereby generating even more fragments. The danger arises, first, from the extremely high orbital velocities involved—between 6.9 and 7.9 km/s in LEO—meaning that even debris as small as 1-10 centimeters can completely destroy a satellite and create additional debris. Second, the persistence of debris in orbit exacerbates the problem: an object at an altitude of 1,200 kilometers may take up to 4,300 years to naturally re-enter the atmosphere, whereas debris at 600 kilometers typically deorbits within about 10 years.

In response to this growing threat, major space agencies, governments, and some private operators have sought since the late 1980s to establish guidelines and standards aimed at mitigating

<footnote n="1">¹Liou (2025, p. 5) notes that more than half of the total orbital mass is located in LEO.</footnote>

<footnote n="2">²According to IADC (2025a, p. 8), a constellation is “a set of spacecraft operating in a coordinated manner to create a single large system, often based on a recurring spacecraft design and operating in similar orbits.”</footnote>

<footnote n="3">³For a detailed overview of large constellations, see GAO (2022) and Académie des Sciences (2024)</footnote>

<footnote n="4">⁴It should be noted that, for the first time, the IADC [2025a] has included a dedicated section on large constellations.</footnote>
or eliminating debris generation, e.g., the ESA's PSS-01-40 policy introduced in 1988. Yet, despite these initiatives, progress has been limited: the overall volume of debris continues to rise, particularly in LEO, as illustrated in Figure 2.

<caption>Figure 2: Changes in LEO environment (Liou (2025, p. 6))</caption>

We are interested in this paper by the management of space debris and wish to address two research questions:

1. Under what conditions, the economic exploitation of space remains profitable?

2. Can satellite operators benefit from coordinating their launching activities?

The first question is related to the so-called Kessler syndrome (Kessler and Cour-Palais (1978), Kessler (1991)), which predicts that the multiplication of debris beyond a certain level would render space physically unexploitable. As pointed out by Adilov et al. (2018) an economic Kessler syndrome can preclude the physical one, that is, the cost induced by these debris would leave firms with no positive profits. Now, suppose that an international agency is created and responsible for (i) actively removing debris to keep the stock constant over time, and (ii) levying taxes on satellite operators to finance its operations. This raises the question of whether the level of taxation required to implement the policy is compatible with operators maintaining strictly positive profits.

By the virtue of joint optimization, coordinating of strategies (cooperation) leads to a higher total collective payoff than the sum of noncooperative individual gains. Our second question is meant to check if the difference is sufficiently high to justify (probably lengthy and costly) negotiations for designing an agreement between the operators.

Our paper is related to Bernhard et al. (2023), with two notable differences. First, we consider any number of constellations, rather than at most two operating at approximately the same altitude. Second, we note that when constellations are launched into orbits with significant eccentricity, they may interfere with a larger number of similar systems than would occur with quasi-circular orbits. To capture this possibility, we modify the original model to allow for multiple orbit types to represent the relative collision frequencies between two constellations.

Our main finding is that the core conclusion of Bernhard et al. (2023) still holds: the difference between jointly optimal and Nash equilibrium payoffs remains small, as do the differences in the
associated strategies. Consequently, the incentive to deviate from the collectively optimal solution is small.

Although the economic analysis of outer space has roots dating back fifty years, with early contributions by Snow (1975), O'Neill (1977), Sandler and Schulze (1981), and Wihlborg and Wijkman (1981), it is, as Bongers et al. (2024) observe, still premature to speak of a distinct and mature field of “space economics.” At present, the most comprehensive overview of the subject is provided by Weinzierl (2018), while Oltrogge and Christensen (2020) offer a broader perspective on space governance, and Bongers et al. (2024) present the most detailed literature review. Our article contributes to the more focused emerging body of work dealing with the financing of space debris removal through a tax on satellite launches, within a framework of imperfect competition modeled as a dynamic game. In this regard, beyond Bernhard et al. (2023), the study most closely related to ours is Guyot and Rouillon (2023) (extending Rouillon (2020)), which analyzes how the combined use of an ad valorem tax, a launch tax, and a certificate trading system can generate appropriate incentives to ensure the long-term sustainability of low Earth orbits.

The rest of the paper is organized as follows. Section 2 introduces the model, and Section 3 characterizes and compares the noncooperative and cooperative solutions. Section 4 concludes. The proofs of the propositions are in the appendix, as the notation of our variables and their orders of magnitude.

# 2 Model

## 2.1 Notation

Our model is meant to be flexible and adaptable to various situations, which comes with the cost of having a large number of parameters. And to keep our calculations manageable, we need to introduce further notation as the combinations of native coefficients.

**Vector notation** For two $n$-vectors $q$ and $r$, the notation $\langle q, r \rangle$ stands for their inner product:

$$ \langle q, r \rangle = \sum_{i=1}^{n} q_i r_i = q^t r. $$

The notation $\sum_i$, with no explicit bounds on $i$ will stand for $\sum_{i=1}^n$, and likewise for $\sum_j$.

$\mathbb{1}$ stands for the $n$-vector whose coordinates are all equal to 1 while $\mathbb{1}^i$ stands for the $n$-vector whose only non-zero coordinate is the $i$-th one: $\mathbb{1}^i = 1$ and $\forall j \neq i$, $\mathbb{1}^j = 0$. ($\mathbb{1}^i$ is therefore the $i$-th basis vector of $\mathbb{R}^n$.)

For each scalar quantity $q_i$ indexed by $i \in \{1, \dots, n\}$, the letter $q$ stands for the $n$-vector of the $q_i$, $q^i$ stands for the $n$-vector whose only nonzero entry is the $i$-th: $q_i^i = q_i$, $q_j^i = 0$ for $i \neq j$, and $q^t$ and $(q^t)^t$ respectively denote their transposed.

We will introduce coefficients $z_{i,0}$ and $z_{i,j}$ related to collision risks. We denote by $z^0$ the $n$-vector of the $z_{i0}$ and $z^i$ the $n$-vector of the $z_{ij}$, $j$ ranging from 1 to $n$, whose $j$-th coordinate is $z_{ij}$.
**Matrix notation** Let Z be the symmetric n × n-matrix of the zij, with zero diagonal elements.
We will also use the diagonal matrices

$$
A = \operatorname{diag}(a_i), \quad \Psi = \operatorname{diag}(\psi_i), \quad M = \operatorname{diag}(\mu_i).
$$

Finally, we introduce two sets of symmetric n × n matrices: Y^i whose only non-zero entry is
Y_{ii}^i = 1, and Z^i with z^i in column i and (z^i)^t in line i:

$$
Z^i = z^i(\mathbb{1}^i)^t + \mathbb{1}^i(z^i)^t =
\begin{pmatrix}
0 & \cdots & 0 & z_{i1} & 0 & \cdots & 0 \\
\vdots & & \vdots & \vdots & \vdots & & \vdots \\
0 & \cdots & 0 & z_{ii-1} & 0 & \cdots & 0 \\
z_{i1} & & z_{ii-1} & 0 & z_{ii+1} & & z_{in} \\
0 & \cdots & 0 & z_{ii+1} & 0 & \cdots & 0 \\
\vdots & & \vdots & \vdots & \vdots & & \vdots \\
0 & \cdots & 0 & z_{in} & 0 & \cdots & 0
\end{pmatrix}
$$

**Parameters** Table 1 gives the list of variables, functions, and parameters

<table><caption>Table 1: Variables, functions, and parameters</caption>
  <tr>
    <td>s<sub>i</sub>(t)</td>
    <td>Number of satellites launched by operator <i>i</i> at year <i>t</i>.</td>
  </tr>
  <tr>
    <td>x<sub>i</sub>(t)</td>
    <td>Number of satellites of operator <i>i</i> aloft at a time <i>t</i>.</td>
  </tr>
  <tr>
    <td>y(t)</td>
    <td>Number of debris aloft at time <i>t</i>, with <i>y</i>(0) = <i>y</i><sub>0</sub></td>
  </tr>
  <tr>
    <td>h(t)</td>
    <td>Number of actively removed debris via ADR the year <i>t</i>.</td>
  </tr>
  <tr>
    <td>R<sub>i</sub>(x<sub>i</sub>)</td>
    <td>= &phi;<sub>i</sub>x<sub>i</sub> + &psi;<sub>i</sub>x<sub>i</sub><sup>2</sup>/2, revenue function.</td>
  </tr>
  <tr>
    <td>C<sub>i</sub>(s<sub>i</sub>)</td>
    <td>= c<sub>i</sub>s<sub>i</sub>, cost of launching s<sub>i</sub>, with c<sub>i</sub> &gt; 0.</td>
  </tr>
  <tr>
    <td>T<sub>i</sub>(x<sub>i</sub>)</td>
    <td>= &eta;<sub>i</sub>x<sub>i</sub> + &zeta;<sub>i</sub>&omega;, tax charged to operator <i>i</i>.</td>
  </tr>
  <tr>
    <td>&Pi;<sub>i</sub></td>
    <td>Operator <i>i</i>'s profit.</td>
  </tr>
  <tr>
    <td>&eta;<sub>i</sub></td>
    <td>Tax rate</td>
  </tr>
  <tr>
    <td>&zeta;<sub>i</sub></td>
    <td>Share of the non-attributable costs &omega; borne by player <i>i</i>, with &Sigma;<sub>i</sub> &zeta;<sub>i</sub> = 1.</td>
  </tr>
  <tr>
    <td>p</td>
    <td>Unit cost of active debris removal (ADR).</td>
  </tr>
  <tr>
    <td>&alpha;<sub>i</sub></td>
    <td>Relative rate of deorbiting either by natural decay or active deorbiting.</td>
  </tr>
  <tr>
    <td>&alpha;<sub>0</sub></td>
    <td>Relative rate of debris deorbiting including arrivals from higher orbits.</td>
  </tr>
  <tr>
    <td>&beta;<sub>i</sub></td>
    <td>Relative rate of on-orbit satellite death.</td>
  </tr>
  <tr>
    <td>&nu;</td>
    <td>Average number of debris caused by a collision.</td>
  </tr>
  <tr>
    <td>&mu;<sub>i</sub></td>
    <td>Half the cost of an evasive maneuver for operator <i>i</i>.</td>
  </tr>
  <tr>
    <td>&rho;</td>
    <td>Discount factor.</td>
  </tr>
  <tr>
    <td>&sigma;</td>
    <td>Expected number of satellites launched outside the agreement.</td>
  </tr>
  <tr>
    <td>&tau;</td>
    <td>“Congestion parameter”:<br/>&tau;z = probability of hitting one among <i>z</i> objects per unit time.</td>
  </tr>
  <tr>
    <td>&phi;<sub>i</sub></td>
    <td>Unit yearly value of a satellite aloft in revenue function R(x<sub>i</sub>).</td>
  </tr>
  <tr>
    <td>&psi;<sub>i</sub></td>
    <td>Coefficient of the term in x<sub>i</sub><sup>2</sup> in the revenue function R(x<sub>i</sub>).</td>
  </tr>
</table>
## 2.2 Dynamic model

We consider *n* large constellations sharing the same altitudes, each with a number $x_i$ of satellites, and *y* debris in the same altitude interval. These numbers will vary with time according to equations to be described below.

**Collision risk and maneuvers** Space debris and collision risks are our main concerns. The frequency of “potential collisions” of *n* orbital objects of a certain type with *m* of another type is modelled by a term à la Lotka-Volterra $\tilde{\tau}nm$, and $(1/2)\tilde{\tau}n^2$ for objects of the same type, where the $\tilde{\tau}$ are parameters to be described now.

The number of collisions of debris with debris will be denoted $(1/2)\tilde{\tau}y^2$, where $\tilde{\tau}$ is a small parameter (typically of order $10^{-6}$ to $10^{-7}$). The rate of collisions of $x_i$ satellites with $y$ space debris is $\tilde{\tau}_{i0}xy$. Concerning two constellations of $x_i$ and $x_j$ satellites respectively, we will assume that collisions are avoided by Just-in-time Collision Avoidance (JCA) maneuvers. The frequency of these maneuvers depends on the rate of “potential collisions”, and given by $\tilde{\tau}_{ij}x_ix_j$. The parameter $\tilde{\tau}_{ij}$ depends on the relative orbital parameters of the two constellations, assuming they share the same altitudes. These maneuvers are costly, because they expand orbit-keeping propellant, thus limiting the usual life time of the satellites. If the “cost” of a maneuver for operator *i* is $2\mu_i$, each collision risk detected, will, on average, cost $\mu_i$, because only one of the two satellites involved in the collision risk will maneuver.

We set $\tilde{\tau}_{ij} = z_{ij}\tilde{\tau}$, and we expect that $z_{ij} := \tilde{\tau}_{ij}/\tilde{\tau}$ is never large. (See below.) Obviously, $\forall i \neq 0$, $z_{ii} = 0$, and we take the $z_{ij} = z_{ji}$ as given. But determining a reasonable number for each of them is a research topic in itself (see, e.g. [12]), which we do not address here. We note however that we mean here collision risk. Such events are much more frequent than actual collisions in the absence of maneuvering. If we assume that a maneuver is initiated each time the perceived probability of collision is larger than a given threshold $P_0$, then the $z_{ij}$ for $i \neq 0$ should be of the order of $1/P_0$. Assuming $P_0$ is the same for all combinations, we have rather absorb this coefficient into the costs $\mu_i$, to keep $z_{ij}$ of the order of a few units.

### 2.2.1 State dynamics

We adopt an infinite-horizon discrete-time dynamic model. The time step is one year, which is what it takes to decide, build, and launch a set of satellites. For $i \in \{1, 2, \dots, n\}$, let $x_i(t)$ be the number of satellites in constellation *i*, $s_i(t)$ the number of launched satellites, and $y(t)$ the number of space debris in the same region of space at time period *t*. At each time period, we assume the following:

* A loss of $\alpha_i x_i(t)$ satellites by (passive and active) deorbitation.

* A number of $\beta_i x_i(t)$ satellites are transformed into space debris by accidental breakdown.

* A loss of $\alpha_0 y$ debris by natural decay and replenishment by debris decaying from higher orbits.
• A loss of $\tau z_{i0}x_i(t)y(t)$ satellites by collisions with space debris, each collision producing $\nu$ new space debris.

• Collisions between satellites of two of the constellations are systematically avoided by JCA.

• A launch of $s_i(t)$ satellites in that constellation.

• A mean launch of $\sigma$ satellites by operators not party in the agreement.

Consequently, we have

$$x_i(t+1) = (1 - \alpha_i - \beta_i - \tau z_{i0}y(t))x_i(t) + s_i(t), \quad x_i(0) = 0, \quad (1)$$

$$y(t+1) = (1 - \alpha_0)y(t) + \langle \beta, x \rangle + \tau \left[ (\nu - 1)\langle z^0, x \rangle + \frac{y(t)}{2}(\nu - 2) \right] y(t) + \sigma - h(t), \quad (2)$$

$$y(0) = y_0. \quad (3)$$

As in Bernhard et al. (2023), we assume that the number of debris must remain constant. More specifically, we impose $y(t) = y_0$, which then leads to

$$h(t) = -\alpha_0 y_0 + \tau(\nu - 2)\frac{y_0^2}{2} + \sigma + \langle \beta + (\nu - 1)y_0\tau z^0, x \rangle,$$

Letting

$$a_i = 1 - \alpha_i - \beta_i - \tau z_{i0}y_0, \quad (4)$$

$$\omega = -\alpha_0 y_0 + \tau(\nu - 2)\frac{y_0^2}{2} + \sigma, \quad (5)$$

$$\eta_i = \beta_i + \tau z_{i0}(\nu - 1)y_0, \quad (6)$$

where $\omega$ is the yearly debris creation that is not attributable to the $n$ players, the state dynamics become

$$x_i(t+1) = a_i x_i(t) + s_i(t), \quad \text{or, equivalently,} \quad x(t+1) = Ax(t) + s(t), \quad (7)$$

and

$$h(t) = \omega + \langle \eta, x \rangle. \quad (8)$$

For the dynamics of $x$ to make sense, the parameter values in (4) must be such that $a_i > 0$ for all $i$.

### 2.2.2 Optimization problems

Denote by $R_i(x_i)$ Player $i$'s revenues from the exploitation of a fleet of $x_i$ satellites. We assume that $R_i(x_i)$ is given by the following concave function:

$$R_i(x_i) = \phi_i x_i - \frac{\psi_i}{2} x_i^2,$$
where $\phi_i$ and $\psi_i$ are positive parameters. The cost of maneuvering to avoid collisions is defined by $\tau\mu_i \sum_j z_{ij}x_j$, where $\tau$ and $\mu_i$ are positive parameters. Let $c_i s_i$ be the cost of launching $s_i$ satellites during one period of time, and $T_i(x_i)$ be the contribution of player $i$ to ADR cost given by

$$T_i(x_i) = \eta_i x_i + \zeta_i \omega,$$

where the individual shares $\zeta_i$ of the non-attributable costs satisfy $\langle 1, \zeta \rangle = 1$. How the $\zeta_i$ are chosen is a classical problem in public good economics, and will not be discussed at this point. Denoting by $p$ the unit cost of ADR, then Player $i$'s total discounted profit is then given by

$$\Pi_i = \sum_{t=0}^{\infty} \rho^t \left[ R_i(x(t)) - c_i s_i - \tau \mu_i \sum_j z_{ij} x_j - p T_i(x_i(t)) \right], \quad (9)$$

where $\rho \in (0, 1)$ is the common discount factor.

By (7) and (9), we have defined a $n$-player infinite-horizon discrete-time dynamic game, with state variable $x$ and one control variable $s_i$ for each player. To further save on notation, let

$$b_i = \phi_i - p\eta_i = \phi_i - p\beta_i - \tau z_{i0} p(\nu - 1)y_0, \quad (10)$$

$$g_i = \rho\phi_i - (1 - \rho a_i)c_i, \quad (11)$$

$$f_i = \rho b_i - c_i(1 - \rho a_i) = g_i - p\rho\eta_i. \quad (12)$$

Letting $b^i = b_i 1^i$, $c^i = c_i 1^i$, and $Q^i = \psi_i Y^i + \tau \mu_i Z^i$, Player $i$'s profit can be rewritten as:

$$\Pi_i = \sum_{t=0}^{\infty} \rho^t \left[ \langle b^i, x \rangle - \frac{1}{2} \langle x, Q^i x \rangle - \langle c^i, s \rangle - p \zeta_i \omega \right]. \quad (13)$$

Further, let

$$\bar{U} = \sum \mu_i Z^i = MZ + ZM, \quad \bar{Q} = \sum_i Q^i = \Psi + \tau \bar{U},$$

where $\bar{U}$ is the matrix of the $u_{ij} = (\mu_i + \mu_j)z_{ij}$. Using the product $\bar{\psi} = \prod_i \psi_i$, we note that

$$\det \bar{Q} = \bar{\psi} - \tau^2 \sum_{ij} \frac{\bar{\psi}}{\psi_i \psi_j} (\mu_i + \mu_j)^2 z_{ij}^2 + O(\tau^3).$$

As a consequence, and since $\tau$ is of the order of $10^{-6}$ to $10^{-7}$, while the $\psi_i$ are of the order of $10^{-3}$, the $\mu_{ij}$ and $z_{ij}$ of a few units, this determinant is not null, and $\bar{Q}$ is invertible.

The matrix $\bar{U}$ has all its entries positive. According to the Perron-Frobenius theorem, the spectral radius of the matrix $\tau \bar{U}$ is of the order of $\tau$; hence (much) smaller than one. Therefore, the inverse $\bar{Q}^{-1}$ may be represented as a series in the powers of $\tau$, allowing one to write:

$$\bar{Q} = \Psi + \tau \bar{U} = (I + \tau \bar{U} \Psi^{-1}) \Psi \quad \text{hence} \quad \bar{Q}^{-1} = \Psi^{-1} - \tau \Psi^{-1} \bar{U} \Psi^{-1} + O(\tau^2). \quad (14)$$

To illustrate, for $n = 3$ we have

$$\delta = \det \bar{Q} = \psi_1 \psi_2 \psi_3 - \tau^2 (\psi_1 u_{23}^2 + \psi_2 u_{31}^2 + \psi_3 u_{12}^2) + 2\tau^3 u_{12} u_{23} u_{31}, \quad (15)$$
and

$$ \bar{Q}^{-1} = \frac{1}{\delta} \left[ \begin{pmatrix} \psi_2 \psi_3 & 0 & 0 \\ 0 & \psi_3 \psi_1 & 0 \\ 0 & 0 & \psi_1 \psi_2 \end{pmatrix} - \tau \begin{pmatrix} 0 & \psi_3 u_{12} & \psi_2 u_{31} \\ \psi_3 u_{12} & 0 & \psi_1 u_{23} \\ \psi_2 u_{31} & \psi_1 u_{23} & 0 \end{pmatrix} + \tau^2 \begin{pmatrix} -u_{23}^2 & u_{23} u_{31} & u_{12} u_{23} \\ u_{23} u_{31} & -u_{31}^2 & u_{12} u_{31} \\ u_{12} u_{23} & u_{12} u_{31} & -u_{12}^2 \end{pmatrix} \right]. $$

# 3 Solutions

In this section, we characterize and compare the jointly optimal solution and Nash equilibrium.

## 3.1 Joint optimal solution

Let

$$ \xi^O = \frac{1}{\rho} \bar{Q}^{-1} f. $$

The superscript O stands for (jointly) optimal solution.

**Proposition 1** *Assuming that all $\xi_i^O$ are non-negative, the jointly optimal satellite launching policy for $i = 1, \dots, n$ is given by*

$$ s_i^O(x_i) = -a_i x_i + \xi_i^O, \qquad (16) $$

and the steady-state value by

$$ x_{i\infty}^O = \xi_i^O. \qquad (17) $$

The value function is given by

$$ V^O(x) = -\frac{1}{2}\langle x, \bar{Q}x \rangle + \langle b + A^t c, x \rangle + \frac{1}{1-\rho} \left[ \frac{1}{2\rho} \langle f, \bar{Q}^{-1}f \rangle - p\omega \right]. \qquad (18) $$

**Proof.** See Appendix B. ■

Let

$$ p^N := \min_i \frac{g_i}{\rho \eta_i} = \min_i \frac{\rho \phi_i - (1 - \rho a_i) c_i}{\rho \eta_i}. \qquad (19) $$

**Proposition 2** *A necessary condition for $s_i^O(x_i)$, $x_{i\infty}^O$, and $V^O(x)$ to be nonnegative is to have $p \le p^N$.*

**Proof.** Condition (19) is equivalent to have $f_i \ge 0$ for all $i$. Note that all entries of $\bar{Q}$ are positive. Therefore, if $\xi^O$ has also all its coordinates nonnegative, so has the product $f = \rho \bar{Q} \xi^O$. Hence $f_i \ge 0$ for all $i$ is a necessary condition. We may add that if $\xi^O$ is with non-negative entries and is not the zero vector, then all the $f_i$ are strictly positive. ■

To give a hint about the satisfaction of the condition $(\bar{Q}^{-1}f)_i \ge 0$, let us consider a case with 3 symmetric players and assume that $f \ge 0$. According to the formulas given above, we obtain the condition

$$ \psi^2 - 4\tau\psi\mu + 4\tau^2\mu^2 \ge 0 \Leftrightarrow \psi(\psi - 4\tau\mu) + 4\tau^2\mu^2 \ge 0. $$
A sufficient (not necessary) condition for the above inequality to hold is to have $\psi \ge 4\tau\mu$. As the order of magnitude<fnref n="5" /> of $\psi$ and $\tau$ are $10^{-3}$ and $10^{-7}$, respectively, and $\mu < 1$, the condition $\psi \ge 4\tau\mu$ is clearly always satisfied with a wide margin. Therefore, in this case, the only constraining condition is $p \le p^N$.

**Profitability** To investigate the conditions under which space exploitation is profitable under joint optimization, we consider the case with $n$ symmetric players. Then, the solution is profitable when the value function evaluated at $x = 0$ is positive, i.e.,

$$V^O(0) = \frac{1}{1-\rho} \left[ \frac{1}{2\rho} \langle f, \bar{Q}^{-1}f \rangle - p\omega \right] > 0, \quad (20)$$

which is equivalent to

$$\mathcal{P}_n^O(p) = \rho^2 \eta^2 \langle \mathbf{1}, \bar{Q}^{-1} \mathbf{1} \rangle (p^N - p)^2 - 2\rho p \omega > 0.$$

A necessary (not sufficient) condition for the above condition to hold is

$$\langle \mathbf{1}, \bar{Q}^{-1} \mathbf{1} \rangle > 0. \quad (21)$$

We note that $\mathcal{P}_n^O(p^N) = -2\rho p\omega$ is negative. Therefore, condition (21) suffices to ensure that there exists a positive $p_n^O$ such that, for $p < p_n^O$, which implies $p < p^N$, the industry is profitable.

Furthermore, using (14), we see that (21) may be written

$$\frac{n}{\psi} \left( n - 2\tau \frac{n-1}{\psi} \right) + O(\tau^2) > 0,$$

which is always true given the orders of magnitude provided in Appendix A. Hence, in practice, as long as the condition $g_i = \rho\phi_i - (1 - \rho a_i)c_i > 0$ is satisfied, the added condition $p < p_n^O$ suffices to ensure both the existence and the profitability of the optimal strategies.

## 3.2 Nash equilibrium

We seek a Nash equilibrium for the profits (13). We need to introduce the matrix $\tilde{Q}$ whose line $i$ is the line $i$ of $Q^i$. Let

$$\tilde{U} = MZ, \quad \tilde{Q} = \Psi + \tau \tilde{U}.$$

We note that

$$\Delta := \det \tilde{Q} = \bar{\psi}[1 + O(\tau^2/\psi^2)],$$

so that with the orders of magnitude of Appendix A, $\Delta > 0$ and $\tilde{Q}$ is invertible. In a similar fashion as we did in the joint optimization case, we may write $\tilde{Q} = (I + \tau \tilde{U} \Psi^{-1}) \Psi$ and hence $\tilde{Q}^{-1} = \Psi^{-1} - \tau \Psi^{-1} \tilde{U} \Psi^{-1} + O(\tau^2)$.

Let

$$\xi_i^N = \frac{1}{\rho} \tilde{Q}^{-1} f.$$

<footnote n="5">⁵See Appendix A.</footnote>
**Proposition 3** Assuming that all $\xi_i^N$ are non-negative, the unique feedback-Nash satellite launching equilibrium strategy for $i = 1, \dots, n$ is given by

$$s_i^N(x_i) = -ax_i + \xi_i^N, \qquad (22)$$

and the steady-state values by

$$x_{i\infty}^N = \xi_i^N. \qquad (23)$$

The value functions are as follows:

$$V_i(x) = -\frac{1}{2}\langle x, Q^i x \rangle + (b_i + a_i c_i)x_i + \frac{1}{\rho(1-\rho)}\left[\langle f^i, \tilde{Q}^{-1}f \rangle - \frac{1}{2}\langle \tilde{Q}^{-1}f, Q^i \tilde{Q}^{-1}f \rangle - \rho p \zeta_i \omega\right]. \qquad (24)$$

**Proof.** See Appendix B. ■

As in the joint optimization case, the equilibrium strategy is at each step an “impulse” bringing the state to its equilibrium value $\xi^N$, which is a strong form of turnpike.<fnref n="6" />

As in the joint optimization case, we have the following

**Proposition 4** A necessary condition for $s_i^N(x_i)$, $x_{i\infty}^N$, and $V^N(x)$ to be nonnegative is to have $p \le p^N$.

**Proof.** Similar to the proof of Proposition 2. ■

**Profitability** Again, we consider the profitability in a symmetric setup, i.e., with all $z_{ij} = 1$. The condition $V_i(0) > 0$ translates into

$$\mathcal{P}_n^N(p) := \rho^2 \eta^2 \left[ \langle \mathbb{1}^i, \tilde{Q}^{-1} \mathbb{1} \rangle - \frac{1}{2} \langle \tilde{Q}^{-1} \mathbb{1}, Q^i \tilde{Q}^{-1} \mathbb{1} \rangle \right] (p^N - p)^2 - \rho \frac{1}{n} \omega p > 0.$$

We can see that the condition $\mathcal{P}_n^N(0) > 0$ is more restrictive than condition (21) of the joint optimization problem. In the same fashion, if it is satisfied, there exists a bound $p_n^N < p^N$ and $p < p_n^N$ that ensure existence and profitability of the equilibrium strategies.

## 3.3 Comparison

Now, we compare the steady-state values and the payoffs obtained in the two solutions. Recall that the numbers of satellites in the steady state are given by

$$x_{\infty}^{O} = \bar{Q}^{-1}f, \quad x_{\infty}^{N} = \tilde{Q}^{-1}f.$$

Using the first order expansion in powers of $\tau$, we have

$$x_{\infty}^{N} - x_{\infty}^{O} = \tau \Psi^{-1}(\bar{U} - \tilde{U}) \Psi^{-1}f = \tau \Psi^{-1} Z M \Psi^{-1}f + O(\tau^2).$$

<footnote n="6">⁶The turnpike principle in optimal control theory states that the optimal trajectory spends most of its time near a steady-state solution. The analogy comes from the idea that a traveller takes a turnpike to a highway as soon as possible after starting their journey and leaving it at the latest moment to reach their final destination.</footnote>
The first term is made of positive coefficients, hence has all its coordinates strictly positive, of the order of $\tau\mu/\psi^2$, i.e., $10^{-1}$ to a few units. The remainder is of the order of $\tau^2/\psi^3$, i.e., $10^{-3}$ at most. Therefore, the number of satellites per constellation in the long term regime is larger in the case of the Nash equilibrium, but may be in a insignificant way (if this “theoretical” difference is less than one).

The dividend of cooperation (DC), measured by the difference between the jointly optimal solution and Nash equilibrium is defined as follows:

$$DC = V^O(0) - \sum_{i=1}^{n} V_i^N(0).$$

**Proposition 5** In the symmetric case, if $p \le p^N$ so that the Nash equilibrium yields a positive profit for the players, then the jointly optimal solution is also profitable.

In the general case, the dividend of cooperation is positive and is small (in the order of $\tau^2$).

**Proof.** See Appendix B.

<figure><img src="image_4.png" /></figure>

The two comparative results were also obtained in Bernhard et al. (2023) with two constellations. The reason is due to the fact that the players’ payoffs are lightly coupled.

# 4 Conclusion

This paper develops a model to study the economic trade-off between expanding the size of mega-constellations and the resulting loss of profitability due to space congestion, whether arising from collisions with debris or from interference among coexisting constellations. The model is deliberately simple and, as a result, necessarily incomplete. This is the cost of gaining an analytical framework that can explore policy choices more transparently and more efficiently than repeated large-scale simulations.

Our goal, however, is a model with enough adjustable parameters to be calibrated for a wide range of scenarios. In particular, it must accommodate all relevant combinations of orbital regimes. We also recognize that constellations can interfere in additional ways, notably through radio-bandwidth spillover. These effects are diverse and technically intricate, and we set them aside in order to keep the focus on debris and collision risks.

Several limitations follow from these modeling choices. We mention at least five. First, we treat the $x_i$, $s_i$ and $y$ as real numbers even though they are integers; conducting the type of analysis we seek using integer quantities is infeasible, and the approximation seems reasonable when the populations are large. Second, although our framework allows constellations to differ in their orbits, we assume that each constellation is internally homogeneous, consisting of a single orbital population. Modern constellations are typically more heterogeneous. Third, operators generally deorbit satellites after a prescribed lifetime, but a model tracking the deorbiting of thousands of individual satellites is intractable. We therefore rely on a constant rate capturing both natural decay and active disposal. Fourth, the number of debris fragments produced in a collision depends heavily on the geometry and physical composition of the objects involved. Since capturing that variability is out of reach, we employ a fixed debris yield $\nu$. Fifth, we adopt a concave economic value function for constellation size, even though a more realistic representation would be sigmoidal: a
single satellite, or even a small handful, yields virtually no economic value. Because the system's dynamics are naturally linear, we forgo this realism and keep a non-homogeneous quadratic form for the sake of analytical tractability.

Despite these simplifications, we hope the model isolates key qualitative insights, and that the conclusions may guide practical policy choices by operators and strengthen the case for establishing an international institution responsible for space debris removal.
# 5 Appendix A: Orders of magnitude

As we used first order expansions in the small parameter $\tau$ in our derivations, we need to have a (rough) estimate of the order of magnitude for the other parameters. Table 2 provides the information.

<table><caption>Table 2: Order of magnitude of the various parameters.</caption><thead><tr><th></th><th>Nature</th><th>Order</th><th>Comments</th></tr></thead><tbody><tr><td>τ</td><td>Congestion coef.</td><td>10<sup>-7</sup></td><td>Based on literature, e.g., [12], and empirical data.</td></tr><tr><td>α<sub>i</sub></td><td>Satellite decay rate</td><td>.3</td><td>depends on the altitude of the constellations considered. If we consider that after 5 years, 80% of the satellites have left the orbital altitude used (ignoring the possibility of periodic raising of it) this leads to α<sub>i</sub> = .275.</td></tr><tr><td>α<sub>0</sub></td><td>Debris decay</td><td>10<sup>-1</sup></td><td>Of the same order of α<sub>i</sub>. May be either positive or negative due to debris decaying from higher orbits.</td></tr><tr><td>β<sub>i</sub></td><td>Natural breakdown</td><td>10<sup>-1</sup></td><td>A function of the technology involved. Arbitrarily estimated as a fraction of α<sub>i</sub>.</td></tr><tr><td>a<sub>i</sub></td><td>Eqn (4)</td><td>10<sup>-1</sup>.</td><td></td></tr><tr><td>y<sub>0</sub></td><td>Number of debris</td><td>10<sup>3</sup></td><td>Dangerous debris at the altitude considered. Kept constant via ADR per an international agreement.</td></tr><tr><td>ρ</td><td>Discount</td><td>.9</td><td></td></tr><tr><td>φ<sub>i</sub></td><td>Satellite revenue</td><td>1</td><td>Following [3].</td></tr><tr><td>ψ<sub>i</sub></td><td>Saturation coefficient</td><td>10<sup>-3</sup> – 10<sup>-4</sup></td><td>If adding a satellite to the constellation i beyond the N-th brings no new value, then ψ<sub>i</sub> = φ<sub>i</sub>/N.</td></tr><tr><td>μ<sub>i</sub></td><td>Half cost maneuver</td><td>&lt; 5</td><td>It is the cost of a maneuver times 1/P<sub>0</sub>. by expanding the orbit keeping propellant.</td></tr><tr><td>z<sub>ij</sub></td><td>Relative collision frequency</td><td>.1 – 10</td><td>We choose z<sub>00</sub> = 1 by definition, and expect that the z<sub>ij</sub> are never more than a few units. (Operators would not use exceedingly crowded orbits.)</td></tr><tr><td>c<sub>i</sub></td><td>Unit launch cost</td><td>2 – 4</td><td>Economic viability of the constellation requires that c<sub>i</sub> &lt; ρφ<sub>i</sub>/(1 – ρa<sub>i</sub>). Another consideration is the time of return on investment in years.</td></tr><tr><td>η<sub>i</sub></td><td>Eqn (6)</td><td>10<sup>-2</sup></td><td>Consequence of above estimations.</td></tr><tr><td>ν</td><td># debris per collision</td><td>&gt; 50</td><td>Is the least defined of our parameters. We take it as a fixed number for lack of a better model. At least 50, possibly several thousands.</td></tr><tr><td>ω</td><td>Eqn (5)</td><td>10</td><td>Or a few tens. Consequence of other estimations.</td></tr></tbody></table>
# 6 Appendix B: Proofs

## 6.1 Proof of proposition 1

The criterion to be optimized is $\Pi = \sum_i \Pi_i$, i.e.:

$$ \Pi = \sum_{t=0}^{\infty} \rho^t \left[ -\frac{1}{2} \langle x, \bar{Q}x \rangle - \langle c, s \rangle + \langle b, x \rangle - p\omega \right]. $$

We make the informed guess that Bellman's return function will be of the form

$$ V^O(x) = -\frac{1}{2} \langle x, \bar{Q}x \rangle + \langle \ell, x \rangle + m. $$

The fact that we succeed in solving Bellman's equation with such a function will prove our guess right. It reads

$$ V^O(x) = \max_{s \ge 0} \left\{ -\frac{1}{2} \langle x, Qx \rangle - \langle c, s \rangle + \langle b, x \rangle - p\omega \right. \\ \left. + \rho \left[ -\frac{1}{2} \langle Ax + s, \bar{Q}(Ax + s) \rangle + \langle \ell, Ax + s \rangle + m \right] \right\}. $$

The r.h.s. is strictly concave in $s$. Differentiating with respect to $s$ and equating to zero, we obtain

$$ s^O = -Ax + \frac{1}{\rho} \bar{Q}^{-1} (\rho \ell - c). $$

Substituting in $V^O(x)$ is made easy by the fact that $Ax + s$ is simple. We find that the square term in $x$ is as “guessed”, and identifying other coefficients, we find

$$ \ell = b + A^t c \implies \rho \ell - c = \rho b - (I - \rho A^t)c = f, $$

and thus

$$ s^O = -Ax + \xi^O \quad \text{with} \quad \xi^O = \frac{1}{\rho} \bar{Q}^{-1} f, $$

and

$$ m = \frac{1}{1 - \rho} \left[ \frac{1}{2\rho} \langle f, \bar{Q}^{-1} f \rangle - p\omega \right]. $$

These formulas coincide with the formulas expanded componentwise in the proposition.

## 6.2 Proof of proposition 3

We seek an Isaacs Value function of the form

$$ V_i(x) = -\frac{1}{2} \langle x, Q^i x \rangle + \langle \ell^i, x \rangle + m^i. $$
We write Isaacs' equation in terms of $b^i = b_i \mathbb{1}^i$ and $c^i = c_i \mathbb{1}^i$:

$$V_i(x) = \max_{s_i} \left\{ \langle b^i, x \rangle - \frac{1}{2} \langle x, Q^i x \rangle - \langle c^i, s \rangle - p \zeta_i \omega \\ + \rho \left[ -\frac{1}{2} \langle Ax + s, Q^i(Ax + s) \rangle + \langle \ell^i, Ax + s \rangle + m^i \right] \right\}.$$

The matrix $Q^i$ is not positive definite, but the right-hand side above is strictly concave in $s_i$, with a quadratic term $-(1/2)\psi_i s_i^2$. We will find the maximum over $\mathbb{R}$ by differentiating. Let $Q_i^i$ stand for the line $i$ of $Q^i$. Note that the derivative of $\langle s, Q^i s \rangle$ with respect to $s_i$ is $2Q_i^i s$, and that of $\langle Q^i A x, s \rangle$ is $Q_i^i A x$. Differentiating all $V_i^N(x)$ with respect to $s_i$ each and equating to zero, we obtain

$$\frac{\partial V_i}{s_i} = -c_i^i - \rho Q_i^i s - \rho Q_i^i A x + \rho \ell_i^i = 0.$$

Let $\tilde{Q}$ denote the matrix whose line $i$ is $Q_i^i$. Regrouping these equations in a vector equation yields

$$\rho \tilde{Q} s = -\rho \tilde{Q} A x + \rho \ell^i - c^i.$$

Hence, with $\rho \ell^i - c^i = \tilde{f}^i$, we have

$$s = s^N = -Ax + \frac{1}{\rho} \tilde{Q}^{-1} (\rho \ell^i - c^i) = -Ax + \frac{1}{\rho} \tilde{Q}^{-1} \tilde{f}^i.$$

Substitute in $V_i$, we get

$$-\frac{1}{2}\langle x, Q^i x \rangle + \langle \ell^i, x \rangle + m^i = \langle b^i, x \rangle - \frac{1}{2}\langle x, Q^i x \rangle - \langle c^i, -Ax + \frac{1}{\rho}\tilde{Q}^{-1}\tilde{f}^i \rangle - p\zeta_i\omega \\ - \frac{1}{2\rho}\langle \tilde{Q}^{-1}\tilde{f}^i, Q^i\tilde{Q}^{-1}\tilde{f}^i \rangle + \langle \ell^i, \tilde{Q}^{-1}\tilde{f}^i \rangle + \rho m^i.$$

Identifying like powers of $x$, we get the right term in $x^2$, and

$$\ell^i = b^i + A^t c^i \implies \ell^i = (b_i + a_i c_i) \mathbb{1}^i.$$

It follows that $\tilde{f}^i = f_i \mathbb{1}^i = f^i$.

Finally, the terms without $x$ read

$$(1 - \rho)m^i = \frac{1}{\rho}\langle f^i, \tilde{Q}^{-1}f \rangle - \frac{1}{2\rho}\langle \tilde{Q}^{-1}f, Q^i\tilde{Q}^{-1}f \rangle - p\zeta_i\omega.$$

## 6.3 Proof of proposition 5

In the symmetric case, we have seen that $p_n^N < p^N$, which proves the first statement.

The dividend of cooperation in the general case is necessarily positive, the total profit under co-operation being the maximum of the sum of the profits of the players over all $n$-uples of admissible controls, it is larger than the same sum under the Nash equilibrium strategies. To get an order of magnitude of the difference, we investigate the *dividend of cooperation* DC:

$$\text{DC} = V^O(0) - \sum_{i=1}^{n} V_i^N(0).$$
We compare the joint optimum profit $V^O(0)$ and the sum of the individual profits in the Nash equilibrium. Let

$$V^O(0) = \frac{1}{2\rho(1-\rho)}W^O - 2\rho\omega p,$$
$$\sum_i V_i^N(0) = \frac{1}{2\rho(1-\rho)}W^N - 2\rho\omega p,$$

with

$$W^O = \langle f, \bar{Q}^{-1}f \rangle,$$

and, remembering that $\sum_i f^i = f$ and $\sum_i Q^i = \bar{Q}$, we obtain

$$W^N = 2\langle f, \tilde{Q}^{-1}f \rangle - \langle \tilde{Q}^{-1}f, \bar{Q}\tilde{Q}^{-1}f \rangle = 2\langle \tilde{Q}^{-1}f, \tilde{Q}\tilde{Q}^{-1}f \rangle - \langle \tilde{Q}^{-1}f, \bar{Q}\tilde{Q}^{-1}f \rangle.$$

Using the fact that a quadratic form only depends on the symmetric part of the coefficient matrix leads to

$$W^N = \langle \tilde{Q}^{-1}f, (\tilde{Q} + \tilde{Q}^t)\tilde{Q}^{-1}f \rangle - \langle \tilde{Q}^{-1}f, \bar{Q}\tilde{Q}^{-1}f \rangle = \langle \tilde{Q}^{-1}f, (\tilde{Q} + \tilde{Q}^t - \bar{Q})\tilde{Q}^{-1}f \rangle.$$

Note then that $\tilde{Q} + \tilde{Q}^t - \bar{Q} = \Psi$. Hence

$$W^N = \langle \tilde{Q}^{-1}f, \Psi\tilde{Q}^{-1}f \rangle$$

We expand both $W^O$ and $W^N$ to first order:

$$W^O = \langle f, \Psi^{-1}f \rangle - \tau \langle \Psi^{-1}f, \bar{U}\Psi^{-1}f \rangle + O(\tau^2),$$

and

$$W^N = \langle \Psi^{-1}f - \tau\Psi^{-1}\tilde{U}\Psi^{-1}f, \Psi[\Psi^{-1}f - \tau\Psi^{-1}\tilde{U}\Psi^{-1}f] \rangle + O(\tau^2).$$

Keeping terms of first order in $\tau$, and noticing that $\tilde{U} + \tilde{U}^t = \bar{U}$, we get

$$\begin{aligned}
W^N &= \langle f, \Psi^{-1}f \rangle - 2\tau \langle f, \Psi^{-1}\tilde{U}\Psi^{-1}f \rangle + O(\tau^2) \\
&= \langle f, \Psi^{-1}f \rangle - \tau \langle f, \Psi^{-1}(\tilde{U} + \tilde{U}^t)\Psi^{-1}f \rangle + O(\tau^2), \\
&= \langle f, \Psi^{-1}f \rangle - \tau \langle \Psi^{-1}f, \bar{U}\Psi^{-1}f \rangle + O(\tau^2).
\end{aligned}$$

Therefore, $W^O = W^N + O(\tau^2)$, i.e., the benefit of cooperation, which is by definition positive, is of the order of $\tau^2$, and therefore very small.

## References

[1] Académie des Sciences Large satellites constellations: challenges and impact. https://academie-sciences.hal.science/hal-04626677v1, 1–61, 2024.

[2] Adilov, N., Alexander, P., and Cunningham, B. (2018). An economic Kessler Syndrome: a dynamic model of earth orbit debris. Economic Letters, 166:79–82.
[3] Béal, S., Deschamps, M. and Moulin, H. Taxing congestion of the space commons. *Acta Astronautica*, 177:313-319, 2020.

[4] Bernhard, P., Deschamps, M. and Zaccour, G. Large satellite constellations and space debris: exploratory analysis of strategic management of the space commons. *European Journal of Operational Research*, 304:1140-1157, 2023.

[5] Bongers, A., Molinari, B., Rouillon, S. and Torres, J. The foundations of the economic of the outer space: a premier overview. *Space Economics Working Papers*, 2024.

[6] GAO (US Government Accountability Office) Large satellites constellations. Mitigating environmental and other effects. https://www.gao.gov/products/gao-22-105166, 1-80, 2022.

[7] Guyot, J. and Rouillon, S. Sustainable management of space activity in low Earth orbit. Journal of Environmental Economics and Policy, 13:188-212, 2023.

[8] IADC Space debris mitigation guidelines. *IADC Steering Group and Working Group 4*, 1-16, 2025a.

[9] IADC Report on the status of the space debris environment. *IADC Working Group 2/Steering Group*, 1-27, 2025b.

[10] Kessler, D. Collision cascading: the limits of population growth in low earth orbit. *Advances in Space Research*, 11-12:63-66, 1991.

[11] Kessler, D. and Cour-Palais, B. Collision frequency of artificial satellites: the creation of a debris belt. *Journal of Geophysical Research*, 83(2):2637-2646, 1978.

[12] Le May, S., Gehly, S., Carter, B. and Flegel, S. Space debris collision probability analysis for proposed global broadband constellations. *Acta Astronautica*, 151:445-455, 2018.

[13] Liou, J-C. U.S space debris environment and activity updates. *62nd Session of the Scientific and Technical Subcommittee, Committee on the Peaceful Uses of Outer Space, United Nations*, 1-11, 2025.

[14] Oltrogge, D. and Christensen, I. Space governance in the new space era. *Journal of Space Safety Engineering*, 7:432-438, 2020.

[15] O'Neill, G. The high frontier: human colonies in space. Space Studies Institute Press, 1977.

[16] Rouillon, S. A physico-economic model of low Earth orbit management. *Environmental and Resource Economics*, 77:695-723, 2020.

[17] Sandler, T. and Shulze, W. The economics of outer space. *Natural Resources Journal*, 21:371-393, 1981.

[18] Snow, M. Investment cost minimization for communication satellite capacity: refinement and application of the Chenery-Manne-Srinivasan model. *Bell Journal of Economics*, 6:621-643, 1975.
[19] Weinzierl, M. Space, the final economic frontier. Journal of Economic Perspectives, 32:173–192, 2018.

[20] Wihlborg, G. and Wijkman, P. Outer space resources in efficient and equitable use: new frontiers for old principles. Journal of Law and Economics, 24:23–43, 1981.