# Downlink Interference Analysis of UAV-based mmWave Fronthaul for Small Cell
Networks

Mohammad Taghi Dabiri, Mazen Hasna, Senior Member, IEEE, and Walid Saad,
Fellow, IEEE Mohammad Taghi Dabiri and Mazen Hasna are with the Department of
Electrical Engineering, Qatar University, Doha, Qatar (E-mail:
<EMAIL_ADDRESS>and hasna@qu.edu.qa).Walid Saad is with the Wireless$@$VT,
Bradley Department of Electrical and Computer Engineering, Virginia Tech,
Arlington, VA 22203 USA.

###### Abstract

In this paper, an unmanned aerial vehicles (UAV)-based heterogeneous network
is studied to solve the problem of transferring massive traffic of distributed
small cells to the core network. First, a detailed three-dimensional (3D)
model of the downlink channel is characterized by taking into account the real
antenna pattern, UAVs’ vibrations, random distribution of small cell base
stations (SBSs), and the position of UAVs in 3D space. Then, a rigorous
analysis of interference is performed for two types pf interference: intra-
cell interference and inter-cell interference. The interference analysis
results are then used to derive an upper bound of outage probability on the
considered system. Using numerical results show that the analytical and
simulation results match one another. The results show that, in the presence
of UAV’s fluctuations, optimizing radiation pattern shape requires balancing
an inherent tradeoff between increasing pattern gain to reduce the
interference as well as to compensate large path loss at mmWave frequencies
and decreasing it to alleviate the adverse effect of a UAV’s vibrations. The
analytical derivations enable the derivation of the optimal antenna pattern
for any condition in a short time instead of using time-consuming extensive
simulations.

###### Index Terms:

Antenna pattern, backhaul/fronthaul links, interference, mmWave communication,
small cell networks, unmanned aerial vehicles (UAVs).

##  I Introduction

Driven by the development of mobile Internet and smart phones, the need for
data traffic increases by a factor of ten every five years and with the
emerging Internet of Things (IoT) predicted to wirelessly connect trillions of
users/devices in ultra-dense networks. [1, 2, 3, 4, 5]. With the dense
deployment of small cell base stations (SBSs) which are vital for the ongoing
network densification, high capacity backhaul/fronthaul links are required
with a low latency of around 100 $\micro$s or less [6]. For the ongoing
network densification, use of dense deployment of small cell base stations
(SBSs) is vital which requires high capacity backhaul/fronthaul links [6].
These demands can be traditionally fulfilled by coaxial cable or optical fiber
in terrestrial networks. However, for massive deployment of small cells, such
solutions will not be flexible and easy to deploy, and they won’t be as cost-
effective as compared to wireless backhaul/fronthaul links [7]. Wireless
backhaul/fronthaul connectivity can be realized using microwave bands for non-
line-of-sight (NLoS) cases or millimeter wave (mmWave) and free-space optical
(FSO) links for line-of-sight (LoS) cases [8]. Microwave links can cover a
wide area but suffer from low data rates. FSO and mmWave-based
backhaul/fronthaul links can meet the capacity requirements of next generation
communication networks, and are lightweight and easy to install. However,
mmWave/FSO links suffer from susceptibility to weather conditions [9] and
require a LoS connection, which is a major challenge in urban regions due to
the presence of small cells in hard-to-reach, near-street-level locations.
Recently, a promising alternative solution was presented in [9] that utilizes
unmanned aerial vehicle (UAVs) as a wireless backhaul/fronthaul hub point
between small cells and the core network. These UAV-hubs acting as networked
flying platforms (NFPs) provide a possibility of wireless LoS fronthaul link.
One of the important challenges related to the considered NFP-based topology
is the design of the fronthaul downlink using directional mmWave antennas
which is the subject of this work. Although there are several works in this
field [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28], which will be reviewed in detail below, a comprehensive analysis of the
considered downlink has not been provided in any of these studies by taking
into account the effects of the real antenna pattern, the UAVs’ fluctuations,
the random distribution of SBSs, and the effects of inter-cell and intra-cell
interference caused by side-lobes of real antenna patterns. For example, in
most of these studies [14, 15, 16, 17, 18, 19, 20, 21, 22, 23], the effect of
the UAV’s vibrations has been neglected, while we will show that by increasing
the gain of the millimeter wave antenna, small fluctuations of less than one
degree also affect the performance of the considered system. Moreover, in most
of these studies, approximate patterns with fixed gain are used, which cannot
obtain the exact results of the performance.

###  I-A Literature Review and Motivation

Although UAV-assisted channel modeling was investigated in recent studies [29,
30, 31, 32], these works are mainly limited to sub-6 GHz frequency bands which
cannot be directly employed to UAV-based mmWave communication systems.
Meanwhile, most of the prior studies on mmWave channel modeling [10, 11, 12]
do not address the presence of UAVs, with the exception of a few recent works
in [13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. For
instance, a new LoS probability model for UAV-to-ground base station (BS)
links is proposed in [14] over realistic urban grid deployments in which the
effects of building height distributions along with their densities and
dimensions are considered. In [15] and [16], the authors provided a new
channel characterization for UAV-based mmWave links by using ray tracing
simulations at two different frequency bands: 60 GHz and 28 GHz. In [17], a
novel ground-to-UAV channel prediction method is provided by using ray tracing
simulations based on the minimum Euclidean distance. However, the results of
these works [14, 15, 16, 17] are limited for omnidirectional mmWave antenna
pattern. High gain antennas can point the signals in a preferred direction in
order to efficiently enhance the received signal power [33]. Moreover, the
narrow beamwidth of a directional antenna concentrates the signal power at the
target point, which results in less interference and higher security in the
physical layer.

Due to these advantages, the use of directional antennas for UAV-based mmWave
communications was studied in [18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28]. A
dynamic algorithm is presented in [18] to get a better performance by
considering dynamic backhauling and mobility of users. A UAV-assisted
configuration is introduced in [19, 20, 21] wherein a non-orthogonal multiple
access technique was employed to provide wireless connectivity for a large
number of users. The use of aerial relay nodes was investigated in [22] for
both wireless access and backhaul mmWave links to allow dynamic routing. A
sectoring method is proposed in [23] to guarantee ground user downlink
coverage using a mmWave-enabled directional UAV link. However, the results of
the works in [18, 19, 20, 21, 22, 23] assume that the UAVs are highly stable
and not affected by vibrations. However, in practice, as the antenna gain
increases, the sensitivity to the small UAV’s vibration increases. This leads
to an unreliable communication due to antenna gain mismatch [34, 35, 36].

More recently, in [25, 26, 27], the authors studied the relationship between
UAV’s vibrations and mmWave antenna gain. The results of these work clearly
show that the performance of a UAV-based link with directional mmWave antenna
is highly dependent on the strength of UAV’s vibrations. However, the results
of [25] and [26] are limited to a special case that does not account for the
effects of interference caused by the side lobe gain of antenna array. More
recently, a UAV-assisted mmWave channel model was proposed in [28] for a UAV-
based network by taking into consideration real parameters such as UAV’s
vibrations and interference caused by antenna side lobes. Generally, the
dominant network traffic is for downloading high-quality videos and watching
live sports matches or concerts. Therefore, the downlink capacity is often
several times greater than that of the uplink and thus, its design is more
important. However, the results of [28] are limited to the uplink and cannot
be directly used for the design and analysis of a downlink. Due to the
asymmetric nature of UAV-based network, as we show in this study, the
interference model of the uplink and downlink will be radically different. In
addition, in [28], the questions regarding the effect of random distribution
on SBSs have not been answered.

###  I-B Contributions

The main contribution of this paper is a holistic interference analysis of
fronthaul downlink transmissions for a UAV-based network in the presence of
the realistic channel parameters. To the best of the authors’ knowledge, there
exists no prior work on UAV-based mmWave downlink communications that jointly
takes into account the actual antenna pattern, the UAVs’ vibrations, and the
random effect of the SBSs. In summary, the paper makes the following key
contributions:

  * •

First, we characterize a detailed 3D modeling of the downlink channel, by
taking into account the real antenna pattern, UAVs’ vibrations, random
distribution of SBSs, and the position of NFPs in 3D space. Using this
characterization, we derive the end-to-end signal-to-noise plus interference
ratio (SINR) by considering all the mentioned parameters.

  * •

We then rigorously analyze the interference in the system. Interference is
divided into two important categories: intra-cell interference and inter-cell
interference. We model both types of interference by considering all channel
parameters. Then, using simulations, the effect of the channel parameters on
the interference is studied and important results are extracted regarding the
reduction of the interference effect.

  * •

Subsequently, we use the obtained results about interference in order to
examine the system performance in terms of outage probability. To this end,
first, the effect of the random distribution of SBSs on the outage probability
is studied. We will show that the important parameter affecting the
interference is the spatial angle between the mmWave directional antennas,
which is a function of the random distribution of the SBSs. Therefore, we
derive the probability density function (PDF) and cumulative distribution
function (CDF) of spatial angle between the mmWave directional antennas. Using
those results, we derive an upper bound on the average outage probability of
the considered system by taking into account all the channel parameters.

  * •

Finally, using extensive simulation and numerical results, we evaluate the
performance of the considered system in terms of outage probability. Our
simulation results show that the derived analytical expression is close to the
simulated system with acceptable accuracy. We also show that in the presence
of UAV’s fluctuations, it is very necessary to choose the optimal antenna
pattern to manage the interference and at the same time to reduce the effect
of UAV’s fluctuations. Our analytical derivations allow us to obtain the
optimal antenna pattern for any condition in a short time instead of using
time-consuming extensive simulations.

The rest of this paper is organized as follows. We characterize the considered
downlink channel in Section II. Interference modeling and analysis are done in
Section III, and outage probability analysis along with related simulations
are provided in Section IV. Finally, conclusions are drawn in Section V.

(a)

(b)

(c)

Figure 1: (a) An illustration of a UAV-assisted HetNet as an alternative solution for backhaul/fronthaul links which uses UAVs and higher directional frequencies (such as mmWave, THz, and optical frequencies) to transfer traffic from the distributed SBSs to the core network. In this topology, each UAV is able to establish fronthaul links for several SBSs. (b) A 3D illustration of cell $i$ for which several SBSs (denoted by $\mathcal{S}_{i,j}$ for $j\in\\{1,...,J\\}$) are connected to $U_{i}$. The spatial angle between the downlink transmissions of $\mathcal{S}_{ij}$ and $\mathcal{S}_{i1}$ is denoted by $\theta_{ij}$. Here, higher values for $\theta_{ij}$ or lower values for antenna beamwidth lead to a lower interference between $\mathcal{S}_{ij}$ and $\mathcal{S}_{i1}$. The position of each $\mathcal{S}_{ij}$ relative to the $U_{i}$ is determined using angles $\theta_{x_{ij}}$ and $\theta_{y_{ij}}$. TABLE I: The list of main notations. Parameter | Description
---|---
Subscript $i$ | Denote the cell $i$
Subscript $j$ | Denote the SBS $j$ in each cell
Subscript $q$ |  $\in\\{s,u\\}$ determines the SBS and UAV nodes
Subscript $w$ |  $\in\\{x,y\\}$ determines $x$ and $y$ axes
$\sigma_{\theta}$ | Characterize the strength of instability of UAVs
$J$ | Number of SBSs in each cell
$I$ | Number of UAVs
$\alpha_{c}$ | Defined in Remark 5
$U_{i}$ | Denote $i$th NFP where $i\in\\{1,...,N_{D}\\}$
$\mathcal{S}_{ij}$ | Denote $j$th SBS in the $i$th cell
$T_{ij}$ | Denote the link between $U_{i}$ and $\mathcal{S}_{ij}$
$A_{\mathcal{S}_{ij}}$ | Denote the array antenna of $\mathcal{S}_{ij}$
$A_{u_{i}}$ | Denote the array antenna of $U_{i}$
$P_{t}$ | Transmitted power
$L_{ij}$ | Link length of $T_{ij}$ link
$h_{L_{ij}}$ | Path loss of a $T_{ij}$ link
$N_{u}$ | Denote square array antenna of $U_{i}$ with
|  $N_{u}\times N_{u}$ elements
$N_{s}$ | Denote square array antenna of $\mathcal{S}_{ij}$ with
|  $N_{s}\times N_{s}$ elements
$h_{i}$ | The height of $i$th NFP
$[x;y;z]$ | Cartesian coordinate system that the position of each $\mathcal{S}_{ij}$
| and $U_{1}$ are characterized as $[x_{ij},y_{ij},0]$ and $[0,0,h_{1}]$
$\lambda$ and $f_{c}$ | Wavelength and carrier frequency
$\alpha$ nad $\beta$ | constants whose values depend on the
| propagation environment
$\sigma^{2}_{N}$ | The thermal noise power
$\Theta_{i}$ |  $=[\Theta_{x_{i}},\Theta_{y_{i}}]$ that indicate the orientation fluctuations
| of $U_{i}$ in $x-z$, and $y-z$ planes. See Fig. 2
$\theta_{w_{ij}}$ | Defined in (9) and Fig. 1c where $w\in\\{x,y\\}$
$\theta_{ij}$ | It is calculated from $\theta_{x_{ij}}$ and $\theta_{y_{ij}}$, similar to (II-C)
$\theta_{wd_{1j}}$ |  $=(\theta_{w_{11}}-\theta_{w_{1j}})$. See Fig. 2.
$\theta_{d_{1j}}$ | It is calculated from $\theta_{xd_{1j}}$ and $\theta_{yd_{1j}}$, similar to (II-C)
$\theta_{\textrm{elev},i}$ | Elevation angle of $\mathcal{S}_{11}$ compared to $U_{i}$
$f_{x}(x)$ | The PDF of RV $x$
$F_{x}(x)$ | The CDF of RV $x$

##  II The System Model

As shown in Fig. 1a, we consider a UAV-assisted HetNet that is offered as an
affordable and easy to deploy option in [9] to solve the concern related to
transferring traffic of the distributed SBSs to the core network. More
precisely, as an alternative for backhaul/fronthaul links, the UAV-based
mmWave or FSO links are suggested in [9] to transfer the SBSs traffic,
particularly in ultra-dense urban areas. However, under cloudy, foggy, and
raining conditions, a FSO communication link experiences high attenuation and
thus, it becomes difficult to guarantee the required backhaul/fronthaul
quality-of-service (QoS) [9]. Unlike FSO links, mmWave links are more robust
to foggy weather. Therefore, we use mmWave for backhaul/fronthaul links of the
considered system. As seen in Fig. 1a, the considered HetNet consists of
several different wireless backhaul/fronthaul links: i) SBS-to-NFP (S2N)
links, ii) inter NFP links, and iii) inter mother NFP links. Our focus will be
on the modeling and design of the S2N links. For an S2N fronthaul link, the
uplink consists of a transmission from a ground remote radio head (RRH) to an
NFP, and the downlink transmission will be essentially a link from an NFP to
an RRH. Here, we focus on the downlink.

