# Path Planning in a dynamic environment using Spherical Particle Swarm
Optimization

Mohssen E. Elshaar1,∗, Mohammed R. Elbalshy1,∗, A. Hussien 2,† and Mohammed
Abido2,†,‡
1Department of Aerospace Engineering, 2Department of Electrical Engineering
†Interdisciplinary Research Center for Sustainable Energy Systems (IRC-SES)
‡SDAIA-KFUPM Joint Research Center for Artificial Intelligence
King Fahd University of Petroleum and Minerals (KFUPM), Dharan, Saudi Arabia

###### Abstract

Efficiently planning an Unmanned Aerial Vehicle (UAV) path is crucial,
especially in dynamic settings where potential threats are prevalent. A
Dynamic Path Planner (DPP) for UAV using the Spherical Vector-based Particle
Swarm Optimisation (SPSO) technique is proposed in this study. The UAV is
supposed to go from a starting point to an end point through an optimal path
according to some flight criteria. Path length, Safety, Attitude and Path
Smoothness are all taken into account upon deciding how an optimal path should
be. The path is constructed as a set of way-points that stands as re-planning
checkpoints. At each path way-point, threats are allowed some constrained
random motion, where their exact positions are updated and fed to the SPSO-
solver. Four test scenarios are carried out using real digital elevation
models. Each test gives different priorities to path length and safety, in
order to show how well the SPSO-DPP is cabable of generating a safe yet
efficient path segments. Finally, a comparison is made to reveal the
persistent overall superior performance of the SPSO, in a dynamic environment,
over both the Particle Swarm Optimisation (PSO) and the Genetic Algorithm
(GA). The methods are compared directly, by averaging costs over multiple
runs, and by considering different challenging levels of obstacle motion. SPSO
outperformed both PSO and GA, showcasing cost reductions ranging from 330% to
675% compared to both algorithms.

###### Index Terms:

UAV, Path Planning, Dynamic Obstacles, Meta-heuristic Optimization.

††∗Equal Contribution.
<EMAIL_ADDRESS>
##  I Introduction

Unmanned aerial vehicles (UAVs), commonly referred to as drones, are
progressively gaining prevalence and proving their utility across diverse
domains and industries, ranging from agriculture and surveillance to disaster
response and aerial photography. One important factor that is crucial to the
efficient use of UAVs during deployment is path planning which involves the
systematic determination of the optimal route a drone or a UAV should follow
to achieve its mission objectives while adhering to various constraints and
minimizing risks [1, 2]. Key considerations in UAV path planning include
obstacle avoidance, energy efficiency, real-time decision-making, and
adaptation to changing conditions. Researchers and engineers have developed a
plethora of approaches to tackle UAV path planning challenges, from
traditional methods like $A^{*}$ search and Dijkstra’s algorithm to advanced
techniques like Rapidly-exploring Random Trees (RRT) and Genetic Algorithms
(GA).

The criteria for optimal path planning extend to diverse objectives depending
on the application, involving the maximization of detection probability,
flight time minimization or Pareto solutions achievement for the navigation of
multi-objective problems [3, 4, 5]. Furthermore, the intended trajectory must
adhere to relevant safety and feasibility limits within the operational
environment of the UAV. Safety considerations include ensuring that the path
can safely steer the UAV past environmental risks such as barriers. The path
must be feasible in terms of altitude, flight time, fuel consumption,
ascending angle and turning rate. As a result, path planning aimed at
increasing safety by permitting collision-free and practicable UAV motion
remains a constantly difficult task. Various path planning strategies have
been proposed including approaches such as graph search, cell decomposition,
potential field, and nature-inspired algorithms [6, 7, 8, 9, 10, 11, 12, 13].

Lately, nature-inspired approaches are increasingly favored in UAV path
planning due to their efficiency in accommodating dynamic UAV constraints and
seeking global optima in intricate scenarios. Various algorithms, drawing
inspiration from nature, are tailored for UAV path planning, including cuckoo
search [14], genetic algorithm (GA) [15, 16], differential evolution (DE)
[17], artificial bee colony (ABC) [18], ant colony optimization (ACO) [19],
and particle swarm optimization (PSO) [1, 15, 19, 20]. PSO, in particular,
finds widespread use and has seen several variants developed. PSO is a
population-based algorithm inspired by bird flocking and fish schooling
behavior that embodies swarm intelligence through cognitive and social
coherence [21]. These qualities allow each particle to navigate the solution
space based on individual and collective experiences, which differs from
standard evolutionary operators such as mutation and crossover. As a result,
PSO efficiently identifies global solutions with stable convergence,
outperforming alternative nature-inspired algorithms [22]. PSO variants
proposed in [23, 24, 25, 26, 27] share a population-based structure but differ
in search space representation and particle-encoded solutions.

Real-world problems are often dynamic, involving changing global and local
optimal solutions that vary over time. Dynamic Path Planners were designed
based on improved Simulated Annealing (SA) [28], Neural Networks (NN) [29],
Reinforcement Learning [30] and Fuzzy Logic [31]. Some studies argue that
Evolutionary Algorithms (EAs) are well-equipped to handle these challenges,
evidenced by their success in dynamic problems [32, 33, 34]. From these
algorithms, PSO has been employed in a variety of path planning issues in
dynamical contexts, including [35, 36, 37]. A new algorithm, known as
Spherical Vector-based Particle Swarm Optimization (SPSO), is introduced in
[38] to tackle the problem of path planning for UAVs in complex and static
threat-filled environments. It begins by creating a cost function that turns
the path planning problem into an optimization task, considering UAV
operational requirements and safety constraints. SPSO is then applied to
efficiently search the UAV’s configuration space to find the optimal path.
This is done by linking particle positions to critical flight parameters,
including speed, turning angle, and climb/dive angle. SPSO was tested
exclusively in a static environment, and its performance in a dynamic
environment remains unexplored. In essence, this work focuses on:

  1. 1.

Introducing Spherical Particle Swarm Optimization Dynamic Path Planner (SPSO-
DPP).

  2. 2.

Investigating the efficiency of SPSO-solver within a dynamic setting
characterized by the presence of randomly moving obstacles.

  3. 3.

Comparing the SPSO-DPP performance with PSO and GA performances in a dynamic
environment.

##  II Methodology

###  II-A Problem formulation

An optimisation problem is used to formulate the UAV path planning problem.
The constraints are obtained from the UAV flying constraints, and the path
objective function is defined as the path length [38].

####  II-A1 Path Cost

For many aerial picture and mapping applications, the primary optimum
criterion is to minimise the path length. The set of $n$ way-points that the
UAV must travel through represents the flight path. Each co-ordinate
$(x_{i},y_{i},z_{i})$ describes a way-point $\vec{P}_{i}$ in detail. The
separation between every two way-points is:

$d_{i}=\|\vec{P}_{i+1}-\vec{P}_{i}\|$

The cost $F_{1}$ that needs to be minimized is the summation of all distances
$d_{i}$, where $i\in\\{1,2,3,..,n-1\\}$.

$F_{1}=\sum_{i=1}^{n-1}d_{i}$ (1)

####  II-A2 Safety Constrain

In order to protect the UAV from any threats resulting from impediments in its
operating region, the planned route needs to be carefully developed. This
means that a safety constraint must be included in the planning phase. A set
that represents all possible threats, $K$, is defined. It is assumed that
every danger is confined within a cylindrical container, and that each
threat’s projection has a radius of $R_{k}$ and a central coordinate $C_{k}$.
The cost associated with potential threats for a given path segment
$\Delta\vec{P}_{i}=\vec{P}_{i+1}-\vec{P}_{i}$ is proportional to the distance,
$d_{k}$, between the segment and the central coordinate $C_{k}$ of each
danger. The threat cost $F_{2}$ is calculated along the way-points $P_{i}$
with regard to the obstacle set $K$ in the following method to account for the
physical dimensions of the UAV, which are defined by its diameter $D$ and the
critical safety threshold $S$ that defines the danger zone for probable
collisions:

$F_{2}=\sum_{i=1}^{n-1}\sum_{k=1}^{K}T_{k}\Delta\vec{P}_{i},$ (2)

where

$T_{k}=\begin{cases}0&\text{if }d_{k}>S+D+R_{k}\\\ (S+D+R_{k})-d_{k}&\text{if
}D+R_{k}<d_{k}\leq S+D+R_{k}\\\ \infty&\text{if}d_{k}\leq D+R_{k}.\end{cases}$
(3)

####  II-A3 Altitude Constrain

The flying altitude is frequently restricted during operation between the
minimum and maximum heights, the two designated extrema. For instance, the
camera must capture visual data at a particular resolution and field of view
for surveying and search applications, which limits the flying altitude. Let
$h_{\text{min}}$ and $h_{\text{max}}$ stand for the lowest and highest
heights, respectively. For each way-point $P_{i}$, the altitude cost is
calculated as:

$H_{i}=\begin{cases}\left|h_{i}-\frac{h_{\text{max}}+h_{\text{min}}}{2}\right|,&\text{if
}h_{\text{min}}\leq h_{i}\leq h_{\text{max}}\\\
\infty,&\text{otherwise}\end{cases}$ (4)

In this representation, the flying height above the ground is shown by
$h_{i}$. It is made so that out-of-range values are penalised and the average
height is maintained by $H_{i}$. The altitude cost is obtained by adding
together all of the way points’ $H_{i}$:

$F_{3}=\sum_{i=1}^{n}H_{i}$ (5)

where $n$ is the total number of way-points.

####  II-A4 Smoothness Constrain

The smoothness cost function assesses the rates of turning and climb, both of
which are vital for creating workable flight trajectories. The turning angle,
represented as $\phi_{i}$, signifies the angle between two successive segments
of the path. This angle is determined by projecting the vectors originating
from point $P^{\prime}_{i}$ to $P^{\prime}_{i+1}$ and from point
$P^{\prime}_{i+1}$ to $P^{\prime}_{i+2}$ onto the horizontal plane denoted as
$O_{xy}$.

$\phi_{ij}=\arctan\left(\frac{\|\overrightarrow{P^{\prime}_{i}P^{\prime}_{i+1}}\times\overrightarrow{P^{\prime}_{i+1}P^{\prime}_{i+2}}\|}{\overrightarrow{P^{\prime}_{i}P^{\prime}_{i+1}}\cdot\overrightarrow{P^{\prime}_{i+1}P^{\prime}_{i+2}}}\right)$
(6)

The angle between the path section from $P_{i}$ to $P_{i}+1$ and its
projection into the horizontal plane - that is, the line from $P^{\prime}_{i}$
to $P^{\prime}_{i}+1$ — is known as the angle of climb, indicated by the
symbol $psi_{i}$. The following formula can be used to compute this angle:

