22 A6.IPB 73 Paris, France 18-22 September 2022 2022the authors

Celina Pasiecznik1 Andrea D’Ambrosio2 Daniel Jang3 Richard Linares4

SM Candidate, Department of Aeronautics and Astronautics, Massachusetts
Institute of Technology, MA 02139<EMAIL_ADDRESS>Postdoctoral
Associate, Department of Aeronautics and Astronautics, Massachusetts Institute
of Technology, MA 02139<EMAIL_ADDRESS>PhD Candidate, Department of
Aeronautics and Astronautics, Massachusetts Institute of Technology, MA 02139,
<EMAIL_ADDRESS>Rockwell International Career Development Professor,
Associate Professor of Aeronautics and Astronautics, Department of Aeronautics
and Astronautics, Massachusetts Institute of Technology, MA 02139,
<EMAIL_ADDRESS>
# A Dynamical Systems Analysis of the Effects of the Launch Rate Distribution
on the Stability of a Source-Sink Orbital Debris Model

###### Abstract

Future launches are projected to significantly increase both the number of
active satellites and aggregate collision risk in Low Earth Orbit (LEO).
Ensuring the long-term sustainability of the space environment demands an
accurate model to understand and predict the effect of launch rate
distribution as a major driver of the evolution of the LEO orbital population.
In this paper, a dynamical systems theory approach is used to analyze the
effect of launch rate distribution on the stability of the LEO environment. A
multi-shell, three-species source-sink model of the LEO environment referred
to as MOCAT-3 for MIT Orbital Capacity Assessment Tool - 3 Species, is used to
study the evolution of the species populations. The three species included in
the model are active satellites, derelict satellites, and debris. Each shell
is modeled by a system of three equations, representing each species, that are
coupled through coefficients related to atmospheric drag, collision rate, mean
satellite lifetime, post-mission disposal probability, and active debris
removal rate. The major sink in the model is atmospheric drag, whereas the
only source apart from collision fragments is the launch rate, making it the
critical manageable factor impacting the orbital capacity. Numerical solutions
of the system of differential equations are computed, and an analysis of the
stability of the equilibrium points is conducted for numerous launch rate
distributions. The stability of the equilibrium points is used to test the
sensitivity of the environment to run-away debris growth, known as Kessler
syndrome, that occurs at the instability threshold. Various bounding cases are
studied from business-as-usual launch rates based on historic launch data, to
high launch rates wherein a fraction of the satellite proposals filed with the
International Telecommunication Union (ITU) are launched. An analysis of the
environment’s response to perturbations in launch rate and debris population
is conducted. The maximum perturbation in the debris population from the
equilibrium state, for which the system remains in a stable configuration, is
calculated. Plots of the phase space about the equilibrium points are
generated. The results will help to better understand the orbital capacity of
LEO and the stability of the space environment, as well as provide improved
guidelines on future launch plans to avoid detrimental congestion of LEO.

Source-Sink, Launch Rate, System Dynamics, Kessler Syndrome, Debris
Evolutionary Model

##  1 Introduction

###  1.1 Motivation

The unprecedented launch rate of satellites into LEO may have severe
consequences on the stability of the orbital environment for decades to come.
In the $200-900$ km altitude range alone, the number of satellites launched
per year has increased over the last decade as shown in Figure 1. Companies
such as Amazon and SpaceX, have announced plans to launch constellations of
thousands of satellites into LEO over the next decade [1] with hundreds of
satellites already launched. A rising launch rate will increase orbital
congestion which in turn raises the chance of debris-generating collisions. A
growing debris population can have catastrophic consequences on space
missions. It is important to study how increased launch activities affect the
evolution of the orbital environment and the production of debris, as well as
how increased debris populations affect the stability of LEO.

Figure 1: Number of objects launched to $200-900$ km altitudes over the past
decade [2].

###  1.2 Literature Review

Various analytic models of the orbital environment have been proposed in the
literature that make use of differential equations to represent the evolution
of the number of objects in space and their interactions. Some of these models
have also been used to test the LEO environment’s sensitivity to run-away
debris growth. This run-away debris growth is known as Kessler syndrome
wherein the congestion of the orbital environment is large enough to cause a
chain reaction of debris generation. To study such debris growth caused by
collisions, Kessler and Cour-Palais developed a source-sink model to predict
detrimental debris population growth [3]. Furthermore, Talent [4] used one ODE
to represent the total number of objects in space and studied various
evolutionary cases for different launch rates to look for catastrophic
behaviour. Zhang et al. [5] developed a model using partial differential
equations and solved the equations numerically to study the long-term
evolution of the space debris environment. A dynamical systems analysis was
conducted by Drmola and Hubik [6] in which three different classes of debris
were used to study various debris accumulation scenarios and whether they lead
to Kessler syndrome.

The MOCAT model was inspired the three population model JASON given in
reference [7]. The MOCAT-3 model has been used to calculate the intrinsic
capacity of LEO in reference [8]. The model was extended to included a
differentiation between slotted and unslotted satellites in reference [9].
Here, a brief overview of the MOCAT-3 model is given, and various numerical
analysis are conducted using the model for the study of launch rates and
debris generation.

###  1.3 Paper objectives

The objective of this paper is to study the evolution of the LEO environment
for various launch rate distributions and to determine the stability of the
environment in response to perturbations in the debris population. The intent
of finding the stability boundary is to determine how large the debris
population can grow before run-away debris growth, referred to as Kessler
syndrome, occurs.

###  1.4 Paper structure

The remainder of the paper is arranged as follows: Section 2 gives an overview
of the source-sink model and methods used for conducting the stability
analysis; Section 3 introduces various launch rate distribution cases; Section
4 presents the equilibrium solutions and the results of perturbations to these
solutions; Section 5 shows the analysis of the instability threshold where
Kessler syndrome occurs; and finally Section 6 remarks on the meaning of the
analysis and draws conclusions about the future of the orbital environment.

##  2 Methods

###  2.1 MOCAT-3 Model

The MOCAT-3 Model was developed in reference [8] wherein a detailed
description of each parameter is given. Here we provide a brief overview of
the model and the key equations describing the evolution of each species.

MOCAT-3 is a probabilistic source-sink model with three species: active
satellites (S), derelict satellites (D) and debris (N). The orbital
environment within the altitude range of $200-900$ km is divided into $20$
spherical orbital shells with a shell thickness of $35$ km, represented by the
variable $d$. The evolution of each species is represented by a set of
differential equations per shell $\\{\dot{S}(h)$, $\dot{D}(h)$,
$\dot{N}(h)\\}$, where $h$ is a value from $1$ to $20$ indicating the shell
number. Shell with $h=1$ is the lowest shell for the altitude range $200-235$
km and shell with $h=20$ is the highest shell for the altitude range $865-900$
km. We assume each object has a near-circular orbit. The launch rate per year
is represented by $\lambda(h)$ and only appears as a source in $\dot{S}$. New
active satellites appear instantly in their orbital shell $h$ and do not cross
through lower shells. Dropping the explicit dependence on shell number, the
set of coupled ordinary differential equations representing the evolution of
each species is given by equations (1), (2), (3).

$\dot{S}=\lambda-S/\Delta
t-\phi_{SN}(\delta+\alpha)NS-\phi_{SD}(\delta+\alpha)DS-\alpha_{a}\phi_{SS}S^{2}\quad$
(1) $\dot{D}=\frac{(1-P)S}{\Delta t}+\phi_{SD}\delta DS+\phi_{SN}\delta
NS-\phi_{DN}ND-\phi_{DD}D^{2}+\frac{D_{+}v_{D+}}{d}-\frac{Dv_{D}}{d}$ (2)
$\dot{N}=K0_{SN}\phi_{SN}\alpha NS+K0_{SD}\phi_{SD}\alpha
DS+K0_{DN}\phi_{DN}ND+K0_{DD}\phi_{DD}D^{2}+\alpha_{a}K0_{SS}\phi_{SS}S^{2}+K0_{NN}\phi_{NN}N^{2}+\frac{N_{+}v_{N+}}{d}-\frac{Nv_{N}}{d}$
(3)

Here $D_{+}$ and $N_{+}$ refer to the populations of $D$ and $N$ in the shell
directly above the current shell, namely $D_{+}=D(h+1)$ and $N_{+}=N(h+1)$.
For the highest shells, we assume these parameters represent the current
shell: $D_{+}=D(h)$ and $N_{+}=N(h)$, for reasons given in section 2.2. In
general, once an object is in orbit it can only flow into lower altitude
shells and not into upper shells. The de-orbiting of objects from higher
shells to lower shells is dictated by a static exponential model for the
atmospheric density described in [8]. Only derelict and debris objects are
assumed to de-orbit from atmospheric drag effects as active satellites are
assumed to have station-keeping capabilities that counter-act these drag
effects. These drag effects are accounted for in the terms containing the
variable $v$ which represent the change in the semi-major axis. Active
satellites $S$ can become derelict $D$ or debris $N$ through collisions but no
species can become an active satellite $S$ for which the only source is
$\lambda$. Furthermore, active satellites directly exit the environment at a
rate of $1/\Delta t$ with a success probability of $P$, the rest becoming
derelict satellites. The lifetime of each active satellite is taken as $\Delta
t=5$, whereas the probability of successful post-mission disposal is $P=0.95$.
The number of fragments created during collisions between the species is
determined by the NASA standard break-up model [10] which determines the
values of $K0$ and $\phi$. These variables are specified for each type of
collision between all species combinations. For the collision model, the
average mass, area, and diameter values used for each species was taken from
reference [11], and are shown in Table 1. The variables $\delta$, $\alpha$,
$\alpha_{a}$ set the proportionality of collisions that become debris objects.
Specifically, $\delta=10$ gives the ratio of collisions that produce disabling
versus lethal debris, $\alpha=0.2$ is the fraction of derelict and debris
objects that an active satellite fails to avoid, and $\alpha_{a}=0.01$ is the
fraction of active satellites that another active satellite fails to avoid.
This section has summarized the parameters of the MOCAT-3 model which was used
for the numerical analysis given in the rest of the paper.

| Active | Derelict | Debris
---|---|---|---
Mass (kg) | 223 | 223 | 0.640
Area ($m^{2}$) | 1.741 | 1.741 | 0.020
Diameter (m) | 1.490 | 1.490 | 0.180
Table 1: Physical characteristics of each species.

###  2.2 Equilibrium Solutions and Stability

The equilibrium points for the set of coupled differential equations (1), (2),
(3) were solved for each shell by finding the population of each type of
species $\\{S_{eq}(h),D_{eq}(h),N_{eq}(h)\\}$ for which the differential
equations equal zero:

$\dot{S}=0,\hskip 4.0pt\dot{D}=0,\hskip 4.0pt\dot{N}=0.$