###  II-A The 3D Antenna Pattern

In order to reduce the effect of interference and also to reduce the negative
effect of channel attenuation at higher frequencies, the use of high gain
antenna is essential, particularly for ultra high data rate backhaul/fronthaul
links. Advances in the fabrication of antenna array technology at higher
frequency allow the creation of large antenna arrays with high gain in a cost
effective and compact form [37, 38]. In addition, as we will show later,
employing a directional mmWave antenna pattern allows us to reuse frequency
bands thus improving the spectral efficiency of the considered system.

We assume that the antennas installed on the NFPs have the same pattern and
the antennas installed on the SBSs have the same pattern as well. We consider
a uniform square array antennas for both NFPs and SBSs, comprising
$N_{q}\times N_{q}$ antenna elements with the same spacing $d_{a}$ between
elements, where the subscript $q=s$ determines the SBS node and the subscript
$q=u$ determines UAV node. The array radiation gain is mainly formulated in
the direction of $\theta$ and $\phi$. By taking into account the effect of all
elements, the array radiation gain will be:

$\displaystyle
G_{q}(N_{q},\theta,\phi)=G_{0}(N_{q})\,\underbrace{G_{e}(\theta,\phi)\,G_{a}(N_{q},\theta,\phi)}_{G^{\prime}_{q}(N_{q},\theta,\phi)},$
(1)

where $G_{a}$ is an array factor, $G_{e}$ is single element radiation pattern
and $G_{0}$ is defined in (8). From 3GPP single element radiation pattern,
$G_{e,\textrm{3dB}}=10\times\log_{10}(G_{e})$ of each single antenna element
will be given by [39]

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&G_{e\textrm{3dB}}=G_{\textrm{max}}-\min\left\\{-(G_{e\textrm{3dB,1}}+G_{e\textrm{3dB,2}}),F_{m}\right\\},\\\
&G_{e\textrm{3dB,1}}=-\min\left\\{-12\left(\frac{\theta_{e}-90}{\theta_{e\textrm{3dB}}}\right)^{2},G_{\textrm{SL}}\right\\},\\\
&G_{e\textrm{3dB,2}}=-\min\left\\{-12\left(\frac{\theta_{x}}{\phi_{e\textrm{3dB}}}\right)^{2},F_{m}\right\\},\\\
&\theta_{e}=\tan^{-1}\left(\frac{\sqrt{1+\sin^{2}(\theta_{x})}}{\sin(\theta_{y^{\prime}})}\right),\end{array}\right.$
(6)

where $\theta_{e\textrm{3dB}}=65^{\circ}$ and
$\phi_{e\textrm{3dB}}=65^{\circ}$ are the vertical and horizontal 3D
beamwidths, respectively, $G_{\textrm{max}}$ is the maximum directional gain
of the antenna element, $F_{m}=30$ dB is the front-back ratio, and
$G_{\textrm{SL}}=30$ dB is the side-lobe level limit. If the amplitude
excitation of the entire array is uniform, then the array factor
$G_{a}(N_{q},\theta,\phi)$ for a square array of $N_{q}\times N_{q}$ elements
can be obtained as [40, eqs. (6.89) and (6.91)]:

$\displaystyle G_{a}(N_{q},\theta,\phi)$
$\displaystyle=\left(\frac{\sin\left(\frac{N_{q}(kd_{a}\sin(\theta)\cos(\phi)+\beta_{x})}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin(\theta)\cos(\phi)+\beta_{x}}{2}\right)}\right.$
$\displaystyle\times\left.\frac{\sin\left(\frac{N_{q}(kd_{a}\sin(\theta)\sin(\phi)+\beta_{y})}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin(\theta)\sin(\phi)+\beta_{y}}{2}\right)}\right)^{2},$
(7)

where $d_{a}=\frac{\lambda}{2}$ and $\beta_{w}$ are the spacing and
progressive phase shift between the elements, respectively.
$k=\frac{2\pi}{\lambda}$ is the wave number, $\lambda=\frac{c}{f_{c}}$ is the
wavelength, $f_{c}$ is the carrier frequency, and $c$ is the speed of light.
Also, in order to guarantee that the total radiated power of antennas with
different $N$ are the same, the coefficient $G_{0}$ is defined as

$\displaystyle
G_{0}(N_{q})=\frac{4\pi}{\int_{0}^{\pi}\int_{0}^{2\pi}G^{\prime}_{q}(\theta,\phi)\sin(\theta)d\theta
d\phi}.$ (8)

Based on (II-A), the maximum value of the antenna gain is equal to
$G_{0}(N_{q})$, which is obtained when $\theta=0$.

###  II-B Topology Description

Let $U_{i}$ represent the $i$th NFP, and let $\mathcal{S}_{i,j}$ be the set of
SBSs connected to $U_{i}$. Also, $A_{\mathcal{S}_{ij}}$ indicates the array
antenna of $\mathcal{S}_{i,j}$, which is directed towards the $U_{i}$. As
shown in Fig. 1b, the set of $\mathcal{S}_{i,j}$ connected to $U_{i}$ forms
the cell $i$. We assume that the SBSs are distributed in the $x-y$ plane and
the $U_{i}$ are flying at height $h_{i}$. Therefore, the position of each
$\mathcal{S}_{ij}$ is characterized as $[x_{ij},y_{ij},0]$ in a Cartesian
coordinate system where $[x,y,z]\in\mathbb{R}^{1\times 3}$. As shown in Fig.
1c, the position of each $\mathcal{S}_{ij}$ relative to the $U_{i}$ is
uniquely determined using angles $\theta_{x_{ij}}$ and $\theta_{y_{ij}}$ as
follows:

$\displaystyle
x_{ij}=h_{i}\tan(\theta_{x_{ij}}),~{}~{}~{}y_{ij}=h_{i}\tan(\theta_{y_{ij}}).$
(9)

Each $U_{i}$ is equipped with $J$ directional mmWave array antennas, which are
denoted by $A_{u_{ij}}$. The direction of $A_{u_{ij}}$ can be electrically or
using a simple stepper motor, aligned towards $A_{\mathcal{S}_{ij}}$ [41].
Therefore, according to (9), UAV $U_{i}$ sets the direction of each of the
antennas on the angle $\theta_{ij}=[\theta_{x_{ij}},\theta_{y_{ij}}]$. We
define $T_{ij}$ as the communication link between $A_{\mathcal{S}_{ij}}$ and
$A_{u_{ij}}$.

In order to analyze the considered system, we need a comprehensive channel
model in which the effect of the spatial distribution of SBSs, the spatial
distribution of NFPs, antenna patterns, and UAVs’ instabilities are
considered. For notational simplicity, we obtain the SINR for link $T_{11}$ in
general, which can be easily used to represent the rest of the links. Note
that we have two types of interference, one is intra-cell interference due to
$T_{1j}$ links and the other is inter-cell interference due to $T_{ij}$ links
where $i\neq 1$. Therefore, the SINR at the receiver (Rx) of
$\mathcal{S}_{11}$ will be:

$\displaystyle\gamma_{11}=\frac{P_{r_{11}}}{\underbrace{\sum_{j=2}^{J}\mathbb{I}_{1j}}_{\text{intra-
cell}}+\underbrace{\sum_{i=2}^{I}\sum_{j=1}^{J}\mathbb{I}_{ij}}_{\text{inter-
cell interference}}+\sigma^{2}_{N}}$ (10)

where $\sigma^{2}_{N}$ is the thermal noise power. $P_{r_{11}}$ is the
received power transmitted by $A_{\mathcal{S}_{11}}$ given by:

$\displaystyle
P_{r_{11}}=P_{t}|h_{L_{11}}|^{2}\underbrace{G_{u_{11}}(N_{u_{11}},\Theta_{u_{11}},\phi_{u_{11}})}_{\text{Pattern
Gain of
$A_{u_{ij}}$}}\underbrace{G_{\mathcal{S}_{11}}(N_{\mathcal{S}_{11}},\Theta_{\mathcal{S}_{11}},\phi_{\mathcal{S}_{11}})}_{\text{Pattern
Gain of $A_{\mathcal{S}_{11}}$}}.$ (11)

In (11), $P_{t}$ is the transmit power, $h_{L_{ij}}=h_{Lf_{ij}}h_{Lm_{ij}}$ is
the channel path loss, $h_{Lf_{ij}}=\left(\frac{\lambda}{4\pi
L_{ij}}\right)^{2}$ is the free-space path loss,
$h_{Lm_{ij}}=e^{-\frac{\mathcal{K}(f)}{2}L_{ij}}$ represents the molecular
absorption loss where $\mathcal{K}(f)$ is the frequency dependent absorption
coefficient, and $L_{ij}$ is the link length of $T_{ij}$. Also, random
variables (RVs) $\Theta_{u_{ij}}$ and $\Theta_{\mathcal{S}_{ij}}$ capture the
orientation fluctuations of $A_{u_{ij}}$ and $A_{\mathcal{S}_{ij}}$, and they
are modeled in the next section. For simplicity, the simplified symbol
$\Theta_{i}=[\Theta_{x_{i}},\Theta_{y_{i}}]$ is used to refer to RV
$\Theta_{u_{ij}}$.

