# Team Coordination on Graphs: Problem, Analysis, and Algorithms

Manshi Limbu, Yanlin Zhou, Gregory Stein, Xuan Wang, Daigo Shishika, and Xuesu
Xiao George Mason University {klimbu2, yzhou30, gjstein, xwang64, dshishik,
<EMAIL_ADDRESS>
###### Abstract

Team Coordination on Graphs with Risky Edges (tcgre) is a recently emerged
problem, in which a robot team collectively reduces graph traversal cost
through support from one robot to another when the latter traverses a risky
edge. Resembling the traditional Multi-Agent Path Finding (mapf) problem, both
classical and learning-based methods have been proposed to solve tcgre,
however, they lacked either computation efficiency or optimality assurance. In
this paper, we reformulate tcgre as a constrained optimization and perform
rigorous mathematical analysis. Our theoretical analysis shows the NP-hardness
of tcgre by reduction from the Maximum 3D Matching problem and that efficient
decomposition is a key to tackle this combinatorial optimization problem.
Further more, we design three classes of algorithms to solve tcgre, i.e.,
Joint State Graph (jsg) based, coordination based, and receding-horizon sub-
team based solutions. Each of these proposed algorithms enjoy different
provable optimality and efficiency characteristics that are demonstrated in
our extensive experiments.

##  I Introduction

Multi-Agent Path Finding (mapf) is a trending problem in robotics [1, 2, 3, 4,
5, 6], as it lies in the core of many robotic applications, e.g., drone swarm
control [2], autonomous delivery [7], and public transportation scheduling
[1]. Adding the possibility of team coordination [8] between robots to mapf
makes the problem more difficult. The need of coordination behaviors on large-
scale multi-robot planning problems may exceed the computation capability of a
centralized planner [9], giving rise to decentralized planning that
distributes the computation into each robot [10]. However, despite the
efficiency, dexterity, and responsiveness, the distribution itself may induce
certain performance degradation and lead to suboptimal team coordination
behaviors. Therefore, some centralized pre-planning is still crucial for
multi-robot planning, especially for large-scale problems.

Team Coordination on Graphs with Risky Edges (tcgre) [11] is such a
centralized planning problem in an environment represented as a graph—multiple
robots travel from their start to goal nodes with possible support from some
nodes to reduce the cost of traversing certain risky edges, requiring team
coordination behaviors to reduce the total cost of team graph traversal. By
converting the environment graph to a Joint State Graph (jsg), optimal
coordination can be solved using Dijkstra’s search algorithm on the jsg [11].
However, the conversion to jsg does not scale well with large environment
graphs and number of robots. To address the curse of dimensionality, a
Critical Joint State Graph (cjsg) approach has been proposed for large graphs
with a small amount of support with up to two robots, still assuring solution
optimality. Reinforcement Learning (RL) has been utilized [12] to reduce the
time complexity and scale the solution to a large group of robots and the size
of the graph, but at the cost of sacrificing optimality.

To acquire theoretical insights into this problem, we reformulate tcgre in a
constrained optimization framework and present rigorous mathematical analysis
of this reformulated problem. We prove the NP-hardness of tcgre by reduction
from the Maximum 3D Matching problem. We further show that such a difficult
combinatorial optimization problem can be effectively addressed by efficient
decomposition. In addition to providing a theoretical explanation for previous
algorithms, we further introduce three distinct classes of methods to solve
the tcgre problem: (1) Based on the idea of jsg, we introduce new search
algorithms that do not need to fully construct the jsg in advance and can be
guided by a new admissible heuristic, while guaranteeing optimal solutions;
(2) Inspired by the Conflict-Based Search [13], the second class of algorithms
is based on coordination and we design a Coordination-Exhaustive Search (ces)
algorithm. ces starts with individual optimal paths and finds the lowest cost
among every possible coordination for every robot to achieve the optimal
solution within polynomial time with respect to the number of robots, under
the assumption that the coordination between every pair of support node and
risky edge is only necessary for a limited number of times; (3) Motivated by
ces’s assumption, we also propose a class of receding-horizon sub-team
solutions that further decomposes the order of coordination by only looking at
sub-team coordinations in a local region. We design a Receding-Horizon
Optimistic Cooperative A* (rhoc-a*) search, in order to reduce the time
complexity without much performance loss. Extensive experiment results are
presented and discussed to inform best ways to solve different tcgre problems.

##  II Related Work

We first review related work on the classical mapf problem and common classes
of algorithms. We then review previous approaches to solve the tcgre problem.

###  II-A mapf and Classes of Algorithms

mapf is a specific type of multi-agent planning problem with a key constraint
that no agents can collide with one another [14]. A feasible solution to the
problem is a joint plan that allows all agents to reach their goals from their
starts. Two common objectives are makespan and total cost. Classical mapf
problem may include extra assumptions, such as no vertex conflict, no edge
conflict, no cycle conflict, and no swapping conflict [15, 16].

Algorithms to solve mapf include A*-based search with exponential space and
time compelexity [17, 15], conflict-based search [13] by decomposing into many
constrained single-agent problems, reduction-based approaches to SAT [18, 19],
ILP [20], ASP [21], or CSP [22, 23], rule-based algorithms based on
Kornhauser’s algorithm [24], Push-and-Rotate [25], or BIBOX [26], and
suboptimal solutions [27, 28] to sacrifice optimality for efficiency.

mapf is NP-hard [29], and no optimal solutions can be found in polynomial
time. The time complexity of all above optimal algorithms [17, 15, 18, 19, 13,
20, 21, 22, 23] is exponential to the number of agents. Similarly, we prove in
this paper that our tcgre problem that utilizes, instead of avoiding,
interactions between agents in the form of support is also NP-hard.

###  II-B Team Coordination on Graphs with Risky Edges (tcgre)

tcgre [11] is a recently proposed problem, in which a team of robots traverse
a graph from their starts to goals and also support each other while
traversing certain risky (high-cost) edges to reduce overall cost. Instead of
focusing on collision-free paths in the traditional mapf, the tcgre problem
pursues team coordination. To solve tcgre, Limbu et al. [11] have proposed jsg
and cjsg, both of which construct a single-agent joint-state graph. After the
construction, the original team coordination problem can be solved using
Dijkstra’s algorithm to solve a shortest path problem with optimality
guarantee. The cjsg construction deals with the team coordination problem more
efficiently, although it can only solve problems with two agents. To scale up
tcgre, RL [11] has been utilized to handle many nodes and robots, but at the
cost of optimality.

In this work, we reformulate the tcgre in a constrained optimization framework
and conduct mathematical analysis of this problem. We prove its NP-hardness
and point out the necessity of efficient decomposition to effectively solve
this problem. We further present three classes of algorithms to solve tcgre
from different perspectives.

##  III Problem Formulation

Assuming a team of $N$ homogeneous robots traverse an undirected graph
$\mathbb{G}=(\mathbb{V},\mathbb{E})$, where $\mathbb{V}$ is the set of nodes
the robots can traverse to and $\mathbb{E}$ is the set of edges connecting the
nodes, i.e., $\mathbb{E}\subset\mathbb{V}\times\mathbb{V}$. The team of robots
traverse in the graph from their start nodes $\mathbb{V}_{0}\subset\mathbb{V}$
to goal nodes $\mathbb{V}_{g}\subset\mathbb{V}$ via edges in $\mathbb{E}$.
Each edge $e_{ij}=(V_{i},V_{j})\in\mathbb{E}$ is associated with a cost
$c_{ij}$, depending on its length, condition, traffic, obstacles, etc.
Specially, some edges with high costs are difficult to traverse through,
denoted as risky edges $\mathbb{E}^{\prime}\subset\mathbb{E}$, but with the
support from a teammate from a supporting node, their costs can be
significantly reduced to $\tilde{c}_{ij}$. In this problem, we only consider
such coordination behaviors between two robots. In one coordination behavior,
one receiving robot receives support while traversing a risky edge, and
another supporting robot offers support from some (nearby) location, called
support node. Note that each risky edge $e_{ij}\in\mathbb{E}^{\prime}$
corresponds to certain support node(s) $\mathbb{S}_{e_{ij}}\subset\mathbb{V}$
($\mathbb{S}_{e_{ij}}=\emptyset$ if $e_{ij}\notin\mathbb{E}^{\prime}$).
Additionally, the coordination also induces some cost for the supporter,
denoted by $c^{\prime}$. A central planner needs to schedule the paths of all
agents and coordination on their ways.

###  III-A Action & Cost Model

Without coordination, at each time step $t$, a robot $n$ can choose to stay
where it is, or move to its neighbor ($V_{i}$ is the neighbor of $V_{j}$ if
$e_{ij}\in\mathbb{E}$). Its movement can be denoted by
$M_{n}^{t}=(l_{n}^{t},l_{n}^{t+1})\in\mathbb{E}$, where
$l_{n}^{t},l_{n}^{t+1}\in\mathbb{V}$ indicate its current and next location
and $l_{n}^{t+1}$ is a neighbor of $l_{n}^{t}$. Specially, the robot stays at
its current location if $l_{n}^{t+1}=l_{n}^{t}$ with zero cost, i.e.,
$c_{ii}=0,\forall i$. The movement set can thus be denoted by
$\mathcal{M}=\\{M_{n}^{t}|\forall n,\forall t\\}$. Moreover, the movement
decision $M_{n}^{t}$ can be rewritten as an 0/1 variable $M_{ij}^{nt}$, where
$M_{ij}^{nt}=1$ represents edge $e_{ij}$ is selected by robot $n$ at time $t$,
and 0 otherwise. A robot can only move once at each time step to a neighbor
node or not at all, i.e., $\sum_{\forall
e_{ij}\in\mathcal{N}_{l_{n}^{t}}}M_{ij}^{nt}=1~{}\&~{}\sum_{\forall
e_{ij}\notin\mathcal{N}_{l_{n}^{t}}}M_{ij}^{nt}=0$, where
$\mathcal{N}_{l_{n}^{t}}=\\{(l_{n}^{t},l_{n}^{t+1})|\forall(l_{n}^{t},l_{n}^{t+1})\in\mathbb{E}\\}$.
The movement set can thus be denoted by $\mathcal{M}=\\{M_{ij}^{nt}|\forall
i,j,\forall n,\forall t\\}$.

When a coordination behavior is available—when robot $n$ is going to traverse
a risky edge, another robot $m$ happens to be at one of the support nodes of
the risky edge or vice versa, i.e., $M_{n}^{t}\in\mathbb{E}^{\prime}$ and
$l_{m}^{t}\in\mathbb{S}_{M_{n}^{t}}$, or $M_{m}^{t}\in\mathbb{E}^{\prime}$ and
$l_{n}^{t}\in\mathbb{S}_{M_{m}^{t}}$—the robot pair needs to decide whether to
provide/receive support. Denote the coordination decision of agent $n$ at time
$t$ as $s_{nm}^{t}$. It is clear that agent $n$’s coordination decision is
dependent on its movement decision, so the cost is twofold: (1) When agent $n$
has no coordination opportunity (the above coordination behavior is not
available for any other robot $m$), i.e., $\forall m$,
$l_{m}^{t}\notin\mathbb{S}_{M_{n}^{t}}$ and
$l_{n}^{t}\notin\mathbb{S}_{M_{m}^{t}}$, its cost $C_{n}^{t}$ is only decided
by its movement, i.e., $C_{n}^{t}=c_{ij}$, where $M_{ij}^{nt}=1$.

(2) When coordination is possible for agent $n$, i.e., $\exists m$,
$l_{m}^{t}\in\mathbb{S}_{M_{n}^{t}}$ or $l_{n}^{t}\in\mathbb{S}_{M_{m}^{t}}$,
the cost $C_{n}^{t}$ can be represented as

$\displaystyle
C_{n}^{t}=\begin{cases}c_{ij},~{}~{}\textit{if}~{}~{}s_{nm}^{t}=~{}~{}0;\\\
\tilde{c}_{ij}~{}~{}~{}\textit{if}~{}~{}s_{nm}^{t}=~{}~{}1;\\\
c^{\prime},~{}~{}~{}\textit{if}~{}~{}s_{nm}^{t}=-1.\end{cases}$ (1)

where $s_{nm}^{t}=1$ means agent $n$ decides to receive support from $m$,
$s_{nm}^{t}=-1$ indicates agent $n$ decides to offer support to $m$, and
$s_{nm}^{t}=0$ stands for no coordination between the robot pair $n$ and $m$
at $t$. Specially, no coordination happens for one single robot, i.e.,
$s_{nn}^{t}=0,\forall n,\forall t$, or when $n$ and $m$ cannot support each
other, i.e., $s_{nm}^{t}=0$ if $\forall m$,
$l_{m}^{t}\notin\mathbb{S}_{M_{n}^{t}}$ and
$l_{n}^{t}\notin\mathbb{S}_{M_{m}^{t}}$. The coordination decision set can be
written as $\mathcal{S}=\\{s_{nm}^{t}|\forall n,m,\forall t\\}$. In addition,
a coordination decision is made for a pair, so $s_{nm}^{t}+s_{mn}^{t}=0$ for
every pair of robots. Furthermore, the robots can wait now (no movement) for
future coordination, but there is no point for all robots to stay still at the
same time, i.e., $\sum_{\forall n}\sum_{\forall i\neq j}M_{ij}^{nt}\neq 0$.

###  III-B Problem Definition

Given the node set $\mathbb{V}$, the edge set $\mathbb{E}$, support nodes for
each edge $\mathbb{S}_{e_{ij}}$, cost of each edge without and with
coordination $c_{ij},\tilde{c}_{ij}$, $N$ robots with their starts
$\mathbb{V}_{0}$ and goals $\mathbb{V}_{g}$, optimize the movement and
coordination decisions $\mathcal{M}$ and $\mathcal{S}$, in order to minimize
the total cost for each agent to traverse from its start to its goal within a
time limit $T$. Formally, the problem can be represented as

$\displaystyle\min_{{\mathcal{M}},{\mathcal{S}}}\sum_{t=0}^{T-1}\sum_{n=1}^{N}C_{n}^{t}.$
(2) $\displaystyle s.t.~{}~{}\sum_{\forall
m\in\\{1,2,...,N\\}}|s_{nm}^{t}|\leq 1,$ $\displaystyle\forall
n\in\\{1,2,...,N\\},\forall t\in\\{0,1,...,T-1\\}.$ (3) $\displaystyle
s_{nm}^{t},s_{mn}^{t}\in\\{-1,0,1\\},s_{nm}^{t}+s_{mn}^{t}=0,$
$\displaystyle\forall n,m\in\\{1,2,...,N\\},\forall t\in\\{0,1,...,T-1\\}.$
(4) $\displaystyle l_{n}^{0}=\mathbb{V}_{0}(n),\forall n\in\\{1,2,...,N\\}.$
(5) $\displaystyle l_{n}^{T}=\mathbb{V}_{g}(n),\forall n\in\\{1,2,...,N\\}.$
(6) $\displaystyle\sum_{\forall
e_{ij}\in{\mathcal{N}_{l_{n}^{t}}}}M_{ij}^{nt}=1,\sum_{\forall
e_{ij}\notin{\mathcal{N}_{l_{n}^{t}}}}M_{ij}^{nt}=0,$ $\displaystyle\forall
n\in\\{1,2,...,N\\},\forall t\in\\{0,1,...,T-1\\}.$ (7)
$\displaystyle\sum_{\forall n}\sum_{\forall i\neq j}M_{ij}^{nt}\neq 0,$
$\displaystyle\forall t\in\\{0,1,...,T-1\\}.$ (8)

Eqn. (2) suggests the goal of the problem is to minimize the total cost of all
agents across all time steps with two decision variables, movement set
$\mathcal{M}$ and coordination set $\mathcal{S}$. Eqn. (3) indicates, at each
time step, each robot can participate in at most one coordination behavior.
Eqn. (4) regulates that, at each time step, one coordination behavior only
occurs between one robot pair. Eqn. (5) and Eqn. (6) set the start and the
goal for each robot. Eqn. (7) guarantees, at each time step, a robot can only
move to a neighbor node or stay still. Eqn. (8) assures no unnecessary
stagnation.

##  IV Mathematical Analysis

Figure 1: Reduction from Maximum 3D Matching (Middle) to tcgre (Left and
Right) and Inspiration for ces.

In this section, we prove our tcgre problem reduces from the Maximum 3D
Matching problem, a NP-hard problem. Then, we start a rigorous analysis on the
mathematical problem, which suggests decomposition is a promising solution to
this combinatorial optimization problem.

###  IV-A NP-Hardness

######  Definition 1.

Maximum 3D Matching: ${\bf{X}},{\bf{Y}},{\bf{Z}}$ are 3 finite sets.
${\bf{T}}$ is the subset of ${\bf{X}}*{\bf{Y}}*{\bf{Z}}$, with triples
$(x,y,z)$, where $x\in{\bf{X}},y\in{\bf{Y}},z\in{\bf{Z}}$.
${\bf{M}}\subset{\bf{T}}$ is a 3D matching if for any two distinct triples
$(x_{1},y_{1},z_{1})$ and $(x_{2},y_{2},z_{2})\in{\bf{M}}$, we have $x_{1}\neq
x_{2}$, $y_{1}\neq y_{2}$, and $z_{1}\neq z_{2}$; each triple has a weight
$w(x_{i},y_{j},z_{k})$. Maximum 3D Matching problem is to find a 3D matching
with maximum total weight.

######  Theorem 1.

Maximum 3D Matching reduces to tcgre.

Proof sketch: Without coordination, we can generate individual optimal paths
for all robots; the cost is $C_{0}$ (Fig. 1 left). With coordination, we let
some robot pairs take detours to the support pairs (i.e., one risky edge and
one support node); the new cost is $C^{\prime}$, where the first robot in the
pair traverses the risky edge while the second visits and supports from the
support node (Fig. 1 right). Minimizing total cost is equivalent to maximizing
the cost reduction $C_{0}-C^{\prime}$.

Consider X contains all robot pairs; ${\bf{Y}}$ is the set of all support
pairs plus an empty set element; Z is a list of time orders of events (time
steps are not necessarily needed, because robots can stay still and wait, Fig.
1 middle). Consider the weight $w(x_{i},y_{j},z_{k})$ as the sum of cost
reduction through coordination for robot pair $x_{i}$ to detour to support
pair $y_{j}$ with time order $z_{k}$. The weight of a 3D matching M is the
total cost reduction of all robots, $C_{0}-C^{\prime}$. Maximum 3D Matching is
to find the paths for all robot pairs to achieve maximum total cost reduction
from their original without-coordination costs, which is a specific case of
the general tcgre problem that only needs to use each robot and support pair
once with a smaller solution space than tcgre.

Since tcgre reduces from Maximum 3D Matching, a classical NP-hard problem
[30], tcgre is also NP-hard. We cannot find an optimal solution in polynomial
time. ∎

###  IV-B Problem Analysis

Because we only care about the total cost of all robots, and every
coordination requires a robot pair, we can reassign the coordination cost
$c^{\prime}$ from the supporter to the receiver in addition to the reduced
edge cost $\tilde{c}_{ij}$, without changing the problem. So, the cost with
coordination in Eqn. (1) becomes

$\displaystyle
C_{n}^{t}=\begin{cases}c_{ij},~{}~{}\textit{if}~{}~{}s_{nm}^{t}=~{}~{}0;\\\
\hat{c}_{ij},~{}~{}\textit{if}~{}~{}s_{nm}^{t}=~{}~{}1;\\\
0,~{}~{}~{}~{}\textit{if}~{}~{}s_{nm}^{t}=-1.\end{cases}$ (9)

where $\hat{c}_{ij}=\tilde{c}_{ij}+c^{\prime}$. Therefore, the original
objective function (Eqn. (2)) can be rewritten as

$\displaystyle\min_{{\mathcal{M}},{\mathcal{S}}}\sum_{t=0}^{T-1}\sum_{n=1}^{N}\sum_{\forall
i\neq j}M_{ij}^{nt}[(1-s_{nm}^{t})c_{ij}+s_{nm}^{t}\hat{c}_{ij}],$ (10)
$\displaystyle
s.t.~{}~{}\eqref{constraint1:supportnum},\eqref{constraint2:supportpair},\eqref{constraint3:starts},\eqref{constraint4:goals},\eqref{constraint5:movenum},\eqref{constraint6:correlation}.$

Notice that because when $s_{nm}^{t}=-1$, $M_{ij}^{nt}=0,\forall i\neq j$, the
last if condition in Eqn. (9) does not need to be considered in Eqn. (10). To
solve the combinatorial optimization problem, a typical approach is dynamic
programming [31], by decoupling the interdependency among the decision
variables. The ideal case is to break down the problem into two sub-problems:
one optimizing the movement decisions ${\mathcal{M}}$ and the other optimizing
the coordination decisions ${\mathcal{S}}$.

Based on such a motivation, if we can find a way to eliminate $n$ from
$M_{ij}^{nt}$, Eqn. (10) can be rewritten as

$\displaystyle\underbrace{\min_{\mathcal{M}}\sum_{t=0}^{T-1}\sum_{\forall
i\neq j}M_{ij}^{nt}}_{\textrm{Sub-Problem
2}}\underbrace{\min_{\mathcal{S}}\sum_{n=1}^{N}[(1-s_{nm}^{t})c_{ij}+s_{nm}^{t}\hat{c}_{ij}]}_{\textrm{Sub-
Problem 1}}.$ (11)

The first half of the function contains only the movement decisions
$\mathcal{M}$, while the second half only has the coordination decisions
$\mathcal{S}$. Then, decomposition is possible and the NP-hard tcgre problem
can be solved with significantly reduced complexity.

##  V Solutions

Based on the mathematical analysis, we propose three classes of algorithms to
solve tcgre from different perspectives with different optimality and
efficiency characteristics. We first propose a class of jsg-based solutions
that utilizes the decomposition of the original problem and provides optimal
solutions to the two sub-problems (Eqn. 11). Second, to reduce the time
complexity, we focus on coordination (i.e., $\mathcal{S}$) and propose a class
of coordination-based solutions that decomposes the problem differently. To be
specific, we introduce Coordination-Exhaustive Search (ces), which can achieve
optimal solutions under a reasonable assumption that each coordination
behavior (i.e., support pair composed of support node and risky edge) is only
needed for a constant time in the optimal solution. Finally, when the same
support may need to be repeated many times and the global optimality does not
need to be guaranteed, a class of algorithms that focuses only on local sub-
team coordination behaviors are introduced, for which we develop Receding-
Horizon Optimistic Cooperative A* (rhoc-a*).

###  V-A jsg-Based Solutions

jsg-based solutions perform the decomposition into the two sub-problems by 1)
implicitly solving $\mathcal{S}$ by calculating the minimum edge cost for each
edge in the jsg; and 2) explicitly solving $\mathcal{M}$ by solving a single-
robot shortest path problem with $\mathcal{S}$ implicitly encoded. As
mentioned in Sec. IV-B, $n$ is effectively eliminated from $M_{ij}^{nt}$ by
building the jsg.