For this set of values $\\{S_{eq}(h),D_{eq}(h),N_{eq}(h)\\}$, the sources and
sinks of the environment balance each other and the system is in equilibrium.
A change in the launch rate generates a new set of equilibrium solutions
because the launch rate is a major source in the active satellites population.
Since each differential equation (1), (2), (3) has degree $2$, the set of $3$
coupled equations has $2^{3}=8$ equilibrium solutions per shell. We eliminate
solutions that are purely imaginary or that contain a real negative part, as
these are considered non-physical solutions for the species populations. For
the launch rate cases we studied, we found that each shell has two sets of
positive, real-valued equilibrium solutions $\\{S_{1}(h),D_{1}(h),N_{1}(h)\\}$
and $\\{S_{2}(h),D_{2}(h),N_{2}(h)\\}$, with one solution set having a larger
number of active satellites than the other: $S_{1}>S_{2}$. We used this
solution set $\\{S_{1}(h),D_{1}(h),N_{1}(h)\\}$ as the influx populations to
the next lower shell $\\{S_{1}(h-1),D_{1}(h-1),N_{1}(h-1)\\}$. For the highest
shell, we assumed the influx of objects from higher altitudes was equaled to
the outflow of objects from that shell. This assumption may differ from
reality as there are many objects located above $900$ km but it is difficult
to measure how many of these objects would de-orbit due to atmospheric drag at
such high altitudes per year. By finding the equilibrium solutions
$\\{S,D,N\\}$ per shell starting with the highest altitude shell and ending
with the lowest altitude shell, we guarantee the equilibrium of the entire
orbital environment in the 200km-900km altitude range because each shell
depends only on the species population within that shell and directly above
it. Our model does not contain any flows from lower shells into higher shells
so the lower shells do not affect the equilibrium of the higher shells.

For a given launch rate and set of initial conditions
$\\{S_{i},D_{i},N_{i}\\}$, we integrate the differential equations (1), (2),
(3) with respect to time to find the number of years required for the source-
sink model to reach equilibrium. The initial conditions we used was based on
Two Line Element data from space-track.org 111The TLE catalog was downloaded
from space-track.org, accessed on August 18th, 2022.. We used the process of
reference [9] that classified each object as an active satellite, derelict
satellite, or as debris according to its mass, diameter and area. We note that
unlike reference [9], we do not distinguish between slotted and unslotted
satellites, thus we combined these two populations for the active satellite
population. The initial populations of each species are displayed per altitude
shell in Figure 2.

Figure 2: Population of each species from TLE data as of August 2022.

###  2.3 Basin of Attraction

_Basin of Attraction Definition:_ “The set of points in the space of system
variables such that initial conditions chosen in this set dynamically evolve
to a particular attractor” [12].

In our model, an attractor is a stable equilibrium point. Once the source-sink
model is at a steady equilibrium state, we studied how the model reacts to
perturbations in $\\{S,D,N\\}$ from equilibrium. In particular, we studied how
a sudden increase in debris $N$ alters the orbital environment by looking at
how the system evolved after the perturbation and whether or not it tended
back toward equilibrium. We found the basin of attraction about the stable
equilibrium point for various initial conditions of debris $N$. One of the
reasons we decided to focus on analyzing the stability of the orbital
environment with respect to changes in the species of debris rather than
changes in active or derelict satellites, is the exact amount of debris
currently in LEO is unknown. Another reason is we wanted to study the
instability threshold for debris creation known as ‘Kessler syndrome’ in which
the amount of debris in the orbital environment is numerous enough that it
continuous to generate more and more debris, creating a chain reaction. To
visualize the basin of attraction, we plotted phase space diagrams.

##  3 Launch Rate Distributions

_Launch Rate Distribution Definition:_ The number of active satellites
launched into orbit per year per altitude shell.

We used various launch rate distributions to study the stability of the LEO
environment. The two types of launch rate distributions that we studied are
_static_ and _dynamic_ launch rates. A static launch rate represents a
constant influx of active satellites per altitude shell per year for a given
number of years. A dynamic launch rate represents a variable launch rate per
altitude shell per year. For each launch rate, a unique set of equilibrium
solutions is found. Thus, a static launch rate is used to study the behaviour
of the system of equations with respect to the equilibrium solutions. However,
a varying launch rate represents a more realistic scenario since the number of
satellites launched per year has changed drastically over the past few
decades. A dynamic launch rate was studied as a separate case wherein the
equilibrium solutions changed with variations in the launch rate.

###  3.1 Static Launch Rate

We studied two cases of static launch rates. For the first case we used the
maximum number of satellites launched in one year per altitude shell over the
past ten years. The second case, we used ‘As Received’ filings database from
the International Communications Union (ITU).

####  3.1.1 _Case 1: Past Launch Rates_

We used the maximum number of satellites launched within one year over the
past ten years for each altitude shell. We use the maximum historic launch
rate per shell instead of using the launch rate from a specific year because
particular years have a low number of launches to certain shells that doesn’t
represent the general behaviour for launch cadence. This launch rate
distribution allowed us to analyze the stability of the current LEO
environment to see if current launch activities are sustainable or if we are
already in danger of run-away debris growth. The number of satellites launched
into each altitude shell was taken from the Union of Concerned Scientists
database [2], and is displayed in Figure 3.

Figure 3: Maximum launch rate over the last ten years per altitude shell [2].

####  3.1.2 _Case 2: ITU Filings_

The ‘As Received’ ITU filings database [13], is a list of satellite notices
filed with the ITU that have not yet been reviewed or published by the ITU. It
should be noted that the ITU states this database is not regulated. However,
this database allows for some forecasting of satellite launches over the next
few years. Each filing includes the altitude and number of satellites that an
organization intends to launch. We have filtered through this database to
eliminate duplicate filings made by the same organization. We used the average
of the apogee and perigee altitude to bin the satellites into altitude shells.
The number of satellites forecast to be launched into each shell is shown in
Figure 4.

Figure 4: ‘As Received’ ITU filings of satellite notices [13].

These filings are valid for several years and the deployment of a satellite or
constellation of satellites into orbit can take over a year. Thus, in using
ITU filings to estimate a launch rate per year, we divided the total number of
satellites in each shell by a number of years $n$.

###  3.2 Dynamic Launch Rate

We studied a dynamic launch rate as a third launch rate distribution. As shown
in Figure 1, the number of objects launched into orbit each year has not
remained constant. By using a dynamic launch rate, we can represent such a
change in launch rate per year.

####  3.2.1 Case 3: Varying Launch Rate per Year

We followed a similar approach to [4] in modeling a dynamic launch rate. We
took the launch rate displayed in Figure 3 as the base case and then increased
this launch rate by 0%,1%,3%,5%, and 7% each year for $50$ years. Then we set
the launch rate to be constant at the rate calculated at the end of the $50$
years. We let the environment evolve for another $800$ years at this constant
launch rate. The total number of objects per year for each incremental launch
rate is displayed in Figure 5. The total number of launched satellites and
species populations at the end of the 800 years for each incremental launch
rate is given in Table 2. As can be seen in the table, the various percentage
increases in launch rate for the first 50 years produces drastically different
total number of objects in the environment at the end of the simulation. For
example, a 7% increase in launch rate produces a total number of objects at
the end of the 800 year simulation that is about two orders of magnitude
larger than the total produced by a 1% increase in launch rate. This numerical
analysis shows how a few percent difference in increasing launch rate per year
creates large differences in the population of each species when propagated
over time.

Figure 5: Total number of objects in orbit for launch rates growing by 0-7% per year for 50 years and then remaining constant for 800 years. % Increase |  Launch $\lambda$ |  Active S |  Derelict D |  Debris N |  Total
---|---|---|---|---|---
0% |  1965 |  9820 |  1383 |  1402 |  12606
1% |  3232 |  16146 |  1984 |  1743 |  19873
3% |  8614 |  42976 |  4621 |  4317 |  51914
5% |  22533 |  111830 |  12387 |  22061 |  146270
7% |  57883 |  270980 |  60571 |  741220 |  1072800
Table 2: Population of each species and total number of objects in orbit at
the end of a period of constant growth rate in launch as shown in Figure 5.

##  4 Results and Discussion

###  4.1 Equilibrium Solutions

For each launch rate, we determined the set of equilibrium solutions. We
computed the equilibrium solutions per shell for the ‘business as usual’ case
3.1.1; the results for which are presented in Figure 6. Since equilibrium
solutions exist, we can conclude that our current launch activity is
sustainable for the $200-900$ km altitude range and that run-away debris
growth will not occur if launch rates remain at these levels.

Figure 6: Equilibrium solutions per species for the constant launch rate case
3.1.1.

We also computed the equilibrium solutions for case 3.1.2 where
$\lambda_{itu}/n$ with $n=7$ years and found that no positive, real-valued
solutions existed. This may seem alarming but realistically we should not
assume that all the satellite notices filed will be used. Only a fraction of
the satellite notice filings will actually be launched. We increased $n$ to
$n=21$, and found a set of equilibrium points for each shell as shown in
Figure 7.

Figure 7: Equilibrium solutions per species for a constant launch rate
proportional to ITU filings given in case 3.1.2.

For the dynamic launch case 3.2.1, positive, real-valued equilibrium solutions
could only be found for the 1% increase in launch rate per year. A summary of
total populations of each species at equilibrium for each launch case is given
in Table 3.

Rate | Total Launch $\lambda$ | Total Active S | Total Derelict D | Total Debris N | Total All
---|---|---|---|---|---
Case 1 | 1965 | 9820 | 1383 | 1402 | 12606
Case 2 | 56410 | 278847 | 22704 | 157870 | 459421
Case 3: 1% | 3232 | 16145 | 2292 | 4281 | 22719
Table 3: Total population of each species at equilibrium for various launch
rate distributions.

Overall, the amount of active satellites $S$ at equilibrium is proportional to
the number of satellites launched. The amount of debris $N$ and derelict
satellites $D$, however, vary with respect to the amount of active satellites,
with higher shells acquiring a larger number of derelict and debris objects
because the sink caused by atmospheric drag removes less objects per year at
higher altitudes due to the lower atmospheric density.

###  4.2 Stability Analysis

The stability of the equilibrium solutions found for Case 1 (3.1.1), Case 2
(3.1.2) and Case 3 (3.2.1) at 1% increment only, was determined by computing
the eigenvalues of the solutions. All eigenvalues were found to be negative
indicating all of these equilibrium solutions are stable. As an example, Table
4 displays the eigenvalues for Case 2 launch rate of section 3.1.2.