$\psi_{ij}=\arctan\left(\frac{z_{i+1}-z_{i}}{\|\overrightarrow{P^{\prime}_{i}P^{\prime}_{i+1}}\|}\right)$
(7)

The calculation of the smoothness cost is as follows:

$F_{4}=a_{1}\sum_{i=1}^{n-2}\phi_{i}+a_{2}\sum_{i=1}^{n-1}\left(\psi_{i}-\psi_{i-1}\right)$
(8)

Here, $a_{1}$ and $a_{2}$ represent the penalty coefficients associated with
the turning and climbing angles, respectively.

####  II-A5 Accumulative Cost Function

Taking into account considerations related to optimizing the path, ensuring
safety, and meeting altitude and smoothness constraints, we can formulate the
overall cost function as follows:

$F(X_{j})=\sum_{k=1}^{4}b_{k}F_{k}(X_{j})$ (9)

Here, each $b_{k}$ serves as a weight coefficient, and $F_{1}$ through $F_{4}$
represent the costs associated with different aspects: path length (1), safety
cost (2), flight altitude (5), and smoothness (8). The decision variable,
denoted as $X_{j}$, encompasses a list of $n$ way-points from the starting
point to the end point, $P_{ij}=(x_{ij},y_{ij},z_{ij})$, satisfying the
condition $P_{ij}\in O$, where $O$ signifies the operational space for the
UAV. With these definitions, the cost function $F$ becomes fully defined, and
the optimization problem (9) is aimed to be solved for the optimum $X^{\ast}$
in which the $F(X^{\ast})$ is minimum.

In the SPSO approach, each flight path is translated into a collection of
vectors, each specifying how the UAV moves from one waypoint to another. These
vectors are described in a spherical coordinate system, consisting of three
key elements: the magnitude $\rho\in[0,\text{path length}]$, the elevation
angle $\psi\in[-\pi/2,\pi/2]$, and the azimuth angle $\phi\in[-\pi,\pi]$. A
flight path indexed as $j$ and comprising of $N$ nodes can be represented as a
hyper-spherical vector with $3N$ dimensions:

$\displaystyle\Omega_{j}$
$\displaystyle=(\rho_{1j},\theta_{1j},\phi_{1j},r_{2j},\theta_{2j},\phi_{2j},$
(10) $\displaystyle\dots,r_{Nj},\theta_{Nj},\phi_{Nj});\text{where }N=n-2$

When we characterize a particle’s position as $\Omega_{i}$, the particle’s
velocity is defined by an incremental vector:

$\displaystyle\Delta\Omega_{j}$ $\displaystyle=(\Delta
r_{1j},\Delta\theta_{1j},\Delta\phi_{1j},\Delta
r_{2j},\Delta\theta_{2j},\Delta\phi_{2j},$ (11) $\displaystyle\dots,\Delta
r_{Nj},\Delta\theta_{Nj},\Delta\phi_{Nj})$

####  II-A6 Spherical vector-based PSO algorithm

The purpose of utilising spherical vectors in SPSO is to improve navigation
safety by establishing a relationship between the speed, turning angle, and
ascending angle of the UAV and the magnitude, elevation, and azimuth
components of these vectors. As a result, the possibility of discovering
excellent solutions is increased since the particles in SPSO investigate
solutions inside the configuration space as opposed to the Cartesian space.
Furthermore, the search space can be significantly reduced by explicitly
incorporating limitations related to turning and ascending angles by adjusting
the elevation and azimuth angles of the spherical vector. Fixing the magnitude
can help further restrict the search space while increasing the search
capacity in some circumstances, such as when the UAV stays at a fixed speed.
The Pseudo-code for Static Optimal Path Planning using SPSO is in [38].

###  II-B Environment Construction

The evaluation scenarios are constructed using actual digital elevation model
(DEM) maps obtained from LiDAR sensors. Regions of Christmas Island in
Australia, characterized by diverse terrain structures, are chosen [38].
Subsequently, these areas are enhanced to create different bench-marking
scenarios, illustrated in this section. Within these scenarios, threats are
depicted by green cylinders when seen in an isometric view, and as concentric
circles when seen in a top view (see figure LABEL:fig:EnvViews).

###  II-C Dynamic obstacles

In addressing the dynamic path planning problem for the Unmanned Aerial
Vehicle (UAV), way-points play a pivotal role as crucial checkpoints,
discretizing the UAV’s path through an environment teeming with randomly
moving obstacles. In this section, the path planning challenge is dynamically
tackled using an iterative solver operating online. This solver consistently
updates the UAV’s path in response to the dynamic nature of the environment,
characterized by obstacles exhibiting random motion, necessitating real-time
adjustments for effective UAV navigation. The threat dynamics are simulated by
random positional shifts, with each obstacle restricted to move in a circular
area with a specified radius, and centered at its current position. All
threats are assumed not to instantly jump on the UAV, making sure that the UAV
does not fall inside its dead-zone. No threat is allowed to overlap with
another threat as well.

The difficulty of the environment could be adjusted through changing both the
threat size and the radii of the motion restricted areas for each threat.

###  II-D SPSO-DPP

The SPSO dynamic path planner utilizes the SPSO-solver introduced by [38]. At
each way-point, the solver reevaluates the environment, incorporating changes
in obstacle positions. Notably, the SPSO-solver is employed at each way-point
to ensure the path’s coherence, optimizing it as the UAV progresses from one
way-point to the next. This adaptive process unfolds dynamically, ensuring the
UAV successfully navigates evolving scenarios, as visually depicted in figures
3, 4, 5 and 6.

Figure 2: SPSO-DPP Flow Chart

###  II-E Path representation

It is crucial to note that the UAV path is represented by a curve that is
split into three sections: 1) a dashed-line section that resembles the
temporary planned path from the next way-point to the destination; 2) a solid-
line portion that represents the finalized path from the UAV’s current way-
point to the next way-point; and 3) a dashed-dotted-line portion that
represents the path taken from the start point to the current point. The dash-
dotted-line portion not only serves as a representation of the finalized path
in space but also signifies the path that the UAV has previously traversed and
is no longer subject to updates. Similarly, the solid-line portion denotes the
current path taken by the UAV in the ongoing step or iteration and is not
revisited for re-planning. In contrast, the dashed-line portion outlines the
prospective path for all upcoming way-points, representing the segment that
necessitates optimization and re-planning. This distinction in path portions
ensures a focused approach, directing computational efforts towards refining
the UAV’s future path while preserving the integrity of its past and present
routes. Within the overall path from the start to the destination, 10 way-
points are positioned. Way-points 1 to 8 act as checkpoints for dynamic
updates, where the solver re-calibrates the UAV’s path. The last 2 way-points
represent a final path segment, requiring no further re-planning as it
constitutes a unique solution.

##  III Results

Figure 3: Base Case: Planned path at each way-point using the same weights as
in [38] Figure 4: Safer Path Case: Planned path at each way-point with
increasing the value of ${b_{2}}$ to 100 for ensuring a safer path for the UAV
Figure 5: Shorter Path Case: Planned path at each way-point with decreasing
${b_{2}}$ to 1 for prioritizing the distance cost to achieve more dominance in
the overall cost function Figure 6: Safer Than Ever Case: Planned path at each
way-point with increasing the value of ${b_{2}}$ to be 500 prioritizing safety
over distance minimization for ensuring a higher level of safety for the UAV

###  III-A Case Scenarios

Four case scenarios are presented to test the effect of cost weights,
specifically the safety and distance weights.

  * •

Case 1 - ”Base Case”: The weights are taken to be similar to [38]. The
parameters are set as follows: ${b_{1}}$ = 1, ${b_{2}}$ = 5, ${b_{3}}$ = 10,
and ${b_{4}}$ = 1.

  * •

Case 2 - ”Safer Path Case”: A modification is made specifically to ${b_{2}}$,
increasing its value to 100. This adjustment is implemented as a means of
ensuring a safer path for the UAV by elevating the threat cost coefficient.
The impact of this change is evident in the planned paths showcased in Figures
3 and 4, illustrating the discernible effect of altering the weights or
importance assigned to the distance and threat objectives within the problem
formulation.

  * •

Case 3 - ”Shorter Path Case”: ${b_{2}}$ is decreased to 1; strategically
prioritizing the distance cost to achieve more dominance in the overall cost
function.

  * •

Case 4 - ”Safer Than Ever Case”: ${b_{2}}$ value is increased to be 500
prioritizing safety over distance minimization. As a consequence, the
algorithm naturally generates longer paths, reflecting a deliberate trade-off.
This adjustment aims to ensure a higher level of safety for the UAV, even at
the expense of producing relatively more extended and resource-intensive
routes.

Each case is represented as 8 sub-figures resembling the past, present and
future path portion at each way-point.

It is obvious from the comparison of the case 1 and case 2 that planned path
tend have safer segments when giving the threat cost a higher weight. This can
be explicitly seen in Figures 3-a and 4-a.

Despite of the closely position obstacles in case 3, the planned path have
segments that go between these obstacles without much consideration of the
high threat the UAV may face. Figure 5-c is a resemblance of a relatively low
threat cost effect with a high distance cost. In Figure 5-d, e, and f, another
intriguing observation is that, at certain way-points, the solver opted for
nearly non-progressive segments of very small lengths instead of choosing a
longer route. This in a way shows the strength of the dynamic planner, where
it is not just capable of determining the best route according to the
momentary obstacles positions, but also, in future, may be further developed
to account for the obstacles direction and speed of motion.

Aggressively increasing the threat cost resulted in a relatively longer
planned route as seen in Figure 6. The solver did not go for a chance margin
of collision, not even a small one. Undoubtedly, these observations strongly
indicate a significant trade-off between safety and the distance of the path.

###  III-B Performance of SPSO-DPP against GA & PSO

The first three case scenarios in the previous section we re-run using SPSO-
DPP in comparison to GA and PSO. Two tests were conducted: one where the
random movement of the threats were held constant across all methods and
scenarios, and another where the the randomness of the threats movements
varied for each run. In the first test, the performance is presented as the
direct comparative costs of each method, while in the second test, it is
reported as the averaged costs from many runs for each scenario. Finally, an
additional test is conducted to evaluate methods average performance across
three different settings of threat dynamics.

The cost summery of all cases are tabulated and shown in Tables I, II and
Figure 7. The highlighted columns indicate that SPSO provide better results
for the cost compared to PSO and GA. Cost 1 signifies the expense from the
first way-point to the end of the path. Similarly, Cost 2 denotes the cost
from the second way-point, where the path is recalculated, to the end of the
path. This pattern continues for the subsequent costs, up to Cost 8, each
reflecting the recalculated path from the corresponding way-point to the
destination.

####  III-B1 First Test