####  V-A1 jsg Construction

In the action model (Sec. III-A), we use $l_{n}^{t}\in\mathbb{V}$ to represent
robot $n$’s location at time $t$. In a joint state graph, however, one state
is the set of all robots’ locations
${\mathbf{L}}^{t}=\\{l_{1}^{t},l_{2}^{t},...,l_{N}^{t}\\}$. The new node set
is $\mathbb{L}=\mathbb{V}^{N}$ and each node ${\mathbf{L}_{p}}\in\mathbb{L}$
correlates to $N$ nodes $(V_{i_{1}},V_{i_{2}},...,V_{i_{n}})$ in $\mathbb{V}$.
By checking the constraints in Eqns. (3), (4), (7), and (8), for each pair of
joint-states, we can form the new edge set $\mathbb{M}\subset\mathbb{L}^{2}$.
Specially, Eqn. (8) assures no self-loops in the jsg. An action is the move
from current state to next state
${\mathbf{M}}^{t}=({\mathbf{L}}^{t},{\mathbf{L}}^{t+1})$, which can also be
written as a 0/1 variable, i.e.,
${\mathbf{M}}_{pq}^{t}=\prod_{n=1}^{N}M_{ij}^{nt}$, where ${\mathbf{e}}_{pq}$
is any edge in the new edge set, and the movement decision set becomes
$\mathcal{M}=\\{{\mathbf{M}}^{t}|\forall t\\}$.