Active S | Derelict D | Debris N
---|---|---
-0.200 | -1.028 | -0.212
-221.167 | -0.591 | -0.211
-885.271 | -0.432 | -0.211
-345.598 | -0.350 | -0.200
-149.086 | -0.253 | -0.200
-86.341 | -0.119 | -0.209
-37.246 | -0.089 | -0.208
-70.472 | -0.077 | -0.208
-36.780 | -0.056 | -0.205
-19.147 | -0.043 | -0.200
-17.605 | -0.001 | -0.203
-10.008 | -0.003 | -0.205
-9.180 | -0.004 | -0.201
-5.522 | -0.007 | -0.201
-4.784 | -0.017 | -0.201
-3.129 | -0.024 | -0.204
-2.350 | -0.026 | -0.203
-1.306 | -0.027 | -0.203
-1.776 | -0.027 | -0.201
-0.777 | -0.147 | -0.202
Table 4: The eigenvalues of each population for the equilibrium solutions
displayed in Figure 7.

Additionally, we used phase portraits to depict the phase space about the
stable equilibrium point per shell. The phase portraits for the equilibrium
points of the Case 1 launch rate, shown in Figure 6, are displayed for two
altitude shells in Figure 8.

Figure 8: Phase portraits about the stable equilibrium state at different
altitudes.

We analyzed how many years were needed for the orbital environment to settle
into its equilibrium state for the launch rate $\lambda_{itu}/21$ assuming the
initial conditions are given by Figure 2. We integrated the set of
differential equations using these initial conditions. The results are shown
for active satellites, derelict satellites, and debris in Figures 9, 10, 11,
respectively. The population of active satellites settles into equilibrium
across all shells within 10 years. The populations of derelict satellites and
debris requires greater than 50 years to reach equilibrium for particular
altitude shells.

Figure 9: Number of active satellites in orbit over time for the launch rate
given in Section 3.1.1. Figure 10: Number of derelict satellites in orbit over
time for the launch rate given in Section 3.1.1. Figure 11: Amount of debris
in orbit over time for the launch rate given in Section 3.1.1.

###  4.3 Perturbing Launch Rate for One Year

We studied how the orbital environment would react to a one-time drastic
increase in launch. This perturbation in launch rate allows for all ITU
filings displayed in 4 to be launched in one year. The utility of this
approach is we can study how the orbital environment reacts to a launch rate
for which no equilibrium solutions exist as stated in section 2.2. We started
with initial conditions given by Figure 2 and a launch rate of
$\lambda_{itu}/21$ as shown in Figure 7. We allowed the system of equations to
evolve for 20 years at this constant launch rate and then we increased the
launch by a factor of twenty for one year. These two launch rates are
displayed in Figure 12.

Figure 12: Launch forecast proportional to ITU ’As received’ filings.

After this one year increase in launch activity, we decreased the launch rate
back to the original rate of $\lambda_{itu}/21$ and allowed the system to
evolve for another 20 years. The results are shown for each species in Figures
13, 14, 15. Overall, the system evolved back toward the stable equilibrium
solution given in Figure 7 within the 20 years following the perturbation in
launch rate. A change in the launch rate in effect perturbs the population of
each species away from equilibrium.

Figure 13: Number of active satellites over time with a one time increase in
launch at 20 years. Figure 14: Number of derelict satellites over time with a
one time increase in launch at 20 years.  Figure 15: Number of debris objects
over time with a one time increase in launch at 20 years.

###  4.4 Perturbations of Equilibrium Solutions

We studied how the orbital environment reacts to the event of a sudden
increase in debris. We analyzed two cases: the first case depicts the effect
of an uniform increase in debris across all shells, and the second case
depicts the effect of debris increase in one shell and compares a perturbation
in debris at a high altitude vs. at a low altitude.

####  4.4.1 Equal Perturbation Across all Shells

Starting with the system at equilibrium depicted in Figure 7 for the constant
launch rate given by case 3.1.2, we perturbed the amount of debris in each
shell by $10,000$ objects at $t=20$ years and allowed the system to evolve for
200 years. The change in the amount of debris is shown in Figure 16. Table 5
summarizes these results by displaying the total amount of each species before
the increase in debris, at the time of the event, and 200 years after. From
Figure 16, we note that such an event creates a minimal effect in lower
altitude shells with each species returning back to its equilibrium state
within a few years. However, for higher altitude shells the scenario is
drastically different with the system remaining out of equilibrium for at
least $200$ years. This analysis shows how an ‘explosion’ type of event that
produces a large amount of debris across all shells, greatly affects the
amount of debris present at higher altitudes for many years following the
event.

Species |  Initial |  At Event |  After 200 Years |  $\Delta$
---|---|---|---|---
S |  292830 |  - |  292629 |  -201
D |  21290 |  - |  22597 |  +1308
N |  107712 |  307712 |  145379 |  +37667
Table 5: Population of each species before and after a sudden increase in
debris across all shells. Figure 16: Debris population over time with an
impulsive increase in debris by 10,000 fragments at 20 years.

####  4.4.2 Impulsive Debris Perturbations in High Shell vs Low Shell

From the equilibrium values displayed in Figure 7, the amount of debris was
increased by 10,000 objects in the second-highest shell (830km-865km) and in
the shell with altitude range (410km-445km). We chose the second highest
altitude shell rather than the highest shell since it contained significantly
more active satellites. For the lower altitude shell, we chose the shell that
contained the greatest amount of active satellites overall since this shell is
the most sensitive to collisions between debris and active satellites. We set
the orbital environment to equilibrium and then we added a perturbation in
debris at $t=20$ years. After this perturbation, we allowed the system to
evolve for 200 years. The results are displayed in Figure 17. Increasing the
amount of debris by 10,000 objects in a high altitude shell has a much more
significant impact on the overall LEO environment than increasing the debris
in a lower altitude shell. In the left side of Figure 17 we see that the
environment quickly returns to its near-equilibrium state after the
perturbation occurs, whereas the environment does not recover to equilibrium
if such a perturbation occurs at a high altitude shell as shown in the right
side of the figure. Rather than returning to its equilibrium state, a large
perturbation of the debris population at a high altitude causes the population
of debris to keep growing across multiple shells over the course of 200 years.
It could be that for a longer period of time $t>200$ years, the system will
return to its equilibrium state or it could be the case that the population of
debris has reached an amount that causes collisions to continuously occur and
debris to grow without end. We studied this behaviour in more detail in the
next section 5.

Figure 17: A comparison of the evolution of the debris population after a
sudden increase in debris at $t=20$ years in two different shells.

##  5 Instability Threshold: Kessler Syndrome

Using the launch rate given in case 3.1.2, we calculated the instability
threshold as the maximum perturbation in debris away from equilibrium for
which the population of debris continues to increase without bound over 1,000
years. The two types of perturbations we used were perturbing debris in all
shells simultaneously as done in section 4.4.1 and perturbing debris in each
shell individually similar to the approach used in section 4.4.2. The system
does not necessarily need to return to its equilibrium solution within 1,000
years of the perturbation but for the perturbation to be considered apart of
the stable region, the amount of debris and derelict satellites need to be
decreasing at the end of the 1,000 years. In other words, at $t=1000$ years
the set of $\\{S,D,N\\}$ must satisfy:

$\hskip 4.0pt\dot{D}\leq 0,\hskip 4.0pt\dot{N}\leq 0$

for all shells. The number of active satellites can be decreasing or
increasing at the end time. In this way, we calculated the threshold at which
run-away debris growth occurs, referred to as Kessler syndrome.

###  5.1 Perturbing All Shells Simultaneously

To the nearest thousandth, the largest perturbation to the debris population
for which the system reverted toward the equilibrium state after 1,000 years
was found to be 29,000 debris objects. Perturbing the debris population by
30,000 objects was found to cause runaway debris growth. These two cases are
displayed in Figure 18. In the left of Figure 18 it is clear that the
perturbation in debris causes debris growth for about 400 years but the system
begins to return to equilibrium as it evolves for the remaining 600 years.
This is not true for a perturbation of 30,000 debris objects as displayed in
the right side of the Figure 18, wherein the population of debris continues to
grow reaching a couple quintillion before the integration fails at $t=600$
years. Such run-away debris growth displays Kessler syndrome as the orbital
environment is unstable and continues to diverge away from its equilibrium
state. Through collisions with active satellites, a perturbation in debris
also causes a change in the population of derelict satellites. The evolution
of the derelict population for a perturbation of 29,000 and 30,000 debris
objects is displayed in Figure 19. Overall, the runaway debris growth occurred
in higher altitude shells, but we are also interested in the instability
threshold of lower altitude shells. Thus, rather than simultaneously
perturbing all shells away from equilibrium, we studied the instability
threshold of each altitude shell in the next section.

Figure 18: A comparison of the evolution of the debris population after a
perturbation in debris occurs across all shells. The stable regime is
displayed on the left and the unstable regime is displayed on the right.
Figure 19: A comparison of evolution of the derelict population after a
perturbation in debris across all shells.

###  5.2 Perturbing Shells Individually

By perturbing the amount of debris in each shell individually, we were able to
analyze how sensitive each altitude shell is to a sudden increase in debris.
We calculated the maximum perturbation in the debris population away from
equilibrium to the nearest thousandth for which the system evolved back toward
equilibrium within 1000 years. The instability threshold hence exists at this
maximum perturbation amount. The results are presented in Table 6. Altitudes
below 410 km are not included in the table because these shells could
withstand a perturbation of debris equal to $10^{8}$ objects. We conclude that
the stability of the orbital environment at lower altitude shells is much more
resilient to perturbations in debris than higher altitude shells. This result
concurs with the result of section 4.4.2. The reasoning for this behaviour is
the sink of the model, namely atmospheric drag, is much greater at lower
altitude shells, which removes debris from the environment preventing
collisions with active and derelict satellites that would create more debris.
We would like to note that this analysis was done for a particular launch rate
taken as a fraction of the ITU filings $\lambda_{itu}/21$ as shown in Figure
7. Debris creation is directly affected by the launch rate since launch
activity is the source of active satellites per shell, and a higher density of
active satellites per shell creates a greater likelihood of collision with
debris. Thus, different evolution of the debris population would occur for a
different launch rate.

Altitude Shell (km) |  Debris at Equilibrium |  Max. Perturbation in Debris
---|---|---
410-445 |  2462 |  53442000
445-480 |  1445 |  28478000
480-515 |  1145 |  17331000
515-550 |  1953 |  9861000
550-585 |  3077 |  5775000
585-620 |  5085 |  3348000
620-655 |  7983 |  2001000
655-690 |  12046 |  1216000
690-725 |  17185 |  610000
725-760 |  14725 |  414000
760-795 |  18644 |  263000
795-830 |  18621 |  91000
830-865 |  703 |  71000
865-900 |  25 |  63000
Table 6: Maximum perturbation in debris per shell before Kessler Syndrome
occurs.

##  6 Remarks and Conclusions

