# Energy Efficiency Analysis of Charging Pads-powered UAV-enabled Wireless
Networks

Yujie Qin, Mustafa A. Kishk, Member, IEEE, and Mohamed-Slim Alouini, Fellow,
IEEE Yujie Qin, Mustafa A. Kishk, and Mohamed-Slim Alouini are with Computer,
Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King
Abdullah University of Science and Technology (KAUST), Thuwal, 23955-6900,
Saudi Arabia (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>slim.alouini@kaust.edu.sa).

###### Abstract

This paper analyzes the energy efficiency of a novel system model where
unmanned aerial vehicles (UAVs) are used to provide coverage for user hotspots
(user clusters) and are deployed on charging pads to enhance the flight time.
We introduce a new notion of ”cluster pairs” to capture the dynamic nature of
the users spatial distribution in order to exploit one of the top advantages
of UAVs, which is the mobility and relocation flexibility. Using tools from
stochastic geometry, we first derive a new distance distribution that is vital
for the energy efficiency analysis. Next, we compute the coverage probability
under two deployment strategies: (i) one UAV per cluster pair, and (ii) one
UAV per cluster. Finally, we compute the energy efficiency for both
strategies. Our numerical results reveal which of the two strategies is better
for different system parameters. Our work investigates some new aspects of the
UAV-enabled communication system such as the dynamic density of users and the
advantages or disadvantages of one- or two-UAV deployment strategies per
cluster pair. By considering the relationships between the densities of user
cluster pairs and the charging pads, it is shown that an optimal cluster pair
density exists to maximize the energy efficiency.

###### Index Terms:

Stochastic geometry, Poisson Point Process, Poisson Cluster Process, Unmanned
Aerial Vehicles, Charging Pads.

##  I Introduction

As a result of rapid technological progress, unmanned aerial vehicles (UAVs,
also known as drones) play an increasingly essential role in enhancing the
performance of current and future wireless networks [1, 2, 3]. The primary
benefits of UAV-enabled networks include the improvement of the coverage
probability by establishing line-of-sight (LoS) links with ground users [4],
[5], deployment in dangerous environments or in natural disasters, such as
severe snow storms or earthquakes, and more services available to remote
Internet of Things (IoT) users in rural areas [6].

In addition, unlike the fixed terrestrial base stations (TBSs), UAVs can
easily satisfy the dynamic traffic demands by functioning as aerial BSs, owing
to their flexibility, mobility, and 3D-deployment ability [7]. At places where
the spatial distributions of active users continuously changing over time,
UAVs are more suitable and eligible as they are flexible enough to optimize
their locations in real time. Under the circumstance in which the users’
locations exhibit certain degree of spatial clustering and dynamic variants,
UAVs can help offload TBSs by offering an additional dynamic capacity [8]. For
instance, in some scenarios where the densities of active users vary with time
rapidly (e.g., commercial area and office area), UAVs will travel to the user
clusters with higher user densities and, therefore, supply the users with
stable connectivity, making UAVs, a feasible and practical alternative.

Generally, one of the most popular tools for analyzing cellular networks is
stochastic geometry [9], and the most widely used cluster model is the Poisson
cluster process (PCP) [10], such as Matern cluster process (MCP) and Thomas
cluster process (TCP) [11]. However, there is always an underlying assumption
of homogeneity, which assumes the user density keeps a constant value and does
not capture the dynamic variants of the system. Generally, analysis becomes
more complex when considering the dynamic variants of the system due to the
correlation on time and locations. This paper develops a novel framework that
takes the dynamic changes of user densities in a day into account and studies
the system’s dynamic performance based on stochastic geometry.

###  I-A Related Work

Literature has many studies related to this work which can be categorized
into: (i) optimization of energy efficiency of a cellular network, (ii) energy
efficiency of UAV-enabled wireless network, and (iii) stochastic geometry-
based analysis of UAV-enabled wireless network. For each of these areas, a
brief introduction is discussed in the following lines.

Optimization of energy efficiency of a cellular network. Authors in [12]
comprehensively surveyed recent advances in system architectures and key
techniques related to the impact of the constrained fronthaul on spectral
efficiency, energy efficiency, and open issues on software-defined networking.
Based on stochastic geometry and the deployment of sleeping strategies,
authors in [13] studied the associated trade-offs of energy-efficient cellular
networks. They derived the success probability and formulated the minimization
of power consumption problem and maximization of energy efficiency problem.
Authors in [14] analyzed the minimization of the network energy cost by
optimizing the BS density for both homogeneous and heterogeneous cellular
networks. The results revealed which type of BSs should be deployed or slept
first when the traffic load is low. In [15], the authors focused on how small
cell access points can be used to offload the traffic and how they can enhance
energy efficiency. They also evaluated a distributed sleep-mode strategy and
the trade-off between traffic offloading and energy consumption based on
stochastic geometry.

Energy efficiency of UAV-enabled wireless network. A comprehensive survey of
UAV-enabled networks was given in [16], which included in-depth discussions on
UAV-based communication systems, classification based on application,
altitude, and network, and summarized some open issues and future research
directions. Authors in [17] studied the energy efficiency of UAV-assisted
network with wireless charging. While UAVs aimed to serve multiple users,
energy harvesting time, transmission power, mission complete time and UAV
trajectories were jointly considered to minimize the total energy consumption.
Authors in [18] and [19] considered a single fixed-wing UAV behaving as a
relay between two ground terminals and maximized the system throughput by
optimizing the transmit power and trajectory of the UAV. In [20], the authors
proposed a novel method to increase the spectral efficiency of UAV small cell
by autonomous reposition according to the mobility and activities of users. A
self-organized approach with full three-dimensional spatial flexibility was
studied in [21] to provide the highest possible average spectral energy
efficiency in a disaster scenario, subject to interference from the remaining
UAVs. In [22] a non-stationary fixed-wing UAV-based relay was analyzed, the
authors proposed to maximize the energy efficiency by optimizing the circular
maneuvering and communication scheme. In [23], BS sleeping strategy was
proposed to maximize the energy efficiency of a two-tier downlink cellular
network.

Stochastic geometry-based analysis of UAV-enabled wireless network. A tutorial
on the stochastic geometry-based analysis of cellular networks was provided in
[24] and [25], which briefly introduced point processes, characterized the
interference, and exact error performance. In addition, they presented the
outage and ergodic rate analysis and showed that stochastic geometry is
capable of capturing realistic network performance. A heterogeneous network
composed of terrestrial and aerial BSs was considered in [5] and [26] where
the locations were modeled as two independent Poisson point processes (PPPs).
For this setup, the authors first derived accurate expressions of coverage
probability and average data rate and then obtained the approximated equations
by using the upper bound of the Gamma distribution. Author in [27, 28, 29]
studied the availability probability and revised coverage probability of UAVs
subject to the limited capacity of charging stations, where the locations of
UAVs and charging stations were modeled as two independent PPPs, and users
were modeled as MCP. In [30, 31, 8], authors studied a system with tethered
UAV, which denotes a UAV connected to a ground station to get a stable power
supply and a reliable data link but with restricted mobility.

While the existing literature focus on the spatial distribution of UAVs and
TBSs, and some strategies such as sleep-mode to optimize energy efficiency,
there is no work to capture the influence of the dynamic variation of user
spatial density on UAV deployment and system energy efficiency.

###  I-B Contribution

This paper investigates the energy efficiency of a UAV-enabled wireless
network while capturing the dynamic changes in the user spatial densities. In
addition, to avoid the flight time limitation of UAVs, we consider a scenario
in which the UAV deploys itself on top of a carefully chosen charging pad. To
achieve that objective, we introduce user cluster pairs where the user density
in the clusters in each pair vary with time. Modeling the locations of the
user cluster pairs as a Poisson bipolar network and the locations of the UAV
charging pads as an independent PPP, we investigate the performance of two UAV
deployment approaches. The contributions of this paper can be summarized as
follows:

Novel Framework and Performance Metrics: Since the user densities are
different according to the time of a day, we introduce a novel framework, the
user cluster pair, which consists of two user clusters with varying densities
of the user. In this way, we are able to analyze the influence of the dynamic
user densities on the system performance. We first define the switching
frequency as the number of changes of the user density per day. We then
compare the energy efficiency under different switching frequencies and two
types of UAV deployment, one UAV deployed for each cluster pair and one UAV
deployed for each user cluster, respectively.

Evaluation of Traveling Distance: Using tools from stochastic geometry, we
derive the expression for the traveling distance, which is defined as the
distance between two charging pads that are the nearest to two user clusters
in a user cluster pair. The resulting expression and its approximation match
well with the simulation and are readily used in the performance analyses.

Energy Efficiency: Given the aforementioned expression of traveling distance,
we are able to study the influence of dynamic user density on the system
performance. We first derive the expressions for coverage probability as
functions of traveling distance. Next, we obtain the energy efficiency of the
above two scenarios, respectively, under the approximation of traveling
distance distribution.

System-Level Insights: By comparing the energy efficiency of the
aforementioned two scenarios, our numerical results reveal various useful
system insights. We show that user densities do have a significate impact on
system performance. For instance, deploying two UAVs is not always better than
only one UAV deployed in some scenarios, owing to the interference and
relationships between the densities of users. More details are provided in
Section VI.

##  II System Model

In wireless networks, the locations of users tend to be clustered and the
density of users in a cluster varies with time. For example, users are
typically concentrated in an office area in the morning or commercial area at
night. These two user clusters exist at the same time, but the user densities
vary with the time of a day. For instance, the user density of the office
cluster is higher than the commercial cluster in the morning and lower at
night. Such two user clusters behave like cluster pairs, while the density of
users of one cluster is higher than the other one but their relationship
changes with time. In our system, we assume that each user cluster pair
consists of two user clusters, say $x_{m}$ and $x_{n}$, and they are located
at a fixed distance $d_{\rm nm}$ in a random orientation. In this case, the
user densities of $x_{m}$ and $x_{n}$ form the density pairs ($\lambda_{\rm
m,h}$,$\lambda_{\rm n,l}$) and ($\lambda_{\rm m,l}$,$\lambda_{\rm n,h}$),
where the subscript $h$ and $l$ denote the higher and lower densities of
users. That is, $\lambda_{\rm m,h}>\lambda_{\rm n,l}$ and $\lambda_{\rm
n,h}>\lambda_{\rm m,l}$, and we are interested in analyzing the influence of
switching frequency on the system performance. Here, switching frequency
denotes the number of changing of density pairs between ($\lambda_{\rm
m,h}$,$\lambda_{\rm n,l}$) and ($\lambda_{\rm m,l}$,$\lambda_{\rm n,h}$) per
day.

Different from the previous work [27], we use charging pads in this work to
prolong the service time of UAVs. By doing so, UAVs can provide service
without being interrupted. Assume that the charging pads are located on the
top of buildings at a fixed altitude $h$. Particularly, we consider that these
buildings are able to connect to the network via fiber cables, hence, no need
for UAVs to wirelessly backhaul. UAVs are deployed on charging pads to get a
stable power supply and provide service to users in user cluster to offload
TBSs. Note that deploying UAVs on the charging pads is different from roof-top
BSs, which is flexible, has low maintenance and installation cost and high
energy efficiency.

We consider a dynamic multi-access cellular network that consists of the
aforementioned user cluster pairs, UAVs, charging pads, and TBSs, where the
switching frequency of user cluster pairs is $n_{t}$. The locations of
charging pads and TBSs are modeled as two independent PPPs, $\Phi_{\rm c}$ and
$\Phi_{\rm t}$, with densities $\lambda_{\rm c}$ and $\lambda_{\rm t}$,
respectively. The locations of user cluster pairs are modeled as a Poisson
bipolar network model $\Phi_{\rm user}$ with density $\lambda_{\rm user}$, as
shown in Fig. 1. The locations of UAVs are conditioned on user cluster pairs
and denoted by $\Phi_{\rm u}$. The locations of users are uniformly
distributed in user clusters with radii $r_{c}$, modeled as MCP, and the
densities of users are $\lambda_{\rm\\{m,n\\},\\{h,l\\}}$ as defined. The
locations and densities of UAVs will be explained in details in the following
part of this paper.

Figure 1: System model. TABLE I: Table of Notations