####  V-A2 Sub-Problem 1

After the construction of jsg, the second half of Eqn. (11) is calculating the
minimum edge cost for each $\mathbf{e}_{pq}\in\mathbb{M}$ by optimizing
$\mathcal{S}$, which is a 0/1 Integer Linear Programming (ILP) problem and can
be solved by classical methods, such as Branch and Bound [32].

####  V-A3 Sub-Problem 2

After solving sub-problem 1, we have the cost $C_{pq}$ for each edge
${\mathbf{e}}_{pq}\in\mathbb{M}$. Sup-problem 2 is to optimize movement
decisions $\mathcal{M}$ to minimize the total cost:

$\displaystyle\min_{{\mathcal{M}}}~{}~{}\sum_{t=0}^{T-1}\sum_{\forall{\mathbf{e}}_{pq}\in\mathbb{M}}$
$\displaystyle{\mathbf{M}}_{pq}^{t}C_{pq}^{t},$ (12)

Now it is a single-robot shortest path problem with non-negative costs and no
self-loops solvable by any shortest path algorithms. Instead of using
Dijkstra’s algorithm on a fully constructed jsg beforehand [11], we present
results with Uniform Cost Search (ucs) and A* (guided by an admissible
heuristics assuming all future risky edges will be supported by a teammate)
while constructing jsg on the fly, i.e., interleaving partial solutions of the
two sub-problems.

###  V-B Coordination-Based Solutions

Some jsg edge costs are simply the sum of individual edge costs of all robots,
suggesting possible total cost separation into the costs with and without
coordination. Thus, while minimizing costs without coordination can be simply
solved for individual robots, the second class of algorithms focuses on
coordination. Specifically, we present a Coordination-Exhaustive Search (ces)
algorithm based on a slightly different and interleaving decomposition of Eqn.
(10):

$\displaystyle\min_{{\mathcal{M}},{\mathcal{S}}}\overbrace{\sum_{t=0}^{T-1}\sum_{n=1}^{N}\sum_{\forall
e_{ij}\in\mathbb{E}}M_{ij}^{nt}c_{ij}}^{\textrm{Cost without Coordination}}-$
(13) $\displaystyle\underbrace{\sum_{t=0}^{T-1}\sum_{n=1}^{N}\sum_{\forall
e_{ij}\in\mathbb{E}^{\prime}or}^{\begin{subarray}{c}l_{n}^{t}=l_{n}^{t+1}\in\\\
\bigcup_{m=1}^{N}\mathbb{S}_{M_{m}^{t}}\end{subarray}}\frac{1}{2}M_{ij}^{nt}[(1-s_{nm}^{t})c_{ij}+(1+s_{nm}^{t})\Delta{c}_{ij}]}_{\textrm{Cost
Reduction due to Coordination}}.$

where $\Delta{c}_{ij}={c}_{ij}-\hat{c}_{ij}$. When $s_{nm}^{t}=1$ (receiving
support) the second part reduces the cost in the first part by
$\Delta{c}_{ij}$; When $s_{nm}^{t}=-1$ (providing support), the cost is
reduced to zero. Inspired by Conflict-Based Search, we can start with finding
the individual shortest path for each robot without coordination. Then, we
find the coordination behaviors (with some detours) that can cause the maximum
cost reduction.

While the first half can be solved individually for each robot (Fig. 1 left),
ces uses an exhaustive search for the second half. If an optimal solution
requires a coordination behavior between a robot pair, it is equivalent to
make the robot pair detour to the support pair (risky edge and support node)
from their original individual shortest paths, while other robots remain on
their individual shortest paths (Fig. 1 right). This robot pair’s shortest
paths are a combination of two path segments—their shortest paths from their
starts to the support pair, and their shortest paths from the support pair to
their goals. Given certain coordination behaviors, we can solve the shortest
path and minimum total cost of all robots for each path segment (Fig. 1
right), as shown in Alg. 1.

ces is a coordination-based method through an exhaustive search. Because one
support pair may be assigned to the same/different robot pair(s) for infinite
times, to conduct an exhaustive search, we assume each support pair can only
occur for a fixed number of times. In our implementation, we assume a support
pair can only happen once, but it can be easily expanded to a more general
case, by repetitively adding the same support pairs to the coordination set.

The ces algorithm is shown in Alg. 2: In lines 1-2, using any shortest path
algorithm, it generates an individual optimal path and cost for each robot,
with original edge costs $\mathbf{c}=\\{c_{ij}|\forall e_{ij}\in\mathbb{E}\\}$
and reduced edge costs $\hat{\mathbf{c}}=\\{\hat{c}_{ij}|\forall
e_{ij}\in\mathbb{E}\\}$, called pessimistic/optimistic paths
$\mathbb{P}_{1}/\mathbb{P}_{2}$ and total cost $C_{1}/C_{2}$. If
$C_{1}=C_{2}$, which means no coordination is needed in the optimal solution,
then simply return $\mathbb{P}_{1}$ and $C_{1}$ (lines 3-4). Else, it starts
the scheduling process. In line 5, it generates a coordination set,
$\mathbb{CS}$ that contains all coordination behaviors, and $\mathbb{SCS}$
that contains all subsets of $\mathbb{CS}$, to decide which coordination
behaviors/support pairs are needed. In line 6, it generates a set of all robot
pairs $\mathbb{RP}$ (order in the pair matters since we need to decide which
robot moves to the support node/risky edge) to determine which robot pair
should be assigned to each support pair. Now, it looks like a Maximum 3D
Matching problem as stated in Definition 1, except that the matching problem
has one more assumption that one robot pair can be only assigned once as shown
in Fig. 1 middle. For all possible sets of support pairs (line 8), it
generates all possible time orders for this support pair set using permutation
(line 9). Notice that robots can wait for one another, so the order of each
coordination behavior, not necessarily the exact time step, is sufficient.
Then, it iterates through every possible support pair order (line 10). There
could be $N(N-1)$ possible robot pairs assigned to each support pair, so a
N-Fold Cartesian Product is applied to generate all possible sets of support
robot pairs (line 11). Then, it explores every set (line 12), where each robot
pair in ${\mathbf{SRP}}$ is assigned to each support pair in ${\mathbf{PSCS}}$
with the same index. Thus, in lines 13-18, it adds the risky edge of each
support pair to the individual coordination set of the first robot of the
robot pair, and the support node of each support pair to the individual
coordination set of the second robot of the robot pair. ${\mathbf{ICS}}$ then
contains individual coordination set of all robots. With ${\mathbf{ICS}}$, we
use Alg. 1 to calculate the shortest paths and minimum total cost of this
assignment in line 19, which is one solution. Last, lines 20-21 records the
best solution with minimum total cost. To sum up, the loop in line 8 decides a
subset of support pairs we need for cost reduction. The loop in line 10
decides an order for the subset. The loop in line 12 selects a robot pair for
each support pair in the subset. Iterating through the three loops explores
every possible solution under the assumption that each support pair can be
applied for a constant number of times, making ces optimal.