###  II-C Modeling UAV Vibrations

In an ideal scenario, the directional antennas are perfectly aligned and the
parameters $\Theta_{q_{ij}}$ in (11) tend to zero and the received signal
power tends to the maximum value as

$\displaystyle P_{r_{11}}\simeq
P_{t}|h_{L_{11}}|^{2}G_{0}(N_{u})G_{0}(N_{s}).$ (12)

Based on the ideal results in (12), a higher value of $N_{q}$, leads to a
higher value for the received power. Also, in Section III, we will show that a
larger $N_{q}$ can lead to a lower interference, and hence, the SINR always
improves by increasing $N_{q}$. However, in real situations, the SINR does not
necessarily improve by increasing $N_{q}$. In practical situations, an error
in mechanical control system, mechanical noise, position estimation errors,
air pressure, and wind speed can affect the UAV’s angular and position
stability. In practice, the instantaneous orientation of a UAV can randomly
deviate from its means denoted by $\Theta_{q}$. This, in turn, leads to
deviations in the angle of departure (AoD) of the transmitter (Tx)’s antenna
pattern. Therefore, as $N_{q}$ increases, the system becomes more sensitive to
the UAV’s vibrations. Therefore, on the one hand, we must increase $N_{q}$ to
reduce interference as well as to increase the desired received power. On the
other hand, we must be careful not to choose a very large value for $N_{q}$
such that even with small fluctuations in antenna direction, the probability
of missing the main lobe and being on the side-lobes increases.

The orientation fluctuations of antennas mounted on $U_{i}$ are captured by
$\Theta_{i}=[\Theta_{x_{i}},\Theta_{y_{i}}]$ where $\Theta_{x_{i}}$ and
$\Theta_{y_{i}}$ are the instantaneous fluctuations of $U_{i}$ in the $x-z$
and $y-z$ Cartesian coordinates, respectively. Therefore, the overall
instantaneous deviation angle of $A_{u_{ij}}$ is obtained as:

$\displaystyle\Theta_{i}$
$\displaystyle=\tan^{-1}\left(\sqrt{\tan^{2}(\Theta_{x_{i}})+\tan^{2}(\Theta_{y_{i}})}\right),$
$\displaystyle\phi_{i}$
$\displaystyle=\tan^{-1}\left(\frac{\tan(\Theta_{y_{i}})}{\tan(\Theta_{x_{i}})}\right).$
(13)

Based on the central limit theorem, the UAV’s orientation fluctuations are
considered to be Gaussian distributed [42, 43, 44]. Therefore, we have
$\Theta_{x_{i}}\sim\mathcal{N}(0,\sigma^{2}_{\theta_{x}})$, and
$\Theta_{y_{i}}\sim\mathcal{N}(0,\sigma^{2}_{\theta_{y}})$. For the ground
SBSs, the estimation error of the exact position of the flying NFPs and the
insufficient speed to track NFPs will cause an angular error. Unlike the
flying NFPs, we assume that the ground SBSs do not face weight and power
consumption limitations and, hence, they can align their antennas with the
considered NFP with a negligible angle error. Therefore, (11) can be rewritten
as follows:

$\displaystyle
P_{r_{11}}=P_{t}|h_{L_{11}}|^{2}G_{0}(N_{s})G_{u}(N_{u},\Theta_{1},\phi_{1}).$
(14)

###  II-D Probability of LoS

In addition to the high propagation loss, mmWave communication systems are
very sensitive to blockages [45]. Therefore, the probability of LoS is an
important factor and can be described as a function of the elevation angle and
environment as follows [46, 47]:

$\displaystyle
P_{\textrm{LoS}}(\theta_{\text{elev}})=\frac{1}{1+\alpha\exp\left(-\beta(\frac{180}{\pi}\theta_{\text{elev}}-\alpha)\right)}$
(15)

where $\alpha$ and $\beta$ are constants whose values depend on the
propagation environment, e.g., rural, urban, or dense urban, and
$\theta_{\textrm{elev}}$ is the elevation angle.

Figure 2: Graphical illustration of intra-cell interference in 2D $x-z$ plane
caused by antenna pattern of $A_{u_{1j}}$. The instantaneous orientation of
$U_{1}$ randomly deviates from its means $\Theta_{x_{1}}$ and $\Theta_{y_{1}}$
in the $x-z$ and $y-z$ planes.

(a)

(b)

Figure 3: (a) Real antenna pattern of a $20\times 20$ standard array antenna
versus $\theta$ for different values of $\phi$. Comparing the accuracy of the
model obtained from (42) with the real antenna pattern model for different
values of $\phi=0:0.33:45$.

##  III Interference Analysis

In this section, we first model interference and then investigate the effect
of interference under different conditions.

###  III-A Interference Modeling

For the considered system, the interference is divided into two categories,
intra-cell and inter-cell interference, as modeled next.

####  III-A1 Intra-cell Interference

For notation simplicity, here, we provide the general channel model for
$T_{11}$ link which connects $U_{1}$ to $\mathcal{S}_{11}$ in presence of
interference from other $\mathcal{S}_{ij}$ and $U_{i}$. The obtained model can
be easily used for other $T_{ij}$ links. For link $T_{11}$, the intra-cell
interference is caused by the interference from other links $T_{1j}$ for
$j\in\\{2,...,J\\}$. Therefore, as shown in Fig. 2, the intera-cell
interference will be:

$\displaystyle\mathbb{I}_{\text{intra}}=\sum_{j=2}^{J}\mathbb{I}_{1j},$ (16)

where $\mathbb{I}_{1j}$ is the interference caused by $T_{1j}$ link. As
graphically illustrated in Fig. 2, $\mathbb{I}_{1j}$ is obtained in (III-A1)
where we have $\theta_{xd_{1j}}=(\theta_{x_{11}}-\theta_{x_{1j}})$, and
$\theta_{yd_{1j}}=(\theta_{y_{11}}-\theta_{y_{1j}})$.

$\displaystyle\mathbb{I}_{1j}=$
$\displaystyle~{}P_{t}|h_{L_{11}}|^{2}G_{0}(N_{s})G_{0}(N_{u})G_{e}\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{1j}}-\Theta_{x_{1}})+\tan^{2}(\theta_{yd_{1j}}-\Theta_{y_{1}})}\right)\right)$
$\displaystyle\times\left(\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{1j}}-\Theta_{x_{1}})+\tan^{2}(\theta_{yd_{1j}}-\Theta_{y_{1}})}\right)\right)\cos(\phi_{u_{1j}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{1j}}-\Theta_{x_{1}})+\tan^{2}(\theta_{yd_{1j}}-\Theta_{y_{1}})}\right)\right)\cos(\phi_{u_{1j}})}{2}\right)}\right.$
$\displaystyle\times\left.\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{1j}}-\Theta_{x_{1}})+\tan^{2}(\theta_{yd_{1j}}-\Theta_{y_{1}})}\right)\right)\sin(\phi_{u_{1j}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{1j}}-\Theta_{x_{1}})+\tan^{2}(\theta_{yd_{1j}}-\Theta_{y_{1}})}\right)\right)\sin(\phi_{u_{1j}})}{2}\right)}\right)^{2}$
(17)

Also parameter $\phi_{u_{1j}}$ is the Roll angle of antenna $A_{u_{1j}}$
relative to $A_{\mathcal{S}_{11}}$, which is important in interference
analysis. For both intra-cell and inter-cell interference, the side-lobes
generally cause interference. Unlike the main lobe, as shown in Fig. 3a, the
side-lobes are functions of $\phi$.

Remark 1. Equation (III-A1) is a function of RVs $\Theta_{x_{1}}$,
$\Theta_{y_{1}}$, $\theta_{yd_{1j}}$, and $\theta_{xd_{1j}}$. However, the RVs
$\theta_{x_{1}}$, and $\theta_{y_{1}}$ are caused by the UAV’s vibrations and
have a coherence time of less than a second. Moreover, the parameters
$\theta_{yd_{1j}}$, and $\theta_{xd_{1j}}$ are caused by the distribution of
SBSs and NFPs in the 3D space, which usually change for a longer period of
time.

(a)

(b)

Figure 4: Characterization the interference caused by $T_{ij}$ link on the
considered $T_{11}$ link. This interference depends on two spatial angles
$\theta^{\prime}_{i,j1}=[\theta_{xd_{i,j1}}+\Theta_{x_{i}},\theta_{yd_{i,j1}}+\Theta_{y_{i}}]$,
$\theta^{\prime\prime}_{i1}=[\theta_{xd_{i1}},\theta_{yd_{i1}}]$, where
$\theta_{xd_{i,j1}}=(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}})$,
$\theta_{yd_{i,j1}}=(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}})$,
$\theta_{xd_{i1}}=(\theta^{\prime}_{x_{i1}}-\theta_{x_{11}})$, and
$\theta_{yd_{i1}}=(\theta^{\prime}_{y_{i1}}-\theta_{y_{11}})$. (a)
Characterization of parameter $\theta^{\prime}_{i,j1}$, and (b)
characterization of parameter $\theta^{\prime\prime}_{i1}$.

####  III-A2 Inter-cell Interference

The inter-cell interference is caused by the interference of $T_{ij}$ links on
the considered $T_{11}$ link when $i\neq 1$ and $j\in\\{1,...,J\\}$.
Therefore, the inter-cell interference will be given by:

$\displaystyle\mathbb{I}_{\text{inter}}=\sum_{i=2}^{I}\sum_{j=1}^{J}\mathbb{I}_{ij},$
(18)

where $I$ is the number of NFPs and $\mathbb{I}_{ij}$ is the interference
caused by a link $T_{ij}$. As shown in Fig. 4, $\mathbb{I}_{ij}$ depends on
the angle between the direction of $A_{u_{ij}}$ with respect to
$\mathcal{S}_{11}$ denoted by
$\theta^{\prime}_{i,j1}=[\theta_{xd_{i,j1}}+\Theta_{x_{i}},\theta_{yd_{i,j1}}+\Theta_{y_{i}}]$,
the angle between the direction of $A_{\mathcal{S}_{11}}$ with respect to
$U_{i}$ denoted by
$\theta^{\prime\prime}_{i,j}=[\theta_{xd_{i,j}},\theta_{yd_{i,j}}]$, and also
the distance between $U_{i}$ and $\mathcal{S}_{11}$ denoted by $L_{i1}$. The
position of $U_{1}$ is characterized as $[0,0,h_{1}]$ in the considered
Cartesian coordinate system $[x,y,z]\in\mathbb{R}^{1\times 3}$, and the
position of $U_{i}$ for $i\in\\{2,...,I\\}$ are characterized as
$[x_{i},y_{i},h_{i}]$. Therefore, $L_{i1}$ is obtained as

$\displaystyle
L_{i1}=\sqrt{(x_{i}^{2}+h_{1}\tan(\theta_{x_{11}}))^{2}+(y_{i}^{2}+h_{1}\tan(\theta_{y_{11}}))^{2}+h_{i}^{2}}.$
(19)

Also it can be easily shown that
$\theta_{xd_{i,j1}}=(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}})$,
$\theta_{yd_{i,j1}}=(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}})$,
$\theta_{xd_{i1}}=(\theta^{\prime}_{x_{i1}}-\theta_{x_{11}})$, and
$\theta_{yd_{i1}}=(\theta^{\prime}_{y_{i1}}-\theta_{y_{11}})$, where