Notation | Description
---|---
$\Phi_{\rm c}$, $\Phi_{\rm t}$; $\lambda_{\rm c}$, $\lambda_{\rm t}$ | PPP of charging pads, PPP of TBSs; density of the charging pads, density of TBSs
$x_{m}$, $x_{n}$, $d_{\rm nm}$ | The reference user cluster pair, distance between centers of $x_{m}$ and $x_{n}$
$\lambda_{\rm m,\\{h,l\\}}$, $\lambda_{\rm n,\\{h,l\\}}$ | User densities of $x_{m}$ and $x_{n}$ at different times, respectively
$\Phi_{\rm user}$,$\Phi_{u}$ | Poisson bipolar network of user cluster pairs, set of locations of UAVs
$\lambda_{\rm user}$, $\lambda_{\rm u,1}$, $\lambda_{\rm u,2}$ | Densities of user cluster pairs, UAVs of the frist and the second scenarios, respectively
$C_{m}$, $C_{n}$; $L$ | The nearest charging pads to the user clusters $x_{m}$ and $x_{n}$, respectively; distance between $C_{m}$ and $C_{n}$
$R_{\rm mm}$, $R_{\rm mn}$; $\theta_{1}$, $\theta_{3}$ | Distances between the $C_{m}$ and the cluster center of $x_{m}$ and $x_{n}$, respectively; related angles
$R_{\rm um}$, $R_{\rm un}$ | Distances between $C_{m}$ and the reference user in $x_{m}$ and $x_{n}$, respectively
$R_{\rm t}$, $R_{\rm nn}$, $\theta_{2}$ | Distance between the reference user and its nearest TBS, $C_{n}$ and cluster center of $x_{n}$, and related angle
$n_{t}$, $b_{w}$, $\gamma$ | Switching frequency of the cluster UAV, bandwidth, SINR threshold
$p_{\rm s}$, $p_{\rm m}$; $p_{\rm tba}$ | Service-, travel-related power consumption of UAVs; total power consumption of TBSs
$\mathcal{A}_{\rm LoS}$, $\mathcal{A}_{\rm NLoS}$, $\mathcal{A}_{\rm TBS}$ | Probability of associating with LoS/NLoS UAV, event of associating with the nearest TBS
${\rm EE}$, ${\rm SE}$, $P_{\rm tot}$ | Energy efficiency, spectrum efficiency, total energy consumption

We compare the energy efficiency of UAV deployments in two scenarios: one UAV
for each cluster pair, or two UAVs for each cluster pair. Without loss of
generality, we perform our analysis at a typical user cluster pair in which
the center of $x_{m}$ is located at the origin and the center of $x_{n}$ is
located at $(d_{\rm nm},0)$. Let $C_{m}$ and $C_{n}$ denote the nearest
charging pads to the user cluster $x_{m}$ and $x_{n}$, respectively.

######  Definition 1 (The First Scenario).

In the case of the first scenario, each user cluster pair has one UAV deployed
at the charging pad which is the nearest to the user cluster that has a higher
user density. In other words, UAV locates at $C_{m}$ when user density pair is
$(\lambda_{\rm m,h},\lambda_{\rm n,l})$ and $C_{n}$ when $(\lambda_{\rm
m,l},\lambda_{\rm n,h})$, as shown in Fig. 2. Therefore, the density of UAVs
is equal to the density of user cluster pairs.

######  Definition 2 (The Second Scenario).

In the case of the second scenario, each user cluster pair has two UAVs
deployed, one for $C_{m}$ and one for $C_{n}$ if they are different.
Otherwise, only one UAV deployed. That is, we only deployed one UAV when the
same charging pad is the nearest to both user clusters.

######  Remark 1.

In the frist scenario, the density of UAVs $\lambda_{\rm u,s_{1}}$ equals to
the density of user cluster pairs $\lambda_{\rm user}$ for sure. However, in
the second scenario, the density of UAVs is a function of $d_{\rm nm}$ and
$\lambda_{\rm c}$. When $\lambda_{\rm c}$ is high and $d_{\rm nm}$ is large,
$\lambda_{\rm u,s_{2}}$ has a high probability of being equivalent to
$2\lambda_{\rm user}$. This is because there will be a very low probability of
having the same charging pad being closest to both clusters.

Figure 2: System model of the first scenario (left) and the second scenario
(right).

###  II-A Distance Distribution

As mentioned in the system model, we assume that the cluster UAV locates at
the charging pad which is the nearest to the user cluster that has a higher
user density in the first scenario. It implies that the cluster UAV needs to
change its location if the density pair of the user clusters switches (e.g.,
UAV needs to travel from $C_{m}$ to $C_{n}$, if the density pair changes from
($\lambda_{\rm m,h}$,$\lambda_{\rm n,l}$) to ($\lambda_{\rm
m,l}$,$\lambda_{\rm n,h}$)). While we consider that UAVs cannot provide
service when traveling, the traveling time is negligible compared with the
time spent on charging pads.

######  Definition 3 (Traveling Distance).

Let $L$ be the distance between two charging pads $C_{m}$ and $C_{n}$, which
is conditioned on $R_{\rm mm}$ and $\theta_{1}$, the distance between $C_{m}$
and $x_{m}$, and its relative angle. To simplify the complexity of the
analysis, we assume that $R_{\rm mm}$ is shorter than $d_{nm}$.

Let $\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$ denote the probability that
$L$ equal to zero in the first scenario. Note that $\mathbb{P}_{\rm(0|R_{\rm
mm},\theta_{1})}$ is also the probability that only one UAV deployed in the
second scenario.

Therefore, the densities of UAVs in the above two scenarios are, respectively,
given by

$\displaystyle\lambda_{\rm u,1}$ $\displaystyle=\lambda_{\rm user},$
$\displaystyle\lambda_{\rm u,2}$ $\displaystyle=(2-\mathbb{E}_{\rm
R_{mm},\theta_{1}}[\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}])\lambda_{\rm
user},$ (1)

where $\lambda_{\rm user}$ is the density of user clusters, $\lambda_{\rm
u,1}$ and $\lambda_{\rm u,2}$ are the densities of UAVs in the first and the
sceond scenario, respectively.

###  II-B User Association

We assume that users associate with the cluster UAV or the nearest TBS that
provides the strongest average received power. Note that the access priority
in this work is limited for each UAV to only its own cluster users and this is
commonly known as closed access [32]. Without loss of generality, we focus on
the reference user randomly selected from the user cluster pair $x_{m}$ and
$x_{n}$. Let $q_{m,\\{h,l\\}}$ and $q_{n,\\{l,h\\}}$ be the probabilities of
the reference user being located in the user clusters $x_{m}$ and $x_{n}$ at
different time of a day, respectively,

$\displaystyle q_{m,\\{h,l\\}}$ $\displaystyle=\frac{\lambda_{\rm
m,\\{h,l\\}}}{\lambda_{\rm m,\\{h,l\\}}+\lambda_{\rm n,\\{l,h\\}}},$
$\displaystyle q_{n,\\{l,h\\}}$ $\displaystyle=\frac{\lambda_{\rm
n,\\{l,h\\}}}{\lambda_{\rm m,\\{h,l\\}}+\lambda_{\rm n,\\{l,h\\}}}.$ (2)

When the reference user associates with the cluster UAV, the received power is
given by,

