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

HAL
open science

Large satellite constellations and space debris:
exploratory analysis of strategic management of the
space commons

Pierre Bernhard, Marc Deschamps, Georges Zaccour

► To cite this version:

Pierre Bernhard, Marc Deschamps, Georges Zaccour. Large satellite constellations and space debris: exploratory analysis of strategic management of the space commons. European Journal of Operational Research, 2023, 304 (3), pp.1140-1157. 10.1016/j.ejor.2022.04.030. hal-03909662

HAL Id: hal-03909662

https://inria.hal.science/hal-03909662

Submitted on 29 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 Attribution 4.0 International License
Large satellite constellations and space debris: exploratory
analysis of strategic management of the space commons*

Pierre Bernhard<sup>a,*</sup>, Marc Deschamps<sup>b,1</sup>, Georges Zaccour<sup>c</sup>

**Abstract**

The use of space through satellites is more and more important for nations, companies, and individuals. However, since the first satellite was sent up in 1957, mankind has been polluting space with debris (i.e., artificial objects with no function), especially in low orbits (between 100 and 2000 km). The current situation is such that: 1/ space agencies send on average several collision risk alerts every day, and 2/ satellites as well as the International Space Station regularly perform avoidance maneuvers to escape being damaged or simply destroyed.

In addition, in the last few years, these problems have become more worrisome and
may permanently change dimension with the advent of mega-constellations of satel-
lites. Indeed, in order to develop telecommunications and high-speed Internet, several
companies (e.g., Starlink, Kuiper, OneWeb, Hongyan, Hongyun, Leosat, Athena) are
planning to send several tens of thousands of satellites into low orbits, which are al-
ready the most polluted.

The purpose of this paper is to provide an economic analysis in terms of dynamic
games of the trade-off between constellation size and cost of preserving the space en-
vironment. Our goal is to contribute to provide a framework for a sustainable develop-
ment of a space economy.

Keywords: Game theory, Satellites constellations, Space congestion, Active Debris Removal, Dynamic games.

<footnote><sup>a</sup>BIOCORE team, Centre INRIA de l'Université de la Côte d'Azur, B.P. 93, F-06902 Sophia Antipolis Cedex, France.</footnote>

<footnote><sup>b</sup>CRESE EA3190, Univ. Bourgogne Franche-Comté F-25000 Besançon, France and BETA CNRS, GREDEG CNRS and OFCE Sciences Po.</footnote>

<footnote><sup>c</sup>GERAD, HEC Montréal, Canada</footnote>

<footnote>*The authors would like to thank the three anonymous Reviewers for their very helpful comments.</footnote>

<footnote><sup>*</sup>Corresponding author</footnote>

<footnote>Email addresses: Pierre.Bernhard@inria.fr (Pierre Bernhard),
Marc.Deschamps@univ-fcomte.fr (Marc Deschamps), Georges.Zaccour@gerad.ca (Georges
Zaccour)</footnote>

<footnote n="1"><sup>1</sup>Marc Deschamps wishes to acknowledge usefull discussions with S. Béal, Ch. Bonnal, M. Diss, H. Moulin and L. Thomas. He was partially supported by grandt "Chrysalide UFC" (Université de Franche Comté).</footnote>
# 1. Introduction

The successful launch of the satellite Sputnik 1 on October 4th, 1957 inaugurated the era of mankind exploitation of the outer space, defined as the domain above the Karman line usually set at 100 km above the Earth (McDowel (2018)).<fnref n="2" /> Since then, the use of space for commercial and military purposes has been steadily growing, and this use is expected to accelerate with the arrival of mega constellations operated by private firms (e.g., Starlink, Kuiper, OneWeb, Hongyan, Hongyun, Leosat, Athena). In 2020, the space economy reached an estimated value of 371 billion dollars, of which 271 for the satellite industry alone (Bryce (2021)).

Satellites play a crucial role in our daily life, and we can no longer do without their services (OECD (2019)). To illustrate, we recall that: (i) Weather forecast entirely depends on measurements and pictures from space, and 26 of the 50 climate variables essential to assess the state of the Earth, are only fully accessible from space; (ii) Satellite data are systematically used to manage major disasters (earthquakes, volcanic eruptions, tsunamis, hurricanes, fires, etc.); and (iii) Navigation satellites such as GPS, Glonass, or Galileo are used to help us travel. On the negative side, the conquest of space has been accompanied by a production of space debris (essentially dead satellites, rocket upper stages, operational and fragmentation debris), which can potentially render space economically useless if their number exceeds a certain threshold.

The objective of this paper is to attempt to reconcile the growth in space exploitation with the preservation of the global commons, that is, to achieve a sustainable development of space.

## 1.1. *Space debris: Some background*

Sputnik 1, a small silver ball of 60 cm in diameter and 84 kg in mass, split into two pieces of debris, i.e., the last stage of the rocket that placed the satellite into orbit (a metal cylinder of 28 m in length, 3 m in diameter and a mass of 6500 kg) and a conical cap (of a mass of approximately 100 kg). Although the numbers may differ from one source to another, there is a broad consensus on the orders of magnitude of space debris. For instance, the former head of space debris activities at ESA (European Space Agency) estimated in 2015 that there were 5,400 objects larger than 1 m, 29,200 between 1 m and 10 cm,<fnref n="3" /> (see Figure 1) 740,000 between 10 cm and 1 cm, 170 million between 1 cm and 1 mm, and 360 billion between 1 mm and 0.1 mm (Klinkrad (2015)). Moreover, these objects are not homogeneously distributed throughout space, and the most densely populated area is in low orbits (between 100 and 2000 km). This region of space, used in particular by observation and meteorological satellites, is not itself uniformly affected by this pollution (see Figure 2).

Even though space is very large, the presence of these space debris and their concentration in certain areas, represent some risks for the exploitation of space due to

<footnote n="2">²The UNOOSA (United Nations Office of Outer Space Affairs) proposes a catalog of all the objects sent into space. Available at https://www.unoosa.org/oosa/osoiindex/index.jspx?lf_id. See AstriaGraph (http://astria.tacc.utexas.edu/AstriaGraph/) for a graphical view</footnote>

<footnote n="3">³The catalog of the objects of more than 10 cm is incomplete because some of them are difficult to detect, and because it excludes the military objects.</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Number of catalogued objects (Liou (2020))</figcaption></figure>

the possible accidental collisions between satellites and debris or collisions between debris. Currently, the probability of loss of a satellite due to debris is estimated at between 3 and 5% over its operational life, while the risk of loss of a satellite during its launch is 1 to 2%, and the risk of loss due to poor satellite design is around 1%. The risk of loss due to a collision with a piece of debris is related to the mass but especially to the size and the hypervelocity of the debris. Indeed, a piece of space debris located at 1,000 km of altitude will remain in orbit for a millennium while moving at 26 500 km/h. Moreover, the bigger it is, the higher the probability of collision with another space object. A piece of debris smaller than 1 mm colliding with a satellite generates essentially a wear of materials, while a piece of debris larger than 1 cm can lead to a fragmentation (with a probable deterioration of its functions) or to a total destruction of a satellite.<fnref n="4" />

To date four accidental collisions happened in low orbit between a satellite and a piece of debris, namely, the collision involving: (i) Cosmos 1934 satellite and a piece of catalogued debris (December 23, 1991);(ii) Cerise satellite and a piece of catalogued debris (July 24, 1996); (iii) Thor Burner 2A satellite and a piece of catalogued debris (July 17, 2005); and (iv) Iridium 33 satellite and the former Cosmos 2251 dead satellite (which had become a piece of debris) (Feb. 10, 2009).<fnref n="5" /> In addition, it should be

<footnote n="4">⁴As an indication, a 1 mm radius debris has a kinetic energy equivalent to a bowling ball thrown at 100 km/h, while a 1 cm radius debris has a kinetic energy equivalent to a big car running at 130 km/h.</footnote>

<footnote n="5">⁵The jump observed in 2009 in Figure 1 corresponds to this accidental collision, while that of 2007 is the deliberate result of an anti-satellite missile firing by the Chinese military on January 11 on one of its dead satellites (Feng-Yun 1-C). On March 27, 2019, the Indian government also carried out an anti-satellite firing on one of its satellites.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: Density in low orbit (Pardini and Anselmo (2021))</figcaption></figure>

mentioned that, as early as 1983, the windshield of the space shuttle Challenger was damaged by a collision with a piece of orbital debris. And on November 15, 2021, an antisatellite test by Russia created at least one thousand debris at the orbital altitude of the ISS, forcing the crew to rush for shelter in their reentry space vehicles.

In a recent paper Foreman et al. (2017) investigate the question of large satellite constellations on the orbital debris environment and use, as case studies, OneWeb and SpaceX. They conclude that: “even the introduction of a single constellation could permanently contaminate the LEO environment” and underline the risk of Kessler syndrome.

According to Pardini and Anselmo (2014), the risk of collision in low earth orbit was multiplied by 4.5 between 1980 and 2010, in a way almost parallel to the increase in the number of debris. To this, we must add that the increase in the intensity of debris on certain orbits could one day lead to a diverging chain reaction on these orbits, whereby the collisions produce more new debris than natural decay eliminates, eventually making these orbits physically unusable: this is the Kessler syndrome (Kessler and Cour-Palais (1978); Kessler (1991)).

The space debris congestion in low earth orbit is likely to undergo a major change of scale with the planned projects<fnref n="6" /> for mega-constellations of telecommunication satellites<fnref n="7" /> placed in low orbits to reduce their latency (i.e., their response time). These projects, with some of them already being implemented, consist in launching several

<footnote n="6">⁶For an overview of all these projects see https://www.newspace.im/.</footnote>

<footnote n="7">⁷See in particular Le May et al. (2018) for an assessment of collision probabilities concerning these mega-constellations.</footnote>
tens of thousands of satellites (e.g., SpaceX alone plans to send 42,000 satellites for its Starlink constellation) in order to allow the development of broadband Internet in every part of the planet. To realize the magnitude of these projects, we note that there are currently only about 2,100 active satellites in all orbits (Undseth et al. (2020)), and the congestion problem is already here. For instance, ESA was forced on September 2, 2019, to maneuver its Aeolus satellite to avoid a collision with a satellite of the SpaceX Starlink constellation. To make things worse, the detailed information about the shielding of the OneWeb and Starlink satellites is not public.

## 1.2. Managing space debris

Considered as a critical problem for low-Earth orbits, space agencies started to deal with space debris more than thirty years ago. Their approach, shared with all the stakeholders of the space sector, is articulated along three complementary axes. The first axis aims at avoiding collisions by improving the surveillance systems (for instance, USSPACECOM operates a network of 29 observation means<fnref n="8" />) and by carrying out avoidance maneuvers (on average, there is one maneuver per year and per satellite, and the international space station carries out about two maneuvers per month) (Bonnal et al. (2020)). The second axis seeks to reduce the production of debris, by (i) protecting the satellites from the effects of collisions (in particular by shielding<fnref n="9" /> and a reflection on the architecture and materials used in the satellites), and (ii) designing space operations so that they do not produce new debris (e.g., prohibition of explosions in space, controlled atmospheric re-entry of the satellite at the end of its life, or placing it in a graveyard orbit).<fnref n="10" /> Finally, the third axis is active debris removal, i.e., the collection and elimination of space debris. This is the most exploratory axis, with a set of techniques (non-contact, with rigid contact, or with soft contact) being still under study (Priyant and Kamath (2019)). The first missions (RemoveDebris, Clean Space-1, AnDROiD) are underway.

Whatever the approach followed to mitigate debris pollution; one unavoidable question is how to finance such effort.

## 1.3. Literature Review

In a review of operations research (OR) in the space industry, Fliege et al. (2012) regrouped the contributions under the following headings: (i) Parameter identification; (ii) Trajectory optimization; (iii) Motion and path planning; (iv) System design and system management; and (v) Satellite optimization. As one can easily imagine, a large literature in engineering, optimization, and control, covers each of these streams, and all will eventually contribute to the above outlined directions for reducing space congestion with debris. Reviewing the literature, even only after Fliege et al. (2012), is