$\displaystyle\theta^{\prime}_{x_{i1}}$
$\displaystyle=\tan^{-1}\left(\frac{h_{i}}{x_{i}+h_{1}\tan(\theta_{x_{11}})}\right)$
$\displaystyle\theta^{\prime}_{y_{i1}}$
$\displaystyle=\tan^{-1}\left(\frac{h_{i}}{y_{i}+h_{1}\tan(\theta_{y_{11}})}\right)$
(20)

Based on this, $\mathbb{I}_{ij}$ is defined in (III-A2).

$\displaystyle\mathbb{I}_{ij}=P_{t}|h_{L_{i1}}|^{2}G_{0}(N_{s})G_{0}(N_{u})$
$\displaystyle~{}~{}~{}~{}\times
G_{e}\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{i,j}})+\tan^{2}(\theta_{yd_{i,j}})}\right)\right)G_{e}\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}}+\Theta_{x_{i}})+\tan^{2}(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}}+\Theta_{y_{i}})}\right)\right)$
$\displaystyle~{}~{}~{}~{}\times\left(\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{i,j}})+\tan^{2}(\theta_{yd_{i,j}})}\right)\right)\cos(\phi_{\mathcal{S}_{11}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{i,j}})+\tan^{2}(\theta_{yd_{i,j}})}\right)\right)\cos(\phi_{\mathcal{S}_{11}})}{2}\right)}\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{i,j}})+\tan^{2}(\theta_{yd_{i,j}})}\right)\right)\sin(\phi_{\mathcal{S}_{11}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{xd_{i,j}})+\tan^{2}(\theta_{yd_{i,j}})}\right)\right)\sin(\phi_{\mathcal{S}_{11}})}{2}\right)}\right)^{2}$
$\displaystyle~{}~{}~{}~{}\times\left(\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}}+\Theta_{x_{i}})+\tan^{2}(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}}+\Theta_{y_{i}})}\right)\right)\cos(\phi_{u_{ij}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}}+\Theta_{x_{i}})+\tan^{2}(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}}+\Theta_{y_{i}})}\right)\right)\cos(\phi_{u_{ij}})}{2}\right)}\right.$
$\displaystyle~{}~{}~{}~{}\times\left.\frac{\sin\left(\frac{N_{q}(kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}}+\Theta_{x_{i}})+\tan^{2}(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}}+\Theta_{y_{i}})}\right)\right)\sin(\phi_{u_{ij}}))}{2}\right)}{N_{q}\sin\left(\frac{kd_{a}\sin\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta^{\prime}_{x_{i1}}-\theta_{x_{ij}}+\Theta_{x_{i}})+\tan^{2}(\theta^{\prime}_{y_{i1}}-\theta_{y_{ij}}+\Theta_{y_{i}})}\right)\right)\sin(\phi_{u_{ij}})}{2}\right)}\right)^{2}\times
P_{\text{LoS}}(\theta_{\text{elev},i}),$ (21) where $\displaystyle
P_{\textrm{LoS}}(\theta_{\text{elev},i})=\frac{1}{1+\alpha\exp\left(-\beta\left[\frac{180}{\pi}\tan^{-1}\left(\sqrt{\left(\frac{h_{i}}{x_{i}+h_{1}\tan(\theta_{x_{11}})}\right)^{2}+\left(\frac{h_{i}}{y_{i}+h_{1}\tan(\theta_{y_{11}})}\right)^{2}}\right)-\alpha\right]\right)}.$
(22)

Remark 2. In order to compute inter-cell interference from (III-A2), for any
distribution of SBSs and $U_{i}$, we calculate parameters
$\theta^{\prime}_{w_{ij}}$, $\theta_{w_{ij}}$, and $L_{i1}$ for
$w\in\\{x,y\\}$, $I\in\\{1,...,I\\}$, and $j\in\\{1,...,J\\}$ based on the
positions of SBSs and $U_{i}$ in the considered Cartesian coordinate system.
For any change in the distribution of SBSs and/or the position of hovering
$U_{i}$, we must recalculate and update the indicated parameters. However, the
parameters $\Theta_{x_{i}}$ and $\Theta_{y_{i}}$ are caused by the UAVs’
angular vibrations and have a coherence time of less than a second, and act as
rapid RVs in the analyses compared to the RVs $\theta^{\prime}_{w_{ij}}$,
$\theta_{w_{ij}}$, and $L_{i1}L_{i1}$.

By substituting (12), (III-A1) and (III-A2), in (10), the instantaneous value
of $\gamma_{11}$ is calculated, which is a complex expression of the channel
parameters. Before analytical derivations which is provided in Section IV, it
is necessary to have a preliminary view of the effect of the parameters on the
$\gamma_{11}$ and model it more simple. For this aim, in the next subsection,
by providing simulation, we will study the inter-cell and intra-cell
interference, which leads to important results. The obtained results will let
us to simplify $\gamma_{11}$ and prepare the groundwork for the comprehensive
analysis of the system in the next section.

TABLE II: Parameter values for simulations of single relay system. Parameters | Values | Parameters | Values
---|---|---|---
$P_{t}$ | 10 mW | $\gamma_{\text{th}}$ | 9 dB
$N_{s}$ | 15 | $N_{u}$ | 5-25
$f_{c}$ | 95 GHz | $J$ | 10
$R$ |  $400$ m | $T$ |  $20^{o}$C
BW | 3 GHz | $\sigma_{\theta}$ | $0.5^{o}-1.5^{o}$
$h_{i}$ | 50-100 m | $d_{x}=d_{y}$ | $\lambda/2$

###  III-B Simulation Results

The values used for the simulations are provided in Table II. In our
simulation setting, we randomly deploy (position of each SBS is a uniform
random variable on $x-y$ plane) 100 SBSs are randomly distributed in an area
on $x-y$ plane with a radius of 1.5 km. Then, ten UAVs are placed at a height
of 80-110 meter from the ground level in such a way that each UAV is connected
to the nearest ten SBSs. An SBS is chosen randomly and we call it
$\mathcal{S}_{11}$ and the UAV connected to it is called $U_{1}$. According to
the position of the other SBSs and $U_{i}$, we calculate the average value of
the parameters related to angle and distance. Then, $5\times 10^{7}$ random
independent values of RVs $\Theta_{x_{i}}$ and $\Theta_{y_{i}}$ are generated
related to the UAVs’ vibrations. After that, using (16) and (III-A1), $5\times
10^{7}$ independent values of intra-cell interference and using (18) and
(III-A2), $5\times 10^{7}$ independent values of inter-cell interference are
produced. Finally, by averaging over $5\times 10^{7}$ independent runs, the
average values of inter-cell and intra-cell interference are obtained.

The papremter $\phi$ is one of the important parameters for the interference
analysis. To show its importance, in Fig. 5, the inter-cell and intra-cell
interference are plotted for $N_{u}=15$. The results of Fig. 5 are provided
for three independent random distributions of SBSs. For each independent
distribution, similar to the steps mentioned above, the inter-cell and intra-
cell interference are plotted, separately, and indicated by “Run $n$” in Fig.
6 for $n\in\\{1,2,3\\}$. As can be observed from the results of Fig. 5, with
any change in the distribution of SBSs and NFPs, the inter-cell and intra-cell
interference changes. A very important point is that with the change of
$\phi_{\mathcal{S}_{11}}$, the inter-cell interference changes significantly.
Since $\phi_{\mathcal{S}_{11}}$ is a tunable parameter, for any distribution
of SBSs and NFPs, the parameter $\phi_{\mathcal{S}_{11}}$ can be adjusted in
such a way that the inter-cell interference falls below the noise level and
becomes almost negligible. For example, based on the results of Fig. 5, for
Run 1, the inter-cell interference falls bellow $10^{-15}$ for
$40^{o}<\phi_{\mathcal{S}_{11}}<45$, while for Run 2, the inter-cell
interference reaches less than $10^{-16}$ for
$\phi_{\mathcal{S}_{11}}=35^{o}$, which is much lower than the noise level
and/or intra-cell interference. However, changing $\phi_{\mathcal{S}_{11}}$
has almost no effect on the intra-cell interference. To find out the reason
for this, it should be noted that $\phi_{\mathcal{S}_{11}}$ is related to the
antenna of $\mathcal{S}_{11}$, which is directed towards $U_{1}$. Therefore,
the main lobe of $A_{\mathcal{S}_{11}}$ is towards $A_{u_{1j}}$, and as can be
observed from the results of Fig. 3, changing parameter
$\phi_{\mathcal{S}_{11}}$ has almost no effect on the main lobe antenna
pattern. As a result, the intra-cell interference is constant with respect to
the change of $\phi_{\mathcal{S}_{11}}$. Due to the importance of the obtained
results, in the following remarks, a summary of the obtained results is
provided, which will be used in the next section for the analysis and
parameter system design.

Remark 3. For any distribution of SBSs and NFPs, the tunable parameter
$\phi_{\mathcal{S}_{11}}$ can be set in such a way that the inter-cell
interference falls below the noise level and thus, inter-cell interference has
a negligible destructive effect on the considered system performance.

Remark 4. Since the antenna pattern in the main lobe is constant with respect
to $\phi$, the parameter $\phi_{\mathcal{S}_{11}}$ has almost no effect on the
intra-cell interference.

Figure 5: Comparison of inter- and intra-cell interference versus
$\phi_{\mathcal{S}_{11}}$ for three independent random distributions of SBSs.
Figure 6: Comparison of intra-cell interference due to $T_{1j}$ link versus
$\theta_{d_{1j}}$ for different values of $\phi_{u_{1j}}$.

Unlike inter-cell interference, intra-cell interference is the main factor in
reducing the quality of the considered system, and cannot be ignored due to
following reasons: (i) as discussed, the main lobe of $A_{\mathcal{S}_{11}}$
is set toward $A_{u_{11}}$ antenna mounted on $U_{1}$, and unintentionally, it
is also set toward the $A_{u_{1j}}$ antennas. (ii) The intra-cell interference
caused by $A_{u_{1j}}$ has a smaller distance $A_{\mathcal{S}_{11}}$ than the
inter-cell interference caused by $A_{u_{ij}}$, and thus, it has less
attenuation. (iii) Most importantly, compared to $A_{u_{ij}}$ antennas, with
more probability, $A_{u_{1j}}$ antennas are placed at a smaller spatial angle
than $A_{\mathcal{S}_{11}}$. For better understanding, according to (III-A2),
the intra-cell interference is caused by the interference between the main
lobe of $A_{\mathcal{S}_{11}}$ and the side-lobe of the $A_{u_{1j}}$ antennas.
For interference limited scenarios (i.e., SINR tends to signal-to-interference
ratio (SIR)), since the main lobe of $A_{\mathcal{S}_{11}}$ antenna
simultaneously boosts the main signal and intra-cell interference, changing
the antenna pattern of $A_{\mathcal{S}_{11}}$ has no effect on the SIR.
Therefore, the only option to reduce the intra-cell interference is to reduce
the interference caused by the side-lobes of $A_{u_{1j}}$, which is a function
of the antenna pattern as well as the spatial angle $\theta_{d_{1j}}$ and
$\phi_{u_{1j}}$a. In particular, the spatial angle $\theta_{d_{1j}}$ is a
function of the distribution of $\mathcal{S}_{1j}$ with respect to the
$U_{1}$, which is a random distribution. For better understanding, in Fig. 6,
the intra-cell interference is drawn versus $\theta_{d_{1j}}$ and different
values of $\phi_{u_{1j}}$. As we expected, by increasing $\theta_{d_{1j}}$,
the antenna side-lobes decreases and, as a result, the intra-cell interference
decreases. However, as can be seen, the interference does not decrease
linearly by increasing $\theta_{d_{1j}}$ and depends on the parameter
$\phi_{u_{1j}}$. An important point is that for small values of
$\theta_{d_{1j}}$, the intra-cell interference is approximately caused by the
main lobe of the antennas, and in this case, the communication link will
definitely be interrupted. This is related to the case when the distance of
one of the $\mathcal{S}_{1j}$ is very close to $\mathcal{S}_{11}$ (in the
order of several meters). In this case, as predicted in 5G, two close
frounthaul links should be integrated together, which is easily possible due
to their very small distance.