SPSO demonstrated superior performance compared to both PSO and GA. In
comparison to PSO, SPSO exhibited cost reductions of up to 330%, 340%, and
632% in cases 1, 2, and 3, respectively. When compared to GA, SPSO
demonstrated cost reductions of up to 390%, 380%, and 675% in cases 1, 2, and
3, respectively. Furthermore, it is noted that the application of the SPSO-
solver leads to a consistently non-increasing trend in path cost as the UAV
progresses through the way-points. This suggests that the SPSO-solver
consistently identifies an optimal way-point for the current movement. In
contrast, the GA exhibited fluctuations at various instances (e.g., Case 1 -
costs 4 and 5, Case 2 - costs 1 and 2, and Case 3 - costs 3 and 4), indicating
that the selected way-point by the algorithm may not necessarily be optimal.

####  III-B2 Second Test

Another evaluation metric for the SPSO-DPP algorithm is to get averaged
performance statistics for the algorithm compared to the state of arts GA and
PSO by running it 10 times with random motion of the obstacles in each run and
calculate the average cost over the total number of runs. Figure 7 shows the
averaged costs of each method across all scenarios and runs. Even though all
methods find a solution path with an almost non-increasing cost, SPSO
outperforms GA and PSO significantly in a dynamic environment on average.

| Scenario 1 | Scenario 2 | Scenario 3
---|---|---|---
| SPSO | PSO | GA | SPSO | PSO | GA | SPSO | PSO | GA
Cost 1 | 1058 | 3269 | 4424 | 946 | 4173 | 3205 | 944 | 4016 | 2837
Cost 2 | 1013 | 2866 | 3127 | 783 | 2307 | 3504 | 637 | 3493 | 3612
Cost 3 | 790 | 2626 | 2650 | 708 | 2440 | 3405 | 477 | 3247 | 2231
Cost 4 | 655 | 2512 | 1861 | 753 | 2223 | 2704 | 368 | 2697 | 2574
Cost 5 | 492 | 1695 | 2418 | 455 | 1589 | 1793 | 237 | 1673 | 1839
Cost 6 | 410 | 1273 | 1478 | 362 | 1221 | 1135 | 226 | 1375 | 1117
Cost 7 | 346 | 1038 | 933 | 281 | 1029 | 800 | 226 | 1098 | 1038
Cost 8 | 184 | 789 | 819 | 271 | 843 | 787 | 226 | 937 | 908

TABLE I: SPSO vs PSO vs GA costs in Fixed Obstacles Motion Figure 7: SPSO vs
PSO vs GA Averaged Costs across all Scenarios

####  III-B3 Third Test

Using the base case Scenario (Case 1), SPSO-DPP, PSO and GA were tested across
various environment difficulties, simulated by varying the obstacles allowed
motion radii. The radii are specifically set to 50, 100, and 200 in this test.
Table II shows the average costs of 10 runs for each radius value. Again,
SPSO-DPP showed obvious domination over the other methods.

| R = 50 | R = 100 | R = 200
---|---|---|---
| SPSO | PSO | GA | SPSO | PSO | GA | SPSO | PSO | GA
Cost 1 | 1032.0 | 4226.8 | 2678.4 | 967.1 | 3460.0 | 4266.3 | 952.0 | 3490.4 | 2796.2
Cost 2 | 849.2 | 3764.4 | 3390.1 | 784.4 | 2780.0 | 3573.5 | 752.5 | 3070.0 | 4227.8
Cost 3 | 676.6 | 2803.6 | 2515.0 | 665.3 | 2432.5 | 2468.6 | 640.6 | 2120.2 | 2946.8
Cost 4 | 562.4 | 2041.1 | 2234.9 | 567.6 | 2516.3 | 2336.3 | 612.5 | 1790.5 | 1814.4
Cost 5 | 507.1 | 1687.2 | 2007.5 | 491.9 | 1892.4 | 1695.6 | 652.9 | 1482.3 | 1758.3
Cost 6 | 383.9 | 1232.9 | 1142.3 | 411.3 | 1225.2 | 1263.6 | 555.2 | 1382.0 | 1918.2
Cost 7 | 343.0 | 1123.7 | 1097.7 | 348.9 | 1108.7 | 892.6 | 375.7 | 973.9 | 1484.6
Cost 8 | 246.4 | 781.3 | 889.7 | 360.0 | 889.8 | 860.6 | 308.2 | 883.8 | 803.3

TABLE II: SPSO vs PSO vs GA Averaged costs Across Varying Threat allowed
Motion Radii

##  IV Limitations

The presented work depends on the instantaneous relative threats positions,
and does not take the velocity estimation of the dynamic threats into
consideration. Thus, the SPSO-DPP does not account for the threats velocity
vectors and their expected movement directions. Such estimation would enable
further investigation of how the performance of SPSO-DPP is influenced by
variations in environmental complexity.

##  V Conclusion

Finally, the paper introduces the SPSO-DPP. The approach uses an online solver
to guide the UAV across areas with randomly moving obstacles by using way-
points as checkpoints. At each way-point, the SPSO-solver is used to
adaptively optimize the path, provide coherence, and make adjustments in real-
time for effective UAV navigation. Four scenarios have been presented, each of
which modifies the parameters to reduce distance or prioritize safety. The
results demonstrate that SPSO consistently produces superior cost values than
both GA and PSO, demonstrating the algorithm’s efficacy in generating optimal
paths. The effect of prioritizing distance reduction above hazard avoidance is
best shown in Case 3, where the decision leads to a more direct path that
circumvents obstacle constraints. This realization emphasizes the necessity of
carefully weighing objectives against constraints in situations including
dynamic UAV path planning.

In various conducted tests, the SPSO-DPP algorithm consistently performs
better in path cost optimization than GA and PSO. Whether the approaches are
directly compared, the costs are averaged over several runs, or various levels
of obstacle motion are taken into account, SPSO-DPP performs better and always
yields optimal solutions. Furthermore, we intend to validate our algorithm’s
output with real-world simulations, to confirm that the algorithm operates as
efficiently and robustly as possible under demanding and dynamic UAV
conditions.

##  VI Acknowledgments

The authors would like to acknowledge the support of IRC for Sustainalbe
Energy Systems through the funded project no. INRE2328. Dr. Abido acknowledges
the support recieved from Saudi Data and AI Authority (SDAIA) and KFUPM under
SDAIA-KFUPM Joint Research Center for Artificial Intelligence no. JRC-
AIRFP-09.