Given a launch rate distribution that is based on historic launch activities
(3.1.1), the current orbital environment will evolve to a stable equilibrium
state. In such a state, the sources of the model, namely due to launch and
collisions, balance the sinks of the model, namely post-mission disposal and
atmospheric drag. Thus if launch activities remain at current levels, Kessler
syndrome will not occur in the $200-900$ km altitude range of the orbital
environment, given the assumptions of our model. This is also true for an
increased but constant launch rate (3.1.2) for which a stable equilibrium
state exists. We note that the evolution of the environment from the current
populations of active, derelict, and debris objects to this equilibrium state
would take decades. Given a dynamic launch rate distribution, an ever-
increasing launch rate entails the system will not reach an equilibrium state.
However, if the launch rate becomes constant after a period of continuous
growth then the system may evolve towards the equilibrium state if the growth
rate was small enough. In our analysis, only a growth rate of 1% in launch
rate per year over 50 years lead to a stable equilibrium state. Larger growth
rates in launch rate entailed no equilibrium state would be reached with each
species population ever-increasing.

Perturbations in the debris population away from equilibrium showed how
sensitive the environment is to an increase in debris. In general,
perturbations in the debris populations in higher altitude shells had more
drastic consequences than perturbations in lower altitude shells due to
increased atmospheric drag forces existing at lower altitudes. Run-away debris
growth is more common at high altitudes, with Kessler syndrome resulting from
significantly smaller perturbations in debris than at lower altitudes. Thus a
debris-generating event occurring at a high altitude is more dangerous than
one occurring at a low altitude, since at higher altitudes such an event can
trigger Kessler syndrome to occur.

## Acknowledgments

The authors would like to thank Thomas Roberts for sharing processed ITU
filing data. The authors wish to acknowledge the support of this work by the
Defense Advanced Research Projects Agency (Grant N66001-20-1-4028). The
content of the information does not necessarily reflect the position or the
policy of the Government. No official endorsement should be inferred.
Distribution statement A: Approved for public release; distribution is
unlimited.

## References

  * [1] C. Davenport, “Thousands more satellites could soon be launched into space. can the federal government keep up?,” Jul 2020.
  * [2] UCS-Satellite-Database-5-1-2022. http://www.ucsusa.org/satellite_database.
  * [3] D. J. Kessler and B. G. Cour-Palais, “Collision frequency of artificial satellites: The creation of a debris belt,” Journal of Geophysical Research: Space Physics, vol. 83, no. A6, pp. 2637–2646, 1978.
  * [4] D. L. Talent, “Analytic model for orbital debris environmental management,” Journal of Spacecraft and Rockets, vol. 29, no. 4, pp. 508–513, 1992.
  * [5] B. Zhang, Z. Wang, and Y. Zhang, “Discrete evolution model based on mean spatial density for space debris environment,” Astrophysics and Space Science, vol. 364, 04 2019.
  * [6] J. Drmola and T. Hubik, “Kessler syndrome: System dynamics model,” Space Policy, vol. 44-45, pp. 29–39, 2018.
  * [7] H. G. Lewis, G. G. Swinerd, R. J. Newland, and A. Saunders, “The fast debris evolution model,” Advances in Space Research, vol. 44, pp. 568–578, 9 2009\.
  * [8] A. D’Ambrosio, M. Lifson, and R. Linares, “The capacity of low earth orbit computed using source-sink modeling,” 2022.
  * [9] M. Lifson, A. D’Ambrosio, D. Arnas, and R. Linares, “How many satellites can we fit in low earth orbit?: Capacity integrating risk-based and intrinsic methods,” Astrodynamics Specialist Conference, Preprint.
  * [10] P. Krisko, “Proper implementation of the 1998 nasa breakup model,” Orbital Debris Quarterly News, vol. 15, no. 4, pp. 1–10, 2011.
  * [11] G. L. Somma, Adaptive remediation of the space debris environment using feedback control.  PhD thesis, University of Southampton, 2019.
  * [12] M. Crisan, “Convergence towards a dynamic theory of linguistics and semantics,” in Convergence and Hybrid Information Technologies (M. Crisan, ed.), ch. 3, Rijeka: IntechOpen, 2010.
  * [13] ITU e-Submission of Satellite Network Filings ‘As Received’. https://www.itu.int/ITU-R/space/asreceived/Publication/AsReceived.


$\displaystyle\deg_{\tilde{m}_{0}}(\tilde{m}^{\prime}):=\sum_{i=1}^{r}a_{i}.$
(0.7.2)

Then, by (a), we have

$\displaystyle\deg_{\tilde{m}_{0}}(\tilde{m}_{0})=0<\deg_{\tilde{m}_{0}}(\tilde{m}_{1})<\cdots<\deg_{\tilde{m}_{0}}(\tilde{m}_{J})=\deg_{\tilde{m}_{0}}(\tilde{m}_{\gamma}).$
(0.7.3)

(d). For the segment $\gamma_{j}$ in the condition (5), suppose that
$\gamma_{j}$ crosses walls $[\mathfrak{d}_{\lambda},f_{{\lambda}}]_{n_{j}}$
$(\lambda\in\Lambda_{j}$) at $t_{j}$. Note that $\gamma^{\prime}=-m_{j-1}$
just before $\gamma$ crosses these walls. Then, the intersection signs in
(0.2.6) are given by

$\displaystyle\epsilon_{j}=\begin{cases}1&\langle n_{j},m_{j-1}\rangle>0,\\\
-1&\langle n_{j},m_{j-1}\rangle<0.\\\ \end{cases}$ (0.7.4)

Thus, $\epsilon_{j}\langle n_{j},m_{j-1}\rangle=|\langle
n_{j},m_{j-1}\rangle|=|\langle n_{j},\tilde{m}_{j-1}\rangle_{1}|$. Therefore,
by Definition 0.2.6, we have

$\displaystyle\mathfrak{p}_{\gamma_{j},\mathfrak{D}}(c_{j-1}x^{\tilde{m}_{j-1}})=c_{j-1}x^{\tilde{m}_{j-1}}\prod_{\lambda\in\Lambda_{j}}f_{{\lambda}}^{|\langle\delta(n_{j})n_{j},\tilde{m}_{j-1}\rangle_{1}|}.$
(0.7.5)

It is crucial that there is _no division_ by $f_{{\lambda}}$ in this
expression for the forthcoming positivity of theta functions.

######  Definition 0.7.2 (General position).

We say that $Q\in M_{\mathbb{R}}\setminus\mathrm{Supp}(\mathfrak{D})$ is _in
general position_ if for each $\ell>0$ there is some neighborhood $U_{\ell}$
of $Q$ such that, for the reduction $\mathfrak{D}_{\ell}$, any broken line for
$Q^{\prime}\in U_{\ell}$ converges to a broken line for $Q$ in the limit
$Q^{\prime}\rightarrow Q$. (Namely, no broken line crosses
$\mathrm{Sing}(\mathfrak{D}_{\ell})$ in the limit.)

Below we always and implicitly assume that $Q$ is in general position.

######  Definition 0.7.3 (Theta function).

Under the same assumption and notations in Definition 0.7.1, the _theta
function $\vartheta_{Q,\tilde{m}_{0}}$_ for $\tilde{m}_{0}$ with endpoint $Q$
is defined by

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}:=\sum_{\gamma\in
B(Q,\tilde{m}_{0})}c_{\gamma}x^{\tilde{m}_{\gamma}},$ (0.7.6)

where $B(Q,\tilde{m}_{0})$ is the set of all broken lines for $\tilde{m}_{0}$
with endpoint $Q$. We also set

$\displaystyle\vartheta_{Q,0}=1.$ (0.7.7)

The function $\vartheta_{Q,\tilde{m}_{0}}$ only depends on the equivalence
class of $\mathfrak{D}$, though the condition for $Q$ to be in general
position depends on the choice of a representative of $\mathfrak{D}$.

######  Lemma 0.7.4 ([GHKK18, Proposition 3.4]).

For any scattering diagram $\mathfrak{D}$, we have

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}\in
x^{\tilde{m}_{0}}\mathbbm{k}[[P_{1}]].$ (0.7.8)

###### Proof.

By Remark (b), each term $c_{\gamma}x^{\tilde{m}_{\gamma}}$ belongs to
$x^{\tilde{m}_{0}}\mathbbm{k}[[P_{1}]]$. Thus, it is enough to show that, for
a given $\tilde{m}^{\prime}\in\tilde{m}_{0}+P_{1}$, there are only finitely
many $\gamma$ such that $\tilde{m}_{\gamma}=\tilde{m}^{\prime}$. We note the
following facts:

  * •

By the constraint (0.7.3), such $\gamma$ bends at most
$\deg_{\tilde{m}_{0}}(\tilde{m}^{\prime})$ times.

  * •

The walls of $\mathfrak{D}$ that contribute to $\tilde{m}^{\prime}$ belongs to
the reduction $\mathfrak{D}_{\ell}$ at
$\ell=\deg_{\tilde{m}_{0}}(\tilde{m}^{\prime})$, which has only finitely many
walls.

  * •

At each bending there are only finitely many possibilities of bending.

Therefore, there are only finitely many possibilities of such $\gamma$. ∎

######  Proposition 0.7.5 ([GHKK18, Proposition 3.8]).

Let $\mathfrak{D}$ be any scattering diagram for a given seed $\mathfrak{s}$,
and let $\mathcal{C}^{+}_{\mathfrak{s}}$ be the positive orthant in (0.2.40).
For $\tilde{m}_{0}\in\tilde{M}^{\circ}$, suppose that
$m_{0}\in\mathcal{C}^{+}_{\mathfrak{s}}\cap M^{\circ}$. Also, suppose that
$Q\in\mathrm{Int}(\mathcal{C}^{+}_{\mathfrak{s}})$. Then, we have

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}=x^{\tilde{m}_{0}}.$ (0.7.9)

###### Proof.

We may assume that $\tilde{m}_{0}\neq 0$. It is clear that $Q+\mathbb{R}_{\geq
0}m_{0}$ is a broken line for $\tilde{m}_{0}$ with end point $Q$, and it does
not intersect any walls of $\mathfrak{D}_{\mathfrak{s}}$; therefore, the
associated monomial is $x^{\tilde{m}_{0}}$. We claim that under the condition
there is no broken line $\gamma$ with bending. Suppose that such $\gamma$
exists. We use the notation in Definition 0.7.1. For each break point $t_{j}$
therein, let $n_{j}\in N^{+}$ be a unique (not necessarily primitive) element
such that $\gamma$ intersects walls in $n_{j}^{\perp}$ at $t_{j}$ and that the
velocity of $\gamma$ is shifted by $-p^{*}(n_{j})$. Let $L_{j}$ be the segment
of $\gamma$ for the interval $I_{j}$. Then, we claim the following:

$\displaystyle L_{j}\subset H_{j}^{-}:=\biggl{\\{}z\in
M_{\mathbb{R}}\biggm{|}\biggl{\langle}\sum_{i=1}^{j}n_{i},z\biggr{\rangle}\leq
0\biggr{\\}}\quad(j=1,\dots,J).$ (0.7.10)