Remark 5. According to the above results, a minimum angle $\alpha_{c}$ is
defined for the $\mathcal{S}_{1j}$ that are close to each other. If the
spatial angle $\theta_{d_{1j}}$ between $T_{11}$ and $T_{1j}$ links is less
than $\alpha_{c}$, both frounthaul links should be integrated according to the
protocols defined in 5G.

##  IV Outage Probability Analysis

In this section, we use the results obtained in Section III, and we analyze
the performance of the considered system in terms of outage probability.

###  IV-A Outage Probability Analysis

Outage probability of $T_{11}$ link (i.e., the probability that $\gamma_{11}$
falls below a threshold $\gamma_{\text{th}}$) is obtained as follows:

$\displaystyle\mathbb{P}_{\text{out}}=\int_{0}^{\gamma_{\text{th}}}f_{\gamma_{11}}(\gamma_{11})\text{d}\gamma_{11}.$
(23)

Therefore, to analyze the considered system, the first step is to know the
distribution function of SINR $\gamma_{11}$. In the following theorem, the PDF
of $\gamma_{11}$ is provided.

Theorem 1. For any random distribution of SBSs, the PDF of $\gamma_{11}$ for
the worst case scenario is obtained as:

$\displaystyle
f_{\gamma_{11}|\theta_{d_{1j}}}(\gamma_{11})=\frac{1}{\sqrt{\beta_{qx}\beta_{qy}}}$
(24)
$\displaystyle\left[\sum_{j=1}^{J}\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}e^{-\frac{(\theta_{d_{1j}}-A^{\prime}_{m})^{2}}{w_{m}^{2}}}+N^{\prime}_{0}\right]^{\frac{\beta_{qx}+\beta_{qy}}{2\beta_{qx}\beta_{qy}}}h_{p}^{\frac{\beta_{qx}+\beta_{qy}}{2\beta_{qx}\beta_{qy}}-1}$
$\displaystyle
I_{0}\left(\beta_{qxy}\ln\left(\left[\sum_{j=1}^{J}\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}e^{-\frac{(\theta_{d_{1j}}-A^{\prime}_{m})^{2}}{w_{m}^{2}}}+N^{\prime}_{0}\right]\gamma_{11}\right)\right)$

where $\beta_{uw}=\frac{\sqrt{2}\sigma_{\theta_{w}}^{2}}{w_{0}}$,
$\beta_{qxy}=\frac{\beta_{qy}-\beta_{qx}}{2\beta_{qx}\beta_{qy}}$,
$N^{\prime}_{0}=\frac{\sigma^{2}_{N}}{P_{t}|h_{L_{11}}|^{2}G_{0}(N_{s})G_{0}(N_{u})}$
is the normalized thermal noise,

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&w_{m}=\frac{1}{N_{s}}~{}~{}\&~{}A^{\prime}_{m}=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}m=0\\\
&w_{m}=\frac{1}{2N_{s}}~{}\&~{}A^{\prime}_{m}=\sin^{-1}\left(\frac{2m+1}{N_{u}}\right),~{}~{}~{}~{}m\neq
0\\\ \end{array}\right.$ (27)
$\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&A^{\prime\prime}_{m}=1,~{}~{}\text{for}~{}~{}\frac{m}{N_{u}}\leq
1\\\ &A^{\prime\prime}_{m}=0,~{}~{}\text{for}~{}~{}\frac{m}{N_{u}}>1\\\
\end{array}\right.$

and $M=6$,
$A_{0}=1,A_{1}=0.05,A_{2}=0.020,A_{3}=0.011,A_{4}=0.006,A_{5}=0.004$,
$A_{6}=0.003$.

###### Proof:

Please refer to Appendix A. ∎

Notice that, as can be seen from (III-A1) and (III-A2), $\gamma_{11}$ is a
complex function of channel RVs, that include the location of SBSs relative to
each other as well as relative to the UAV connected to it, and the orientation
fluctuations of the directional antennas mounted on hovering UAVs. Therefore,
it seems impossible to obtain the exact distribution function for
$\gamma_{11}$. Using the results obtained in the Section III, and using some
derivations, the distribution function of $\gamma_{11}$ for the worst-case is
provided in Theorem 1. In network design, operators usually design the system
for the worst-case scenario or consider a guard band for analysis. In order to
derive the PDF of $\gamma_{11}$ in (24), we provide and use the approximated
antenna pattern in (42) as a function of the sum of several Gaussian
functions. As we see from Fig. 3b, the approximated model provided in (42)
models well the main lobe along with the maximum amount of side-lobes of real
antenna pattern. We use (42) to model the interference, and since the actual
value of the interference is equal to or less than the interference obtained
from (42), the interference obtained from (42) is for the worst scenario.

Theorem 2. For any random distribution of SBSs, the outage probability for the
worst case scenario is derived (IV-A).

$\displaystyle\leavevmode\resizebox{390.25534pt}{}{$+Q\left(\mathbb{B}_{2}\sqrt{\frac{-2\ln\left(\left[\sum_{j=1}^{J}\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}e^{-\frac{(\theta_{d_{1j}}-A^{\prime}_{m})^{2}}{w_{m}^{2}}}+N^{\prime}_{0}\right]\gamma_{th}\right)}{\beta_{ux}+\beta_{uy}}},\mathbb{B}_{1}\sqrt{\frac{-2\ln\left(\left[\sum_{j=1}^{J}\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}e^{-\frac{(\theta_{d_{1j}}-A^{\prime}_{m})^{2}}{w_{m}^{2}}}+N^{\prime}_{0}\right]\gamma_{th}\right)}{\beta_{ux}+\beta_{uy}}}\right)$}.$
(28)

In (IV-A), $Q(a,b)$ is the Marcum Q-function, $\gamma_{th}$ is the SINR
threshold, and

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\mathbb{B}_{1}=\frac{\sqrt{1-T_{q}^{4}}}{2T_{q}}\sqrt{\frac{1+T_{q}}{1-T_{q}}},\\\
&\mathbb{B}_{2}=\mathbb{B}_{1}\frac{1-T_{q}}{1+T_{q}},~{}~{}~{}T_{q}=\frac{\sigma_{\theta_{\text{min}}}}{\sigma_{\theta_{\text{max}}}},\\\
&\sigma_{\theta_{\text{max}}}=\text{max}\\{\sigma_{\theta_{x}},\sigma_{\theta_{y}}\\},~{}~{}~{}~{}\sigma_{\theta_{\text{min}}}=\text{min}\\{\sigma_{\theta_{x}},\sigma_{\theta_{y}}\\}\end{array}\right.$
(32)

###### Proof:

Please refer to Appendix B. ∎

The outage probability obtained in Theorem 2 is based on the distribution
function obtained in Theorem 1 and therefore calculates the upper bound of
outage probability. Another point is that the outage probability provided in
Theorem 2 is conditional on knowing $\theta_{d_{1j}}$. The parameter
$\theta_{d_{1j}}$ are function of the position of the $\mathcal{S}_{1j}$
relative to the connected UAV, which are random variables. Therefore, in order
to obtain the average outage probability, it is necessary to have the
distribution function of $\theta_{d_{1j}}$, which is analyzed in the following
two propositions.

Proposition 1. The PDF of RV $\theta_{11}$ is obtained as:

$\displaystyle
f_{\theta_{{11}}}(\theta_{11})=\frac{h_{1}}{R\cos^{2}(\theta_{{11}})}~{}~{}\text{for}~{}~{}0<\theta_{11}<\tan^{-1}\left(\frac{R}{h_{1}}\right).$
(33)

###### Proof:

Please refer to Appendix C. ∎

The distribution of $\theta_{x_{11}}$ which determines the position of
$\mathcal{S}_{11}$ is very important for performance analysis because the
intra-cell interference depends on the position $\mathcal{S}_{11}$. In
Proposition 1, the distribution function of $\theta_{x_{11}}$ is provided,
which is a function of the UAV’s flight height $h_{1}$ and the UAV’s coverage
radius $R$. In Fig. 7a, the distribution function of $\theta_{x_{11}}$ is
plotted for different flight heights and $R=400$m. The results of Fig. 7a show
that $f_{\theta_{11}}(\theta_{11})$ is an ascending function of $\theta_{11}$.
It means that with a low probability, $S_{11}$ is located at a small spacial
angle with respect to the $U_{1}$ and as the height of $U1$ decreases,
$f_{\theta_{11}}$ decreases for lower values of $\theta_{11}$, which is not
desirable.

(a)

(b)

Figure 7: (a) The distribution function of $\theta_{x_{11}}$ for different
$U_{1}$’s flight heights and $R=400$m. (b) The distribution function of
$\theta_{d_{1j}}$ for different values of $\theta_{x_{11}}$.

Proposition 2. The PDF and CDF of $\theta_{d_{1j}}$ are given by:

$\displaystyle f_{\theta_{d_{1j}}}(\theta_{d_{1j}})=\int
f_{\theta_{d_{1j}}|\theta_{x_{1j}}}(\theta_{d_{1j}})f_{\theta_{x_{1j}}}(\theta_{x_{1j}})\text{d}\theta_{x_{1j}}.$
(34) $\displaystyle F_{\theta_{d_{1j}}}(\theta_{d_{1j}})=\int
F_{\theta_{d_{1j}}|\theta_{x_{1j}}}(\theta_{d_{1j}})f_{\theta_{x_{1j}}}(\theta_{x_{1j}})\text{d}\theta_{x_{1j}},$
(35)

where

$\displaystyle
f_{\theta_{d_{1j}}|\theta_{x_{1j}}}(\theta_{d_{1j}})=\frac{h_{1}}{\sqrt{R^{2}-h_{1}^{2}\tan^{2}(\theta_{x_{1j}})}}$
$\displaystyle\times\frac{\sin(\theta_{d_{1j}})}{\cos^{3}(\theta_{d_{1j}})\sqrt{\tan^{2}(\theta_{d_{1j}})-\tan^{2}(\theta_{x_{1j}}-\theta_{x_{11}})}}$
(36) $\displaystyle F_{\theta_{d_{1j}}|\theta_{x_{1j}}}(\theta_{d_{1j}})=$
(37)
$\displaystyle\frac{h_{1}}{\sqrt{R^{2}-h_{1}^{2}\tan^{2}(\theta_{x_{1j}})}}\sqrt{\tan^{2}(\theta_{d_{1j}})-\tan^{2}(\theta_{x_{1j}}-\theta_{x_{11}})},$
(38)

###### Proof:

Please refer to Appendix D. ∎

Note that $\theta_{d_{1j}}$ is a function of $\theta_{x_{1j}}$ and
$\theta_{y_{1j}}$. In order to calculate the average outage probability, we
need to know the PDF and CDF of $\theta_{d_{1j}}$, which are derived in
Proposition 2. The important point is that $\theta_{d_{1j}}$ is also a
function of $\theta_{11}$ which is modeled in proposition 1. In Fig. 7b, the
distribution function of $\theta_{d_{1j}}$ is plotted for different values of
$\theta_{11}$. As mentioned in the previous section, the larger
$\theta_{d_{1j}}$ is, the less interference it creates, and as a result, the
system performance improves. The results of Fig. 7b show that the smaller
$\theta_{11}$ is, the probability of having large values of $\theta_{d_{1j}}$
increases, and as a result, we expect that the interference decreases and the
system performance improves. On the other hand, according to the Proposition
1, as the flight height increases, the probability of having a smaller
$\theta_{11}$ increases, which leads to a decrease in the average outage
probability. The obtained results are summarized in the remark below.