## References

  * [1] M. D. Phung, C. H. Quach, T. H. Dinh, and Q. Ha, “Enhanced discrete particle swarm optimization path planning for uav vision-based surface inspection,” Automation in Construction, vol. 81, pp. 25–33, 2017.
  * [2] M. D. Phung, T. H. Dinh, Q. P. Ha, et al., “System architecture for real-time surface inspection using multiple uavs,” IEEE Systems Journal, vol. 14, no. 2, pp. 2925–2936, 2019.
  * [3] M. D. Phung and Q. P. Ha, “Motion-encoded particle swarm optimization for moving target search using uavs,” Applied Soft Computing, vol. 97, p. 106705, 2020.
  * [4] L. Lin and M. A. Goodrich, “Uav intelligent path planning for wilderness search and rescue,” in 2009 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 709–714, IEEE, 2009.
  * [5] C. Yin, Z. Xiao, X. Cao, X. Xi, P. Yang, and D. Wu, “Offline and online search: UAV multiobjective path planning under dynamic urban environment,” IEEE Internet of Things Journal, vol. 5, no. 2, pp. 546–558, 2018.
  * [6] R. W. Beard, T. W. McLain, M. A. Goodrich, and E. P. Anderson, “Coordinated target assignment and intercept for unmanned air vehicles,” IEEE transactions on robotics and automation, vol. 18, no. 6, pp. 911–922, 2002.
  * [7] T. W. McLain and R. W. Beard, “Coordination variables, coordination functions, and cooperative timing missions,” Journal of Guidance, Control, and Dynamics, vol. 28, no. 1, pp. 150–161, 2005.
  * [8] D. Eppstein, “Finding the k shortest paths,” SIAM Journal on computing, vol. 28, no. 2, pp. 652–673, 1998.
  * [9] P. O. Pettersson and P. Doherty, “Probabilistic roadmap based path planning for an autonomous unmanned helicopter,” Journal of Intelligent & Fuzzy Systems, vol. 17, no. 4, pp. 395–405, 2006.
  * [10] Y. Lin and S. Saripalli, “Sampling-based path planning for uav collision avoidance,” IEEE Transactions on Intelligent Transportation Systems, vol. 18, no. 11, pp. 3179–3192, 2017.
  * [11] J. Barraquand, B. Langlois, and J.-C. Latombe, “Numerical potential field techniques for robot path planning,” IEEE transactions on systems, man, and cybernetics, vol. 22, no. 2, pp. 224–241, 1992.
  * [12] Y.-b. Chen, G.-c. Luo, Y.-s. Mei, J.-q. Yu, and X.-l. Su, “Uav path planning using artificial potential field method updated by optimal control theory,” International Journal of Systems Science, vol. 47, no. 6, pp. 1407–1420, 2016.
  * [13] B. Di, R. Zhou, and H. Duan, “Potential field based receding horizon motion planning for centrality-aware multiple uav cooperative surveillance,” Aerospace Science and Technology, vol. 46, pp. 386–397, 2015.
  * [14] P.-C. Song, J.-S. Pan, and S.-C. Chu, “A parallel compact cuckoo search algorithm for three-dimensional path planning,” Applied Soft Computing, vol. 94, p. 106443, 2020.
  * [15] V. Roberge, M. Tarbouchi, and G. Labonté, “Comparison of parallel genetic algorithm and particle swarm optimization for real-time uav path planning,” IEEE Transactions on industrial informatics, vol. 9, no. 1, pp. 132–141, 2012.
  * [16] V. Roberge, M. Tarbouchi, and G. Labonté, “Fast genetic algorithm path planner for fixed-wing military uav using gpu,” IEEE Transactions on Aerospace and Electronic Systems, vol. 54, no. 5, pp. 2105–2117, 2018.
  * [17] Y. Fu, M. Ding, C. Zhou, and H. Hu, “Route planning for unmanned aerial vehicle (uav) on the sea using hybrid differential evolution and quantum-behaved particle swarm optimization,” IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 43, no. 6, pp. 1451–1465, 2013.
  * [18] C. Xu, H. Duan, and F. Liu, “Chaotic artificial bee colony approach to uninhabited combat air vehicle (ucav) path planning,” Aerospace science and technology, vol. 14, no. 8, pp. 535–541, 2010.
  * [19] X. Yu, W.-N. Chen, T. Gu, H. Yuan, H. Zhang, and J. Zhang, “Aco-a*: Ant colony optimization plus a* for 3-d traveling in environments with dense obstacles,” IEEE Transactions on Evolutionary Computation, vol. 23, no. 4, pp. 617–631, 2018.
  * [20] Y. Fu, M. Ding, and C. Zhou, “Phase angle-encoded and quantum-behaved particle swarm optimization applied to three-dimensional route planning for uav,” IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans, vol. 42, no. 2, pp. 511–526, 2011.
  * [21] J. Kennedy, “Swarm intelligence,” in Handbook of nature-inspired and innovative computing: integrating classical models with emerging technologies, pp. 187–219, Springer, 2006.
  * [22] Z.-L. Gaing, “Particle swarm optimization to solving the economic dispatch considering the generator constraints,” IEEE transactions on power systems, vol. 18, no. 3, pp. 1187–1195, 2003.
  * [23] P. Das and P. K. Jena, “Multi-robot path planning using improved particle swarm optimization algorithm through novel evolutionary operators,” Applied Soft Computing, vol. 92, p. 106312, 2020.
  * [24] Y. Zhang, D.-w. Gong, and J.-h. Zhang, “Robot path planning in uncertain environment using multi-objective particle swarm optimization,” Neurocomputing, vol. 103, pp. 172–185, 2013.
  * [25] Z. Wei-Min, L. Shao-Jun, and Q. Feng, “$\theta$-pso: a new strategy of particle swarm optimization,” Journal of Zhejiang University-SCIENCE A, vol. 9, no. 6, pp. 786–790, 2008.
  * [26] V. Hoang, M. D. Phung, T. H. Dinh, and Q. P. Ha, “Angle-encoded swarm optimization for uav formation path planning,” in 2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 5239–5244, IEEE, 2018.
  * [27] Y. Zhang, D.-w. Gong, X.-y. Sun, and N. Geng, “Adaptive bare-bones particle swarm optimization algorithm and its convergence analysis,” Soft Computing, vol. 18, pp. 1337–1352, 2014.
  * [28] K. Shi, Z. Wu, B. Jiang, and H. R. Karimi, “Dynamic path planning of mobile robot based on improved simulated annealing algorithm,” Journal of the Franklin Institute, January 2023.
  * [29] G. G. R. de Castro, M. F. Pinto, I. Z. Biundini, A. G. Melo, A. L. M. Marcato, and D. B. Haddad, “Dynamic path planning based on neural networks for aerial inspection,” Journal of Control, Automation and Electrical Systems, vol. 34, pp. 85–105, 2023.
  * [30] G. Fu, Y. Gao, L. Liu, M. Yang, and X. Zhu, “Uav mission path planning based on reinforcement learning in dynamic environment,” Journal of Function Spaces, 2023.
  * [31] A. Hentout, A. Maoudj, and M. Aouache, “A review of the literature on fuzzy-logic approaches for collision-free path planning of manipulator robots,” Artificial Intelligence Review, vol. 56, no. 4, pp. 3369–3444, 2023.
  * [32] J. Branke, “Evolutionary approaches to dynamic optimization problems-updated survey,” in GECCO Workshop on evolutionary algorithms for dynamic optimization problems, pp. 27–30, 2001.
  * [33] J. Branke and H. Schmeck, “Designing evolutionary algorithms for dynamic optimization problems,” Advances in evolutionary computing: theory and applications, pp. 239–262, 2003.
  * [34] J. Branke, Evolutionary optimization in dynamic environments, vol. 3.  Springer Science & Business Media, 2012.
  * [35] T. Blackwell, “Particle swarm optimization in dynamic environments,” Evolutionary computation in dynamic and uncertain environments, pp. 29–49, 2007.
  * [36] S. Alaliyat, R. Oucheikh, and I. Hameed, “Path planning in dynamic environment using particle swarm optimization algorithm,” in 2019 8th International Conference on Modeling Simulation and Applied Optimization (ICMSAO), pp. 1–5, IEEE, 2019.
  * [37] D. Parrott and X. Li, “A particle swarm model for tracking multiple peaks in a dynamic environment using speciation,” in Proceedings of the 2004 congress on evolutionary computation (IEEE Cat. No. 04TH8753), vol. 1, pp. 98–103, IEEE, 2004.
  * [38] Q. P. H. Manh Duong Phung, “Safety-enhanced uav path planning with spherical vector-based particle swarm optimization,” Applied Soft Computing 107376, vol. 107, 2021.


# Epidemic spreading and herd immunity in a driven non-equilibrium system of
strongly-interacting atoms

Dong-Sheng Ding<EMAIL_ADDRESS>Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, Anhui 230026, China.
Synergetic Innovation Center of Quantum Information and Quantum Physics,
University of Science and Technology of China, Hefei, Anhui 230026, China.
Zong-Kai Liu  Key Laboratory of Quantum Information, University of Science and
Technology of China, Hefei, Anhui 230026, China.  Synergetic Innovation Center
of Quantum Information and Quantum Physics, University of Science and
Technology of China, Hefei, Anhui 230026, China.  Hannes Busche  Department of
Physics, Chemistry and Pharmacy, Physics@SDU, University of Southern Denmark,
5230 Odense M, Denmark  Bao-Sen Shi<EMAIL_ADDRESS>Key Laboratory of
Quantum Information, University of Science and Technology of China, Hefei,
Anhui 230026, China.  Synergetic Innovation Center of Quantum Information and
Quantum Physics, University of Science and Technology of China, Hefei, Anhui
230026, China.  Guang-Can Guo  Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei, Anhui 230026, China.
Synergetic Innovation Center of Quantum Information and Quantum Physics,
University of Science and Technology of China, Hefei, Anhui 230026, China.
Charles S. Adams<EMAIL_ADDRESS>Department of Physics, Joint Quantum
Centre (JQC) Durham-Newcastle, Durham University, South Road, Durham DH1 3LE,
United Kingdom.  Franco Nori<EMAIL_ADDRESS>Theoretical Quantum Physics
Laboratory, RIKEN Cluster for Pioneering Research, Wako-shi, Saitama 351-0198,
Japan.  RIKEN Center for Quantum Computing (RQC), Wako-shi, Saitama 351-0198,
Japan.  Physics Department, The University of Michigan, Ann Arbor, MI
48109-1040,USA.

(August 30, 2024)

It is increasingly important to understand the spatial dynamics of epidemics.
While there are numerous mathematical models of epidemics, there is a scarcity
of physical systems with sufficiently well-controlled parameters to allow
quantitative model testing. It is also challenging to replicate the macro non-
equilibrium effects of complex models in microscopic systems. In this work, we
demonstrate experimentally a physics analog of epidemic spreading using
optically-driven non-equilibrium phase transitions in strongly interacting
Rydberg atoms. Using multiple laser beams we can impose any desired spatial
structure. We observe spatially localized phase transitions and their
interplay in different parts of the sample. These phase transitions simulate
the outbreak of an infectious disease in multiple locations, as well as the
dynamics towards “herd immunity” and “endemic state” in different regimes. The
reported results indicate that Rydberg systems are versatile enough to model
complex spatial-temporal dynamics.

Self-organization and non-equilibrium dynamics of complex systems emerge not
only in physics, but also in other fields such as earth science, biology, and
economics (haken2006information, ). In these cases, many interacting degrees
of freedom are repelling, such as sand and rice grains (per1987self, ),
vortices (field1995superconducting, ; olson1997superconducting, ), etc. and
they can be temporally stuck in metastable states, due to local traps or due
to static friction. When slowly driven towards marginal stability, these
extended systems exhibit avalanches at various length scales (per1987self, ;
field1995superconducting, ; olson1997superconducting, ). Rydberg atoms allow
to model spatial dynamics with emerging complexity (schauss2012observation, ;
labuhn2016tunable, ; bernien2017probing, ; de2019observation, ;
bluvstein2021controlling, ). These atoms combine the precision of atomic
physics with strong interactions with neighboring atoms, more common in
condensed matter systems (eisert2015quantum, ; browaeys2020many, ). Rydberg
atoms are not only ideal for applications in quantum technology and computing
(saffman2010quantum, ; Adams_2019, ; morgado2020quantum, ), but also as a
potential model system for simulations (bloch2012quantum, ;
georgescu2014quantum, ; morgado2020quantum, ). Precise control of the
excitation probability using narrowband lasers (gallagher2005rydberg, ;
firstenberg2016nonlinear, ) makes Rydberg atoms an excellent candidate to
study non-equilibrium physics, as instabilities of equilibrium states are
ubiquitous processes occurring in a variety of driven Rydberg systems. For
example, Rydberg atoms in different driven configurations exhibit many
fascinating characteristics of complex systems, including aggregate formation
(urvoy2015strongly, ), non-equilibrium phase transitions
(carr2013nonequilibrium, ; malossi2014full, ), critical points
(marcuzzi2014universal, ), self-organized criticality (helmrich2020signatures,
; ding2020phase, ), epidemic growth (wintermantel2021epidemic, ), and
hydrodynamics (klocke2021hydrodynamic, ).

There is considerable interest in simulating the outbreak of infectious
diseases. The accuracy of epidemic models has been considerably improved over
time (eubank2004modelling, ; longini2005containing, ). It is desirable to
study and replicate the macro non-equilibrium effects of complex dynamics in
well-controlled microscopic atomic systems. Moreover, many non-equilibrium
phenomena cannot be broken down to just two global stable states; and numerous
real systems are spatially inhomogeneous, e.g. the outbreak of an infectious
disease can occur in multiple locations (pastor2015epidemic, ). The added
interplay between multiple spatial domains results in additional complexity
that manifests as both local and non-local dynamics. Consequently, the study
of phase transitions in engineered spatially inhomogeneous systems is an
important tool to improve our understanding of non-equilibrium dynamics in
extended systems.

Figure 1: Epidemic model simulation with Rydberg avalanche excitations. (a)
The model consists of three states: the interacting Rydberg state (red), the
non-interacting Rydberg state (green) and the ground state (black). The right
diagram shows the energy levels of atoms under a laser drive. Here, $\Delta$
is defined as the effective detuning between the Rydberg state $\lvert
r\rangle$ and the laser, and $\Gamma_{r}$ is the decay rate of $\lvert
r\rangle$. (b) The algorithm of our 2D cellular automaton. The rule of
updating cells in the susceptible-infected-recovered model is shown in
(b1–b3): the cells are updated in the next step according to its current state
and the state of its eight nearest neighbors. The susceptible infected-
susceptible model satisfies the condition $\beta\gg\gamma$, $\gamma\rightarrow
0,\mu=0.01$; while the susceptible-infected-recovered process obeys
$\beta/\gamma>1$ and $\mu\rightarrow 0$. The interacting Rydberg atoms spread
out in the entire interacting region because the population-dependent energy
shift $\Delta_{\rm shift}$ induces a nonlinear ‘facilitated’ excitation, in
which the effective detuning $\Delta$ between the laser and the Rydberg energy
level of the excited state $\lvert r\rangle$ is shifted. (c) A 2D snapshot of
interacting Rydberg atoms (red) spreading in a $N=100\times 100$ inhomogeneous
Rydberg density system (with fractions $f_{R}=0$–$0.9$) after $i=200$
iterations. There is an isolated boundary of contagion (marked in white), a
domain wall separating the inhomogeneous distributions of Rydberg density.

