# Nearly-frustration-free ground state preparation

Matthew Thibodeau<EMAIL_ADDRESS>Bryan K. Clark<EMAIL_ADDRESS>Department of Physics, University of Illinois at Urbana-Champaign, IL 61801,
USA  IQUIST and Institute for Condensed Matter Theory and NCSA Center for
Artificial Intelligence Innovation, University of Illinois at Urbana-
Champaign, IL 61801, USA

###### Abstract

Solving for quantum ground states is important for understanding the
properties of quantum many-body systems, and quantum computers are potentially
well-suited for solving for quantum ground states. Recent work [14] has
presented a nearly optimal scheme that prepares ground states on a quantum
computer for completely generic Hamiltonians, whose query complexity scales as
$\delta^{-1}$, i.e. inversely with their normalized gap. Here we consider
instead the ground state preparation problem restricted to a special subset of
Hamiltonians, which includes those which we term “nearly-frustration-free”:
the class of Hamiltonians for which the ground state energy of their block-
encoded and hence normalized Hamiltonian $\alpha^{-1}H$ is within $\delta^{y}$
of -1, where $\delta$ is the spectral gap of $\alpha^{-1}H$ and $0\leq y\leq
1$. For this subclass, we describe an algorithm whose dependence on the gap is
asymptotically better, scaling as $\delta^{y/2-1}$, and show that this new
dependence is optimal up to factors of $\log\delta$. In addition, we give
examples of physically motivated Hamiltonians which live in this subclass.
Finally, we describe an extension of this method which allows the preparation
of excited states both for generic Hamiltonians as well as, at a similar
speedup as the ground state case, for those which are nearly frustration-free.

##  1 Introduction

Steady progress in the development of quantum computing hardware has spurred
interest in possible applications of these machines. One natural candidate is
the ground state preparation problem: given a many-body Hamiltonian $H$,
prepare its ground state in a quantum register. Given a quantum many-body
system, physicists are often interested in the properties of its ground state
because this state predominates at low temperature; preparing ground states
using a quantum computer is an important first step to studying them
synthetically, with full control over their generating Hamiltonian.

To prepare the ground state of $H$ in a quantum register, one must distinguish
it from all other states in the spectrum. Since, by definition, the ground
state is separated in energy from the rest of the spectrum by the spectral gap
$\Delta$, it should be no great surprise that $\Delta$ is a key parameter
controlling the complexity of quantum algorithms that prepare ground states.
For example, both phase estimation [13] combined with amplitude amplification
[2], as well as quantum singular value tranformation (QSVT) [9], can find the
ground state in time $O(1/\Delta)$. The latter is more efficient with respect
to the error $\epsilon$, with a cost which goes as $\log(\epsilon)$ instead of
$(1-\epsilon)$ [7]. In fact, the QSVT is known to be optimal with respect to
query complexity for preparing the ground state in a quantum register for
generic Hamiltonians [14]. The query complexity measures the number of times
that $H$ must be accessed and roughly corresponds to the circuit depth. The
optimality of this procedure comes from a reduction to Grover’s algorithm
[10].

For QSVT, beyond $\Delta$ and $\epsilon$, there is an additional dependence on
the overlap $\gamma=\braket{\phi}{\mu}$ with the starting ansatz, and an
additional term $\alpha$ related to the cost of loading (i.e. block-encoding)
$H$ into a quantum computer; roughly, for generic local lattice models, a
block-encoding with $\alpha\propto N_{\textrm{sites}}$ is easily obtained
[17]. In fact, both in our algorithms and those previously given, $\Delta$
always appears in the complexity normalized to $\Delta/\alpha\sim\Delta/N$;
Section 3 covers the details of this normalization.

In many situations, one might know some special properties that $H$ may
satisfy which then may evade this general optimality result. Examples of
potential properties could include $k$-locality, frustration-free ground
states, or stoquastic Hamiltonians [3, 19]. It is useful, then, to find
quantum circuits that prepare the ground state of $H$ after assuming that $H$
is a member of some restricted subclass of all hermitian operators, or to
bound the minimum depth of such circuits.

One of the most important of these subclasses is that of frustration-free
Hamiltonians, where the ground state of $H$ is simultaneously a ground state
of all the operators which sum up to make $H$. A naive application of the
algorithm from ref. [14] still scales as $1/\Delta$ even in the frustration
free case. Nonetheless, there is a straightforward way to do better for many
frustration free Hamiltonians using spectral gap amplification [19] or uniform
spectral amplification [15]. These two approaches take Hamiltonians which are
frustration free when decomposed into operator terms which are projectors or
linear combination of unitaries respectively and generates a new Hamiltonian
with the same ground state which has an $O(1)$ gap. Given a standard (block-
encoded) description of our original Hamiltonian, one query of the new
Hamiltonian takes $O(1/\sqrt{\Delta})$ queries of the original Hamiltonian.
Using the optimal ground-state generation on the new Hamiltonian then takes
$O(1)$ time (given that the gap is $O(1)$) leading to a total number of
queries (and circuit depth) which scales as $O(1/\sqrt{\Delta})$,
significantly improving over the generic result of $1/\Delta$. For the case of
a local Hamiltonian, this also gives a speedup from $O(N)$ to $O(\sqrt{N})$
via the $\alpha$-dependence of the complexity. The apparent violation of the
optimality relation is not surprising, as the optimality theorem of [14]
relies on a specific Hamiltonian $H_{G}$ that is extremely frustrated.

The improvement of the frustration free case from $1/\Delta$ to
$1/\sqrt{\Delta}$ motivates asking whether there is some intermediate class of
$H$ between fully frustration-free and the generic case where a speedup might
be achieved, as well as a corresponding proof of optimality. Further, while
this all applies to ground state problems, there are analogous questions one
can ask about excited states.

##  2 Summary of Results

In this work, we address the question of finding a method that prepares a
ground state of $H$ for a class of Hamiltonians that are near frustration-
free, giving a scaling which continuously tunes between the $1/\Delta$ of the
generic case to the $1/\sqrt{\Delta}$ for the frustration-free case. While our
primary goal will be to show this for ground state problems, we will also
extend the work of ref. [14] to preparing low-lying excited states and
furthermore show how they can similarly be accelerated in nearly frustration-
free situations. To accomplish these goals, the property of nearly
frustration-free Hamiltonians that we use is that the ground state energy of
the Hamiltonian is close to the frustration-free energy. We then provide a
lower-bound on the complexity of preparing ground states of such Hamiltonians,
showing that the scaling with $\Delta$ of our algorithm is optimal up to
logarithmic factors.

To be more precise, we need to use the language of block encoding [4, 17]. A
block encoding of $H$ is a unitary operator $U$ that has a scaled copy of the
Hamiltonian $\alpha^{-1}H$ as its upper-left corner. The existence of a block-
encoding $U$ of $H$ is a very powerful tool because it allows the application
quantum signal processing [17, 9, 11, 16, 5], which lets us transform $H$ by a
suitable polynomial $p(x)$, which acts on the spectrum of $H$ by
$p(H)=\sum_{\lambda}p(\lambda)\ket{\lambda}\bra{\lambda}$.

Now let us define a parameter $y$ in terms of the other parameters of our
problem:

######  Definition 1.

Let $H_{0}$ be a Hermitian operator on $n$ qubits with ground state energy
$\mu_{0}$ and gap $\Delta$, and let $U$ be a block-encoding of
$\alpha^{-1}H_{0}$ (with $\alpha\geq||H_{0}||$, the spectral norm of $H_{0}$)
using $m$ extra qubits. Set $H=\alpha^{-1}H_{0}$, and let
$\delta=\alpha^{-1}\Delta$ and $\mu=\alpha^{-1}\mu_{0}$; note that
$-1\leq\mu\leq 1$. Given these parameters, there is a $y$ such that $0\leq
y\leq 1$ and $\mu\leq-1+2\delta^{y}$.

The relationship between $\alpha$, spectral norm $||H_{0}||$, and the ground
state energy $\mu$ will play an important role in the efficiency of our
algorithm. In particular, the actual promise we require of $H$ is that there
exists a block-encoding where $\alpha$ is close to the magnitude of the ground
state energy $|\mu|$. This closeness is encoded in the value of the parameter
$y$ which specifies the distance between $\mu$ and its minimum possible value
of $-1$, and we develop an algorithm for the ground state preparation problem
such that the query complexity to $U$ is controlled by $y$.

######  Proposition 1.

Let $H_{0}$ and $U$ be as specified in Definition 1, and let $\ket{\phi}$
satisfy $|\braket{\phi}{\mu_{0}}|=\gamma$. Then the query complexity of
preparing the ground state of $H_{0}$ to precision $\epsilon$ starting with
the ansatz $\ket{\phi}$, as measured by queries to $U$ and $U^{\dagger}$, is

$\displaystyle
O\left(\frac{\delta^{y/2-1}}{\gamma}\log^{5/2}\frac{1}{\epsilon}\log\frac{1}{\delta}\right)$
(1)

To accomplish this, as in [14], we will also apply a polynomial of $H$,
$p(H)\approx\ket{\mu}\bra{\mu}$, which is approximately the projection onto
the ground state $\ket{\mu}$. Note that the polynomial we apply and that of
ref [14], though closely related, differs in an important way which is
necessary to achieve this improved algorithm. In particular, $p(x)$ is
necessarily a quickly-oscillating function of $x$ near the ground state energy
so as to filter out the first excited state. The key improvement for our
polynomial is to shift the location of this oscillation from $x=0$ to outer
edges of the domain $[-1,1]$ which allows a reduction in the required degree
of $p$ and therefore a reduction in the query complexity.

As an example, in the case that $H_{0}$ is frustration-free when decomposed
into a linear combination of unitaries, there exists a block-encoding $U$ such
that $\alpha=|\mu_{0}|$, so that $\mu=-1$ and $y=1$. Then, our complexity
results agree with that achieved by combining uniform spectral amplification
[15] with the standard approach of ref. [14], giving that the ground state can
be prepared with complexity $O(1/\sqrt{\delta})$. In addition, when $\mu$ is
close enough to $-1$ and $y>0$, we find an intermediate scaling between
$1/\sqrt{\delta}$ and $1/\delta$.

In addition to providing an improved algorithm, we also show that the
asymptotic dependence of the query complexity on $\delta$ is optimal, up to a
factor of $\log\delta$, for any Hamiltonians which satisfy only the properties
of Definition 1.

######  Proposition 2.

Let $H_{0}$ be as specified in Definition 1. If $\gamma=\Omega(1)$ and
$\delta\to 0+$, then the query complexity of preparing its ground state is
$\Omega(\delta^{y/2-1})$.