Remark 6. Based on the results obtained from propositions 1 and 2, as the
flight height increases, the probability of smaller values of $\theta_{d11}$
decreases, and as a result, the interference decreases and the system
performance improves. Theorem 3. The average outage probability for the worst
case scenario is derived as:

$\displaystyle\mathbb{P}_{\text{out}}=\int\mathbb{P}_{\text{out}|\theta_{x_{1j}}}f_{x_{1j}}(x_{1j})\text{d}\theta_{x_{1j}}$
(39)

where $\mathbb{P}_{\text{out}|\theta_{x_{1j}}}$ is derived in (IV-A). In
(IV-A), we have
$\mathbb{G}_{1}=\max\\{\alpha_{c},|\theta_{x_{1j}}-\theta_{x_{11}}|\\}$, and
$\mathbb{G}_{2}=(\theta_{d\text{max}}-\mathbb{G}_{1})$.

$\displaystyle\mathbb{P}_{\text{out}|\theta_{d_{1j}}}=\mathbb{P}_{\text{out}|\theta_{x_{1j}}}=\sum_{d=0}^{D-1}\sum_{j=1}^{J-1}\binom{J}{j}\left(\frac{h_{1}}{\sqrt{R^{2}-h_{1}^{2}\tan^{2}(\theta_{x_{1j}})}}\right)^{J-1}$
(40)

###### Proof:

Please refer to Appendix E. ∎

The number of random parameters of the considered system is at least $2(J+1)$,
therefore it requires to solve a $2(J+1)$-dimensional integral, which will be
very time-consuming. To overcome this problem, in Theorem 3, the upper bound
of the average outage probability of the considered system is derived, which
only requires solving a one-dimensional integral and reduces the computational
load, significantly.

###  IV-B Numerical and Simulation Results

Next we provide extensive simulation and numerical results in order to assess
the performance of the considered system in term of outage probability. In our
simulations, we follow the same path presented in the previous section for
simulating interference. After obtaining $5\times 10^{7}$ independent values
of intra- and inter-cell interference, we calculate $5\times 10^{7}$
independent values of the SINR and finally calculate the outage probability.
This outage probability is obtained for only a random distribution of SBSs
whose simulation time is in the order of minutes. To calculate the average
outage probability, we simulate the outage probability for 500 random
distributions of SBSs and finally calculate the average outage probability,
which requires several hours of simulation time. However, using the results of
Theorem 3, the upper bound of outage probability will be computed in the order
of one second.

In Fig. 8, the analytical results obtained from Theorem 3 are compared with
the outage probability results obtained from the simulations. As we observe,
the analytical results compute the upper bound of real outage probability with
acceptable accuracy while it reduce the simulation time, significantly. In
Fig. 8, the results are obtained for different values of $\alpha_{c}$. As we
expect, for the larger $\alpha_{c}$, the probability of having
$\theta_{d_{1j}}<\alpha_{c}$ becomes zero and thus, the outage probability of
the considered system decreases. Another important point is related to the
antenna pattern, which is controlled by parameter $N_{u}$. The results of Fig.
8 clearly show that the parameter $N_{u}$ has a significant effect on the
system performance. First, by increasing $N_{u}$, the main received power
increases and also the side-lobes become smaller and the interference
decreases. Therefore, by increasing $N_{u}$, the outage probability decreases.
However, by increasing $N_{u}$, the beam-width decreases and the sensitivity
of the antenna alignment to the UAV’s vibrations increases. Therefore, the
optimal selection $N_{u}$ is of great importance, and the provided analytical
expressions help us to analyze the effect of various parameters on the system
performance in a short time instead of time consuming Monte-Carlo simulations.

In addition to $\alpha_{c}$, there are other important parameters that affect
the performance of the considered system, the most important of which are
$\sigma_{\theta}$ and $\theta_{11}$. To get a better view, in Fig. 9, the
outage probability is provided for two different values of the UAV’s
instability ($\sigma_{\theta}=1^{o}$ and $1.4^{o}$) and also two different
values of $\theta_{11}=20^{o}$ and $60^{o}$. As we expect, the system
performance decreases by increasing $\sigma_{\theta}$. In addition, the
simulation results show that by reducing $\theta_{11}$ from $60^{o}$ to
$20^{o}$, the performance of the considered system improves. The reason for
this is that by reducing $\theta_{11}$, the probability of having large
$\theta_{d_{1j}}$ increases, and thus, the interference decreases and the
system performance improves.

Figure 8:  Outage probability versus $N_{u}$ for different values of
$\alpha_{c}$.  Figure 9:  Outage probability versus $N_{u}$ for two different
valus of $\sigma_{\theta}$ and two different values of $\theta_{11}$.

##  V Conclusion

In this paper, we have studied mmWave downlink transmissions for an NFP-based
network in the presence of real channel parameters. Our aim was to design a
mmWave downlink and reuse frequency using directional antennas. In particular,
we have performed a comprehensive modeling of inter-cell and intra-cell
interference and by studying the effect of parameters, valuable solutions were
provided to reduce the interference. From the obtained results, the upper
bound of the outage probability was derived. Then, by providing numerical and
simulation results, the effect of important channel parameters on the system
performance was investigated and precious results were obtained, including the
effect of integrating SBSs close to each other, how to choose the optimal
antenna pattern to reduce interference and as well as reduce the negative
effect of UAV fluctuations.

## Appendix A

Based on the results of [48, 49], the main-lobe of the antenna pattern can be
well approximated by Gaussian distribution. Therefore, (14) can be rewritten
as follows:

$\displaystyle P_{r_{11}}=P_{t}|h_{L_{11}}|^{2}G_{0}(N_{s})G_{0}(N_{u})$
$\displaystyle\times\exp\left(-\frac{\left(\tan^{-1}\left(\sqrt{\tan^{2}(\Theta_{x_{1}})+\tan^{2}(\Theta_{y_{1}})}\right)\right)^{2}}{w_{0}^{2}}\right),$
(41)

where $w_{0}=1/N_{u}$. Also, after a comprehensive search, we approximate the
antenna pattern along with the side-lobes as a function of the sum of several
Gaussian functions as:

$\displaystyle G_{u}(N_{u})\simeq
G_{0}(N_{u})\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}\exp\left(-\frac{(\theta-A^{\prime}_{m})^{2}}{w_{m}^{2}}\right)$
(42)

where, $M=6$, and

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&w_{m}=\frac{1}{N_{s}}~{}~{}\&~{}A^{\prime}_{m}=0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}m=0\\\
&w_{m}=\frac{1}{2N_{s}}~{}\&~{}A^{\prime}_{m}=\sin^{-1}\left(\frac{2m+1}{N_{u}}\right),~{}~{}~{}~{}m\neq
0\\\ \end{array}\right.$ (45)

and $A_{0}=1,A_{1}=0.05,A_{2}=0.020,A_{3}=0.011,A_{4}=0.006,A_{5}=0.004$,
$A_{6}=0.003$, and

$\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&A^{\prime\prime}_{m}=1,~{}~{}\text{for}~{}~{}\frac{m}{N_{u}}\leq
1\\\ &A^{\prime\prime}_{m}=0,~{}~{}\text{for}~{}~{}\frac{m}{N_{u}}>1\\\
\end{array}\right.$

In Fig. 3b, for different values of $\phi=0:0.33:45$, the accuracy of the
model obtained from (42) is compared with the real antenna pattern model. As
we see, the approximated model provided in (42) models well the main lobe
along with the maximum amount of side-lobes. We use (42) to model the
interference, and since the actual value of the interference is equal to or
less than the interference obtained from (42) (the actual value of the
interference depends on RV $\phi_{u_{i}}$), the interference obtained from
(42) is for the worst scenario. Therefore, the outage probability obtained
from (42) will be the upper bound of the actual outage probability. Now, from
(10), (III-A1), (A), (42) and using the results of Remarks 3 and 4, the
instantaneous SINR can be written as

$\displaystyle\gamma_{11}\simeq\frac{\exp\left(-\frac{\left(\tan^{-1}\left(\sqrt{\tan^{2}(\Theta_{x_{1}})+\tan^{2}(\Theta_{y_{1}})}\right)\right)^{2}}{w_{0}^{2}}\right)}{\sum_{j=1}^{J}\sum_{m=0}^{M}A_{m}A^{\prime\prime}_{m}\exp\left(-\frac{(\theta_{d_{1j}}-A^{\prime}_{m})^{2}}{w_{m}^{2}}\right)+N^{\prime}_{0}}$
(46)

where
$N^{\prime}_{0}=\frac{\sigma^{2}_{N}}{P_{t}|h_{L_{11}}|^{2}G_{0}(N_{s})G_{0}(N_{u})}$
is the normalized thermal noise. Let us define RV
$\Theta^{\prime}_{1}=\frac{\left(\tan^{-1}\left(\sqrt{\tan^{2}(\Theta_{x_{1}})+\tan^{2}(\Theta_{y_{1}})}\right)\right)^{2}}{w_{0}^{2}}$.
Note that $\Theta^{\prime}_{1}$ is a random variable caused by the instability
and vibrations of $U_{1}$, which is in the order of one degree (depending on
the weather conditions, wind speed and the type of UAV, it can be a little
more or less) [27]. Therefore, since for small values of $x$ we have
$\tan(x)\simeq x$, then $\Theta_{1}$ can be approximated as follows:

$\displaystyle\Theta^{\prime}_{1}\simeq\frac{\Theta_{x_{1}}^{2}+\Theta_{y_{1}}^{2}}{w_{0}^{2}}.$
(47)

Using [50, eqs. (3) and (5)], the distribution of $\Theta$ is obtained as

$\displaystyle f_{\Theta^{\prime}_{1}}(\Theta_{1})=$ (48)
$\displaystyle\frac{1}{\sqrt{\beta_{ux}\beta_{uy}}}\exp\left(-\frac{\beta_{ux}+\beta_{uy}}{2\beta_{ux}\beta_{uy}}\Theta^{\prime}_{1}\right)I_{0}\left(\frac{\beta_{ux}-\beta_{uy}}{2\beta_{ux}\beta_{uy}}\Theta^{\prime}_{1}\right)$

where $\beta_{uw}=\frac{\sqrt{2}\sigma_{\theta_{w}}^{2}}{w_{0}}$, and
$I_{0}(\cdot)$ is the modified Bessel function of the first kind with order
zero. Using (46) and (48), and after some derivations, the PDF of
$\gamma_{11}$ is derived in (24).

## Appendix B

To calculate the outage probability, we need the CDF of $\Theta^{\prime}_{1}$.
Based on (24) and using [51], the CDF of $\Theta^{\prime}_{1}$ is derived as

$\displaystyle F_{\Theta^{\prime}_{1}}(\Theta^{\prime}_{1})$
$\displaystyle=Q\left(\mathbb{B}_{1}\sqrt{\frac{2\Theta^{\prime}_{1}}{\beta_{ux}+\beta_{uy}}},\mathbb{B}_{2}\sqrt{\frac{2\Theta^{\prime}_{1}}{\beta_{ux}+\beta_{uy}}}\right)$
$\displaystyle~{}~{}~{}-Q\left(\mathbb{B}_{2}\sqrt{\frac{2\Theta^{\prime}_{1}}{\beta_{ux}+\beta_{uy}}},\mathbb{B}_{1}\sqrt{\frac{2\Theta^{\prime}_{1}}{\beta_{ux}+\beta_{uy}}}\right)$
(49)

where $Q(a,b)$ is the Marcum Q-function, and