1 ${\mathbb{P}}=[\emptyset]*N$;

2 $totalcost=0$;

3 for  _$n=1$ to $N$_ do

4 if  _${\mathbf{ICS}}[n]=\emptyset$_ then

5
${\mathbf{P}},C=ShortestPath(\mathbb{G},\mathbb{V}_{0}[n],\mathbb{V}_{g}[n],\mathbf{c})$;

6 ${\mathbb{P}}[n]={\mathbf{P}}$;

7 $totalcost=totalcost+C$;

8 else

9 for  _$item\in{\mathbf{ICS}}[n]$_ do

10 $start=\mathbb{V}_{0}[n]$;

11 if  _$item$ is a risky edge_ then

12 ${\mathbf{P}},C=ShortestPath(\mathbb{G},start,item[0],\mathbf{c})$;

13 ${\mathbb{P}}[n].extend({\mathbf{P}}\bigcup item)$;

14 $totalcost=totalcost+C+\hat{c}_{item}$;

15 $start=item[1]$;

16 else if  _$item$ is a support node_ then

17 ${\mathbf{P}},C=ShortestPath(\mathbb{G},start,item,\mathbf{c})$;

18 ${\mathbb{P}}[n].extend({\mathbf{P}})$;

19 $totalcost=totalcost+C$;

20 $start=item$;

21${\mathbf{P}},C=ShortestPath(\mathbb{G},start,\mathbb{V}_{g}[n],\mathbf{c})$;

22 ${\mathbb{P}}[n].extend({\mathbf{P}})$;

23 $totalcost=totalcost+C$;

24return _$\mathbb{P},totalcost$_

Algorithm 1 $CostCalculation$
($\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g},{\mathbf{ICS}}$)

1
$\mathbb{P}_{1},C_{1}=MultipleShortestPath(\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g},\mathbf{c})$;

2
$\mathbb{P}_{2},C_{2}=MultipleShortestPath(\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g},\hat{\mathbf{c}})$;

3 if  _$C_{1}\leq C_{2}$_ then

4 return _$\mathbb{P}_{1},C_{1}$_ ;

5Generate subsets
$\mathbb{SCS}=\\{{\mathbf{SCS}}|\forall{\mathbf{SCS}}\subset\mathbb{CS}\\}$ of
the coordination set $\mathbb{CS}$;

6 Generate a set of all robot pairs $\mathbb{RP}$;

7 $\mathbb{P}_{min},C_{min}\leftarrow\mathbb{P}_{1},C_{1}$;

8 for  _${\mathbf{SCS}}\in\mathbb{SCS}$_ do

9 ${\mathbb{PSCS}}=AllPermutations({\mathbf{SCS}},len({\mathbf{SCS}}))$;

10 for  _${\mathbf{PSCS}}\in\mathbb{PSCS}$_ do

11 ${\mathbb{SRP}}=CartesianProduct(\mathbb{RP},len({\mathbf{PSCS}}))$;

12 for  _${\mathbf{SRP}}\in{\mathbb{SRP}}$_ do

13 ${\mathbf{ICS}}=[\emptyset]*N$;

14 for  _$n=1$ to $len(\bf{PSCS})$_ do

15$SP={\bf{PSCS}}[n]$;

16 $RP={\bf{SRP}}[n]$;

17 ${\mathbf{ICS}}[RP[0]].append(SP[0])$;

18 ${\mathbf{ICS}}[RP[1]].append(SP[1])$;

19$\mathbb{P},C=CostCalculation(\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g},{\mathbf{ICS}})$;

20 if  _$C <C_{min}$_ then

21 $\mathbb{P}_{min},C_{min}\leftarrow\mathbb{P},C$;

22return _$\mathbb{P}_{min},C_{min}$_.

Algorithm 2 ces ($\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g}$)

There are total $O(2^{|\mathbb{CS}|})$ subsets in $\mathbb{SCS}$ . For each
subset, there are $O(|\mathbb{CS}|!)$ permutations. For every permutation,
there are $O(N^{2})$ possible robot pairs for each support pair, totally
$O(N^{2|\mathbb{CS}|})$ assignments. Therefore, the number of possible
solutions is $O(2^{|\mathbb{CS}|}\cdot
N^{2|\mathbb{CS}|}\cdot|\mathbb{CS}|!)=O((2N^{2})^{|\mathbb{CS}|}\cdot|\mathbb{CS}|!)$.
For each solution, we run the shortest path algorithm for $O(|\mathbb{CS}|)$
times, each run costs $O(|\mathbb{E}|log(|\mathbb{V}|))$. Therefore, the time
complexity of ces is
$O((2N^{2})^{|\mathbb{CS}|}\cdot|\mathbb{CS}|!\cdot|\mathbb{E}|log(|\mathbb{V}|))$,
which is not exponential to the number of robots $N$ anymore. Note that the
above algorithm is for directed graphs. For an undirected graph, each edge is
actually two directed edges in a directed graph, so after line 10, there
should be an additional loop that iterates through all possible directions of
selected risky edges, which won’t change the time complexity class.

###  V-C Receding-Horizon Sub-Team Solutions

tcgre’s computation complexity arises from two fronts, the large size of the
graph and the team. Therefore, the third class of algorithms reduces the
complexity from both fronts by planning with a limited horizon and for a sub-
team of all robots at a time, efficiently facilitating local coordination
while compromising global optimality. Such sub-team local coordination within
the receding horizon prioritizes actions that yield the best short-term
outcomes and potentially allows dynamic adaptation to changing circumstances,
e.g., updated graph structure from robot perception. One specific algorithm is
Receding-Horizon Optimistic Cooperative A* (rhoc-a*), which provides flexible
and efficient solution by assuring optimal robot pair coordination within the
horizon while assuming optimistic cooperation beyond the horizon (Fig. 2).

Figure 2: Receding-Horizon Optimistic Cooperative A*.

Alg. 3 presents rhoc-a*. All robots are initially not at their goals and
therefore on duty (lines 1-2). We compute the heuristic for all nodes assuming
always-available support (line 3). rhoc-a* iterates until all robots arrive at
their goals (lines 4). If there is still at least one pair of robots on duty
(line 5), rhoc-a* sequentially plans for each robot pair (lines 6-7). The
computation efficiency is enabled by only looking at a small jsg with only two
robots, $n$ and $m$, within $K$ steps (lines 8-12). Notice that optimality is
assured on this small jsg while the cost-to-go on the horizon $K$ is estimated
by the always-available support heuristic. Lines 13-18 address the situation
where only one robot is on duty, who has to traverse to its goal alone.

For $N$ robots, generating all possible pairs of robots can be done in
$O(N^{2})$ time. rhoc-a*’s time complexity for each pair’s search of $K$ steps
can be approximated as $O(b^{K})$, where
$b=O(\frac{|\mathbb{E}|}{|\mathbb{V}|})$ is the effective branching factor in
the joint action space. Thus, in $K$ steps, running rhoc-a* for all pairs
results in complexity of $O(N^{2}(\frac{|\mathbb{E}|}{|\mathbb{V}|})^{K})$.
There will be no cycles for any robot, due to graph search. As a result, there
will be at most $O(|V|)$ steps for each robot, so we need to run the $K$-step
A* search for $\frac{|V|}{K}$ times for one robot pair. There will be
$\frac{N}{2}$ runs for $N$ robots. Therefore, the time complexity of the
algorithm is
$O(N^{3}\cdot\frac{|V|}{K}\cdot(\frac{|\mathbb{E}|}{|\mathbb{V}|})^{K})$.

1 Initialize $atGoal_{n}\leftarrow False,\forall n\in\\{1,2,...,N\\}$

2 ${\bf{OnDuty}}=\\{n|\forall n,atGoal_{n}=False\\}$

3 Compute optimistic heuristic $h(\cdot)$ for each node in $\mathbb{G}$

4 while  _$len({\bf{OnDuty}})\neq 0$_ do

5 if  _$len({\bf{OnDuty}})\geq 2$_ then

6 $\mathbb{RP}=\\{\\{n,m\\}|\forall n,m\in\bf{OnDuty}\\}$

7 for  _each pair $\\{n,m\\}\in\mathbb{RP}$_ do

8 Initialize a $K$-step jsg $\mathbb{G}_{nm}$ for the pair $\\{n,m\\}$ with
start $\\{l_{n}^{t},l_{m}^{t}\\}$;

9 if  _not $atGoal_{n,m}$_ then

10 A* on $\mathbb{G}_{nm}$ for $K$ steps using $h(\cdot)$;

11 Update $l_{n}^{t},l_{m}^{t}$, $atGoal_{n,m}$, and ${\bf{OnDuty}}$;

12 Update individual and total costs;

13else

14 $n={\bf{OnDuty}}.pop()$;

15 Initialize $K$-step graph $\mathbb{G}_{n}$ with start $l_{n}^{t}$;

16 A* on $\mathbb{G}_{n}$ for $K$ steps using $h(\cdot)$;

17 Update $l_{n}^{t}$, $atGoal_{n}$, and ${\bf{OnDuty}}$;

18 Update individual and total costs;

19return Paths and costs for all robots.

Algorithm 3 rhoc-a*($\mathbb{G},\mathbb{V}_{0},\mathbb{V}_{g}$, $K$)

##  VI Results