<footnote n="8">⁸In 2021 a consortium of European countries has created the European Union Space Surveillance and Tracking support framework (https://www.eusst.eu/)</footnote>

<footnote n="9">⁹This leads to additional costs: an additional kilogram at launch costs between 15,000 and 25,000 € depending on the orbit of the satellite concerned.</footnote>

<footnote n="10">¹⁰See in particular ISO 24113 and associated standards, as well as the French law of April 9, 2008 on space operations.</footnote>
well beyond the objective of this paper and would take the readers along paths very re-
mote from our paper's focus. To give a hint about some of the optimization techniques
used in the literature, we give few examples mainly in system design and system man-
agement, one of the favorite topics in OR.

The design of an optimal production mix of launch vehicles by an aerospace com-
pany has been considered in Morgan et al. (2006). The authors formulate the prob-
lem as a bin packing one and proposed a heuristic to obtain efficiently the sought
mix. The optimal management of a satellite system has been studied in a series of
papers. In Chen et al. (2019), a multi-satellite scheduling problem is formulated as
a mixed-integer linear-programming (MILP) model, and a detailed analysis approach
of possible conflicts between constraints is developed, which allowed to significantly
speed up the solution process. Tangpattanakul et al. (2015) model the selection and
scheduling of observations for an agile Earth observing satellite as a multi-objective
optimization problem and developed an indicator-based multi-objective local search
(IBMOLS) to solve it. Typically, the users of satellites made requests for photographs
of areas of interest, in a time window, but not all of them can be satisfied. Bianchessi
et al. (2007) consider the problem of selecting the requests to be fulfilled by several
satellites performing multiple orbits over a given planning horizon, and solve it using
some heuristics coupled with a column generation algorithm to accelerate the search.
Nagendra et al. (2020) show how satellite big data analytics helped achieving better
coordination and collaboration between rescue teams for humanitarian relief efforts in
the case of the 2018 Kerala floods.

In parallel to the OR and engineering literature, a stream of papers focused on the
economics of space exploitation, with the first paper seems to be Sandler and Schulze
(1981). Weinzierl (2018) provides a general perspective on space economics, and
Oltrogge and Christensen (2020) gives an overview of the governance system. It is only
very recently that the problem of space debris attracted the attention of researchers, de-
spite the fact that the Kessler syndrome concept has been around for a while. To the best
of our knowledge, Adilov et al. (2015) are the first to theoretically analyze space debris
as an externality. They propose a model designed to couple the incentives of satellite
launching operators and the cleanup of space debris. Macauley (2015) provides an em-
pirical estimation of different policies that could be implemented to mitigate the effect
of space debris. These policies are regrouped into three categories, namely, launch
taxes, technology standards, and penalties for debris generation. The empirical formu-
las representing these policies are applied in a two-period model to assess the short-
and long-term costs and benefits of each of them. In the same vein, Grzelka and Wag-
ner (2019) offer an economic analysis of a series of policies that combine some ex-ante
launching incentives, e.g., increasing satellite quality, and ex-post take back debris
measures, similar to environmental cleaning in terrestrial context. (Also, see Muller
et al. (2011) for an economic evaluation of space debris removal.) Grzelka and Wag-
ner (2019) derive some specific recommendations for the US, taken as a case study.
Rouillon (2020) compares different approaches, i.e., maximum carrying capacity, open
access and optimal policy, to manage the low Earth orbit. These ways are typically used
in the analysis of (terrestrial) resource problems (e.g., high sea fisheries). The author
derives some insights based on realistic numerical illustrations. Rao et al. (2020) recall
that space pollution is nothing, but another tragedy of the commons triggered by agents
who do not fully internalize the damage cost imposed on all when making their satellite launching decisions. The authors argue that an internationally harmonized orbital-use fee can correct this negative externality and even increase the value of space industry.

The optimal decisions derived in the above contributions assume a competitive space sector, while Klima et al. (2016) and Béal et al. (2020) assume imperfect competition. Klima et al. (2016) propose an empirically calibrated game model to discuss why space users do not contribute enough to debris cleanup. The reason is similar to the one repeatedly obtained in other environmental contexts, namely, free riding is the best individual response, with however the long-term effect of destroying the commons. Note that Adilov et al. (2018) show that it is possible that the physical impossibility to use an orbit can be preceded by an economic impossibility, i.e., an economic Kessler syndrome can predate the physical Kessler syndrome. By launching satellites, space-faring agencies create space congestion, and in turn space pollution (debris). Béal et al. (2020) compare a tax applied to each new launch under two modes of play, i.e., noncooperative and cooperative (joint optimization) games. The authors obtain that under centralization, it is twice as easy to finance the cleanup cost than in the noncooperative equilibrium. The qualitative result that centralization is more efficient than decentralization is true by definition. The large size of the improvement brought by centralization provides a strong incentive to the actors to act in this direction.

Some of the above cited papers and ours have some connections to the literature dealing with the tragedy of the commons. (See Hardin (1968) for the seminal paper and Frischmann et al. (2019) for an overview.) Indeed, similar to high seas and the environment, the outer space does not belong to any juridiction and no firm can be a priori excluded from exploiting space. However, there are some differences that are worth mentioning. A first particularity is that the whole setup is dynamic, with elapsing time being a critical ingredient. The “resource” is free space, but the way it is exploited is quite special. Indeed, placing satellites into orbit leads to a time stream of revenues, and there is a specific cost for joining the commons, before exploiting it. These features have different types of implications on players’ payoffs then, say, on fishers’ outcomes. Further, when a satellite leaves the commons, the resource is restored completely and instantly. As a consequence, each player has an incentive to de-orbit its own satellites when they get old to leave room for newer ones.

## 1.4. *Research questions and contribution*

The developments of mega constellations by few companies require an analysis of space congestion that accounts for the strategic interactions between these big players. To capture this feature we assume that the space industry is formed of two players whose main business is to launch and exploit satellites in low-Earth orbits. The number of players is restricted to two because it seems impossible to plan for more than two mega constellations operating at the same orbital altitude.

We aim at answering the following research questions:

1. Is it possible to reconcile continuing growth in space exploitation, while preserving the global commons?

2. How do noncooperative and joint optimization satellites launching strategies compare?
3. Is it feasible to set a tax that fully covers the cost of an active debris removal program, while allowing firms to be profitable?

In a nutshell, we show that sustainable exploitation of space is feasible, and provide the conditions under which a tax can fully finance an active debris removal (ADR) program. By definition a joint optimal solution leads to higher total payoff than a Nash (and any noncooperative) equilibrium. Surprisingly, we obtain that the difference is rather small.

To the best of our knowledge, our paper is the first analysis of space debris issues in a fully dynamic and strategic context. As for any other form of pollution emissions, it is the accumulation of space debris that causes the environmental damage. Consequently, only a full dynamic description of the evolution of space debris can shed a light on their potential harm, including on the occurrence of the Kessler syndrome. Moreover, we enlarge to a new area, space management, the literature applying dynamic games to environmental problems; see, e.g., De Frutos and Martín-Herrán (2019) and El Ouardighi et al. (2020) for recent contributions. We hope that scholars in operations research, applied mathematics and economics will push further the theory and analysis developed here.

The remainder of this paper is organized as follows. Section 2 introduces our model by defining the state dynamics as well as the optimization problems faced by firms managing satellite constellations. Sections 3 and 4 define and compare the case where the optimization is joint with the case where the game between the two constellations is noncooperative à la Nash, respectively under linear and then strictly concave revenue functions. Section 5 proposes two taxation schemes (one on created debris and the other on satellite launches) to ensure the neutrality of the presence of these mega-constellations with respect to the ex-ante situation with respect to space debris. Our last section concludes the article by proposing avenues of development and by exposing the limits of our analysis.

## 2. The model

We propose an infinite-horizon discrete-time dynamic model, with a one-year time step. We start by introducing the state dynamics that specify the evolution of the number of satellites and other objects in the space over time, and next the players’ optimization problems.

### 2.0. Notation summary

Our model has a large amount of notation, in order to be adaptable to many different contexts: LEO, GEO, sun-synchronous, more or less coplanar orbits, and also to several revenue structures and taxing schemes. Beyond the “native” parameters, we also introduce in the sequel compound notation to simplify the calculations and the resulting formulas.

To help the reader we summarize here the notation introduced so far, and to be introduced in the sequel.
**Native parameters**

$C_i(s_i)$ Cost of launching $s_i$ satellites the same year.

$c_i, d_i$ $C_i(s_i) = c_i s_i + d_i s_i^2 / 2$.

$h(t)$ Number of actively removed debris via ADR the year $t$.

$k_i$ Approximation of $\pi_i x_i$ in the linear-cost model.

$p$ Unit cost of active debris removal (ADR).

$s_i(t)$ Number of satellites launched by operator $i$ at year $t$.

$T_i(x_i)$ $= \eta_i x_i + \zeta_i \omega$ Number of active satellites removal charged to operator $i$.

$x_i(t)$ Number of satellites of operator $i$ aloft at a time $t$.

$y(t)$ Number of debris aloft at time $t$.

$y_0$ $y(0)$, the number of debris sought to be kept.

$\alpha_i$ Relative rate of deorbiting either by natural decay or active deorbiting.

$\alpha_0$ Relative rate of deorbiting including arrivals from higher orbits.

$\beta_i$ Relative rate of on-orbit satellite death.

$\gamma$ Tax coefficient in the launching taxing scheme.

$\zeta_i$ Share of the non-attributable costs $\omega$ born by player $i$. $\zeta_1 + \zeta_2 = 1$.

$\kappa$ Relative per-ADR tax margin for robustness.

$\nu$ Avarage number of debris caused by a collision.

$\Pi_i$ Operator $i$'s profit.

$\pi_i$ Half the cost of an evasive maneuver for operator $i$.

$\rho$ Discount factor.

$\zeta$ Number of satellites launched on the orbit considered by other operators.

$\sigma$ $\mathbb{E}\zeta$.

$\tau$ “Congestion parameter”:

$\tau z$ probability of hitting one among $z$ objects per unit time.

$\varphi_i$ Unit yearly value of a satellite aloft in the linear-cost model.

$\psi_i$ Coefficient of the term in $x_i^2$ in the strictly concave cost problem.

**Further notation**

$a_i = 1 - \alpha_i - \beta_i - \tau y_0$.

$b_i = \varphi_i - p\eta_i = \varphi_i - p\beta_i - \tau b(\nu - 1)y_0$.

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

$g_i = \rho\varphi_i - (1 - \rho a_i)c_i$.

N: Superscript N: Nash equilibrium values.

NS: Superscript NS: Nash equilibrium values in the symmetric problems.

O: Superscript O: Optimal values in the joint optimisation problems,

OS: Superscript OS: Optimal in the symmetric joint optimisation problems.

~: A tilde accent for the strictly concave revenue function case.

$R_i(x_1, x_2)$: Revenue function.

$V_i(x_1, x_2)$: Bellman or Isaacs Value function.

$\Delta = \psi_1\psi_2 - \tau^2\pi_1\pi_2$.

$\delta = \psi_1\psi_2 - \tau^2(\pi_1 + \pi_2)^2$.

$\gamma^\dagger, \gamma^\ddagger, \gamma^-, \gamma^+, \gamma^*$ locally used in subsection 5.2.

$\eta_i = \beta_i + \tau(\nu - 1)y_0$.

$\theta = \rho(\psi + 2\tau\pi)\omega/2$ except in subsubsection 5.2.2. (See equation (30).)

$\xi^N = (\psi_i f_j - \tau \pi_i f_j) / (\rho \Delta)$.
$$ \xi_i^O = [\psi_j f_i - \tau(\pi_1 + \pi_2) f_j]/(\rho\delta). $$

$$ \omega = -\alpha_0 y_0 + \tau(v - 2)y_0^2/2 + \sigma. \text{ Non attributable yearly debris creation.} $$

## 2.1. State dynamics

For clarity, we build our space dynamics in three steps. First, we model the effect of space congestion on the survival rate of satellites of a focal constellation, assuming away additional launches and other factors such as natural failures and active deorbiting of satellites at the end of their useful life. Next, we develop the full model while still considering a single constellation. Finally, we extend it to two interacting constellations.

Denote by $x(t)$ the number of satellites of the focal constellation aloft and operative at time $t$, and by $y(t)$ the number of other objects aloft at time $t$ at roughly the same altitude.

### 2.1.1. Evolution of number of active satellites

As in Béal et al. (2020), we assume that all satellites are of the same size, and are on essentially circular orbits. We adopt the model of probability of collisions of their section 5: the probability of a single satellite colliding with any of $z$ objects is taken to be $P = \tau z$ per time period, for a small factor $\tau$ depending on the relative geometries of the orbits considered.

To estimate the expected number of living satellites at the end of a period, let $P_k$ be the probability of exactly $k$ collisions occur in that period. Therefore, we have

$$ P_0 = (1 - \tau y)^x \approx 1 - \tau xy, $$

$$ P_1 = x\tau y(1 - \tau y)^{x-1} \approx \tau xy[1 - (x - 1)\tau y]. $$

Neglecting terms in $\tau^2$, the expectation of $x(t + 1)$, knowing $x(t) = x$, is then given by

$$ \mathbb{E}x(t + 1) = xP_0 + (x - 1)P_1 \approx x(1 - \tau xy) + (x - 1)\tau xy[1 - (x - 1)\tau y] \approx x(1 - \tau y). $$

For the sake of simplicity, we will be content with a deterministic model

$$ x(t + 1) = x(t)(1 - \tau y(t)) \quad (1) $$

allowing for a non integer $x$.

The above dynamics can be interpreted as a discrete-time description of the evolution of a species $x$ in the presence of a predator $y$. This model is commonly used in population dynamics, e.g., in fisheries. In fact, (1) is one of the two Lotka-Voltera equations, the other one describing the evolution of $y$. In the next subsection, the dynamics of both state variables $x$ and $y$ are fully specified.

### 2.1.2. Dynamics of one constellation

Satellites leave the constellation either by natural decay, or by natural failures at a rate $\beta$, or through collisions with debris or other objects at a rate $\tau y(t)$, or by active deorbiting at the end of their useful life, typically 20 to 25 years. For mathematical tractability, we do not keep track of the age of each object (there may be hundreds or
thousands aloft), and assume that both the vintage effect and the natural decay can be well-approximated by one parameter, namely, a deorbiting rate $\alpha$. The rates $\alpha, \beta$ and $\tau$ are given. A natural extension would be to let them be described by some stochastic process, which would, however, complicate considerably the analysis and render any effort to gain a qualitative insight into the problem of space congestion completely vain.

Denote by $s(t)$ the number of satellites that the decision maker chooses to add to the constellation during the period $[t, t+1]$. (The delay between the decision time and the actual launch is due to, e.g., necessary planning and fabrication.) Consequently, the dynamics of $x$ become

$$x(t+1) = (1 - \alpha - \beta)x(t) - \tau x(t)y(t) + s(t), \quad x(0) = 0. \quad (2)$$

Now, we turn to the dynamics of $y$. There are a number $y_0$ of other objects at the start of creation of the constellation. Being at the same altitude, they have the same rate of natural decay as constellation satellites, but decaying debris from higher orbits enter the orbital level considered. The balance will be taken as $\alpha_0 y(t)$, where $\alpha_0 \in (0, 1)$. Failing constellation satellites become debris. Collisions between two objects, satellites or debris, engender new debris, whose number depends on a series of factors, e.g., sizes of the two objects and their velocities. The determination of this number requires the development of a statistical model of breakup. Instead, we make the simplifying assumption that a collision between two objects generates $\nu > 2$ of new debris. Further, there are an average of $\tau xy$ collisions between satellites and debris and $\tau y^2/2$ collisions within the set of debris. Collisions between active satellites of the same constellation are ruled out by design.

A small number of lone satellites may be launched at the same altitude. The impact of their presence on $y$ can be represented by a random process $\zeta(t)$. However, to keep it simple, we assume that this effect can be well-approximated by the positive mean $\sigma$ of the process. Finally, at each time period, a number $h(t)$ of debris are removed via active debris removal (ADR). Consequently, we end up with the following dynamics:

$$y(t+1) = (1 - \alpha_0)y(t) + \beta x(t) + \tau \left[ x(t)(\nu - 1) + \frac{y(t)}{2}(\nu - 2) \right] y(t) + \sigma - h(t), \\ y(0) = y_0.$$

**Remark 1.** Incidentally, we also have a model of the so-called Kessler syndrome, when the debris are so numerous that their collisions generate more debris and the process diverges. In our model, this happens if

$$y > \frac{2}{\tau(\nu - 2)} \alpha_0.$$

### 2.1.3. Dynamics of two constellations

We extend our state dynamics to account for the presence of two constellations, controlled by two players (firms or consortia), at roughly the same altitude. Denote
with an index $i \in \{1, 2\}$ the same quantities as above without index, but now pertaining to player $i$. We will assume that collisions between satellites of the two constellations are avoided through just-in-time collision avoidance (JCA) maneuvers. (See below). The dynamics become

$$x_i(t + 1) = (1 - \alpha_i - \beta_i - \tau y(t)) x_i(t) + s_i(t), \quad (3)$$

$$y(t + 1) = (1 - \alpha_0)y(t) + \beta_1 x_1(t) + \beta_2 x_2(t) \\ + \tau \left[ (x_1(t) + x_2(t))(v-1) + \frac{y(t)}{2}(v-2) \right] y(t) + \sigma - h(t), \quad (4)$$

$$x_i(0) = 0, \quad i = 1, 2, \quad y(0) = y_0. \quad (5)$$

Up to now, we have not discussed how the number of debris removal is determined. We assume that, following an international agreement, it has been decided to keep the number of debris constant at $y_0$,<fnref n="11" /> and an international “ADR agency” is in charge of actively removing excess debris. This translates into having

$$h(t) = -\alpha_0 y_0 + \tau(v-2) \frac{y_0^2}{2} + \sigma + [\beta_1 + \tau(v-1)y_0] x_1(t) + [\beta_2 + \tau(v-1)y_0] x_2(t),$$

and

$$x_i(t + 1) = (1 - \alpha_i - \beta_i - \tau y_0) x_i(t) + s_i(t), \quad i = 1, 2.$$

To save on notation, let

$$a_i = 1 - \alpha_i - \beta_i - \tau y_0, \quad i = 1, 2, \quad (6)$$

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

$$\eta_i = \beta_i + \tau(v-1)y_0, \quad i = 1, 2. \quad (8)$$

we then have

$$x_i(t + 1) = a_i x_i(t) + s_i(t), \quad i = 1, 2, \quad (9)$$

and

$$h(t) = \omega + \eta_1 x_1(t) + \eta_2 x_2(t). \quad (10)$$

Clearly, the parameter values must be such that $a_i > 0, i = 1, 2$; otherwise, the dynamics of $x$ would not make much sense.

## 2.2. Optimization problems

We will develop two models of the optimization problems faced by the players, and for each model we will compare two possible strategic behaviors: cooperative optimum versus noncooperative Nash equilibrium. Each player's profit-maximization problem involves:

<footnote n="11">¹¹Thus avoiding the quandary of determining who to charge for old debris cleaning.</footnote>
* A revenue function that depends on the total number of operative satellites and is denoted by $R_i(x_1, x_2)$. It will result that

$$ \frac{\partial R_i}{\partial x_i} > 0, \quad \frac{\partial^2 R_i}{\partial x_i^2} \le 0, \quad \frac{\partial R_i}{\partial x_{3-i}} < 0, \quad \frac{\partial^2 R_i}{\partial x_i \partial x_{3-i}} \le 0. \quad (11) $$

The first two derivatives mean that revenues are concave increasing in own fleet of satellites. Concavity of revenue function is a standard assumption in economics. The last two derivatives imply that the two “products” are strategic substitutes. To illustrate, a Cournot duopoly model satisfies the properties in (11).

Such continuous revenue functions do not fit all types of constellations. Typically, a more realistic revenue structure for telecommunication constellations would display a discontinuity at low satellite numbers.<fnref n="12" /> Taking such a feature into account would require a completely different set of mathematical tools, and does not seem feasible in the context of the current article.

* A cost function of launching new satellites $C_i(s_i)$, assumed to be convex increasing, with $C_i(0) = 0$, that is, there is no fixed cost.

* A tax payment $T_i(x_i)$ as contribution to ADR cost. We want this payment to reflect the contribution of each player to these costs, i.e., its contribution to the creation of debris $\eta_i x_i$, plus a share $\zeta_i$ of the non directly attributable part $\omega$, with $\zeta_1 + \zeta_2 = 1$. Hence we let

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

How the $\zeta_i$ are chosen is a classical problem of public good economics, and will not be discussed at this point. The unit cost of ADR will be denoted by $p$.

*Just-in-case collision avoidance (JCA).* The precise definition of $R_i$ will depend on how we take into account the risk of collisions between satellites of the two constellations. According to our model, the mean frequency of such events is $\tau x_1 x_2$. Collisions are avoided by letting one of the two satellites involved maneuver. This has a cost in terms of expenditure of propellant, which, in turn, translates into time in orbit, as the same propellant is used for orbit keeping. We assume that each constellation is in charge of maneuvering half of the time, and let the cost of one such maneuver be $2\pi_i$. Hence the cost incurred on the average by player $i$ is $\tau \pi_i x_1 x_2$ per year.

**Remark 2.** One might object that each collision risk causes a maneuver and the mere risk is much more frequent than actual collisions in the absence of evasive maneuver. Say ten times more frequent if a 10% collision risk identified causes such an action. This is easily accounted for in our model by simply multiplying both coefficients $\pi_i$ by the same factor.

<footnote n="12">¹²We thank an anonymous Reviewer for this remark.</footnote>
Denoting by $\rho \in (0, 1)$ the common discount factor, Player i's profit is as follows:

$$ \Pi_i = \sum_{t=0}^{\infty} \rho^t [R_i(x_1(t), x_2(t)) - C_i(s_i(t)) - pT_i(x_i(t))] . \quad (12) $$

By (9-12), we have defined a two-player infinite-horizon discrete-time dynamic game, with two state variables ($x_1$ and $x_2$) and one control variable $s_i$ for each player.

### 3. Linear revenue functions

In this section, we keep $R_i$ linear to get as simple a model as possible, but still satisfying the conditions (11). We introduce for each player a positive coefficient $\varphi_i$ representing the economic value of one satellite in orbit for one year, and a positive coefficient $k_i$ approximating $\pi_i x_i$ in the JCA term, to obtain

$$ R_i(x_1, x_2) = \varphi_i x_i - \tau k_i x_{3-i}, \quad i = 1, 2. $$

Simultaneously, we assume that the $C_i$ are quadratic, defined by two positive coefficients $c_i$ and $d_i$, and of the form

$$ C_i(s_i) = c_i s_i + \frac{d_i}{2} s_i^2. $$

We chose a convex increasing function to account for the possible additional complexity and cost involved when there is multiple launching during the same year because of, e.g., the saturation of the launch facilities, and the increasing cost of personnel. (In the next section, we will consider the case where $d_i = 0, i = 1, 2$.)

Under these assumptions, the objective functional and the dynamics in (9)-(12) are linear in the state variables. For this class of linear-state dynamic games, it is well-known that the value functions are linear and that feedback and open-loop Nash equilibria coincide (see, e.g., Başar and Olsder (1999), Engwerda (2005) and Haurie et al. (2012)). Further, as the dynamic optimization problem in the coordinated case and the dynamic game problem in the noncooperative game are autonomous, i.e., they do not depend explicitly on time, we look for stationary solutions.

Unless it causes an ambiguity, we shall omit from now on the time argument.

To further save on notation, we will, from now on, let for $i \in \{1, 2\}$:

$$ b_i = \varphi_i - p\eta_i = \varphi_i - p\beta_i - \tau p(\nu - 1)y_0, $$

$$ g_i = \rho\varphi_i - (1 - \rho a_i)c_i, $$

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

Note that $b_i$ and $f_i$ are increasing with the marginal revenue coefficient $\varphi_i$, but decreasing or invariant with the failure rate coefficient $\beta_i$, the “price” $p$ at which ADR is charged to the players, and the number $y_0$ of debris.

Notice that sending a satellite at time $t = 0$ has a cost of $c_i$ or more (depending on whether $d_i$ is null or positive) and, even if the constellation is alone, brings from time $t = 1$ on a revenue at most $a_i^t \varphi_i$ per time step. Hence, the operation can only yield a profit under the following hypothesis:
**Hypothesis 1.** The coefficients $\rho$, $a_i$, and $c_i$ are, for $i \in \{1, 2\}$, such that

$$ \frac{\rho\varphi_i}{1 - \rho a_i} > c_i, $$

or, equivalently, $g_i > 0$.

### 3.1. Joint optimal solution

Suppose that the two players agree to coordinate their launching strategies to maximize their joint profit. Then, they solve a standard infinite-horizon discrete-time dynamic optimization problem. The control and state variables are superscripted with $O$ (for Optimal solution).

Let $V(x_1, x_2)$ be the value function, that is, the maximal total discounted future profit of the players over an infinite horizon when the number of satellites aloft is $x = (x_1, x_2)$.

**Proposition 1.** If $f_i - \rho\tau k_j \ge 0$, then the optimal number of yearly launches is given for $i \in \{1, 2\}$ and $j = 3 - i$ by

$$ s_i^O = \frac{f_i - \rho\tau k_j}{d_i(1 - \rho a_i)}, $$

and the steady-state value by

$$ x_{i\infty}^O = \frac{f_i - \rho\tau k_j}{d_i(1 - \rho a_i)(1 - a_i)} = \frac{s_i^O}{1 - a_i}. $$

The value function is as follows:

$$ V^O(x_1, x_2) = \sum_{i=1}^{2} \left\{ \frac{b_i - \tau k_j}{1 - \rho a_i} x_i + \frac{1}{1 - \rho} \left[ \frac{1}{2d_i} \left( \frac{f_i - \rho\tau k_j}{1 - \rho a_i} \right)^2 - \zeta_i p\omega \right] \right\}. $$

**Proof.** See Appendix.

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

We make the following remarks: First, the result that the optimal launching policy is constant and the value function is linear in the state is, as alluded to it before, a by-product of the linear-state structure of the dynamic optimization problem.

Under this strategy, each firm launches at every time step (year) the expected number of lost satellites when at full stationary constellation number. Because at the beginning, the number of satellites in the constellation is smaller, the number of satellites actually lost is smaller, and the excess in launches causes the actual number of satellites to grow toward its stationary value. (Reached in finite time in the real process because it happens in integer numbers.) This is, as expected, a “turnpike” behavior.

Second, the total discounted payoff that the players obtain over the infinite planning horizon is given by

$$ V^O(0, 0) = \frac{1}{1 - \rho} \sum_{i=1}^{2} \left[ \frac{1}{2d_i} \left( \frac{f_i - \rho\tau k_j}{1 - \rho a_i} \right)^2 - \zeta_i p\omega \right]. \quad (13) $$

It is easy to verify that this total payoff is decreasing in cost parameters, that is, in $c_i$, $d_i$, $k_i$ and $p$, and increasing in the revenue parameters $\varphi_i$. (Recall that $\zeta_1 + \zeta_2 = 1$.)
Third, the launching strategy is dictated by the familiar rule of marginal cost equal marginal benefit. To see it, we write the optimality condition as follows:

$$d_i s_i + c_i = \rho \frac{\partial V}{\partial x_i}, \quad i = 1, 2.$$

The left-hand side is the marginal cost and the right-hand side is the discounted marginal value (given by the derivative of the value function) of an additional launch.

Fourth, the proposition is stated under the assumption that $i \in \{1, 2\}$, $j = 3 - i$, $s_i^O \ge 0$, which is the case if

$$p \le \min_i \frac{\rho\varphi_i - (1 - \rho a_i)c_i - \rho\tau k_j}{\rho\eta_i} = \frac{g_i - \rho\tau k_j}{\rho\eta_i} \triangleq p^O, \quad i = 1, 2. \quad (14)$$

This inequality puts an upper bound on the tax rate that an agency can collect. To ensure that this upper bound is positive, we need to assume that, for $i \in \{1, 2\}$:

$$g_i - \rho\tau k_j > 0 \iff c_i < \frac{\rho(\varphi_i - \tau k_j)}{1 - \rho a_i}, \quad (15)$$

which is a strengthening of Hypothesis 1 to take into account the negative externality each player exerts on the other one.

Further, if $s_i^O \ge 0$, $i = 1, 2$, then the steady-state value $\hat{x}_i^O$ is also positive. It is decreasing in $c_i$, $d_i$, $k_j$, $\eta_i$ and $p$, and increasing in $\varphi_i$, which is quite intuitive. May be less intuitive is the fact that the optimal launching rate $s_i^O$ is decreasing with the satellite failure rate $\beta_i$.

Using (10), we conclude that the debris removal at any period $t$ is constant and given by

$$h^O(t) = h^O = \omega + \sum_{i=1}^{2} \eta_i \frac{f_i - \rho\tau k_j}{d_i(1 - \rho a_i)(1 - a_i)} = \omega + \sum_{i=1}^{2} \frac{\eta_i s_i^O}{(1 - a_i)}.$$

*Profitability.* Clearly, if the price $p$ of ADR is too high, the operation will not be profitable for the firms who will therefore exit this industry. We write the condition for profitability $V(0, 0) \ge 0$ in the symmetric case where $a_1 = a_2 = a$, $c_1 = c_2 = c$, $k_1 = k_2 = k$. It reads:

$$\mathcal{P}^O(p) \triangleq \rho^2 \eta^2 p^2 - 2 \left[ (g - \rho \tau k) \rho \eta + d(1 - \rho a)^2 \frac{\omega}{2} \right] p + (g - \rho \tau k)^2 \ge 0.$$

Notice that $\mathcal{P}^O(p^O) = -p\omega/2 < 0$ implies that the polynomial $\mathcal{P}^0$ has two real roots, the smaller one being less than $p^O$, while Hypothesis 1 and therefore $\mathcal{P}(0) > 0$ ensures that this smaller root is positive. We conclude that in the symmetric coordinated problem, there exists a pair of strategies that yields a joint profit if and only if condition (15) holds and further

$$p \le p_-^O \triangleq p^O - \frac{1}{\rho^2 \eta^2} \left[ \sqrt{[(g - \rho \tau k) \rho \eta + d(1 - \rho a)^2 \omega/2]^2 - \rho^2 \eta^2 (g - \rho \tau k)^2} - d(1 - \rho a)^2 \omega/2 \right].$$
### 3.2. Nash equilibrium

If the game is played noncooperatively, then the players seek a Nash equilibrium. Denote by $V_i^N(x_1, x_2)$ the value function of Player i. To save on notation, let, for $i \in \{1, 2\}$ and $j = 3 - i$:

$$m^i = \frac{1}{1 - \rho} \left[ \frac{f_i^2}{2d_i(1 - \rho a_i)^2} - \frac{\tau \rho k_i f_j}{d_j(1 - \rho a_j)^2} - p \zeta_i \omega \right].$$

**Proposition 2.** If

$$p \le p^N \triangleq \min_i \frac{g_i}{\rho \eta_i} = \min_i \frac{\rho \varphi_i - (1 - \rho a_i) c_i}{\rho \eta_i}, \quad (16)$$

the unique Nash equilibrium of satellite launches is given by

$$s_i^N = \frac{f_i}{d_i(1 - \rho a_i)}, \quad i = 1, 2, \quad (17)$$

and the steady-state values by

$$x_{i\infty}^N = \frac{f_i}{d_i(1 - \rho a_i)(1 - a_i)} = \frac{s_i^N}{1 - a_i}, \quad i = 1, 2.$$

The value functions are as follows, still for $i \in \{1, 2\}$, $j = 3 - i$:

$$V_i^N(x_1, x_2) = \frac{b_i}{1 - \rho a_i} x_i - \frac{\tau k_i}{1 - \rho a_j} x_j + m^i.$$

**Proof.** See Appendix.

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

Hypothesis 1 suffices to ensure that $p^N$ is positive.

Each player's total payoff is increasing in its own number of satellites aloft and decreasing in the rival's one. Although we do not have in this setup a demand law involving the two stocks $x_1$ and $x_2$, still each player suffers from the presence of the other through space congestion.

The results carry similar messages to the ones in the joint optimization case. Indeed, again we see that the launching strategy follows the basic rule of marginal cost equal marginal revenues. The condition for having strictly positive number of launches is that $p < p^N$ as given by (16).

Under the above condition, the steady state will also be positive. Further, it is easy to verify that Player i's strategy is decreasing in $c_i$, $d_i$, $\beta_i$ and $p$, and increasing in its marginal revenue parameter $\varphi_i$. In this noncooperative behavior, this strategy is independent of its competitor's parameters. Finally, the number of debris that must be removed per period is constant and given by

$$h^N(t) = h^N = \omega + \eta_1 \frac{s_1^N}{1 - a_1} + \eta_2 \frac{s_2^N}{1 - a_2}.$$

The total discounted payoff that the players obtain over the infinite planning horizon is given by $m_i$, i.e.,

$$V^N(0, 0) = \frac{1}{1 - \rho} \left[ \frac{f_i^2}{2d_i(1 - \rho a_i)^2} - \frac{\tau \rho k_i f_j}{d_j(1 - \rho a_j)^2} - p \zeta_i \omega \right]. \quad (18)$$
*Profitability.* We may investigate under which condition the operation will be profitable for the players, i.e., $V^N(0, 0) \ge 0$. As previously, we write this condition in the symmetric case:

$$ \mathcal{P}^N \triangleq \rho^2 \eta^2 p^2 - 2 \left[ (g - \rho \tau k) \rho \eta + d(1 - \rho a)^2 \frac{\omega}{2} \right] p + g(g - 2 \rho \tau k) \ge 0. $$

Again, $\mathcal{P}^N(p^N) = -2d(1-\rho a)^2\omega/2 < 0$. The implication is that in the Nash equilibrium symmetric case, the pair of strategies $s_i^N$ as given by equation (17) yield a net profit for the players if and only if, on the one hand

$$ g - 2\tau\rho k > 0, \qquad (19) $$

and on the other hand

$$ p \le p_{-}^N \triangleq p^O - \frac{1}{\rho^2 \eta^2} \left[ \sqrt{ \left[ (g - \rho \tau k) \rho \eta + d(1 - \rho a)^2 \omega / 2 \right]^2 - \rho^2 \eta^2 g(g - 2 \rho \tau k) - d(1 - \rho a)^2 \omega / 2 } \right]. $$

### 3.3. Comparison

We wish to compare the jointly optimal solution and the Nash equilibrium. This will be done assuming that both conditions (14) and (16) be met, but we notice that the first is always the most demanding. Assuming, thus, that $p \le p^O$, we have the following results:

**Proposition 3.** *The players launch more satellites in the Nash equilibrium than in the jointly optimal solution.*

**Proof.** Computing the differences between the strategies in the two modes of play, we get

$$ s_i^N - s_i^O = \frac{\rho \tau k_j}{d_i(1 - \rho a_i)} > 0. $$

Although the difference is positive, it seems to be negligible, because $\tau$ is a “small” parameter estimated to be between $10^{-6}$ and $10^{-7}$ (see Section 2.1.1). The fact that the players launch more satellites in the noncooperative equilibrium than in the coordinated solution is not surprising. In the jointly optimal solution, the decision maker in charge of managing both firms, internalizes the impact of a launch by either of them on both of them. (Note that in $s_i^O$ we have the term $-\tau k_i$, absent in $s_i^N$.) In the noncooperative solution, the “cross” impact (or the negative externality) is ignored, which leads to a larger number of launches. This over exploitation of the commons when the agents are selfish is a result that has been repeatedly obtained in the pollution control literature, as well as in the renewable resources (fisheries, forests) literature; see, e.g., the surveys in Jørgensen et al. (2010) and Long (2011).

Concerning the constraints to impose to have profitable solutions, we have:
**Proposition 4.** The constraints on the parameters of the problem and on the maximum p allowable to have profitable operations are both more restrictive for the Nash equilibrium than for the coordinated solution.

The total profit under cooperation is higher than the sum of individual noncooperative profits. The difference, or dividend of cooperation DC is given by

$$DC = V^O(0, 0) - V_1^N(0, 0) - V_2^N(0, 0) = \tau^2 \frac{\rho^2}{1-\rho} \sum_{i=1}^{2} \frac{k_j^2}{2d_i(1-\rho a_i)^2}. \quad (20)$$

**Proof.** Clearly, inequality (19) implies inequality (15). And we observe that $P^O(p) = P^N(p) + \rho^2\tau^2k^2$. Hence $P^N(p) \ge 0$ implies $P^O(p) > 0$. Equality (20) is obtained by a straightforward difference using equations (13) and (18). ■

This difference is indeed positive, but it is of the order of $\tau^2$, hence possibly very small. One should bear in mind however that $k_i$ may be rather large (several hundreds time $\varphi_i$ if the constellations are large) and that $1 - \rho a_i$ may be rather small.

So, clearly, the players should prefer coordination over noncooperation, if they can agree on a sharing of the total cooperative payoff that is mutually beneficial. But it is debatable, and depends on more refined numerical values, whether the quest for this additional profit is worth a huge effort.

What about the preference of the other party, namely, the space? In our framework, the quality of the space environment is the same in both mode of plays, because we constrain the system to reach the same target $y_0$. However, the economic implications vary with the mode of play as the cost associated with ADR is not the same. The difference in debris removals per period is given by

$$h^N - h^O = \tau \rho \sum_{i=1}^{2} \frac{\eta_i k_i}{d_i(1 - \rho a_i)(1 - a_i)} > 0.$$

Recalling that $\rho a_i < a_i < 1$ and $\eta_i = \beta_i + \tau(v-1)y_0$, clearly this difference is increasing in $a_i$, $k_i$ and $\eta_i$, and decreasing in $d_i$. However, $a_i$ being decreasing with $\beta_i$ while $\eta_i$ is increasing, the dependence on the failure rate $\beta_i$ is ambiguous.

Note that this difference is independent of the other cost parameter $c_i$, which cancels out. Finally, the larger the initial congestion $y_0$, the larger the environmental benefit of cooperation.

**Remark 3.** If the players agree to cooperate, then the question of how to divide the collective gain among them arises. One option would be to adopt the Nash bargaining solution (NBS) that gives to each firm its noncooperative payoff plus half of the dividend of cooperation. Using the NBS, Player i, i = 1, 2 ends up getting

$$NBS_i = m_i + \frac{\tau^2}{1-\rho} \sum_{i=1}^{2} \frac{k_{3-i}^2}{4d_i(1-\rho a_i)^2}.$$

**Remark 4.** To make a link between the tax policy given by

$$T_i(x_i) = \eta_i x_i + \zeta_i \omega.$$
and the model’s parameter, let us substitute for $\eta_i$ and consider the steady-state values of $x_i$ at Nash equilibrium and cooperation, respectively, to get

$$T_i^N = \frac{\beta_i + \tau(\nu - 1)y_0}{1 - a_i} s_i^N + \zeta_i \omega = \frac{(\beta_i + \tau(\nu - 1)y_0)f_i}{(1 - a_i)d_i(1 - \rho a_i)} + \zeta_i \omega, \quad \text{for } i = 1, 2,$$

$$T_i^C = \frac{\beta_i + \tau(\nu - 1)y_0}{1 - a_i} s_i^O + \zeta_i \omega = \frac{(\beta_i + \tau(\nu - 1)y_0)(f_i - \rho \tau k_j)}{(1 - a_i)d_i(1 - \rho a_i)} + \zeta_i \omega, \quad \text{for } i \neq j = 1, 2.$$

Clearly, the tax is increasing in the number of launches, which are behind space con- gestion. As a consequence of $s_i^N > s_i^O$, the tax is higher when the players behave noncooperatively than when they cooperate. Further, the tax is increasing in $y_0$, which is the taget level of ADR activities. Finally, when the two players cooperate, the tax paid by Player i involves $k_j$, a parameter that is part of the just-in-case collision avoid- ance cost. Once more, we see that cooperation internalizes each other’s cost, while it does not in a noncooperative mode of play.

## 4. Strictly concave revenue functions

In this section, we assume that the economic value of x satellites aloft is strictly concave, of the form $\varphi_i x_i^2 - (\psi_i/2)x_i^2$. Furthermore, we keep the nonlinear JCA term in its original form $\tau\pi_i x_1 x_2$, so that the overall revenue function becomes

$$R_i(x_1, x_2) = \varphi_i x_i - \frac{1}{2} \psi_i x_i^2 - \tau \pi_i x_1 x_2, \quad i = 1, 2, \qquad (21)$$

where $\varphi_i$, $\psi_i$ and $\pi_i$ are strictly positive. Clearly, this function satisfies the properties in (11). We can interpret it in classical terms of economic analysis. Let us rewrite it as

$$R_i(x_1, x_2) = x_i \left( \varphi_i - \frac{1}{2} \psi_i x_i - \tau \pi_i x_{3-i} \right), \quad i = 1, 2.$$

Then, the revenues as given by (21) are interpreted as production capacity, times a de- creasing function $D(x_1, x_2) = \varphi_i - \frac{1}{2}\psi_i x_i - \tau\pi_i x_{3-i}$, which plays the role of an inverse demand function. Such a revenue function is common in an imperfectly competitive market à la Cournot, assuming that the firms produce at their full capacity. The pa- rameter $\varphi_i$ is the maximum willingness to pay when $x_1 = x_2 = 0$. Alternatively, it can be interpreted as the marginal revenue at $(x_1, x_2) = (0, 0)$. The fact that $\tau\pi_i$ is positive means that the two goods (constellations) are strategic substitutes. In (21), the two con- stellations are imperfect substitutes. (If $\psi_1 = \psi_2 = \psi$, $\pi_1 = \pi_2 = \pi$ and $\psi = 2\tau\pi$, then the two products would be homogeneous.) In Adilov et al. (2018), the authors also introduce an inverse demand function in, however, a competitive market (open access).

Assuming that the launching cost is linear ($C_i(s_i) = c_i s_i$), Player i’s profit that it seeks to maximize, is then as follows:

$$\tilde{\Pi}_i = \sum_{t=0}^{\infty} \rho^t \left[ \varphi_i x_i(t) - \frac{1}{2} \psi_i x_i^2(t) - \tau \pi_i x_1(t) x_2(t) - c_i s_i(t) - p(\eta_i x_i(t) + \zeta_i \omega) \right]. \quad (22)$$
In the next two propositions, we state the results in the coordinated and nonco-operative case, respectively. As the results can be interpreted in the same way as the corresponding ones in the linear revenue function scenario, we shall avoid repeating (to some extent) the same arguments and focus on the new features that are induced by the strict concavity of the revenue functions.

### 4.1. Joint optimal solution

Let

$$ \delta = \psi_1\psi_2 - \tau^2(\pi_1 + \pi_2)^2, $$

and, for $i \in \{1, 2\}$ and $j = 3 - i$:

$$ \xi_i^O = \frac{1}{\rho\delta}[\psi_j f_i - \tau(\pi_1 + \pi_2)f_j]. $$

Recall that the superscript $O$ stands for optimal solution. To distinguish with the linear revenue function setup, we tilded the results.

**Proposition 5.** *Assuming the following expression is non-negative, the optimal policy is given by*

$$ \tilde{s}_i^O(x_i) = -a_i x_i + \xi_i^O, \quad i = 1, 2, $$

and steady-state values by

$$ \tilde{x}_{i\infty}^O = \xi_i^O, \quad i = 1, 2. $$

The value function is given by

$$ \tilde{V}^O(x_1, x_2) = \frac{1}{1-\rho} \left[ \frac{1}{2\rho\delta} (\psi_2 f_1^2 - 2\tau(\pi_1 + \pi_2)f_1 f_2 + \psi_1 f_2^2) - p\omega \right] \\ - \frac{1}{2} (\psi_1 x_1^2 + 2\tau(\pi_1 + \pi_2)x_1 x_2 + \psi_2 x_2^2) + (a_1 c_1 + b_1)x_1 + (a_2 c_2 + b_2)x_2. $$

**Proof.** See Appendix.

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

The optimal strategy of Player $i$ is linear and decreasing in its constellation size $x_i$, and independent of the other player's constellation. The linear form of the policy is due to the linear-quadratic structure of the dynamic optimization problem. The result that $\tilde{s}_i^O(x_i)$ is decreasing in $x_i$ is reminiscent to the concavity of the revenue function.

The main takeaway of Proposition 5 is that the optimal launching policy consists in keeping constant over time the number of operative satellites aloft. Indeed, substituting for $\tilde{s}_i^O(x_i)$ in the dynamics, we obtain

$$ x_i(t+1) = a_i x_i(t) + \tilde{s}_i^O(x_i) = \xi_i^O, $$

for all $t$. This an extreme form of turnpike property, with the state being brought back to the turnpike at each time step. This is due to the fact that the cost of launching satellites has been taken linear.

Consequently, the number of launches at any $t$ can be expressed as

$$ \tilde{s}_i^O(x_i) = (1 - a_i) \tilde{x}_{i\infty}^O, \quad i = 1, 2. $$

Recalling that $a_i = 1 - \alpha_i - \beta_i - \tau y_0$, then the optimal strategy is to replace the satellites that left the constellation either by natural failures at a rate $\beta_i$, or through collisions with debris or other objects at a rate $\tau y_0$, or by passive or active deorbiting at a rate $\alpha_i$.
*Symmetric case.* To get some additional qualitative hints into the results, we consider a symmetric setup, that is, we let

$$a_1 = a_2 = a, \quad c_1 = c_2 = c, \quad \psi_1 = \psi_2 = \psi, \quad \pi_1 = \pi_2 = \pi,$$

and hence $f_1 = f_2 = f$. The results are superscripted with OS for Optimal Symmetric solution. It turns out that the condition for the $x_i$ and $s_i$ to be positive is the same as in the Nash equilibrium of the linear revenue function. Hence we have

**Corollary 1.** If

$$p \le p^N = \frac{\rho\varphi - (1 - \rho a)c}{\rho\eta} = \frac{g}{\rho\eta},$$

then the unique symmetric optimal number of launches by Player i is given by

$$\tilde{s}_i^{OS}(x_i) = -ax_i + \frac{f}{\rho(\psi + 2\tau\pi)}, \quad i = 1, 2,$$

and the steady-state values by

$$\tilde{x}_{i\infty}^{OS} = \frac{f}{\rho(\psi + 2\tau\pi)}, \quad i = 1, 2.$$

The value function is quadratic and given by:

$$\tilde{V}^{OS}(x_1, x_2) = (ac + b)(x_1 + x_2) - \frac{1}{2}\psi(x_1^2 + x_2^2) + 2\tau\pi x_1 x_2 + \frac{1}{1 - \rho}\left(\frac{f^2}{\rho(\psi + 2\tau\pi)} - p\omega\right).$$

**Proof.** By simple substitution. ■

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

The optimal policy $\tilde{s}_i^{OS}(x_i)$ is as before, and for the same reasons, linear and decreasing in own constellation size. As expected, it is increasing in $\varphi$ and decreasing in $\psi$. Further, the larger $\tau\pi$, i.e., the higher the level of substitutability between the two constellations, the lower the number of launches, a standard result in any quantity model à la Cournot. As in the simpler linear case analyzed before, the strategy is decreasing in the cost parameters $c$ and $p$, which also is intuitive.

The upper bound $p^N$ on $p$ is increasing in $\varphi$, $\rho$ and $a$, decreasing in $c$ and $\eta$, hence decreasing in $\beta$ and independent of $\pi$ and $\psi$.

The total discounted payoff over the infinite planning horizon is given by

$$\tilde{V}^{OS}(0, 0) = \frac{1}{1 - \rho} \left( \frac{f^2}{\rho(\psi + 2\tau\pi)} - p\omega \right) = \frac{\rho(\psi + 2\tau\pi)(\tilde{x}_{i\infty}^{OS})^2 - p\omega}{1 - \rho}.$$

We shall compare this value to the corresponding one in the noncooperative game.

What is the impact of product substitutability on the results? When the collision risk between the two constellations is ignored or maneuvering is free ($\pi = 0$), each constellation revenue is independent of the other player's satellites aloft, the interaction between the two players reduces to the congestion effect. When the two products (constellations) are substitutes, which is the case for a positive $\tau\pi$, then they launch less
satellites, and the steady-state values $\tilde{x}_{i\infty}^{OS}$, $i = 1, 2$ are lower, and so are the payoffs. Indeed, we have

$$ \frac{\partial \tilde{s}_i^{OS}}{\partial \pi} < 0, \quad \frac{\partial \tilde{x}_{i\infty}^{OS}}{\partial \pi} < 0, \quad \frac{\partial \tilde{V}^{OS}(0, 0)}{\partial \pi} < 0. $$

Therefore, the negative economic externality brought by strategic product substitutability, comes with a positive one in the form of less space congestion and lower number of debris to be removed.

*Profitability.* We investigate the profitability of the optimal strategies in the symmetric case: it requires that $\tilde{V}^{OS}(0, 0) \ge 0$, i.e.,

$$ \tilde{P}^{O}(p) \triangleq (g - \rho\eta p)^2 - p\rho(\psi + 2\tau\pi)\omega \ge 0, $$

Clearly, $\tilde{P}^{O}(0) = g^2 \ge 0$, and $\tilde{P}^{O}(p^N) = -p\rho(\psi + 2\tau\pi)\omega < 0$. So $\tilde{P}^{O}$ has two positive real roots, the smaller one being less than $p^N$. Let

$$ 2\theta \triangleq \rho(\psi + 2\tau\pi)\omega. \qquad (23) $$

In the symmetric case, there exists a pair of strategies ensuring a positive joint profit if and only if

$$ p \le \tilde{p}_{-}^{O} \triangleq p^{N} - \frac{1}{\rho^{2}\eta^{2}} \left[ \sqrt{(\rho\eta g + \theta)^{2} - g^{2}} - \theta \right]. \qquad (24) $$

## 4.2. Nash equilibrium

Let

$$ \Delta = \psi_1\psi_2 - \tau^2\pi_1\pi_2. $$

For $i \in \{1, 2\}$ and $j = 3 - i$, denote by

$$ \xi_i^N = \frac{\psi_j f_i - \tau \pi_i f_j}{\rho \Delta}. $$

To distinguish with the linear revenue function setup, we tilded the results.

**Proposition 6.** Assuming that the following expression is non-negative, the unique feedback-Nash equilibrium is given by

$$ \tilde{s}_i^N(x_i) = -ax_i + \xi_i^N, \quad i = 1, 2, $$

and the steady-state values by

$$ \tilde{x}_{i\infty}^{N} = \xi_{i}^{N}, \quad i = 1, 2. $$

The value functions are as follows, for $i \in \{1, 2\}$:

$$ \tilde{V}_i^N(x_1, x_2) = -\frac{1}{2}\psi_i x_i^2 - \tau\pi_i x_1 x_2 + (b_i + a_i c_i)x_i + \frac{\psi_i(\psi_j f_i - \tau\pi_i f_j)^2}{2\rho(1-\rho)\Delta^2} - \frac{p\zeta_i\omega}{1-\rho}. $$
**Proof.** See Appendix.

The feedback-Nash equilibrium shares a series of features with the optimal solution. In particular, the strategy is decreasing in each player's own constellation size and independent of the competitor's one. Further, $x_i$ is constant over time and given by

$$ \tilde{x}_i^N(t) = \xi_i^N = \tilde{x}_{i\infty}^N, \quad i = 1, 2. $$

The number of launches is then

$$ \tilde{s}_i^N(x_i) = (1 - a_i) \tilde{x}_{i\infty}^N, \quad i = 1, 2. $$

*Symmetric case.* Again, to get additional insight, we consider a symmetric setup. The results are superscripted with NS for Nash Symmetric equilibrium.

**Corollary 2.** *Assuming that,*

$$ p < p^N = \frac{g}{\rho\eta} = \frac{\rho\varphi - (1 - \rho a)c}{\rho\eta}, $$

the unique symmetric feedback-Nash equilibrium is given by

$$ \tilde{s}_i^{NS}(x_i) = -ax_i + \frac{f}{\rho(\psi + \tau\pi)}, \quad i = 1, 2, $$

and the steady-state values by

$$ \tilde{x}_{i\infty}^{NS} = \frac{f}{\rho(\psi + \tau\pi)}, \quad i = 1, 2. $$

The value functions are as follows, for $i \in \{1, 2\}$, $j = 3 - i$:

$$ \tilde{V}_i^{NS}(x_1, x_2) = (b - ac - \frac{1}{2}\psi x_i - \tau\pi_i x_j)x_i + \frac{\psi f^2}{2\rho(1 - \rho)(\psi + \tau\pi)^2} - \frac{p\omega}{2(1 - \rho)}. $$

**Proof.** By direct substitution.

As before, each player's strategy is decreasing in its own constellation size and independent of the rival's one. Further, the state remains constant over time and is given by

$$ \tilde{x}_i^{NS}(t) = \frac{f}{\rho(\psi + \tau\pi)} = \tilde{x}_{i\infty}^{NS}, \quad i = 1, 2. $$

The number of launches is then

$$ \tilde{s}_i^{NS}(x_i) = \frac{(1 - a)f}{\rho(\psi + \tau\pi)}, \quad i = 1, 2. $$

*Profitability.* Finally, we investigate the profitability of the Nash equilibrium in the symmetric case. The condition $\tilde{V}^N(0, 0) \ge 0$ reads

$$ \tilde{\mathcal{P}}^N(p) \triangleq \psi(g - p\rho\eta)^2 - p\rho(\psi + \tau\pi)^2\omega \ge 0. $$

Clearly, $\tilde{\mathcal{P}}^N(0) = g^2 > 0$ and $\tilde{\mathcal{P}}^N(p^N) = -p^N\rho(\psi + \tau\pi)^2\omega < 0$, so we may again ascertain that $\tilde{\mathcal{P}}^N$ has two real positive roots, the smaller one being less than $p^N$. We conclude that in the symmetric case, there exists a unique Nash equilibrium that is profitable if and only if

$$ p \le \tilde{p}^N \triangleq p^N - \frac{1}{2\psi\eta} \left[ \sqrt{(\psi + \tau\pi)^2\omega[(\psi + \tau\pi)\omega + 4\psi\eta p^N]} - (\psi + \tau\pi)^2\omega \right]. $$
### 4.3. Comparison

Comparing the coordinated and noncooperative solutions, we first notice that the conditions $\xi_i^O \ge 0$ and $\xi_i^N \ge 0$ respectively read

$$ \psi_j f_i - \tau(\pi_i + \pi_j) f_j \ge 0 \quad \text{and} \quad \psi_j f_i - \tau \pi_i f_j \ge 0, $$

so that, as previously, the first one is always the more demanding. More precisely, we evaluate

$$ \xi_i^N - \xi_i^O = \frac{\tau}{\rho\Delta} \left[ \pi_j f_j - \frac{\tau}{\delta}(\pi_1^2 + \pi_1\pi_2 + \pi_2^2)(\psi_j f_i - \tau\pi_i\psi_j) \right], $$

which, given that $\tau$ (and even $\tau\pi$) is a small parameter, is positive but small. This leads to the following conclusions:

1. The players launch more satellites in the Nash equilibrium than in the coordinated solution, but the difference is small. Indeed, we have

$$ \tilde{s}_i^N - \tilde{s}_i^O = \xi_i^N - \xi_i^O > 0. $$

This positive difference is due to the fact that, in the coordinated solution the launching policy by each player internalizes its impact on the other constellation, whereas it does not in the selfish Nash equilibrium.

2. The implication is that the size of the two constellations is larger in the Nash equilibrium than in the coordinated solution, with

$$ \tilde{x}_i^N - \tilde{x}_i^O = \xi_i^N - \xi_i^O > 0. $$

3. As the difference in debris removal at each period is given by:

$$ \tilde{h}^N - \tilde{h}^O = \sum_{i=1}^{2} \eta_i (\tilde{x}_i^N - \tilde{x}_i^O) > 0, $$

we conclude that the ADR cost is higher in the Nash equilibrium than in the coordinated solution.

4. It is easy to verify that the differences $\tilde{s}_i^N - \tilde{s}_i^O$, $\tilde{x}_i^N - \tilde{x}_i^O$, and $\tilde{h}^N - \tilde{h}^O$ are decreasing in $\pi$. In particular, if the two constellations are strategically independent, that is, $\pi_i = \pi_j = 0$, then a noncooperation mode of play exacerbates less space congestion. The ADR cost differential is also lower.

The net result is similar to the linear revenue function case:

**Proposition 7.** In the symmetric case, if $p \le \tilde{p}^N$ so that the Nash equilibrium yields a profitable value for the players, then the optimum coordinated solution is profitable.

In every cases, the total profit under cooperation is larger than the sum of individual noncooperative profits, but possibly by a small amount. The difference $\overline{DC}$ is given by

$$ \overline{DC} = \tilde{V}^O(0, 0) - \sum_{i=1}^{2} \tilde{V}_i^N(0, 0) = \frac{\tau^2}{2\delta\Delta^2\rho(1-\rho)} [\psi_1\psi_2(\psi_2\pi_1^2f_1^2 + \psi_1\pi_2^2f_2^2) - \tau P(\tau)], $$

where $P$ is a second-degree polynomial in $\tau$ with positive constant and leading coefficients.
**Proof.** See appendix.

As a consequence, we conclude, as in the linear revenue function case, that the players should prefer coordination over noncooperation, if they can agree on a sharing of the total cooperative payoff that is mutually beneficial. But again, the benefit might be small.

**5. Active debris removal and taxes**

In this section, we discuss the taxation schemes that aim at financing ADR cost, proposing an alternative to the scheme hypothesized so far.

**5.1. Robustifying the debris taxation scheme**

Until now, we have assumed that the supra-national agency in charge of ADR charges the exact cost of actively removing the expected number of debris created each year. This assumes that some rule has allowed the stakeholders to decide how to share the unattributable part of this debris production, due to collisions among the stock of debris already there at the start of the agreement and to the supposedly small number of extra satellites launched by operators not part of the agreement and taxing scheme.

Actually, this average-cost tax based approach requires that any fixed cost be also shared between the players. But beyond, this scheme may lead to a deficit if random disturbances cause a significant deviation over the predicted mean of debris production, or alternatively to an increased number of debris. To hedge against this case, a provision may be added. One way of “robustifying” the tax would be to implement a cost-plus taxing scheme defined by a tax level $p'$:

$$p' = (1 + \kappa)p, \qquad (25)$$

where $\kappa \in [0, 1]$ is a margin to cover all unforeseen in an accounting sense. Under this tax scheme, the firm is viable only if its expected cumulated discounted profit $V(0, 0)$ is positive. All previous calculations hold substituting $p'$ to $p$. The bounds on $p$ for our proposed solutions to hold now bear on $p'$.

**5.2. Taxing satellites launch**

In the above approach, the players pay a tax to cover the ADR cost. To make a parallel with environmental economics, it is as if the regulator collects taxes from firms to clean the environment rather than taxing their pollutant emissions. An alternative would be to tax satellites launch, which could be possibly easier to implement.

Now, the viability of such a scheme must be investigated for all stakeholders: the profitability for the firms exploiting satellite constellations and the budget equilibrium for the agency in charge of ADR, whose income is the proceeds of the taxes and the expense the cost $ph$ of ADR.

We investigate the feasibility of that approach on the symmetric cases of the previous models, assuming a linear tax $\gamma s$ on launches.

All calculations performed above hold with $c$ replaced by $c^+ = c + \gamma$, and $p$ set to zero, hence $f$ replaced by $g^+$:

$$g^+ = g - (1 - \rho a_i)\gamma \qquad (26)$$
And p now stands for the price the ADR agency pays per debris removal.

The condition for $s^N$, $\xi^0$ and $\xi^N$ to be non-negative is now $g^+ \ge 0$, i.e.,

$$ \gamma \le \gamma^* \triangleq \frac{g}{1 - \rho a} = \frac{\rho \varphi - (1 - \rho a)c}{1 - \rho a}. \qquad (27) $$

while the condition for $s^0$ to be non-negative is the stronger requirement that $g^+ \ge \rho \tau k$, i.e.

$$ \gamma \le \gamma^\dagger = \frac{g - \rho \tau k}{1 - \rho a}. \qquad (28) $$

We will also need the yet stronger condition $g^+ \ge 2\rho\tau k$, i.e.,

$$ \gamma \le \gamma^\ddagger = \frac{g - 2\rho\tau k}{1 - \rho a}. \qquad (29) $$

### 5.2.1. Linear revenue function

Concerning this new taxing scheme, for the symmetric linear revenue function model, we state the following fact:

**Proposition 8.** For the symmetric model with linear revenue function, the launching taxation scheme is profitable for the players under condition (28) in the joint optimization problem and (29) in the Nash equilibrium problem.

Budget equilibrium of the ADR agency is obtained if

$$ \gamma \ge \frac{p\eta\rho}{1 - \rho a}, $$

which is compatible with the preceding conditions provided that p be less than $p^0$ in the joint optimization problem, and less than $p^N - 2\rho\tau k/(1 - \rho a)$ in the Nash equilibrium.

**Proof.** See appendix.

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

### 5.2.2. Strictly concave revenue function

In this section, we set a modification of notation (23) as follows. The parameter $\theta$ will have a different definition for the joint optimization problem than for the Nash equilibrium problem:

$$ 2\theta = \begin{cases} \rho(\psi + 2\tau\pi)\omega & \text{for the joint optimization problem} \\ \rho(\psi + \tau\pi)\omega & \text{for the Nash equilibrium problem} \end{cases} \qquad (30) $$

We state the following result:

**Proposition 9.** For the symmetric model with strictly concave revenue function, the launching taxation scheme is both profitable for the players and budget balancing for the ADR agency if, on the one hand

$$ p \le p^N - \frac{1}{\rho^2 \eta^2} \left[ \sqrt{4\theta(\rho\eta g + \theta)} - \theta \right] < p^N, \qquad (31) $$
and, on the other hand, $\gamma \in [\gamma^{-}, \gamma^{+}]$ where, using $\varepsilon \in \{-, +\}$:

$$ \gamma^{\varepsilon} = \frac{1}{1 - \rho a} \left[ g + \rho \rho \eta + \varepsilon \sqrt{(g - \rho \rho \eta)^2 - 4\theta} \right]. $$

And it holds that $\gamma^{-} \le \gamma^{+} \le \gamma^{*}$.

**Proof.** See appendix.

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

We may expect that the ADR agency would in fact choose a $\gamma$ close to, but not less than, $\gamma^{-}$.

## 6. Conclusion

This paper uses dynamic game theory to analyze how to preserve the common good of outer space while not stopping the development of mega-constellation space activities. Two main conclusions emerge from our modeling: (a) compared to the case where firms coordinate their strategies, competition increases congestion and reduces individual and collective payments, but since the difference between these alternatives may be small, the systematic superiority of coordination over competition is debatable, and (b) it is possible to finance the cleanup of space debris created by these firms through an international agency that would tax them according to the amount of debris they produce or, if the cost of cleanup is not too high, by taxing the satellite launches performed by these firms.

These conclusions are qualitatively robust to changes in the magnitude of the parameters and we believe that they point to research perspectives that could be of interest to other researchers in operational research, economics, applied mathematics or space engineering. Among the possible extensions, we would like to highlight three. First, while we focus on an infinite horizon modeling, it would be interesting to have a short-term analysis to understand and regulate the constitution of these mega-constellations. Then, because of our hypotheses, our analysis is only relevant for large constellations and it would be interesting to analyze the case of smaller constellations (in tens of satellites), which are currently the most numerous. Finally, a more refined, stochastic, analysis of collisions leading to a stochastic modeling of the problems seems to us to be an important point to better understand and regulate the issue of space debris.

On the other hand, our analysis is only a step in our understanding of these issues and it naturally has some limitations. We emphasize some of them and discuss possible extensions.

First, our approach ensures that mega-constellations are neutral in terms of space debris compared to the situation where they would not exist, but it does not solve the problem of financing the further cleanup of space debris. Indeed, various estimates, notably from the IADC (Inter agency Space Debris Coordination Committee), indicate that it is already necessary to remove annually 5 to 10 space debris of more than 10 cm to contain the increase of debris in low earth orbit.

Secondly, nothing is said in this article about uncontrolled entries into the atmosphere and the damage they could cause to the Earth, as for instance the entry on January 24 in 1978 of the Kosmos 954 satellite whose nuclear reactor irradiated a 600 km2 area.
Third, our framework assumes the existence of an international agency, which some specialists believe to be illusory in the short term (Pelton (2020)), and it does not address the question of the possibility of bankruptcy of these firms (such as OneWeb or Leosat), whereas this is what happened in the late 1990s to all of the first generation of large constellation projects (Iridium, Globalstar, ICO, Orbcomm, and Teledesic) (Pelton and Madry (2020)).

Fourth, by focusing on the issue of space debris, our analysis omits the fact that these mega-constellations can also lead to other nuisances. As an illustration of these, we think of the fact that these thousands of satellites placed in low orbit produce a light pollution that is already beginning to disturb the work of astronomers as shown by the call<fnref n="13" /> of the International Astronomical Union of June 3 in 2019, or its communiqué of February 12 in 2020, which clearly warn that both optical astronomy and radio astronomy could be very strongly impacted.

Fifth, the ADR agency is playing a passive role in our framework. An interesting extension would be to consider the agency as a strategic player defining Markovian tax that depends on the stock of debris. Similar work has been done in the literature on emissions control; see, e.g., Benchekroun and Long (1998).

Sixth, our approach of taxing debris production is not the only possible one. It has been suggested that the international community could mimic for low earth orbits the slot allocation process already in force for geo-stationary satellites. This is a more complex undertaking, because “slots” are now more complex orbits. However, how to do that has been investigated in Arnas et al. (2021).

Finally, our model does not consider in details the altitude chosen by each constellation, and treats the 100 km - 2,000 km altitude bands as essentially uniform. Clearly, this is a simplifying assumption that should be dropped in future investigation by considering different altitudes and possibly orbits.<fnref n="14" />

In view of these elements, it seems necessary that the issue of the sustainable development of space be the subject of much greater attention and research.

# 7. References

## References

Adilov, N., Alexander, P., and Cunningham, B. (2015). An economic analysis of earth orbit pollution. *Environmental and Resource Economics*, 60:81–98.

Adilov, N., Alexander, P., and Cunningham, B. (2018). An economic “Kessler Syndrome”: a dynamic model of earth orbit debris. *Economic Letters*, 166:79–82.

Arnas, D., Lifson, M., Linares, R., and Avendaño, M. E. (2021). Definition of low earth orbit slotting architectures using 2d lattice flower constellations. *Advances in Space Research*, 67:3696–3711.

<footnote n="13">¹³Available at https://www.iau.org/news/announcements/detail/ann19035/.</footnote>

<footnote n="14">¹⁴We would like to thank anonymous Reviewers for these last three suggestions.</footnote>
Başar, T. and Olsder, G.-J. (1999). *Dynamic Noncooperative Game Theory*. SIAM, second edition edition.

Béal, S., Deschamps, M., and Moulin, H. (2020). Taxing congestion of the space commons. *Acta Astronautica*, 177:313–319.

Benchekroun, H. and Long, N. V. (1998). Efficiency inducing taxation for polluting oligopolists. *Journal of Public Economics*, 70:325–342.

Bianchessi, N., Cordeau, J. F., Desrosiers, J., Laporte, G., and Raymond, V. (2007). A heuristic for the multi-satellite, multi-orbit and multi-user management of earth observation satellites. *European Journal of Operational Research*, 177(2):750–762.

Bonnal, C., McKnhight, D., Phipps, C., Dupont, C., Missonier, S., Lequette, L., Merle, M., and Rommelaere, S. (2020). Just in time collision avoidance –a review. *Acta Astronautica*, 170:637–651.

Bryce (2021). State of satellite industry report. Technical report, Bryce Tech, Alexandria, VA, USA.

Chen, X., Reinelt, G., Dai, G., and Spitz, A. (2019). A mixed integer linear programming model for multi-satellite scheduling. *European Journal of Operational Research*, 275(2):694–707.

De Frutos, J. and Martín-Herrán, G. (2019). Spatial vs. non-spatial transboundary pollution control in a class of cooperative and non-cooperative dynamic games. *European Journal of Operational Research*, 276(1):379–394.

El Ouardighi, F., Kogan, K., Gnecco, G., and Sanguineti, M. (2020). Transboundary pollution control and environmental absorption efficiency management. *Annals of Operations Research*, 287(2):653–681.

Engwerda, J. (2005). *LQ Dynamic Optimization and Differential Games*. John Wiley & Sons.

Fliege, J., Kaparis, K., and Khosravi, B. (2012). Operations research in the space industry. *European Journal of Operational Research*, 217(2):233–240.

Foreman, V. L., Siddiqi, A., and de Weck, O. L. (2017). Large satellite constellation orbital debris impacts: Case studies of OneWeb and SpaceX proposals. In *AIAA SPACE and Astronautics Forum and Exposition*, pages 5200–52017.

Frischmann, B., Marciano, A., and Ramello, G. (2019). Tragedy of the commons after 50 years. *Journal of Economic Perspectives*, 33:211–228.

Grzelka, Z. and Wagner, J. (2019). Managing satellite debris in low-earth orbit: incentivizing ex ante satellite quality and ex post take-back programs. *Environmental and Resource Economics*, 74:319–336.

Hardin, G. (1968). The tragedy of the commons. *Science*, 162:1243–1248.
Haurie, A., Krawczyk, J. B., and Zaccour, G. (2012). *Games and Dynamic Games*. Scientific World, Singapore.

Jørgensen, S., Martin-Herrán, G., and Zaccour, G. (2010). Dynamic games in the economics and management of pollution. *Environmental Modelling and Assessment*, 15:433–467.

Kessler, D. (1991). Collision cascading: the limits of population growth in low earth orbit. *Advances in Space Research*, 11–12:63–66.

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

Klima, R., Bloembergen, D., Savani, R., Tuyls, K., Hennes, D., and Izzo, D. (2016). Space debris removal: a game theoretic analysis. *Games*, 7:18 p.

Klinkrad, D. (2015). The hazard of orbital debris. In Pelton, J. and Allahdadi, F., editors, *Handbook of cosmic Hazard and Planetary Defense*, pages 807–834. Springer.

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

Liou, J.-C. (2020). Risk from orbital debris. Technical report, Brookings Institution.

Long, N.-V. (2011). Dynamic games in the economics of natural resources: A survey. *Dynamic Games and Applications*, 1:115–148.

Macauley, M. (2015). The economics of space debris: estimating the costs and benefits of debris mitigation. *Acta Astronautica*, 115:160–164.

McDowel, J. (2018). The edge of space: revisiting the karman line. *Acta Astronautica*, 151:668–677.

Morgan, L., Morton, A., and Daniels, R. (2006). Simultaneously determining the mix of space launch vehicles and the assignment of satellites to rockets. *European Journal of Operational Research*, 172(3):747–760.

Muller, C., Rozanova, O., and Urdanoz, M. (2011). Economic valuation of debris removal. In *Limiting future collision risk to spacecraft: an assessment of NASA's meteoroid and orbital debris programs*, Washington DC, USA. 68th international astronautical congress, National research council, The National Academies Press.

Nagendra, N., Narayanamurthy, G., and Moser, R. (2020). Management of humanitarian relief operations using satellite big data analytics: the case of kerala floods. *Annals of Operations Research*.

OECD (2019). The space economy in figures: how space contributes to the global economy. Technical report, OECD.
Oltrogge, D. and Christensen, I. (2020). Space governance in the new space era. *Journal of Space Safety Engineering*, 7:432–438.

Pardini, C. and Anselmo, L. (2014). Review of past on-orbit collisions among cataloged objects and examination of the catastrophic fragmentation concept. *Acta Astronautica*, 100:30–39.

Pardini, C. and Anselmo, L. (2021). Evaluating the impact of space activities in low earth orbit. *Acta Astronautica*, 184:11–22.

Pelton, J. (2020). Global insurance fund for orbital space debris removal. *Journal of Space Safety Engineering*, 7:99–100.

Pelton, J. and Madry, S. (2020). Introduction to the small satellite revolution and its many implications. In Pelton, J. and Madry, S., editors, *Handbook of small satellites*. Springer.

Priyant, M. C. and Kamath, S. (2019). Review on active space debris removal methods. *Space Policy*, 47:194–206.

Rao, A., Burgess, M., and Kaffine, D. (2020). Orbital-use fees could more than quadruple the value of the space industry. *Proceedings of the National Academy of Sciences*, 117:12756–12762.

Rouillon, S. (2020). A physico-economic model of low earth orbit management. *Environmental and Resource Economics*, 77:695–723.

Sandler, T. and Schulze, W. (1981). The economics of outer space. *Natural Resources Journal*, 21:371–393.

Tangpattanakul, P., Nicolas, J., and Pierre, L. (2015). A multi-objective local search heuristic for scheduling earth observations taken by an agile satellite. *European Journal of Operational Research*, 245(2):542–554.

Undseth, M., Jolly, C., and Olivari, M. (2020). Space sustainability: the economics of space debris in perspective. Technical Report 87, OECD Science, technology and industry, Policy Papers.

Weinzierl, M. (2018). Space, the final economic frontier. *Journal of Economic Perspectives*, 32:173–192.
# Appendix A. Proofs

## Appendix A.1. Notation

In these proofs, we will make use of vectorial notation as follows. The letters $x$, $s$, $b$, $c$, $f$, $g$, $\varphi$, $\eta$, $\zeta$ without indices stand for 2-vectors whose coordinates are the same letters indexed. The index $j$ always means $3-i$. The notation $\langle \cdot, \cdot \rangle$ stands for the scalar product:

$$ \langle \varphi, x \rangle = \varphi_1 x_1 + \varphi_2 x_2 . $$

We also let

$$ A = \operatorname{diag}\{a_1, a_2\} = \begin{pmatrix} a_1 & 0 \\ 0 & a_2 \end{pmatrix}, \quad D = \operatorname{diag}\{d_1, d_2\} = \begin{pmatrix} d_1 & 0 \\ 0 & d_2 \end{pmatrix}, $$

$$ \varphi^1 = \begin{pmatrix} \varphi_1 \\ 0 \end{pmatrix}, \quad \varphi^2 = \begin{pmatrix} 0 \\ \varphi_2 \end{pmatrix}, $$

and likewise for $c^1, c^2, \eta^1, \eta^2, b^1, b^2, f^1, f^2$, also $D_1 = \operatorname{diag}\{d_1, 0\}, D_2 = \operatorname{diag}\{0, d_2\}$. An exception is

$$ k^1 = \begin{pmatrix} 0 \\ k_1 \end{pmatrix}, \quad k^2 = \begin{pmatrix} k_2 \\ 0 \end{pmatrix}, \quad k = k^1 + k^2 = \begin{pmatrix} k_2 \\ k_1 \end{pmatrix}. $$

Observe that $\varphi = \varphi^1 + \varphi^2$, and likewise for $c, \eta, b, f$, and $D$. We also set

$$ Q^1 = \begin{pmatrix} \psi_1 & \tau\pi_1 \\ \tau\pi_1 & 0 \end{pmatrix}, \quad Q^2 = \begin{pmatrix} 0 & \tau\pi_2 \\ \tau\pi_2 & \psi_2 \end{pmatrix}, $$

and $Q = Q^1 + Q^2$. Using these notation, and $b^i = \varphi^i - p\eta^i$, the individual profits write, in the linear revenue function model

$$ \Pi_i = \sum_{t=0}^{\infty} \rho^t \left[ \langle b^i - \tau k^i, x \rangle - \langle c^i, s \rangle - \frac{1}{2} \langle s, D_i s \rangle - p \zeta_i \omega \right], $$

and in the strictly concave revenue function model

$$ \tilde{\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]. $$

*Dynamics.* For all the problems considered, the dynamics are

$$ x(t+1) = Ax(t) + s(t), \quad x(0) = x_0 = 0. $$

## Appendix A.2. Proof of proposition 1

For the joint optimization problem, the unique criterion is

$$ \Pi^O = \sum_{t=0}^{\infty} \rho^t \left[ \langle b - \tau k, x \rangle - \langle c, s \rangle - \frac{1}{2} \langle s, Ds \rangle - p\omega \right]. $$
Consider the problem with an arbitrary initial state $x_0$, and let

$$V^O(x_0) = \max_{s(\cdot)\ge 0} \Pi^0.$$

We make the informed guess that $V^O$ is of the form

$$V^O(x) = \langle \ell^O, x \rangle + m^O,$$

and write Bellman's equation:

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

i.e., using $\langle \ell^O, Ax \rangle = \langle A^t \ell^O, x \rangle$ as we will do in all the sequel<fnref n="15" />

$$V^O(x) = \max_{s \ge 0} \left\{ -\frac{1}{2} \left[ \langle s, Ds \rangle + 2 \langle c - \rho \ell^O, s \rangle \right] + \langle b + \rho A^t \ell^O - \tau k, x \rangle - p\omega + \rho m^O \right\}.$$

Completing the square, this may also be written

$$V^O = \max_{s \ge 0} \left\{ -\frac{1}{2} \langle s - s^O, D(s - s^O) \rangle \right\} + \frac{1}{2} \langle s^O, Ds^O \rangle + \langle b + \rho A^t \ell^O - \tau k, x \rangle - p\omega + \rho m^O,$$

with

$$s^O = D^{-1}(\rho \ell^O - c).$$

The matrix $D$ is positive definite. Hence the maximum is obviously reached at $s = s^O$, so that

$$\langle \ell^O, x \rangle + m^O = \langle b + \rho A^t \ell^O - \tau k, x \rangle + \frac{1}{2} \langle \rho \ell^O - c, D^{-1}(\rho \ell^O - c) \rangle - p\omega + \rho m^O.$$

Identifying terms, it follows that

$$\ell^O = b + \rho A^t \ell^O - \tau k \implies \ell_i^O = \frac{b_i - \tau k_j}{1 - \rho a_i},$$

$$m^O = \frac{1}{2} \langle \rho \ell^O - c, D^{-1}(\rho \ell^O - c) \rangle - p\omega + \rho m^O.$$

Using

$$\rho \ell^O - c = (I - \rho A^t)^{-1} [\rho b - (I - \rho A^t)c - \rho \tau k] = (I - \rho A^t)^{-1} (f - \rho \tau k),$$

we get

$$s^O = D^{-1}(I - \rho A^t)^{-1}(f - \rho \tau k) \implies s_i^O = \frac{f_i - \rho \tau k_j}{d_i(1 - \rho a_i)},$$

$$\begin{aligned} m^O &= \frac{1}{1-\rho} \left[ \frac{1}{2} \left\langle (I - \rho A^t)^{-1}(f - \rho \tau k), D^{-1}(I - \rho A^t)^{-1}(f - \rho \tau k) \right\rangle - p\omega \right] \\ &= \frac{1}{1-\rho} \left[ \sum_{i=1}^{2} \frac{1}{2d_i} \left( \frac{f_i - \rho \tau k_j}{1 - \rho a_i} \right)^2 - p\omega \right]. \end{aligned}$$

<footnote n="15">¹⁵We keep the upper index $t$ on $A$ meaning $A$ transpose although our matrix $A$ is symmetric, because the calculation using the vector and matrix notation holds for more general dynamics and criterion.</footnote>
Appendix A.3. *Proof of proposition 2*

We set, for an arbitrary initial state $x_0$, player j playing its equilibrium strategy $s_j^N$:

$$V_i^N(x_0) = \max_{s_i \ge 0} \Pi_i(s_i, s_j^N).$$

Proceeding as previously, we guess a form $V_i^N(x) = \langle \ell^i, x \rangle + m^i$ and write Isaacs' equation:

$$V_i^N(x) = \max_{s_i \ge 0} \left[ -\frac{d_i}{2} s_i^2 + (\rho \ell_i^i - c_i) s_i \right] + \langle b^i - \tau k^i + \rho A^t \ell^i, x \rangle + \rho \ell_j^i s_j^N - p \zeta_i \omega + \rho m^i.$$

or, completing the square

$$V_i^N(x) = \max_{s_i \ge 0} \left[ -\frac{d_i}{2}(s_i - s_i^N)^2 \right] + \frac{d_i}{2}(s_i^N)^2 + \langle b^i - \tau k^i + \rho A^t \ell^i, x \rangle + \rho \ell_j^i s_j^N - p \zeta_i \omega + \rho m^i,$$

with

$$s_i^N = \frac{\rho \ell_i^i - c_i}{d_i},$$

which provides the maximum sought. The same reasoning applies to $s_j^N$. We may substitute these in $V_i^N$:

$$V_i^N = \langle b^i - \tau k^i + \rho A^t \ell^i, x \rangle + \frac{(\rho \ell_i^i - c_i)^2}{2d_i} + \rho \ell_j^i \frac{\rho \ell_j^j - c_j}{d_j} - \zeta_i p \omega + \rho m^i.$$

Identifying with $\langle \ell^i, x \rangle + m^i$, we obtain

$$\ell^i = b^i - \tau k^i + \rho A^t \ell^i \implies \ell^i = (I - \rho A^t)^{-1}(b^i - \tau k^i) = \begin{pmatrix} \frac{b_i}{1 - \rho a_i} \\ \frac{-\tau k_i}{1 - \rho a_j} \end{pmatrix},$$

and hence

$$s_i^N = \frac{\rho b_i - (1 - \rho a_i)c_i}{d_i(1 - \rho a_i)} = \frac{f_i}{d_i(1 - \rho a_i)},$$

and

$$m^i = \frac{1}{1 - \rho} \left[ \frac{1}{2d_i} \left( \frac{f_i}{1 - \rho a_i} \right)^2 - \frac{\rho \tau k_i f_j}{d_j(1 - \rho a_j)^2} - \zeta_i p \omega \right].$$

Appendix A.4. *Proof of proposition 5*

The criterion to be optimized is $\Pi = \Pi_1 + \Pi_2$, i.e.:

$$\Pi = \sum_{t=0}^{\infty} \rho^t \left[ -\frac{1}{2} \langle x, Qx \rangle - \langle c, s \rangle + \langle b, x \rangle - p\omega \right].$$
Again, we make the informed guess that Bellman's return function will be of the form

$$ \tilde{V}^0(x) = -\frac{1}{2} \langle x, Qx \rangle + \langle \tilde{l}, x \rangle + \tilde{m}. $$

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

$$ \tilde{V}^0(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, Q(Ax + s) \rangle + \langle \tilde{l}, Ax + s \rangle + \tilde{m} \right] \right\}. $$