$\displaystyle\left\\{\\!\\!\\!\\!\\!\\!\\!\begin{array}[]{rl}&\mathbb{B}_{1}=\frac{\sqrt{1-T_{q}^{4}}}{2T_{q}}\sqrt{\frac{1+T_{q}}{1-T_{q}}},\\\
&\mathbb{B}_{2}=\mathbb{B}_{1}\frac{1-T_{q}}{1+T_{q}},~{}~{}~{}T_{q}=\frac{\sigma_{\theta_{\text{min}}}}{\sigma_{\theta_{\text{max}}}},\\\
&\sigma_{\theta_{\text{max}}}=\text{max}\\{\sigma_{\theta_{x}},\sigma_{\theta_{y}}\\},~{}~{}~{}~{}\sigma_{\theta_{\text{min}}}=\text{min}\\{\sigma_{\theta_{x}},\sigma_{\theta_{y}}\\}\end{array}\right.$
(53)

Using (46) and (B), and after some derivations, the CDF of $\gamma_{11}$ is
derived in (IV-A).

## Appendix C

Here, we want to derive the distribution of
$\theta_{11}=[\theta_{x_{11}},\theta_{y_{11}}]$. Using (II-C), it can be
formulated as

$\displaystyle\theta_{11}$
$\displaystyle=\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{x_{11}})+\tan^{2}(\theta_{y_{11}})}\right).$
(54)

Considering that each UAV can cover SBSs up to a maximum distance of
$L_{\text{max}}$, so for $U_{1}$ placed at $[0,0,h_{1}]$ in a Cartesian
coordinate system $[x,y,z]\in\mathbb{R}^{1\times 3}$, the SBSs located in a
circular plane with a radius of $R$ can be connected to $U_{1}$ where
$R=\sqrt{L_{\text{max}}^{2}-h_{1}^{2}}$. As mentioned earlier, the position of
$\mathcal{S}_{1j}$ is characterized by $[x_{1j},y_{1j},0]$ in the considered
Cartesian coordinate system where $x_{1j}$ and $y_{1j}$ are uniform RV that

$\displaystyle
f_{x_{1j}}(x_{1j})=\frac{1}{2R},~{}~{}~{}\text{for}~{}~{}~{}-R<x_{1j}<R,$ (55)

and

$\displaystyle
f_{y_{1j}|x_{1j}}(y_{1j})=\frac{1}{2Y_{\text{max}}}~{}~{}~{}\text{for}~{}~{}~{}-Y_{\text{max}}<y_{1j}<Y_{\text{max}}.$
(56)

In (56), we have $Y_{\text{max}}=\sqrt{R^{2}-x_{1j}^{2}}$. Based on (9), we
have

$\displaystyle\theta_{x_{1j}}=\tan^{-1}\left(\frac{x_{1j}}{h_{1}}\right),~{}~{}~{}~{}\&~{}~{}~{}~{}\theta_{y_{1j}}=\tan^{-1}\left(\frac{y_{1j}}{h_{1}}\right).$
(57)

Using (55) and (57), the PDF of $\theta_{x_{1j}}$ is obtained as

$\displaystyle F_{\theta_{x_{11}}}(\theta_{x_{11}})$
$\displaystyle=\text{Prob}\left\\{\tan^{-1}\left(\frac{x_{11}}{h_{1}}\right)<\theta_{x_{11}}\right\\}$
$\displaystyle=F_{x_{11}}\left(h_{1}\tan\left(\theta_{x_{11}}\right)\right).$
(58)

Without loss of generality, we assume that $\mathcal{S}_{11}$ is placed on the
$x$-axis. Therefore, based on (54), we have $\theta_{11}=|\theta_{x_{11}}|$
and thus:

$\displaystyle
f_{\theta_{{11}}}(\theta_{11}=|\theta_{x_{11}}|)=\frac{\text{d}F_{|\theta_{x_{11}}|}(|\theta_{x_{11}}|)}{\text{d}|\theta_{x_{11}}|},~{}~{}~{}0\\!<\\!\theta_{11}\\!<\\!\tan^{-1}\left(\\!\frac{R}{h_{1}}\\!\right)$
(59)

Finally, using (C) and (59), the PDF of $\theta_{11}$ is derived in (33).

## Appendix D

In this appendix, we derive the distribution of $\theta_{d_{1j}}$. Based on
the assumption of Appendix C, $\theta_{d_{1j}}$ can be formulated as

$\displaystyle\theta_{d_{1j}}=\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{x_{1j}}-\theta_{x_{11}})+\tan^{2}(\theta_{y_{1j}})}\right).$
(60)

Based on (60), the CDF of $\theta_{d_{1j}}$ conditioned on $\theta_{x_{1j}}$
is obtained as

$\displaystyle F_{\theta_{d_{1j}}|\theta_{x_{1j}}}$ (61)
$\displaystyle\\!\\!=F_{\theta_{y_{1j}}|\theta_{x_{1j}}}\left(\tan^{-1}\left(\sqrt{\tan^{2}(\theta_{d_{1j}})-\tan^{2}(\theta_{x_{1j}}-\theta_{x_{11}})}\right)\right).$

The CDF of $\theta_{y_{1j}}$ is obtained as

$\displaystyle F_{\theta_{y_{1j}}}(\theta_{y_{1j}})=\int
F_{\theta_{y_{1j}}|\theta_{x_{1j}}}(\theta_{y_{1j}})f_{\theta_{x_{1j}}}(\theta_{x_{1j}})\text{d}\theta_{x_{1j}}.$
(62)

Following the results of Appendix C, the PDF of $\theta_{x_{1j}}$ and the CDF
of $\theta_{y_{1j}}$ conditioned on $\theta_{x_{1j}}$ are derived respectively
as

(63) $\displaystyle
F_{\theta_{y_{1j}}|\theta_{x_{1j}}}(\theta_{y_{1j}})=\frac{h_{1}}{\sqrt{R^{2}-h_{1}^{2}\tan^{2}(\theta_{x_{1j}})}}\tan(\theta_{y_{1j}}).$
(64)

From (61) and (64), we obtain

$\displaystyle F_{\theta_{d_{1j}}|\theta_{x_{1j}}}(\theta_{d_{1j}})=$ (65)
$\displaystyle\frac{h_{1}}{\sqrt{R^{2}-h_{1}^{2}\tan^{2}(\theta_{x_{1j}})}}\sqrt{\tan^{2}(\theta_{d_{1j}})-\tan^{2}(\theta_{x_{1j}}-\theta_{x_{11}})}.$

By derivation from (65), the PDF of $\theta_{d_{1j}}$ conditioned on
$\theta_{x_{1j}}$ is derived in (IV-A).

## Appendix E

We have previously obtained outage probability condition on $\theta_{d_{1j}}$
in Theorem 2. Using (IV-A), the average outage probability of the considered
system is obtained as:

$\displaystyle\mathbb{P}_{\text{out}}=\int...\int\mathbb{P}_{\text{out}|\theta_{d_{1j}}}f_{\theta_{d_{1j}}}(\theta_{d_{1j}})\text{d}\theta_{d_{12}}...\text{d}\theta_{d_{1J}}.$
(66)

Outage probability of the considered system is computed by substituting (IV-A)
and (34) in (66), which requires to solve a $2(J-1)$ dimensional numerical
integration. Since (IV-A) and (34) are approximately complex, solving the
above $2(J-1)$ dimensional integral will be very complicated if not
impossible. In the following, we use the results of Theorem 2 as well as
Propositions 1 and 2, and by removing ineffective terms, we simplify (66).

Based on the results of Theorem 2, by increasing $\theta_{d_{1j}}$, the
average outage probability decreases. In other words, to calculate the average
outage probability, we can only consider the effect of $\mathcal{S}_{1j}$
close to $\mathcal{S}_{11}$ that have the smallest $\theta_{d_{1j}}$. For
better understanding, a numerical example is provided. For example, suppose
that for 4 random $\mathcal{S}_{1j}$, the random value of spatial angle
$\theta_{d_{1j}}$ are equal to $20^{o}$, $40^{o}$, $100^{o}$, and $140^{o}$,
which respectively causes average outage probability of $10^{-3}$, $10^{-5}$,
$10^{-7}$, and $10^{-8}$, respectively. It can be easily seen that the final
outage probability is dominated to the closest $\mathcal{S}_{1j}$, because we
have $(10^{-3}+10^{-5}+10^{-7}+10^{-8})\simeq 10^{-3}$. To find the SBSs with
the smallest $\theta_{d_{1j}}$, we sectorize the parameter $\theta_{d_{1j}}$
as:

$\displaystyle\theta_{d_{1j}}\simeq\frac{d\mathbb{G}_{2}}{D}+\mathbb{G}_{1},~{}~{}~{}\text{for}~{}~{}~{}\frac{d\mathbb{G}_{2}}{D}<\theta_{d_{1j}}-\mathbb{G}_{1}<\frac{(d+1)\mathbb{G}_{2}}{D},$

where $d\in\\{0,1,...,D\\}$,
$\mathbb{G}_{1}=\max\\{\alpha_{c},|\theta_{x_{1j}}-\theta_{x_{11}}|\\}$, and
$\mathbb{G}_{2}=(\theta_{d\text{max}}-\mathbb{G}_{1})$. For obtaining outage
probability, we consider the effect of the $\mathcal{S}_{1j}$ located in the
nearest angular sectors. It is clear that as $D$ increases, the accuracy of
the derivations increases. In the simulations, we show that $D=80$ provides
acceptable accuracy. According to the mentioned points, the
$2(J-1)$-dimensional integration (66) is simplified to the following one-
dimensional integral:

$\displaystyle\mathbb{P}_{\text{out}}=\int\mathbb{P}_{\text{out}|\theta_{x_{1j}}}f_{\theta_{x_{1j}}}(\theta_{x_{1j}})\text{d}\theta_{x_{1j}}.$
(67)

where

$\displaystyle\mathbb{P}_{\text{out}|\theta_{x_{1j}}}=\sum_{d=0}^{D-1}\sum_{j=1}^{J-1}\binom{J}{j}\mathbb{P}_{\text{out}|\frac{d\mathbb{G}_{2}}{D}}$
(68)
$\displaystyle\times\left(\text{Prob}\left\\{\frac{d}{D}<\frac{\theta_{d_{1j}}-\mathbb{G}_{1}}{\mathbb{G}_{2}}<\frac{d+1}{D}\Big{|}\theta_{x_{1j}}\right\\}\right)^{j}$
(69)
$\displaystyle\times\left(\text{Prob}\left\\{\frac{d+1}{D}<\frac{\theta_{d_{1j}}-\mathbb{G}_{1}}{\mathbb{G}_{2}}\Big{|}\theta_{x_{1j}}\right\\}\right)^{J-j-1}.$
(70)

In (68), $\binom{n}{m}$ is the combinations of $m$ and $n$. Also, for each
$d$, the term
$\left(\text{Prob}\left\\{\frac{d}{D}<\frac{\theta_{d_{1j}}-\mathbb{G}_{1}}{\mathbb{G}_{2}}<\frac{d+1}{D}\Big{|}\theta_{x_{1j}}\right\\}\right)^{j}$
is proportional to the probability of placing $j$ of $J-1$ SBSs in the $d$th
sector, and the term
$\left(\text{Prob}\left\\{\frac{d+1}{D}<\frac{\theta_{d_{1j}}-\mathbb{G}_{1}}{\mathbb{G}_{2}}\Big{|}\theta_{x_{1j}}\right\\}\right)^{J-j-1}$
is proportional to the probability of placing the rest of the SBSs in sectors
with
$\theta_{d_{1j}}>\left(\frac{(d+1)\mathbb{G}_{2}}{D}+\mathbb{G}_{1}\right)$.
Finally, using (37), (38), (IV-A), and (68), and after a series of
derivations, the average outage probability of the considered system is
derived in (39) and (IV-A).