$\displaystyle p_{\rm u}=\left\\{\begin{aligned} p_{\rm l}&=\eta_{\rm
l}\rho_{\rm u}G_{\rm l}D_{\rm\\{um,un\\}}^{-\alpha_{\rm l}},\text{in the case
of LoS},\\\ p_{\rm n}&=\eta_{\rm n}\rho_{\rm u}G_{\rm
n}D_{\rm\\{um,un\\}}^{-\alpha_{\rm n}},\text{in the case of NLoS},\\\
\end{aligned}\right.$ (3)

in which $\rho_{\rm u}$ is the transmit power of UAVs, $\alpha_{\rm l}$ and
$\alpha_{\rm n}$ are path-loss exponent, $G_{\rm l}$ and $G_{\rm n}$ present
the fading gains that follow gamma distribution with shape and scale
parameters ($m_{\rm l},\frac{1}{m_{\rm l}}$) and ($m_{\rm n},\frac{1}{m_{\rm
n}}$), $\eta_{\rm l}$ and $\eta_{\rm n}$ denote the mean additional losses for
LoS and NLoS transmissions, respectively,
$D_{\rm\\{um,un\\}}=\sqrt{h^{2}+R_{\rm\\{um,un\\}}^{2}}$, where
$R_{\rm\\{um,un\\}}$ are the horizontal between the reference user in $x_{m}$
or $x_{n}$ and the cluster UAV, respectively. The occurrence probability of a
LoS link established by the reference user and the cluster UAV is given in
[33] as

$\displaystyle P_{\rm l}(r)$
$\displaystyle=\frac{1}{1+a\exp(-b(\frac{180}{\pi}\arctan(\frac{h}{r})-a))},$
(4)

where $a$ and $b$ are two environment variables (e.g., urban, dense urban, and
highrise urban), and $h$ is the altitude of the UAV. Consequently, the
probability of NLoS link is $P_{\rm n}(r)=1-P_{\rm l}(r)$.

When the user associates with the nearest TBS, the received power is

$\displaystyle p_{\rm t}$ $\displaystyle=\rho_{\rm t}HR_{\rm t}^{-\alpha_{\rm
t}},$ (5)

where $H$ is the fading gain that follows exponential distribution with
average power of unity, $\rho_{\rm t}$ is the transmit power of TBSs and
$R_{\rm t}$ is the distance between the user and the nearest TBS.

We then define the association probability with UAVs, and an association event
that denotes the reference user being served by the nearest TBS.

######  Definition 4 (Association Probability and Event).

Conditioned on the serving base station is at distance $r$ away, let
$\mathcal{A}_{\rm NLoS}(r)$ and $\mathcal{A}_{\rm LoS}(r)$ be the
probabilities that the reference user associates with the cluster N/LoS UAV,
and let $\mathcal{A}_{\rm TBS}(r)$ be the event that the user associates with
the nearest TBS,

$\displaystyle\mathcal{A}_{\rm NLoS}(r)=\mathbb{P}(p_{n}(r)>p_{t}),$
$\displaystyle\mathcal{A}_{\rm LoS}(r)=\mathbb{P}(p_{l}(r)>p_{t}),$
$\displaystyle\mathcal{A}_{\rm TBS}(r)=\mathbbm{1}(p_{t}(r)>p_{u}).$ (6)

The reason of defining $\mathcal{A}_{\rm TBS}(r)$ as an association event will
be explained in detail in Remark 3.

If the signal-to-interference-plus-noise ratio (SINR) of the established
channel between the reference user and its serving base station is beyond the
designated threshold, this user is said to be successfully served. We define
the probability of successfully serving as coverage probability.

######  Definition 5 (Coverage Probability).

In the first scenario, the coverage probability $P_{\rm cov,s_{1}}$ is given
by

$\displaystyle P_{\rm cov,s_{1}}=$ $\displaystyle\alpha P_{\rm
cov,s_{1}|C_{m}}+(1-\alpha)P_{\rm cov,s_{1}|C_{n}}$ $\displaystyle=$
$\displaystyle\alpha(q_{m,h}P_{\rm cov,s_{1}|C_{m},q_{m,h}}+q_{n,l}P_{\rm
cov,s_{1}|C_{m},q_{n,l}})+(1-\alpha)(q_{m,l}P_{\rm
cov,s_{1}|C_{n},q_{m,l}}+q_{n,h}P_{\rm cov,s_{1}|C_{n},q_{n,h}}),$ (7)

where $\alpha$ is ratio of the time of cluster UAV spent at $C_{m}$ and the
total time of a day. $P_{\rm cov,s_{1}|c_{\\{m,n\\}},q_{m,\\{h,l\\}}}$ and
$P_{\rm cov,s_{1}|c_{\\{m,n\\}},q_{n,\\{l,h\\}}}$ are the coverage
probabilities of the reference user in $x_{m}$ and $x_{n}$ when UAV locates at
$C_{m}$ and $C_{n}$ are, respectively, given by

$\displaystyle P_{\rm
cov,s_{1}|C_{\\{m,n\\},qm,\\{h,l\\}}}=\sum_{bs}\mathbb{E}\left[\mathcal{A}_{bs}(r)\mathbb{P}(\operatorname{SINR}\geq\gamma\mid
r,bs)\right],\quad bs\in{\rm\\{N/LoS,TBS\\}},$ (8)

where $\gamma$ is the SINR threshold.

In the case of the second scenario, the coverage probability $P_{\rm
cov,s_{2}}$ is given by

$\displaystyle P_{\rm cov,s_{2}}$ $\displaystyle=P_{\rm cov,s_{2},U}+P_{\rm
cov,s_{2},T},$ (9)

where $P_{\rm cov,s_{2},U}$ and $P_{\rm cov,s_{2},T}$ are the coverage
probabilities of the cluster UAVs and the nearest TBS,

$\displaystyle P_{\rm cov,s_{2},U}$ $\displaystyle=\sum_{\rm
N/LoS}\mathbb{E}[\mathcal{A}_{\rm N/LoS}(r)\mathbb{P}({\rm
SINR}(r)\geq\gamma\mid r,{\rm N/LoS})],$ $\displaystyle P_{\rm cov,s_{2},T}$
$\displaystyle=\mathbb{E}[\mathcal{A}_{\rm TBS}(r)\mathbb{P}({\rm
SINR}(r)\geq\gamma\mid r,{\rm T})].$ (10)

Let $\Phi_{\rm u_{l}}$ and $\Phi_{\rm u_{n}}$ be subsets of $\Phi_{\rm u}$,
which denote the locations of LoS UAVs and NLoS UAVs, respectively.
Conditioned on the serving BS $b_{s}$, the SINR and aggregate interference
power are defined as

$\displaystyle{\rm SINR}=$
$\displaystyle\frac{\max(p_{u},p_{t})}{I+\sigma^{2}},$ $\displaystyle I=$
$\displaystyle\sum_{n_{i}\in\Phi_{\rm u_{n}}/b_{s}}\eta_{\rm n}\rho_{\rm
u}G_{\rm n_{i}}D_{\rm n_{i}}^{-\alpha_{\rm n}}+\sum_{l_{j}\in\Phi_{\rm
u_{l}}/b_{s}}\eta_{\rm l}\rho_{\rm u}G_{\rm l_{j}}D_{\rm l_{j}}^{-\alpha_{\rm
l}}+\sum_{t_{k}\in\Phi_{\rm t}/b_{s}}\rho_{\rm t}H_{\rm t_{k}}D_{\rm
t_{k}}^{-\alpha_{\rm t}},$ (11)

in which $\sigma^{2}$ is the additive white Gaussian noise (AWGN) power,
$D_{\rm n_{i}}$, $D_{\rm l_{j}}$ and $D_{\rm t_{k}}$ are the distances between
the reference user and the interfering NLoS, LoS UAVs, and TBSs, respectively.

###  II-C Energy Efficiency

We consider the UAV’s power consumption composed of travel-related power
$p_{\rm m}$ and service-related power $p_{\rm s}$. Based on [34], $p_{\rm s}$
is a constant and ${p}_{\rm m}$ is

$\displaystyle p_{\rm m}=p_{\rm 0}\left(1+\frac{3v^{2}}{u_{\rm
tip}^{2}}\right)+\frac{p_{\rm i}v_{\rm 0}}{v}+\frac{1}{2}d_{\rm 0}\rho
sa_{1}v^{3},$ (12)

where $p_{\rm 0}$ and $p_{\rm i}$ present the blade profile power and induced
power, $U_{\rm tip}$ is the tip speed of the rotor blade, $v_{\rm 0}$ is the
mean rotor induced velocity in hover, $\rho$ is the air density, $a_{1}$ is
the rotor disc area, $d_{\rm 0}$ is fuselage drag ratio, $s$ is rotor
solidity, and $v$ is the traveling velocity.

The energy efficiency (EE) performance metric in this paper is defined as the
ratio of the network throughput to the total power consumption, and we are
interested in investigating the influence of switching frequency $n_{t}$ on
the system performance.

######  Definition 6 (Energy Efficiency).

EE of a cellular network is defined in [23] as,

$\displaystyle{\rm EE}=\frac{{\rm SE}}{{\rm P_{\rm tot}}},$ (13)

where ${\rm SE}$ is spectral efficiency and $P_{\rm tot}$ is the total energy
consumption of UAVs and TBSs [35],

$\displaystyle{\rm SE}$ $\displaystyle=\text{b}_{\rm
w}\log_{2}{(1+\gamma)}(\lambda_{\rm u}P_{\rm cov,U}+\lambda_{\rm t}P_{\rm
cov,T}),$ (14) $\displaystyle P_{\rm tot}$ $\displaystyle=\lambda_{\rm
u}(\beta_{\rm tra}(n_{t})p_{m}+\beta_{\rm ser}(n_{t})p_{s})+\lambda_{\rm
t}p_{\rm tbs},$ (15)

in which $\text{b}_{\rm w}$ is the transmission bandwidth, $p_{\rm tbs}$ is
the overall power consumption of TBSs, $\beta_{\rm ser}(n_{t})$ and
$\beta_{\rm tra}(n_{t})$ are fractions of serving and traveling time of UAVs,
respectively.

##  III Distance Distribution

To capture the distribution of traveling distance $L$, which is the most
critical distance in our work, we first derive some prerequisite distance
distributions.

In the first scenario, since the user cluster pair is symmetrical, it is
enough to analyze the case where the cluster UAV locates at $C_{m}$ . With
that been said, even though we have two user clusters and two nearest charging
pads and cluster UAV needs to travel between them, the system performances
(both coverage probability and energy efficiency) of UAV locating at $C_{m}$
and $C_{n}$ are similar and influenced by the densities of user cluster
$x_{m}$ and $x_{n}$ only. Therefore, in the following parts of this paper, we
assume that the cluster UAV locates at $C_{m}$, and all the distances are
obtained based on $R_{\rm mm}$ and $\theta_{1}$, which is the distance between
$C_{m}$ and the center of $x_{m}$ and related angle.

######  Lemma 1 (Distribution of $R_{\rm mm}$ and $R_{\rm nn}$).

As mentioned in Definition 3, $R_{\rm mm}$ is assumed to be less than
$d_{nm}$. Thus, the probability density function (PDF) of $R_{\rm mm}$ is,

$\displaystyle f_{\rm R_{\rm mm}}(r)$ $\displaystyle=\frac{2\pi\lambda_{\rm
c}re^{-\pi\lambda_{\rm c}r^{2}}}{1-e^{-\pi\lambda_{\rm c}d_{\rm
nm}^{2}}},\quad 0\leq r\leq d_{\rm nm}.$ (16)

Let $R_{\rm nn}$ be the distance between $C_{n}$ and the center of $x_{n}$,
conditioned on $R_{\rm mm}$ and $\theta_{1}$, the PDF of $R_{\rm nn}$ is

$\displaystyle f_{\rm R_{\rm nn}}(r|R_{\rm
mm},\theta_{1})=\left\\{\begin{aligned} &2\lambda_{\rm
c}r\pi\exp(-\lambda_{\rm c}\pi r^{2}),\quad\text{\rm if}\quad
0<r<d_{nm}-R_{\rm mm},\\\ &2\lambda_{\rm
c}r(\pi-\theta_{d}(r))\exp\bigg{(}-\lambda_{\rm c}\bigg{(}\pi
r^{2}-\int_{d_{nm}-R_{\rm mm}}^{r}2\theta_{d}(z)z{\rm d}z\bigg{)}\bigg{)},\\\
&\text{\rm if}\quad d_{\rm nm}-R_{\rm mm}<r<R_{\rm mn},\\\
\end{aligned}\right.$ (17)

where $\theta_{d}$ is a function of $d_{nm}$ and $R_{\rm mm}$,

$\displaystyle\theta_{d}(r,R_{\rm mm})$ $\displaystyle=\arccos(\frac{d_{\rm
nm}^{2}+r^{2}-R_{\rm mm}^{2}}{2d_{\rm nm}r}).$ (18)

We then derive the distributions of the distances between the cluster UAV and
the reference user in $x_{m}$ and $x_{n}$, $R_{\rm um}$ and $R_{\rm un}$,
respectively, conditioned on $R_{\rm mm}$.

######  Lemma 2 (Distribution of $R_{\rm um}$ and $R_{\rm un}$).

In the first scenario, let $R_{\rm um}$ and $R_{\rm un}$ be the horizontal
distances between the cluster UAV and the reference user in user clusters
$x_{m}$ and $x_{n}$, respectively. Conditioned on $R_{\rm\\{mm,mn\\}}<r_{c}$,
the PDF of $R_{\rm um}$ and $R_{\rm un}$ are given by,

$\displaystyle f_{\rm R_{\\{um,un\\}}}(r|R_{\rm mm})=\left\\{\begin{aligned}
&\frac{2r}{r_{c}^{2}},\quad\text{\rm if}\quad 0<r<r_{c}-R_{\rm\\{mm,mn\\}},\\\
&\frac{2r}{\pi
r_{c}^{2}}\arccos\bigg{(}\frac{r^{2}-r_{c}^{2}+R_{\rm\\{mm,mn\\}}^{2}}{2R_{\rm\\{mm,mn\\}}r}\bigg{)},\text{\rm
if}\quad r_{c}-R_{\rm\\{mm,mn\\}}<r<R_{\rm\\{mm,mn\\}}+r_{c},\\\
\end{aligned}\right.$ (19)

otherwise, if $R_{\rm\\{mm,mn\\}}>r_{c}$, we have

$\displaystyle f_{\rm R_{\rm\\{um,un\\}}}(r|R_{\rm mm})$
$\displaystyle=\frac{2r}{\pi
r_{c}^{2}}\arccos(\frac{r^{2}-r_{c}^{2}+R_{\rm\\{mm,mn\\}}^{2}}{2R_{\rm\\{mm,mn\\}}r}),\quad\text{\rm
if}\quad R_{\rm\\{mm,mn\\}}-r_{c}<r<R_{\rm\\{mm,mn\\}}+r_{c},$ (20)

where the subscript $\\{um,un\\}$ corresponding to $\\{mm,mn\\}$,
respectively, and $R_{\rm mn}$ is the distance between $C_{m}$ and the center
of $x_{n}$, which can be represented as function of $R_{\rm mm}$ and
$\theta_{1}$

$\displaystyle R_{\rm mn}$ $\displaystyle=\sqrt{R_{\rm mm}^{2}+d_{\rm
nm}^{2}-2d_{\rm nm}R_{\rm mm}\cos(\theta_{1})}.$ (21)

Now that we have obtained all the required distance distributions, we can
derive the distribution of traveling distance $L$.

######  Theorem 1 (Distribution of $L$).

We first compute the probability of $L=0$, denoted by
$\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$,

$\displaystyle\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}=\exp(-\lambda_{\rm
c}({\rm Area}(R_{\rm mm},\theta_{1})),$ (22)

in which,

$\displaystyle{\rm Area}(R_{\rm mm},\theta_{1})=$ $\displaystyle R_{\rm
mm}^{2}(\pi-\theta_{1})+R_{\rm mn}^{2}(\pi-\theta_{3}(R_{\rm
mm},\theta_{1}))+R_{\rm mm}d_{\rm nm}\sin(\theta_{1}),$
$\displaystyle\theta_{3}(R_{\rm mm},\theta_{1})=$
$\displaystyle\arccos(\frac{d_{\rm nm}^{2}+R_{\rm mn}(R_{\rm
mm},\theta_{1})^{2}-R_{\rm mm}^{2}}{2d_{\rm nm}R_{\rm mn}(R_{\rm
mm},\theta_{1})}).$ (23)

Conditioned on $R_{\rm mm}$ and $\theta_{1}$, the CDF of L is given in $F_{\rm
L}(l)$, please refer to Appendix -A for more detail.

###### Proof:

See Appendix -B. ∎

######  Remark 2.

Note that the above distribution of L is conditioned on $R_{\rm mm}$ and
$\theta_{1}$, and we need to take the expectation over $R_{\rm mm}$, which is
a random variable. Hence, computing the unconditioned distribution of $L$ is
very complex and takes a long time. Therefore, in the following paper, we use
the expectation of $R_{\rm mm}$ to approximate the distribution of $L$. More
numerical detail will be explained in Section VI.

##  IV Coverage Probability

To capture the energy efficiency of the aforementioned two scenarios, we need
to obtain the coverage probabilities first. In the following theorem, we aim
to characterize the association probability and event. Recall the association
policy in Sec. II, the reference user associates with the cluster UAV or the
nearest TBS, which provides the strongest average received power.

######  Theorem 2 (Association Probability and Event).

Recall that $\mathcal{A}_{\rm NLoS}(r)$ and $\mathcal{A}_{\rm LoS}(r)$ are the
probabilities that the user associates with the cluster N/LoS UAV. Conditioned
on the cluster UAV is at distance $r$ away with channel type N/LoS, the
association probabilities are given as,

$\displaystyle\mathcal{A}_{\rm LoS}(r)$ $\displaystyle=\exp(-\pi\lambda_{\rm
t}d_{\rm lt}^{2}(r)),$ $\displaystyle\mathcal{A}_{\rm NLoS}(r)$
$\displaystyle=\exp(-\pi\lambda_{\rm t}d_{\rm nt}^{2}(r)),$ (24)

in which,

$\displaystyle d_{\rm lt}(r)$ $\displaystyle=(\rho_{\rm t}/\rho_{\rm
u})^{\frac{1}{\alpha_{\rm t}}}\eta_{\rm l}^{-\frac{1}{\alpha_{\rm
t}}}r^{\frac{\alpha_{\rm l}}{\alpha_{\rm t}}},$ $\displaystyle d_{\rm nt}(r)$
$\displaystyle=(\rho_{\rm t}/\rho_{\rm u})^{\frac{1}{\alpha_{\rm t}}}\eta_{\rm
n}^{-\frac{1}{\alpha_{\rm t}}}r^{\frac{\alpha_{\rm n}}{\alpha_{\rm t}}}.$ (25)

Recall $\mathcal{A}_{\rm TBS-N/LoS,\\{m,n\\}}(r)$ is the event that the user
associates with the TBS. Conditioned on the serving TBS is a distance $r$ away
and the UAV has the channel type N/LoS at distance $R_{\rm\\{um,un\\}}$ away,

$\displaystyle\mathcal{A}_{\rm TBS-
LoS,\\{m,n\\}}(r,R_{\rm\\{um,un\\}})=\mathbbm{1}\bigg{(}{\rm
LoS(R_{\rm\\{um,un\\}})}\bigg{)}\mathbbm{1}\bigg{(}R_{\rm\\{um,un\\}}>\sqrt{d_{\rm
tl}^{2}(r)-h^{2}}\bigg{)},$ $\displaystyle\mathcal{A}_{\rm TBS-
NLoS,\\{m,n\\}}(r,R_{\rm\\{um,un\\}})=\mathbbm{1}\bigg{(}{\rm
NLoS(R_{\rm\\{um,un\\}})}\bigg{)}\mathbbm{1}\bigg{(}R_{\rm\\{um,un\\}}>\sqrt{d_{\rm
tn}^{2}(r)-h^{2}}\bigg{)},$ (26)

where the subscript $\\{m,n\\}$ denote the location of the reference user, in
$x_{m}$ or $x_{n}$ and corresponding distance $R_{\rm um}$ or $R_{\rm un}$,
respectively,

$\displaystyle d_{\rm tl}(r)$ $\displaystyle=\max\bigg{(}(\rho_{\rm
u}/\rho_{\rm t})^{\frac{1}{\alpha_{\rm l}}}\eta_{\rm l}^{\frac{1}{\alpha_{\rm
l}}}r^{\frac{\alpha_{\rm t}}{\alpha_{\rm l}}},h\bigg{)},$ $\displaystyle
d_{\rm tn}(r)$ $\displaystyle=\max\bigg{(}(\rho_{\rm u}/\rho_{\rm
t})^{\frac{1}{\alpha_{\rm n}}}\eta_{\rm n}^{\frac{1}{\alpha_{\rm
n}}}r^{\frac{\alpha_{\rm t}}{\alpha_{\rm n}}},h\bigg{)}.$ (27)

###### Proof:

See Appendix -C. ∎

######  Remark 3.

Note that in Theorem 2, $\mathcal{A}_{\rm NLoS}(r)$ and $\mathcal{A}_{\rm
LoS}(r)$ are the association probabilities but $\mathcal{A}_{\rm
TBS-N/LoS,\\{m,n\\}}(r)$ is an association event. The reason is that the
association with a TBS is conditioned on $R_{\rm\\{um,un\\}}$, and should take
the expectation over $R_{\rm\\{um,un\\}}$. However, the Laplace transform of
the interference is also conditioned on $R_{\rm\\{um,un\\}}$. Hence, we write
this event as indicator function and take the expectation over
$R_{\rm\\{um,un\\}}$ together. The Laplace transform of the interference will
be introduced in the following lemma.

Laplace transform of the aggregate interference is the final requirement to
derive the coverage probability, which is provided in the following lemma.

######  Lemma 3 (Laplace Transform).

Conditioned on associating with the cluster UAV $bs$ at distance $r$, the
Laplace transform of the aggregate interference power is given by

$\displaystyle\mathcal{L}_{\rm I_{u},\\{s_{1},s_{2}\\}}(s,r)=$
$\displaystyle\exp\bigg{\\{}-2\pi\lambda_{\rm
u,\\{1,2\\}}\biggl{(}\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm
n}}{m_{\rm n}+s\eta_{\rm n}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm
n}}{2}}}\bigg{)}^{m_{\rm n}}\bigg{]}xP_{\rm n}(x){\rm d}x$
$\displaystyle\quad+\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm l}}{m_{\rm
l}+s\eta_{\rm l}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm
l}}{2}}}\bigg{)}^{m_{\rm l}}\bigg{]}xP_{\rm l}(x){\rm d}x\biggl{)}\bigg{\\}}$
$\displaystyle\times\exp\biggl{(}-2\pi\lambda_{\rm
t}\int_{t(r)}^{\infty}\bigg{[}1-\bigg{(}\frac{1}{1+s\rho_{\rm
t}x^{-\alpha_{\rm t}}}\bigg{)}\bigg{]}x{\rm d}x\biggl{)},$ (28)