This is proved through a reduction to Grover’s algorithm, using a suitably
transformed version of the Hamiltonian $H(1/2)$ in Theorem 10 of [14].

##  3 Block Encodings, Polynomial Approximation and Degree Reduction

As stated above, to prepare the ground state we will form a quantum circuit
that encodes a projector onto $\ket{\mu_{0}}$, accurate up to some error
$\epsilon$. To do so we will peform quantum signal processing on a block-
encoding of $\alpha^{-1}H_{0}$ for some normalization $\alpha$. We say that
$U$ is an $(\alpha,m,\epsilon)$-block-encoding of $H_{0}$ if

$\displaystyle||(\bra{0^{m}}\otimes I)U(\ket{0^{m}}\otimes
I)-\alpha^{-1}H_{0}||_{2}\leq\epsilon$ (2)

i.e. that $U$ is a unitary operator using $m$ extra qubits beyond those used
by $H=\alpha^{-1}H_{0}$, and the restriction of $U$ to the
$\ket{0^{m}}$-flagged subspace is $\epsilon$-close in operator norm to $H$.
Given a block-encoding, we can form a polynomial transformations of $H$ using
the following theorem:

######  Theorem 1 ([9], Theorem 2).

Let U be an $(\alpha,m,0)$-block-encoding of a Hermitian matrix $H_{0}$. Let
$p\in\mathbf{R}[x]$ be a degree-$l$ even or odd real polynomial with $l>0$ so
that $|p(x)|\leq 1$ for any $x\in[-1,1]$. Then there exists a
$(1,m+1,0)$-block-encoding $\widetilde{U}$ of $p(\alpha^{-1}H_{0})$ using $l$
queries of $U$, $U^{\dagger}$, and $O((m+1)l)$ other primitive quantum gates.

The theorem makes clear that the degree $l$ of the polynomial transformation
$p$ is the key driver of circuit complexity, and the idea of this work is to
show that for some class of $H$, the required $l$ for projection into the
ground state can be systematically reduced.

The basic tools of quantum signal processing are these polynomial
transformations of the input $H$, most often polynomial approximations of
analytic functions. That holds true for our application, where the function at
hand is a regulated version of the unit step function which we will take to be
$\Theta_{k}(x)=\tfrac{1}{2}(1+\operatorname{erf}kx)$ for $k>0$.

The key to achieving a speedup over the general case is by exploiting a
feature of Chebyshev approximations on the interval $x\in[-1,1]$: quickly-
oscillating functions can be approximated (to a given precision) more
efficiently when their oscillation occurs near the endpoints $x=\pm 1$. In
particular, the Chebyshev approximation of the shifted step function
$\Theta_{k}(x+1)$ can be truncated at a lower degree than that of
$\Theta_{k}(x)$ while maintaining the same desired precision. We will choose
$k$ such that, given an error $\epsilon$, the error function is within
$\epsilon$ of the step function outside of a region of width $\delta$:

$\displaystyle|\Theta_{k}(x)-\Theta(x)|$
$\displaystyle\leq\epsilon\quad\quad|x|\geq\delta/2$ (3)

The details of the approximation are made precise by the following Lemma,
which adapts Lemma 16 from [15].

######  Lemma 1.

Fix $\epsilon,\delta>0$ and $\eta$ such that $|\eta|<1-\delta$, and let
$k=\frac{\sqrt{2}}{\delta}\log^{1/2}\frac{1}{2\pi\epsilon^{2}}$. Then the
shifted error function $f(x)=\operatorname{erf}k(x-\eta)$ can be approximated
to error $\epsilon$ on the interval $x\in[-1,1]$ with a polynomial $p(x)$ such
that $|p(x)|\leq 1$ for $|x|\leq 1$ and $p$ is of order
$n=O(\frac{\sqrt{1-|\eta|}}{\delta}\log\frac{\sqrt{1-|\eta|}}{\delta}\log^{3/2}\frac{1}{\epsilon})=\widetilde{O}(\frac{\sqrt{1-|\eta|}}{\delta})$,
where the $\widetilde{O}$ notation hides polylogarithmic factors in
$\epsilon,\delta$ and $\eta$.

###### Proof.

Expand $f(x)=\sum_{j}a_{j}T_{j}(x)$ in the Chebyshev basis; we will form
$p(x)$ by truncating this series. As stated in [15], the truncation error for
$p_{n}(x)=\sum_{j=0}^{n}a_{j}T_{j}(x)$ is

$\displaystyle\max_{x\in[-1,1]}|f(x)-p(x)|\leq\frac{2M\rho^{-n}}{\rho-1}=O(\epsilon)\quad\quad
M=\max_{z\in E_{\rho}}|f(x)|$ (4)