We show it by the induction on $j$. This is true for $j=1$, because the
support of the wall crossed by $\gamma$ at $t_{1}$ is in $n_{1}^{\perp}$, and
$\gamma$ is in $\mathcal{C}^{+}_{\mathfrak{s}}$ for $t\rightarrow-\infty$. Now
suppose that (0.7.10) hold up to $j$. Then, we have
$\langle\sum_{i=1}^{j}n_{i},\gamma(t_{j+1})\rangle\leq 0$. Also, $\langle
n_{j+1},\gamma(t_{j+1})\rangle=0$. Thus, $\gamma(t_{j+1})\in H_{j+1}^{-}$.
Moreover, the new velocity vector
$\gamma^{\prime}_{j+1}=-m_{0}-\sum_{i=1}^{j+1}p^{*}(n_{i})$ satisfies the
inequality

$\displaystyle\begin{split}\biggl{\langle}\sum_{i=1}^{j+1}n_{i},\gamma^{\prime}_{j+1}\biggr{\rangle}&=\biggl{\langle}\sum_{i=1}^{j+1}n_{i},-m_{0}\biggr{\rangle}+\biggl{\\{}\sum_{i=1}^{j+1}n_{i},\sum_{i=1}^{j+1}n_{i}\biggr{\\}}\\\
&=\biggl{\langle}\sum_{i=1}^{j+1}n_{i},-m_{0}\biggr{\rangle}\leq
0.\end{split}$ (0.7.11)

Thus, the fact (0.7.10) holds for $j+1$. Since
$\mathrm{Int}(\mathcal{C}^{+}_{\mathfrak{s}})\cap H_{j}^{-}=\emptyset$,
$\gamma$ cannot reach $Q\in\mathrm{Int}(\mathcal{C}^{+}_{\mathfrak{s}})$. ∎

###  0.7.2 Transitivity on theta functions

Now we assume that $\mathfrak{D}$ is consistent.

######  Proposition 0.7.6 ([CPS, Lemma 4.9]).

Let $\mathfrak{D}$ be a consistent scattering diagram. Let
$\tilde{m}_{0}\in\tilde{M}^{\circ}$ and $Q,\,Q^{\prime}\in
M_{\mathbb{R}}\setminus\mathrm{Supp}(\mathfrak{D})$. Then, for any admissible
curve $\gamma$ from $Q$ to $Q^{\prime}$, we have

$\displaystyle\vartheta_{Q^{\prime},\tilde{m}_{0}}=\mathfrak{p}_{\gamma,\mathfrak{D}}(\vartheta_{Q,\tilde{m}_{0}}).$
(0.7.12)

The rest of this subsection is devoted to prove this proposition. We may
assume that $\tilde{m}_{0}\neq 0$.

Let $J$ be the maximal ideal of $\mathbbm{k}[[P_{1}]]$ generated by
$P_{1}\setminus\\{0\\}$.

######  Lemma 0.7.7 ([CPS, Lemma 4.7]).

Let $\mathfrak{D}_{\ell}$ be the reduction of $\mathfrak{D}$ at $\ell$.
Suppose that $Q$ and $Q^{\prime}$ are in the same chamber of
$\mathrm{Supp}(\mathfrak{D}_{\ell})$. Then, we have the equality

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}\equiv\vartheta_{Q^{\prime},\tilde{m}_{0}}\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.13)

###### Proof.

Since we consider $\vartheta_{Q,\tilde{m}_{0}}$ only modulo
$x^{\tilde{m}_{0}}J^{\ell+1}$, one can replace $\mathfrak{D}$ with
$\mathfrak{D}_{\ell}$ to calculate $\vartheta_{Q,\tilde{m}_{0}}$. First,
suppose that $Q$ and $Q^{\prime}$ are connected by a curve $\beta$ in the same
chamber such that any point on $\beta$ is in general position for
$\mathfrak{D}_{\ell}$. Then, each broken line with endpoint $Q$ continuously
deformed to a broken line with endpoint $Q^{\prime}$. Therefore, the theta
function is unchanged.

Suppose that there is no such curve between $Q$ and $Q^{\prime}$. Let us
consider the special case where $Q$ and $Q^{\prime}$ are connected by a curve
$\beta$ in the same chamber such that there is a unique point $Q_{0}$ on
$\beta$ that is not in general position. (The general case reduces to this
case.) Then, there is at least one broken line $\gamma$ with endpoint $Q$ such
that, in the limit $Q\rightarrow Q_{0}$, the broken line intersects with some
boundary $\partial\mathfrak{d}$ or joint
$\mathfrak{j}=\mathfrak{d}_{1}\cap\mathfrak{d}_{2}$ of $\mathfrak{D}_{\ell}$.

(a)$\gamma_{1}$$\gamma_{2}$$\gamma^{\prime}_{1}$$\gamma_{A}$$\gamma_{B}$$\mathfrak{j}$$\textstyle{\bullet}$$\textstyle{\bullet}$$\textstyle{\bullet}$

(b)$\gamma$$\gamma^{\prime}$$\gamma_{A}$$\gamma_{B}$$\mathfrak{c}$$\textstyle{\bullet}$$\textstyle{\bullet}$$\textstyle{\bullet}$

Figure 17: Broken lines in the proof of Lemma 0.7.7.

First, let us consider the case $\mathfrak{j}$. See Figure 17 (a). Suppose
that $\gamma$ crosses the walls that contain the joint $\mathfrak{j}$ with the
incoming attached monomial $c_{\mathrm{in}}x^{\tilde{m}_{\mathrm{in}}}$ and
the outgoing attached monomial $c_{\mathrm{out}}x^{\tilde{m}_{\mathrm{out}}}$.
There may be several broken lines $\gamma_{1}=\gamma$, …, $\gamma_{p}$ with
endpoint $Q$ such that their bending patterns differ from $\gamma$ only inside
these walls. They have the common incoming attached monomial
$c_{\mathrm{in}}x^{\tilde{m}_{\mathrm{in}}}$ and the common exponent
$\tilde{m}_{\mathrm{out}}$ for the outgoing attached monomials. Then, all
these broken lines have a common limit for $Q\rightarrow Q_{0}$. Let
$\gamma_{A}$ be the half of a small loop around $\mathfrak{j}$, which is on
the same side of $\mathfrak{j}$ with $\gamma_{i}$’s as in Figure 17 (a). Then,
the sum of the outgoing attached monomials for $\gamma_{1}$, …, $\gamma_{p}$
is exactly the terms of the monomials in
$\mathfrak{p}_{\gamma_{A},\mathfrak{D}_{\ell}}(c_{\mathrm{in}}x^{\tilde{m}_{\mathrm{in}}})$
with exponents $x^{\tilde{m}_{\mathrm{out}}}$. We repeat the same argument to
the opposite side of $\mathfrak{j}$ with the opposite half $\gamma_{B}$ of the
small loop around $\mathfrak{j}$. Then, by the consistency of $\mathfrak{D}$,
we have

$\displaystyle\mathfrak{p}_{\gamma_{A},\mathfrak{D}_{\ell}}(c_{\mathrm{in}}x^{\tilde{m}_{\mathrm{in}}})\equiv\mathfrak{p}_{\gamma_{B},\mathfrak{D}_{\ell}}(c_{\mathrm{in}}x^{\tilde{m}_{\mathrm{in}}})\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.14)

Thus, passing through $Q_{0}$ does not change $\vartheta_{Q,\tilde{m}_{0}}$
modulo $x^{\tilde{m}_{0}}J^{\ell+1}$. Therefore, we obtain the equality
(0.7.13).

Next, we consider the remaining case $\partial\mathfrak{d}$. See Figure 17
(b). We may concentrate on a face $\mathfrak{c}$ of $\mathfrak{d}$ of
codimension 2 such that $\mathfrak{c}$ is not a joint. Then, $\mathfrak{c}$ is
contained in some walls of $\mathfrak{D}_{\ell}$ with the common normal vector
$n$. In this case, for each broken line for $Q$, there is a unique broken line
for $Q^{\prime}$ with the same attached incoming and outgoing monomial modulo
$x^{\tilde{m}_{0}}J^{\ell+1}$, thanks to the consistency of $\mathfrak{D}$.
Thus, passing through $Q_{0}$ does not change $\vartheta_{Q,\tilde{m}_{0}}$
modulo $x^{\tilde{m}_{0}}J^{\ell+1}$. ∎

######  Lemma 0.7.8 ([CPS, Lemma 4.8]).

Suppose that $Q$ and $Q^{\prime}$ are in different chambers of
$\mathrm{Supp}(\mathfrak{D}_{\ell})$. Then, for any admissible curve $\gamma$
from $Q$ to $Q^{\prime}$, we have

$\displaystyle\vartheta_{Q^{\prime},\tilde{m}_{0}}\equiv\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}(\vartheta_{Q,\tilde{m}_{0}})\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.15)

###### Proof.

We may assume that $Q$ and $Q^{\prime}$ are in adjacent chambers of
$\mathrm{Supp}(\mathfrak{D}_{\ell})$, since the general case follows from this
case by composing the automorphisms
$\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}$. Let $n\in N_{\mathrm{pr}}^{+}$ be
the normal vector of the walls separating $Q$ and $Q^{\prime}$. We decompose

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}=\sum_{\tilde{m}\in\tilde{M}^{\circ}}c_{\tilde{m}}x^{\tilde{m}}=\vartheta_{+}+\vartheta_{0}+\vartheta_{-},$
(0.7.16)
$\displaystyle\vartheta_{+}=\sum_{\scriptstyle\tilde{m}\in\tilde{M}^{\circ}\atop\scriptstyle\langle
n,m\rangle>0}c_{\tilde{m}}x^{\tilde{m}},\quad\vartheta_{0}=\sum_{\scriptstyle\tilde{m}\in\tilde{M}^{\circ}\atop\scriptstyle\langle
n,m\rangle=0}c_{\tilde{m}}x^{\tilde{m}},\quad\vartheta_{-}=\sum_{\scriptstyle\tilde{m}\in\tilde{M}^{\circ}\atop\scriptstyle\langle
n,m\rangle<0}c_{\tilde{m}}x^{\tilde{m}}.$ (0.7.17)

Also, we decompose
$\vartheta_{Q^{\prime},\tilde{m}_{0}}=\vartheta^{\prime}_{+}+\vartheta^{\prime}_{0}+\vartheta^{\prime}_{-}$
in the same way. Thanks to Lemma 0.7.7, we may take $Q$ and $Q^{\prime}$ to be
very close to the wall and also move parallel to the wall inside each chamber
as we like. Without loosing generality, one can assume that $\langle
n,Q\rangle>0$ and $\langle n,Q^{\prime}\rangle<0$.