Regroup terms with s:

$$ \tilde{V}^0(x) = \max_{s \ge 0} \left\{ -\frac{1}{2} \langle s, \rho Qs \rangle + \langle \rho \tilde{l} - \rho QAx - c, s \rangle \right\} \\ - \frac{1}{2} \langle x, (Q + \rho A'QA)x \rangle + \langle b + \rho A'\tilde{l}, x \rangle - p\omega + \rho \tilde{m}. $$

Completing the square, this reads

$$ \tilde{V}^0(x) = \max_{s \ge 0} \left\{ -\frac{1}{2} \langle s - \tilde{s}^0, \rho Q(s - \tilde{s}^0) \rangle \right\} + \frac{1}{2} \langle \tilde{s}^0, \rho Q \tilde{s}^0 \rangle \\ - \frac{1}{2} \langle x, Qx \rangle - \frac{1}{2} \langle Ax, \rho QAx \rangle + \langle b + \rho A' \tilde{l}, x \rangle - p\omega + \rho \tilde{m}, $$

with

$$ \tilde{s}^0 = -Ax + \frac{1}{\rho} Q^{-1} (\rho \tilde{l} - c). $$

Expanding the square term in $\tilde{s}^0$, it follows

$$ \tilde{V}^0(x) = \frac{1}{2} \langle Ax, \rho QAx \rangle - \langle \rho \tilde{l} - c, Ax \rangle + \frac{1}{2\rho} \langle \rho \tilde{l} - c, Q^{-1}(\rho \tilde{l} - c) \rangle \\ - \frac{1}{2} \langle x, Qx \rangle - \frac{1}{2} \langle Ax, \rho QAx \rangle + \langle b + \rho A' \tilde{l}, x \rangle - p\omega + \rho \tilde{m}, $$

hence,

$$ \tilde{V}^0(x) = -\frac{1}{2} \langle x, Qx \rangle + \langle b + A'c, x \rangle + \frac{1}{2\rho} \langle \rho \tilde{l} - c, Q^{-1}(\rho \tilde{l} - c) \rangle - p\omega + \rho \tilde{m}. $$

The square term in $x$ is as “guessed”, and identifying other coefficients, we find

$$ \tilde{l} = b + A'c \implies \rho \tilde{l} - c = \rho b - (I - \rho A')c = f, $$