Here we model epidemic spreading dynamics using a laser-driven thermal Rydberg
gas. We have observed the analog of typical transmission profiles for
susceptible infected-susceptible (SIS) and susceptible-infected-recovered
(SIR) processes. The intensity pattern of the laser beams creates a spatially
inhomogeneous Rydberg excitation rate, and thus spatial domains with different
Rydberg atom densities. Within each domain, the system can be found in one of
two phases. Below a critical threshold density, interactions between Rydberg
atoms are negligible; We refer to this as the ‘non-interacting’ NI phase.
Above the threshold, interactions between Rydberg atoms induce a shift and
broadening of the Rydberg lines that facilitate Rydberg excitation in adjacent
regions; This triggers a localized excitation avalanche and transitions to a
‘strongly-interacting’ I phase. We observe separate phase jumps and hysteresis
loops associated with different spatial domains within our medium. The
multiple hysteresis loops combine to produce multistability in the optical
transmission. The experimental atomic system shows an initial exponential
spread of the (microscopic) “epidemic”, until the system reaches two cases:
one is “herd immunity”, which stops the spread in the SIR process; another is
“endemic state”, which is a stationary state in the SIS case.

In order to gain insight into the spatial dynamics of Rydberg non-equilibrium
phase transitions and, in particular, the spreading of the interacting phase,
we build an avalanche model via the SIS and SIR models introduced in Refs.
(pastor2015epidemic, ). The elements of the effective model are schematically
shown in Fig. 1(a); where $\beta$ is the spreading rate from the non-
interacting Rydberg state to the interacting Rydberg state, $\mu$ is the rate
of the reverse process, and $\gamma$ is the decay rate of the interacting
Rydberg state to the ground state. The spreading process is simulated by a
cellular automaton with the rules shown in Fig. 1(b). The model uses
$N=m\times m$ 2D cells. Each cell updates according to its current state and
the state of its eight nearest neighbors. First, the cells are randomly filled
with fractions $f_{R}$ of Rydberg atoms. A cell is in the NI phase if
$f_{R}<f_{R,c}$ [green cells in Fig. 1(b)], in the I phase if $f_{R}>f_{R,c}$
(red), or in a depleted phase (black) without Rydberg atoms in the cell. These
situations correspond to a susceptible person, an infected person, or an
immune site, in the original epidemic model (pastor2015epidemic, ).

The procedure schematically shown in Fig. 1(b) describes the following steps:
if the fraction $f_{R}$ of Rydberg atoms in an arbitrary cell [red cells in
figure 1(b1)] exceeds the critical fraction $f_{R,c}$, and is thus in the
interacting phase, it increases the excitation probability of neighboring
Rydberg atoms [green cells in Fig. 1(b1)] $f_{R}$ to be above $f_{R,c}$, such
that the neighboring cells transition to the interacting phase as well in the
next iteration [red cells in Fig. 1(b2)]. This occurs due to level shifts
$\Delta_{\rm shift}$ by strong interactions that result in a ‘facilitated’
excitation process, which triggers an avalanche process throughout the
interaction region, corresponding to the epidemic spreading to adjacent sites;
see the energy diagrams in Fig. 1(b1–b3). We consider the condition
$\beta\gg\gamma,\mu$ for the following SIS simulation, which corresponds to
fast scanning either the probe intensity or $\Delta_{c}$, and the system
changes from the NI to the I phase. Otherwise, the system would oscillate near
the critical point, and the oscillations between the phases display a bimodal
distribution of transmission levels as demonstrated in Ref. (ding2020phase, ),
see also (supplement, ). The boundary condition in our model is such that
interacting Rydberg atoms at the edges disappear as they would move out of the
excitation volume defined by the laser beams.

Figure 2: Epidemic nonlinear spreading (a) SIS model by simulating avalanche
behavior of interacting Rydberg atoms via a cellular automaton in an
$N=100\times 100$ cells system. (b) Measured phase transition versus probe
intensity. The solid lines in (a) and (b) are fit by the hyperbolic tangent
function $A+B\tanh[(x-C)/\omega]$, with parameters $A=B=0.47$, $C=0.6$ and
$\omega=0.05$ for (a), and $A=0.12\text{, }B=0.11$, $C=0.45$ and
$\omega=0.035$ for (b). (c) SIR model showing a peak versus time. (d) The
measured transmission in time without scanning laser detuning. The solid lines
in (c) and (d) are fit by the Gaussian function $Be^{[-\omega(t-C)^{2}]}$ with
parameters $B=0.98$, $C=27.6$ and $\omega=0.0041$ for (c), and $B=0.034$,
$C=0.8$ and $\omega=46$ for (b). The dynamics is divided in three regimes:
initial, exponential growth, and “herd immunity” or “endemic state”. The gray
areas in (a-d) are the exponential growth regions (or nonlinear outbreak) in
the epidemic process.

We have simulated the dynamics of ‘facilitated’ Rydberg excitations according
to the cellular automaton, see Fig. 1(c). After $i=200$ iterations in a 2D
inhomogeneous Rydberg density system from fractions $f_{R}=0$ to $f_{R}=0.9$,
respectively, with $m=100$, there is an isolated boundary edge [shown in white
in Fig. 1(c)] or “domain wall” that splits up the NI and I phases. On the
contrary, there are no obvious domain walls if using a uniform atomic density.
In the latter case, we record the Rydberg atoms density in the I phase after
$200$ iterations, and the simulated result is shown in Fig. 2(a). The
interacting Rydberg fraction versus the Rydberg density $f_{R}$ shows a
critical point near a Rydberg fraction $f_{R,c}=0.55$–$0.65$. In the
experiment, we measure the transmission of the probe beam against the Rydberg
atoms fraction by increasing the probe intensity in a two-photon
electromagnetically induced transparency (EIT) scheme, see more details on the
experimental configuration in (supplement, ). There is a sudden jump in the
probe transmission spectrum given in Fig. 2(b), which corresponds to the
system transition from the NI to the I phase within the exponential growth
regime. The threshold effect in the spreading of the interacting Rydberg atoms
is consistent with the nonlinear spreading characteristics in the SIS model,
as predicted in Fig. 2(a). In the case, the system is scanned in a fast rate
(see more details in (supplement, )), the NI and I phases tends to dynamical
equilibrium, which predicts a stationary “endemic state”, a main character of
SIS model. When we consider the SIR process, the parameters $\beta$ and $\mu$
satisfy $\beta/\gamma>1$ and $\mu\rightarrow 0$, corresponding to the case of
measuring the Rydberg excitations without scanning the laser detuning in the
experiment. We have simulated the results in Fig. 2(c) with $\beta=0.95$ and
$\gamma=0.2$. The interacting Rydberg atoms fraction first increases due to
the contagion effect in the exponential regime, and then suddenly decreases
because the system reaches “herd immunity” via the “infection recovery
channel” in the SIR process. We have also measured the probe transmission
spike near the critical point in the time domain, as shown in Fig. 2(d). The
experimental Rydberg system shows the exponential spread of the (microscopic)
epidemic until it reaches “herd immunity”.

Figure 3: Multi-domain bistability. Comparison of the (a) theoretically, and
(b) experimentally observed probe transmission. The probe transmission versus
the coupling detuning $\Delta_{c}$ are shown in (c) theory and (d)
experiments. The two Rydberg fractions $f_{R1}$ and $f_{R2}$ are fixed in (c)
and the perturbing field Rabi frequency $\Omega_{\textrm{pert}}$ is fixed in
(d). The areas indicated by A and B show their individual bistability regions,
which add up to multi-stability where they overlap. The hysteresis loops of
multistability are shown in (c) theory and (d) experiment. Positive ($+$) and
negative ($-$) scan directions are shown in red and black, respectively.
Figure 4: Multiple phase jumps. Panels (a)-(d) show, respectively, spectra
with single-, double- and triple-jumps in the transmission. The different
spectra are obtained by adding the additional perturbing fields 1 and 2 to
create two and three domains. In these cases, the angles between the probe and
the two perturbing fields are 7∘ and 14∘, respectively. The dashed lines given
in (a) and (b) are the slopes of the jumps. Here,
$\Omega_{\textrm{pert,1}}=2\pi\times 3.7$ MHz and
$\Omega_{\textrm{pert,2}}=2\pi\times 2.2$ MHz. Panels (e)-(h) show our
simulations using the SIS cellular automaton (drossel1992self, ;
ding2020phase, ) with one, two, and three spatial domains, respectively. The
solid lines in these simulations are fit by the function $A+\sum
B_{i}\tanh((x-C_{i})/\omega_{i})$ including a series of tanh steps with
different position $C_{i}$ and width $\omega_{i}$. We initialize the domain
with atomic fractions $f_{R1}=f_{R}$ in (e) and $f_{R1}=f_{R}+0.2$,
$f_{R2}=f_{R}+0.3$ in (f), $f_{R1}=f_{R}+0.2$ and $f_{R3}=f_{R}$ in (g) and
$f_{R2}=f_{R}+0.3$, $f_{R1}=f_{R}+0.15$, $f_{R3}=f_{R}$ in (h). The insets
show the spread of interacting atoms (red) from left to right following 100
iterations.

The nonlinear spreading of the I phase is strongly dependent on the initial
Rydberg density $f_{R}$($t$=0). This reveals that one uniform spatial domain
leads to one jump, and multiple spatial domains lead to multiple jumps. This
corresponds to the scenario of an outbreak of an infectious disease in
different locations. In the experiment, a spatial inhomogeneity is introduced
by subjecting part of the atomic vapor to localized perturbing beams. The
additional perturbing fields have the same frequency as the optical probe
field, but intersect the latter at an angle, creating regions with different
intensities, as explained in (supplement, ). Note that, a domain boundary can
also arise from the attenuation of the probe light along with the medium
(carr2013nonequilibrium, ). For any spatially varying intensity distribution,
the distinct phases are characterized by a spatial boundary between the states
of high and low population. The occurrence of optical multistability and non-
equilibrium phase transitions in a multi-domain system can be simulated using
an EIT master equation with a mean-field model, see (supplement, ). The
multiple discrete transmission levels arise as the transition threshold is
reached for different detunings in each domain.