$Q^{\prime}$$Q$$\textstyle{\bullet}$$\textstyle{\bullet}$$\textstyle{\bullet}$$\textstyle{\bullet}$
Figure 18: Broken lines in the proof of Lemma 0.7.8.

(i). $\vartheta_{+}$ and $\vartheta^{\prime}_{+}$. For any exponent
$\tilde{m}$ of $\vartheta^{\prime}_{+}$, $\langle n,m\rangle>0$ holds. Then,
any broken line for $\vartheta^{\prime}_{+}$ is obtained by extending a broken
line for $\vartheta_{+}$. See Figure 18. Thus, we have

$\displaystyle\vartheta_{+}^{\prime}\equiv\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}(\vartheta_{+})\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.18)

(ii). $\vartheta_{-}$ and $\vartheta^{\prime}_{-}$. The situation is opposite,
and we have

$\displaystyle\vartheta_{-}\equiv\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}^{-1}(\vartheta^{\prime}_{-})\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.19)

(iii). $\vartheta_{0}$ and $\vartheta^{\prime}_{0}$. The broken lines for both
do not cross the walls separating $Q$ and $Q^{\prime}$ in the final step.
Then, the proof of Lemma 0.7.7 applies to $\vartheta_{0}$ and
$\vartheta^{\prime}_{0}$, and we have

$\displaystyle\vartheta_{0}^{\prime}\equiv\vartheta_{0}\equiv\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}(\vartheta_{0})\mod
x^{\tilde{m}_{0}}J^{\ell+1}.$ (0.7.20)

We can unify the above three equalities into a single one

$\displaystyle\vartheta_{Q^{\prime},\tilde{m}_{0}}\equiv\mathfrak{p}_{\gamma,\mathfrak{D}_{\ell}}(\vartheta_{Q,\tilde{m}_{0}})\mod
x^{\tilde{m}_{0}}J^{\ell+1}$ (0.7.21)

as desired. ∎

Since $\ell$ is arbitrary in (0.7.15), we have the desired equality (0.7.12).
This completes the proof of Proposition 0.7.6.

###  0.7.3 Positivity and mutation invariance of theta functions

Now let us specialize to a CSD $\mathfrak{D}_{\mathfrak{s}}$ for a given seed
$\mathfrak{s}$.

We have the following fundamental result on the positivity of theta functions
for a CSD.

######  Theorem 0.7.9 (Positivity of theta functions [GHKK18, Theorem 1.13 &
Remark 3.2]).

For a CSD $\mathfrak{D}_{\mathfrak{s}}$, every theta function
$\vartheta_{Q,\tilde{m}_{0}}\in x^{\tilde{m}_{0}}\mathbbm{k}[[P_{1}]]$ has
only positive integer coefficients.

###### Proof.

We choose a positive realization of $\mathfrak{D}_{\mathfrak{s}}$ in Theorem
0.5.2, so that every wall function $f[tn]^{s}$ of
$\mathfrak{D}_{\mathfrak{s}}$ has only positive integer coefficients. Then,
the product in the right hand side of (0.7.5) has also only positive integer
coefficients. Thus, $c_{\gamma}$ in (0.7.1) is a positive integer. Therefore,
the theta function defined by (0.7.6) has only positive integer coefficients.
∎

We conclude the section by showing the mutation invariance of theta functions
of a CSD and its consequences. Let us define a piecewise-linear transformation
in parallel to $T_{k}$ in (0.6.25),

$\displaystyle\begin{matrix}\tilde{T}_{k}=\tilde{T}_{k,\mathfrak{s}}:&\tilde{M}^{\circ}&\rightarrow&\tilde{M}^{\circ}\\\
&\tilde{m}&\mapsto&\tilde{m}+[\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}]_{+}p_{1}^{*}(e_{k})=\begin{cases}\tilde{S}_{k}(\tilde{m})&m\in\mathcal{H}_{k}^{+},\\\
\tilde{m}&m\in\mathcal{H}_{k}^{-},\end{cases}\end{matrix}$ (0.7.22)

where $\tilde{S}_{k}(\tilde{m})$ is the one in (0.6.10).

######  Proposition 0.7.10 ([GHKK18, Proposition 3.6]).

For a given seed $\mathfrak{s}$, let
$\mathfrak{s}^{\prime}=\mu_{k}(\mathfrak{s})$ be the mutation of
$\mathfrak{s}$ in direction $k$.

(a). There is a one-to-one correspondence between the broken lines for
$\tilde{m}_{0}$ with endpoint $Q$ with respect to a CSD
$\mathfrak{D}_{\mathfrak{s}}$ and the broken lines for
$\tilde{T}_{k}(\tilde{m}_{0})$ with endpoint $T_{k}(Q)$ with respect to a CSD
$\mathfrak{D}_{\mathfrak{s}^{\prime}}=T_{k}(\mathfrak{D}_{\mathfrak{s}})$. The
locus of the line is transformed by the piecewise-linear map $T_{k}$ in
(0.6.25).

(b). (Mutation invariance of theta functions.) The theta function
$\vartheta^{\mathfrak{s}}_{Q,\tilde{m}_{0}}$ with respect to
$\mathfrak{D}_{\mathfrak{s}}$ and the one
$\vartheta^{\mathfrak{s}^{\prime}}_{Q,\tilde{m}_{0}}$ with respect to
$\mathfrak{D}_{\mathfrak{s}^{\prime}}$ are related by

$\displaystyle\vartheta^{\mathfrak{s}^{\prime}}_{T_{k}(Q),\tilde{T}_{k}(\tilde{m}_{0})}=\begin{cases}\tilde{S}_{k}(\vartheta_{Q,\tilde{m}_{0}}^{\mathfrak{s}})&Q\in\mathcal{H}_{k}^{+},\\\
\vartheta_{Q,\tilde{m}_{0}}^{\mathfrak{s}}&Q\in\mathcal{H}_{k}^{-}.\end{cases}$
(0.7.23)

where $\mathcal{H}_{k}^{\pm}$ is the one in (0.6.23), and $\tilde{S}_{k}$ acts
on the exponents of $x$ in $\vartheta_{Q,\tilde{m}_{0}}^{\mathfrak{s}}$.

###### Proof.

(a). In $\mathcal{H}_{k}^{-}$ the walls and broken lines are unchanged, while
in $\mathcal{H}_{k}^{+}$ the walls and broken lines are transformed by the
linear map $S_{k}$. Thus, the rule for broken lines are preserved in each
$\mathcal{H}_{k}^{\pm}$. Therefore, we only need to check the rule when a
broken line crosses the wall
$\mathbf{w}_{e_{k}}=[e_{k}^{\perp},1+x^{p_{1}^{*}(e_{k})}]_{e_{k}}$.

Case (1). Suppose that a broken line $\gamma$ with respect to
$\mathfrak{D}_{\mathfrak{s}}$ crosses $e_{k}^{\perp}$ from
$\mathcal{H}_{k}^{-}$ with attached monomial $cx^{\tilde{m}}$. Then, $\langle
e_{k},-\tilde{m}\rangle_{1}>0$. After crossing $e_{k}^{\perp}$, the new
attached monomial is a term in
$cx^{\tilde{m}}(1+x^{p_{1}^{*}(e_{k})})^{|\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}|}$.
Applying $\tilde{S}_{k}$, it transforms to a term in

$\displaystyle\begin{split}&\quad\
cx^{\tilde{m}+\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}p_{1}^{*}(e_{k})}(1+x^{p_{1}^{*}(e_{k})})^{|\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}|}\\\
&=cx^{\tilde{m}-|\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}|p_{1}^{*}(e_{k})}(1+x^{p_{1}^{*}(e_{k})})^{|\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}|}\\\
&=cx^{\tilde{m}}(1+x^{p_{1}^{*}(e^{\prime}_{k})})^{|\langle\delta_{k}e^{\prime}_{k},\tilde{m}\rangle_{1}|}.\end{split}$
(0.7.24)

This is the rule for a broken line with respect to
$T_{k}(\mathfrak{D}_{\mathfrak{s}})$.

Case (2). Suppose that a broken line $\gamma$ with respect to
$\mathfrak{D}_{\mathfrak{s}}$ crosses $e_{k}^{\perp}$ from
$\mathcal{H}_{k}^{+}$ with attached monomial $cx^{\tilde{m}}$. Then, $\langle
e_{k},-\tilde{m}\rangle_{1}<0$. After crossing $e_{k}^{\perp}$, the new
attached monomial is a term in
$cx^{\tilde{m}}(1+x^{p_{1}^{*}(e_{k})})^{|\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}|}$,
which is rewritten as

$\displaystyle
cx^{\tilde{m}+\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}p_{1}^{*}(e_{k})}(1+x^{p_{1}^{*}(e^{\prime}_{k})})^{|\langle\delta_{k}e^{\prime}_{k},\tilde{m}\rangle_{1}|}.$
(0.7.25)

Since the attached monomial for the incoming segment is transformed from
$cx^{\tilde{m}}$ to
$cx^{\tilde{m}+\langle\delta_{k}e_{k},\tilde{m}\rangle_{1}p_{1}^{*}(e_{k})}$,
this is the correct rule for a broken line with respect to
$T_{k}(\mathfrak{D}_{\mathfrak{s}})$.

(b). Let $\gamma$ be a broken line with respect to
$\mathfrak{D}_{\mathfrak{s}}$, and let $\gamma^{\prime}=T_{k}(\gamma)$ be the
corresponding broken line with respect to
$\mathfrak{D}_{\mathfrak{s}^{\prime}}$. Let $c_{\gamma}x^{\tilde{m}_{\gamma}}$
and $c^{\prime}_{\gamma}x^{\tilde{m}^{\prime}_{\gamma}}$ be the monomials
attached with $\gamma$ and $\gamma^{\prime}$, respectively. Then, by (a), we
have

$\displaystyle
c^{\prime}_{\gamma}x^{\tilde{m}^{\prime}_{\gamma}}=\begin{cases}c_{\gamma}x^{\tilde{S}_{k}(\tilde{m}_{\gamma})}&Q\in\mathcal{H}_{k}^{+},\\\
c_{\gamma}x^{\tilde{m}_{\gamma}}&Q\in\mathcal{H}_{k}^{-}.\end{cases}$ (0.7.26)

Thus, we have (0.7.23). ∎

We have the following corollary of Propositions 0.7.5 and 0.7.10.

######  Corollary 0.7.11 ([GHKK18, Corollary 3.9]).

Let $\mathfrak{D}_{\mathfrak{s}}$ be a CSD, and let
$\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}$ be the $G$-cone in (0.6.75).
For $\tilde{m}_{0}\in\tilde{M}^{\circ}$, suppose that
$m_{0}\in\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}\cap M^{\circ}$. Also,
suppose that
$Q\in\mathrm{Int}(\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}})$. Then, we
have

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}=x^{\tilde{m}_{0}}.$ (0.7.27)