and thus

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

and

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

These formulas coincide with the formulas expanded componentwise in the proposition.
Appendix A.5. *Proof of proposition 6*

We seek an Isaacs Value function of the form

$$ \tilde{V}_i^N(x) = \max_{s_i \ge 0} \tilde{\Pi}_i(s_i, \tilde{s}_j^N) = -\frac{1}{2} \langle x, Q^i x \rangle + \langle \tilde{l}^i, x \rangle + \tilde{m}^i, $$

with moreover

$$ \tilde{l}^1 = \begin{pmatrix} \tilde{l}_1 \\ 0 \end{pmatrix}, \quad \tilde{l}^2 = \begin{pmatrix} 0 \\ \tilde{l}_2 \end{pmatrix}, \quad \tilde{l}^N = \begin{pmatrix} \tilde{l}_1 \\ \tilde{l}_2 \end{pmatrix}. $$

so that $\langle \tilde{l}^i, x \rangle = \tilde{l}_i x_i$. Again, the fact we do solve Isaacs' equation with this "guess" will prove it correct. Isaacs' equation is now, for $i \in \{1, 2\}$:

$$ \tilde{V}_i^N(x) = \max_{s_i \ge 0} \left\{ -\frac{1}{2} \langle x, Q^i x \rangle + \langle b^i, x \rangle - \langle c^i, s \rangle - \zeta_i p \omega \right. \\ \left. + \rho \left[ -\frac{1}{2} \langle Ax + s, Q^i(Ax + s) \rangle + \langle \tilde{l}^i, Ax + s \rangle + \tilde{m}^i \right] \right\}. $$