Figure 3: True and Naive Optimality vs. Time with jsg-ucs, jsg-a*, ces,
rhoc-a*, and Naive. Each data point denotes the result for the experiment with
# of Robots–# of Nodes. For visibility, cluttered areas are magnified in the
dashed boxes.

We conduct experiments on a variety of graphs to evaluate the optimality and
efficiency of the three classes of proposed algorithms. To be specific, we
implement ucs and a* for jsg-based solutions, ces for coordination-based
solutions, and rhoc-a* for receding-horizon sub-team solutions. We first
conduct a set of large-scale, method-agnostic experiments on a variety of
randomly generated graphs and then present focused experiments to study the
pros and cons of specific methods.

###  VI-A Large-Scale Method-Agnostic Experiments

To evaluate each method in an objective manner, we generate a set of graphs
with randomly generated support pairs, including sparse, moderate, and dense
connectivities and five different numbers of nodes
($|\mathbb{V}|\in\\{10,15,20,25,30\\}$), three graphs each type, i.e., a total
of 45 distinct graphs. A total of 900 trials are conducted with five different
team sizes ($N\in\\{3,4,5,6,7\\}$) and four methods.

We evaluate the optimality and runtime of all methods along with a naive
approach, in which each robot executes its individual optimal path without
coordination. While the True Optimality value is defined as the optimal cost
divided by the actual cost, for scenarios where the optimal cost cannot be
found due to excessive computation, we define Naive Optimality to be the naive
cost divided by actual cost. If a data point does not exist in Fig. 3, the
corresponding method cannot produce a solution for the robot and node number.
As shown in Fig. 3, the jsg-based solutions achieve optimal solutions, but
require significant runtime even in small graphs with only a few robots and
fail to produce a solution when the problem becomes larger; ces has better
runtime but loses some performance because we assume each support pair can be
applied only once; with a fine-tuned $K$, surprisingly, rhoc-a* in many cases
achieves better results than ces with less runtime.

Figure 4: ces Planning Time on Graphs with Two Support Pairs

###  VI-B Focused Experiments

####  VI-B1 ces’s Insensitivity to Robot and Node Numbers

Fig. 4 showcases that, when there are not many support pairs, ces works well
with different numbers of robots and different sizes of graphs (polynomial
time to both $N$ and $|\mathbb{V}|$). However, its runtime increases
drastically with the number of support pairs, as shown in our method-agnostic
experiments, which verifies our time complexity analysis (Sec. V-B).

####  VI-B2 rhoc-a*’s Sensitivity to Planning Horizon

Fig. 5 showcases how rhoc-a*’s computation time scales with different planning
horizons. A large horizon $K$ comes closer to solving the original tcgre
problem with multiple robot pairs, which significantly increases the solution
time. While the total cost can be reduced with a longer horizon, it is
necessary to strike a balance between horizon and efficiency.

Figure 5: rhoc-a* Planning Time.

##  VII Conclusions and Discussions

We present a systematic problem formulation and mathematical analysis of
tcgre, which proves its NP-hardness and shows efficient decomposition is the
key to solving this problem. We propose three classes of solutions with a set
of implementations and present their experiment results.

As given by the analysis in Sec. IV-B, all of the proposed solutions are
trying to solve a form of a decomposed problem. For example, jsg-based
solutions solve a 0/1 ILP problem [20] and a single-agent shortest path
problem, after constructing a jsg; coordination-based solutions, like ces,
deal with a 3D matching problem embedded with multiple single-agent shortest
path problems. By applying some approximation methods to the subproblems—for
the former, forming only a few edges instead of all feasible edges and
calculating approximate edge costs; for the latter, omitting unpromising
matchings—we can significantly reduce runtime without sacrificing too much
performance. rhoc-a*, though efficient, does not consider the order of the
robot pair selection, with which its performance can improve while still
maintaining coordination efficiency.

## References

  * [1] J. L. Adler and V. J. Blue, “A cooperative multi-agent transportation management and route guidance system,” Transportation Research Part C: Emerging Technologies, vol. 10, no. 5-6, pp. 433–454, 2002.
  * [2] Y. Hu, M. Chen, W. Saad, H. V. Poor, and S. Cui, “Distributed multi-agent meta learning for trajectory design in wireless drone networks,” IEEE Journal on Selected Areas in Communications, vol. 39, no. 10, pp. 3177–3192, 2021.
  * [3] C. Yu, A. Velu, E. Vinitsky, J. Gao, Y. Wang, A. Bayen, and Y. Wu, “The surprising effectiveness of ppo in cooperative, multi-agent games,” 2022.
  * [4] P. Surynek, “An optimization variant of multi-robot path planning is intractable,” in Proceedings of the AAAI conference on artificial intelligence, vol. 24-1, pp. 1261–1263, 2010.
  * [5] J. Foerster, N. Nardelli, G. Farquhar, T. Afouras, P. H. Torr, P. Kohli, and S. Whiteson, “Stabilising experience replay for deep multi-agent reinforcement learning,” in International conference on machine learning, pp. 1146–1155, PMLR, 2017.
  * [6] J. K. Gupta, M. Egorov, and M. Kochenderfer, “Cooperative multi-agent control using deep reinforcement learning,” in Autonomous Agents and Multiagent Systems: AAMAS 2017 Workshops, Best Papers, São Paulo, Brazil, May 8-12, 2017, Revised Selected Papers 16, pp. 66–83, Springer, 2017.
  * [7] M. Liu, H. Ma, J. Li, and S. Koenig, “Task and path planning for multi-agent pickup and delivery,” in Proceedings of the International Joint Conference on Autonomous Agents and Multiagent Systems (AAMAS), 2019.
  * [8] Z. Yan, N. Jouandeau, and A. A. Cherif, “A survey and analysis of multi-robot coordination,” International Journal of Advanced Robotic Systems, vol. 10, no. 12, p. 399, 2013.
  * [9] R. Luna and K. E. Bekris, “Efficient and complete centralized multi-robot path planning,” in 2011 IEEE/RSJ International Conference on Intelligent Robots and Systems, pp. 3268–3275, IEEE, 2011.
  * [10] J. Capitan, M. T. Spaan, L. Merino, and A. Ollero, “Decentralized multi-robot cooperation with auctioned pomdps,” The International Journal of Robotics Research, vol. 32, no. 6, pp. 650–671, 2013.
  * [11] M. Limbu, Z. Hu, S. Oughourli, X. Wang, X. Xiao, and D. Shishika, “Team coordination on graphs with state-dependent edge costs,” in 2023 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), pp. 679–684, IEEE, 2023.
  * [12] M. Limbu, Z. Hu, X. Wang, D. Shishika, and X. Xiao, “Team coordination on graphs with reinforcement learning,” in 2024 IEEE International Conference on Robotics and Automation (ICRA), IEEE, 2024.
  * [13] G. Sharon, R. Stern, A. Felner, and N. R. Sturtevant, “Conflict-based search for optimal multi-agent pathfinding,” Artificial Intelligence, vol. 219, pp. 40–66, 2015.
  * [14] R. Stern, “Multi-agent path finding–an overview,” Artificial Intelligence: 5th RAAI Summer School, Dolgoprudny, Russia, July 4–7, 2019, Tutorial Lectures, pp. 96–115, 2019.
  * [15] T. Standley, “Finding optimal solutions to cooperative pathfinding problems,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 24-1, pp. 173–178, 2010.
  * [16] A. Felner, R. Stern, S. Shimony, E. Boyarski, M. Goldenberg, G. Sharon, N. Sturtevant, G. Wagner, and P. Surynek, “Search-based optimal solvers for the multi-agent pathfinding problem: Summary and challenges,” in Proceedings of the International Symposium on Combinatorial Search, vol. 8-1, pp. 29–37, 2017.
  * [17] M. R. Ryan, “Exploiting subgraph structure in multi-robot path planning,” Journal of Artificial Intelligence Research, vol. 31, pp. 497–542, 2008.
  * [18] P. Surynek, “Towards optimal cooperative path planning in hard setups through satisfiability solving,” in Pacific Rim international conference on artificial intelligence, pp. 564–576, Springer, 2012.
  * [19] P. Surynek, “Makespan optimal solving of cooperative path-finding via reductions to propositional satisfiability,” arXiv preprint arXiv:1610.05452, 2016.
  * [20] J. Yu and S. M. LaValle, “Planning optimal paths for multiple robots on graphs,” in 2013 IEEE International Conference on Robotics and Automation, pp. 3612–3617, IEEE, 2013.
  * [21] E. Erdem, D. Kisa, U. Oztok, and P. Schüller, “A general formal framework for pathfinding problems with multiple agents,” in Proceedings of the AAAI Conference on Artificial Intelligence, vol. 27-1, pp. 290–296, 2013.
  * [22] P. Surynek, A. Felner, R. Stern, and E. Boyarski, “Efficient sat approach to multi-agent path finding under the sum of costs objective,” in Proceedings of the twenty-second european conference on artificial intelligence, pp. 810–818, 2016.
  * [23] R. Barták, N.-F. Zhou, R. Stern, E. Boyarski, and P. Surynek, “Modeling and solving the multi-agent pathfinding problem in picat,” in 2017 IEEE 29th International Conference on Tools with Artificial Intelligence (ICTAI), pp. 959–966, IEEE, 2017.
  * [24] D. Kornhauser, G. Miller, and P. Spirakis, “Coordinating pebble motion on graphs, the diameter of permutation groups, and applications,” in 25th Annual Symposium onFoundations of Computer Science, 1984., pp. 241–250, 1984.
  * [25] B. De Wilde, A. W. Ter Mors, and C. Witteveen, “Push and rotate: a complete multi-agent pathfinding algorithm,” Journal of Artificial Intelligence Research, vol. 51, pp. 443–492, 2014.
  * [26] P. Surynek, “A novel approach to path planning for multiple robots in bi-connected graphs,” in 2009 IEEE International Conference on Robotics and Automation, pp. 3613–3619, IEEE, 2009.
  * [27] R. C. Holte, M. B. Perez, R. M. Zimmer, and A. J. MacDonald, “Hierarchical a*: Searching abstraction hierarchies efficiently,” in AAAI/IAAI, Vol. 1, pp. 530–535, 1996.
  * [28] R. E. Korf, “Real-time heuristic search,” Artificial intelligence, vol. 42, no. 2-3, pp. 189–211, 1990.
  * [29] O. Goldreich, Finding the shortest move-sequence in the graph-generalized 15-puzzle is NP-hard.  Springer, 2011.
  * [30] V. Kann, “Maximum bounded 3-dimensional matching is max snp-complete,” Information Processing Letters, vol. 37, no. 1, pp. 27–35, 1991.
  * [31] R. Bellman, “Dynamic programming,” Science, vol. 153, no. 3731, pp. 34–37, 1966.
  * [32] E. L. Lawler and D. E. Wood, “Branch-and-bound methods: A survey,” Operations research, vol. 14, no. 4, pp. 699–719, 1966.