###### Proof.

By (0.6.75),

$\displaystyle\mathcal{C}^{+}_{\mathfrak{s}_{i}}=(T_{k_{i-1},\mathfrak{s}_{i-1}}\circ\cdots\circ
T_{k_{0},\mathfrak{s}_{0}})(\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}).$
(0.7.28)

Let $T:=T_{k_{i-1},\mathfrak{s}_{i-1}}\circ\cdots\circ
T_{k_{0},\mathfrak{s}_{0}}$ and
$\tilde{T}:=\tilde{T}_{k_{i-1},\mathfrak{s}_{i-1}}\circ\cdots\circ\tilde{T}_{k_{0},\mathfrak{s}_{0}}$.
Let $S:M_{\mathbb{R}}\rightarrow M_{\mathbb{R}}$ be the linear map that
locally coincides with $T$ on the cone
$\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}$. Let $\tilde{T}$ and
$\tilde{S}$ be the corresponding maps for
$\tilde{M}^{\circ}\rightarrow\tilde{M}^{\circ}$. Then, applying (0.7.23)
repeatedly, we obtain

$\displaystyle\vartheta^{\mathfrak{s}_{i}}_{T(Q),\tilde{T}(\tilde{m}_{0})}=\tilde{S}(\vartheta_{Q,\tilde{m}_{0}}),$
(0.7.29)

where $\tilde{S}$ acts on the exponents of $x$. Meanwhile,
$\vartheta^{\mathfrak{s}_{i}}_{T(Q),\tilde{T}(\tilde{m}_{0})}=x^{\tilde{T}(\tilde{m}_{0})}$
by Proposition 0.7.5. Thus, we have
$\vartheta_{Q,\tilde{m}_{0}}=x^{\tilde{m}_{0}}$. ∎

Combining Proposition 0.7.6 and Corollary 0.7.11, we have another fundamental
theorem on theta functions for a CSD, which provides the identification of
theta functions with cluster monomials or cluster variables for the
corresponding cluster pattern.

######  Theorem 0.7.12 ([GHKK18, Theorem 4.9]).

Let $\mathfrak{D}_{\mathfrak{s}}$ be a CSD, and let
$\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}$ be the $G$-cone in (0.6.75).
For $\tilde{m}_{0}\in\tilde{M}^{\circ}$, suppose that
$m_{0}\in\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}}\cap M^{\circ}$. Also,
suppose that $Q\in\mathrm{Int}(\mathcal{C}_{\mathfrak{s}}^{+})$. Then, for any
admissible curve $\gamma$ from any point in
$\mathrm{Int}(\mathcal{C}^{\mathfrak{s}}_{\mathfrak{s}_{i}})$ to $Q$, we have

$\displaystyle\vartheta_{Q,\tilde{m}_{0}}=\mathfrak{p}_{\gamma,\mathfrak{D}_{\mathfrak{s}}}(x^{\tilde{m}_{0}}).$
(0.7.30)

Notes

All results are taken from [GHKK18, §3, §4] and [CPS, §4] adapted to the
present formulation of scattering diagrams.

##  .8 Source code of Ordering Algorithm

For the reader’s convenience we present a source code of Ordering Algorithm
(Algorithm 0.5.7) written for SageMath 9.4 111Sage Mathematics Software, The
Sage Development Team, https://www.sagemath.org.. It is an almost direct
translation of Algorithm 0.5.7, and no effort is made for the efficiency and
sophistication. (If any minor or major bug is found, please let us know.)

The following examples explain how to use it on Sage Notebook.

######  Example .8.1.

(a). Let us order the product $\Psi[e_{2}]^{5}\Psi[e_{1}]^{3}$ modulo
$G^{>5}$.

    In:
    order(5,[[0,1,5],[1,0,3]])
    Out:
    [[1, 0, 3],
     [3, 1, 5],
     [2, 1, 15],
     [3, 2, 125],
     [1, 1, 15],
     [2, 2, 60],
     [2, 3, 270],
     [1, 2, 30],
     [1, 3, 30],
     [1, 4, 15],
     [0, 1, 5]]

(b). One can order any (not necessarily anti-ordered) product with factors of
the form $\Psi[tn]^{s/t}$ $(n\in N_{\mathrm{pr}}^{+},\ s,t\in\mathbb{Z}_{>0})$
as specified in Algorithm 0.5.7.

    In:
    order(6,[[0,1,1],[1,2,1],[2,1,1],[2,2,1/2],[1,0,2]])
    Out:
    [[1, 0, 2],
     [4, 1, 1],
     [3, 1, 2],
     [2, 1, 2],
     [4, 2, 19],
     [3, 2, 16],
     [1, 1, 2],
     [2, 2, 13/2],
     [3, 3, 33],
     [2, 3, 10],
     [1, 2, 1],
     [2, 4, 9/2],
     [1, 3, 1],
     [0, 1, 1]]

Enjoy exploring the Badlands!

Source code of Ordering Algorithm for SageMath 9.4

⬇

1 # Ordering Algorithm for dilogarithm products

2 # written by Tomoki Nakanishi for SageMath 9.4

3 # ver.2021.10.30.b

4 # usage: order(L,C)

5

6def bf(a,b):

7 # bilinear form {a,b} such that {e_2,e_1}=1

8 # a, b: rank 2 integer vectors

9 res = a[1]*b[0]-a[0]*b[1]

10 return res

11

12def check_ordered(C):

13 # check dilog product C is ordered completely

14 # C has factors of [n]^{c}, c>0 any, len(C)>1

15 l = len(C) # l>1

16 for i in range(l-1):

17 if bf(C[i][:2],C[i+1][:2]) > 0: # if {n’,n}>0

18 return False # not ordered

19 elif bf(C[i][:2],C[i+1][:2]) == 0: # {n’,n}=0

20 if C[i][0]+C[i][1] > C[i+1][0]+C[i+1][1]: # deg(n’)>deg(n)

21 return False # not ordered

22 return True # ordered

23

24def decompose_initial(C):

25 # decompose dilog product C for the main routine (1)

26 # C has factors [tn]^{s/t}

27 res = []

28 l = len(C)

29 for i in range(l):

30 if C[i][2].is_integer() and C[i][2] > 1: # if the power is integer >1

31 for j in range(C[i][2]):

32 res = res+[[C[i][0],C[i][1],1]] # decompose to [tn]’s

33 else:

34 res = res+[C[i]]

35 return res

36

37def decompose(p,C):

38 # decompose dilog product C for subroutine (p), p>1

39 # p>0: integer

40 # C=[n’]^{c’}[n]^{c}) is an p-exchangeable pair

41 res = []

42 c0 = C[0][2] # multiple of 1/p

43 c1 = C[1][2] # multiple of 1/p

44 for i in range(p*c0):

45 res = res+[C[0][:2]+[1/p]] # decompose to [n]^{1/p}

46 for i in range(p*c1):

47 res = res+[C[1][:2]+[1/p]] # decompose to [n]^{1/p}

48 return res

49

50def join(C):

51 # join factors with common [n] in dilog product C

52 # C has factors of [n]^{c}, len(C)>1

53 res = C

54 l = len(res) # l>1

55 flag_joined = False

56 while not flag_joined:

57 l = len(res)

58 for i in range(l-1):

59 if res[i][:2] == res[i+1][:2]:

60 res = res[:i]+[res[i][:2]+[res[i][2]+res[i+1][2]]]+res[i+2:]

61 break

62 if i == l-2: # if this is the last pair, it is joined

63 flag_joined = True

64 return res

65

66def order_p_partial(L,p,C):

67 # order dilog product C up to G^{>L} by p-pentagon (i.e.,{n’,n}=p)

68 # L,p>0: integer

69 # C has factors of [n]^{1/p}, len(C)>1

70 res = C

71 flag_ordered_p = False

72 while not flag_ordered_p: # finish if p-ordered

73 l = len(res) # l>1

74 for i in range(l-1):

75 if bf(res[i][:2],res[i+1][:2]) < 0: # if ordered, do nothing

76 pass

77 elif bf(res[i][:2],res[i+1][:2]) == 0: # if parallel

78 if res[i][0]+res[i][1] > res[i+1][0]+res[i+1][1]:

79 # if deg(left)>deg(right), commute it

80 res = res[:i]+[res[i+1]]+[res[i]]+res[i+2:]

81 break

82 elif bf(res[i][:2],res[i+1][:2]) > 0: # if anti-ordered

83 if res[i][0]+res[i][1]+res[i+1][0]+res[i+1][1] > L:

84 # if deg(n+n’)>L, commute it

85 res = res[:i]+[res[i+1]]+[res[i]]+res[i+2:]

86 break

87 elif bf(res[i][:2],res[i+1][:2]) == p:

88 # else if {n’,n}=p, apply p-pentagon relation

89 if res[i][2] == 1/p and res[i+1][2]== 1/p:

90 res = res[:i]+[res[i+1]]+[[res[i][0]+res[i+1][0],\

91 res[i][1]+res[i+1][1],1/p],res[i]]+res[i+2:]

92 break

93 if i == l-2: # if we reach the last pair, it is p-ordered

94 flag_ordered_p = True

95 return res

96

97def order_p(L,p,C):

98 # subroutine (p); recursively defined

99 # L,p>0: integer

100 # for p=1, C has factors [tn]^{s/t} (n:primitive vector, s,t>0 integers)

101 # for p>1, C has factors of [n]^{1/p}

102 # len(C)>1

103 res = C

104 while not check_ordered(res): # do until completely ordered

105 if p == 1: # only in the main routine (1)

106 res=decompose_initial(res) # decompose to [n]

107 res = order_p_partial(L,p,res) # p-order the product

108 l = len(res) # l>1

109 for i in range(l-1):

110 if bf(res[i][:2],res[i+1][:2]) > 0: # find q-admissible pair

111 q = bf(res[i][:2],res[i+1][:2])

112 Cin = decompose(q,[res[i],res[i+1]]) # decompose to [n]^{1/q}

113 Cout = order_p(L,q,Cin) # go to subroutine (q) (recursively)

114 res = res[:i]+Cout+res[i+2:]

115 break

116 res = join(res)

117 return res

118

119def order(L,C):

120 """␣␣␣␣L:␣integer␣>␣0,␣C:␣list␣of␣[n1,n2,c]

121␣␣␣␣n1,n2:␣integers>=0,␣c:␣rational>0␣such␣that␣c*gcd(n1,n2)=integer>0␣"""

122 # main routine (1)

123 # L>0: integer

124 # C has factors [tn]^{s/t} (n:primitive positive vector, s,t>0 integers)

125 if not (L > 0 and L.is_integer()): # check L

126 return print("the␣first␣argument␣is␣illegal")

127 l = len(C) # l>=1

128 for i in range(l):

129 n1 = C[i][0]

130 n2 = C[i][1]