or, reordering terms

$$ \tilde{V}_i^N(x) = \max_{s_i \ge 0} \left\{ -\frac{1}{2} \langle s, \rho Q^i s \rangle + \langle \rho \tilde{l}^i - c^i - \rho Q^i Ax, s \rangle \right\} \\ - \frac{1}{2} \langle x, Q^i x \rangle - \frac{1}{2} \langle Ax, \rho Q^i Ax \rangle + \langle \rho A^t \tilde{l}^i + b^i, x \rangle - \zeta_i p \omega + \rho \tilde{m}^i. $$

The matrix $Q^i$ is not positive definite, but the term in curly braces 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^t$ 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^t s_i$, and that of $\langle Q^i Ax, s \rangle$ is $Q_i^t Ax$. To exploit this fact, we introduce the (non-symmetric) matrix

$$ \tilde{Q} = \begin{pmatrix} Q_1^1 \\ Q_2^1 \end{pmatrix} = \begin{pmatrix} \psi_1 & \tau \pi_1 \\ \tau \pi_2 & \psi_2 \end{pmatrix}, $$

whose determinant is $\Delta = \psi_1 \psi_2 - \tau^2 \pi_1 \pi_2$. Differentiating in both $\tilde{V}_i^N(x)$ for $i \in \{1, 2\}$ we get