where the subscript $\\{s_{1},s_{2}\\}$ denote the first or the second
scenario,

$\displaystyle t(r)=\left\\{\begin{aligned} d_{\rm nt}(r),&\quad\text{\rm
if}\quad bs\in\Phi_{\rm U_{n}},\\\ d_{\rm lt}(r),&\quad\text{\rm if}\quad
bs\in\Phi_{\rm U_{l}}.\\\ \end{aligned}\right.$ (29)

If the reference user associates with the nearest TBS $bs$, which is at
distance $r$, conditioned on the hotspot UAV having a LoS or NLoS link with
the user, the Laplace transforms of the aggregate interference are,
respectively, given by

$\displaystyle\mathcal{L}_{\rm I_{\rm
tl,\\{m,n\\}},\\{s_{1},s_{2}\\}}(s,r,R_{\rm\\{um,un\\}})=$
$\displaystyle\exp\bigg{\\{}-2\pi\lambda_{\rm
u,\\{1,2\\}}\biggl{(}\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm
n}}{m_{\rm n}+s\eta_{\rm n}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm
n}}{2}}}\bigg{)}^{m_{\rm n}}\bigg{]}$ $\displaystyle xP_{\rm n}(x){\rm
d}x+\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm l}}{m_{\rm l}+s\eta_{\rm
l}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm l}}{2}}}\bigg{)}^{m_{\rm
l}}\bigg{]}xP_{\rm l}(x){\rm d}x\biggr{)}\bigg{\\}}$
$\displaystyle\times\exp\biggl{(}-2\pi\lambda_{\rm
t}\int_{r}^{\infty}\bigg{[}1-\bigg{(}\frac{1}{1+s\rho_{\rm t}x^{-\alpha_{\rm
t}}}\bigg{)}\bigg{]}x{\rm d}x\biggl{)}$
$\displaystyle\times\bigg{(}\frac{m_{\rm l}}{m_{\rm l}+s\eta_{\rm l}\rho_{\rm
u}(R_{\rm\\{um,un\\}}^{2}+h^{2})^{-\frac{\alpha_{\rm l}}{2}}}\bigg{)}^{m_{\rm
l}},$ $\displaystyle\mathcal{L}_{\rm I_{\rm
tn,\\{m,n\\}},\\{s_{1},s_{2}\\}}(s,r,R_{\rm\\{um,un\\}})=$
$\displaystyle\exp\bigg{\\{}-2\pi\lambda_{\rm
u,\\{1,2\\}}\biggl{(}\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm
n}}{m_{\rm n}+s\eta_{\rm n}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm
n}}{2}}}\bigg{)}^{m_{\rm n}}\bigg{]}$ $\displaystyle xP_{\rm n}(x){\rm
d}x+\int_{0}^{\infty}\bigg{[}1-\bigg{(}\frac{m_{\rm l}}{m_{\rm l}+s\eta_{\rm
l}\rho_{\rm u}(x^{2}+h^{2})^{-\frac{\alpha_{\rm l}}{2}}}\bigg{)}^{m_{\rm
l}}\bigg{]}xP_{\rm l}(x){\rm d}x\biggr{)}\bigg{\\}}$
$\displaystyle\times\exp\biggl{(}-2\pi\lambda_{\rm
t}\int_{r}^{\infty}\bigg{[}1-\bigg{(}\frac{1}{1+s\rho_{\rm t}x^{-\alpha_{\rm
t}}}\bigg{)}\bigg{]}x{\rm d}x\biggl{)}$
$\displaystyle\times\bigg{(}\frac{m_{\rm n}}{m_{\rm n}+s\eta_{\rm n}\rho_{\rm
u}(R_{\rm\\{um,un\\}}^{2}+h^{2})^{-\frac{\alpha_{\rm l}}{2}}}\bigg{)}^{m_{\rm
n}}.$ (30)

Note that the subscripts $\\{s_{\rm 1},s_{\rm 2}\\}$ denote the first and the
second scenario, where the densities of UAVs are different, $\lambda_{\rm
u,1}$ and $\lambda_{\rm u,2}$, respectively,

###### Proof:

The Laplace transform is computed using MGF of Gamma distribution, PGFL of
inhomogeneous PPP, similar to [27] Lemma 7. ∎

Having developed all the expressions for the relevant distances and the
association probabilities, we study the coverage probability as explained in
Definition 5.

######  Theorem 3 (Coverage Probability of the First Scenario).

Recall that $q_{m,\\{h,l\\}}$ and $q_{n,\\{l,h\\}}$ are fractions of user
densities and denote the probabilities of the reference user being located in
the user cluster $x_{m}$ and $x_{n}$ at different times. Let $P_{\rm
cov,s_{1}|C_{m}}$ be the coverage probability when the cluster UAV is located
at $C_{m}$,

$\displaystyle P_{\rm cov,s_{1}|C_{m}}=q_{m,h}P_{\rm
cov,s_{1}|C_{m},q_{m,h}}+q_{n,l}P_{\rm cov,s_{1}|C_{m},q_{n,l}}$
$\displaystyle=q_{m,h}(P_{\rm cov,s_{1},U|C_{m},q_{m,h}}+P_{\rm
cov,s_{1},T|C_{m},q_{m,h}})+q_{n,l}(P_{\rm cov,s_{1},U|C_{m},q_{n,l}}+P_{\rm
cov,s_{1},T|C_{m},q_{n,l}}),$ (31)

in which, $P_{\rm cov,s_{1},U|C_{m},q_{m,h}}$ and $P_{\rm
cov,s_{1},U|C_{m},q_{n,l}}$ are the conditional coverage probabilities when
the reference user locates at $x_{m}$ or $x_{n}$, respectively,

$\displaystyle P_{\rm cov,s_{1},U|C_{m},q_{m,h}}$ $\displaystyle=P_{\rm
cov,s_{1},U_{l}|C_{m},q_{m,h}}+P_{\rm cov,s_{1},U_{n}|C_{m},q_{m,h}},$
$\displaystyle P_{\rm cov,s_{1},U|C_{m},q_{n,l}}$ $\displaystyle=P_{\rm
cov,s_{1},U_{l}|C_{m},q_{n,l}}+P_{\rm cov,s_{1},U_{n}|C_{m},q_{n,l}},$ (32)

where,