131 if not (n1.is_integer() and n2.is_integer()\

132 and n1 >= 0 and n2 >= 0 and n1+n2 > 0): # check (n1,n2)

133 return print("the␣second␣argument␣is␣illegal")

134 s = C[i][2]*math.gcd(n1,n2)

135 if not (s > 0 and s.is_integer()): # check c

136 return print("the␣second␣argument␣is␣illegal")

137 if l == 1: # for a single factor, just return

138 return C

139 res = order_p(L,1,C) # enter the main routine (1)

140 return res \end{lstlisting}

141

142

143%%%%%%%%%%%% part 3 %%%%%%%%%%%%%%

144%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

## References

  * [Bri17] T. Bridgeland, _Scattering diagrams, Hall algebras and stability conditions_ , Algebraic Geometry 4 (2017), 523–561; arXiv:1603.00416 [math.AG].
  * [CI12] G Cerulli Irelli, _Cluster algebras of type $A_{2}^{(1)}$_, Algebras and Representation Theory 15 (2012), 977–1021; arXiv:0904.2543 [math.RA].
  * [CPS] M. Carl, M. Pumperla, and B. Siebert, _A tropical view of Landau-Ginzburg models_ , preprint, 2010, available at https://www.math.uni-hamburg.de/home/siebert/preprints/LGtrop.pdf.
  * [DM21] B. Davison and T. Mandel, _Strong positivity for quantum theta bases of quantum cluster algebras_ , Invent. Math. (2021), published online; arXiv:1910.12915 [math.RT].
  * [FG09] V. V. Fock and A. B. Goncharov, _Cluster ensembles, quantization and the dilogarithm_ , Annales Sci. de l’École Norm. Sup. 42 (2009), 865–930; arXiv:math/0311245 [math.AG].
  * [FG16] V. Fock and A. Goncharov, _Cluster poisson varieties at infinity_ , Selecta Math. (N.S.) 22 (2016), 2569–2589; arXiv:1104.0407 [math.AG].
  * [FZ02] S. Fomin and A. Zelevinsky, _Cluster algebras I. Foundations_ , J. Amer. Math. Soc. 15 (2002), 497–529 (electronic); arXiv:math/0104151 [math.RT].
  * [FZ03a] , _Cluster algebras II. Finite type classification_ , Invent. Math. 154 (2003), 63–121; arXiv:math/0208229 [math.RA].
  * [FZ03b] , _Y-systems and generalized associahedra_ , Ann. of Math. 158 (2003), 977–1018; arXiv:hep–th/0111053.
  * [FZ07] , _Cluster algebras IV. Coefficients_ , Compositio Mathematica 143 (2007), 112–164; arXiv:math/0602259 [math.RT].
  * [GHK15] M. Gross, P. Hacking, and S. Keel, _Birational geometry of cluster algebras_ , Algebr. Geom. 2 (2015), 137–175; arXiv:math/1309.2573.
  * [GHKK18] M. Gross, P. Hacking, S. Keel, and M. Kontsevich, _Canonical bases for cluster algebras_ , J. Amer. Math. Soc. 31 (2018), 497–608; arXiv:1411.1394 [math.AG].
  * [GPS10] M. Gross, R. Pandharipande, and B. Siebert, _The tropical vertex_ , Duke Math. J. (2010), no. 153, 297 – 362; arXiv:0902.0779 [math.AG].
  * [Gro11] M. Gross, _Tropical geometry and mirror symmetry_ , CBMS Regional Conf. Ser. in Math., no. 114, Amer. Math. Soc., 2011.
  * [GS11] M. Gross and B. Siebert, _From affine geometry to complex geometry_ , Annals of Math. 174 (2011), 95–138; arXiv:math/0703822.
  * [IOTW] K. Igusa, K. Orr, G. Todorov, and J. Weyman, _Modulated semi-invariants_ , arXiv:1507.03051 [math.RT].
  * [IT21] K Igusa and G. Todorov, _Picutre groups and maximal green sequences_ , Electronic Research Archive 29 (2021), 3031–3068; arXiv:2007.14584 [math.RT].
  * [ITW] K. Igusa, G. Todorov, and J. Weyman, _Picture groups of finite type and cohomology in type $A_{n}$_, arXiv:1609.02636 [math.RT].
  * [Jac79] N. Jacobson, _Lie algebras_ , Dover Publications, New York, 1979.
  * [KS06] M. Kontsevich and Y. Soibelman, _Affine structures and non-Archimedean analytic spaces_ , Prog. Math. 244 (2006), 321–385; arXiv:math/0406564.
  * [KS08] , _Stability structures, motivic Donaldson-Thomas invariants and cluster transformations_ , 2008, arXiv:0811.2435 [math.AG].
  * [KS14] , _Wall-crossing strunctures in Donaldson-Thomas invariants, integrable systems and mirror symmetry_ , Homological mirror symmetry and tropical geometry, Lect. Notes Unione Ital., vol. 15, Springer, 2014, pp. 197–308; arXiv:1303.3253 [math.AG].
  * [Lew81] L. Lewin, _Polylogarithms and associated functions_ , North-Holland, Amsterdam, 1981.
  * [Mag54] W. Magnus, _On the exponential solution of differential equations for a linear operator_ , Commun. Pure Appl. Math. 7 (1954), 649–673.
  * [Mat21] K. Matsushita, _Consistency relations of rank 2 cluster scattering diagrams of affine type and the pentagon relation_ , 2021, arXiv:2112.04743 [math.QA].
  * [Mul16] G. Muller, _The existence of a maximal green sequence is not invariant under quiver mutation_ , Electron. J. Combin. 23 (2016), #P2.47, 23 pages; arXiv:1503.04675.
  * [Nak21] T. Nakanishi, _Synchronicity phenomenon in cluster patterns_ , J. London Math. Soc. 103 (2021), 1120–1152; arXiv:1906.12036 [math.RA].
  * [NC13] A. Nájera Chávez, _On the c-vectors and g-vectors of the Markov cluster algebra_ , Séminaire Lontharingien de Combinatoire 69 (2013), B29d;arXiv:1112.2357.
  * [Rea] N. Reading, _Dominance phenomena: mutation, scattering, and cluster algebras_ , arXiv:1802.10107 [math.CO].
  * [Rea14] , _Universal geometric cluster algebras_ , Mathematische Zeitschrift 277 (2014), 499–547; arXiv:1209.3987.
  * [Rea20a] , _A combinatorial approach to scattering diagrams_ , Algebraic Combinatorics 3 (2020), 603–636; arXiv:1806.05094 [math.CO].
  * [Rea20b] , _Scattering fans_ , Int. Math. Res. Notices 23 (2020), 9640–9673; arXiv:1712.06968.
  * [Rei09] M. Reineke, _Poisson automorphisms and quiver moduli_ , J. Inst. Math. Jussieu 9 (2009), 653–667, arXiv:0804.3214 [math.RT].

## Index

  * admissible curve Definition 0.2.5
  * admissible region §0.6.5
  * $B$-equivalent Definition 0.1.22
  * Badlands §0.3.5, 5th item
  * Baker-Campbell-Hausdorff (BCH) formula §0.1.2
  * broken line Definition 0.7.1
  * chamber item (a).
    * $G$- Definition 0.6.23
  * cluster element Definition 0.3.13
  * cluster scattering diagram (CSD) Definition 0.3.11
  * cone Definition 0.2.1
    * convex rational polyhedral Definition 0.2.1
    * $G$- Definition 0.6.23
    * positivie/negative §0.2.4
    * strongly convex Definition 0.2.1
  * decomposition of $G$
    * at $z$ §0.2.3
    * by $e_{i}$ §0.3.4
    * by $n$ §0.3.1
  * degree function §0.1.1
  * dilogarithm, see Euler dilogarithm
  * dilogarithm element Definition 0.1.10
  * Euler dilogarithm §0.1.4
  * fan
    * $G$- §0.6.6
  * finiteness condition Definition 0.2.3
  * fixed data Definition 0.1.1
    * nondegenerate Definition 0.1.1
  * $G$-chamber Definition 0.6.23
  * $G$-cone Definition 0.6.23
  * $G$-fan §0.6.6
  * general point (in $M_{\mathbb{R}}$) Definition 0.2.8
  * general position Definition 0.7.2
  * intersection sign Definition 0.2.6
  * joint Definition 0.5.13
    * parallel Definition 0.5.13
    * perpendicular Definition 0.5.13
  * matrix
    * exchange §0.1.1
    * skew-symmetrizable §0.1.1
  * mutation
    * of a cluster scattering diagram Definition 0.6.7
    * seed (for a fixed data) Definition 0.6.1
  * mutation invariance
    * of a cluster scattering diagram Theorem 0.6.9
    * of theta functions Proposition 0.7.10
  * normal vector (of a wall) Definition 0.2.2
  * normalization factor Definition 0.1.18
  * normalized automorphism Definition 0.4.9
  * normalized form §0.1.5
  * ordered/anti-ordered §0.3.5, §0.5.2
  * Ordering Algorithm Algorithm 0.5.7
    * source code 8.§
  * Ordering Lemma Proposition 0.5.4
  * $p$-exchangeable pair Definition 0.5.5
  * parallel subgroup Definition 0.1.5
  * path-ordered product Definition 0.2.6
  * pentagon relation Proposition 0.1.14
  * positive element §0.1.1
  * positive realization §0.5.1
  * positivity (of theta functions) Theorem 0.7.9
  * primitive §0.1.1
  * principal extension
    * of a fixed data Definition 0.4.1
  * principal $x$-representation §0.4.2
  * reduction (of a scattering diagram) Definition 0.2.3
  * rescaling (of a fixed data) Definition 0.1.16
  * scattering diagram Definition 0.2.3
    * cluster Definition 0.3.11
    * consistent Definition 0.2.13
    * equivalent Definition 0.2.7
    * extended Definition 0.6.14
    * trivial Definition 0.2.12
    * with minimal support Definition 0.2.23
  * seed
    * for a fixed data Definition 0.1.2
  * singular locus (of a scattering diagram) Definition 0.2.4
  * structure group §0.1.2
    * extended §0.6.2
  * support
    * of a scattering diagram Definition 0.2.4
    * of a wall Definition 0.2.2
  * theta function Definition 0.7.3
  * total wall element Definition 0.2.10
  * wall Definition 0.2.2
    * attached to a joint Construction 0.5.14
    * incoming Definition 0.3.9
    * outgoing Definition 0.3.9
    * unreachable Example 0.6.25
  * wall element Definition 0.2.2
  * wall function Definition 0.4.10
  * $x$-representation §0.1.3
    * principal §0.4.2
  * $x$-variable (cluster variable) §0.1.1
  * $y$-representation §0.1.3
  * $y$-variable (coefficient) §0.1.1
  * Zassenhaus formula §0.1.2