$$ \tilde{Q}s = \rho(\tilde{l}^N - \tilde{Q}Ax) - c, $$

hence, using the notation

$$ \rho \tilde{l}^N - c = \tilde{f} = \begin{pmatrix} \tilde{f}_1 \\ \tilde{f}_2 \end{pmatrix}, \quad \rho \tilde{l}^1 - c^1 = \tilde{f}^1 = \begin{pmatrix} \tilde{f}_1 \\ 0 \end{pmatrix}, \quad \rho \tilde{l}^2 - c^2 = \tilde{f}^2 = \begin{pmatrix} 0 \\ \tilde{f}_2 \end{pmatrix}, $$

$$ s = \tilde{s}^N = -Ax + \frac{1}{\rho} \tilde{Q}^{-1} (\rho \tilde{l}^N - c) = -Ax + \frac{1}{\rho} \tilde{Q}^{-1} \tilde{f}. $$
Substitute in $\tilde{V}^N$:

$$ 
\begin{aligned}
\tilde{V}^N(x) = & -\frac{1}{2} \langle Ax, \rho Q^i Ax \rangle + \langle Ax, Q^i \tilde{Q}^{-1} \tilde{f} \rangle - \frac{1}{2\rho} \langle \tilde{Q}^{-1} \tilde{f}, Q^i \tilde{Q}^{-1} \tilde{f} \rangle \\
& + \langle \rho Q^i Ax, Ax \rangle - \langle \tilde{f}^i, Ax \rangle + \frac{1}{\rho} \langle \tilde{f}^i, \tilde{Q}^{-1} \tilde{f} \rangle - \langle Q^i Ax, \tilde{Q}^{-1} \tilde{f} \rangle \\
& - \frac{1}{2} \langle x, Q^i x \rangle - \frac{1}{2} \langle Ax, \rho Q^i Ax \rangle + \langle \rho A^t \tilde{\ell}^i + b^i, x \rangle - \zeta_i p\omega + \rho \tilde{m}^i.
\end{aligned}
$$