In order to simulate the non-equilibrium behavior of the spatial domains in
our phenomenological model, we consider distinct fractions $f_{R1}$ and
$f_{R2}$ to each domain and their overall optical response. We then plot a
color map of the transmission difference for the $+/-$ scan directions against
$f_{\text{R2}}$ and $\Delta_{c}$ by setting $f_{\text{r1}}=0.33$, as in Fig.
3(a). For the experiments, we extract the transmission difference for the
$+/-$ scan directions and plot the probe transmission versus
$\Omega_{\textrm{pert}}$ and $\Delta_{c}$ in Fig. 3(b). In both theory and
experiments, we observe two bistable regimes, and a multistable regime with
overlapping regions. We also plot the transmission for the $+/-$ scan
directions for theory and experiments, as shown in Figs. 3(c) and (d). The
optical responses for multiple spatial domains show a double optical
bistability, as described in (supplement, ).

It is insightful to increase the number of spatial domains in the experiments,
and hence the number of phase jumps, by adding a few perturbing laser beams.
Figures 4(a-d) show the transmission spectra with 1, 2, and 3 discrete jumps,
obtained by scanning $\Delta_{c}$ from red- to blue-detuning. Figure 4(a)
refers to the single domain created by the probe field, showing a sudden jump
in the transmission of light [gray region (1), where $\Delta_{c}/2\pi\sim-19$
MHz]. As we add the perturbing field 1 overlapping with the probe, the jump in
Fig. 4(a) is split into two jumps in Fig. 4(b), these located at the detuning
$\Delta_{c}/2\pi\sim-19$ MHz and $\Delta_{c}/2\pi\sim-24$ MHz, respectively.
The three jumps in Fig. 4(d) correspond to three regions of higher Rydberg
density associated with the probe and the perturbing fields 1 and 2. This
added interplay between multiple spatial domains results in additional
complexity that manifests in the additional jumps in Fig. 4(b). Meanwhile, the
phase transition is weakened, as seen from the decreased height of each jump
and the susceptibility of the phase transition defined as
$d\textrm{T}/d\Delta_{c}$ ($d\textrm{T}$ is the transmission difference of the
jump) is reduced near the transition point $\Delta_{c}/2\pi\sim-19$ MHz.
Figure 4(c) shows the corresponding situation when applying a perturbing field
2. In these two cases, as $\Delta_{c}$ approaches resonance, $f_{R}$ increases
and eventually one domain reaches the threshold $f_{R}^{(1)}$ at position 1
and another domain reaches the threshold $f_{R}^{(2)}$, the system undergoes
two non-equilibrium phase transitions accompanied by two sudden jumps in the
optical transmission. Figure 4(d) shows three jumps created by both the probe
and perturbing fields 1 and 2. Changing the alignment determines whether these
regions overlap or are separate. Here, we extend the epidemic model to
multiple spatial domains according to the 2D SIS cellular automaton. The
qualitative agreement between the predicted phase jumps, Figs. 4(e-h), and the
experimental observations, Figs. 4(a-d), confirms that we can associate each
jump with an individual domain.

Non-equilibrium phase transitions within multi-domains often correspond to
optical multistability, as demonstrated here, which is the extension of
optical bistability to three or more stable output states. Before this work,
multistability has been predicted and investigated in coupled atom-cavity
systems (kitano1981optical, ; cecchi1982observation, ; joshi2003optical, ;
sheng2012realization, ) and semiconductor microcavities
(gippius2007polarization, ; paraiso2010multistability, ;
sarkar2010polarization, ; cerna2013ultrafast, ; goblot2019nonlinear, ).
Optical bistability in Rydberg atoms has previously been studied both in
theory (lee2012collective, ; marcuzzi2014universal, ; weimer2015variational, ;
vsibalic2016driven, ; levi2016quantum, ; vsibalic2016driven, ) and experiments
(carr2013nonequilibrium, ; malossi2014full, ; de2016intrinsic, ;
weller2016charge, ; ding2020phase, ). However, optical multistability in
Rydberg atoms is more challenging, and has not been explored before. A mean-
field master Bloch equation can be used to qualitatively simulate
multistability, because the individual Rydberg atoms interact with each other
at random, in such a way that each Rydberg atom in a compartment is treated
similarly and indistinguishably. More precise simulations in large-scale
networks must consider the subtle differences between various states of the
Rydberg atoms.

In summary, we have studied optically-driven multi-domain non-equilibrium
dynamics in strongly interacting Rydberg atom gases. The measured hysteresis
and phase jumps can be understood well and reproduced qualitatively using the
epidemic model and using a mean-field master equation. Our results highlight
the rich range of non-equilibrium phenomena that are accessible even in a
relatively simple experiment and provide observational data to benchmark
theoretical models of non-equilibrium dynamics under arbitrary spatial
structures. Specifically, the observed dynamics of Rydberg atoms in different
time and space scales could predict “herd immunity”, “endemic state” and the
outbreak of a disease in multiple location. The reported multi-domain
dependent exotic phases could help build a Rydberg-based simulate platform for
more complex phenomena in real-world, such as in epidemiology
(pastor2015epidemic, ), ecosystems and climate (haken2006information, ) or
other complex systems (morgado2020quantum, ).

###### Acknowledgements.

We thank for Christopher G. Wade, Kevin J. Weatherill and Igor Lesanovsky for
helpful discussions on the phase boundaries and non-equilibrium dynamics. We
acknowledge funding from National Key R&D Program of China (2017YFA0304800),
NSFC funding (Grant Nos. U20A20218, 61525504, 61722510, 61435011), the Youth
Innovation Promotion Association of CAS Grant No. 2018490, EPSRC through grant
agreements EP/M014398/1, EP/R002061/1, EP/L023024/1, EP/P012000/1,
EP/R035482/1, EP/S015973/1, as well as, DSTL, and Durham University. The
European Union’s Horizon 2020 Research and Innovation Program under Grant No.
845218 (Marie Curie Fellowship to H. B.). F.N. was supported in part by NTT
Research, JSPS Grant No. JP20H00134, ARO Grant No. W911NF-18-1-0358, AOARD
Grant No. FA2386-20-1-4069, and FQXi Grant No. FQXi-IAF19-06.

## References

  * (1) H. Haken, Information and self-organization: A macroscopic approach to complex systems (Springer Science & Business Media, 2006).
  * (2) P. Bak, C. Tang, K. Wiesenfeld, Self-organized criticality: and explanation of 1/f noise.  Phys. Rev. Lett. 59, 381–384 (1987).
  * (3) S. Field, J. Witt, F. Nori, X. Ling, Superconducting vortex avalanches.  Phys. Rev. Lett. 74, 1206 (1995).
  * (4) C. Olson, C. Reichhardt, F. Nori, Superconducting vortex avalanches, voltage bursts, and vortex plastic flow.  Phys. Rev. B. 56, 6175 (1997).
  * (5) P. Schauß, et al., Observation of spatially ordered structures in a two-dimensional Rydberg gas.  Nature 491, 87–91 (2012).
  * (6) H. Labuhn, et al., Tunable two-dimensional arrays of single Rydberg atoms for realizing quantum Ising models.  Nature 534, 667–684 (2016).
  * (7) H. Bernien, et al., Probing many-body dynamics on a 51-atom quantum simulator.  Nature 551, 579–584 (2017).
  * (8) S. de Léséleuc, et al., Observation of a symmetry-protected topological phase of interacting bosons with Rydberg atoms.  Science 365, 775–780 (2019).
  * (9) D. Bluvstein, et al., Controlling quantum many-body dynamics in driven Rydberg atom arrays.  Science 371, 1355–1359 (2021).
  * (10) J. Eisert, M. Friesdorf, C. Gogolin, Quantum many-body systems out of equilibrium.  Nature Physics 11, 124–130 (2015).
  * (11) A. Browaeys, T. Lahaye, Many-body physics with individually controlled Rydberg atoms.  Nature Physics 16, 132–142 (2020).
  * (12) M. Saffman, T. Walker, K. Mølmer, Quantum information with Rydberg atoms.  Reviews of Modern Physics 82, 2313 (2010).
  * (13) C. S. Adams, J. D. Pritchard, J. P. Shaffer, Rydberg atom quantum technologies.  Journal of Physics B 53, 012002 (2019).
  * (14) M. Morgado, S. Whitlock, Quantum simulation and computing with Rydberg-interacting qubits.  AVS Quantum Sci 3, 023501 (2021).
  * (15) I. Bloch, J. Dalibard, S. Nascimbene, Quantum simulations with ultracold quantum gases.  Nature Physics 8, 267–276 (2012).
  * (16) I. M. Georgescu, S. Ashhab, F. Nori, Quantum simulation.  Reviews of Modern Physics 86, 153 (2014).
  * (17) T. F. Gallagher, Rydberg atoms, vol. 3 (Cambridge University Press, 2005).
  * (18) O. Firstenberg, C. S. Adams, S. Hofferberth, Nonlinear quantum optics mediated by Rydberg interactions.  Journal of Physics B 49, 152003 (2016).
  * (19) A. Urvoy, et al., Strongly correlated growth of Rydberg aggregates in a vapor cell.  Phys. Rev. Lett. 114, 203002 (2015).
  * (20) C. Carr, R. Ritter, C. Wade, C. S. Adams, K. J. Weatherill, Nonequilibrium phase transition in a dilute Rydberg ensemble.  Phys. Rev. Lett. 111, 113901 (2013).
  * (21) N. Malossi, et al., Full counting statistics and phase diagram of a dissipative Rydberg gas.  Phys. Rev. Lett. 113, 023006 (2014).
  * (22) M. Marcuzzi, E. Levi, S. Diehl, J. P. Garrahan, I. Lesanovsky, Universal nonequilibrium properties of dissipative Rydberg gases.  Phys. Rev. Lett. 113, 210401 (2014).
  * (23) S. Helmrich, et al., Signatures of self-organized criticality in an ultracold atomic gas.  Nature 577, 481–486 (2020).
  * (24) D.-S. Ding, H. Busche, B.-S. Shi, G.-C. Guo, C. S. Adams, Phase diagram and self-organizing dynamics in a thermal ensemble of strongly interacting Rydberg atoms.  Phys. Rev. X. 10, 021023 (2020).
  * (25) T. Wintermantel, et al., Epidemic growth and Griffiths effects on an emergent network of excited atoms.  Nature Communications 12, 1–6 (2021).
  * (26) K. Klocke, T. Wintermantel, G. Lochead, S. Whitlock, M. Buchhold, Hydrodynamic stabilization of self-organized criticality in a driven Rydberg gas.  Phys. Rev. Lett. 126, 123401 (2021).
  * (27) S. Eubank, et al., Modelling disease outbreaks in realistic urban social networks.  Nature 429, 180–184 (2004).
  * (28) I. M. Longini, et al., Containing pandemic influenza at the source.  Science 309, 1083–1087 (2005).
  * (29) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, A. Vespignani, Epidemic processes in complex networks.  Reviews of modern physics 87, 925 (2015).
  * (30) Supplementary materials, including experimental setup, optical multistability, theoretical analysis and fast and slow scan effect.
  * (31) B. Drossel, F. Schwabl, Self-organized critical forest-fire model.  Phys. Rev. Lett. 69, 1629 (1992).
  * (32) M. Kitano, T. Yabuzaki, T. Ogawa, Optical tristability.  Phys. Rev. Lett. 46, 926 (1981).
  * (33) S. Cecchi, G. Giusfredi, E. Petriella, P. Salieri, Observation of optical tristability in sodium vapors.  Phys. Rev. Lett. 49, 1928 (1982).
  * (34) A. Joshi, M. Xiao, Optical multistability in three-level atoms inside an optical ring cavity.  Phys. Rev. Lett. 91, 143904 (2003).
  * (35) J. Sheng, U. Khadka, M. Xiao, Realization of all-optical multistate switching in an atomic coherent medium.  Phys. Rev. Lett. 109, 223906 (2012).
  * (36) N. Gippius, et al., Polarization multistability of cavity polaritons.  Phys. Rev. Lett. 98, 236401 (2007).
  * (37) T. Paraïso, M. Wouters, Y. Léger, F. Morier-Genoud, B. Deveaud-Plédran, Multistability of a coherent spin ensemble in a semiconductor microcavity.  Nature Materials 9, 655–660 (2010).
  * (38) D. Sarkar, et al., Polarization bistability and resultant spin rings in semiconductor microcavities.  Phys. Rev. Lett. 105, 216402 (2010).
  * (39) R. Cerna, et al., Ultrafast tristable spin memory of a coherent polariton gas.  Nature Communications 4, 2008 (2013).
  * (40) V. Goblot, et al., Nonlinear polariton fluids in a flatband reveal discrete gap solitons.  Phys. Rev. Lett. 123, 113901 (2019).
  * (41) T. E. Lee, H. Haeffner, M. Cross, Collective quantum jumps of Rydberg atoms.  Phys. Rev. Lett. 108, 023602 (2012).
  * (42) H. Weimer, Variational principle for steady states of dissipative quantum many-body systems.  Phys. Rev. Lett. 114, 040402 (2015).
  * (43) N. Šibalić, C. G. Wade, C. S. Adams, K. J. Weatherill, T. Pohl, Driven-dissipative many-body systems with mixed power-law interactions: Bistabilities and temperature-driven nonequilibrium phase transitions.  Phys. Rev. A. 94, 011401 (2016).
  * (44) E. Levi, R. Gutiérrez, I. Lesanovsky, Quantum non-equilibrium dynamics of Rydberg gases in the presence of dephasing noise of different strengths.  Journal of Physics B 49, 184003 (2016).
  * (45) N. R. de Melo, et al., Intrinsic optical bistability in a strongly driven Rydberg ensemble.  Phys. Rev. A. 93, 063863 (2016).
  * (46) D. Weller, A. Urvoy, A. Rico, R. Löw, H. Kübler, Charge-induced optical bistability in thermal Rydberg vapor.  Phys. Rev. A. 94, 063820 (2016).