# Spinel LiGa5O8 prospects as ultra-wideband-gap semiconductor: band
structure, optical properties and doping.

Walter R. L. Lambrecht<EMAIL_ADDRESS>Department of Physics, Case
Western Reserve University, 10900 Euclid Avenue, Cleveland, Ohio 44106-7079,
USA

###### Abstract

LiGa5O8 in the spinel type structure is investigated as a potential ultra-
wide-band-gap semiconductor. The band structure is determined using the
quasiparticle self-consistent $GW$ method and the optical properties are
calculated at the Bethe Salpeter Equation level including electron-hole
interaction effects. The optical gap including exciton effects and an estimate
of the zero-point motion electron phonon coupling renormalizations is
estimated to be about 5.2$\pm 0.1$ eV with an exciton binding energy of about
0.4 eV. Si doping as potential $n$-type dopant is investigated and found to be
a promising shallow donor.

##  I Introduction

There is a need for ultra-wide-band-gap semiconductors both for optoelectronic
and high-power applications. The high-power applications benefit from a large
gap because it is found that the breakdown fields scale with the band gap.
Recently, there has been a lot of interest in $\beta$-Ga2O3, because of its
gap of $\sim$4.9 eV.[1, 2] In spite of such a large gap it can be doped with
Si, Ge or Sn to give semiconducting properties. This finding has sparked a
world-wide development of this material in epitaxial form. The development of
heterojunction devices involving even higher gap (AlxGa1-x)2O3 is making rapid
progress. For the higher Al-concentrations $n$-type doping is difficult but
$p$-type doping appears to be elusive so far even for Ga2O3. This is
attributed to formation of self-trapped hole polarons, related in turn to high
valence band effective masses. While there are some reports of low-level
$p$-type doping by means of defect complexes, these remain controversial and
the hole mobility is very low. [3, 4] $\beta$-Ga2O3 also has led to a renewed
interest in fundamental properties because of the unusual monoclinic structure
with both tetrahedral and octahedral Ga sites.

On the other hand, LiGaO2 (lithium gallate) has an even wider band gap and was
recently also predicted to be $n$-type dopable by Si or Ge. [5, 6] Like Ga2O3
bulk crystals can be grown and offer the possibility of homo-epitaxial growth.
It has a much simpler wurtzite based crystal structure with all Ga and Li in
tetrahedral coordinations and can be thought of as derived from ZnO by
replacing the group-II ion Zn by alternating group-I (Li) and group III Ga
ions. Further band gap tuning is possible by replacing Ga by Al and Li by Na
in the same $Pna2_{1}$ crystal structure or closely related structures.[7, 8]
The anisotropy of the material leads to interesting exciton splittings with
large exciton binding energies and recent calculations [9] were found to be in
excellent agreement with experiment[10, 11, 12] and give an exciton gap of 6
eV. It thus appears useful to develop epitaxial growth and doping of this
material to pursue it as an active semiconductor material and not merely an
optical material.

Figure 1: Crystal structure of LiGa5O8, blue spheres Li, light green,
octahedral Ga, dark green tetrahedral Ga, red, O

In the Li-Ga-O ternary system, another stable compound exists with formula
LiGa5O8. One can think of it as taking 4 LiGaO2 units and replacing 3 Li with
one Ga to maintain the charge balance. However, LiGa5O8 is known to have a
different crystal structure. It has a cubic spinel-type structure with Ga
occurring in both tetrahedral and octahedral sites and Li in octahedral
sites.[13] Its structure and related information can be found in [14, 15]. We
used the experimental crystal structure from the International
Crystallographic Database (ICSD),[16] with cubic lattice constant of 8.203 Å
and space group (No. 212) $P4_{3}32$ or $O^{6}$ as shown in Fig. 1.
Interestingly, a defective spinel phase also occurs in Ga2O3 as a metastable
phase and has been labeled the $\gamma$-phase [17]. It has recently been
observed to occur during certain processing steps in doped $\beta$-Ga2O3 and a
model for how the transition from $\beta$ to $\gamma$ phase can proceed via a
succession of defect formations was proposed by Hsien-Lien et al. [18, 19]
Since LiGa5O8 is a compound with composition in between Ga2O3 and LiGaO2 it is
of interest to determine its band structure and optical properties to assess
its potential as UWBG semiconductor. It has previously been studied as a
phosphorescent material when doped with Cr.[20, 21, 22] Very recently, LiGa5O8
was grown epitaxially and reported to exhibit $p$-type conduction.[23] The
relation to the $\gamma$-phase of Ga2O3, which is known to have a higher band
gap than the $\beta$-phase, provides additional interest. It suggests that
adding a small amount of Li may stabilize the $\gamma$-phase.

This paper presents results on the electronic band structure using the state-
of-the-art quasiparticle self-consistent $GW$ method and also calculates
optical properties, including excitonic effects. Since the unit cell already
has 56 atoms, (8 formula units) we can already obtain an idea of doping by
simply replacing one Ga by a Si. We thus also study a Li4Ga19SiO32 cell with
Si placed on either a tetrahedral or octahedral Ga site.

##  II Computational Methods

The calculations were performed using the Questaal suite,[24] which uses the
full-potential linearized muffin-tin orbital method(FP-LMTO) [25] to implement
both density functional theory (DFT) and many-body perturbation theory (MBPT)
approaches such as Hedin’s $GW$ method for quasiparticle band structures ($G$
for Green’s function and $W$ for screened Coulomb interaction)[26, 27] in a
quasiparticle self-consistent version called QS$GW$, [28, 29] and the Bethe
Salpeter Equation (BSE) approach for optical response calculations.[30, 31,
32] Specifically, we here use an extension of the QS$GW$ method in which the
screened Coulomb interaction is calculated including ladder diagrams, which
has been dubbed QS$G\hat{W}$,[32] and involves solving a BSE equation not just
in the ${\bf q}\rightarrow 0$ limit for optical response but for the mesh of
${\bf q}$ for which $W({\bf q},\omega)$ is needed. The muffin-tin-orbital
basis functions have atom centered spherical harmonic times smoothed Hankel
functions as envelope functions [33, 24], which are then replaced in its own
and every other sphere by an expansion in spherical harmonics times radial
solutions of the Schrödinger equation at a linearization energy and its energy
derivative which match in value and slope to the envelope functions. This
process is called augmentation. The expansion of the basis function centered
on one site about another is carried out by means of structure constants. The
basis set angular momentum cut-offs for first and second set of smoothed
Hankel function basis sets were $\ell_{max}$=3,3 for Ga and O and 3,2 for Li.
Augmentations inside the spheres were carried out up to a higher cut-off of
$\ell_{max}\leq 4$. In the QS$G\hat{W}$ calculations 32 valence bands and 4
conduction bands were included in the active space in which the ladder
diagrams are evaluated. The PBEsol functional was used for the initial DFT
calculations [34]. However, we did not minimize the structure within this
functional but used the experimental structural parameters and used the PBEsol
only as starting band structure for the subsequent $GW$ calculations. In the
QS$GW$ method, the energy dependent self-energy matrix $\Sigma_{ij}({\bf
k},\omega)$ in the basis of the Kohn-Sham eigenstates $\psi_{i}$ is replaced
by a Hermitian $\tilde{\Sigma}_{ij}({\bf
k})=\scriptstyle\frac{1}{2}\textstyle\Re{[\Sigma_{ij}({\bf
k},\epsilon_{i})+\Sigma_{ij}({\bf k},\epsilon_{j})]}$ and this acts as a non-
local exchange-correlation potential whose difference from the DFT exchange-
correlation potential is added to the Kohn-Sham Hamiltonian from which the
$\Sigma$ is calculated in the next iteration until convergence is reached. The
self-energy is calculated explicitly up to some maximum energy
$\omega_{max}=3.5$ Ry. The self-energy matrix $\Sigma_{ij}({\bf k},\omega)$ is
approximated by its diagonal above 3.0 Ry with an average value evaluated over
the range $3.0<E<3.5$ Ry. This procedure is used to facilitate interpolation
of the $GW$ bands to other k-points than the mesh on which the self-energy is
calculated explicitly. This is similar to a Wannier function interpolation
procedure where the muffin-tin orbitals themselves serve as Wannier functions.
For a full description of these technical aspects of the method we refer the
reader to Ref. [29].

Figure 2: Band structure of LiGa5O8 in GGA (red dashed) and QS$G\hat{W}$ (blue
solid) method. The Brillouin zone high-symmetry points follow the convention
for a simple cubic structure, $\Gamma=(0,0,0)$, $X=(\pi/a,0,0)$,
$S=(\pi/a,\pi,a,0)$, (labeled $M$ in [35, 36]), $R=(\pi/a,\pi/a,\pi/a)$ Figure
3: Zoom in near the valence band maximum of LiGa5O8 in QS$G\hat{W}$ method.

##  III Results

The band structure of LiGa5O8 is shown in Fig. 2 as obtained in the GGA and
QS$G\hat{W}$ methods. The conduction band minimum (CBM) occurs at $\Gamma$ and
the top valence band is very flat but with valence band maximum (VBM) at $S$
as shown in Fig. 3. The gap is thus slightly indirect with indirect gap of
5.72 eV and direct gap at $\Gamma$ at 5.84 eV. The gap is severely
underestimated by GGA. Further details on the gap are given in Table 1. The
difference between the QS$G\hat{W}$ and GGA gap is 3.24 eV and results from a
downward shift of the VBM by 1.43 eV and 1.81 upward shift of the CBM. These
individual band edge shifts are obtained by directly comparing the DFT and the
QS$GW$ band edges relative to the common reference potential, which in the
LMTO method is set by an average of the potentials at the muffin-tin radii. It
assumes that the self-consistent charge density and the potential apart from
the $GW$-self-energy terms are the same in both cases, which is an excellent
approximation.