Cancelling terms that cancel and substituting for $\tilde{f}^i$, we end up with

$$ \tilde{V}_i^N(x) = -\frac{1}{2} \langle x, Q^i x \rangle + \langle b^i + A^t c^i, x \rangle - \frac{1}{2\rho} \langle \tilde{Q}^{-1} \tilde{f}, Q^i \tilde{Q}^{-1} \tilde{f} \rangle + \frac{1}{\rho} \langle \tilde{f}^i, \tilde{Q}^{-1} \tilde{f} \rangle - \zeta_i p\omega + \rho \tilde{m}^i. $$

Identifying coefficients, we find that the square term in $x$ is as “guessed”, and also

$$ \tilde{\ell}^i = b^i + A^t c^i \implies \tilde{\ell}^1 = \begin{pmatrix} b_1 + a_1 c_1 \\ 0 \end{pmatrix}, \quad \tilde{\ell}^2 = \begin{pmatrix} 0 \\ b_2 + a_2 c_2 \end{pmatrix}. $$

Substituting in $\tilde{f}$, it follows that $\tilde{f} = f$ according to our standard definition. Hence also

$$ \tilde{m}^i = \frac{1}{1-\rho} \left[ \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 - \zeta_i p\omega \right]. $$

Finally, using

$$ \tilde{Q}^{-1} f = \frac{1}{\Delta} \begin{pmatrix} \psi_2 f_1 - \tau \pi_1 f_2 \\ -\tau \pi_2 f_1 + \psi_1 f_2 \end{pmatrix}, $$