# Explicit Constructions of Finite Groups as Monodromy Groups

Ra-Zakee Muhammad, Javier Santiago and Eyob Tsegaye  Department of
Mathematics, Pomona College, 333 N. College Way Claremont, CA 91711
<EMAIL_ADDRESS>Department of Computer Science, University of
Puerto Rico, Río Piedras, 17 Ave Universidad, Ste 1701, San Juan, Puerto Rico
00925-2537<EMAIL_ADDRESS>Deparment of Mathematics, Stanford
University, Building 380, Stanford, California 94305<EMAIL_ADDRESS>
###### Abstract.

In 1963, Greenberg proved that every finite group appears as the monodromy
group of some morphism of Riemann surfaces. In this paper, we give two
constructive proofs of Greenberg’s result. First, we utilize free groups,
which given with the universal property and their construction as discrete
subgroups of $\operatorname{PSL}_{2}(\mathbb{R})$, yield a very natural
realization of finite groups as monodromy groups. We also give a proof of
Greenberg’s result based on triangle groups $\Delta(m,n,k)$. Given any finite
group $G$, we make use of subgroups of $\Delta(m,n,k)$ in order to explicitly
find a morphism $\pi$ such that $G\simeq\operatorname{Mon}(\pi)$.

##  1\. Introduction

Greenberg’s celebrated result [8] from 1963 shows than any finite group can be
realized as the monodromy group of some morphism of Riemann surfaces. Over
$40$ years later, in [13], Wolfart provided a method to make Greenberg’s ideas
explicit for finite cyclic groups. Specifically, he utilizes surjections from
triangle groups in order to construct the cyclic group as a monodromy group.
Inspired by Wolfart, in this paper we seek to make Greenberg’s result explicit
for all finite groups. That is, given any finite group $G$, our goal is to
find an explicit morphism $\pi$ with monodromy group $\operatorname{Mon}(\pi)$
such that $G\simeq\operatorname{Mon}(\pi)$. On that account, we will be
presenting two constructive proofs of Greenberg’s classical result.

We first present a natural construction arising from the theory of free
groups, where we make use of the appearance of all groups $G$ as quotient
groups of free groups [5], as well as the identification of the free group of
rank $2$ with the congruence subgroup $\Gamma(2)$ [7]. Then we give a proof of
our main result, which is a construction which utilizes surjections from
appropriate subgroups of triangle groups. The formal construction is as
follows:

######  Theorem 1.1.

Let $S_{n}$ be the symmetric group of degree $n$ and let
$\varphi\colon\Delta(2,n,n-1)\to S_{n}$ be the homomorphism given by

$\gamma_{0}\mapsto(1\;2),\quad\gamma_{1}\mapsto(1\;2\;\cdots\;n),\quad\gamma_{\infty}\mapsto(n\;\cdots\;3\;2).$

Then,

  1. (i)

$\varphi$ is a surjective group homomorphism. Moreover, if $n<5$, then
$\varphi$ is an isomorphism.

  2. (ii)

If $G$ is a subgroup of $S_{n}$, then $\Gamma=\varphi^{-1}(G)$ and
$\Gamma^{\prime}=\operatorname{ker}(\varphi)$ are both Fuchsian groups.

  3. (iii)

$G\simeq\operatorname{Mon}(\pi)$, where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$ is
the canonical projection and $G$, $\Gamma$, and $\Gamma^{\prime}$ are as
stated in (ii).

Note that by Cayley’s theorem, we may identify any finite group $G$ as a
subgroup of $S_{n}$ for appropriate $n$. The organization of the paper is the
following. In Section 2, we define the concepts of Fuchsian groups and
triangle groups, in addition to giving a description of monodromy in simple
algebraic terms. Moreover, we also go into detail on the explicit
constructions of free groups $F_{r}$ as subgroups of $\Gamma(2)$. In Section
3, we then give both constructive proofs of Greenberg’s result.

##  2\. Preliminaries

Recall $\operatorname{SL}_{2}(\mathbb{R})$, as a topological space, can be
naturally identified with a subset of $\mathbb{R}^{4}$. Consequently,
$\operatorname{PSL}_{2}(\mathbb{R})$ is endowed with a quotient topology. In
this paper, we will be interested in subgroups
$\Gamma\leq\operatorname{PSL}_{2}(\mathbb{R})$ such that the topology induced
in $\Gamma$ is discrete.

######  Definition 2.1.

[7] A Fuchsian group is a subgroup
$\Gamma\leq\operatorname{PSL}_{2}(\mathbb{R})$ such that for all
$\gamma\in\Gamma$ there is a neighbourhood $V$ of $\gamma$ in
$\operatorname{PSL}_{2}(\mathbb{R})$ with $V\cap\Gamma=\\{\gamma\\}$.

By definition, it is clear that if $\Gamma$ is a Fuchsian group and
$\Gamma^{\prime}\leq\Gamma$, then $\Gamma^{\prime}$ is also a Fuchsian group.
This will be relevant when working with subgroups of
$\operatorname{PSL}_{2}(\mathbb{Z})$, since
$\operatorname{PSL}_{2}(\mathbb{Z})$ is a discrete subgroup of
$\operatorname{PSL}_{2}(\mathbb{R})$. For more details regarding Fuchsian
groups, see [10, 9, 6].

As mentioned before, the goal of this paper is to explicitly realize any
finite group $G$ as the monodromy group of some branched cover. Now, on one
hand, given a branched cover of Riemann surfaces, we may realize the monodromy
group by lifting paths around the branch points. However, we may also take an
algebraic point of view. In [7], we see that monodromy groups can also be
explicitly constructed from Fuchsian groups. Namely, given two Fuchsian groups
$\Gamma^{\prime},\Gamma$ with $\Gamma^{\prime}\leq\Gamma$, we have
$\displaystyle{\frac{\Gamma}{\bigcap_{\gamma\in\Gamma}\gamma^{-1}\Gamma^{\prime}\gamma}\simeq\operatorname{Mon}(\pi)}$,
where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$.
Furthermore,
$\displaystyle{\bigcap_{\gamma\in\Gamma}\gamma^{-1}\Gamma^{\prime}\gamma}$ is
the largest normal subgroup of $\Gamma$ that is contained in $\Gamma^{\prime}$
(also called the normal core). Therefore, if $\Gamma^{\prime}$ is a normal
subgroup of $\Gamma$, we have the following simple description of
$\operatorname{Mon}(\pi)$.

######  Corollary 2.2.

[7, Corollary 2.59] Let $\Gamma^{\prime},\Gamma$ be Fuchsian groups. If
$\Gamma^{\prime}\trianglelefteq\Gamma$, then
$\Gamma/\Gamma^{\prime}\simeq\operatorname{Mon}(\pi)$, where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$.

Making Greenberg’s result explicit is then a purely algebraic question: Given
any finite group $G$, can we find Fuchsian groups $\Gamma^{\prime},\Gamma$
with $\Gamma^{\prime}\trianglelefteq\Gamma$ such that
$G\simeq\Gamma/\Gamma^{\prime}$? Theorem 3.1 and Theorem 1.1 provide two
separate constructions that answer this question.