From Table 1 we can see that the differences between the QS$GW$ and
QS$G\hat{W}$ gaps in this material is very small. This is somewhat unusual. In
most materials studied so far,[32] including the ladder diagrams reduces the
$GW$ self-energy by about 10-20 %. For LiGaO2, we already found a rather small
reduction (5 %) of the self-energy when including ladder diagrams but here the
reduction is only 1 %. However, the small effect of ladder diagrams found here
may be a result of including an insufficient number of bands in the active
space, i.e. the space spanned by the $N_{v}\times N_{c}\times N_{k}$ vertical
excitation basis states in the two-particle Hamiltonian solved in the BSE
step. $N_{v}$ is the number of occupied states included, $N_{c}$ the number of
empty conduction bands included and $N_{k}$ the number of k-points. We
estimate the reduction of the self-energy by

$1-\frac{E_{g}^{QSG\hat{W}}-E_{g}^{DFT}}{E_{g}^{QSGW}-E_{g}^{DFT}}.$ (1)

From the study of standard semiconductors, the guideline is to include the
upper set of anion-$p$-like valence bands and a number of conduction bands up
to an equivalent energy above the CBM. Unfortunately, for LiGa5O8 including
all O-$2p$ bands would require 96 valence bands and including at least one
band per Ga in the conduction band would require 20 conduction bands. This
presently exceeds the memory allocation requirements we can afford. To check
the importance of this convergence issue, we revisited our previous
calculations of LiGaO2 but with a higher number of conduction bands. Instead
of using all 52 occupied bands (20 Ga-$d$ like, 8 O-$2s$ and 24 O-$2p$) and
only 4 conduction bands which was the default used in [9] we now use 24
valence bands and 12 conduction bands. This gives a gap of 6.737 eV instead of
7.016 eV in [9]. With a QS$GW$ gap of 7.218 eV and GGA gap of 3.314 eV, the
$\Sigma_{QSG\hat{W}}$ is then estimated to be reduced relative to
$\Sigma_{QSGW}$ by 12 % instead of only 5 %. It is important to note that
while this reduces the quasiparticle gap of LiGaO2 by about 0.3 eV, the same
reduction of the screened Coulomb interaction also affects the exciton binding
energy. Calculating the exciton gap accurately requires a denser k mesh.
However, for these low lying excitations we can fortunately use a smaller
$N_{v}$ and $N_{c}$ but need a larger $N_{k}$ and an extrapolation as function
of $N_{k}$. This is different from the calculation of the ladder diagrams
which involves the overall effect of the screening and requires a higher
$N_{c}$. Using a similar extrapolation procedure as in [9] for the excitons we
then find an extrapolated exciton gap of 6.4 eV. This is in excellent
agreement with the result obtained in [9]. However, to compare with
experiment, we also need to include a zero-point motion electron-phonon
coupling correction, which was there estimated to be $-0.4$ eV. This finally
gave excellent agreement with the experimental value of 6.0 eV for the lowest
exciton gap.[12] In conclusion, the main result of that paper of an exciton
gap at 6.0 eV and the qualitative analysis of the exciton states presented in
that paper holds when we here use a better converged QS$G\hat{W}$
calculations. However, the quasiparticle gap and the exciton binding energy
are both reduced by about 0.3 eV. So, it does represent a correction to those
previous results but one which is not so easy to test experimentally as no
independent measurements of the exciton binding energy or the quasiparticle
gap are available but only the resulting optical gap or exciton energy. This
re-analysis of the LiGaO2 case suggests that the QS$G\hat{W}$ gap in Table 1
is still overestimated. Using a reduction factor due to ladders of about 12 %
as suggested by the re-evaluation of LiGaO2 we would obtain a indirect gap of
5.4$\pm 0.1$ eV and direct gap at $\Gamma$ of 5.5 eV.

Table 1: Band gaps of LiGa5O8 in eV | GGA | QS$GW$ | QS$G\hat{W}$
---|---|---|---
indirect ($S-\Gamma$) | 2.48 | 5.75 | 5.72
direct $\Gamma$ | 2.58 | 5.87 | 5.84
Figure 4: Real and imaginary part of the dielectric function. Figure 5:
Exciton gap as function of k-mesh.

The dielectric function is shown in Fig. 4 in the independent particle
approximation (IPA) and the BSE approximation. We can see a significant effect
from including the electron-hole interaction effects and excitons with a large
exciton binding energy. We caution that the dielectric constant
$\varepsilon(\omega=0)$ obtained here is likely overestimated as is the
intensity of $\varepsilon_{2}(\omega)$ because of difficulties in calculating
the contribution of the self-energy to the velocity matrix elements but the
peak positions should be accurate. In this calculation we only used $N_{k}=2$
in a $N_{k}\times N_{k}\times N_{k}$ mesh for the BSE equation. To obtain the
excitons more accurately, we calculated the lowest exciton for finer k-meshes
with $N_{k}=4,6$ but with smaller number of bands 10 valence bands instead of
32 while keeping 4 conduction bands. We can see that reducing the number of
valence bands makes a difference of only 0.01 eV with slightly lower gaps
obtained with more valence bands included. Extrapolating the gap as function
of $1/N_{k}^{3}$, i.e. the total number of mesh points in the BZ, as we did in
[9] we find an extrapolated lowest exciton at 5.48 eV, close to our calculated
gap for $N_{k}=6$. This would amount to a exciton binding energy of $\sim$0.4
eV. The BSE calculations here only include direct transitions. Therefore we
deduce the exciton binding energy with respect to the direct gap at $\Gamma$.
One may expect a similar binding energy for an exciton related to the indirect
gap which is $\sim$0.1 eV smaller. Neither the quasiparticle gap nor the
exciton gap here include electron-phonon coupling renormalization. In LiGaO2,
the latter was found to be $-0.3\pm 0.1$ eV. A similar value is expected here
and would reduce our direct gap of Table 1 from 5.87 to 5.6$\pm$0.1 eV
(quasiparticle ) and 5.2.$\pm$0.1 eV (exciton) and another 0.1 eV lower for
the indirect quasiparticle and exciton gaps. Both the quasiparticle gap and
exciton binding energy could be somewhat lower if we assume a larger reduction
of $\hat{W}$ relative to $W$ due to ladder diagrams but as mentioned earlier
for LiGaO2, we expect these two errors to compensate each other and give a
similar optical gap. This gap is slightly higher than the gap of 5.0 eV
reported [17] for $\gamma$-Ga2O3.

(a)  (b)

Figure 6: Band structure of SiLi4Ga19O32 with Si in tetrahedral site (a),
octahedral site (b). Red indicates the Si-dopant contribution to the bands,
bands without Si are light grey.

(a) (b)

Figure 7: Bond lengths (Å) near Si relaxed structure for: (a) tetrahedral, (b)
octahedral site. The blue sphere is Si, the green Ga and the teal one is Li.
The longest Si-O bond is indicated by the crosses in (a).

To evaluate the $n$-type dopability of the system, we replaced one Ga atom in
the 56 atom cell with Si. Fig.6(a) shows the band structure for the
tetrahedral Si site after relaxation of the atomic positions keeping the
volume of the cell fixed and in the GGA. One can see that the extra valence
electron essentially starts filling the conduction band but no levels occur in
the gap. The gap in fact slightly increased from 2.58 eV in pure LiGa5O8 in
the GGA to 2.74 eV. This is partly due to the relaxation. Without relaxation
the gap would have decreased slightly to 2.39 eV. The Si contribution to the
bands in this figure is indicated by the red color while the background color
of the bands without Si contribution is grey. The relaxation is significant.
As can be seen in Fig. 7, the nearest neighbor O move inward toward the Si but
with also a distortion and symmetry lowering making one Si-O bond longer than
the other three. The bond lengths change from 1.776 Å for the Ga-O tetrahedral
bond to 1.675 Å (3 bonds) and 1.695 Å (one Si-O) bond. The Ga atoms bonded to
these nearest neighbor O atoms then have a larger bond length again by about
0.1 Å than the original octahedral Ga-O bond length of 2.05 Å. In spite of
these relaxations there clearly is no deep level formation. This indicates
that Si would be a shallow donor and could lead to $n$-type doping.

We have also calculated Si on an octahedral site with similar results. First,
we find that the relaxed total energies were within error bar
indistinguishable. So Si is likely to occupy either octahedral or tetrahedral
sites with equal likelihood. Secondly, we find no states in the gap in either
case, indicating that Si is a shallow donor on both sites as seen in Fig.
6(b). Similar to the tetrahedral case, we find a slight increase in gap due to
the Si replacement to 2.688 eV in GGA. Two of the octahedral bonds around Si
are 1.814 Å, two are 1.864 Å, and two are 1.784Å. The Ga-O octahedral bonds in
the host are 1.982 Å. Thus, there is again an inward relaxation of nearby O
toward the Si.

While we have here only studied the Si doping via the GGA band structure, we
expect that the shallow donor character will stay valid at the QS$GW$ level. A
similar situation was found recently for LiAlO2.[7]

Of course, a full evaluation of doping possibilities will also require a study
of native defects and compensation issues. Work along those lines is in
progress and will be published elsewhere.

##  IV Conclusion

In summary, LiGa5O8 in the spinel type cubic structure is a good candidate
UWBG semiconductor. It is predicted to have slightly higher gap than
$\gamma$-Ga2O3 of about 5.3 eV but not as high as LiGaO2 and appears to be
$n$-type dopable by Si on either its tetrahedral or octahedral Ga site. Its
cubic structure may have advantages.

###### Acknowledgements.

This work made use of the High Performance Computing Resource in the Core
Facility for Advanced Research Computing at Case Western Reserve University
and was supported by the U.S. Air Force Office of Scientific Research (AFOSR)
under grant no. FA9550-22-1-0201. I thank Hongping Zhao for communicating her
results on epitaxial growth of LiGa5O8 prior to publication, which initiated
my interest in this material. Data Availability: The data that supports the
findings of this study are available within the article