it is a straightforward calculation, though cumbersome, to check that according to the statement to prove

$$ \tilde{m}^i = \frac{1}{1-\rho} \left[ \frac{\psi_i (\psi_j f_i - \tau \pi_i f_j)^2}{2\rho \Delta^2} - \zeta_i p\omega \right]. $$

## Appendix A.6. Proof of proposition 7

Concerning the symmetric case, we notice that

$$ \psi \tilde{P}^0(p) = \tilde{P}^N(p) + p \rho \tau^2 \pi^2 \omega. $$

Therefore, if $\tilde{P}^N(p)$ is positive, so is $\tilde{P}^0(p)$. Thus $\tilde{p}_-^N < \tilde{p}_-^O$.

In the general case, the total profit under cooperation being the maximum of the sum of the profits of the players over all pairs 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:

$$ \overline{DC} = \tilde{V}^O(0) - \sum_{i=1}^{2} \tilde{V}_i^N(0). $$
Let us calculate the sum $\sum_{i=1}^{2} \tilde{V}_{i}^{N}(0)$, using the fact that $\sum_{i=1}^{2} f^{i} = f$ and $\sum_{i=1}^{2} Q^{i} = Q$:

$$ \sum_{i=1}^{2} \tilde{V}_{i}^{N}(0) = \frac{1}{\rho(1-\rho)} \left[ \langle f, \tilde{Q}^{-1}f \rangle - \frac{1}{2} \langle \tilde{Q}^{-1}f, Q\tilde{Q}^{-1}f \rangle \right]. $$

We write the first scalar product as

$$ \langle f, \tilde{Q}^{-1}f \rangle = \langle \tilde{Q}^{-1}f, \tilde{Q}\tilde{Q}^{-1}f \rangle = \langle \tilde{Q}^{-1}f, \tilde{Q}^{s}\tilde{Q}^{-1}f \rangle, $$

where $\tilde{Q}^s$ is the symmetric part of $\tilde{Q}$:

$$ \tilde{Q}^s = \begin{pmatrix} \psi_1 & \tau \frac{\pi_1+\pi_2}{2} \\ \tau \frac{\pi_1+\pi_2}{2} & \psi_2 \end{pmatrix}. $$

Then,

$$ \tilde{Q}^s - \frac{1}{2}Q = \frac{1}{2} \begin{pmatrix} \psi_1 & 0 \\ 0 & \psi_2 \end{pmatrix} \triangleq \frac{1}{2}\Psi. $$

Therefore,

$$ \sum_{i=1}^{2} \tilde{V}_{i}^{N}(0) + \frac{p\omega}{1-\rho} = \frac{1}{2\rho(1-\rho)} \langle \tilde{Q}^{-1}f, \Psi\tilde{Q}^{-1}f \rangle = \frac{1}{2\rho(1-\rho)} \langle f, \tilde{Q}^{-t}\Psi\tilde{Q}^{-1}f \rangle, $$

and

$$ \overline{DC} = \frac{1}{2\rho(1-\rho)} \langle f, (Q^{-1} - \tilde{Q}^{-t}\Psi\tilde{Q}^{-1})f \rangle. $$

There remains to perform the straightforward computation

$$ \tilde{Q}^{-1} - \tilde{Q}^{-t}\Psi\tilde{Q}^{-1} = \frac{1}{\Delta^2\delta}R, $$

where R is the symmetric matrix with coefficients as follows:

$$ R_{11} = \tau^2 \psi_1 \psi_2^2 \pi_1^2 + \tau^4 \psi_2 \pi_2^2 (2\pi_1^2 + 2\pi_1 \pi_2 + \pi_2^2), $$

$$ R_{12} = R_{21} = -\tau^3 \psi_1 \psi_2 (\pi_1^3 + \pi_1^2 \pi_2 + \pi_1 \pi_2^2 + \pi_2^3) - \tau^5 \pi_1^2 \pi_2^2 (\pi_1 + \pi_2), $$

$$ R_{22} = \tau^2 \psi_1^2 \psi_2 \pi_2^2 + \tau^4 \psi_1 \pi_1^2 (\pi_1^2 + 2\pi_1 \pi_2 + 2\pi_2^2). $$

## Appendix A.7. Proof of proposition 8

In this section, we adopt the notation of section 5. That is, we are in the symmetric case. The un-indexed letters x, s, a, g, η are the scalar common variables or parameters of the two players. Moreover, the notation $s^*$ stands for $s^O$ or $s^N$ depending on whether we consider the joint optimization problem or the Nash equilibrium problem.

**Profitability.** In the joint optimization problem, the condition $g^+ - \rho\tau k \ge 0$ is imposed by the non-negativity of $s^O$. Furthermore, it follows from the equalities

$$ V^O(0, 0) = \frac{(g^+ - \rho\tau k)^2}{2d(1-\rho)(1-\rho a)^2} \quad \text{and} \quad V^N(0, 0) = \frac{g^+(g^+ - 2\rho\tau k)}{2d(1-\rho)(1-\rho a)^2} $$

that the new taxation scheme is profitable for the players if $V(0, 0)$ is non-negative, i.e., under condition (28) for the joint optimization problem, or (29) for the Nash equilibrium.
*Budget Equilibrium.* In the linear revenue function model, the optimal strategies are constant at $s = s^*$. Therefore, the cumulated discounted tax levied is

$$T = 2 \frac{\gamma s^*}{1 - \rho}.$$

The trajectories followed from $x_i(0) = 0$ are, for both players

$$x_i(t) = (1 + a + \cdots + a^{t-1})s^* = \frac{1 - a^t}{1 - a}s^*.$$

Therefore, the expected cumulated discounted ADR expense is

$$\text{ADR} = \frac{2p\eta s^*}{1-a} \sum_{t=0}^{\infty} \rho^t (1-a^t) = \frac{2p\eta s^*}{1-a} \left[ \frac{1}{1-\rho} - \frac{1}{1-\rho a} \right] = \frac{2p\eta\rho s^*}{(1-\rho)(1-\rho a)}.$$

Budget equilibrium is reached provided that

$$[\gamma(1 - \rho a) - p\eta\rho]s^* \geq 0.$$

With the added condition that $s^* \geq 0$, this is ensured by the condition

$$\gamma \geq \frac{p\eta\rho}{1 - \rho a}.$$

Together with condition (28) or (29), depending on which applies, this imposes the limits

$$p \leq \frac{g - \rho \tau k}{\rho \eta} = p^0 \quad \text{or} \quad p \leq \frac{g - 2\rho \tau k}{\rho \eta} = p^N - \frac{2\rho \tau k}{1 - \rho a},$$

respectively.

## Appendix A.8. *Proof of proposition 9*

In this section, the convention on notations is as in the preceding one, concerning the symmetric model.

*Profitability.* We now have

$$\tilde{V}^O(0, 0) = \frac{(g^+)^2}{\rho(1 - \rho)(\psi + 2\tau\pi)}, \qquad \tilde{V}^N(0, 0) := \frac{\psi(g^+)^2}{\rho(1 - \rho)(\psi + \tau\pi)^2},$$

which are always non-negative. Hence the only condition for the profitability of the scheme is that $g^+$ be non-negative, i.e., $\gamma \leq \gamma^*$, given by equation (27).

*Budget equilibrium.* Notation $\xi$ now holds either for $\xi^O$ for the joint optimization problem, or for $\xi^N$ for the Nash equilibrium. With the simple strategy $s_i = -a_i x_i + \xi_i$ coupled with $x_i(0) = 0$, and $\forall t > 0$, $x_i(t) = \xi_i$, the part $\gamma \xi_i$ of the taxes is collected at each time step, while the part $-\gamma_i a_i x_i$ is collected from step one on. Therefore, the cumulated discounted taxes levied in the symmetric case is

$$T = \frac{2\gamma}{1 - \rho}(1 - \rho a)\xi.$$
In a similar fashion, the total ADR expense for the agency is

$$ADR = \frac{2p}{1-\rho} \left( \rho\eta\xi + p\frac{\omega}{2} \right).$$

Hence, budget equilibrium is ensured if

$$\gamma(1 - \rho a)\xi - p(\rho\eta\xi + p\frac{\omega}{2}) \geq 0,$$

where

$$\xi = \frac{g^{+}\omega}{\theta} = \frac{(g - (1 - \rho a)\gamma)\omega}{\theta}.$$

This yields a second degree polynomial $Q(\gamma)$ which must be non-negative:

$$Q(\gamma) \triangleq -(1 - \rho a)^2 \gamma^2 + (1 - \rho a)[g + p\rho\eta]\gamma - p(g\rho\eta + \theta) \geq 0.$$

This is possible only if the discriminant of $Q$ is non-negative, i.e.

$$\mathcal{R}(p) \triangleq \rho^2 \eta^2 p^2 - 2(\rho\eta g + 2\theta)p + g^2 \geq 0.$$

Note that $\mathcal{R}(p^N) = -4\theta < 0$ and $\mathcal{R}(0) = g^2 > 0$. Therefore $\mathcal{R}$ has a positive real root less than $p^N$, and its positivity is ensured by inequality (31):

$$p \leq p^N - \frac{1}{\rho^2 \eta^2} \left[ \sqrt{4\theta(\rho\eta g + \theta)} - \theta \right].$$

It follows from the remark that $\mathcal{R}(p) = \tilde{\mathcal{P}}^0(p) - 2\theta p$ that this new condition is strictly more restrictive than that for a profitable solution in the previous taxing scheme.

Finally, the constraint it creates on $\gamma$ is $\gamma \in [\gamma^-, \gamma^+]$ where, using $\varepsilon \in \{-, +\}$:

$$\gamma^{\varepsilon} = \frac{1}{1 - \rho a} \left[ g + p\rho\eta + \varepsilon \sqrt{(g - p\rho\eta)^2 - 4\theta} \right].$$

Observe that $Q(\gamma^*) = -p\theta < 0$ and $Q'(\gamma^*) = (1 - \rho a)(p\rho\eta - g) < 0$ for $p < p^N$. Therefore $\gamma^* > \gamma^+$, and the constraint (31) and $\gamma^- \leq \gamma \leq \gamma^+$ together ensure both the profitability of the scheme for the players and a positive balance sheet for the ADR agency.