for any $\rho>1$, where $\rho$ is the parameter of the Bernstein ellipse
$E_{\rho}=\\{z:z=\frac{1}{2}(\rho e^{i\theta}+\rho^{-1}e^{-i\theta}\\}$. Using
the same error estimates given in [15], we may bound

$\displaystyle M\leq\max_{z\in
E_{\rho}}\exp\operatorname{Re}\left(-(k(z-\eta))^{2}\right)$ (5)

Over the range $0\leq\theta<2\pi$, the exponent
$\alpha=\operatorname{Re}(k(z+\eta))^{2}$ achieves its maximum value of

$\displaystyle\alpha_{0}=\frac{k^{2}}{4\rho^{2}(1+\rho^{4})}(\rho^{2}-1)^{2}(1-2\eta^{2}\rho^{2}+\rho^{4})$
(6)

If we choose $\rho=1+a$ and expand around small $a$, then we find that
$\alpha_{0}=2k^{2}(1-\eta^{2})a^{2}+O(a^{3})$. Requiring that
$\alpha_{0}=O(1)$ fixes the magnitude of $a$ to

$\displaystyle
a=O\left(\frac{1}{\sqrt{k^{2}(1-\eta^{2})}}\right)=O\left(\frac{\delta}{\sqrt{1-\eta^{2}}\log^{1/2}1/\epsilon}\right)$
(7)

after substituting in $k=O(\delta^{-1}\log^{1/2}(1/\epsilon))$. Note that so
long as $|\eta|<1-\delta$, which is always true, $a$ is indeed asymptotically
small and scales at least as strongly as $O(\sqrt{\delta})$. We are interested
in the case where $|\eta|$ is near 1, so we may set
$\sqrt{1-\eta^{2}}=O(\sqrt{1+\eta}\sqrt{1-\eta})=O(\sqrt{1-|\eta|})$. Thus,
the approximation error from (4) is

$\displaystyle\epsilon=O(a^{-1}e^{-na})\implies n$
$\displaystyle=O\left(\frac{\sqrt{1-|\eta|}}{\delta}\log\frac{\sqrt{1-|\eta|}}{\delta}\log^{3/2}\frac{1}{\epsilon}\right)$
(8) $\displaystyle=\widetilde{O}\left(\frac{\sqrt{1-|\eta|}}{\delta}\right)$
(9)

as desired.

∎

To make use of this result, we will always put the location of the shift
between the ground state energy $\mu$ and the first excited state
$\mu+\delta$, so that $\eta=\mu+\delta/2$. In particular, if
$\mu=-1+O(\delta)$, then the degree is
$n=\widetilde{O}(\frac{1}{\sqrt{\delta}})$. Essentially, if the location $\mu$
of the unit step is close to $-1$ (relative to the width $\delta$) then the
dependence of $n$ on $\delta$ can be reduced by up to a square-root factor
relative to the unshifted case where $n=\widetilde{O}(\frac{1}{\delta})$.

Figure 1: In A, the spectrum of $H$, the Hamiltonian whose ground state
$\ket{\mu_{0}}$ we would like to prepare. To block-encode, $H$ must be
normalized to fit its spectrum in $[-1,1]$ as shown in B and C, along with the
polynomial functions $p_{1}(H/\alpha_{1})$ and $p_{2}(H/\alpha_{2})$ that
increase the effective gap to $\Omega(1)$. In B, the normalization
$\alpha_{1}$ is such that $\mu_{1}=\mu_{0}/\alpha_{1}$ is far from $-1$
compared to the normalized gap $\delta_{1}=\Delta/\alpha_{1}$, while in C
another normalization $\alpha_{2}$ gives us $\mu_{2}<-1+\sqrt{\delta_{2}}$.
Because Chebyshev polynomials oscillate faster near the endpoints of the unit
interval, $p_{2}$ requires a lower degree than does $p_{1}$. D shows the
analogous process for preparing excited states. In E, we show the process of
generating $g(\tfrac{1}{1+|z|}(H_{G}-zI))$, with $z$ defined below, whose
spectrum is shown on the vertical axis. This defines a model
$(\mu,\delta)$-problem that we use to show optimality. Although its ground
state is close to -1, note that the new gap is asymptotically smaller. Diagram
F shows the logical flow of the optimality proof for Proposition 2. The
original Grover Hamiltonian with gap $\delta=\frac{2}{\sqrt{N}}$ must take
$\Omega(1/\delta)$ time to solve. Using this Hamiltonian we show optimality of
our algorithm for all $0\leq y\leq 1$. Optimality of $y=0$ follows directly
from the top arrow as the algorithm solves the $y=0$ problem optimally (this
optimality was already known from previous works [14]). For all other $y$, we
map $H_{G}$ to $H^{\prime}$, with gap $\delta^{1+y/(2-y)}$, which when solved
by our algorithm finds the marked element of Grover’s in time $1/\delta$,
demonstrating optimality for all $y$.

##  4 Projecting to the Ground State

By using the Chebyshev approximation of $f(x)$ from Lemma 1 as our projector,
we can apply this speedup to the ground state preparation problem, identifying
$\mu$ with the (normalized) ground state energy and $\delta$ with the spectral
gap. Because QSP works best with polynomials of definite parity, we need to do
the projection in two steps. First, using the odd component of $p(x)$, we
generate a new “effective gap” in the spectrum of size
$\delta^{\prime}=1+O(\epsilon)$; this has query complexity
$\widetilde{O}(\frac{\sqrt{1+\mu+\delta/2}}{\delta})$. Second, we use the
method of [14] to project into the ground state; this has complexity
$\widetilde{O}(\frac{1}{\delta^{\prime}})=O(1)$. We compose these two quantum
circuits and thus multiply their complexities, for a final query complexity
that has a dependence on $\mu$ and $\delta$ of
$\widetilde{O}(\frac{\sqrt{1+\mu+\delta/2}}{\delta})$. This algorithm
furnishes a proof of Proposition 1 when we have $\mu=-1+\delta^{y}$ for some
$0\leq y\leq 1$.

###### Proof.

Proof of Proposition 1. Fix $H_{0}$ and its $(\alpha,m,0)$ block-encoding
$U_{1}$, and denote $p_{\text{odd}}=\frac{1}{2}(p(x)-p(-x))$ for the order-$n$
polynomial $p$ defined in Lemma 1. Given $U_{1}$, there is a $(1,m+1,0)$
block-encoding $U_{2}$ of $p_{\text{odd}}(H)$ that uses $n$ queries to $U_{1}$
and $U_{1}^{\dagger}$ as well as $O((m+1)n)$ other primitive quantum gates
[9]. By construction, $p(\mu)\leq-1+O(\epsilon)$ and
$|p(\lambda)|=O(\epsilon)$ for any eigenvalue with $\lambda>\mu$, so the
spectral gap of $p_{\text{odd}}(H)$ is now $\Omega(1)$. Thus, we may use the
algorithm of [14] to project into the ground state $\ket{\mu}$ using
$O(\frac{1}{\gamma}\log\frac{1}{\epsilon})$ queries to $U_{2}$, or

$\displaystyle
O\left(\frac{n}{\gamma}\log\frac{1}{\epsilon}\right)=O\left(\frac{1}{\gamma}\delta^{y/2-1}\log^{5/2}\frac{1}{\epsilon}\log\frac{1}{\delta}\right)=\widetilde{O}\left(\frac{1}{\gamma}\delta^{y/2-1}\right)$
(10)

queries to $U_{1}$ and $U_{1}^{\dagger}$, as desired. ∎

##  5 Optimality of the Algorithm

The algorithm of Proposition 1 achieves the optimal dependence of its query
complexity on $\delta$ and $\mu$, up to a factor of $\log\delta$, in the case
where the starting ansatz is good i.e. $\gamma=\Omega(1)$. This lower bound is
formalized by Proposition 2, and we prove it by adapting the optimality
argument presented in [14] and thus reducing unstructured search to the ground
state preparation problem. We will focus on a “Grover Hamiltonian” $H_{G}$
defined as follows: given a Hilbert space of dimension $N$, define the usual
marker oracle $U_{t}$ and Grover diffusion operator $D$ as

$\displaystyle U_{t}=I-2\ket{t}\bra{t}\quad\quad
D=I-2\ket{u}\bra{u}\quad\quad\ket{u}=\frac{1}{\sqrt{N}}\sum_{s}\ket{s}$ (11)

The Hamiltonian $H_{G}=\frac{1}{2}D+\frac{1}{2}U_{t}$ encodes in its ground
state $\ket{\mu}\propto\ket{u}+\ket{t}$ the solution $\ket{t}$ of the
corresponding unstructured search problem; the spectral gap is
$\delta=\frac{2}{\sqrt{N}}$. Thus, for any chosen (small) $\delta>0$, one can
choose a large enough $N$ such that $H_{G}$ has gap $O(\delta)$; this fixes
the ground state energy of $H_{G}$ as $\mu=-\delta/2$. We take the starting
ansatz to be $\ket{u}$, which gives us
$\gamma=\braket{u}{\mu}\approx\frac{1}{\sqrt{2}}+O(1/\sqrt{N})=\Omega(1)$.

Our algorithm from Proposition 1 works with ground state energies that are
near -1 and, hence, a ground state energy very close to 0 doesn’t achieve any
speedup relative to the generic case. One first naive attempt to prove this
lower bound is to simply shift the Hamiltonian so the ground state is closer
to -1. This will improve the bound over the generic case but doesn’t saturate
the upper bound of Proposition 1. Instead, we perform another transformation
of $H_{G}$ such that the algorithm of Proposition 1 on this transformed
Hamiltonian solves Grover’s in time $O(\sqrt{N})$; any algorithm faster than
one in Proposition 1 then would violate the known lower bound of Grover’s. In
the rest of this section, we describe in detail how this works.

It is useful to measure the ground state energy relative to the endpoint $-1$,
so we define the notation $\eta=O_{-1}(\mu)$ to mean that the number $\eta$ is
asymptotically closer to $-1$ than is $\mu$, or more precisely that
$\eta+1<c(\mu+1)$ in the $\delta\to 0$ limit for constant $c$. In this
section, we refer to preparing the ground state of a (normalized) Hamiltonian
$H$ with ground state energy $O_{-1}(\mu)$ and spectral gap $O(\delta)$ (we
will ignore small constant factors) as a $(\mu,\delta)$-problem.

The lower bound of ref.[14] for generic Hamiltonians can be represented in the
following way. For any $\delta$ there is a corresponding $(0,\delta)$-problem
that encodes the solution to the unstructured search, which has a query
complexity lower bound of $\Omega(1/\delta)$, thus proving that the generic
$(0,\delta)$-problem is $\Omega(1/\delta)$. We need to formulate a lower bound
for the general $(\mu,\delta)$-problem when $\mu$ is asymptotically close to
$-1$.

Given a Grover Hamiltonian $H_{G}$ with gap $\delta^{1-y}$ for some chosen
$0\leq y<1$, we can shift it to $H_{\textrm{shifted}}=H_{G}-zI$. This shift
affects both $\mu$ and $\delta$ because it modifies the block-encoding
normalization $\alpha$. Requiring that the $H_{\textrm{shifted}}$ has gap
$\delta$ fixes $x$ and results in a $(-1+\delta^{y},\delta)$-problem. Since
the ground state of $H_{\textrm{shifted}}$ can’t be prepared in time faster
then $\Omega(1/\delta^{1-y})$, the generic $(-1+\delta^{y},\delta)$-problem
has a lower bound of $\Omega(\delta^{y-1})$ queries when $y<1$.

On the other hand, the algorithm of Proposition 1 when $\mu=-1+\delta^{y}$
solves any $(-1+\delta^{y},\delta)$-problem with query complexity
$O(\delta^{y/2-1}\operatorname{polylog}\tfrac{1}{\delta})$, which can be quite
a bit larger than $\delta^{y-1}$ in the $\delta\to 0$ limit that we are
considering. We will need an additional step to raise the Grover-derived lower
bound and prove that the $y/2-1$ exponent of our algorithm is indeed optimal.

###### Proof.

Proof of Proposition 2. The extra ingredient required is another polynomial
transformation of the input $H_{G}$, denoted $g(x)$. We’ll use $g(x)$ to map
$\mu$ to approximately $-1$, and we’ll see that this mapping can be made to
shrink the gap $\delta$ less than the simple shift. This polynomial must have
two properties: it must have minimum value $\min_{x}g(x)=-1$ on the interval
$x\in[-1,1]$, and it must be one-to-one in a sufficiently large region around
the ground state. Specifically, we will use the explicit polynomial
$g(x)=ax+bx^{3}+cx^{5}$ with $a=19/6,b=-16/3,c=8/3$, which achieves its
minimum of $-1$ at $x=-\frac{1}{2}$. Now, we create our model
$(\mu^{\prime},\delta^{\prime})$-problem in two steps: first, we introduce the
$(-\tfrac{1}{2}+x_{0},\delta)$-problem $H$, which has the same spectrum as
$H_{G}$ but is shifted by a constant and scaled such that its ground state is
exactly $-\frac{1}{2}+x_{0}$ for some chosen $x_{0}>0$, with $|x_{0}|<1/2$:

$\displaystyle H=\frac{1}{1+|z|}\left(H_{G}-zI\right)\quad\quad
z=\frac{1-2x_{0}-\delta}{1+2x_{0}}$ (12)

Using the linear-combination-of-unitaries approach [17], $H$ can be block-
encoded with $\alpha=1$. Then, we map this $H$ through the polynomial $g$, and
$g(H)$ will define a particular $(-1+\delta^{\prime
y},\delta^{\prime})$-problem with an associated lower bound in terms of
$\delta^{\prime}$.

We choose the offset as $x_{0}=\delta^{\nu}$ for some $0<\nu\leq 1$; note that
$x_{0}$ is small in the limit $\delta\to 0$ that we are concerned with. (This
shift does not asymptotically affect the gap size.) Since
$g^{\prime}(-\tfrac{1}{2})=0$ and $g^{\prime\prime}(-\tfrac{1}{2})=O(1)$, we
have

$\displaystyle g(\mu)$ $\displaystyle=g(\tfrac{1}{2}+x_{0})\simeq-1+x_{0}^{2}$
(13) $\displaystyle g(\mu+\delta)$
$\displaystyle\simeq-1+x_{0}^{2}+2x_{0}\delta+\delta^{2}$ (14)

Since the new gap of $g(H)$ is $g(\mu+\delta)-g(\mu)$, we have used $g$ to map
the $(0,\delta)$-problem to a
$(-1+x_{0}^{2},2x_{0}\delta+\delta^{2})$-problem. Since $x_{0}=\delta^{\nu}$,
to leading order in $\delta$, we’ve reduced the Grover problem to a
$(-1+\delta^{2\nu},\delta^{\nu+1})$-problem, or, rescaling $\delta$ to
$\delta^{\prime}=\delta^{\nu+1}$, a
$(-1+\delta^{\prime\frac{2\nu}{\nu+1}},\delta^{\prime})=(-1+\delta^{\prime
y},\delta^{\prime})$-problem for any $0<y\leq 1$.

Because $g(H)$ and $H$ have the same ground state, and because querying $g(H)$
once requires $O(1)$ queries of $H$, the original $\Omega(1/\delta)$ bound for
preparing the ground state of $H$ also applies to $g(H)$. This implies that
the solution to the generic $(-1+\delta^{\prime y},\delta^{\prime})$-problem
for any $y\leq 1$ has query complexity bounded by

$\displaystyle\Omega(\tfrac{1}{\delta})=\Omega\left(\delta^{\prime\tfrac{-1}{1+\nu}}\right)=\Omega\left(\delta^{\prime\tfrac{y}{2}-1}\right)$
(15)

as desired.

∎

##  6 Preparing Excited States with Projectors

A straightforward generalization of the projection polynomial $p(x)$ from
Section 3 allows the preparation of excited states. This applies both to any
general Hamiltonian, extending the work of ref.[14], to excited states as well
as to nearly frustration free Hamiltonians where we again get an additional
speedup even for low-lying excited states. More precisely, given a normalized
excited-state energy $\eta$ (with, as before, $|\eta|\leq 1$), and a radius
$\delta_{1}$ which separates the excitation from both the ground state and
other excitations, we can form a new polynomial $p_{\eta}$ such that
$p_{\eta}(H)$ selects the subspace of eigenstates with energy
$\lambda\in(\eta-\delta_{1},\eta+\delta_{1})$. Moreover, the degree of
$p_{\eta}$ scales exactly the same way as that of the ground state projector
$p$, as $n=\widetilde{O}\left(\frac{\sqrt{1+\eta}}{\delta_{1}}\right)$.

######  Corollary 1.

If a Hamiltonian $H$, with block-encoding $U$, has a unique first excited
state $\ket{\eta}$ separated by an energy gap of at least $\delta_{1}$ from
the rest of the spectrum, it can be prepared to error $\epsilon$ with
$\widetilde{O}\left(\frac{\sqrt{1+\eta}}{\gamma\delta_{1}}\right)$ queries to
$U$ and $U^{\dagger}$, where $\gamma=|\braket{\phi_{0}}{\eta}|$ is the overlap
with the starting ansatz $\ket{\phi_{0}}$.

###### Proof.

As in Lemma 1, we use shifted error functions to accomplish this. Note that
the function
$f_{1}(x)=\tfrac{1}{2}(\operatorname{erf}k(x-\eta-\delta_{1}/2)-\operatorname{erf}k(x-\eta+\delta_{1}/2))$,
with $k=O(\frac{1}{\delta_{1}}\log^{1/2}\frac{1}{\epsilon})$, satisfies
$f_{1}(x)\leq-1+O(\epsilon)$ when $|x-\eta|\leq\delta_{1}/2$ and
$|f_{1}(x)|\leq O(\epsilon)$ when $|x-\eta|\geq\delta_{1}$. By the estimates
of Lemma 1, $f_{1}(x)$ can be approximated to error $\epsilon$ with a
degree-$n_{1}$ polynomial $p_{1}$, where

$\displaystyle
n_{1}=O\left(\frac{\sqrt{1+\eta}}{\delta_{1}}\log\frac{\sqrt{1+\eta}}{\delta_{1}}\log^{3/2}\frac{1}{\epsilon}\right)$
(16)

Again form the antisymmetric $p_{1,\text{odd}}=(p_{1}(x)-p_{1}(-x)$; as
before, the spectral gap of $p_{1,\text{odd}}(H)$ is now $\Omega(1)$, but the
ground state is now $\ket{\eta}$ with eigenvalue
$\lambda_{\eta}=-1+O(\epsilon)$. We can then use the algorithm of [14] to
project to $\ket{\eta}$ with the same query complexity as in Proposition 1, or

$\displaystyle
O\left(\frac{\sqrt{1+\eta}}{\delta_{1}\gamma}\log^{5/2}\frac{1}{\epsilon}\log\frac{1}{\delta}\right)=\widetilde{O}\left(\frac{\sqrt{1+\eta}}{\delta_{1}\gamma}\right)$
(17)

queries to $U$ and $U^{\dagger}$, as desired. ∎

Although we can form $p_{\eta}$ for any $|\eta|\leq 1$, as with ground state
preparation an additional speedup is possible if $\eta$ is close to the
endpoints $-1$ or $1$. This could be the case for the low-lying excited states
of a nearly-frustration-free $H_{0}$, for the following reason: Suppose the
Hamiltonian $H$ has a unique first excited state $\ket{\eta}$. Then, if the
ground state energy is $\mu=-1+\delta^{y}$ and $y\leq 1$, then we have
$\eta=-1+\delta^{y}+\delta=O_{-1}(\mu)$. This implies that the query
complexity of this first-excited-state preparation problem has gap dependence
$\widetilde{O}(\delta^{y/2}/\delta_{1})$, where
$\delta_{1}=\min(\delta,\delta^{\prime})$ and $\delta^{\prime}$ separates
$\ket{\eta}$ from the next excited state.

##  7 Example Hamiltonians with Speedup

In this section, we consider a class of Hamiltonians for which our algorithms
achieve a speedup over the standard approaches. In particular, we will look at
systems for which we have sub-volumetric perturbations to frustration-free
Hamiltonians. Note that for frustration free Hamiltonians, our algorithm gives
the expected $1/\sqrt{\Delta}$ dependence as one would achieve from uniform or
spectral gap amplification.

Sub-volumetric perturbations show up in a number of interesting physical
scenarios; we now enumerate some examples. To begin with, many physical
systems have defects (i.e. line or point defects) which intercolate amongst
the bulk degrees of freedom. For example, a frustration free $N\times N$
system perturbed by a density of $N$ impurities (i.e. approximately one per
row), or a constant number of such line-defects, gets an efficiency gain from
our algorithm. Note that there is a long history of interesting physics
related to impurity problems with relevance to both physics examples, such as
disordered solids [6], as well as algorithmic methods such as DMFT [8] and NRG
[18].

In addition to defect problems, there is interesting physics that often
happens at interfaces [12]. One can consider a system with two frustration-
free Hamiltonians which are connected at an interface. While the 2-d bulk of
both of these models will be frustration free, the 1-d interface will have
Hamiltonian terms that generically are frustrating. Our algorithm is more
efficient at finding ground states in this case. As a final example, quantum
computers are designed to “replace” classical simulation techniques for
finding ground states. A very common technique in numerical methods such as
DMRG is to introduce an additional pinning field [20, 1] to promote certain
phases or test robustness of phase diagrams. The introduction of this pinning
field on the boundary is sub-volumetric and hence has a scaling with $\Delta$
for which our algorithm applies.

Each of these scenarios follow from the decomposition of the full system
Hamiltonian into a frustration-free piece $H_{0}$ and a perturbation $V$. Here
we explicitly show that such Hamiltonians generically have the key property
that $\mu<-1+\delta^{y}$ for $y>0$. Consider a frustration-free $H_{0}$
defined on a 2D lattice of size $N\times N$ such as
$H_{0}=\sum_{j=1}^{N^{2}}\alpha_{j}U_{j}$, where the $U_{j}$ may for instance
be 2-body local interactions, and suppose for simplicity that $|\alpha_{j}|=1$
for all $j$. Without loss of generality, we may assume that $-1$ is an
eigenvalue of each $U_{j}$ by absorbing any phases into the coefficients
$\alpha_{j}$, so that the ground state is $\mu=-N^{2}$. Say that we add a 1-D
defect of $N$ perturbations
$V=\sum_{k=1}^{N}V_{k}=\sum_{k=1}^{N}vU^{\prime}_{k}$, where $v$ is the
magnitude of each perturbation. We also assume that the the $U^{\prime}_{k}$
are unitary and mutually commute with each other but do not commute with, and
hence frustrate, $H_{0}$. Deducing the $N$ dependence of this system, we have
$\alpha=N^{2}+vN=O(N^{2})$ when $v$ is held fixed. So long as the system is
gapped, so that $\Delta$ does not vanish in the large $N$ limit, we have
$\mu_{0}\leq-N^{2}+Nv$ or $\mu\leq-1+\frac{v}{N}$. Thus, because
$\delta=\frac{\Delta}{N^{2}}$, we have the situation described in Proposition
1 with $y=1/2$. The algorithm of that Proposition can then prepare the ground
state of $H_{0}+\sum_{k=1}^{N}V_{k}$ with a query complexity dependence on $N$
of

$\displaystyle O\left(\delta^{y/2-1}\right)=O(\delta^{-3/4})=O(N^{3/2})$ (18)

which is better than $O(N^{2})$ obtained without considering that
$H_{0}+\sum_{k=1}^{N}V_{k}$ is nearly frustration free. This argument
generalizes to $d$-dimensional systems where perturbations are added sub-
volumetrically.

## Acknowledgements

We acknowledge support from the Department of Energy grant DOE DESC0020165.
M.T. thanks Samuel Goldman for fruitful discussions.

## References

  * Assaad and Herbut [2013] Fakher F. Assaad and Igor F. Herbut.  Pinning the order: The nature of quantum criticality in the hubbard model on honeycomb lattice.  _Phys. Rev. X_ , 3:031010, Aug 2013.  doi: 10.1103/PhysRevX.3.031010.  URL https://link.aps.org/doi/10.1103/PhysRevX.3.031010.
  * Brassard et al. [2002] Gilles Brassard, Peter Høyer, Michele Mosca, and Alain Tapp.  Quantum amplitude amplification and estimation.  _Quantum Computation and Information_ , page 53–74, 2002.  ISSN 0271-4132.  doi: 10.1090/conm/305/05215.  URL http://dx.doi.org/10.1090/conm/305/05215.
  * Bravyi et al. [2008] Sergey Bravyi, David P. Divincenzo, Roberto Oliveira, and Barbara M. Terhal.  The complexity of stoquastic local Hamiltonian problems.  _Quantum Information and Computation_ , 8(5):0361–0385, 2008.  ISSN 15337146.  doi: 10.26421/qic8.5-1.
  * Chakraborty et al. [2019] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The power of block-encoded matrix powers: Improved regression techniques via faster Hamiltonian simulation.  In _Leibniz International Proceedings in Informatics, LIPIcs_ , volume 132, pages 1–58, 2019.  ISBN 9783959771092.  doi: 10.4230/LIPIcs.ICALP.2019.33.
  * Dong et al. [2021] Yulong Dong, Xiang Meng, K. Birgitta Whaley, and Lin Lin.  Efficient phase-factor evaluation in quantum signal processing.  _Physical Review A_ , 103(4):1–26, 2021.  ISSN 24699934.  doi: 10.1103/PhysRevA.103.042419.
  * Dyre and Schrøder [2000] Jeppe C Dyre and Thomas B Schrøder.  Universality of ac conduction in disordered solids.  _Reviews of Modern Physics_ , 72(3):873, 2000\.
  * Ge et al. [2019] Yimin Ge, Jordi Tura, and J. Ignacio Cirac.  Faster ground state preparation and high-precision ground energy estimation with fewer qubits.  _Journal of Mathematical Physics_ , 60(2):1–25, 2019.  ISSN 00222488.  doi: 10.1063/1.5027484.
  * Georges et al. [1996] Antoine Georges, Gabriel Kotliar, Werner Krauth, and Marcelo J. Rozenberg.  Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions.  _Rev. Mod. Phys._ , 68:13–125, Jan 1996.  doi: 10.1103/RevModPhys.68.13.  URL https://link.aps.org/doi/10.1103/RevModPhys.68.13.
  * Gilyén et al. [2019] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: Exponential improvements for quantum matrix arithmetics.  _Proceedings of the Annual ACM Symposium on Theory of Computing_ , pages 193–204, 2019.  ISSN 07378017.  doi: 10.1145/3313276.3316366.
  * Grover [1996] Lov K. Grover.  A fast quantum mechanical algorithm for database search.  In _Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing_ , STOC ’96, page 212–219, New York, NY, USA, 1996. Association for Computing Machinery.  ISBN 0897917855.  doi: 10.1145/237814.237866.  URL https://doi.org/10.1145/237814.237866.
  * Haah [2019] Jeongwan Haah.  Product decomposition of periodic functions in quantum signal processing.  _Quantum_ , 3(September):1–22, 2019.  ISSN 2521327X.  doi: 10.22331/q-2019-10-07-190.
  * Hellman et al. [2017] Frances Hellman, Axel Hoffmann, Yaroslav Tserkovnyak, Geoffrey S. D. Beach, Eric E. Fullerton, Chris Leighton, Allan H. MacDonald, Daniel C. Ralph, Dario A. Arena, Hermann A. Dürr, Peter Fischer, Julie Grollier, Joseph P. Heremans, Tomas Jungwirth, Alexey V. Kimel, Bert Koopmans, Ilya N. Krivorotov, Steven J. May, Amanda K. Petford-Long, James M. Rondinelli, Nitin Samarth, Ivan K. Schuller, Andrei N. Slavin, Mark D. Stiles, Oleg Tchernyshyov, André Thiaville, and Barry L. Zink.  Interface-induced phenomena in magnetism.  _Rev. Mod. Phys._ , 89:025006, Jun 2017.  doi: 10.1103/RevModPhys.89.025006.  URL https://link.aps.org/doi/10.1103/RevModPhys.89.025006.
  * Kitaev [1995] A. Yu. Kitaev.  Quantum measurements and the Abelian Stabilizer Problem.  pages 1–22, 1995.  URL http://arxiv.org/abs/quant-ph/9511026.
  * Lin and Tong [2020] Lin Lin and Yu Tong.  Near-optimal ground state preparation.  _Quantum_ , 4:1–22, 2020.  ISSN 2521327X.  doi: 10.22331/Q-2020-12-14-372.
  * Low and Chuang [2017a] Guang Hao Low and Isaac L. Chuang.  Hamiltonian simulation by uniform spectral amplification.  _arXiv_ , pages 1–32, 2017a.  ISSN 23318422.
  * Low and Chuang [2017b] Guang Hao Low and Isaac L. Chuang.  Optimal Hamiltonian Simulation by Quantum Signal Processing.  _Physical Review Letters_ , 118(1):1–6, 2017b.  ISSN 10797114.  doi: 10.1103/PhysRevLett.118.010501.
  * Low and Chuang [2019] Guang Hao Low and Isaac L. Chuang.  Hamiltonian Simulation by Qubitization.  _Quantum_ , 3:163, 2019.  ISSN 2521-327X.  doi: 10.22331/q-2019-07-12-163.
  * Pižorn and Verstraete [2012] Iztok Pižorn and Frank Verstraete.  Variational numerical renormalization group: Bridging the gap between NRG and density matrix renormalization group.  _Physical Review Letters_ , 108(6):1–8, 2012\.  ISSN 00319007.  doi: 10.1103/PhysRevLett.108.067202.
  * Somma and Boixo [2013] R. D. Somma and S. Boixo.  Spectral gap amplification.  _SIAM Journal on Computing_ , 42(2):593–610, 2013\.  ISSN 00975397.  doi: 10.1137/120871997.
  * Stoudenmire and White [2012] E. M. Stoudenmire and Steven R. White.  Studying two-dimensional systems with the density matrix renormalization group.  _Annual Review of Condensed Matter Physics_ , 3(1):111–128, 2012.  ISSN 19475454.  doi: 10.1146/annurev-conmatphys-020911-125018.


# Where to Deploy an Airborne Relay in Unknown Environments: Feasible
Locations for Throughput and LoS Enhancement

Juan David Pabon, Matthew C. Valenti, and Xi Yu
West Virginia University, Morgantown, WV, USA.

###### Abstract

The deployment of heterogeneous teams of both air and ground mobile assets
combines the advantages of mobility, sensing capability, and operational
duration when performing complex tasks. Air assets in such teams act to relay
information between ground assets but must maintain unblocked paths to enable
high-capacity communication modes. Obstacles in the operational environment
may block the line of sight (LoS) between air assets and ground assets
depending on their locations and heights. In this paper, we analyze the
probability of spanning a two-hop communication between a pair of ground
assets deployed in an environment with obstacles at random locations and with
random heights (i.e. a Poisson Forest) using an air asset at any location near
the ground assets. We provide a closed-form expression of the LoS probability
based on the 3-dimensional locations of the air asset. We then compute a 3-D
manifold of the air asset locations that satisfy a given LoS probability
constraint. We further consider throughput as a measure of communication
quality, and use it as an optimization objective.

##  I Introduction

The use of coordinated autonomous mobile agents, consisting of both ground and
air assets, has become prevalent in various applications such as long-term
monitoring and post-disaster rescue in large and complex operational
environments like urban areas and forests [1, 2, 3, 4, 5, 6]. The deployment
of such a heterogeneous set of mobile agents offers several advantages. The
air assets provide swift mobility and have large sensing areas, while the
ground assets offer high sensing accuracy and long operational duration [7].
In different scenarios, air assets may act as communication relays [8] or may
provide remote sensing of the battleground [9]. The locations of these air
assets should be carefully planned such that they maintain connectivity with
the proper sets of ground assets.

Streaming a large volume of data between assets, especially between air assets
and ground assets, is essential for their successful coordination [10]. Modern
communication technologies tend to boost the data transmission rate by using
very short wavelengths (e.g. mmWave frequencies or even visible light), which
have compromised penetration capability and are, therefore, easily blocked
[11]. In real-world operational environments, it is common to encounter
obstacles of various locations, shapes, and heights that can block the line of
sight (LoS) between different assets. Such LoS blockage may also hinder
localization and mapping techniques that rely on cameras or lidar [12, 13].

A benefit of using mobile agents is that they can reposition themselves to
recover the necessary LoS paths. Assuming that a detailed map of the obstacles
is provided, mobile assets can carefully plan their routes through the
environment to maintain LoS [14]. However, such an assumption is usually a
luxury for various reasons. Civilian activities, military operations, and
disasters may change the landscape of a previously known area. New obstacles
can emerge, existing structures may be destroyed or altered, and the terrain
can be transformed. Certain environments, like forests, can pose additional
difficulties. They may contain a large number of obstacles, such as trees and
dense vegetation, which can be challenging to catalog and represent accurately
on a map. Relying solely on pre-existing maps for path planning and obstacle
avoidance is imprudent in light of these factors.

An alternative approach to address the limitations of fixed, deterministic
maps is to model obstacles in urban areas or forests as distributed randomly.
For instance, the authors of [15, 16, 17] model the locations of buildings in
an urban area using a Manhattan Poisson line process with randomly distributed
building heights. Forested areas were modeled as having randomly located trees
in [18, 19, 20] and with random heights in [21, 22, 23]. A forest-like
cluttered environment can be modeled via a Poisson point process (PPP) and is
referred to as a Poisson forest.

In environments where the locations of the obstacles and their heights can be
modeled as a stochastic process, the probability that the path between two
assets crosses an obstacle, thereby blocking the LoS, can be computed
accordingly. The so-called _LoS probability_ is then the complement of the
probability that the path is blocked. Our previous work [24] demonstrated the
benefits of deploying an air asset as a communication relay for ground assets
by deriving the LoS probability and throughput in a Poisson forest. However,
in that work, we did not solve the problem of where to deploy the air assets,
but instead assumed that the air asset was located precisely halfway between
the two ground assets.

Deploying an air asset precisely in the middle of a pair of ground assets
presents challenges, especially when dealing with moving ground assets or a
constantly moving air asset such as a fixed-wing aircraft. Consequently, this
paper focuses on the issue of identifying optimal locations for the air asset
_away_ from the midpoint. Our criterion for finding these locations is based
on achieving a constant LoS probability. The nexus of these locations form a
curve in 2-D space or a surface in 3-D space, defining a volume within which
the LoS probability can be guaranteed. By constraining the air asset’s
location or path to remain within this volume, the desired performance target
is maintained.

In addition to considering locations of constant LoS probability, our study
also takes into account _capacity_ and _throughput_ metrics. More precisely,
the capacity is the _Shannon_ capacity of the LoS link, while the throughput
is defined as the expected value of capacity with respect to the LoS
probability. These metrics compensate for distance and discourage infeasible
solutions that suggest placing the air asset at very high altitudes. Similar
to LoS probability, the metrics of capacity and throughput can be utilized to
define a surface in 3-D space where the values remain constant. As will be
shown, the surface of constant throughput is highly nonlinear yet relatively
easy to compute using the theory developed in this paper. Consequently, the
contributions presented in this paper significantly simplify the challenging
task of planning the locations of air assets. The proposed relay placement
strategies in this paper offer a moderate level of complexity, which can
contribute to improved real-time decision-making and response capabilities.

The paper is organized as follows. Sec. II formulates the problem. Sec. III
provides an analysis that provides closed-form solutions to determining the
positions of an air asset that result in a constant value of LoS probability.
Sec. IV continues the analysis by considering capacity and throughput,
similarly finding positions for the air asset that result in constant values
of these quantities. Sec. V provides numerical results that illustrate the
analysis and allow general trends to be noted. Finally, Sec. VI concludes the
paper.

##  II Problem Formulation

Consider a pair of ground assets equipped with communication devices whose
heights are equal to $h_{g}$. These assets are deployed in a operational space
($e.g.$ a forest) with randomly distributed obstacles ($e.g.$ trees). The
location of these obstacles can be described by a two-dimensional PPP with a
fixed density $\lambda_{f}$. This density represents the expected number of
_potential_ obstacles per unit area111In this paper, we distinguish between
_potential_ obstacles and _blockages_. A potential obstacle is an obstacle of
height $H$ that may or may not block the LoS path. On the other hand, a
blockage is a potential obstacle that is tall enough to block the path.
Whether the LoS is blocked or not depends on the heights of the obstacle and
the communicating devices involved.. Let $N$ be the number of potential
obstacles in an operational environment of area $A_{f}$, which is referred to
as a Poison forest. From the properties of a PPP, the probability that the
number of potential obstacles in the area is equal to $n$ is

$\mathbb{P}\\{N=n\\}=\frac{({\lambda_{f}A_{f})}^{n}}{n!}e^{-\lambda_{f}A_{f}}.$

In a Poison forest, the height of any obstacle is represented by a non-
negative random variable $H$. The distribution of $H$ may vary [21, 22, 23].
Let $F_{H}(h)$ denote the cumulative distribution function (CDF) of $H$. While
the analysis in this paper is not limited to any particular distribution, to
illustrate our method and provide closed-form expressions, this paper focuses
on the specific realistic cases where $H$ assumes a uniform distribution or a
truncated Gaussian distribution. For the uniform distribution,
$H\in[0,h_{\max}]$, where $h_{\max}$ is the maximum height that an obstacle
could have. The truncated Gaussian is created by taking a Gaussian random
variable with mean $\mu$ and standard deviation $\sigma$, and then
conditioning the variable on $h\geq 0$.

Let $(x,y)$ denote a location on the ground plane. The ground assets are
assumed to be distance $g$ apart and placed at locations $(0,0)$ and $(g,0)$
on the ground plane as shown in Fig. 1. Any potential obstacle with a height
above $h_{g}$ located between the ground assets blocks the LoS between the
assets, and hence, becomes a blockage. Consider that the obstacles have a non-
trivial thickness $W$ and that the average thickness of the obstacles is
$\mathbb{E}(W)$. In the Poisson forest, the distribution of potential
obstacles along the straight line joining any two points can be characterized
by a 1-D Poisson process with fixed density
$\lambda_{0}=\mathbb{E}(W)\lambda_{f}$ describing the expected number of
obstacles per unit length.

Figure 1: Top-view showing the location of the air and ground assets in the
ground plane. The ground assets are located at coordinates $(0,0)$ and
$(g,0)$, and the air asset is located at an altitude of $h_{a}$ directly above
the coordinate $(x_{a},y_{a})$. The values $r_{a}$ and $r_{b}$ represent the
horizontal distances from the ground assets to the air asset.

When the ground assets are spaced far apart, it is very likely that at least
one obstacle will block the LoS of the direct path between them. To mitigate
this issue, an air asset can be deployed as a relay to aid in communication.
Relaying follows a two-hop protocol, with the first hop involving a
transmission from the first ground asset to the air asset, and the second hop
involving a transmission from the air asset to the second ground asset. Let
the altitude (or height) of the air asset be $h_{a}$ and its location, as
projected onto the ground plane, be $(x_{a},y_{a})$ as shown in Fig. 1. Define
the _horizontal distance_ between a ground asset and the air asset to be the
distance between the ground asset and the projection of the air asset onto the
ground plane. Let $r_{a}$ be the horizontal distance between the ground asset
at $(0,0)$ and the air asset, and $r_{b}$ be the horizontal distance between
the ground asset at $(g,0)$ and the air asset, as shown in Fig. 1. It follows
that $r_{a}=\sqrt{x_{a}^{2}+y_{a}^{2}}$ and
$r_{b}=\sqrt{(g-x_{a})^{2}+y_{a}^{2}}$.

In our previous work [24], the LoS probability is computed considering a fixed
position of the air asset (in the middle of the ground assets). Maintaining
the air asset at a fixed position can introduce vulnerabilities, as it becomes
more susceptible to attacks. Additionally, certain air assets (_e.g._ fixed-
wing aircraft or blimps) may lack precise control mechanisms to remain in the
same position consistently. Furthermore, the ability to explore a large
terrain or facilitate communication among multiple assets can be compromised
if the air asset is required to remain in a static position at all times.

When the air asset has the flexibility to be deployed throughout the entire
operational environment, it becomes essential to identify the positions that
can deliver the desired performance and fulfill the communication requirements
of the ground assets. The performance can be defined in terms of Los
probability and throughput, which is a metric that trades off between
communication distance and the effect of LoS blockage. In the following
sections, we calculate regions in 3-D space identifying where the air relay
will provide desired values for LoS probability, capacity, and throughput.
Furthermore, the altitude of the air asset is optimized by maximizing the
throughput at each possible position in the horizontal plane.

##  III Positions with Constant LoS probability

In [24], the LoS probability for the hop between one of the ground assets,
$i\in\\{a,b\\}$, and the air asset222Note that the LoS probability is
reciprocal; i.e., the LoS probability from a ground asset to an air asset is
the same as the LoS probability from the air asset to the same ground asset.
was found to be:

$\mathbb{P}_{LoS}^{ga}(r_{i})=\exp\left(-\int_{0}^{r_{i}}\lambda(r)dr\right).$
(1)

$\lambda(r)$ is the density of a PPP that represents blockages; _i.e._ , only
those potential obstacles that exceed a _critical height_ $h_{c}$ and
therefore may block the LoS path. $h_{c}$ is the minimum height of an obstacle
able to block the ground-air LoS at a given distance $r$ from a ground asset.
$h_{c}$ is distance dependent, increasing with the distance $r$. Because
$h_{c}$ is distance dependent, so is the density $\lambda(r)$, and thus the
PPP is inhomogeneous.

More specifically, $\lambda(r)$ depends on the density $\lambda_{0}$
representing the locations of the potential obstacles along with the
probability distribution of the obstacles’ heights evaluated at the critical
height $h_{c}(\cdot)$. In particular, the density $\lambda(r)$ is determined
as follows (see [24] for details):

$\lambda(r)=\lambda_{0}\left[1-F_{H}\left(h_{c}\left(r\right)\right)\right].$
(2)

When the value of $\mathbb{P}^{ga}_{LoS}(r_{i})$ is fixed, $i.e.$ equal to a
constant within the interval $(0,1]$, (1) can be solved for $r_{i}$. The
solution provides the horizontal distance between the ground and the air asset
that results in the fixed value of $\mathbb{P}^{ga}_{LoS}(r_{i})$. Once
solved, the locations of constant LoS probability for a single hop form a
circle around the ground asset. This circle has a radius of $r_{i}$, and by
restricting the movement of the air asset to fly along this circular path, a
constant LoS probability can be maintained.

However, in this paper, we are concerned with using the air asset to provide a
relay connecting two ground assets. For the two ground assets to be connected,
we assume that there must be a clear LoS from the air asset to _both_ of the
ground assets. Define $\mathbb{P}_{LoS}^{gag}(r_{a},r_{b})$ to be the
probability of LoS for two-hop ($i.e.$, ground-air-ground) communication,
which requires a LoS from the first ground asset to the air asset and a LoS
from the air asset to the second ground asset. Assuming that the two paths are
independent (a reasonable assumption unless their angular distance is close,
see [25]), this probability is the product of the LoS probabilities of each of
the two hops, $i.e.$,

$\mathbb{P}_{LoS}^{gag}(r_{a},r_{b})=\mathbb{P}_{LoS}^{ga}(r_{a})\mathbb{P}_{LoS}^{ga}(r_{b}).$
(3)

As with the case for one ground asset, it is possible to fix the value of LoS
probability in (3) and find the location of the air asset that provides a
constant value. However, because there are now two ground assets, the region
of constant LoS probability will become an ellipse rather than a circle, as
will be shown subsequently for height distributions that are uniform or
truncated Gaussian.

###  III-A LoS for Uniform Distribution

When $H$ is uniform over $[0,h_{\max}]$, the end-to-end LoS probability can be
expressed as

$\mathbb{P}_{LoS}^{gag}(r_{a},r_{b})=e^{-\Lambda_{u}(r_{a})}e^{-\Lambda_{u}(r_{b})}$
(4)

where $\Lambda_{u}(\cdot)$ is the result of the integral in (1). Assuming that
the air asset is deployed at $h_{a}>h_{\max}$, $\Lambda_{u}(\cdot)$ is equal
to

$\Lambda_{u}(r_{i})=\lambda_{0}r_{ci}\left[1-\dfrac{h_{g}}{h_{\max}}-\left(\dfrac{h_{a}-h_{g}}{2r_{i}h_{\max}}\right)r_{ci}\right],$

where $r_{ci}$ for $i\in\\{a,b\\}$ is the critical distance beyond which the
critical height is taller than $h_{\max}$. It is found to be:

$r_{ci}=r_{i}\left(\dfrac{h_{\max}-h_{g}}{h_{a}-h_{g}}\right).$

Fixing the value of $\mathbb{P}_{LoS}^{gag}(r_{a},r_{b})=P$ and rewriting (4),
the following equation is obtained:

$r_{a}+r_{b}=c_{u},$ (5)

where $c_{u}$ is

found to be:

$c_{u}=\dfrac{-(h_{a}-h_{g})\ln(P)}{\lambda_{0}(h_{\max}-h_{g})\left(1-\dfrac{h_{g}}{h_{\max}}-\dfrac{h_{\max}-h_{g}}{2h_{\max}}\right)}.$

Using (5), the relation between the distance from the air asset to each ground
asset and the desired value of $P$ can be obtained. For a fixed $h_{a}$, (5)
describes an ellipse since the sum of the distances from the air asset to the
ground assets $a$ and $b$ is constant, _i.e._ , the value $c_{u}$.

In (5), $r_{a}$ and $r_{b}$ can be expressed in terms of the position of the
air asset in the ground plane. Then, the relation between the location of the
air asset and the desired value of $P$ can be obtained. In Fig. 1, it is
observed that the coordinate $(x_{a},y_{a})$ is the location of the air asset
in the ground plane. By replacing $(r_{a},r_{b})$ with the corresponding
values of $(x_{a},y_{a})$, (5) can be written as

$\sqrt{x_{a}^{2}+y_{a}^{2}}+\sqrt{(g-x_{a})^{2}+y_{a}^{2}}=c_{u}.$ (6)

After simplifying (6) and substituting $(x,y)=(x_{a},x_{b})$, the following
equation is obtained

$\dfrac{4}{c_{u}^{2}}\left(x-\dfrac{g}{2}\right)^{2}+\dfrac{4}{c_{u}^{2}-g^{2}}y^{2}=1.$
(7)

When the air asset flies above coordinate $(x,y)$, it will have the desired
two-hop LoS probability $P$. When $(c_{u}^{2}-g^{2})>0$, the locations of
$(x,y)$ allowed by (7) form an ellipse with a major axis along the x-axis and
foci located at $(0,0)$ and $(g,0)$; _i.e.,_ at the positions of the ground
assets.

###  III-B LoS for Truncated Gaussian Distribution

Now consider the case that $H$ is a truncated Gaussian. To find the positions
of the airborne relay producing a constant two-hop LoS probability, start with
(1). For the truncated Gaussian distribution, the integral of (1) is equal to
$r_{i}\left(c+ck\right)$, where $c=\lambda_{0}/(2\Phi(\mu/\sigma))$ and
$\Phi(\cdot)$ is the CDF of the standard normal distribution. The value of $k$
is given by

$k=\left(\dfrac{e^{-a^{2}}+\sqrt{\pi}\left[(d-a)\text{erf}(a-d)+a\text{erf}(a)\right]-e^{-(a-d)^{2}}}{\sqrt{\pi}d}\right)$

where erf($\cdot$) is the error function, $a=(\mu-h_{g})/(\sqrt{2}\sigma$) and
$d=(h_{a}-h_{g})/(\sqrt{2}\sigma)$. Substituting the result of this
integration into (1), the following equation is obtained:

$\mathbb{P}_{LoS}^{ga}(r_{i})=e^{-r_{i}(c+ck)}.$

It follows from (3) that the two-hop LoS probability is equal to

$\mathbb{P}_{LoS}^{gag}(r_{a},r_{b})=e^{-(r_{a}+r_{b})(c+ck)}.$

Fixing the desired value of $\mathbb{P}_{LoS}^{gag}=P$ gives

$P=e^{-(r_{a}+r_{b})(c+ck)}.$

Rewriting this equation gives

$r_{a}+r_{b}=c_{t}$ (8)

where $c_{t}=-\ln(P)/(c+ck)$. As with (5), equation (8) describes an ellipse,
only now the two distances sum to a different value; _i.e.,_ they sum to
$c_{t}$.

By expressing $r_{a}$ and $r_{b}$ in terms of the position of the air asset in
the ground plane $(x,y)=(x_{a},y_{a})$, and assuming $g<c_{t}$, simplifying
(8) yields

$\dfrac{4}{c_{t}^{2}}\left(x-\dfrac{g}{2}\right)^{2}+\dfrac{4}{c_{t}^{2}-g^{2}}y^{2}=1.$
(9)

As with the uniform distribution, air assets flying above coordinate $(x,y)$
will provide the desired two-hop LoS probability $P$, and for fixed $h_{a}$,
the locus of all $(x,y)$ forms an ellipse.

##  IV Link Capacity and Expected Throughput

Sec. III found the positions of the air asset that produce a fixed two-hop LoS
probability. However, the LoS probability does not consider the signal power
loss that occurs as the distance between the air asset and the ground asset
increases. To provide a more comprehensive assessment of communication
performance, this section introduces the computation of capacity and
throughput metrics. These metrics take into account the distance-dependent
signal power loss, allowing for a more suitable objective when determining the
optimal placement of the airborne relay.

###  IV-A Positions with Constant Link Capacity

As a measure of capacity, we specifically consider the _Shannon_ capacity. The
Shannon capacity represents the maximum achievable data rate for an unblocked
link. It quantifies the upper limit of information transmission rate over the
channel, taking into account such factors as the channel propagation model,
distance transmitted, transmitted power, transmission bandwidth, and antenna
gains. By considering the Shannon capacity, we can assess the maximum data
rate that can be achieved under ideal conditions for a given communication
link, such as between the air asset and a ground asset.

The (Shannon) capacity of the link between ground asset $i\in\\{a,b\\}$ and
the air asset can be written as

$C_{i}=B\log_{2}(1+\mathsf{SNR}),$ (10)

where $B$ is the signal bandwidth and $\mathsf{SNR}$ is the signal-to-noise
ratio. When $B$ is in units of Hertz, $C_{i}$ is in units of bits-per-second
(bps). When expressed in dB, the value of $\mathsf{SNR}$ is given by

$\mathsf{SNR}^{\text{dB}}=\mathsf{SNR}^{\text{dB}}_{0}-10\alpha\log_{10}\left(\dfrac{d_{i}}{d_{0}}\right),$
(11)

where $d_{i}$ is the Euclidean distance from the air asset to ground asset
$i$, $\alpha$ is the path-loss exponent, and $d_{0}$ is the reference
distance. $\mathsf{SNR}^{\text{dB}}_{0}$ is the signal-to-noise ratio when the
transmission distance is equal to $d_{0}$ assuming free-space propagation up
to that distance. Note that, while it depends on the communication distance,
the capacity does not depend on the LoS probability.

For the ground asset located at the origin ($i=a$), the Euclidean distance
$d_{a}^{2}=x^{2}+y^{2}+h_{a}^{2}$, where $(x,y,h_{a})$ is the 3-D position of
the air asset. Rearranging (10) yields:

$x^{2}+y^{2}+h_{a}^{2}=d^{2}_{0}\left(\dfrac{10^{\mathsf{SNR}_{0}^{\text{dB}}/10}}{2^{C_{i}/B}-1}\right)^{2/\alpha}.$
(12)

For the other ground asset ($i=b$), which is located at coordinate $(g,0)$,
the term $x$ in (12) must be replaced by $g-x$ per Fig. 1.

As with the LoS probability, (12) shows that the airborne relay should be
located on a circle around the ground asset when there is just a single hop.
However, when two-hop (ground-air-ground) relaying is considered, the end-to-
end capacity is limited by the minimum of the capacities of the two hops.
Moreover the duplexing operation at the relay must be taken into account.
Thus, the two-hop capacity is $C=\min(C_{a},C_{b})/2$, where the
multiplication by $1/2$ accounts for the time-division duplexing operation at
the relay.

It follows that capacity is determined by the shorter of the two links. Hence,
when the air asset is closer to the ground asset located at the origin,
$d_{a}<d_{b}$ and $C_{a}>C_{b}$. Conversely, when $d_{b}<d_{a}$, then
$C_{b}>C_{a}$. Since $d_{a}<d_{b}$ when $x<g/2$, the capacity for the ground-
air-ground link is given by

$C=\frac{1}{2}\min(C_{a},C_{b})=\begin{cases}C_{b}/2&\text{if }x\leq g/2\\\
C_{a}/2&\text{if }x>g/2.\end{cases}$ (13)

###  IV-B Positions with Constant Throughput

The LoS probability serves as a valuable indicator for predicting the presence
or absence of a link between the assets, while the capacity provides an
estimation of the maximum achievable data rate if a link exists. However,
neither of these parameters alone provides a comprehensive assessment of the
overall communication quality. The LoS probability accounts for blockages but
does not consider the impact of transmission distance on signal power loss. On
the other hand, capacity considers the transmission distance but overlooks the
blockage process. An appropriate metric for assessing the quality of the link
is the _throughput_ , which strikes a balance between LoS probability and
capacity [24].

The capacity provided in (10) is the capacity of a given link when it is
unblocked. When the link is blocked, its capacity is zero. Thus, in
environments with random blocking, the capacity of link $i$ is a random
variable that assumes a value of $C_{i}$ with probability
$\mathbb{P}_{LoS}^{ga}$ and a value of zero with probability
$1-\mathbb{P}_{LoS}^{ga}$. The throughput of a link in the presence of
blocking is the expected value of the capacity of the link, where the
expectation is with respect to the LoS probability. In particular, the
throughput of a single hop is

$\displaystyle T_{i}$ $\displaystyle=$
$\displaystyle\mathbb{P}_{LoS}^{ga}C_{i}.$ (14)

For a two-hop communication, the throughput is the expectation of the end-to-
end capacity, where the expectation is with respect to the two-hop LoS
probability; _i.e.,_ it is as follows:

$T=\begin{cases}\mathbb{P}_{LoS}^{gag}C_{b}/2&\text{if }x\leq g/2.\\\
\mathbb{P}_{LoS}^{gag}C_{a}/2&\text{if }x>g/2.\end{cases}$ (15)

Fixing the value of throughput in (15) and allowing the 3-D position
$(x,y,h_{a})$ of the air asset to vary, we can find a surface in 3-D space
that guarantees the desired throughput. Also, the positions of the air asset
that maximize the throughput can be found, for instance, by fixing the
coordinate on the ground plane $(x,y)$ and determining the altitude $h_{a}$
that maximizes the throughput.

##  V Numerical Results

This section presents numerical results to illustrate the theoretical concepts
discussed in Sections III and IV. The target values of LoS probability,
capacity, and throughput are varied to determine the 2-D and 3-D surfaces for
which these values are constant. This provides insight into possible flight
paths for the air asset that provide the necessary performance metrics.

Unless otherwise specified, it is assumed that the values of the key physical
parameters are assumed to be $\lambda_{0}=0.02$, $h_{g}=2$ m, $B=100$ MHz,
$d_{0}=1$ m, $\mathsf{SNR}_{0}^{\text{dB}}=50$ dB, and $\alpha=2.3$. This
path-loss coefficient corresponds to measured LoS pathloss at $38$ GHz [26].
The value of $\mathsf{SNR}_{0}^{\text{dB}}$ corresponds to a carrier frequency
of 38 GHz, a bandwidth of 100 MHz, a transmit power of $0$ dBm, a receiver
noise figure of $11$ dB, and antenna gains of $12.1$ dBi for both the transmit
and receive antennas, which are the gains reported for a compact 6-element
array operating at $38$ GHz in [27]. Both kinds of obstacle height
distributions are considered; for the case that the heights are uniform we use
$h_{\max}=29$ m and for the case that they are truncated Gaussian we use
$\mu=19$ m and $\sigma=10$ m.

###  V-A LoS Probability

This section shows the results obtained when the desired LoS probability $P$
and the distance between ground assets $g$ take different values in (7) and
(9).

Fig. 2 considers the case that $h_{a}=100$ m and $g=60$ m. In this figure, red
dots indicate the location of the two ground assets, while the black ellipses
show the locations for the air asset that provide constant LoS probabilities
equal to $P=0.8$ (inner ellipse), $P=0.65$ (middle ellipse), and $P=0.5$
(outer ellipse). The subplot on the left corresponds to the case that the
obstacle height distribution is uniform while the subplot on the right
corresponds to the case that it is truncated Gaussian. For greater LoS
probability the eccentricity of the ellipses increases and for smaller
probabilities the eccentricity decreases and the major axis of the ellipse
increases its length.

Figure 2: Air asset positions producing fixed LoS probability. Obstacle
heights are (a) uniform distribution and (b) truncated Gaussian distribution.
Red dots indicate the location of the ground assets.

When $h_{a}$ changes and the same LoS probability is required, the elliptic
cone presented in Fig. 3 is obtained. This surface represents the positions
that produce the desired LoS probability $P$. In this case, $P=0.7$ and the
uniform obstacle height distribution is considered. The minimum altitude
$h_{a}$ at which the air asset can be deployed to produce the desired $P$ with
a given distance $g$ between ground assets can be determined by equating
$c_{u}$ to $g$ and solving for $h_{a}$. This minimum altitude determines the
distance from the bottom of the elliptic cone to ground. At the bottom of the
cone $h_{a}=44.29$ m.

Figure 3: Air asset positions producing $P=0.7$. The contour plot of the
surface is shown in the x-y plane. Red dots indicate the location of the
ground assets. Uniform height distribution is considered.

Fig. 4 shows the positions of the air asset that guarantee a fixed LoS
probability for three different values of distance between ground assets
($g=\\{10,30,50\\}$ m) when $h_{a}=100$ m. In this case, the obstacle height
distribution is truncated Gaussian. When the distance between the ground
assets is small, the ellipse has little eccentricity (i.e., it is almost a
circle), while the eccentricity grows with increasing $g$. As $g$ tends to
$c_{t}$ (or, for uniform height distributions, $h_{u}$), the positions for the
air asset will tend towards a straight line connecting the ground assets. That
distance $g$ will be the maximum distance that the assets could separate while
the air asset is able to guarantee the desired LoS probability.

Figure 4: Ellipses representing air asset positions producing $P=0.8$. The
black dot represents the position of ground asset $a$. Dots with different
colors represent positions of ground asset $b$ and each dot has an ellipse
with the same color. Truncated Gaussian height distribution is considered.

###  V-B Capacity and Throughput

In Fig. 5, the value of capacity in (13) is fixed at $C=100$ Mbps, and the 3-D
region of constant capacity is shown. When the air asset’s position is such
that $x\leq g/2$, the capacity of the ground-air-ground link is limited by the
capacity of the link between the air asset and ground asset $b$, and the
surface is obtained evaluating (12) at $r_{b}$. Similarly, for $x>g/2$, the
capacity is determined by the link between ground asset $a$ and the air asset,
and the surface is obtained evaluating (12) at $r_{a}$. Any position inside
the volume covered by this surface will produce a capacity greater than 100
Mbps. Because capacity does not account for the presence of obstacles, the
region does not depend on the obstacle height distribution.

Figure 5: Air asset’s positions producing $C=100$ Mbps. A contour plot is
projected onto the xy-plane.

Next, we consider regions of constant throughput, as it is a metric that
balances capacity with LoS probability. Fig. 6 shows a 3-D surface
representing the positions of the air asset that guarantee a throughput equal
to 80 Mbps when the uniform height distribution is considered. The following
observations can be made about the surface that is shown. As shown in 3, it
can be observed that as the altitude $h_{a}$ of the air asset increases, the
contours of constant LoS probability expand, indicating a larger area covered
by LoS connections. On the other hand, in Figure 5, as $h_{a}$ increases, the
contours of constant capacity shrink. Since the throughput is the product of
LoS probability and capacity, in Fig. 6, both behaviors are observed, with the
cross section areas initially increasing with $h_{a}$ as the regions of
constant LoS probability expand. But then, after a certain value of $h_{a}\geq
53.6$ m, the area of the cross sections decreases as the constant-capacity
contours contract with increasing $h_{a}$. The volume of the region contained
by the surface is inversely proportional to the throughput; i.e., if a smaller
throughput were considered, then the region shown would be larger.

Figure 6: (a) air asset’s positions producing $T=80$ Mbps. (b) contour plot of
the surface shown in (a). The obstacle heights are uniformly distributed.
Figure 7: Air asset positions producing the maximum throughput. The color of
the surface represents the value of throughput indicated by the color bar. The
contours in the xy-plane represent the contours for the same air asset
altitude. The obstacle heights are uniformly distributed.

In addition to identifying regions of constant throughput, it is also possible
to optimize equation (15) with respect to the altitude $h_{a}$ of the air
asset. This optimization process determines the air asset altitude that
maximizes the throughput for a given position in the ground plane. Fig. 7
shows the air asset altitude that maximizes throughput for each position of
the air asset over the ground plane. It is observed that the maximum possible
throughput is obtained when the air asset is located in $(g/2,0,36.3)$ for
$g=60$ m. The surface shown in Fig. 7 allows us to determine the positions
across which the air asset should move if it is required to obtain the maximum
possible throughput for any of the positions. Additionally, the contours for
different heights of the surface in Fig. 7 are shown in the xy-plane.

##  VI Conclusions and future work

In this paper, we have developed a framework for determining feasible
locations to deploy an airborne relay to enhance communication between two
ground assets while satisfying performance metrics including LoS probability,
capacity, and throughput. As a metric, throughput strikes a balance between
deploying the air asset at locations that increases LoS probability while
simultaneously decreases signal power loss due to distance.

A key contribution of this paper is an approach for calculating the throughput
based on the location of the air asset. This methodology takes into account
the heights of obstacles and their locations described by an inhomogeneous
PPP. Closed-form expressions are derived for determining the regions where the
air asset can be deployed to achieve desired values of LoS probability or
capacity. Knowledge of the maximum achievable throughput facilitates the
deployment of the air asset in locations that optimize communication
performance, and can be used to identify the optimal altitude to deploy the
air asset.

Future research directions include investigating optimal trajectory planning
for the air asset when the ground assets are in motion. Additionally,
optimizing trajectory planning to deploy the air asset in scenarios involving
multiple ground assets requiring communication would be a valuable area of
exploration.

## References

  * [1] J. Langelaan and S. Rock, “Towards autonomous UAV flight in forests,” in _Proc. AIAA Guidance, Navigation, and Control Conference and Exhibit_ , 2005, p. 5870.
  * [2] X. Zheng, S. Koenig, D. Kempe, and S. Jain, “Multirobot forest coverage for weighted and unweighted terrain,” _IEEE Transactions on Robotics_ , vol. 26, no. 6, pp. 1018–1031, 2010.
  * [3] K. Harikumar, J. Senthilnath, and S. Sundaram, “Multi-UAV oxyrrhis marina-inspired search and dynamic formation control for forest firefighting,” _IEEE Transactions on Automation Science and Engineering_ , vol. 16, no. 2, pp. 863–873, 2018.
  * [4] M. S. Couceiro, D. Portugal, J. F. Ferreira, and R. P. Rocha, “Semfire: Towards a new generation of forestry maintenance multi-robot systems,” in _Proc. IEEE/SICE International Symposium on System Integration (SII)_ , 2019, pp. 270–276.
  * [5] Y. Tian, K. Liu, K. Ok, L. Tran, D. Allen, N. Roy, and J. P. How, “Search and rescue under the forest canopy using multiple UAVs,” _The International Journal of Robotics Research_ , vol. 39, no. 10-11, pp. 1201–1221, 2020.
  * [6] L. F. Oliveira, A. P. Moreira, and M. F. Silva, “Advances in forest robotics: A state-of-the-art survey,” _Robotics_ , vol. 10, no. 2, p. 53, 2021.
  * [7] B. Grocholsky, J. Keller, V. Kumar, and G. Pappas, “Cooperative air and ground surveillance,” _IEEE Robotics & Automation Magazine_, vol. 13, no. 3, pp. 16–25, 2006.
  * [8] L. Chaimowicz, A. Cowley, D. Gomez-Ibanez, B. Grocholsky, M. Hsieh, H. Hsu, J. Keller, V. Kumar, R. Swaminathan, and C. Taylor, “Deploying air-ground multi-robot teams in urban environments,” in _Multi-Robot Systems. From Swarms to Intelligent Automata Volume III: Proceedings from the 2005 International Workshop on Multi-Robot Systems_. Springer, 2005, pp. 223–234.
  * [9] M. Wei and V. Isler, “Air to ground collaboration for energy-efficient path planning for ground robots,” in _Proc. IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , 2019, pp. 1949–1954.
  * [10] R. Doriya, S. Mishra, and S. Gupta, “A brief survey and analysis of multi-robot communication and coordination,” in _Proc. International Conference on Computing, Communication & Automation_, 2015, pp. 1014–1021.
  * [11] E. Hriba, M. C. Valenti, K. Venugopal, and R. W. Heath, “Accurately accounting for random blockage in device-to-device mmwave networks,” in _Proc. IEEE Global Communications Conference_ , 2017, pp. 1–6.
  * [12] M. F. Mysorewala, D. O. Popa, and F. L. Lewis, “Multi-scale adaptive sampling with mobile agents for mapping of forest fires,” _Journal of Intelligent and Robotic Systems_ , vol. 54, no. 4, pp. 535–565, 2009.
  * [13] B. Benjamin, G. Erinc, and S. Carpin, “Real-time wifi localization of heterogeneous robot teams using an online random forest,” _Autonomous Robots_ , vol. 39, no. 2, pp. 155–167, 2015.
  * [14] A. Gasparetto, P. Boscariol, A. Lanzutti, and R. Vidoni, “Path planning and trajectory planning algorithms: A general overview,” _Motion and Operation Planning of Robotic Systems_ , pp. 3–27, 2015.
  * [15] F. Baccelli and X. Zhang, “A correlated shadowing model for urban wireless networks,” in _Proc. IEEE Conference on Computer Communications (INFOCOM)_ , 2015, pp. 801–809.
  * [16] M. Gapeyenko, D. Moltchanov, S. Dmitri, and R. W. Heath Jr, “Line-of-sight probability for mmwave-based UAV communications in 3D urban grid deployments,” _IEEE Transactions On Wireless Communications_ , vol. 20, no. 10, pp. 6566–6579, 2021.
  * [17] E. Hriba, M. C. Valenti, and R. W. Heath, “Optimization of a millimeter-wave uav-to-ground network in urban deployments,” in _Proc. IEEE Military Communications Conference (MILCOM)_ , 2021, pp. 861–867.
  * [18] S. Karaman and E. Frazzoli, “High-speed flight in an ergodic forest,” in _Proc. IEEE International Conference on Robotics and Automation_ , 2012, pp. 2899–2906.
  * [19] ——, “High-speed motion with limited sensing range in a Poisson forest,” in _Proc. 51st IEEE Conference on Decision and Control (CDC)_ , 2012, pp. 3735–3740.
  * [20] B. Martinez R and G. A. Pereira, “Fast path computation using lattices in the sensor-space for forest navigation,” in _Proc. IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 1117–1123.
  * [21] T. Kohyama and T. Hara, “Frequency distribution of tree growth rate in natural forest stands,” _Annals of Botany_ , vol. 64, no. 1, pp. 47–57, 1989.
  * [22] J. M. Felfili, “Diameter and height distributions in a gallery forest tree community and some of its main species in central Brazil over a six-year period (1985-1991),” _Brazilian Journal of Botany_ , vol. 20, no. 2, pp. 155–162, 1997.
  * [23] F. Mauro, R. Valbuena, J. Manzanera, and A. García-Abril, “Influence of global navigation satellite system errors in positioning inventory plots for tree-height distribution studies,” _Canadian Journal of Forest Research_ , vol. 41, no. 1, pp. 11–23, 2011.
  * [24] J. D. Pabon, S. Alkandari, M. C. Valenti, and X. Yu, “Air-aided communication between ground assets in a Poisson forest,” in _Proc. IEEE Military Communications Conference (MILCOM)_ , 2022, pp. 133–139.
  * [25] E. Hriba and M. C. Valenti, “The impact of correlated blocking on millimeter-wave personal networks,” in _Proc. IEEE Military Communications Conference (MILCOM)_ , 2018, pp. 1–6.
  * [26] T. S. Rappaport, S. Sun, R. Mayzus, H. Zhao, Y. Azar, K. Wang, G. N. Wong, J. K. Schulz, M. Samimi, and F. Gutierrez, “Millimeter wave mobile communications for 5g cellular: It will work!” _IEEE Access_ , vol. 1, pp. 335–349, 2013.
  * [27] Y. Rahayu and M. I. Hidayat, “Design of 28/38 GHz dual-band triangular-shaped slot microstrip antenna array for 5G applications,” in _Proc. International Conference on Telematics and Future Generation Networks (TAFGEN)_ , 2018, pp. 93–97.