## References

  * Matsumoto  _et al._ [1974] T. Matsumoto, M. Aoki, A. Kinoshita, and T. Aono, Absorption and Reflection of Vapor Grown Single Crystal Platelets of $\beta$-Ga 2 O 3, Jap. J. Appl. Phys. 13, 1578 (1974).
  * Sasaki  _et al._ [2013] K. Sasaki, M. Higashiwaki, A. Kuramata, T. Masui, and S. Yamakoshi, ”{MBE} grown ga2o3 and its power device applications ”, J. Cryst. Growth 378, 591 (2013), the 17th International Conference on Molecular Beam Epitaxy.
  * Chi  _et al._ [2023] Z. Chi, C. Sartel, Y. Zheng, S. Modak, L. Chernyak, C. M. Schaefer, J. Padilla, J. Santiso, A. Ruzin, A.-M. Gonçalves, J. von Bardeleben, G. Guillot, Y. Dumont, A. Pérez-Tomás, and E. Chikoidze, Native defects association enabled room-temperature p-type conductivity in $\beta$-Ga2O3, Journal of Alloys and Compounds 969, 172454 (2023).
  * Chikoidze  _et al._ [2019] E. Chikoidze, C. Sartel, H. Mohamed, I. Madaci, T. Tchelidze, M. Modreanu, P. Vales-Castro, C. Rubio, C. Arnold, V. Sallet, Y. Dumont, and A. Pérez-Tomás, Enhancing the intrinsic p-type conductivity of the ultra-wide bandgap Ga2O3 semiconductor, J. Mater. Chem. C 7, 10231 (2019).
  * Boonchun  _et al._ [2019] A. Boonchun, K. Dabsamut, and W. R. L. Lambrecht, First-principles study of point defects in LiGaO2, Journal of Applied Physics 126, 155703 (2019).
  * Dabsamut  _et al._ [2020] K. Dabsamut, A. Boonchun, and W. R. L. Lambrecht, First-principles study of n- and p-type doping opportunities in LiGaO2, Journal of Physics D: Applied Physics (2020).
  * Popp and Lambrecht [2022] P. Popp and W. R. L. Lambrecht, Quasiparticle self-consistent $GW$ band structures and phase transitions of ${\mathrm{LiAlO}}_{2}$ in tetrahedrally and octahedrally coordinated structures, Phys. Rev. Mater. 6, 104605 (2022).
  * Radha  _et al._ [2021] S. K. Radha, A. Ratnaparkhe, and W. R. L. Lambrecht, Quasiparticle self-consistent $GW$ band structures and high-pressure phase transitions of $\mathrm{Li}\mathrm{Ga}{\mathrm{O}}_{2}$ and $\mathrm{Na}\mathrm{Ga}{\mathrm{O}}_{2}$, Phys. Rev. B 103, 045201 (2021).
  * Dadkhah  _et al._ [2023] N. Dadkhah, W. R. L. Lambrecht, D. Pashov, and M. van Schilfgaarde, Improved quasiparticle self-consistent electronic band structure and excitons in $\beta\text{$-$}{\mathrm{LiGaO}}_{2}$, Phys. Rev. B 107, 165201 (2023).
  * Trinkler  _et al._ [2017] L. Trinkler, A. Trukhin, B. Berzina, V. Korsaks, P. Ščajev, R. Nedzinskas, S. Tumėnas, M. Chou, L. Chang, and C.-A. Li, Luminescence properties of LiGaO2 crystal, Optical Materials 69, 449 (2017).
  * Trinkler  _et al._ [2022] L. Trinkler, V. Pankratov, A. Trukhin, B. Berzina, M. Chou, and L. Chang, Anisotropic photoluminescence of $\beta$-LiGaO2 crystal, Optical Materials 132, 112856 (2022).
  * Tuménas  _et al._ [2017] S. Tuménas, P. Mackonis, R. Nedzinskas, L. Trinkler, B. Berzina, V. Korsaks, L. Chang, and M. Chou, Optical properties of lithium gallium oxide, Applied Surface Science 421, 837 (2017).
  * Joubert  _et al._ [1963] J. Joubert, M. Brunel, A. Waintal, and A. Durif, Etude cristallographique du gallate de lithium et de sa solution solide avec l‘aluminate, Comptes Rendus Hebdomadaires des Seances de l’Academie des Sciences (1884 - 1965) 256, 5324–5326 (1963).
  * [14] Data retrieved from Materials Project for LiGa5O8 (mp-28146) from database version v2022.10.28.
  * Jain  _et al._ [2013] A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder, and K. A. Persson, Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Materials 1, 011002 (2013).
  * [16] International Crystallographic Data Base, id=33716.
  * Biswas and Nishinaka [2022] M. Biswas and H. Nishinaka, Thermodynamically metastable $\alpha$-, $\epsilon$\- (or $\kappa$-), and $\gamma$-Ga2O3: From material growth to device applications, APL Materials 10, 060701 (2022), https://pubs.aip.org/aip/apm/article-pdf/doi/10.1063/5.0085360/16490464/060701_1_online.pdf .
  * Huang  _et al._ [2023a] H.-L. Huang, C. Chae, J. M. Johnson, A. Senckowski, S. Sharma, U. Singisetti, M. H. Wong, and J. Hwang, Atomic scale defect formation and phase transformation in Si implanted $\beta$-Ga2O3, APL Materials 11, 061113 (2023a), https://pubs.aip.org/aip/apm/article-pdf/doi/10.1063/5.0134467/17941303/061113_1_5.0134467.pdf .
  * Huang  _et al._ [2023b] H.-L. Huang, J. M. Johnson, C. Chae, A. Senckowski, M. H. Wong, and J. Hwang, Atomic scale mechanism of $\beta$ to $\gamma$ phase transformation in gallium oxide, Applied Physics Letters 122, 251602 (2023b), https://pubs.aip.org/aip/apl/article-pdf/doi/10.1063/5.0156009/18006597/251602_1_5.0156009.pdf .
  * De Clercq  _et al._ [2017] O. Q. De Clercq, L. I. D. J. Martin, K. Korthout, J. Kusakovskij, H. Vrielinck, and D. Poelman, Probing the local structure of the near-infrared emitting persistent phosphor LiGa5O8:Cr3+, J. Mater. Chem. C 5, 10861 (2017).
  * Sousa and Carvalho [2020] O. Sousa and I. Carvalho, Theoretical study of the structural, energetic, electronic and magnetic properties of the host matrix LiGa5O8 doped with Cr3+, Journal of Solid State Chemistry 289, 121472 (2020).
  * Huang  _et al._ [2018] W. Huang, X. Gong, R. Cui, X. Li, L. Li, X. Wang, and C. Deng, Enhanced persistent luminescence of LiGa5O8:Cr3+ near-infrared phosphors by codoping Sn4+, Journal of Materials Science: Materials in Electronics 29, 10535 (2018).
  * [23] K. Zhang, V. G. T. Vangipuram, H.-L. Huang, J. Hwang, and H. Zhao, Discovery of a Robust P-Type Ultrawide Bandgap Oxide Semiconductor: LiGa5O8, Advanced Electronic Materials n/a, 2300550, https://onlinelibrary.wiley.com/doi/pdf/10.1002/aelm.202300550 .
  * Pashov  _et al._ [2019] D. Pashov, S. Acharya, W. R. Lambrecht, J. Jackson, K. D. Belashchenko, A. Chantis, F. Jamet, and M. van Schilfgaarde, Questaal: A package of electronic structure methods based on the linear muffin-tin orbital technique, Computer Physics Communications , 107065 (2019).
  * Kotani and van Schilfgaarde [2010] T. Kotani and M. van Schilfgaarde, Fusion of the LAPW and LMTO methods: The augmented plane wave plus muffin-tin orbital method, Phys. Rev. B 81, 125117 (2010).
  * Hedin [1965] L. Hedin, New method for calculating the one-particle green’s function with application to the electron-gas problem, Phys. Rev. 139, A796 (1965).
  * Hedin and Lundqvist [1969] L. Hedin and S. Lundqvist, Effects of electron-electron and electron-phonon interactions on the one-electron states of solids, in  _Solid State Physics, Advanced in Research and Applications_ , Vol. 23, edited by F. Seitz, D. Turnbull, and H. Ehrenreich (Academic Press, New York, 1969) pp. 1–181.
  * van Schilfgaarde  _et al._ [2006] M. van Schilfgaarde, T. Kotani, and S. Faleev, Quasiparticle Self-Consistent $GW$ Theory, Phys. Rev. Lett. 96, 226402 (2006).
  * Kotani  _et al._ [2007] T. Kotani, M. van Schilfgaarde, and S. V. Faleev, Quasiparticle self-consistent GW method: A basis for the independent-particle approximation, Phys.Rev. B 76, 165106 (2007).
  * Onida  _et al._ [2002] G. Onida, L. Reining, and A. Rubio, Electronic excitations: density-functional versus many-body Green’s-function approaches, Rev. Mod. Phys. 74, 601 (2002).
  * Cunningham  _et al._ [2018] B. Cunningham, M. Grüning, P. Azarhoosh, D. Pashov, and M. van Schilfgaarde, Effect of ladder diagrams on optical absorption spectra in a quasiparticle self-consistent $\mathit{GW}$ framework, Phys. Rev. Materials 2, 034603 (2018).
  * Cunningham  _et al._ [2023] B. Cunningham, M. Gruening, D. Pashov, and M. van Schilfgaarde, QS$G\hat{W}$: Quasiparticle Self consistent GW with ladder diagrams in W (2023), arXiv:2302.06325 [cond-mat.mtrl-sci] .
  * Bott  _et al._ [1998] E. Bott, M. Methfessel, W. Krabs, and P. C. Schmidt, Nonsingular Hankel functions as a new basis for electronic structure calculations, J. Math. Phys. 39, 3393 (1998).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett. 100, 136406 (2008).
  * [35] Bilbao Crystallographic Server, https://www.cryst.ehu.es.
  * Aroyo  _et al._ [2014] M. I. Aroyo, D. Orobengoa, G. de la Flor, E. S. Tasci, J. M. Perez-Mato, and H. Wondratschek, Brillouin-zone database on the Bilbao Crystallographic Server, Acta Crystallographica Section A 70, 126 (2014).