In Theorem 3.1, we take $\Gamma$ to be a free group, whose rank depends on the
number of generators of the finite group $G$. This construction is based on
both the fact that $\Gamma(2)$ is a free group generated by
$A=\pm\begin{bmatrix}1&2\\\ 0&1\end{bmatrix},B=\pm\begin{bmatrix}1&0\\\
2&1\end{bmatrix}\in\operatorname{PSL}_{2}(\mathbb{Z})$ [2], as well as the
classical result that the free group of rank two contains freely generated
subgroups of rank $r$ for any $r\in\mathbb{N}$. The following lemma gives an
explicit construction of free groups of rank $r$ as subgroups of $\Gamma(2)$.

######  Lemma 2.3.

[4, Proposition 7.80] Let $A=\pm\begin{bmatrix}1&2\\\ 0&1\end{bmatrix}$ and
$B=\pm\begin{bmatrix}1&0\\\ 2&1\end{bmatrix}$ and let
$F_{r}=\langle\;B^{-j}AB^{j}\;|\;j=0,1,\dots,r-1\;\rangle$. Then $F_{r}$ is a
free subgroup of rank $r$ in $\Gamma(2)$.

For Theorem 1.1 we choose $\Gamma$ to be a partiular subgroup of a triangle
group.

######  Definition 2.4.

Let $m,n,k\in\mathbb{N}$. Then a group with presentation

$\Delta(m,n,k)=\langle\gamma_{0},\gamma_{1},\gamma_{\infty}\;|\;\gamma_{0}^{m}=\gamma_{1}^{n}=\gamma_{\infty}^{k}=\gamma_{0}\gamma_{1}\gamma_{\infty}=1\rangle$

is called a triangle group.

An important result regarding triangle groups is the following:

######  Theorem 2.5.

[1] Let $m,n,k\in\mathbb{N}$. Then the triangle group $\Delta(m,n,k)$ is a
Fuchsian group.

######  Remark 2.6.

[7] Notice $\Gamma(2)$ can be identified with the triangle group
$\Delta(\infty,\infty,\infty)$.

In both Theorem 3.1 and Theorem 1.1, we will construct a surjective
homomorphism $\varphi$ from $\Gamma$ onto the finite group $G$, and
$\operatorname{ker}(\varphi)$ will be the natural choice for the normal
subgroup $\Gamma^{\prime}$.

##  3\. Constructions

In this section, we will be using Corollary 2.2 to explicitly realize finite
groups $G$ as monodromy groups. As discussed in Section 2, given a finite
group $G$, our approach is to find a Fuchsian group $\Gamma$ and a surjective
homomorphism $\varphi\colon\Gamma\to G$ such that
$G\simeq\Gamma/\Gamma^{\prime}$, where
$\Gamma^{\prime}=\operatorname{ker}(\varphi)$. The kernel of such a
homomorphism $\varphi$ is then a normal Fuchsian subgroup of $\Gamma$. As a
result, we can then use Corollary 2.2 to conclude
$G\simeq\operatorname{Mon}(\pi)$, where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$ is
the canonical map. Two separate constructions will be provided, one with free
groups and one with subgroups of triangle groups.

###  3.1. Free group construction

From group theory, we have that any group $G$ is isomorphic to a quotient
group of some free groups. Here, we will be using the structure of groups as
quotients of free groups in order to explicitly express finite groups as
monodromy groups.

Consider both $\Gamma(2)=\langle A,B\rangle$ and $F_{r}$ as constructed in
Lemma 2.3.

######  Theorem 3.1.

Let $G$ be an $r$-generated finite group, $\Gamma=F_{r}$ and $\varphi$ a
surjective group homomorphism $\varphi\colon\Gamma\to G$ with
$\Gamma^{\prime}=\ker\varphi$. Then,

  1. (i)

$\Gamma$ and $\Gamma^{\prime}$ are both Fuchsian groups. Moreover,
$\Gamma^{\prime}$ is a free group of rank $1+\lvert G\rvert(r-1)$.

  2. (ii)

$G\simeq\operatorname{Mon}(\pi)$, where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$ is
the canonical projection.

###### Proof.

$(i)$ By the construction in Lemma 2.3, we see that
$\Gamma\leq\operatorname{PSL}_{2}(\mathbb{Z})$, and thus, it is a Fuchsian
group. Since $\Gamma^{\prime}=\operatorname{ker}(\varphi)$ is a subgroup of
the Fuchsian group $\Gamma$, it then follows that $\Gamma^{\prime}$ is also a
Fuchsian group. Notice further, by the Schreier index formula [11], the rank
of $\Gamma^{\prime}$ is $1+\lvert G\rvert(r-1)$.
$(ii)$ Since $\Gamma^{\prime},\Gamma$ are Fuchsian groups, we then have the
canonical map between Riemann surfaces
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$.
Consequently, given that $\varphi$ is surjective, we have
$G\simeq\Gamma/\Gamma^{\prime}$, and from Corollary 2.2 we can conclude
$G\simeq\operatorname{Mon}(\pi)$. ∎

######  Remark 3.2.

Given the group presentation for $G$, then $\operatorname{ker}(\varphi)$ can
be seen as smallest normal subgroup containing the relations in $G$.

Recall that the free group $\Gamma(2)$ can be identified with the triangle
group $\Delta(\infty,\infty,\infty)$. Thus, in Theorem 3.1, we can view
$\Gamma,\Gamma^{\prime}$ as subgroups of $\Delta(\infty,\infty,\infty)$. We
now give another proof of Greenberg’s result, this time using subgroups of the
triangle group $\Delta(m,n,k)$ with $m,n,k$ finite.

###  3.2. Triangle group construction

The specific triangle group we will be working with is $\Delta(2,n,n-1)$. As
stated in Theorem 1.1, let

$\varphi\colon\Delta(2,n,n-1)\to S_{n}$

be the homomorphism given by

$\gamma_{0}\mapsto(1\;2),\quad\gamma_{1}\mapsto(1\;2\;\cdots\;n),\quad\gamma_{\infty}\mapsto(n\;\cdots\;3\;2).$

###### Proof of Theorem 1.1.

$(i)$ Since $S_{n}$ is generated by the transposition $(1\textbf{ }2)$ and the
$n$-cycle $(1\textbf{ }2\textbf{ }\cdots\text{ }n)$ and these are mapped to by
$\gamma_{0}$ and $\gamma_{1}$ respectively, $\varphi$ is a surjective map.
Now, notice if $n\geq 5$, then $\frac{1}{2}+\frac{1}{n}+\frac{1}{n-1}<1$. In
that case, $\Delta(2,n,n-1)$ corresponds to a tiling of the hyperbolic plane.
And so, it follows that $\Delta(2,n,n-1)$ is an infinite group [12]. For the
cases with $n<5$ we go one by one. If $n=2$, then $\Delta(2,2,1)$ is the
cyclic group $C_{2}\simeq S_{2}$. If $n=3$, then $\Delta(2,3,2)$ is the
dihedral group $D_{3}\simeq S_{3}$. And finally, if $n=4$, then
$\Delta(2,4,3)$ is the octahedral group, which is isomorphic to $S_{4}$ [3].
$(ii)$ Since $\Gamma=\varphi^{-1}(G)$ and
$\Gamma^{\prime}=\operatorname{ker}(\varphi)$ are subgroups of
$\Delta(2,n,n-1)$, by Theorem 2.5 both $\Gamma,\Gamma^{\prime}$ are Fuchsian
groups.
$(iii)$ Consider the restriction

$\varphi\Big{|}_{\Gamma}\colon\Gamma\to G.$

By the first isomorphism theorem,
$G\simeq\Gamma/\operatorname{ker}\left(\varphi\Big{|}_{\Gamma}\right)$, and
since
$\Gamma^{\prime}=\ker\varphi=\operatorname{ker}\left(\varphi\Big{|}_{\Gamma}\right)$,
we have $G\simeq\Gamma/\Gamma^{\prime}$. Consequently, given that both
$\Gamma,\Gamma^{\prime}$ are Fuchsian groups, by Corollary 2.2, it then
follows that

$G\simeq\Gamma/\Gamma^{\prime}\simeq\operatorname{Mon}(\pi),$

where
$\pi\colon\Gamma^{\prime}\backslash\mathbb{H}\to\Gamma\backslash\mathbb{H}$. ∎

######  Remark 3.3.

One may observe that the homomorphism $\varphi$ in Section 3.2 is really a
projection map $F_{2}/N\to F_{2}/K$, where $N$ and $K$ are the normal closures
of the relations defining $\Delta(2,n,n-1)$ and $S_{n}$, respectively. Hence,
the kernel of this map can be explicitly computed as $K/N$.

##  4\. Acknowledgements

The authors would like to thank Dr. Edray Goins, Dr. Duane Cooper and Dr.
Alexander Barrios for providing support and guidance during our research
project. We would also like to thank the National Science Foundation Grant No.
DMS-1659138, the Alfred P. Sloan Foundation, Grant No. G-2017-9876 and the
Mathematical Sciences Research Institute (MSRI) for supporting our research.

## References

  * [1] Alan F. Beardon.  The geometry of discrete groups, volume 91.  Springer Science & Business Media, 2012.
  * [2] Bryan Birch.  Noncongruence subgroups, covers and drawings.  The Grothendieck Theory of Dessins d’Enfants, 200:25, 1994.
  * [3] Pete Clark and John Voight.  Algebraic curves uniformized by congruence subgroups of triangle groups.  Transactions of the American Mathematical Society, 371(1):33–82, 2019.
  * [4] Cornelia Druţu and Michael Kapovich.  Geometric Group Theory.  American Mathematical Society, 2018.
  * [5] David S. Dummit and Richard M. Foote.  Abstract algebra.  John Wiley & Sons, Inc., Hoboken, NJ, third edition, 2004.
  * [6] Emilio Bujalance Garcia, John William Scott Cassels, Nigel James Hitchin, et al.  Topics on Riemann surfaces and Fuchsian groups, volume 287.  Cambridge University Press, 2001.
  * [7] Ernesto Girondo and Gabino González-Diez.  Introduction to compact Riemann surfaces and dessins d’enfants, volume 79.  Cambridge University Press, 2012.
  * [8] Leon Greenberg.  Maximal fuchsian groups.  Bulletin of the American Mathematical Society, 69(4):569–573, 1963\.
  * [9] Gareth A Jones and David Singerman.  Complex functions: an algebraic and geometric viewpoint.  Cambridge university press, 1987.
  * [10] Svetlana Katok.  Fuchsian groups.  University of Chicago press, 1992.
  * [11] Jean-Pierre Serre.  Trees, 1980.
  * [12] Jerry Shurman.  Geometry of the Quintic.  John Wiley & Sons, 1997.
  * [13] Jürgen Wolfart.  ABC for Polynomials, Dessins d’Enfants, and Uniformization: a Survey.  Universitätsbibliothek Johann Christian Senckenberg, 2005.