## References

  * [1] W. Saad, M. Bennis, and M. Chen, “A vision of 6G wireless systems: Applications, trends, technologies, and open research problems,” _IEEE network_ , vol. 34, no. 3, pp. 134–142, 2019.
  * [2] M. Giordani, M. Polese, M. Mezzavilla, S. Rangan, and M. Zorzi, “Toward 6G networks: Use cases and technologies,” _IEEE Commun. Mag._ , vol. 58, no. 3, pp. 55–61, 2020.
  * [3] Y. Zeng, Q. Wu, and R. Zhang, “Accessing from the sky: A tutorial on UAV communications for 5G and beyond,” _Proceedings of the IEEE_ , vol. 107, no. 12, pp. 2327–2375, 2019.
  * [4] 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.
  * [5] L. U. Khan, W. Saad, D. Niyato, Z. Han, and C. S. Hong, “Digital-twin-enabled 6G: Vision, architectural trends, and future directions,” _IEEE Communications Magazine_ , vol. 60, no. 1, pp. 74–80, 2022.
  * [6] S. Ericsson, Stockholm, “Microwave towards 2020,” _White Paper_ , Sep. 2015, accessed on May 18, 2017.
  * [7] R. Taori and A. Sridharan, “Point-to-multipoint in-band mmwave backhaul for 5G networks,” _IEEE Commun. Mag._ , vol. 53, no. 1, pp. 195–201, 2015.
  * [8] M. T. Dabiri, S. Khankalantary, M. J. Piran, I. S. Ansari, M. Uysal, W. Saad, and C. S. Hong, “UAV-assisted free space optical communication system with amplify-and-forward relaying,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 9, pp. 8926–8936, 2021.
  * [9] M. Alzenad, M. Z. Shakir, H. Yanikomeroglu, and M.-S. Alouini, “FSO-based vertical backhaul/fronthaul framework for 5G+ wireless networks,” _IEEE Commun. Mag._ , vol. 56, no. 1, pp. 218–224, 2018.
  * [10] T. S. Rappaport, G. R. MacCartney, M. K. Samimi, and S. Sun, “Wideband millimeter-wave propagation measurements and channel models for future wireless communication system design,” _IEEE Tran. Commun._ , vol. 63, no. 9, pp. 3029–3056, 2015.
  * [11] S. Ju, Y. Xing, O. Kanhere, and T. S. Rappaport, “Millimeter wave and sub-terahertz spatial statistical channel model for an indoor office building,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 6, pp. 1561–1575, 2021.
  * [12] I. A. Hemadeh, K. Satyanarayana, M. El-Hajjar, and L. Hanzo, “Millimeter-wave communications: Physical channel models, design considerations, antenna constructions, and link-budget,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 2, pp. 870–913, 2017.
  * [13] L. Xu, M. Chen, M. Chen, Z. Yang, C. Chaccour, W. Saad, and C. S. Hong, “Joint location, bandwidth and power optimization for THz-enabled UAV communications,” _IEEE Commun. Let._ , vol. 25, no. 6, pp. 1984–1988, 2021\.
  * [14] M. Gapeyenko, D. Moltchanov, S. Andreev, and R. W. Heath, “Line-of-sight probability for mmwave-based UAV communications in 3D urban grid deployments,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 10, pp. 6566–6579, 2021.
  * [15] W. Khawaja, O. Ozdemir, and I. Guvenc, “UAV air-to-ground channel characterization for mmWave systems,” in _Proc. IEEE 86th Vehicular Technology Conference (VTC-Fall), Toronto, ON, Canada_ , Sep. 2017, pp. 1–5.
  * [16] ——, “Temporal and spatial characteristics of mmWave propagation channels for UAVs,” in _Proc. 2018 11th Global Symposium on Millimeter Waves (GSMM), Boulder, CO, USA, USA_ , May 2018, pp. 1–6.
  * [17] W. Khawaja, O. Ozdemir, and I. Guvenc, “Channel prediction for mmWave ground-to-air propagation under blockage,” _IEEE Antennas Wireless Propag. Lett._ , vol. 20, no. 8, pp. 1364–1368, 2021.
  * [18] N. Tafintsev, D. Moltchanov, M. Gerasimenko, M. Gapeyenko, J. Zhu, S.-p. Yeh, N. Himayat, S. Andreev, Y. Koucheryavy, and M. Valkama, “Aerial access and backhaul in mmwave B5G systems: Performance dynamics and optimization,” _IEEE Commun. Mag._ , vol. 58, no. 2, pp. 93–99, 2020.
  * [19] N. Rupasinghe, Y. Yapıcı, I. Güvenç, and Y. Kakishima, “Non-orthogonal multiple access for mmWave drone networks with limited feedback,” _IEEE Trans. Commun._ , vol. 67, no. 1, pp. 762–777, Jan. 2019\.
  * [20] N. Rupasinghe, Y. Yapıcı, I. Güvenç, M. Ghosh, and Y. Kakishima, “Angle feedback for NOMA transmission in mmWave drone networks,” _IEEE J. Sel. Topics Signal Process._ , vol. 13, no. 3, pp. 628–643, 2019.
  * [21] Y. Yapıcı, N. Rupasinghe, I. Güvenç, H. Dai, and A. Bhuyan, “Physical layer security for noma transmission in mmwave drone networks,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 4, pp. 3568–3582, 2021.
  * [22] M. Gapeyenko, V. Petrov, D. Moltchanov, S. Andreev, N. Himayat, and Y. Koucheryavy, “Flexible and reliable UAV-assisted backhaul operation in 5G mmWave cellular networks,” _IEEE J. Sel, Areas Commun._ , vol. 36, no. 11, pp. 2486–2496, 2018.
  * [23] S. Kumar, S. Suman, and S. De, “Dynamic resource allocation in UAV-enabled mmWave communication networks,” _IEEE Internet Things J._ , vol. 8, no. 12, pp. 9920–9933, 2020.
  * [24] M. T. Dabiri, M. O. Hasna, N. Zorba, and T. Khattab, “Enabling long mmWave aerial backhaul links via fixed-wing UAVs: Performance and design,” _arXiv preprint arXiv:2203.08282_ , 2022.
  * [25] M. T. Dabiri, M. O. Hasna, T. Khattab, and K. Qaraqe, “A study of multihop mmW aerial backhaul links,” _arXiv preprint arXiv:2203.04387_ , 2022.
  * [26] M. T. Dabiri, H. Safi, S. Parsaeefard, and W. Saad, “Analytical channel models for millimeter wave UAV networks under hovering fluctuations,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 4, pp. 2868–2883, 2020.
  * [27] M. T. Dabiri, M. Rezaee, V. Yazdanian, B. Maham, W. Saad, and C. S. Hong, “3D channel characterization and performance analysis of UAV-assisted millimeter wave links,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 110–125, 2020.
  * [28] M. T. Dabiri and M. O. Hasna, “3D uplink channel modeling of UAV-based mmWave fronthaul links for future small cell networks,” _arXiv preprint arXiv:2203.14293_ , 2022, Accepted for IEEE Transactions on Vehicular Technology.
  * [29] A. A. Khuwaja, Y. Chen, N. Zhao, M.-S. Alouini, and P. Dobbins, “A survey of channel modeling for UAV communications,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 4, pp. 2804–2821, 2018.
  * [30] T. M. Nguyen, W. Ajib, and C. Assi, “A novel cooperative NOMA for designing UAV-assisted wireless backhaul networks,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 11, pp. 2497–2507, 2018.
  * [31] D. W. Matolak and R. Sun, “Air–ground channel characterization for unmanned aircraft systems—part I: Methods, measurements, and models for over-water settings,” _IEEE Trans. Vehic. Technol._ , vol. 66, no. 1, pp. 26–44, Jan. 2017.
  * [32] R. Sun and D. W. Matolak, “Air-ground channel characterization for unmanned aircraft systems part II: Hilly and mountainous settings.” _IEEE Trans. Vehic. Technol._ , vol. 66, no. 3, pp. 1913–1925, Mar. 2017.
  * [33] X. Lu, W. Yang, X. Guan, Q. Wu, and Y. Cai, “Robust and secure beamforming for intelligent reflecting surface aided mmwave MISO systems,” _IEEE Wireless Commun. Let._ , vol. 9, no. 12, pp. 2068–2072, 2020.
  * [34] Z. Guan and T. Kulkarni, “On the effects of mobility uncertainties on wireless communications between flying drones in the mmWave/THz bands,” in _Proc. IEEE Conference on Computer Communications Workshops (INFOCOM WKSHPS), Paris, France, France_ , April 2019, pp. 768–773.
  * [35] W. Zhong, L. Xu, X. Liu, Q. Zhu, and J. Zhou, “Adaptive beam design for UAV network with uniform plane array,” _Physical Communication_ , vol. 34, pp. 58–65, 2019.
  * [36] J. Pokorny, A. Ometov, P. Pascual, C. Baquero, P. Masek, A. Pyattaev, A. Garcia, C. Castillo, S. Andreev, J. Hosek _et al._ , “Concept design and performance evaluation of UAV-based backhaul link with antenna steering,” _J. Commun. Netw._ , vol. 20, no. 5, pp. 473–483, 2018.
  * [37] A. S. W. Ghattas, A. A. R. Saad, and E. E. M. Khaled, “Compact patch antenna array for 60 GHz millimeter-wave broadband applications,” _Wireless Personal Communications_ , vol. 114, no. 4, pp. 2821–2839, 2020.
  * [38] M. Asaadi, I. Afifi, and A.-R. Sebak, “High gain and wideband high dense dielectric patch antenna using FSS superstrate for millimeter-wave applications,” _IEEE Access_ , vol. 6, pp. 38 243–38 250, 2018.
  * [39] 3GPP TR 37.840 v12.1.0, “Technical specification group radio access network; study of radio frequency (RF) and electromagnetic compatibility (EMC) requirements for active antenna array system (AAS) base station,” _Tech. Rep._ , 2013.
  * [40] C. A. Balanis, _Antenna theory: analysis and design_. John wiley & sons, 2016.
  * [41] R. W. Heath, N. Gonzalez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” _IEEE journal of selected topics in signal processing_ , vol. 10, no. 3, pp. 436–453, 2016.
  * [42] M. T. Dabiri, S. M. S. Sadough, and M. A. Khalighi, “Channel modeling and parameter optimization for hovering UAV-based free-space optical links,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 2104–2113, 2018.
  * [43] A. Kaadan, H. H. Refai, and P. G. LoPresti, “Multielement FSO transceivers alignment for inter-UAV communications,” _J. Lightw. Technol._ , vol. 32, no. 24, pp. 4785–4795, 2014.
  * [44] M. T. Dabiri, S. M. S. Sadough, and I. S. Ansari, “Tractable optical channel modeling between UAVs,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 12, pp. 11 543–11 550, 2019.
  * [45] T. Bai and R. W. Heath, “Coverage and rate analysis for millimeter-wave cellular networks,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 2, pp. 1100–1114, 2014.
  * [46] A. Al-Hourani, S. Kandeepan, and A. Jamalipour, “Modeling air-to-ground path loss for low altitude platforms in urban environments,” in _2014 IEEE global communications conference_. IEEE, 2014, pp. 2898–2904.
  * [47] A. Al-Hourani, S. Kandeepan, and S. Lardner, “Optimal lap altitude for maximum coverage,” _IEEE Wireless Commun. Let._ , vol. 3, no. 6, pp. 569–572, 2014\.
  * [48] S. Tripathi, N. V. Sabu, A. K. Gupta, and H. S. Dhillon, “Millimeter-wave and terahertz spectrum for 6G wireless,” in _6G Mobile Wireless Networks_. Springer, 2021, pp. 83–121.
  * [49] M. T. Dabiri and M. Hasna, “Pointing error modeling of mmwave to thz high-directional antennas,” _arXiv preprint arXiv:2206.10756_ , 2022.
  * [50] M. Nakagami, “The m-distribution—A general formula of intensity distribution of rapid fading,” in _Statistical methods in radio wave propagation_. Elsevier, 1960, pp. 3–36.
  * [51] (2016) Hoytdistribution. Wolfram Research. [Online]. Available: https://reference.wolfram.com/language/ref/HoytDistribution.html