$\displaystyle P_{\rm cov,s_{1},U_{l}|C_{m},q_{m,h}}$
$\displaystyle=\int_{h}^{\sqrt{(d_{\rm nm}+r_{c})^{2}+h^{2}}}P_{\rm
l}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm LoS}(r)\sum_{k=0}^{m_{\rm
l}-1}\frac{(m_{\rm l}g_{\rm l}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{\rm u,s_{1}}}(s,r)|_{\rm s=m_{\rm l}g_{\rm
l}(r)}f_{\rm R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (33)
$\displaystyle P_{\rm cov,s_{1},U_{n}|C_{m},q_{m,h}}$
$\displaystyle=\int_{h}^{\sqrt{(d_{\rm nm}+r_{c})^{2}+h^{2}}}P_{\rm
n}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm NLoS}(r)\sum_{k=0}^{m_{\rm
n}-1}\frac{(m_{\rm n}g_{\rm n}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{\rm u,s_{1}}}(s,r)|_{\rm s=m_{\rm n}g_{\rm
n}(r)}f_{\rm R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (34)
$\displaystyle P_{\rm cov,s_{1},U_{l}|C_{m},q_{n,l}}$
$\displaystyle=\int_{h}^{\sqrt{(2d_{\rm nm}+r_{c})^{2}+h^{2}}}P_{\rm
l}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm LoS}(r)\sum_{k=0}^{m_{\rm
l}-1}\frac{(m_{\rm l}g_{\rm l}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{\rm u,s_{1}}}(s,r)|_{\rm s=m_{\rm l}g_{\rm
l}(r)}f_{\rm R_{\rm un}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (35)
$\displaystyle P_{\rm cov,s_{1},U_{n}|C_{m},q_{n,l}}$
$\displaystyle=\int_{h}^{\sqrt{(2d_{\rm nm}+r_{c})^{2}+h^{2}}}P_{\rm
n}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm NLoS}(r)\sum_{k=0}^{m_{\rm
n}-1}\frac{(m_{\rm n}g_{\rm n}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{\rm u,s_{1}}}(s,r)|_{\rm s=m_{\rm n}g_{\rm
n}(r)}f_{\rm R_{\rm un}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (36)

in which,

$\displaystyle g_{\rm l}(r)=\gamma(\rho_{\rm u}\eta_{\rm
l})^{-1}r^{\alpha_{\rm l}},$ $\displaystyle g_{\rm n}(r)=\gamma(\rho_{\rm
u}\eta_{\rm n})^{-1}r^{\alpha_{\rm n}},$ (37)

and,

$\displaystyle P_{\rm cov,s_{1},T|C_{m},q_{m,h}}$
$\displaystyle=\int_{0}^{\infty}2\pi\lambda_{\rm t}r\exp(-\pi\lambda_{\rm
t}r^{2})\bigg{(}\int_{\sqrt{d_{\rm tl}^{2}(r)-h^{2}}}^{R_{\rm mm}+r_{c}}f_{\rm
R_{um}}(z|R_{\rm mm})P_{\rm l}(z)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tl,m,s_{1}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r,z){\rm d}z$
$\displaystyle+\int_{\sqrt{d_{\rm tn}^{2}(r)-h^{2}}}^{R_{\rm mm}+r_{c}}f_{\rm
R_{um}}(z|R_{\rm mm})P_{\rm n}(z)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tn,m,s_{1}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r,z){\rm
d}z\bigg{)}{\rm d}r,$ (38) $\displaystyle P_{\rm cov,s_{1},T|C_{m},q_{n,l}}$
$\displaystyle=\int_{0}^{\infty}2\pi\lambda_{\rm t}r\exp(-\pi\lambda_{\rm
t}r^{2})\bigg{(}\int_{\sqrt{d_{\rm tl}^{2}(r)-h^{2}}}^{R_{\rm mn}+r_{c}}f_{\rm
R_{un}}(z|R_{\rm mm})P_{\rm l}(z)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tl,n,s_{1}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r,z){\rm d}z$
$\displaystyle+\int_{\sqrt{d_{\rm tn}^{2}(r)-h^{2}}}^{R_{\rm mn}+r_{c}}f_{\rm
R_{un}}(z|R_{\rm mm})P_{\rm n}(z)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tn,n,s_{1}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r,z){\rm
d}z\bigg{)}{\rm d}r.$ (39)

###### Proof:

Above equations derived by the fact that (i) the definition: $\bar{F}_{\rm
G}(g)=\frac{\Gamma_{u}(m,g)}{\Gamma(m)}$, where
$\Gamma_{u}(m,g)=\int^{\infty}_{mg}t^{m-1}e^{-t}dt$ is the upper incomplete
Gamma function, and (ii) the definition
$\frac{\Gamma_{u}(m,g)}{\Gamma(m)}=\exp(-g)\sum^{m-1}_{k=0}\frac{g^{k}}{k!}$,
similar to [27] Theorem 2. ∎

######  Remark 4.

As mentioned in Sec. III, it is enough to analyze the scenario where the
cluster UAV locates at $C_{m}$, and the difference of the coverage probability
between the cluster UAV located at $C_{m}$ and $C_{n}$ (the difference between
$P_{\rm cov,s_{1}|C_{n}}$ and $P_{\rm cov,s_{1}|C_{m}}$) is the probability of
the location of the reference user. That is, $q_{m,l}$ and $q_{n,h}$ are not
equal to $q_{n,l}$ and $q_{m,h}$, if the user density pairs are different.
Hence, the coverage probability when the cluster UAV located at $C_{n}$ is

$\displaystyle P_{\rm cov,s_{1}|C_{n}}$ $\displaystyle=q_{m,l}(P_{\rm
cov,s_{1},U|C_{n},q_{m,l}}+P_{\rm cov,s_{1},T|C_{n},q_{m,l}})+q_{n,h}(P_{\rm
cov,s_{1},U|C_{n},q_{n,h}}+P_{\rm cov,s_{1},T|s_{1},C_{n},q_{n,h}}),$ (40)

where $P_{\rm cov,s_{1},U|C_{n},q_{m,l}}$, $P_{\rm
cov,s_{1},T|C_{n},q_{m,l}}$, $P_{\rm cov,s_{1},U|C_{n},q_{n,h}}$ and $P_{\rm
cov,s_{1},T|s_{1},C_{n},q_{n,h}}$ are exactly the same as $P_{\rm
cov,s_{1},U|C_{m},q_{m,h}}$, $P_{\rm cov,s_{1},T|C_{m},q_{m,h}}$, $P_{\rm
cov,s_{1},U|C_{m},q_{n,l}}$ and $P_{\rm cov,s_{1},T|s_{1},C_{m},q_{n,l}}$,
respectively.

Now we present an approximation for the above equations by using the upper
bound of the CDF of the Gamma distribution in the case that the expressions of
coverage probability require evaluating higher-order derivatives of the
Laplace transform.

######  Lemma 4 (Approximated Coverage Probability).

$P_{\rm cov,U_{\rm\\{l,n\\}}|C_{m},q_{m,\\{h,l\\}}}$ and $P_{\rm
cov,U_{\\{l,n\\}}|C_{m},q_{n,\\{l,h\\}}}$ can be approximated by using the
upper bound of the CDF of the Gamma distribution as,

$\displaystyle P_{\rm cov,U_{\rm\\{l,n\\}}|C_{m},q_{m,h}}=$
$\displaystyle\int_{h}^{\sqrt{(d_{\rm
nm}+r_{c})^{2}+h^{2}}}P_{\rm\\{l,n\\}}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm\\{L,NL\\}oS}(r)\sum_{k=1}^{m_{\rm\\{l,n\\}}}\binom{m_{\rm\\{l,n\\}}}{k}$
$\displaystyle\quad\times(-1)^{k+1}\bigg{(}\mathcal{L}_{\rm\sigma^{2}+I_{\rm
u,s_{1}}}(k\beta_{2}m_{\rm\\{l,n\\}}g_{\rm\\{l,n\\}}(r)),r\bigg{)}f_{\rm
R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ $\displaystyle P_{\rm
cov,U_{\\{l,n\\}}|C_{m},q_{n,l}}$ $\displaystyle=\int_{h}^{\sqrt{(2d_{\rm
nm}+r_{c})^{2}+h^{2}}}P_{\rm
l}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm\\{L,NL\\}oS}(r)\sum_{k=1}^{m_{\rm\\{l,n\\}}}\binom{m_{\rm\\{l,n\\}}}{k}$
$\displaystyle\quad\times(-1)^{k+1}\bigg{(}\mathcal{L}_{\rm\sigma^{2}+I_{\rm
u,s_{1}}}(k\beta_{2}m_{\rm\\{l,n\\}}g_{\rm\\{l,n\\}}(r)),r\bigg{)}f_{\rm
R_{\rm un}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (41)

in which $\beta_{2}=(m_{\rm{\\{l,n\\}}}!)^{-\frac{1}{m_{\rm{\\{l,n\\}}}}}$.

###### Proof:

Approximate coverage probability is derived by using the definition of CCDF of
Gamma function, upper incomplete Gamma function and Binomial theorem, similar
to [27] Lemma 8. ∎

Following the same way, we obtain the expressions of the coverage probability
of the second scenario.

######  Theorem 4 (Coverage Probability of the Second Scenario).

The coverage probability of the second scenario is,

$\displaystyle P_{\rm cov,s_{2}}$ $\displaystyle=P_{\rm cov,s_{2},U}+P_{\rm
cov,s_{2},T}=P_{\rm cov,s_{2},U_{l}}+P_{\rm cov,s_{2},U_{n}}+P_{\rm
cov,s_{2},T},$ (42)

where,

$\displaystyle P_{\rm cov,s_{2},U_{l}}=$ $\displaystyle\int_{h}^{\sqrt{(d_{\rm
nm}+r_{c})^{2}+h^{2}}}P_{\rm l}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm
LoS}(r)\sum_{k=0}^{m_{\rm l}-1}\frac{(m_{\rm l}g_{\rm l}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{u,s_{2}}}(s)|_{\rm s=m_{\rm l}g_{\rm
l}(r)}f_{\rm R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (43)
$\displaystyle P_{\rm cov,s_{2},U_{n}}=$ $\displaystyle\int_{h}^{\sqrt{(d_{\rm
nm}+r_{c})^{2}+h^{2}}}P_{\rm n}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm
NLoS}(r)\sum_{k=0}^{m_{\rm n}-1}\frac{(m_{\rm n}g_{\rm n}(r))^{k}}{k!}$
$\displaystyle\quad\times(-1)^{k}\frac{\partial^{k}}{\partial
s^{k}}\mathcal{L}_{\rm\sigma^{2}+I_{u,s_{2}}}(s)|_{\rm s=m_{\rm n}g_{\rm
n}(r)}f_{\rm R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r,$ (44)
$\displaystyle P_{\rm cov,s_{2},T}=$
$\displaystyle\int_{0}^{\infty}2\pi\lambda_{\rm t}r\exp(-\pi\lambda_{\rm
t}r^{2})\bigg{(}\int_{\sqrt{d_{\rm tl}(r)^{2}-h^{2}}}^{d_{\rm
nm}+r_{c}}\mathcal{A}_{\rm TBS-LoS}(r)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tl,s_{2}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r|x)f_{\rm R_{um}}(x){\rm
d}x$ $\displaystyle+\int_{\sqrt{d_{\rm tn}(r)^{2}-h^{2}}}^{d_{\rm
nm}+r_{c}}\mathcal{A}_{\rm TBS-NLoS}(r)\mathcal{L}_{\rm\sigma^{2}+I_{\rm
tn,s_{2}}}(\gamma\rho_{\rm t}^{-1}r^{\alpha_{\rm t}},r|x)f_{\rm R_{um}}(x){\rm
d}x\bigg{)}{\rm d}r.$ (45)

Similar to the Theorem 3, we use the upper bound approximation to simlify the
high order derivation terms in the Laplace transform in Theorem 4 of the
second scenario.

######  Lemma 5 (Approximated Coverage Probability).

By using the upper bound, $P_{\rm cov,s_{2},U_{\\{l,n\\}}}$ can be written as,

$\displaystyle P_{\rm cov,s_{2},U_{\\{l,n\\}}}=$
$\displaystyle\int_{h}^{\sqrt{(d_{\rm
nm}+r_{c})^{2}+h^{2}}}P_{\rm\\{l,n\\}}(\sqrt{r^{2}-h^{2}})\mathcal{A}_{\rm\\{L,NL\\}oS}(r)\sum_{k=1}^{m_{\rm\\{l,n\\}}}\binom{m_{\rm\\{l,n\\}}}{k}(-1)^{k+1}$
$\displaystyle\quad\times\bigg{(}\mathcal{L}_{\rm\sigma^{2}+I_{u,s_{2}}}(k\beta_{2}m_{\rm\\{l,n\\}}g_{\rm\\{l,n\\}}(r))\bigg{)}f_{\rm
R_{\rm um}}(\sqrt{r^{2}-h^{2}}|R_{\rm mm}){\rm d}r.$ (46)

##  V Energy Efficiency

Finally, we compute the energy efficiency of these two scenarios. We first
derive the expressions of SE and total power consumption.

######  Lemma 6 (Spectral Efficiency and Total Power Consumption).

The SE (bit/sec/m2) can be formulated as,

$\displaystyle{\rm SE}_{\\{s_{1},s_{2}\\}}$ $\displaystyle=\text{b}_{\rm
w}\log_{2}{(1+\gamma)}(\lambda_{\rm u,\\{1,2\\}}P_{\rm
cov,\\{s_{1},s_{2}\\},U}+\lambda_{\rm t}P_{\rm cov,\\{s_{1},s_{2}\\},T}),$
(47)

in which $\text{b}_{\rm w}$ is the transmission bandwidth, while $P_{\rm
cov,s_{2},U}$ and $P_{\rm cov,s_{2},T}$ are defined in (43), (44) and (45),
$P_{\rm cov,s_{1},U}$ and $P_{\rm cov,s_{1},T}$ are

$\displaystyle P_{\rm cov,s_{1},U}=$ $\displaystyle\alpha(q_{m,h}P_{\rm
cov,s_{1},U|C_{m},q_{m,h}}+q_{n,l}P_{\rm cov,s_{1},U|C_{m},q_{n,l}})$
$\displaystyle+(1-\alpha)(q_{m,l}P_{\rm
cov,s_{1},U|C_{n},q_{m,l}}+q_{n,h}P_{\rm cov,s_{1},U|C_{n},q_{n,h}}),$ (48)
$\displaystyle P_{\rm cov,s_{1},T}=$ $\displaystyle\alpha(q_{m,h}P_{\rm
cov,s_{1},T|C_{m},q_{m,h}}+q_{n,l}P_{\rm cov,s_{1},T|C_{m},q_{n,l}})$
$\displaystyle+(1-\alpha)(q_{m,l}P_{\rm
cov,s_{1},T|C_{n},q_{m,l}}+q_{n,h}P_{\rm cov,s_{1},T|C_{n},q_{n,h}}).$ (49)

Total power consumption ${\rm P_{\rm tot}}$ (Watt/m2) can be formulated as,

$\displaystyle P_{\rm tot,s_{1}}$ $\displaystyle=\lambda_{\rm
u,1}\bigg{(}n_{t}\frac{2L}{24\times 3600v}p_{\rm m}+(1-n_{t}\frac{2L}{24\times
3600v})p_{\rm s}\bigg{)}+\lambda_{\rm t}p_{\rm tbs},$ (50) $\displaystyle
P_{\rm tot,s_{2}}$ $\displaystyle=\lambda_{\rm u,2}p_{\rm s}+\lambda_{\rm
t}p_{tbs},$ (51)

where $n_{t}$ is the frequency of traveling per day, $v$ is the velocity when
traveling and in numerical results section, we use the optimal $v$ which
minimizes the energy consumption, $p_{\rm m}$ and $p_{\rm s}$ are the motion-
related and service-related power, respectively.

###### Proof:

Recall that $\beta_{\rm ser}(n_{t})$ and $\beta_{\rm tra}(n_{t})$ are time
fractions, the ratio of serving, taveling time to total time of a day, given
by

$\displaystyle\beta_{\rm ser}(n_{t})$ $\displaystyle=n_{t}\frac{2L}{24\times
3600v},$ $\displaystyle\beta_{\rm tra}(n_{t})$ $\displaystyle=1-\beta_{\rm
ser}(n_{t}).$ (52)

Proof completes by substituting (52) into (15). ∎

Recall the energy efficiency introduced in Definition 6, which is the ratio of
the network throughput to the total power consumption.

######  Theorem 5 (Energy Efficiency).

EE in both scenarios are, respectively, given by,

$\displaystyle{\rm EE}_{1}$ $\displaystyle=\mathbb{E}_{\rm
L}\bigg{[}\frac{\text{b}_{\rm w}\log_{2}{(1+\gamma)}(\lambda_{\rm u,1}P_{\rm
cov,s_{1},U}+\lambda_{\rm t}P_{\rm cov,s_{1},T})}{\lambda_{\rm
u,1}\bigg{(}n_{t}\frac{2L}{24\times 3600v}p_{\rm m}+(1-n_{t}\frac{2L}{24\times
3600v})p_{\rm s}\bigg{)}+\lambda_{\rm t}p_{\rm tbs}}\bigg{]},$ (53)
$\displaystyle{\rm EE}_{2}$ $\displaystyle=\frac{\text{b}_{\rm
w}\log_{2}{(1+\gamma)}(\lambda_{\rm u,2}P_{\rm cov,s_{1},U}+\lambda_{\rm
t}P_{\rm cov,s_{1},T})}{\lambda_{\rm u,2}p_{\rm s}+\lambda_{\rm t}p_{tbs}}.$
(54)

where $n_{t}$ is the switching frequency, which is the number of times the
cluster UAV travels back and forth between $C_{m}$ and $C_{n}$ per day.

###### Proof:

The above equations follow by substituting (47), (50) and (51) into (13). ∎

######  Remark 5.

Computing (53) requires the PDF of $L$, which needs to take the derivative of
the CDF of $L$ given in Appendix -A. Here we use the expectation to
approximate (53), now the energy efficiency of the first scenario is given as

$\displaystyle{\rm EE}_{1}$ $\displaystyle=\frac{\text{b}_{\rm
w}\log_{2}{(1+\gamma)}(\lambda_{\rm u,1}P_{\rm cov,s_{1},U}+\lambda_{\rm
t}P_{\rm cov,s_{1},T})}{\lambda_{\rm
u,1}\bigg{(}n_{t}\frac{2\mathbb{E}[L]}{24\times 3600v}p_{\rm
m}+(1-n_{t}\frac{2\mathbb{E}[L]}{24\times 3600v})p_{\rm
s}\bigg{)}+\lambda_{\rm t}p_{\rm tbs}}.$ (55)

##  VI Numerical Results

In this section, we validate our analytical results with simulations and
evaluate the impact of various system parameters on the network performance.
To simplify the selection of parameters, we set $\lambda_{\rm
m,l}=\lambda_{\rm n,l}$ and $\lambda_{\rm m,h}=\lambda_{\rm n,h}$.

Unless stated otherwise, we use the simulation parameters as listed herein
Table II.

TABLE II: Table of Parameters

Parameter | Symbol | Simulation Value
---|---|---
Density of TBSs | $\lambda_{\rm t}$ |  $10^{-6}$ m-2
Traveling-related power | $P_{\rm m}$ | 161.8 W
Service-related power | $P_{\rm s}$ | 10 W
UAV altitude | $h$ | 60 m
UAV velocity | $v$ | 18.46 m/s
Radius of MCP disk | $r_{c}$ | 120 m
N/LoS environment variable | $a,b$ | 25.27,0.2
Transmission power of UAVs and TBSs |  $\rho_{\rm u}$, $\rho_{\rm t}$ | 0.4 W, 10 W
Total power comsuption of TBSs | $p_{t}$ | 318 W
SINR threshold | $\gamma$ | 0 dB
Noise power | $\sigma^{2}$ |  $10^{-9}$ W
N/LoS and active charging station path-loss exponent | $\alpha_{\rm n},\alpha_{\rm l},\alpha_{\rm t}$ | $4,2.1,4$
N/LoS fading gain | $m_{\rm n},m_{\rm l}$ | $1,3$
N/LoS additional loss | $\eta_{\rm n},\eta_{\rm l}$ |  $-20,0$ dB

For the simulation of the considered system model, we first compute the
accuracy of the approximation of traveling distance $L$, as mentioned in
Remark 2. We simulate two cases, $d_{\rm nm}=300$ and $d_{\rm nm}=600$. The
traveling distance increases with the increase of $d_{\rm nm}$. In both cases,
the results in Fig. 3 show that using the expectation of $R_{\rm mm}$ provides
a close approximation of the distribution of $L$. We then use the approximate
distribution of $L$ in the following part of this paper.

Figure 3: Using the approximation, the distance distribution of $L$,
conditioned on $\lambda_{c}=10$ charging pad/km2.  Figure 4: Coverage
probabilities under different densities of charging pads, conditioned on
$\lambda_{\rm t}=1$ TBS/km2, $\lambda_{u}=10$ user cluster pair/km2, and
$d_{\rm nm}=300$ m. Figure 5: Coverage probabilities under different densities
of user cluster pairs, conditioned on $\lambda_{\rm c}=100$ charging pad/km2,
$\lambda_{t}=1$ TBS/km2, and $d_{\rm nm}=300$ m.

Next, Fig. 4 and Fig. 5 show the coverage probabilities increase with the
increase of the density of charging pads, and decrease when increasing the
density of user cluster pairs. Both two figures show that the second scenario
is not always better than the first one, especially when the density of one
user cluster is much higher than that of the other one. The system performance
of one UAV per cluster pair is better, if the user density of one user cluster
is much higher than another one. The reason is that in the second scenario,
even though we have one UAV for each user cluster, and the reference user gets
closer to the serving cluster UAV. However, the interference becomes higher
due to the increasing density of UAVs compared with the first scenario. In
addition, the interference has a higher impact on the coverage probabilities
than distance.

Figure 6: Energy efficiency under different densities of charging pads,
conditioned on $\lambda_{\rm u}=10$ user cluster pair/km2, $\lambda_{t}=1$
TBS/km2, and $d_{\rm nm}=300$ m.

Finally, we plot the energy efficiency against $\lambda_{\rm c}$ and
$\lambda_{\rm u}$ in Fig. 6 and Fig. 7, respectively. Note that, the gaps
between simulation and analysis in energy efficiency when $n_{t}$ is large
($n_{t}=200$) are due to the approximation we used in Remark 5. As expected,
energy efficiency increases with the increase of the densities of charging
pads due to the increase of the coverage probability. The system’s energy
efficiency decreases with the increase of $d_{\rm nm}$ since the traveling-
related energy consumption is higher than service-related energy. Besides, if
two user clusters are too far from each other, we cannot ignore the impact of
traveling time on the coverage probability. Fig. 6 shows that the second
scenario has a better performance when $\lambda_{\rm c}$ is high, that is
because of the higher densities of UAVs. With the increase of the density of
charging pads, the probability of $L=0$ decreases dramatically. Hence the
densities of UAV, as well as SE, becomes larger compared with the first
scenario.

Figure 7: Energy efficiency under different densities of user cluster pairs,
conditioned on $\lambda_{\rm c}=100$ charging pad/km2, $\lambda_{t}=1$
TBS/km2, and $d_{\rm nm}=300$ m.

Fig. 7 shows the influence of $\lambda_{\rm u}$. Energy efficiency increases
first due to the increase of coverage probability and densities of UAVs, and
then decreases owing to the decreasing of coverage probability because of the
higher interference.

The results reveal an interesting system insight that two-UAV deployment
strategy does not always achieve a better performance, which also depends on
the density of users and user clusters. For instance, when the density of one
user cluster is much greater than the other, the user cluster density is very
high or the density of charging pads is low, deploying one UAV per cluster
pair is better.

##  VII Conclusion

This paper presented a novel system model that considered the dynamic density
of users and captured the energy efficiency in two scenarios, deploying one
UAV or two UAVs per cluster pair. Firstly, we derived the PDF of traveling
distance. Next, we computed the coverage probability of the above two
scenarios. Finally, we used the aforementioned distance and coverage
probability to compute the energy efficiency of the considered system model.
We have validated that our distance distribution is accurate and shown that
one UAV deployment is better than two UAVs deployment if the density of the
user cluster is high or the density of charging pads is low.

Our work investigated some new aspects of the performance of UAV-enabled
wireless networks, such as the dynamic variations of user clusters and one or
two UAVs deployment, and obtained a new distance distribution. By considering
the relationships between the densities of user cluster pairs and charging
pads, optimal energy efficiency can be achieved by deploying different numbers
of UAVs. A more complicated model can be developed in future work, such as the
optimization of UAV locations.

###  -A Distribution of Traveling Distance L

We are interested in obtaining the distance distribution of L, which is shown
in Fig. 2, the distance between $C_{m}$ and $C_{n}$, which are the nearest
charging pads to the centers of user cluster $x_{m}$ and $x_{n}$,
respectively. Conditioned on $R_{\rm mm}$ and $\theta_{1}$, the CDF of $L$ is
given by,

$\displaystyle F_{\rm L}(l)=\mathbb{P}_{\rm(0|R_{\rm
mm},\theta_{1})},\quad\text{if}\quad l=0,$ (56)

otherwise,

(i) in the case of $\theta_{1}+\theta_{3}(R_{\rm mm},\theta_{1})<\pi/2$,

  * •

if $0<l\leq R_{\rm mn}+d_{\rm nm}-R_{\rm mm}$, $F_{\rm L}(l)$ is given by

$\int_{R_{\rm mn}-l}^{R_{2d}}f_{3}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $R_{\rm mn}+d_{\rm nm}-R_{\rm mm}<l\leq 2R_{\rm
mm}\sin(\frac{\theta_{1}}{2})$, $F_{\rm L}(l)$ is given by

$\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{3}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $2R_{\rm mm}\sin(\frac{\theta_{1}}{2})\leq l\leq\min(2R_{\rm
mm}\sin(\theta_{1}),R_{\rm mn})$, $F_{\rm L}(l)$ is given by

$e^{-\lambda_{\rm c}\pi(R_{\rm mn}-l)^{2}}-\int_{R_{\rm mn}-l}^{d_{\rm
nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm mm},\theta_{1}){\rm d}r-\int_{d_{\rm
nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r-\int_{R_{2d}}^{R3}f_{9}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $2R_{\rm mm}\sin(\theta_{1})\leq l\leq R_{\rm mn}$, $F_{\rm L}(l)$ is given
by

$\exp(-\lambda_{\rm c}\pi(R_{\rm mn}-l)^{2})-\int_{R_{\rm mn}-l}^{d_{\rm
nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm mm},\theta_{1}){\rm d}r-\int_{d_{\rm
nm}-R_{\rm mm}}^{R_{3}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $R_{\rm mn}<l\leq\min(R_{\rm mn}+d_{\rm nm}-R_{\rm mm},2R_{\rm
mm}\sin(\theta_{1}))$, $F_{\rm L}(l)$ is given by

$1-\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{R_{2d}}^{R_{\rm mn}}f_{9}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r$

  * •

if $R_{\rm mn}+d_{\rm nm}-R_{\rm mm}<l\leq 2R_{\rm mm}\sin(\theta_{1})$,
$F_{\rm L}(l)$ is given by

$1-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r-\int_{R_{2d}}^{R_{\rm mn}}f_{9}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $\max(R_{\rm mn},2R_{\rm mm}\sin(\theta_{1}))\leq l\leq R_{\rm mn}+d_{\rm
nm}-R_{\rm mm}$, $F_{\rm L}(l)$ is given by

$1-\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{\rm
mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $\max(2R_{\rm mm}\sin(\theta_{1}),R_{\rm mn}+d_{\rm nm}-R_{\rm mm})<l\leq
2R_{\rm mn}$, $F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{R_{\rm mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

(ii) in the case of $\cos(\pi-\theta_{1}-\theta_{3}(R_{\rm
mm},\theta_{1}))>\sin(\theta_{1})$,

  * •

if $0<l\leq 2R_{\rm mm}\sin(\theta_{1})$, $F_{\rm L}(l)$ is given by

$\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $2R_{\rm mm}\sin(\theta_{1})<l\leq 2R_{\rm
mm}\cos(\pi-\theta_{1}-\theta_{3}(R_{\rm mm},\theta_{1})),$ $F_{\rm L}(l)$ is
given by

$\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{R_{2d}}^{R_{\rm mn}}f_{7}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $2R_{\rm mm}\cos(\pi-\theta_{1}-\theta_{3}(R_{\rm
mm},\theta_{1}))<l\leq\min(R_{\rm mn}+d_{\rm nm}-R_{\rm mm},2R_{\rm mm}),$
$F_{\rm L}(l)$ is given by

$\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{R_{2d2}}^{R_{\rm mn}}f_{7}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{l-R_{\rm mn}}^{d_{\rm nm}-R_{\rm mm}}f_{1}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r$ $+\int_{d_{\rm nm}-R_{\rm
mm}}^{R_{2d}}f_{5}(r,l,R_{\rm mm},\theta_{1}){\rm d}r+\mathbb{P}_{\rm(0|R_{\rm
mm},\theta_{1})}$

  * •

if $R_{\rm mn}+d_{\rm nm}-R_{\rm mm}<l\leq 2R_{\rm mm},$ $F_{\rm L}(l)$ is
given by

$\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{R_{2d2}}^{R_{\rm mn}}f_{7}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\int_{l-R_{\rm mn}}^{R_{2d}}f_{5}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $2R_{\rm mm}<l\leq R_{\rm mn},$ $F_{\rm L}(l)$ is given by

$\exp(-\lambda_{\rm c}\pi(R_{\rm mn}-l)^{2})-\int_{R_{\rm mn}-l}^{d_{\rm
nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm mm},\theta_{1}){\rm d}r-\int_{d_{\rm
nm}-R_{\rm mm}}^{R3}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $R_{\rm mn}<l\leq R_{\rm mn}+d_{\rm nm}-R_{\rm mm},$ $F_{\rm L}(l)$ is
given by

$1-\int_{l-R_{\rm mn}}^{d_{\rm nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{\rm
mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $\max(2R_{\rm mm},R_{\rm mn}+d_{\rm nm}-R_{\rm mm})<l\leq 2R_{\rm mn},$
$F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{R_{\rm mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

(iii) in the case of $\theta_{1}+\theta_{3}(R_{\rm
mm},\theta_{1})\geq\pi/2\quad\text{and}\quad\cos(\pi-\theta_{1}-\theta_{3}(R_{\rm
mm},\theta_{1}))\leq\sin(\theta_{1})$,

  * •

if $0<l\leq 2R_{\rm mm}\cos(\pi-\theta_{1}-\theta_{3})$, $F_{\rm L}(l)$ is
given by

$\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $2R_{\rm mm}\cos(\pi-\theta_{1}-\theta_{3})<l\leq 2R_{\rm
mm}\sin(\frac{\theta_{1}}{2}),$ $F_{\rm L}(l)$ is given by

$\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{1}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{3}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

$2R_{\rm mm}\sin(\frac{\theta_{1}}{2})<l\leq\min(2R_{\rm
mm}\sin(\theta_{1}),R_{\rm mn}),$ $F_{\rm L}(l)$ is given by

$\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{1}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r+\int_{R_{2d}}^{R_{\rm mn}}f_{6}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r\\\ +\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$

  * •

if $R_{\rm mn}<l\leq\min(2R_{\rm mm}\sin(\theta_{1}),R_{\rm mn}+d_{\rm
nm}-R_{\rm mm}),$ $F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{d_{\rm nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{4}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{R_{2d}}^{R_{\rm mn}}f_{9}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r$

  * •

if $R_{\rm mn}+d_{\rm nm}-R_{\rm mm}<l\leq 2R_{\rm mm}\sin(\theta_{1}),$
$F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{R_{2d}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r-\int_{R_{2d}}^{R_{\rm mn}}f_{9}(r,l,R_{\rm mm},\theta_{1}){\rm d}r$

  * •

if $2R_{\rm mm}\sin(\theta_{1})<l\leq R_{\rm mn}$, $F_{\rm L}(l)$ is given by

$\exp(-\lambda_{\rm c}\pi(R_{\rm mn}-l)^{2})-\int_{R_{\rm mn}-l}^{d_{\rm
nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm mm},\theta_{1}){\rm d}r-\int_{d_{\rm
nm}-R_{\rm mm}}^{R3}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r+g(l,R_{\rm
mm},\theta_{1})$

  * •

if $\max(R_{\rm mn},2R_{\rm mm}\sin(\theta_{1}))<l\leq R_{\rm mn}+d_{\rm
nm}-R_{\rm mm},$ $F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{d_{\rm nm}-R_{\rm mm}}f_{2}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r-\int_{d_{\rm nm}-R_{\rm mm}}^{R_{\rm
mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm d}r+g(l,R_{\rm mm},\theta_{1})$

  * •

if $\max(R_{\rm mn}+d_{\rm nm}-R_{\rm mm},2R_{\rm mm}\sin(\theta_{1}))<l\leq
2R_{\rm mn}$, $F_{\rm L}(l)$ is given by

$1-\int_{l-R_{\rm mn}}^{R_{\rm mn}}f_{4}(r,l,R_{\rm mm},\theta_{1}){\rm
d}r+g(l,R_{\rm mm},\theta_{1})$

where $g(l,R_{\rm mm},\theta_{1})=\int_{R_{2d}}^{R_{2d2}}f_{8}(r,l,R_{\rm
mm},\theta_{1}){\rm d}r$ if $\theta_{1}>\frac{\pi}{2}$ and $l<2R_{\rm mm}$,
otherwise, $0$,

TABLE III: Table of Notations

Function | Definition
---|---
$f_{1}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\theta(r,l,R_{\rm mm},\theta_{1})}{\pi}$
$f_{2}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\pi-\theta(r,l,R_{\rm mm},\theta_{1})}{\pi}$
$f_{3}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\theta(r,l,R_{\rm mm},\theta_{1})}{\pi-\theta_{d}(r,R_{\rm mm},\theta_{1})}$
$f_{4}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\pi-\theta(r,l,R_{\rm mm},\theta_{1})}{\pi-\theta_{d}(r,l,R_{\rm mm},\theta_{1})}$
$f_{5}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\theta(r,l,R_{\rm mm},\theta_{1})-\theta_{d}(r,R_{\rm mm},\theta_{1})}{\pi-\theta_{d}(r,R_{\rm mm},\theta_{1})}$
$f_{6}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\theta(r,l,R_{\rm mm},\theta_{1})+\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{d}(r,R_{\rm mm},\theta_{1})}{2\pi-2\theta_{d}(r,R_{\rm mm},\theta_{1})}$
$f_{7}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{\theta(r,l,R_{\rm mm},\theta_{1})-\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{d}(r,R_{\rm mm},\theta_{1})}{2\pi-2\theta_{d}(r,R_{\rm mm},\theta_{1})}$
$f_{8}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{-\theta(r,l,R_{\rm mm},\theta_{1})+\theta_{3}(R_{\rm mm},\theta_{1})+\theta_{d}(r,R_{\rm mm},\theta_{1})}{2\pi-2\theta_{d}}$
$f_{9}(r,l,R_{\rm mm},\theta_{1})$ | $f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{2\pi-\theta(r,l,R_{\rm mm},\theta_{1})-\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{d}(r,R_{\rm mm},\theta_{1})}{2\pi-2\theta_{d}(r,R_{\rm mm},\theta_{1})}$
$\theta_{3}(R_{\rm mm},\theta_{1})$ | $\arccos(\frac{d_{\rm nm}^{2}+R_{\rm mn}(R_{\rm mm},\theta_{1})^{2}-R_{\rm mm}^{2}}{2d_{\rm nm}R_{\rm mn}(R_{\rm mm},\theta_{1})})$
$\theta(r,l,R_{\rm mm},\theta_{1})$ | $\arccos(\frac{R_{\rm mn}(R_{\rm mm},\theta_{1})^{2}+r^{2}-l^{2}}{2d_{\rm nm}R_{\rm mn}(R_{\rm mm},\theta_{1})})$
$\theta_{d}(r,R_{\rm mm},\theta_{1})$ | $\arccos(\frac{d_{\rm nm}^{2}+r^{2}-R_{\rm mm}^{2}}{2d_{\rm nm}r})$
$\theta_{4}(l,R_{\rm mm},\theta_{1})$ | $\pi-\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{1}-\arccos(\frac{l}{2R_{\rm mm}})$
$\theta_{5}(l,R_{\rm mm},\theta_{1})$ | $\pi-\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{1}+\arccos(\frac{l}{2R_{\rm mm}})$
$R_{2d}(l,R_{\rm mm},\theta_{1})$ | $\sqrt{(l\cos(\theta_{4})-R_{\rm mn}(R_{\rm mm},\theta_{1}))^{2}+(l\sin(\theta_{4}))^{2}}$
$R_{2d2}(l,R_{\rm mm},\theta_{1})$ | $\sqrt{(l\cos(\theta_{5})-R_{\rm mn}(R_{\rm mm},\theta_{1}))^{2}+(l\sin(\theta_{5}))^{2}}$

###  -B Proof of the Distribution of Traveling Distance L

In this appendix, we provide proofs two cases: (i)
$\mathbb{P}(\mathcal{N}(A_{0})=0|R_{\rm mm},\theta_{1})$ and (ii) one case in
$F_{\rm L}(l)$, which is condtioned on $\theta_{1}+\theta_{3}(R_{\rm
mm},\theta_{1})<\pi/2$ and $2R_{\rm mm}\sin(\frac{\theta_{1}}{2})\leq
l\leq\min(2R_{\rm mm}\sin(\theta_{1}),R_{\rm mn})$.

####  -B1 Proof of $\mathbb{P}(\mathcal{N}(A_{0})=0|R_{\rm mm},\theta_{1})$

Figure 8: The possible location of $C_{n}$, such that $L$ is greater than 0.

Note that in our system, we assume that $L$ is the distance between the two
charging pads. Therefore, if $C_{m}$ locates at $R_{\rm mm}$ away and
corresponding angle is $\theta_{1}$, $C_{n}$ can only be located in ${\rm
A_{0}}$, which is colored in Fig. 8. Hence, the probability of $L=0$ is the
probability of no point falling in $A_{0}$,

$\displaystyle\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$
$\displaystyle=\mathbb{P}(\mathcal{N}(A_{0})=0|R_{\rm mm},\theta_{1})$ (57)

obviously, the area of $A_{0}$ is given by

$\displaystyle|A_{0}|=\pi R_{\rm mm}^{2}\frac{\pi-\theta_{1}}{\pi}+\pi R_{\rm
mn}^{2}\frac{\pi-\theta_{3}}{\pi}+R_{\rm mm}d_{\rm nm}\sin\theta_{1}-\pi
R_{\rm mm}^{2},$ (58)

$\theta_{3}$ and $R_{\rm mn}$ can be represented by $\theta_{1}$ and $R_{\rm
mm}$ as,

$\displaystyle R_{\rm mn}$ $\displaystyle={\rm abs}\bigg{(}\frac{d_{\rm
nm}-R_{\rm mm}\cos\theta_{1}}{\cos\theta_{3}}\bigg{)},$
$\displaystyle\theta_{3}$ $\displaystyle=\arccos(\frac{d_{\rm nm}^{2}+R_{\rm
mn}^{2}-R_{\rm mm}^{2}}{2d_{\rm nm}R_{\rm mn}(R_{\rm mm},\theta_{1})}),$

therefore, $\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$ now becomes,

$\displaystyle\mathbb{P}_{\rm(0|R_{\rm mm},\theta_{1})}$
$\displaystyle=\exp(-\lambda_{c}|A_{0}(R_{\rm mm},\theta_{1})|),$ (59)

where $|A_{0}|$ denotes the area of $A_{0}$.

####  -B2 Proof of One Case In The Distribution of $L$

In this part, we only prove one case which is $\theta_{1}+\theta_{3}(R_{\rm
mm},\theta_{1})<\pi/2$, $2R_{\rm mm}\sin(\frac{\theta_{1}}{2})\leq
l\leq\min(2R_{\rm mm}\sin(\theta_{1}),R_{\rm mn})$ and $2R_{\rm
mm}\sin(\theta_{1})<R_{\rm mn}$, as shown in Fig. 8. Two red dash lines denote
$2R_{\rm mm}\sin(\frac{\theta_{1}}{2})$ and $2R_{\rm mm}\sin(\theta_{1})$,
respectively.

The probability of $L<l$ is the intersection of the red ball $B(C_{m},l)$ and
$A_{0}$,

$\displaystyle F_{L}(l)$ $\displaystyle=\mathbb{P}(L<l)=1-\mathbb{P}(L>l),$
(60)

where $\mathbb{P}(L>l)$ is the probability of $C_{n}$ falling in $A_{1}$.

Figure 9: Area of L not equal to 0, $\mathbb{P}_{\rm(0|R_{\rm
mm},\theta_{1})}$.

We then compute the probability of $C_{n}$ is located in $A_{1}$. Note that
$A_{1}$ is consisted of $4$ areas, shown in Fig. 9. The probability of $C_{n}$
located in these 4 areas are different, we need to compute them seperately,

$\displaystyle\mathbb{P}(L>l)$
$\displaystyle=\mathbb{P}(A_{2})+\mathbb{P}(A_{3})+\mathbb{P}(A_{4})+\mathbb{P}(A_{5}),$
(61) $\displaystyle\mathbb{P}(A_{2})$ $\displaystyle=1-e^{-\lambda_{\rm
c}\pi(R_{\rm mn}-l)^{2}},$ $\displaystyle\mathbb{P}(A_{3})$
$\displaystyle=\int_{R_{\rm mn}-l}^{d_{\rm nm}-R_{\rm mm}}f_{\rm R_{\rm
nn}}(r|R_{\rm mm},\theta_{1})\frac{\pi-\theta(r,l,R_{\rm
mm},\theta_{1})}{\pi}{\rm d}r,$ $\displaystyle\mathbb{P}(A_{4})$
$\displaystyle=\int_{d_{\rm nm}-R_{\rm mm}}^{R_{2d}}f_{\rm R_{\rm
nn}}(r|R_{\rm mm},\theta_{1})\frac{\pi-\theta(r,l,R_{\rm
mm},\theta_{1})}{\pi-\theta_{d}(r,l,R_{\rm mm},\theta_{1})}{\rm d}r,$
$\displaystyle\mathbb{P}(A_{5})$ $\displaystyle=\int_{R_{2d}}^{R_{\rm
mn}}f_{\rm R_{\rm nn}}(r|R_{\rm mm},\theta_{1})\frac{2\pi-\theta(r,l,R_{\rm
mm},\theta_{1})-\theta_{3}(R_{\rm mm},\theta_{1})-\theta_{d}(r,R_{\rm
mm},\theta_{1})}{2\pi-2\theta_{d}(r,R_{\rm mm},\theta_{1})}{\rm d}r,$ (62)

where $\mathbb{P}(A_{2})$ is simply derived by 1 minus the probability of no
point falling in $A_{1}$. Taking $C_{n}$ falls in $A_{3}$ for example, the
probability of $C_{n}$ located at $r$ away is

$\displaystyle\mathbb{P}(R_{\rm nn}=r)$ $\displaystyle=2\lambda_{\rm
c}r\exp(-r\pi^{2})\frac{2\pi-2\theta}{2\pi}=f_{\rm
R_{nn}}(r)\frac{\pi-\theta}{\pi},$ (63)

where $f_{\rm R_{nn}}(r)$ is the distance distribuiton of $R_{\rm nn}$,
$\frac{\pi-\theta}{\pi}$ denotes the area that $C_{n}$ can be located (except
the uncolored area near $A_{3}$), then taking the integration over $r$ from
$R_{\rm mn}-l$ to $d_{\rm nm}-R_{\rm mm}$ (the borders of $A_{3}$) completes
the proof.

Figure 10: Area of L not equal to 0, $\mathbb{P}_{\rm(0|R_{\rm
mm},\theta_{1})}$.

Note that in $A_{4}$, the minus term $\theta_{d}$ in denominator is owing to
the distance distribution of $R_{\rm nn}$, because we have already considered
that the probability of $C_{n}$ falling in the intersection of $B(x_{m},R_{\rm
mm})$ and $B(x_{n},R_{\rm mn})$ equals to 0.

The proof of the probabilities of $C_{n}$ falling in other areas follows the
same method, and the geometry relationships are shown in Fig. 9.

###  -C Proof of Lemma 2

If the reference user and the cluster UAV establish a LoS link with distance
$r$, the association probability $\mathcal{A}_{\rm LoS}(r)$ is

$\displaystyle\mathcal{A}_{\rm LoS}(r)$ $\displaystyle=\mathbb{P}(p_{\rm
l}>p_{\rm t})=\mathbb{P}\bigg{(}\eta_{\rm l}\rho_{\rm u}r^{-\alpha_{\rm
l}}>\rho_{\rm t}R_{\rm t}^{-\alpha_{\rm t}}\bigg{)}$
$\displaystyle=\mathbb{P}\bigg{(}R_{\rm t}>(\rho_{\rm t}/\rho_{\rm
u})^{\frac{1}{\alpha_{\rm t}}}\eta_{\rm l}^{-\frac{1}{\alpha_{\rm
t}}}r^{\frac{\alpha_{\rm l}}{\alpha_{\rm
t}}}\bigg{)}=\exp\bigg{(}-\pi\lambda_{\rm t}d_{\rm lt}^{2}(r)\bigg{)},$ (64)

in which $d_{\rm lt}(r)=(\rho_{\rm t}/\rho_{\rm u})^{\frac{1}{\alpha_{\rm
t}}}\eta_{\rm l}^{-\frac{1}{\alpha_{\rm t}}}r^{\frac{\alpha_{\rm
l}}{\alpha_{\rm t}}}$.

Recall that $\mathcal{A}_{\rm TBS-LoS,\\{m,n\\}}(r,R_{\rm mm})$ denotes the
event that the average received power from the nearest TBS, which locates at
$r$ away, is stronger than cluster N/LoS UAV that locates at $R_{\rm um}$
away,

$\displaystyle\mathcal{A}_{\rm TBS-LoS,m}(r,R_{\rm um})$
$\displaystyle=\mathbbm{1}\bigg{(}{\rm LoS(R_{\rm
um})}\bigg{)}\mathbbm{1}\bigg{(}(p_{\rm t}>p_{\rm l})\bigg{)}$
$\displaystyle=\mathbbm{1}\bigg{(}{\rm LoS(R_{\rm
um})}\bigg{)}\mathbbm{1}\bigg{(}\rho_{\rm t}r^{-\alpha_{\rm t}}>\eta_{\rm
l}\rho_{\rm u}(R_{\rm um}^{2}+h^{2})^{\frac{-\alpha_{\rm l}}{2}}\bigg{)}$
$\displaystyle=\mathbbm{1}\bigg{(}{\rm LoS(R_{\rm
um})}\bigg{)}\mathbbm{1}\bigg{(}R_{\rm um}>\sqrt{d_{\rm
tl}^{2}(r)-h^{2}}\bigg{)},$ (65)

in which $d_{\rm tl}(r)=\max((\rho_{\rm u}/\rho_{\rm t})^{\frac{1}{\alpha_{\rm
l}}}\eta_{\rm l}^{\frac{1}{\alpha_{\rm l}}}r^{\frac{\alpha_{\rm
t}}{\alpha_{\rm l}}},h)$.

## References

  * [1] S. Sekander, H. Tabassum, and E. Hossain, “Multi-tier drone architecture for 5G/B5G cellular networks: Challenges, trends, and prospects,” _IEEE Communications Magazine_ , vol. 56, no. 3, pp. 96–103, 2018.
  * [2] B. Li, Z. Fei, and Y. Zhang, “UAV communications for 5G and beyond: Recent advances and future trends,” _IEEE Internet of Things Journal_ , vol. 6, no. 2, pp. 2241–2263, 2018.
  * [3] Y. Zeng, R. Zhang, and T. J. Lim, “Wireless communications with unmanned aerial vehicles: Opportunities and challenges,” _IEEE Communications Magazine_ , vol. 54, no. 5, pp. 36–42, 2016.
  * [4] A. Fotouhi, H. Qiang, M. Ding, M. Hassan, L. G. Giordano, A. Garcia-Rodriguez, and J. Yuan, “Survey on UAV cellular communications: Practical aspects, standardization advancements, regulation, and security challenges,” _IEEE Communications Surveys & Tutorials_, vol. 21, no. 4, pp. 3417–3442, 2019.
  * [5] M. Alzenad and H. Yanikomeroglu, “Coverage and rate analysis for vertical heterogeneous networks (VHetNets),” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 12, pp. 5643–5657, Dec. 2019.
  * [6] N. Cheng, F. Lyu, W. Quan, C. Zhou, H. He, W. Shi, and X. Shen, “Space/aerial-assisted computing offloading for IoT applications: A learning-based approach,” _IEEE Journal on Selected Areas in Communications_ , vol. 37, no. 5, pp. 1117–1129, 2019.
  * [7] M. Mozaffari, W. Saad, M. Bennis, Y.-H. Nam, and M. Debbah, “A tutorial on UAVs for wireless networks: Applications, challenges, and open problems,” _IEEE communications surveys & tutorials_, vol. 21, no. 3, pp. 2334–2360, 2019.
  * [8] M. Kishk, A. Bader, and M.-S. Alouini, “Aerial base station deployment in 6G cellular networks using tethered drones: The mobility and endurance tradeoff,” _IEEE Vehicular Technology Magazine_ , vol. 15, no. 4, pp. 103–111, 2020.
  * [9] M. Matracia, M. Kishk, and M.-S. Alouini, “Coverage analysis for UAV-assisted cellular networks in rural areas,” _IEEE Open Journal of Vehicular Technology_ , pp. 1–1, 2021.
  * [10] C. Saha, M. Afshang, and H. S. Dhillon, “Enriched $k$-tier HetNet model to enable the analysis of user-centric small cell deployments,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 3, pp. 1593–1608, March 2017.
  * [11] M. Afshang, C. Saha, and H. S. Dhillon, “Nearest-neighbor and contact distance distributions for Thomas cluster process,” _IEEE Wireless Communications Letters_ , vol. 6, no. 1, pp. 130–133, 2017.
  * [12] M. Peng, C. Wang, V. Lau, and H. V. Poor, “Fronthaul-constrained cloud radio access networks: Insights and challenges,” _IEEE Wireless Communications_ , vol. 22, no. 2, pp. 152–160, 2015.
  * [13] Y. S. Soh, T. Q. Quek, M. Kountouris, and H. Shin, “Energy efficient heterogeneous cellular networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 31, no. 5, pp. 840–850, 2013.
  * [14] D. Cao, S. Zhou, and Z. Niu, “Optimal combination of base station densities for energy-efficient two-tier heterogeneous cellular networks,” _IEEE Transactions on Wireless Communications_ , vol. 12, no. 9, pp. 4350–4362, 2013\.
  * [15] M. Wildemeersch, T. Q. S. Quek, C. H. Slump, and A. Rabbachin, “Cognitive small cell networks: Energy efficiency and trade-offs,” _IEEE Transactions on Communications_ , vol. 61, no. 9, pp. 4016–4029, 2013.
  * [16] S. Shakoor, Z. Kaleem, M. I. Baig, O. Chughtai, T. Q. Duong, and L. D. Nguyen, “Role of UAVs in public safety communications: Energy efficiency perspective,” _IEEE Access_ , vol. 7, pp. 140 665–140 679, 2019.
  * [17] Z. Yang, W. Xu, and M. Shikh-Bahaei, “Energy efficient uav communication with energy harvesting,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 2, pp. 1913–1927, 2020.
  * [18] Y. Zeng, R. Zhang, and T. J. Lim, “Throughput maximization for UAV-enabled mobile relaying systems,” _IEEE Transactions on Communications_ , vol. 64, no. 12, pp. 4983–4996, 2016.
  * [19] Y. Zeng and R. Zhang, “Energy-efficient UAV communication with trajectory optimization,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 6, pp. 3747–3760, 2017.
  * [20] A. Fotouhi, M. Ding, and M. Hassan, “Dynamic base station repositioning to improve spectral efficiency of drone small cells,” in _2017 IEEE 18th International Symposium on A World of Wireless, Mobile and Multimedia Networks (WoWMoM)_ , 2017, pp. 1–9.
  * [21] M. Gruber, “Role of altitude when exploring optimal placement of UAV access points,” in _2016 IEEE Wireless Communications and Networking Conference_ , 2016, pp. 1–5.
  * [22] D. H. Choi, S. H. Kim, and D. K. Sung, “Energy-efficient maneuvering and communication of a single UAV-based relay,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 50, no. 3, pp. 2320–2327, 2014.
  * [23] L. Li, M. Peng, C. Yang, and Y. Wu, “Optimization of base-station density for high energy-efficient cellular networks with sleeping strategies,” _IEEE Transactions on Vehicular Technology_ , vol. 65, no. 9, pp. 7501–7514, 2016.
  * [24] H. ElSawy, A. Sultan-Salem, M.-S. Alouini, and M. Z. Win, “Modeling and analysis of cellular networks using stochastic geometry: A tutorial,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 167–203, 2016.
  * [25] H. ElSawy, E. Hossain, and M. Haenggi, “Stochastic geometry for modeling, analysis, and design of multi-tier and cognitive cellular wireless networks: A survey,” _IEEE Communications Surveys Tutorials_ , vol. 15, no. 3, pp. 996–1019, 2013.
  * [26] B. Galkin, J. Kibida, and L. A. DaSilva, “A stochastic model for UAV networks positioned above demand hotspots in urban environments,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 7, pp. 6985–6996, 2019.
  * [27] Y. Qin, M. A. Kishk, and M.-S. Alouini, “On the influence of charging stations spatial distribution on aerial wireless networks,” _IEEE Transactions on Green Communications and Networking_ , vol. 5, no. 3, pp. 1395–1409, 2021.
  * [28] ——, “Performance evaluation of UAV-enabled cellular networks with battery-limited drones,” _IEEE Communications Letters_ , vol. 24, no. 12, pp. 2664–2668, 2020.
  * [29] ——, “Drone charging stations deployment in rural areas for better wireless coverage: Challenges and solutions,” _IEEE Internet of Things Magazine_ , vol. 5, no. 1, pp. 148–153, 2022.
  * [30] O. M. Bushnaq, M. A. Kishk, A. Celik, M.-S. Alouini, and T. Y. Al-Naffouri, “Optimal deployment of tethered drones for maximum cellular coverage in user clusters,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 3, pp. 2092–2108, 2021.
  * [31] M. A. Kishk, A. Bader, and M.-S. Alouini, “On the 3-D placement of airborne base stations using tethered UAVs,” _IEEE Transactions on Communications_ , vol. 68, no. 8, pp. 5202–5215, 2020.
  * [32] H. S. Dhillon, R. K. Ganti, F. Baccelli, and J. G. Andrews, “Modeling and analysis of k-tier downlink heterogeneous cellular networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 30, no. 3, pp. 550–560, 2012\.
  * [33] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal LAP altitude for maximum coverage,” _IEEE Wireless Communications Letters_ , vol. 3, no. 6, pp. 569–572, Dec 2014.
  * [34] Y. Zeng, J. Xu, and R. Zhang, “Energy minimization for wireless communication with rotary-wing UAV,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 4, pp. 2329–2345, April 2019.
  * [35] M. Di Renzo, A. Zappone, T. T. Lam, and M. Debbah, “System-level modeling and optimization of the energy efficiency in cellular networks - a stochastic geometry framework,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 4, pp. 2539–2556, 2018.