\section{Introduction}
Unmanned aerial vehicles (UAVs) are receiving significant interest to enhance the wireless coverage of cellular networks~\cite{Mozaffari}. UAVs are used when terrestrial cellular systems get damaged, in hard to reach areas, and wherein there is an occasional need for supplementary coverage. One of the main features of UAVs is the existence of line-of-sight (LOS) links, which can improve the signal quality compared to non-line-of-sight (NLOS) links. The LOS transmission is considered as one of the key factors that can boost capacity in beyond $5$G (B$5$G) ultra dense networks. The placement of UAVs at elevated altitudes can effectively avoid obstacles, thus enabling millimeter wave (mmWave) based technologies such as multiple-input-multiple-output (MIMO) transmissions~\cite{Evans}. In addition, UAVs can be further integrated with emerging network architectures such as Internet of things~\cite{Bushnaq,Bushnaq2} and multiple edge computing~\cite{Amer}. In this work, we consider a hybrid aerial-terrestrial cellular network where UAVs are deployed above a built-up dense urban area to assist terrestrial networks and are backhauled through mmWave links.
\vspace{-0.5cm}
\subsection{Related Work}
\vspace{-0.1cm}
The foreseen potential of UAVs’ to support massive wireless connectivity has attracted significant research efforts to model, analyze, and design UAVs networks. Instances of such research efforts include developing air-to-air and air-to-ground communications channel models~\cite{Khawaja,Chang}, optimal placement and resource allocation for UAVs~\cite{Yang,Cai}, trajectory planning~\cite{Wu, Hu,Bushnaq}, charging ground sensors using UAVs radio frequency (RF) signals~\cite{Saab,Benbuk} and evaluating the performance of UAV-assisted cellular networks~\cite{Wang,Galkin,Qi}.

To best analyze the performance of UAV networks, stochastic geometry becomes indispensable to account for the mutual interference between the active aerial-terrestrial links~\cite{Sultan}. In~\cite{Chetlur}, the authors use stochastic geometry to provide an expression of the coverage probability for a UAV network under guaranteed LOS conditions. However, the model in~\cite{Chetlur} overlooks the effect of LOS/NLOS components of the aerial channel. A probabilistic LOS/NLOS propagation model for aerial channels is adopted in~\cite{Wang} for UAV-assisted cellular networks. Such LOS/NLOS assumption leads to more accurate results and design insights at the expense of more involved analysis. In~\cite{Osama}, the authors use stochastic geometry to model a network where tethered UAVs are deployed to assist traffic offloading under practical challenges.

Terrestrial BSs can be backhauled via wired or wireless links. In contrast, wireless backhauling is mandatory for UAVs to support their mobility. The UAV backhaul links can be provided by terrestrial BSs using sub-$6$~GHz technologies, mmWave technologies~\cite{Xiao}, or free space optics (FSO)~\cite{Pappi}. For instance, the authors in~\cite{Alzenad} consider point-to-point FSO backhaul links for UAVs. Laser beams are used in~\cite{Lahmeri} to boost UAVs operations with both power and backhauling capabilities. Across the RF spectrum, backhauling at the mmWave range is shown to be superior to sub-$6$~GHz backhauling due to the larger bandwidth and the active beam steering capabilities. In this regard, 3GPP envisions an integrated access and backhaul architecture for small cell BSs in which the same infrastructure and spectral resources are used by the macro BSs on the access links and the small BSs on the backhaul links~\cite{zhang2}. This new architecture introduces new design challenges such as the proper partitioning of resources between the access and the backhaul links~\cite{Saha2,Saha1}. As for UAV networks, the authors in~\cite{Xiao} utilize mmWave backhauling, where issues related to beam tracking, LOS blockage and UAV discovery are investigated.

Despite UAV communications multiple promises, modeling and analyzing the impact of UAV backhauling remains relatively unexplored. The main focus in the literature is mostly on the UAV access links while assuming guaranteed UAV backhaul. The work in~\cite{Boris} is one of few exceptions that explicitly accounts for the UAV backhaul. However, this work assumes that backhauling BSs cannot directly serve UEs, which underutilizes the role of terrestrial BSs. Furthermore, the UAV network in~\cite{Boris} is modeled as an infinite Poisson Point Process (PPP) and the backhauling impact is only incorporated via simulations. The backhaul effect is also considered in~\cite{Chen}, where a single UAV is introduced to maximize the throughput in a dense urban environment. Hence, the work in~\cite{Chen} does not account for the possible coexistence of multiple UAVs and the subsequent mutual interference between them. Up to the authors' best knowledge, the joint consideration of the access and mmWave backhaul links in a large-scale hybrid aerial-terrestrial cellular network has not been investigated in literature, and so this paper addresses the details and intricacies of this problem from a stochastic geometry perspective.
\vspace{-0.7cm}
\subsection{Contribution and Organization}
\vspace{-0.2cm}
In this paper, a hybrid aerial-terrestrial network is considered, where UAVs are used to assist BSs in a dense urban environment in providing coverage to UEs. The UAVs are backhauled through mmWave links, which are subject to beamforming misalignment errors. Using stochastic geometry, we aim to evaluate the performance of the network and the effect of its parameters, while highlighting the impact of mmWave backhauling on UE coverage. To the best of the authors' knowledge, this work is the first to jointly consider the access and mmWave backhaul links in a large-scale hybrid aerial-terrestrial network, where the UE can be either served by a BS or by a UAV. The main contributions of this paper can be summarized as follows:
\begin{enumerate}
\item \textbf{\textit{System model:}} We consider a hybrid aerial-terrestrial network with a finite number of backhaul-enabled UAVs. The UAV network is modeled as a binomial point process (BPP), which fits more realistic use cases than infinite PPPs. The developed mathematical framework explicitly accounts for the LOS/NLOS mmWave links along with possible beamforming misalignment errors.
\item \textbf{\textit{Comparison between backhaul aware/unaware transmission models:}} The proposed model differentiates between backhaul aware and backhaul unaware transmission models. For the unaware scenario, the serving UAV transmits data to UE without considering its backhaul link quality. In the backhaul aware scheme, the serving UAV is aware of the backhaul status and does not communicate with the UE in case of backhaul link outage. The obtained results show that prior knowledge of the backhaul link quality relieves parts of the network interference, and hence, improves the overall coverage probability.
\item \textbf{\textit{Design guidelines and insights:}} The analytical results, validated using Monte-Carlo simulations, highlight the fundamental impact of the backhaul link quality on the UE coverage. Based on the proposed framework, we investigate the impact of various parameters including the intensity of terrestrial BSs, the fraction of backhaul-enabled BSs, the altitude and number of UAVs, as well as the beamforming misalignment errors. The results reveal that the number and altitude of UAVs must be chosen carefully to optimize the UE experience.
\end{enumerate}
\vspace{-0.7cm}
\begin{table}[htp!]
	\centering
	\caption{Notations Summary.}
	\vspace{-0.3cm}
	\begin{tabular}{m{0.23\linewidth}|m{0.73\linewidth}}
		\hline\hline
		\textbf{Notation} & \textbf{Description} \\\hline
		$D_{u}$, $r_{u}$, $v_0$ & Disk in which the UAVs are distributed, Radius of $D_{u}$, Distance from the UE to origin\\\hline
		$N_{u}$, $h_{u}$, $\lambda_{g}$, $h_{g}$ & Number and height of UAVs, Density and height of terrestrial BSs with reference to ground\\\hline
		$\delta_{b}$, $\lambda_{b}$ & Fraction and density of terrestrial backhaul-enabled BSs\\\hline
		$\Phi_{u}$, $\Phi_{g}$, $\Phi_{g_l}$, $\Phi_{g_n}$ &Set of UAVs, terrestrial BSs, LOS BSs and NLOS BSs respectively \\\hline
		$\kappa_{u_l}(\cdot)$, $\kappa_{u_n}(\cdot)$, $\kappa_{b_l}(\cdot)$, $\kappa_{b_n}(\cdot)$ & LOS and NLOS probabilities of the access link and backhaul link, respectively\\\hline
		$P_{u}$, $P_{g}$, $P_{b}$ & Transmit power of UAVs, BSs on the access link, BSs on the backhaul link, respectively\\\hline
		$\eta_{g}$, $\eta_{l}$, $\eta_{n}$ & Path-loss exponent parameter for terrestrial BS, LOS UAV, or NLOS UAV, respectively \\\hline
		$m_{l}$, $m_{n}$ & Nakagami-m fading parameter for aerial LOS link, or aerial NLOS link, respectively \\\hline
		$\Omega_{g,x_{i}}$, $\Omega_{u_{l},y_{j}}$, $\Omega_{u_{n},y_{j}}$, $\Omega_{b_{l},x_{i},y_{j}}$, $\Omega_{b_{n},x_{i},y_{j}}$ & Small scale fading gain between UE and $i$-th BS, $j$-th LOS or NLOS UAV on the access link and between $j$-UAV and $i$-th LOS or NLOS BS on the backhaul link, respectively\\\hline
		$s_{g,x_{i}}$, $s_{u_{l}, y_{j}}$, $s_{u_{n}, y_{j}}$, $s_{b_{l},x_{i},y_{j}}$, $s_{b_{n},x_{i},y_{j}}$  & Horizontal distance between the projections of the UE and the $i$-th BS or the $j$-th LOS or NLOS UAV and between the $i$-th LOS or NLOS BS and the $j$-th UAV, respectively \\\hline
		$z_{g,x_{i}}$, $z_{u_{l},y_{j}}$, $z_{u_{l},y_{j}}$, $z_{b_{l},x_{i}, y_{j}}$, $z_{b_{n},x_{i}, y_{j}}$ & Distance between the UE and the $i$-th BS or the $j$-th LOS or NLOS UAV and between the $i$-th LOS or NLOS BS and the $i$-th UAV, respectively\\\hline
		$P^{r}_{g,x_{i}}$, $P^{r}_{u_{l},y_{j}}$, $P^{r}_{u_{n},y_{j}}$, $P^{r}_{b_{l},x_{i}, y_{j}}$, $P^{r}_{b_{n},x_{i}, y_{j}}$ & Received power from the $i$-th BS or the $j$-th LOS/NLOS UAV at the UE on the access link, or from $i$-th LOS/NLOS BS to the $j$-th UAV on the backhaul link, respectively \\\hline
		$x_{g}$, $y_{u_l}$, $y_{u_n}$, $x_{b_l}$, $x_{b_n}$ & Distances between the UE and its serving BS, LOS or NLOS UAV on the access link and between the UAV and its serving LOS or NLOS BS on the backhaul link, respectively\\\hline
		$C_{l}$, $C_{n}$, $\sigma_{b}^2$  & LOS and NLOS near-field path loss and noise power for backhaul connection\\\hline
		$G_{s}^{(\mathrm{max})}$, $G_{s}^{(\mathrm{min})}$, $\theta_{s}$ & Antenna parameters for BSs ($s=g$) and UAVs ($s=u$) for backhaul connection \\\hline
		$\varepsilon_{g}$, $\varepsilon_{u}$, $\sigma_{g}^{2}$, $\sigma_{u}^{2}$ & Additive beam-steering errors and corresponding variances of BS and UAV antennas \\\hline
		$\hat{I}_{g}$, $\hat{I}_{u_l}$, $\hat{I}_{u_n}$ & Interference from all BSs or UAVs except serving BS, LOS or NLOS UAV, respectively\\ \hline
		$I_{u}$, $I_{g_l}$, $I_{g_n}$ & Interference from all UAVs or BSs if UE associates to a BS, a LOS or NLOS UAV, respectively\\\hline
		$A_{g}$, $A_{u_l}$, $A_{u_n}$, $A_{b_l}$, $A_{b_n}$ & Association probabilities between the UE and a BS, a LOS or a NLOS UAV for backhaul unaware and between the UAV and a LOS or a NLOS BS for backhaul, respectively\\\hline
		$S_{l}(\cdot)$, $S_{n}(\cdot)$, $S(\cdot)$ & Conditional backhaul probabilities given that the UAV is associated with a LOS BS, or a NLOS BS and overall backhaul probability, respectively\\\hline
		$P_{cov,g}$, $P_{cov,u_l}$, $P_{cov,u_n}$ & Conditional coverage probabilities given that the UE is associated with a terrestrial BS, a LOS UAV, or a NLOS UAV for backhaul unaware transmission, respectively\\\hline
		$\tilde{A}_{g}$, $\tilde{A}_{u_l}$, $\tilde{A}_{u_n}$, $\tilde{A}_{f}$ & Association probabilities between the UE and a terrestrial BS, a LOS UAV, or a NLOS UAV and service failure probability for backhaul aware transmission, respectively \\\hline
		$\tilde{P}_{cov,g}$, $\tilde{P}_{cov,u_l}$, $\tilde{P}_{cov,u_n}$ & Conditional coverage probabilities given that the UE is associated with a terrestrial BS, a LOS UAV, or a NLOS UAV for backhaul aware transmission, respectively\\\hline
		$\tilde{N_u}$, $I_{\tilde{u}}$ & Number of UAVs with successful backhaul links for backhaul aware transmission, Interference from UAVs with successful backhaul links\\\hline
		$P_{cov}$, $\tilde{P}_{cov}$ & Overall coverage probabilities for backhaul unaware and aware transmission, respectively \\\hline
		$\tau_a$, $\tau_b$& $\mathrm{SIR}$ threshold on the access link, $\mathrm{SINR}$ threshold on the backhaul link\\\hline
		\hline
	\end{tabular}
	\label{table:notations_access}
	\vspace{-0.5cm}
\end{table}

\subsection{Paper Notations and Organization}
\vspace{-0.1cm}
The subscripts $\{\cdot\}_u$, $\{\cdot\}_g$, and $\{\cdot\}_b$ refer to UAV, ground BS and backhaul, respectively. The subscripts $\{\cdot\}_{l}$, $\{\cdot\}_{n}$ refer to LOS and NLOS, respectively. $\tilde{\{\cdot\}}$ denotes the backhaul aware transmission scheme. $\mathbb{P}\{\cdot\}$ denotes probability, $\mathbb{E}[\cdot]$ denotes expectation, $\mathcal{L}_x(\cdot)$ denotes the Laplace transform of a random variable $x$, and $_{2}F_{1}(a,b;c;x)=\sum_{n=0}^{\infty}\frac{(a)_n(b)_n}{(c)_n}\frac{x^n}{n!}$ denotes the Gauss hyper geometric function. Other notations are listed in Table~\ref{table:notations_access}.

The rest of the paper is organized as follows. Section~\ref{sec:system_model} describes the system model. The analysis of the backhaul unaware and aware transmission scenarios are presented in Sections~\ref{sec:unaware} and \ref{sec:aware}, receptively. Numerical results and design insights are discussed in Section~\ref{sec:results} and validated using Monte-Carlo simulations. Finally, the paper is concluded in Section~\ref{sec:conclusion}.
\vspace{-0.6cm}
\section{System Model}\label{sec:system_model}
\vspace{-0.5cm}
\subsection{Network Model}
\vspace{-0.2cm}
UAVs are deployed to improve the downlink (DL) coverage of a one-tier cellular network, where UEs can be served by either terrestrial BSs or UAVs. The BSs provide wireless access to the UEs, and mmWave backhaul links to UAVs. Fig.~\ref{fig:system_model} presents the proposed network model. Since we consider a DL scheme, the UAVs transmit to UEs and receive from BSs for backhauling purposes. The UAVs are considered as simple relays with no queuing capabilities. On the other hand, the BSs are assumed to have saturated buffers which always have packets for their served UEs, where UEs can be directly served from the BS or through a UAV relay. The analysis of spatio-temporal models with non-saturated BSs buffers and queuing UAVs is beyond the scope of this paper and is left as a future research direction.
\begin{figure}[t]
	\centering
	\begin{minipage}[t]{0.5\linewidth}
		\includegraphics[width=\linewidth,height=0.21\textheight]{Figures/model_x0.pdf}
		\vspace{-0.8cm}
		\caption{Proposed system model}
		\label{fig:system_model}
	\end{minipage}
	\begin{minipage}[t]{0.4\linewidth}
		\includegraphics[width=\linewidth,height=0.12\textheight]{Figures/antenna_new-2.pdf}
		\vspace{-0.8cm}
		\caption{Sectored-pattern antenna.}
		\label{fig:antenna}
	\end{minipage}
	\vspace{-1cm}
\end{figure}

The ground BSs are spatially distributed according to a PPP $\Phi_{g}=\{\textbf{x}_{i}\}$, where $\textbf{x}_{i}$ is the location of the $i$-th BS. We assume that the density of the ground BSs is $\lambda_g$ and that they are all at the same height $h_g$ and transmit with the same power $P_{g}$ to UEs. A fraction $\delta_{b}$ of the terrestrial BSs is equipped with backhauling capabilities and transmit to the UAVs with the same power $P_{b}$ for backhaul. Thus, the density of backhaul-enabled BSs $\lambda_{b}$ is equal to $\delta_{b}\lambda_{g}$. We assume that $N_u$ UAVs hover at an altitude $h_u$ and are distributed uniformly in a finite disk $D_u$ with radius $r_u$ forming a BPP $\Phi_u=\{\textbf{y}_{j}\}$, where $\textbf{y}_{j}$ refers to the location of the $j$-th UAV. Without loss of generality, we assume that the BPP is centered around $\textbf{o'}=(0,0,h_{u})$ and that all the UAVs transmit with the same power $P_{u}$ to UEs. For the access link, the BSs and UAVs are equipped with single antennas, so as to serve single antenna UEs on the sub-$6$~GHz band. On the other hand, dedicated antenna arrays at the BSs and UAVs are used to align the mmWave transmissions for backhaul. Using separate antennas for mmWave backhaul and access at the terrestrial BSs has two main benefits. Firstly, it allows the BSs to steer their antennas to align with those of the UAVs, and hence, to maximize the directionality gain. Secondly, it reduces the interference on both access and backhaul links as they operate on different frequency bands. The analysis is conducted for a UE positioned at an arbitrary distance $\textbf{v}_{0}=(v_{0},0,0)$ from the origin $\textbf{o}=(0,0,0)$ on the ground. All frequency resources are universally reused across the network.
\vspace{-0.7cm}
\subsection{Channel Model}
\vspace{-0.3cm}
\subsubsection{Terrestrial Access Channel (BS-UE)}
Due to the nature of terrestrial communication, the fading channel of the BS-UE access link consists of a large-scale fading modeled using a distance-dependent path-loss with path-loss exponent $\eta_g$, and a small-scale Rayleigh fading with exponential distribution and unit mean. The signal power received at the UE from the $i$-th BS located at $x_{i}$ can thus be expressed as $P_{g,x_{i}}^{r}=P_{g}(s_{g,x_{i}}^2+h_{g}^2)^{-\eta_g/2}\Omega_{g,x_{i}}$, where $\Omega_{g,x_{i}}$ is the small-scale fading, $s_{g,x_{i}}=\sqrt{z_{g,x_{i}}^{2}-h_{g}^{2}}$ is the horizontal distance separating the UE and the projection of the $i$-th BS on the ground and $z_{g,x_{i}}=||\textbf{x}_{i}-\textbf{v}_{0}||$ is the actual distance.
\subsubsection{Aerial Access (UAV-UE) and Backhaul (BS-UAV) Channels}\label{sec:aerial_channel}
Both the aerial access and backhaul channels are affected by obstacles in the environment which break the LOS links. According to the ITU recommendation report~\cite{ITU}, the probability of a LOS link between a transmitter and a receiver, with heights $h_{\mathrm{TX}}$ and $h_{\mathrm{RX}}$, respectively, is given by~\cite{ITU}:
\begin{equation}\small
P_{\mathrm{LOS}}(r)=\prod_{n=0}^{m}\left[1-\exp\left(-\frac{h_{\mathrm{TX}}-\frac{(n+\frac{1}{2})\left(h_{\mathrm{TX}}-h_{\mathrm{RX}}\right)^{2}}{m+1}}{2\gamma^{2}}\right)\right]
\label{eq:halim}
\end{equation}
where $r$ denotes the horizontal distance between the transmitter and the receiver and $m=\lfloor\frac{r\sqrt{\alpha \beta}}{1000}-1\rfloor$, where $\alpha$, $\beta$ and $\gamma$ are environment related parameters given in Table I in~\cite{Holis}. As we can see from (\ref{eq:halim}), the LOS probability is not a continuous function of the horizontal distance which makes the analysis intractable. As the UE is considered at the ground level, we use the approximation proposed in~\cite{Hourani} to simplify the LOS probability on the access link between a UAV located at a distance $r$ and the UE by a modified Sigmoid function given as:
\begin{equation}\small
\kappa_{u_l}(r)=\frac{1}{1+a \exp (-b[\frac{180}{\pi}\arctan(h_u/\sqrt{r^{2}-h_{u}^{2}})-a])},
\label{eq:LOS_probability}
\end{equation}
where $\sqrt{r^{2}-h_{u}^{2}}$ is the Euclidean horizontal distance separating the projections of the UAV and the UE and $a$ and $b$ are constant values that depend on the environment and are given in Table I of~\cite{Bor}. The NLOS probability is $\kappa_{u_n}=1-\kappa_{u_l}$.

The LOS probability model in (\ref{eq:LOS_probability}) used for the UAV-UE channel was initially developed and approximated for UEs positioned at low height levels from the ground. This model is unsuitable for the BS-UAV aerial channel as both the BSs and UAVs are positioned above the ground level. We use the approximation proposed in~\cite{Cherif} to simplify (\ref{eq:halim}) as a more tractable exponential function, representing the LOS probability between a UAV and a BS separated by a horizontal distance $r$ on the backhaul link as:
\begin{equation}\small
\kappa_{b_l}(r)=-c\exp\left(-d \arctan \left(\frac{|h_u-h_g|}{r}\right)\right)+e,
\label{eq:LOS_probability2}
\end{equation}
where $r$ is the horizontal distance separating the projections of the BS and the UAV, $c$, $d$ and $e$ are parameters that depend on the environment and height of terrestrial BSs and are given in Table II of~\cite{Cherif}. The NLOS probability on the backhaul link is given as $\kappa_{b_n}=1-\kappa_{b_l}$.

The aerial fading channel is characterized using the combination of two components: (a) a large-scale fading modeled using a distance-dependent path-loss and (b) a small-scale Nakagami-$m$ fading modeled using a gamma-distributed random variable. We consider different path-loss exponents ($\eta_{l}$ for LOS and $\eta_{n}$ for NLOS) and fading parameters for the LOS and the NLOS links ($m_{l}$ for LOS and $m_n$ for NLOS). The power received at the reference UE from the $j$-th UAV located at $\textbf{y}_{j}$ is $P^{r}_{u_{\zeta},y_{j}}=P_{u}z_{u_{\zeta},y_{j}}^{-\eta_{\zeta}}\Omega_{u_{\zeta},y_{j}}$, where $z_{u_{\zeta},y_{j}}=||\textbf{y}_{i}-\textbf{v}_{0}||$ is the distance separating the UE from the $i$-th UAV, $\Omega_{u_{\zeta},y_{j}}$ is the small-scale fading, $\eta_{\zeta}$ is the path-loss exponent, and $\zeta\in\{l,n\}$ indicates whether the $j$-th UAV has a LOS link or NLOS link with the UE.

UAVs are connected to ground BSs through mmWave backhaul links. Although mmWave channels share basic propagation characteristics like power law path-loss with sub-$6$~GHz channels, they also have critical distinctions. For instance, beamforming is important to compensate for the high propagation loss in mmWave frequencies as well as the consumed power in analog to digital conversion in large antenna arrays. Thus, the BSs and UAVs must steer their antennas to maximize the directionality gain. However, aligning the antennas is subject to beam-steering errors, which affects the transmission on the backhaul link. We approximate the array patterns of the UAVs and the BSs antennas by the model shown in Fig.~\ref{fig:antenna} and given as
	\vspace{-0.1cm}
	\begin{equation}\small
	G_{s}(\varphi)=\begin{cases}
	G_{s}^{(\text{max})}, & |\varphi|\leq\theta_{s} \\
	G_{s}^{(\text{min})}, & |\varphi|>\theta_{s}
	\end{cases}
	\label{eq:gain}
	\vspace{-0.2cm}
	\end{equation}
	where $\varphi\in[-\pi,\pi)$ is the angle of the boresight direction, $G_{s}^{(\text{max})}$, $G_{s}^{(\text{min})}$ and $\theta_{s}$ are the gains of the main and side lobes and the beamwidth for the BSs and UAVs $(s\in\{g,u\})$, respectively.

	For the desired link, the directivity gain is $G_{b,0}=G_{g}^{(\text{max})}G_{u}^{(\text{min})}$ in the absence of the beamforming misalignment errors modeled explicitly in Section~\ref{sec:misalignment}. For the interfering BSs, the beams are assumed to be randomly oriented with respect to each other, and the steering angles are distributed uniformly in $[-\pi,\pi)$. Since (\ref{eq:gain}) produces either gains $G_{s}^{(\text{max})}$ and $G_{s}^{(\text{min})}$, $s\in\{g,u\}$ over all possible input angles, the resulting gain distribution of a BS or UAV is given as
	\vspace{-0.1cm}
	\begin{equation}\small
	f_{G_{s,I}}(g)=\frac{\theta_{s}}{2\pi}\delta\left(g-G_{s}^{(\text{max})}\right)+\left(1-\frac{\theta_{s}}{2\pi}\right)\delta\left(g-G_{s}^{(\text{min})}\right)
	\label{eq:fsI}
	\vspace{-0.1cm}
	\end{equation}
	where $\delta(\cdot)$ is the Dirac delta function, and $s\in\{g,u\}$.

	Accordingly, the directivity gain of an interfering link is a discrete random variable denoted as $G_{b,I}$, with a probability mass function (PMF) that can be formulated as $f_{G_{b,I}}(g)=\left(f_{G_{g,I}}\otimes f_{G_{u,I}}\right)(g)$, where $f_{G_{g,I}}(\cdot)$ and $f_{G_{u,I}}(\cdot)$ are the PMFs of the antenna gains of the interfering BS and the corresponding UAV. Finally, based on (\ref{eq:gain}) and (\ref{eq:fsI}), the directivity gain of the interfering links takes the values $G_{k}$ with probability $p_{k}$ ($k\in\{1,2,3,4\}$) defined in Table~\ref{table:antenna}, where $c_{g}=\frac{\theta_{g}}{2\pi}$ and  $c_{u}=\frac{\theta_{u}}{2\pi}$. $\theta_{g}$ and $\theta_{u}$ are the beamwidth of the BS and the UAV antennas.
\begin{table}[t]
	\centering
	\caption{Probability mass function of $G_{b,I}$ and $G_{b,0}$.}
	\vspace{-0.3cm}
	\begin{tabular}{c|c|c|c|c}
		\hline
		\textbf{k} & \textbf{1} & \textbf{2} & \textbf{3} & \textbf{4} \\\hline
		$G_{k}$ & $G_{g}^{(\mathrm{max})}G_{u}^{(\mathrm{max})}$&$G_{g}^{(\mathrm{max})}G_{u}^{(\mathrm{min})}$ & $G_{g}^{(\mathrm{min})}G_{u}^{(\mathrm{max})}$ & $G_{g}^{(\mathrm{min})}G_{u}^{(\mathrm{min})}$\\\hline
		$p_{k}$ & $c_{g} c_{u}$ & $c_{g}(1-c_{u})$& $(1-c_{g})c_{u}$ & $(1-c_{g})(1-c_{u})$ \\\hline
		$p_{k,0}$ & $F_{\lvert\varepsilon_{g}\rvert}(\frac{\theta_{g}}{2})F_{\lvert\varepsilon_{u}\rvert}(\frac{\theta_{u}}{2})$ & $F_{\lvert\varepsilon_{g}\rvert}(\frac{\theta_{g}}{2})\bar{F}_{\lvert\varepsilon_{u}\rvert}(\frac{\theta_{u}}{2})$ & $\bar{F}_{\lvert\varepsilon_{g}\rvert}(\frac{\theta_{g}}{2})F_{\lvert\varepsilon_{u}\rvert}(\frac{\theta_{u}}{2})$ & $\bar{F}_{\lvert\varepsilon_{g}\rvert}(\frac{\theta_{g}}{2})\bar{F}_{\lvert\varepsilon_{u}\rvert}(\frac{\theta_{u}}{2})$\small\\\hline
	\end{tabular}
	\label{table:antenna}
	\vspace{-0.9cm}
\end{table}

The received power from the $i$-th BS located at $\textbf{x}_{i}$ at the $j$-th UAV located at $\textbf{y}_{i}$ is given as $P^{r}_{b,x_{i},y_{j}}=P_{b}G_{b}C_{\xi}\left(s_{b_{\xi},x_{i},y_{j}}^{2}+\Delta_{h}^{2}\right)^{-\eta_{\xi}/2}\Omega_{\xi,x_{i},y_{j}}$, where $P_{b}$ is the BS transmit power on the backhaul link, $\Delta_{h}=|h_{u}-h_{g}|$ is the difference between the BSs and the UAVs heights, $\Omega_{b_{\xi},x_{i}, y_{j}}$ is the small-scale fading, $C_{\xi}$ is the path-loss intercept and $\eta_{\xi}$ is the path-loss exponent, where $\xi\in\{l,n\}$ indicates whether the $i$-th BS has a LOS link or a NLOS link with the $j$-th UAV.
\vspace{-0.7cm}
\subsection{Association Strategy and Performance Metrics}
\vspace{-0.2cm}
This paper studies the performance of a hybrid aerial-terrestrial network where both BSs and UAVs are used to serve UEs. That is, each UE connects to the BS or the UAV that offers the maximum average received signal strength. Note that the UE can either connect to a LOS or a NLOS UAV. The coverage probability, defined as the probability that the signal-to-interference ratio ($\mathrm{SIR}$) exceeds a threshold $\tau_a$, is the main performance metric. Each UAV selects its backhauling BS based on the minimum path-loss association rule. Hence, the BS that provides the backhaul is the best in terms of signal strength. Each UAV and its backhauling BS steer their beams to align their transmissions in order to strengthen the intended link and mitigate interference. Such beam steering along with the high bandwidth and severe path-loss of mmWave frequencies relieve the interference dominance and make the backhaul links sensitive to the ambient noise. Thus, it is more adequate to assess the backhaul links through the signal-to-interference-and-noise-ratio ($\mathrm{SINR}$) rather than the $\mathrm{SIR}$. A backhaul transmission is considered successful if and only if the $\mathrm{SINR}$ exceeds a threshold $\tau_b$, which defines the backhaul probability.

We differentiate between two schemes, namely, backhaul unaware and backhaul aware transmission schemes. In the first scheme, the serving UAV transmits directly to the UE regardless of whether the backhaul link is successful or not. Thus, no further post processing is required when the UAV receives packets from the BS and forwards them to the UE. The aerial coverage in the unaware scheme requires two conditions:
(i) $\mathrm{SIR}>\tau_a$: The received $\mathrm{SIR}$ at the UE from its serving UAV must exceed a threshold $\tau_a$, (ii) $\mathrm{SINR}>\tau_b$: The received $\mathrm{SINR}$ at the serving UAV from the BS to which it connects for backhaul support must exceed a threshold $\tau_b$.

In the backhaul aware scenario, the serving UAV is aware of the backhaul status and checks the integrity of the packets before transmitting them to the UE. This will impose a processing burden on the UAV side but will improve the probability of coverage and limits the aerial interference to the subset of UAVs with successful backhaul links only. An active transmission in the aware scheme implies that the condition ($\mathrm{SINR}>\tau_b$) is already satisfied. Otherwise, the UE halts operation and goes to a service failure status in this specific time slot. Thus, the backhaul aware transmission success is only subject to the $\mathrm{SIR}$ condition ($\mathrm{SIR}>\tau_a$) of the access link.
\vspace{-0.6cm}
\section{Backhaul unaware transmission}\label{sec:unaware}
\vspace{-0.15cm}
In this section, we derive the coverage probability for backhaul unaware transmission. In this scenario, the association rule is only dependent on the quality of the UE-UAV/BS connection. However, a UAV can serve the UE only if it has a successful backhaul link with a ground BS.
\vspace{-0.6cm}
\subsection{Association Probabilities}
\vspace{-0.1cm}
To get coverage, the UE can associate to either a terrestrial BS or a UAV. The set of UAVs is divided into two sub-processes according to the LOS/NLOS status of the access links. According to the association rule, each UE connects to the BS or the UAV that offers the maximum average received power. The serving BS/UAV is not necessarily the nearest to the UE due to the difference in path-loss parameters and transmit powers. However, within a particular set, the path-loss parameters and transmit powers are the same for all links. Therefore, for a specific set, the closest BS/UAV provides a larger average received power than that provided by any other in this set. Thus, the UE is served by the closest BS, LOS UAV or NLOS UAV. The association rule implies an exclusion region on the locations of the nearest interfering device and therefore on the locations of all the other interfering devices in each set for each association type. Specifically, when the UE associates with a terrestrial BS located at a horizontal distance $x_{g}$, all the LOS UAVs are further than $E_{g_{l}}(x_{g})$. Similarly, all the NLOS UAVs are further than $E_{g_{n}}(x_{g})$.
\begin{figure}[t]
	\centering
	\begin{minipage}[t]{0.47\linewidth}
		\includegraphics[width=\linewidth,height=0.2\textheight]{Figures/exclusion_area.pdf}
		\vspace{-1.1cm}
		\caption{Exclusion regions $E_{g_{l}}(x_{g})$ and $E_{g_{n}}(x_{g})$ on the LOS and NLOS UAVs.}
		\label{fig:exclusion_region}
	\end{minipage}
	\begin{minipage}[t]{0.47\linewidth}
		\includegraphics[width=\linewidth,height=0.2\textheight]{Figures/exclusion_tareq.pdf}
		\vspace{-1.1cm}
		\caption{Top view snapshot of BSs, LOS/NLOS UAVs and projections of exclusions regions.}
		\label{fig:exclusion_region_horizontal}
	\end{minipage}
	\vspace{-0.8cm}
\end{figure}
Fig.~\ref{fig:exclusion_region} visualizes the exclusion regions on the LOS and NLOS UAVs when the UE associates to a terrestrial BS. Fig.~\ref{fig:exclusion_region_horizontal} presents a top view snapshot on the network showing the exclusion regions on the LOS and NLOS UAVs locations as red dashed circles. Similarly, the exclusion regions $E_{u_l}(x_{u_l})$ and $E_{ln}(x_{u_l})$ are created based on the respective locations of the terrestrial BSs and the NLOS UAVs when the UE associates with a LOS UAV where $x_{u_{l}}$ is the distance to the serving UAV. Finally, the exclusions regions $E_{u_n}(x_{u_n})$ and $E_{nl}(x_{u_n})$ are created based on the respective locations of the BSs and the LOS UAVs, when the UE associates to a NLOS UAV with $x_{u_n}$ being the distance to the serving NLOS UAV. The respective exclusion regions expressions of $E_{g_l}(x)$, $E_{g_n}(x)$, $E_{u_l}(x)$, $E_{u_n}(x)$, $E_{ln}(x)$, and $E_{nl}(x)$ are given as:
\begin{equation}\small
E_{g_\zeta}(x)=\left(\frac{P_{u}}{P_{g}}\right)^{\frac{1}{\eta_{\zeta}}}\left(x^2+h_{g}^2\right)^{\frac{\eta_g}{2\eta_{\zeta}}}, \: E_{u_\zeta}(x)=\sqrt{\left(\frac{P_{g}}{P_{u}}\right)^{\frac{2}{\eta_g}}x^{\frac{2\eta{\zeta}}{\eta_{g}}}-h_{g}^2}, \: E_{\zeta\bar{\zeta}}(x)=x^{\frac{\eta_{\zeta}}{\eta_{\bar{\zeta}}}},
\label{eq:Eg}
\end{equation}
where $\zeta \in \{l,n\}$, $\bar{\zeta}=n$ if $\zeta=l$, and $\bar{\zeta}=l$ if $\zeta=n$.

The association probabilities are defined as the probabilities that the UE associates to either a LOS UAV, a NLOS UAV, or a BS, which are characterized in the following lemma.
\vspace{-0.4cm}
\begin{lemma}
	\emph{Denote by $A_{u_l}$, $A_{u_n}$ and $A_g$ the probabilities that the UE at $\textbf{v}_0$ is served by a LOS UAV, a NLOS UAV or a ground BS, respectively. $A_{u_l}$, $A_{u_n}$ and $A_{g}$ are then given as follows
		\vspace{-0.1cm}
		\begin{equation}\small
		A_{u_\zeta}=N_u\int_{h_u}^{w_{p}}f_{W}(r)\kappa_{u_\zeta}(r)\exp\left(-\pi\lambda_g E^2_{u_\zeta}(r)\right)\left(\int_{r}^{w_{p}}f_{W}(w)\kappa_{u_\zeta}(w)\mathrm{d}w+\int_{E_{\zeta\bar{\zeta}}(r)}^{w_{p}}f_{W}(w)\kappa_{u_{\bar{\zeta}}}(w)\mathrm{d}w\right)^{N_u-1}\mathrm{d}r,
		\label{eq:Aa}
		\end{equation}
		\vspace{-0.8cm}
		\begin{equation}\small
		A_{g}=2\pi\lambda_g\int_{0}^{E_{u}}r \exp\left(-\pi\lambda_g r^2\right)\left(\int_{E_{g_l}(r)}^{w_{p}}f_{W}(w)\kappa_{u_l}(w)\mathrm{d}w+\int_{E_{g_n}(r)}^{w_{p}}f_{W}(w)\kappa_{u_n}(w)\mathrm{d}w\right)^{N_{u}}\mathrm{d}r,
		\label{eq:Ag}
		\end{equation}
		where
		$\zeta \in \{l,n\}$, $\bar{\zeta}=n$ if $\zeta=l$, and $\bar{\zeta}=l$ if $\zeta=n$. $E_u=\mathrm{max}\left(E_{u_l}(w_p), E_{u_n}(w_p)\right)$ and $f_{W}(\cdot)$ is the probability density function (PDF) of the distance $w$ from an arbitrary UAV to the reference UE and is given in Lemma~$2$ of~\cite{Chetlur} as
		\vspace{-0.3cm}
		\begin{equation}\small
		f_{W}(w)=\begin{cases}
		f_{W_{1}}(w)=\frac{2w}{r_{u}^2}, & h_u\leq w\leq w_m\\
		f_{W_{2}}(w)=\frac{2w}{\pi r_{u}^2}\arccos\left(\frac{w^2+v_{0}^2-d^2}{2v_0\sqrt{w^2-h_{u}^2}}\right), & w_m\leq w \leq w_p
		\end{cases}
		\label{eq:fw}
		\end{equation}
		where $r_u$ is the radius of the disk of UAVs, $w_m=\sqrt{(r_u-v_0)^2+h_u^2}$ and $w_p=\sqrt{(r_u+v_0)^2+h_u^2}$.}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:association}.}
	\end{IEEEproof}
	\label{lemma:coverage_association}
	\vspace{-0.2cm}
\end{lemma}
Note that the upper limit of the outer integral in (\ref{eq:Ag}) is not $\infty$ and is equal to $E_{u}=\mathrm{max}\left(E_{u_l}(w_p), E_{u_n}(w_p)\right)$, where $E_{u_l}(\cdot)$ and $E_{u_n}(\cdot)$ are given in (\ref{eq:Eg}). Since $w_p$ is the farthest distance between any UAV and the UE, the signal received from the nearest BS is lower than that of any LOS UAV or NLOS UAV if this BS is located further than $E_{u_l}(w_p)$ or $E_{u_n}(w_p)$.
\vspace{-0.5cm}
\subsection{UE-BS Conditional Coverage Probability}
The conditional coverage probability is the probability that the received $\mathrm{SIR}$ is higher than the threshold $\tau_a$ given the association status. When the UE associates to a BS, this probability is expressed as $P_{cov,g}=\mathbb{P}\left[\mathrm{SIR}\geq \tau_a | s=g \right]$ where $\mathbb{P}\left[s=g\right]=A_g$ is the BS association probability. Due to the universal frequency reuse, the aggregate interference $I_{agg,g}$ includes all non-serving BSs (denoted as $\hat{I}_{g}$) and all UAVs (denoted as $I_{u}$). $P_{cov,g}$  is given in Lemma~\ref{lemma:pcovg}.
\vspace{-0.2cm}
\begin{lemma}
	\emph{The conditional coverage probability $P_{cov,g}$ given that the UE connects to a BS is
		\begin{equation}\small
		P_{cov,g}=\int_{0}^{E_u}\mathcal{L}_{\hat{I}_g}(s_1)\mathcal{L}_{I_u}(s_1)f_{X_g}(x_g)\mathrm{d}x_{g},
		\end{equation}
		where $s_{1}=\frac{\tau_a(x_g^2+h_g^2)^{\frac{\eta_g}{2}}}{P_{g}}$, and $E_u=\mathrm{max}\left(E_{u_l}(w_p), E_{u_n}(w_p)\right)$. $\mathcal{L}_{\hat{I}_g}(s_{1})$ and $\mathcal{L}_{I_u}(s_{1})$ are the Laplace transforms of the aggregate interference of the interfering BSs and all the UAVs. $f_{X_g}(x_g)$ is the PDF of the conditional distance $x_g$ to the serving BS.}
	\begin{IEEEproof}
		\emph{The conditional coverage probability $P_{cov,g}$ is calculated as
			\begin{equation}\small
			\begin{aligned}
			&P_{cov,g}=\mathbb{P}\left[\mathrm{SIR}\geq \tau_a | s=g \right]=\mathbb{P}\left[\frac{P_{g}(x_{g}^2+h_{g}^2)^{-\eta_g/2}\Omega_{g,0}}{I_{agg,g}}\geq \tau_a | s=g \right]\\
			&\stackrel{(a)}{=}\mathbb{E}_{X_g}\left[\mathbb{E}_{I_{agg,g}}\left[\exp\left(-\frac{\tau_a (\hat{I}_g+I_u)}{P_{g}(x_{g}^2+h_{g}^2)^{-\eta_g/2}}\right)\right]\right]\stackrel{(b)}{=}\int_{0}^{E_u}\mathcal{L}_{\hat{I}_g}(s_1)\mathcal{L}_{I_u}(s_1)f_{X_g}(x_g)\mathrm{d}x_{g},
			\end{aligned}
			\label{eq:pcovgproof}
			\end{equation}
			where $s_1=\frac{\tau_a(x_g^2+h_g^2)^{\frac{\eta_g}{2}}}{P_{g}}$, (a) follows from the exponential distribution of $\Omega_{g,0}$ and from the expression of the interference $I_{agg,g}=\hat{I}_{g}+I_u$, and (b) follows from the independence of $\hat{I}_g$ and $I_{u}$ and from the Laplace transform definition. The upper limit of the integral in (\ref{eq:pcovgproof}) $E_{u}=\mathrm{max}\left(E_{u_l}(w_p), E_{u_n}(w_p)\right)$ is the maximum distance from the UE to its serving BS considering the exclusion regions of the LOS and NLOS UAVs. $E_{u_l}(\cdot)$ and $E_{u_n}(\cdot)$ are presented in (\ref{eq:Eg}).}
	\end{IEEEproof}
	\label{lemma:pcovg}
	\vspace{-0.1cm}
\end{lemma}
To obtain the expression of $P_{cov,g}$, Laplace transforms of interference terms and the PDF of the distance to serving BS must be computed. Lemma~\ref{lemma:UE_BS_laplace} and Lemma~\ref{lemma:fxg} present these results.
\vspace{-0.4cm}
\begin{lemma}
	\emph{The Laplace transform of the interference $\hat{I}_g$ of all BSs except the serving BS is
		\begin{equation}\small
		\mathcal{L}_{\hat{I}_{g}}(s_{1})=\exp\left[\frac{-2\pi\lambda_g s_1 P_{g}\left(x_{g}^2+h_{g}^2\right)}{\left(\eta_{g}-2\right)\left(\left(x_{g}^2+h_{g}^2\right)^{\frac{\eta_{g}}{2}}+s_1 P_{g}\right)}{_{2}}F_{1}\left(1,1;2-\frac{2}{\eta_{g}};\frac{1}{1+\frac{\left(x_{g}^2+h_{g}^2\right)^{\frac{\eta_{g}}{2}}}{s_1 P_{g}}}\right)\right].
		\label{eq:Ligb}
		\end{equation}
		where ${_{2}}F_{1}$ denotes the Gauss hyper geometric function.
		The Laplace transform of the interference $I_u$ from the UAVs when the UE associates to a BS is
		\begin{equation}\small
		\begin{aligned}
		\mathcal{L}_{I_{u}}(s_{1})&=\left[\frac{1}{\int_{E_{g_l}(x_{g})}^{w_{p}}f_{W}(w)\kappa_{u_l}(w)\mathrm{d}w+\int_{E_{g_n}(x_{g})}^{w_{p}}f_{W}(w)\kappa_{u_n}(w)\mathrm{d}w}\left(\int_{E_{g_l}(x_{g})}^{w_{p}}\left(1+\frac{s_{1}P_{u}v^{-\eta_{l}}}{m_{l}}\right)^{-m_{l}}\right.\right.\\
		&\times f_{W}(v)\kappa_{u_l}(v)\mathrm{d}v+\int_{E_{g_n}(x_{g})}^{w_{p}}\left(1+\frac{s_{1}P_{u}v^{-\eta_{n}}}{m_{n}}\right)^{-m_{n}}f_{W}(v)\kappa_{u_n}(v)\mathrm{d}v\bigg)\bigg]^{N_u},
		\end{aligned}
		\label{eq:Lia}
		\end{equation}
		where $E_{g_l}(x_g)$ and $E_{g_n}(x_g)$ are the minimum distances at which the LOS and NLOS UAVs are placed when the UE associates to a BS at $x_g$. $f_{W}(\cdot)$ and $\kappa_{u_l}(\cdot)$ are given in (\ref{eq:fw}) and (\ref{eq:LOS_probability}).
	}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:laplace}}
	\end{IEEEproof}
	\label{lemma:UE_BS_laplace}
	\vspace{-0.2cm}
\end{lemma}
The PDF of the conditional distance from the UE to the serving BS $x_g$ is provided in Lemma~\ref{lemma:fxg}.
\vspace{-1cm}
\begin{lemma}
	\emph{The PDF of the distance $x_g$ from the UE to its serving BS is given as
		\begin{equation}\small
		f_{X_{g}}(x_g)=\frac{1}{A_{g}}2\pi\lambda_g x_g \exp\left(-\pi\lambda_g x_{g}^2\right)\left(\int_{E_{g_l}(x_g)}^{w_p}f_{W}(w)\kappa_{u_l}(w)\mathrm{d}w + \int_{E_{g_n}(x_g)}^{w_p}f_{W}(w)\kappa_{u_n}(w)\mathrm{d}w\right)^{N_u},
		\label{eq:fxg}
		\end{equation}
		where $E_{g_l}(\cdot)$ and $E_{g_n}(\cdot)$ are given in (\ref{eq:Eg}). $A_{g}$ and $f_W(\cdot)$ are given in Lemma~\ref{lemma:coverage_association}.
	}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:distance}.}
	\end{IEEEproof}
	\label{lemma:fxg}
\end{lemma}
\vspace{-1cm}
\subsection{UE-UAV Conditional Coverage Probability}
For a UAV to serve a UE, its backhaul link should be successful. The conditional coverage probability given that the UE associates to a LOS/NLOS UAV is the joint probability of two events and can be expressed as $P_{cov,u_{\zeta}}=\mathbb{P}\left[\mathrm{SIR}\geq \tau_a,\mathrm{SINR}\geq\tau_b\right | s=u_{\zeta}]$, where the first term corresponds to the UAV-UE link and the second term to the BS-UAV backhaul. Since the UE and the UAV connect to different BSs, and for mathematical tractability, we approximate $P_{cov,u_{\zeta}}\approx \mathbb{P}\left[\mathrm{SIR}\geq \tau_a | s=u_{\zeta}\right]\times\mathbb{P}\left[\mathrm{SINR}\geq\tau_b | s=u_{\zeta} \right]$. Such approximation is further validated numerically in Section~\ref{sec:results}. Note that, due to the different propagation environments considered for access and backhaul, and the different used probabilistic LOS models, the backhaul probability is independent of the UE association and the access coverage.
\subsubsection{Backhaul Probability}
We start by providing an expression for the backhaul probability $S(\tau_b)=\mathbb{P}\left[\mathrm{SINR}\geq\tau_b\right]$. The analysis is performed for a reference UAV located at height $h_u$ and at the center of the disk $D_u$. The set of backhaul-enabled BSs is divided into two sub-processes: the LOS BSs set $\Phi_{g_l}$ with density $\lambda_b\kappa_{b_l}(r)$ and the NLOS BSs set $\Phi_{g_n}$ with density $\lambda_b \kappa_{b_n}(r)$, where $r$ is the UAV-BS distance and $\kappa_{b_l}(\cdot)$ is given in (\ref{eq:LOS_probability2}). Since the UAV connects to the BS that provides the strongest signal strength, the serving BS can only be either the nearest BS in $\Phi_{g_l}$ or in $\Phi_{g_n}$. Lemma~\ref{lemma:backhaul_association} provides the probabilities of connecting to a LOS or a NLOS BS.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The probability $A_{b_\xi}$ that a UAV connects to a BS with LOS/NLOS backhaul link is
		\begin{equation}\small
		A_{b_\xi}=\int_{0}^{\infty}e^{-2\pi\lambda_b\int_{0}^{E_{b_\xi}(x)}\kappa_{b_{\bar{\xi}}}(t)t\mathrm{d}t}f_{s_{b_\xi}}(x)\mathrm{d}x,
		\label{eq:Aln}
		\vspace{-0.2cm}
		\end{equation}
		where $\bar{\xi}=n$ if $\xi=l$ and $\bar{\xi}=l$ if $\xi=n$. $E_{b_\xi}(x)=\sqrt{\left(\frac{C_{\bar{\xi}}}{C_{\xi}}\right)^{\frac{2}{\eta_{\bar{\xi}}}}(x^2+\Delta_h^2)^{\frac{\eta_{\xi}}{\eta_{\bar{\xi}}}}-\Delta_h^2}$, $\Delta_h=|h_u-h_g|$ and $f_{s_{b_\xi}}(x)=2\pi\lambda_b x \kappa_{b_\xi}(x)e^{-2\pi\lambda_b\int_{0}^{x}\kappa_{b_\xi}(r)r\mathrm{d}r}$. }
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:LOS_association}.}
	\end{IEEEproof}
	\label{lemma:backhaul_association}
	\vspace{-0.3cm}
\end{lemma}
The PDF of the distances $x_{b_l}$ and $x_{b_n}$ from the UAV to the serving BS in $\Phi_{g_l}$ and $\Phi_{g_n}$ are provided in Lemma~\ref{lemma:backhaul_distance}.
\vspace{-0.3cm}
\begin{lemma}
	\emph{Given that a UAV connects to a BS in $\Phi_{g_\xi}$ to get backhaul, the PDF of the distance $x_{b_\xi}$ to the serving BS is
		\vspace{-0.5cm}
		\begin{equation}\small
		f_{X_{b_\xi}}(x)=\frac{f_{s_{b_\xi}}(x)}{A_{b_\xi}}e^{-2\pi\lambda_b\int_{0}^{E_{b_\xi}(x)}\kappa_{b_{\bar{\xi}}}(t)t\mathrm{d}t},
		\label{eq:fsLN}
		\end{equation}
		where $\xi \in \{l,n\}$ indicates if the UAV connects to a LOS or a NLOS BS for backhaul support. $\bar{\xi}=n$ if $\xi=l$ and $\bar{\xi}=l$ if $\xi=n$. $A_{b_\xi}$, $f_{s_{b_\xi}}(\cdot)$ and $E_{b_\xi}(\cdot)$ are given in Lemma~\ref{lemma:backhaul_association}.
	}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:backhaul_distance}.}
	\end{IEEEproof}
	\label{lemma:backhaul_distance}
\end{lemma}
\vspace{-0.3cm}
Finally, the backhaul probability is presented in Theorem~\ref{theorem:backhaul_probability}.
\vspace{-0.4cm}
\begin{theorem}[Backhaul Probability]
	\emph{The backhaul probability $S(\tau_b)$ can be derived as
		\vspace{-0.2cm}
		\begin{equation}\small
		S(\tau_b)=A_{b_l}S_{l}(\tau_b)+A_{b_n}S_{n}(\tau_b),
		\label{eq:Stb}
		\vspace{-0.3cm}
		\end{equation}
		where $S_l(\tau_b)$ and $S_n(\tau_b)$ are the conditional backhaul probabilities given that the UAV is connected to a LOS/NLOS BS. $A_{b_l}$ and $A_{b_n}$ are the association probabilities. $S_{\xi}(\tau_b)$ is given as
		\vspace{-0.4cm}
		\begin{equation}\small
		\begin{aligned}
		&S_{\xi}(\tau_b)\approx\sum_{q=1}^{m_\xi}(-1)^{q+1} {m_\mathrm{\zeta}\choose q}\int_{0}^{\infty}e^{-\frac{q \gamma_\xi\left(x^2+\Delta_h^2\right)^{\frac{\eta_{\xi}}{2}}\tau_b\sigma_{b}^2}{P_{b}C_{\xi}G_{b,0}}-Q_{\xi}(q,x)-V_{\xi}(q,x)}f_{X_{b_\xi}}(x)\mathrm{d}x,
		\end{aligned}
		\label{eq:SLN}
		\end{equation}
		where $f_{X_{b_\xi}}(x)$ is given in (\ref{eq:fsLN}) and
		\vspace{-0.2cm}
		\begin{equation}\small
		\begin{aligned}
		Q_{\xi}(q,x)&=2\pi\lambda_b\sum_{k=1}^{4}p_k\int_{x}^{\infty}H\left(m_\xi,\frac{q\gamma_\xi\bar{G}_k\tau_b(x^2+\Delta_h^2)^{\frac{\eta_{\xi}}{2}}}{m_{\xi}(t^2+\Delta_h^2)^{\frac{\eta_{\xi}}{2}}}\right)\times \kappa_{b_{\xi}}(t)t\mathrm{d}t,
		\end{aligned}
		\end{equation}
		\vspace{-0.2cm}
		\begin{equation}\small
		\begin{aligned}
		V_{\xi}(q,x)&=2\pi\lambda_b\sum_{k=1}^{4}p_k\int_{E_{b_\xi}(x)}^{\infty}H\left(m_{\bar{\xi}},\frac{q C_{\bar{\xi}}\gamma_{\xi}\bar{G}_k\tau_b(x^2+\Delta_h^2)^{\frac{\eta_{\xi}}{2}}}{m_{\bar{\xi}}C_{\xi}(t^2+\Delta_h^2)^{\frac{\eta_{\bar{\xi}}}{2}}}\right)\times \kappa_{b_{\bar{\xi}}}(t)t\mathrm{d}t,
		\end{aligned}
		\end{equation}
		where $H(m,x)=1-1/(1+x)^m$, For $\xi \in \{l,n\}$, $\bar{\xi}=n$ if $\xi=l$ and $\bar{\xi}=l$ if $\xi=n$, $m_\xi$ and $\gamma_{\xi}=m_{\xi}(m_{\xi}!)^{-\frac{1}{m_{\xi}}}$ are the Nakagami-m small scale fading parameters; for $k\in\{1,2,3,4\}$, $\bar{G}_k=\frac{G_k}{G_{b,0}}$, $G_k$ and $p_k$ are defined in Table~\ref{table:antenna}.}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:backhaul_probability}.}
	\end{IEEEproof}
	\label{theorem:backhaul_probability}
\end{theorem}
\vspace{-0.4cm}
\subsubsection{Beamforming Alignment Errors}\label{sec:misalignment}
To optimize the backhaul connectivity between the UAV and the BS, the antennas must be perfectly aligned. However, perfect alignment in mmWave systems is challenging due to the required high directivity that leads to a narrow beamwidth of the main lobe. In the absence of beam-steering errors, the maximum directivity gain offered by directional beamforming is $G_{g}^{(\text{max})} G_{u}^{(\text{min})}$, where $G_{g}^{(\text{max})}$ and $G_{u}^{(\text{min})}$ are the gains of the main lobes of the BS and the UAV antennas, respectively. To account for the beamforming alignment error of the desired link, a beam steering error model similar to~\cite{Wildman} is considered. We denote by $\varphi_{s}^{*}$, $s \in \{g,u\}$, the angles corresponding to error-free beam-steering and by $\varepsilon_s$ the additive beam-steering errors. $\varepsilon_g$ and $\varepsilon_u$ are assumed to be randomly distributed, independent of each other and to have a symmetric distribution around $\varphi_{g}^{*}$ and $\varphi_{u}^{*}$, respectively. Let $F_{|\varepsilon_{s}|}(x) = \mathbb{P}(|\varepsilon_{s}|\leq x)$ be the cumulative distribution function (CDF) of $|\varepsilon_{\mathrm{s}}|$ for $s \in \{g,u\}$, and let $\bar{F}_{|\varepsilon_{s}|}(x)=\left(1-F_{|\varepsilon_{s}|}(x)\right)$. Then, the PMF of the directivity gain of the desired link $G_{b,0}$ is $f_{G_{b,0}}(g)=\left(f_{G_{g,0}}\otimes f_{G_{u,0}}\right)(g)$, where $f_{G_{g,0}}(.)$ and $f_{G_{u,0}}(.)$ are the PMFs of the directivity gains of the serving BS and the corresponding UAV, respectively and $\otimes$ is the convolution operator. $f_{G_{s,0}}(g)$ can be written as:
	\vspace{-0.1cm}
	\begin{equation}\small
	f_{G_{s,0}}(g)=F_{|\varepsilon_{\mathrm{s}}|}\left(\frac{\theta_s}{2}\right)\delta\left(g-G_{s}^{(\text{max})}\right)+\bar{F}_{|\varepsilon_{s}|}\left(\frac{\theta_s}{2}\right)\delta\left(g-G_{s}^{(\text{min})}\right)
	\vspace{-0.2cm}
	\end{equation}
	where $\delta(\cdot)$ is the Dirac delta function, and $G_{s}^{(\text{max})}$, $G_{s}^{(\text{min})}$ and $\theta_{s}$ are the respective gains of the main and side lobes and the beamwidth of the BS and UAV antennas, for $s \in \{g,u\}$, as illustrated earlier in (\ref{eq:gain}). Similarly to the directivity gain of the interfering links, the directivity gain of the desired link $G_{b,0}$ is a discrete random variable that takes the values $G_{k}$ with probability $p_{k,0}$, where $G_{k}$ and $p_{k_0}$ are defined in Table~\ref{table:antenna}. The backhaul probability can be expressed as
	\begin{equation}\small
	\begin{aligned}
	&S(\tau_b)=\mathbb{E}_{G_{b,0}}\left[S(\tau_b,G_{b,0})\right]=\int_{0}^{\infty}S(\tau_b,g)f_{G_{b,0}}(g)\mathrm{d}g\\
	&=F_{|\varepsilon_{g}|}(\theta_{g}/2)F_{|\varepsilon_{u}|}(\theta_{u}/2)S(\tau_b,G_{g}^{(\text{max})}G_{u}^{(\text{max})})+F_{|\varepsilon_{g}|}(\theta_{g}/2)\bar{F}_{|\varepsilon_{u}|}(\theta_{u}/2)S(\tau_b,G_{g}^{(\text{max})}G_{u}^{(\text{min})})\\
	&+\bar{F}_{|\varepsilon_{g}|}(\theta_{g}/2)F_{|\varepsilon_{u}|}(\theta_{u}/2)S(\tau_b,G_{g}^{(\text{min})}G_{u}^{(\text{max})})+\bar{F}_{|\varepsilon_{g}|}(\theta_{g}/2)\bar{F}_{|\varepsilon_{u}|}(\theta_{u}/2)S(\tau_b,G_{g}^{(\text{min})}G_{u}^{(\text{min})}).
	\end{aligned}\label{eq:misalignment}

	\end{equation}
\subsubsection{Conditional Coverage Probability and Distance Distribution}
The UE can connect to a LOS or a NLOS UAV depending on the association rule. The coverage probability $P_{cov,u_{\zeta}}$ when the UE associates to a LOS UAV ($\zeta=l$) or a NLOS UAV ($\zeta=n$) is given in Lemma~\ref{lemma:pcova}.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The conditional coverage probability $P_{cov,u_{\zeta}}$ given that UE connects to a $\zeta$-UAV is
		\begin{equation}\small
		P_{cov,u_{\zeta}}=S(\tau_b) \int_{h_u}^{w_p}\sum_{k=0}^{m_{\zeta}-1}\frac{(-s_{2_{\zeta}})^k}{k!}\left[\frac{\partial^k}{\partial s_{2_{\zeta}}^{k}}\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})\mathcal{L}_{\hat{I}_{u_{\zeta}}}(s_{2_{\zeta}})\right] f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})\mathrm{d}y_{u_{\zeta}},
		\label{eq:pcova}
		\end{equation}
		where $s_{2_{\zeta}}=\frac{m_{\zeta}\tau_a y_{u_{\zeta}}^{\eta_{\zeta}}}{P_{u}}$ and $f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})$ is the PDF of the conditional distance to the serving $\zeta$-UAV with $\zeta\in\{l,n\}$. $\mathcal{L}_{\hat{I}_{u_{\zeta}}}(s_{2_{\zeta}})$ and $\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})$ are the Laplace transforms of the aggregate interference of all the UAVs except the serving UAV $\hat{I}_{u_{\zeta}}$ and of all the BSs $I_{g_{\zeta}}$.}
	\begin{IEEEproof}
		\emph{See Appendix~\ref{app:pcova}.}
	\end{IEEEproof}
	\label{lemma:pcova}
	\vspace{-0.3cm}
\end{lemma}
As noted in (\ref{eq:pcova}), the Laplace transforms of the interference terms $\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})$ and $\mathcal{L}_{\hat{I}_{u_{\zeta}}}(s_{2_{\zeta}})$ and the PDF of the conditional distance to the serving UAV $f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})$ must be computed to obtain the final expression of $P_{cov,u_{\zeta}}$. Lemma~\ref{lemma:UE_UAV_laplace} presents these Laplace transforms as follows.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The Laplace transform of the interference $I_{g_{\zeta}}$ of all BSs when the UE associates to a $\zeta$-UAV with $\zeta\in\{l,n\}$ is given as
		\begin{equation}\small
		\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})=\exp\left[\frac{-2\pi\lambda_g s_{2_{\zeta}} P_{g}\left(E_{u_{\zeta}}^{2}(y_{u_{\zeta}})+h_{g}^2\right)}{\left(\eta_{g}-2\right)\left(\left(E_{u_{\zeta}}^{2}(y_{u_{\zeta}})+h_{g}^2\right)^{\frac{\eta_{g}}{2}}+s_{2_{\zeta}} P_{g}\right)}{_{2}}F_{1}\left(1,1;2-\frac{2}{\eta_{g}};\frac{1}{1+\frac{\left(E_{u_{\zeta}}^{2}(y_{u_{\zeta}})+h_{g}^2\right)^{\frac{\eta_{g}}{2}}}{s_{2_{\zeta}} P_{g}}}\right)\right].
		\label{eq:laplace_Ig}
		\end{equation}
		where $E_{u_{\zeta}}(y_{u_{\zeta}})$ is given in (\ref{eq:Eg}) and ${_{2}}F_{1}(\cdot)$ is the Gauss hyper geometric function. The Laplace transform of the interference $\hat{I}_{u_{\zeta}}$ from all UAVs except the serving UAV is given as
		\begin{equation}\small
		\begin{aligned}
		\mathcal{L}_{\hat{I}_{u_{\zeta}}}(s_{2_{\zeta}})&=\left[\frac{1}{\int_{y_{u_{\zeta}}}^{w_{p}}f_{W}(w)\kappa_{u_{\zeta}}(w)\mathrm{d}w+\int_{E_{\zeta\bar{\zeta}}\left(y_{u_{\zeta}}\right)}^{w_{p}}f_{W}(w)\kappa_{u_{\bar{\zeta}}}(w)\mathrm{d}w}\left(\int_{y_{u_{\zeta}}}^{w_{p}}\left(1+\frac{s_{2_{\zeta}}P_{u}v^{-\eta_{\zeta}}}{m_{\zeta}}\right)^{-m_{\zeta}}\right.\right.\\
		&\times f_{W}(v)\kappa_{u_\zeta}(v)\mathrm{d}v+\int_{E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})}^{w_{p}}\left(1+\frac{s_{2_{\zeta}}P_{u}v^{-\eta_{\bar{\zeta}}}}{m_{\bar{\zeta}}}\right)^{-m_{\bar{\zeta}}}f_{W}(v)\kappa_{u_{\bar{\zeta}}}(v)\mathrm{d}v\bigg)\bigg]^{N_u-1}.
		\end{aligned}
		\end{equation}}
	\vspace{-0.4cm}
	\begin{IEEEproof}
		\emph{When the UE associates to a $\zeta$-UAV, the BSs are located further than $E_{u_{\zeta}}(y_{u_{\zeta}})$, the interfering $\zeta$-UAV are located further than $y_{u_{\zeta}}$ and the $\bar{\zeta}$-UAV exist beyond $E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})$, where $y_{u_{\zeta}}$ is the distance from the UE to the serving UAV. Thus, $\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})$ is obtained by following a similar approach to Lemma~\ref{lemma:UE_BS_laplace} while replacing the lower limit of the integral in (\ref{eq:prooflemma3}) by $E_{u_{\zeta}}(y_{u_{\zeta}})$. The interference from all UAVs except the serving UAV can be expressed as $\sum_{y_i\in\Phi_{u}\backslash y_{u_{\zeta}}}I_{u,y_j}$, where $I_{u,y_j}$ is the interference from the $j$-th UAV and can be obtained by following a similar approach to Lemma~\ref{lemma:UE_BS_laplace} as $\mathcal{L}_{\hat{I}_{u,y_{j}}}(s_{2_{\zeta}})=\left(\mathbb{E}_{I_{u,y_j}}\left[\exp\left(-s_{2_\zeta} I_{u,y_j}\right)\right]\right)^{N_u-1}$. The final expression of $\mathcal{L}_{\hat{I}_{u_\zeta}}(s_{2_{\zeta}})$ is obtained by replacing the lower limits of the integrals by $y_{u_{\zeta}}$ and $E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})$ and by plugging $f_{V_{u_l}}(v,y_{u_\zeta})$ and $f_{V_{u_n}}(v,E_{\zeta\bar{\zeta}}(y_{u_{\zeta}}))$ in $\mathbb{E}_{I_{u,y_j}}\left[\exp\left(-s_{2_{\zeta}} I_{u,y_j}\right)\right]$ in (\ref{eq:Euai}).}
	\end{IEEEproof}
	\label{lemma:UE_UAV_laplace}
\end{lemma}
\vspace{-0.2cm}
Next, we derive the PDF of the distance from the UE to the serving UAV $y_{u_{\zeta}}$ in Lemma~\ref{lemma:fxazeta}.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The PDF of the distance $y_{u_{\zeta}}$ to the serving $\zeta$-UAV ($\zeta \in \{l,n\}$) is given as
		\begin{equation}\small
		\begin{aligned}
		f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})&=\frac{N_u}{A_{u_{\zeta}}}f_{W}(y_{u_{\zeta}})\kappa_{u_\zeta}(y_{u_{\zeta}})\exp\left(-\pi\lambda_g E_{u_{\zeta}}(y_{u_{\zeta}})^2\right)\left(\int_{y_{u_{\zeta}}}^{w_p}f_{W}(w)\kappa_{u_\zeta}(w)\mathrm{d}w\right.\\
		&\left.+\int_{E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})}^{w_p}f_{W}(w)\kappa_{u_{\bar{\zeta}}}(w)\mathrm{d}w\right)^{N_u-1},
		\end{aligned}
		\label{eq:fxa}
		\end{equation}
		where $f_W(\cdot)$ and $A_{u_\zeta}$ are given in Lemma~\ref{lemma:coverage_association}. $E_{u_{\zeta}}(\cdot)$ and $E_{\zeta\bar{\zeta}}(\cdot)$ are given in (\ref{eq:Eg}).}
	\begin{IEEEproof}
		\emph{This proof follows a similar approach to Lemma~\ref{lemma:fxg}.}
	\end{IEEEproof}
	\label{lemma:fxazeta}
\end{lemma}
\vspace{-1cm}
\subsection{Overall Coverage Probability}
\vspace{-0.2cm}
After deriving the association probabilities, the conditional distance distributions and the conditional coverage probabilities, the overall coverage probability is obtained through the law of total probability as
$P_{cov}=A_{u_l} P_{cov,u_l}+A_{u_n} P_{cov,u_n}+A_g P_{cov,g}$.
\vspace{-0.5cm}
\section{Backhaul Aware Transmission}\label{sec:aware}
\vspace{-0.1cm}
For the backhaul aware transmission scenario, the UE activity is dependent on the BS-UAV instantaneous backhaul link in case of aerial coverage. That is, the UAV is aware of the backhaul status as it checks the integrity of packets received from the BS before transmitting them to the UE. If the backhaul $\mathrm{SINR}$ satisfies the threshold, the UAV can serve the UE successfully. When the backhaul is down, the UAV refrains its transmission to the UE. The UE does not change its association,\footnote{Changing association based on instantaneous link status may lead to undesirable excessive handovers known as the ping pong effect.} but rather goes to a service failure event in this time slot. As result, the backhaul aware transmission scheme has the following two main advantages, i) the established DL transmissions are subject to the access link $\mathrm{SIR}$ only, since UAVs have the ability to refrain from transmitting corrupted packets, and ii) the interference is relieved because only successfully backhauled UAVs contribute to it, thereby improving the overall coverage probability. Such advantages come at the cost of extra processing at the UAVs because they check the validity of the packets before transmitting to UE. To characterize the coverage probability of the backhaul aware scheme, we follow the same steps of Section~\ref{sec:unaware}. In particular, we start by deriving the transmission probabilities, then we find the conditional coverage probabilities, then we obtain the distance distributions, and finally we provide the unified expression for the coverage probability.
\vspace{-1.4cm}
\subsection{Transmission Probabilities}\label{sec:aware_association}
\vspace{-0.2cm}
In the context of backhaul aware transmission, we define the transmission probability as the probability that the UE is being served by a BS or a successfully backhauled LOS/NLOS UAV. Note that in the case where the UAV-UE access link is better than the BS-UE link and the UAV backhaul link is not successful, the UE is considered in a service failure state and fails to connect. Hence, the service failure event here refers to the event where the UE is not served in the DL in a time slot not as the event that DL is established and a decoding error occurs. It is important to note that an established DL connection is different from UE to UAV association, where the latter is independent from the backhaul link state. When the UE is served by a terrestrial BS, the transmission probability is expressed as $\tilde{A}_g=\mathbb{P}[s=g]$. Note that the transmission probability when the UE connects to a BS is the same as the association probability to a BS for the unaware transmission scenario and is presented in (\ref{eq:Ag}). For a UAV to transmit to the UE, the backhaul link of this UAV should be successful. The UAV transmission probability is the joint probability of two events and can be expressed as $\tilde{A}_{u_{\zeta}}=\mathbb{P}[s=u_{\zeta}, \mathrm{SINR}>\tau_b]$, where the first term corresponds to the UE-UAV association rule $\left(\zeta \in\{l,n\}\right)$ and the second term to the BS-UAV backhaul link. The independence assumption is also considered and validated in Section~\ref{sec:results}. Thus, the UAV transmission probability can be expressed as $\tilde{A}_{u_{\zeta}}=\mathbb{P}[s=u_{\zeta}]\times\mathbb{P}[\mathrm{SINR}>\tau_b]$, where $\mathbb{P}[s=u_{\zeta}]=A_{u_{\zeta}}$ is the UAV association probability for the backhaul unaware scenario derived in Lemma~\ref{lemma:coverage_association} and $\mathbb{P}[\mathrm{SINR}>\tau_b]=S(\tau_b)$ is the backhaul probability given in (\ref{eq:Stb}). The service failure probability is defined as $\tilde{A}_{f}=\mathbb{P}[s=u_l , \mathrm{SINR}<\tau_b \| s=u_n, \mathrm{SINR}<\tau_b]$ and derived as
\vspace{-0.2cm}
\begin{equation}
\begin{aligned}
\tilde{A}_{f} &=\mathbb{P}\left[s=u_l \right]\mathbb{P}\left[\mathrm{SINR} < \tau_b \right]+\mathbb{P}\left[s=u_n \right]\mathbb{P}\left[\mathrm{SINR} < \tau_b \right]=A_{u_l} \left(1-S(\tau_b)\right)+A_{u_n} \left(1-S(\tau_b)\right),
\end{aligned}
\vspace{-0.5cm}
\end{equation}
where $A_{u_l}$ and $A_{u_n}$ are derived in Lemma~\ref{lemma:coverage_association} and $S(\tau_b)$ in Theorem~\ref{theorem:backhaul_probability}.
\vspace{-0.4cm}
\subsection{Conditional and Overall Coverage Probabilities}\label{sec:aware_conditional}
\vspace{-0.2cm}
The conditional coverage probability given that the UE is served by a BS is expressed as $\tilde{P}_{cov,g}=\mathbb{P}[\mathrm{SIR}\geq\tau_a|s=g]$. Due to the dependency of the UAV transmission on the quality of the backhaul link, the aggregate interference for the backhaul aware scenario $\tilde{I}_{agg,g}$ includes the interference from only the UAVs with successful backhaul links denoted as $I_{\tilde{u}}$ and from all BSs except the serving BS denoted as $\hat{I}_{g}$. An expression for $\tilde{P}_{cov,g}$ is given in Lemma~\ref{lemma:pcovglemma}.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The conditional coverage probability given that the UE is served by a BS for the backhaul aware scenario is given as
		\begin{equation}\small
		\tilde{P}_{cov,g}=\int_{0}^{E_u}\mathcal{L}_{\hat{I}_g}(s_1)\mathcal{L}_{I_{\tilde{u}}}(s_1)f_{X_g}(x_g)\mathrm{d}x_{g},
		\label{eq:pcovg_aware}
		\end{equation}
		where $s_{1}=\frac{\tau_a(x_g^2+h_g^2)^{\frac{\eta_g}{2}}}{P_{g}}$, $E_u = \mathrm{max}\left(E_{u_l}(w_p),E_{u_n}(w_{p})\right)$, $\mathcal{L}_{\hat{I}_g}(s_{1})$ is given in (\ref{eq:Ligb}) and $f_{X_g}(x_g)$ in (\ref{eq:fxg}). $\mathcal{L}_{I_{\tilde{u}}}(s_1)$ is the Laplace transform of the aggregate interference from the UAVs that have successful backhaul links with ground BSs.}
	\begin{IEEEproof}
		\emph{This proof follows a similar approach to Lemma~\ref{lemma:pcovg}.
		}
	\end{IEEEproof}
	\label{lemma:pcovglemma}
	\vspace{-0.3cm}
\end{lemma}
A lower bound for the Laplace transform $\mathcal{L}_{I_{\tilde{u}}}(s_1)$ is derived in Lemma~\ref{lemma:laplace_aware}.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The Laplace transform of the aggregate interference $I_{\tilde{u}}$ from all the UAVs that have a successful backhaul link when the UE is served by a BS can be lower bounded by
		\begin{equation}\small
		\begin{aligned}
		\mathcal{L}_{I_{\tilde{u}}}(s_{1})&\geq\left[\frac{1}{\int_{E_{g_l}(x_{g})}^{w_{p}}f_{W}(w)\kappa_{u_l}(w)\mathrm{d}w+\int_{E_{g_n}(x_{g})}^{w_{p}}f_{W}(w)\kappa_{u_n}(w)\mathrm{d}w}\left(\int_{E_{g_l}(x_{g})}^{w_{p}}\left(1+\frac{s_{1}P_{u}v^{-\eta_l}}{m_l}\right)^{-m_l}\right.\right.\\
		&\times f_{W}(v)\kappa_{u_l}(v)\mathrm{d}v+\int_{E_{g_n}(x_{g})}^{w_{p}}\left(1+\frac{s_{1}P_{u}v^{-\eta_n}}{m_n}\right)^{-m_n} f_{W}(v)\kappa_{u_n}(v)\mathrm{d}v\bigg)\bigg]^{N_{u}\times S(\tau_b)},
		\end{aligned}
		\label{eq:Lia_aware}
		\end{equation}
		where $E_{g_l}(\cdot)$ and $E_{g_n}(\cdot)$ are given in (\ref{eq:Eg}), $f_{W}(\cdot)$ and $w_p$ are given in Lemma~\ref{lemma:coverage_association}. }
	\begin{IEEEproof}
		\emph{We denote $\tilde{N_u}$ the number of UAVs that have successful backhaul links. $I_{\tilde{u}}$ is given as $\sum_{y_{j\in\tilde{\Phi_{u}}}}I_{u,y_j}$ where $I_{u,y_j}$ denotes the interference from the $j$-th UAV and $\tilde{\Phi_{u}}$ is the set of UAVs with successful backhaul links. The Laplace transform of $I_{\tilde{u}}$ can be calculated as
			\begin{equation}\small
			\begin{aligned}
			&\mathcal{L}_{I_{\tilde{u}}}(s_{1})=\mathbb{E}_{I_{\tilde{u}}}\left[e^{-s_{1}I_{\tilde{u}}}\right]=\mathbb{E}_{I_{\tilde{u}}}\left[\exp\left(-s_{1}\sum_{j=1}^{\tilde{N_u}}I_{u,y_j}\right)\right]
			\stackrel{(a)}{=}\mathbb{E}_{\tilde{N_u}}\left[\prod_{j=1}^{\tilde{N_u}} \mathbb{E}_{I_{u,y_j}}\left[\exp\left(-s_{1}I_{u,y_j}\right)\right]\right]\\
			&\stackrel{(b)}{=}\mathbb{E}_{\tilde{N_u}}\left[\left(\mathbb{E}_{I_{u,y_j}}\left[\exp\left(-s_1 I_{u,y_j}\right)\right]\right)^{\tilde{N_u}}\right]\stackrel{(c)}{\geq}\left(\mathbb{E}_{I_{u,y_j}}\left[\exp\left(-s_1 I_{u,y_j}\right)\right]\right)^{\mathbb{E}\left[\tilde{N_u}\right]},
			\end{aligned}
			\end{equation}
			where (a) follows from the independent and identically distributed (iid) fading gains and from their independence of the interferers distances in the expression of the interference $I_{u,y_j}$. (b) follows from the iid distribution of the interferers UAVs. (c) follows from the Jensen's inequality that gives a lower bound for $\mathcal{L}_{I_{\tilde{u}}}(s_{1})$. The number of UAVs that have successful backhaul links $\tilde{N_u}$ follows a binomial distribution with backhaul probability $S(\tau_b)$ and mean $\mathbb{E}[\tilde{N_u}]=N_u S(\tau_b)$. The remaining proof follows the procedure adopted in Lemma~\ref{lemma:UE_BS_laplace}.}
	\end{IEEEproof}
	\label{lemma:laplace_aware}
	\vspace{-0.3cm}
\end{lemma}
Similar to the UE-BS conditional coverage probability, the conditional coverage probability given that the UE is served by a $\zeta$-UAV is expressed as $\tilde{P}_{cov,u_{\zeta}}=\mathbb{P}[\mathrm{SIR}\geq\tau_a|s=u_{\zeta}]$. The aggregate interference $\tilde{I}_{agg,u_{\zeta}}$ includes the interference from the UAVs with successful backhaul links except the serving UAV denoted as $\hat{I}_{\tilde{u_{\zeta}}}$ and from all BSs denoted as $I_{g_{\zeta}}$. An expression for $\tilde{P}_{cov,u_{\zeta}}$ is given in the following lemma.
\vspace{-0.3cm}
\begin{lemma}
	\emph{The coverage probability $P_{cov,u_{\zeta}}$ given that the UE is served by a $\zeta$-UAV is
		\begin{equation}\small
		\tilde{P}_{cov,u_{\zeta}}= \int_{h_u}^{w_p}\sum_{k=0}^{m_{\zeta}-1}\frac{(-s_{2_{\zeta}})^k}{k!}\left[\frac{\partial^k}{\partial s_{2_{\zeta}}^{k}}\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})\mathcal{L}_{\hat{I}_{\tilde{u_{\zeta}}}}(s_{2_{\zeta}})\right] f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})\mathrm{d}y_{u_{\zeta}},
		\label{eq:pcova_aware}
		\end{equation}
		where $\zeta\in\{l,n\}$, $s_{2_{\zeta}}=\frac{m_{u_{\zeta}}\tau_a y_{u_{\zeta}}^{\eta_{u_{\zeta}}}}{P_{u}}$. $f_{Y_{u_{\zeta}}}(y_{u_{\zeta}})$ is given in (\ref{eq:fxa}) and $\mathcal{L}_{I_{g_{\zeta}}}(s_{2_{\zeta}})$ in (\ref{eq:laplace_Ig}). $\mathcal{L}_{\hat{I}_{\tilde{u_{\zeta}}}}(s_{2_{\zeta}})$ is the Laplace transform of the aggregate interference of all the UAVs that have successful backhaul links except the serving UAV and is given as
		\begin{equation}\small
		\begin{aligned}
		\mathcal{L}_{\hat{I}_{\tilde{u_{\zeta}}}}(s_{2_{\zeta}})&\geq \left(\frac{1}{\int_{y_{u_{\zeta}}}^{w_{p}}f_{W}(w)\kappa_{u_\zeta}(w)\mathrm{d}w+\int_{E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})}^{w_{p}}f_{W}(w)\kappa_{u_{\bar{\zeta}}}(w)\mathrm{d}w}\left(\int_{y_{u_{\zeta}}}^{w_{p}}\left(1+\frac{s_{2_{\zeta}}P_{u}v^{-\eta_{\zeta}}}{m_{\zeta}}\right)^{-m_{\zeta}}\right.\right.\\
		&\times f_{W}(v)\kappa_{u_\zeta}(v)\mathrm{d}v+\int_{E_{\zeta\bar{\zeta}}(y_{u_{\zeta}})}^{w_{p}}\left(1+\frac{s_{2_{\zeta}}P_{u}v^{-\eta_{\bar{\zeta}}}}{m_{\bar{\zeta}}}\right)^{-m_{\bar{\zeta}}}f_{W}(v)\kappa_{u_{\bar{\zeta}}}(v)\mathrm{d}v\bigg)\bigg)^{N_u S(\tau_b)-1}.
		\end{aligned}
		\end{equation}
		\vspace{-0.5cm}
	}
	\begin{IEEEproof}
		\emph{The proof of this lemma follows a similar approach to Lemma~\ref{lemma:UE_BS_laplace} and Lemma~\ref{lemma:laplace_aware}.}
	\end{IEEEproof}
	\label{lemma:pcova_aware}
	\vspace{-0.2cm}
\end{lemma}
The overall coverage probability for the backhaul aware transmission scenario can be derived through the law of total probability as $\tilde{P}_{cov}=\tilde{A}_{u_l} \tilde{P}_{cov,u_l}+\tilde{A}_{u_n} \tilde{P}_{cov,u_n}+\tilde{A}_g \tilde{P}_{cov,g}$, where $\tilde{A}_{u_l}$, $\tilde{A}_{u_n}$ and $\tilde{A}_g$ are given in Section~\ref{sec:aware_association} and $\tilde{P}_{cov,u_l}$, $\tilde{P}_{cov,u_n}$ and $\tilde{P}_{cov,g}$ are given in Section~\ref{sec:aware_conditional}.
\vspace{-0.5cm}
\section{Numerical results and Discussions}\label{sec:results}
In this section, we verify the analytical analysis against Monte-Carlo simulations and highlight the impact of different system parameters. Unless it is stated explicitly otherwise, we assume that the UE is at the origin ($v_0=0$) and we employ the simulation parameters listed in Table~\ref{table:parameters}. Assuming dense urban environment, we set the environmental parameters $a = 0.136$ and $b = 11.95$ as in~\cite{Saad} to determine the LOS probability on the UAV-UE access links. The dense urban environment parameters of the LOS probability model adopted for the BS-UAV backhaul links are given in~\cite{Cherif} as $h_{g}=25$~m, $c=1$, $d=0.106$ and $e=1$.
\begin{table*}[t]
	\centering
	\caption{Simulation Parameters.}
	\vspace{-0.3cm}
	\begin{tabular}{c|c||c|c||c|c}
		\hline\hline
		\textbf{Parameter} & \textbf{Value} & \textbf{Parameter} & \textbf{Value}  & \textbf{Parameter} & \textbf{Value} \\\hline
		($P_{g}$, $P_{b}$, $P_{u}$) & ($20$, $10$, $1$)~W & ($\eta_{g}$, $\eta_l$, $\eta_n$) & ($4$ , $2.5$, $4$) & ($\theta_{g}, \theta_{u}$) & ($20^{\circ}$, $20^{\circ}$) \\\hline
		$r_{c}$ & $1000$~m & ($m_l$, $m_n$) &($3$, $2$) & ($a$, $b$) & ($11.95$, $0.136$) \\\hline
		($h_{g}$, $h_{u}$) & ($25$, $100$)~m & ($C_l$, $C_n$) & ($-69.8$, $-69.8$)~dB & ($c$, $d$, $e$) & ($1$, $0.106$, $1$) \\\hline
		$\lambda_{g}$ & $10$~BS/$\mathrm{km}^2$ &($G_{g}^{(\text{max})}$, $G_{g}^{(\text{min})}$) & ($18$, $-2$)~dB & $\sigma_{b}^2$ & $4\cdot10^{-11}$~W \\\hline
		$N_{u}$, $\delta_{b}$ & $10$, $1$ & ($G_{u}^{(\text{max})}$, $G_{u}^{(\text{min})}$) & ($18$, $-2$)~dB & ($\tau_a$, $\tau_b$) & ($0$, $10$)~dB \\\hline\hline
	\end{tabular}
	\label{table:parameters}
	\vspace{-0.5cm}
\end{table*}
      \begin{figure*}[t]
   	\begin{subfigure}[t]{0.33\textwidth}
   	\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
   			\input{Figures/Fig5a.tex}
   			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Backhaul probability.}
		\label{fig:backhaul_ha}
   	\end{subfigure}
   	~
   	\begin{subfigure}[t]{0.33\textwidth}
   		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
   			\input{Figures/Fig5b.tex}
   			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Association probabilities.}
		\label{fig:association}
   	\end{subfigure}
   	~
   	\begin{subfigure}[t]{0.33\textwidth}
   		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
   			\input{Figures/Fig5c.tex}
   			\end{tikzpicture}}
	\vspace{-0.2cm}
	\caption{Coverage probability.}
	\label{fig:aware_ha_beta}
   	\end{subfigure}
	\vspace{-0.4cm}
	\caption{Validation of analytical results and impact of UAVs height for backhaul aware.}
	\vspace{-1cm}
   \end{figure*}

\vspace{-0.8cm}
\subsection{Numerical results}
\vspace{-0.25cm}
We analyze the impact of the UAV height on the capability of the UAV to get a good backhaul link in Fig.~\ref{fig:backhaul_ha}. The solid lines represent the analytical results and the markers represent the simulation results. It is noted that the analytical results match perfectly with the simulations, which validates the proposed model. When a UAV hovers at higher altitude, its distance to the serving BS increases, which causes a degradation in the received signal power. However, a higher altitude means a higher LOS probability. For very low altitudes above the terrestrial BSs, the LOS probability is the dominant factor, and hence, the backhaul probability improves. As the UAV increases its altitude, the path-loss dominates and the backhaul probability degrades. An optimal height can be noted in Fig.~\ref{fig:backhaul_ha}, which maximizes the backhaul probability.

In Fig.~\ref{fig:association}, we show the aerial and terrestrial association probabilities in addition to the service failure probability for the backhaul aware scenario as function of the height of the UAVs. In the low heights region, the UE tends to associate more with UAVs and less with terrestrial BSs as the height of these UAVs increases. This is due to the fact that, as the height increases, more UAVs come within LOS conditions with respect to the UE. However, with a further increase in height, the UAV association probability starts to decrease because of the significant path-loss caused by the increasing distances between the UE and the UAVs. At the same time, the backhaul probability degrades which causes the increase of the service failure probability that corresponds to a good UAV-UE access connection with bad backhauling. The existence of UAVs in this case is harming the performance of the UE as backhauling is not guaranteed at high heights. With a further increase in height, the service failure probability starts to decrease because of the significant path-loss encountered on the UAV-UE access links. In this high height regime, the role of UAV becomes minor and the UE sticks to terrestrial connections.

Fig.~\ref{fig:aware_ha_beta} demonstrates how the increase of the UAVs height affects the coverage probability for different $\mathrm{SIR}$ thresholds in the backhaul aware scenario. For low heights, the coverage probability increases as the UAV height increases. This is due to the fact that, as the height increases, more UAVs will have successful backhaul links. For higher altitudes, the backhaul link deteriorates and the distances from the UE to its serving device and interfering UAVs increase, and so the $\mathrm{SIR}$ and the coverage probability decrease. This corresponds to the increase in the service failure probability as seen in Fig.~\ref{fig:association} which reflects the existence of some UEs with good UAV-UE access links and bad backhaul links. With a further increase in height, the coverage probability starts to increase again until a certain limit after which it stabilizes. For these high heights, there is no point of using UAVs and all the UEs are served with the terrestrial network. Fig.~\ref{fig:aware_ha_beta} reveals that the optimal UAV height is related to the value of the $\mathrm{SIR}$ threshold $\tau_a$. For higher thresholds, the UAVs must be deployed at lower height to maximize the coverage probability.

\begin{figure*}[t]
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig6a.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Number of UAVs.}
		\label{fig:coverage_aware_unaware100}
	\end{subfigure}
	~
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig6b.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Beamforming misalignment error.}
		\label{fig:misalignment_error}
	\end{subfigure}
	~
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig6c.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Density of terrestrial BSs.}
		\label{fig:coverage_lambda}
	\end{subfigure}
	\vspace{-0.4cm}
	\caption{Coverage probability as function of number of UAV (for backhaul unaware and aware transmissions), beamforming misalignment error and density of BSs.}
	\vspace{-1cm}
\end{figure*}
In Fig.~\ref{fig:coverage_aware_unaware100}, we illustrate how the number of UAVs impacts the coverage probability for $\tau_a=0$~dB and $h_u=100$~m for different backhaul thresholds~$\tau_b$. Solid and dashed lines present the results of the backhaul aware and unaware transmission scenarios, respectively. It is clearly seen that increasing the number of UAVs first improves the coverage probability but then deteriorates it for both scenarios. This is because adding more UAVs has both a positive and a negative impact. First, a higher UAVs number increases the average received power from the serving UAV since the distance between the UE and the UAV is lower and the LOS probability is higher. However, the added UAVs impose higher interference levels. Due to the bounded region in which the UAVs are distributed and since the UE-UAV distance is always $\geq h_u$ and the LOS probability $\leq 1$, the power of the received signal reaches a maximum when a sufficient number of UAVs is deployed. On the other side, the interference remains increasing when more UAVs are deployed and the negative effect becomes the dominant factor. Fig.~\ref{fig:coverage_aware_unaware100} also reveals that there exists an optimal number of UAVs that must be deployed to optimize the UE service. Furthermore, Fig.~\ref{fig:coverage_aware_unaware100} shows that the awareness of the backhaul link improves the coverage probability compared to the unaware scenario. This is due to the exclusion of the UAVs that do not have a successful backhaul link from the interfering UAVs. As the number of UAVs in the network increases, the gap between the aware and unaware coverage results becomes larger. It is worth to note here that, even though a service failure event is considered for backhaul aware transmission, its performance remains superior to the unaware scenario for different network parameters.

To evaluate the impact of the beamforming misalignment error on the backhaul link, we introduce the alignment error in the backhaul probability expression in (\ref{eq:misalignment}). We assume that the additive beam steering errors $\varepsilon_{g}$ and $\varepsilon_{u}$ for both BSs and UAVs follow a Gaussian distribution with zero mean and variances $\sigma_{g}^2$ and $\sigma_{u}^2$, respectively. Thus, $\lvert \varepsilon_{s}\rvert$ follows a half-normal distribution with $F_{\lvert\varepsilon_{s}\rvert}=\mathrm{erf}\left(\frac{x}{\sqrt{2}\sigma_{s}}\right)$ and $\bar{F}_{\lvert\varepsilon_{s}\rvert}(x)=1-F_{\lvert\varepsilon_{s}\rvert}(x)$ where $s\in\{g,u\}$ and $\mathrm{erf}(\cdot)$ denotes the error function. Fig.~\ref{fig:misalignment_error} presents the coverage probability as function of the threshold $\tau_b$ for different values of the misalignment error variance ($\sigma_{u}=\sigma_{g}$). We can clearly notice that the quality of the backhaul link has a direct impact on the UE experience. When the BS and UAV antennas are not aligned on the backhaul link, the coverage probability deteriorates. Thus, overlooking the backhaul effect can lead to misleading coverage insights. The obtained results validate the importance of our work compared to the works that consider guaranteed backhaul links.

Fig.~\ref{fig:coverage_lambda} shows the coverage probability $P_{cov}$ as function of $\lambda_{g}$ and the $h_{u}$ considering backhaul aware transmission. When the UAVs hover at low altitudes, increasing the number of BSs decreases $P_{cov}$. Higher densities of BSs result in a higher backhaul probability since the distance separating the UAV from its serving BS decreases, which increases the average received signal power and consequentially enhances the LOS probability. At the same time, at lower heights, the UE tends to associate more to UAVs, and since backhaul aware transmission is considered, more successfully backhauled UAVs means high interference levels and lower coverage probability. At the optimal height, the UEs are striking a balance between good backhauling and association with UAVs. As the height further increases, the coverage probability starts decreasing and reaches a break even point after which adding more terrestrial BSs improves the coverage probability. At this break even point, connecting to a UAV or a terrestrial BS would give the same coverage performance. This analysis is valid until a certain value of the density of BSs, below which insufficient backhauling is encountered. This is clearly noticed for low BSs densities as the UAV has less BSs to align beam with, and initiate backhaul connection. For high heights, the impact of the aerial network is minimal and the UE associates with the terrestrial network. Thus, adding more BSs would decrease the distance separating the UE from its serving BS, which increases the average received signal power and consequentially the coverage probability.

\begin{figure*}[t]
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig7a.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Backhaul probability.}
		\label{fig:backhaul_fraction}
	\end{subfigure}
	~
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig7b.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Coverage probability.}
		\label{fig:coverage_fraction}
	\end{subfigure}
	~
	\begin{subfigure}[t]{0.33\textwidth}
		\resizebox{0.97\columnwidth}{!}{\begin{tikzpicture}[thick,scale=1, every node/.style={scale=1.3},font=\Huge]
			\input{Figures/Fig7c.tex}
			\end{tikzpicture}}
		\vspace{-0.2cm}
		\caption{Aware vs. Instantaneous.}
		\label{fig:outage_no_outage}
	\end{subfigure}
	\vspace{-0.4cm}
	\caption{Backhaul and coverage probabilities as function of the fraction of BSs with backhaul and comparison between backhaul aware transmission and backhaul aware instantaneous association.}
	\vspace{-1cm}
\end{figure*}

Fig.~\ref{fig:backhaul_fraction} and Fig.~\ref{fig:coverage_fraction} show the impact of having a fraction~$\delta_{b}\leq1$ of backhaul-enabled BSs on both the backhaul and the coverage probabilities, respectively. For low heights, adding more backhaul-enabled BSs increases both the backhaul and the coverage probabilities. This is due to the fact that, as more BSs provide backhauling, the distance separating the UAV from its serving BS decreases, which increases the average received power and consequentially enhances the backhaul probability. As this probability increases, and since at low heights, the UE associates more to UAVs, the coverage probability increases. Thus, at low heights, all the terrestrial BSs must be able to provide backhaul capabilities to optimize the coverage. At higher heights, the backhaul and coverage probabilities increase with the increase of $\delta_{b}$ until a certain limit after which the two probabilities stabilize. Thus, adding more backhaul-enabled BSs will not bring any benefit and only a fraction is needed to achieve the optimal coverage. As the UAVs height further increases, the backhaul probability starts to deteriorate with the increase of $\delta_{b}$. This is caused by the increased interference levels added to the high path loss and LOS probability encountered on the backhaul links. The coverage probability, however, is not affected with the increase of $\delta_{b}$ since at high heights, the UE tends to associate more to the terrestrial network. The obtained results suggest that, a good coverage probability can be achieved when a fraction of the terrestrial BSs are backhaul-enabled. The fraction of backhaul-enabled BSs that optimizes the coverage of UEs is highly dependent on the height of the UAVs.

Finally, we compare the proposed backhaul aware transmission scheme with the instantaneous backhaul aware association scenario. In the latter scenario, when the UE-UAV access link quality is better than the UE-BS link but the UAV does not have a successful backhaul link, the UE is not considered in a service failure state. On the contrary, the UE changes its association to the best terrestrial BS alternative. The instantaneous backhaul aware association scenario models the case when UAVs are used to complement existing networks by providing additional capacity to hotspots areas and network coverage in hard to reach areas. In such scenario, when the UAV connection does not have any added value, the UE operates through the terrestrial connection. The backhaul aware transmission scheme is more applied in unexpected natural disasters. In such scenarios, the existing networks can be damaged. Thus, if the connection to UAV fails and the UE-BS access link quality is bad, the UE cannot connect to the network.

Fig.~\ref{fig:outage_no_outage} presents the comparison between backhaul aware transmission and the simulation results of the instantaneous backhaul aware transmission scheme as a function of the number and the height of UAVs. Note here that the analytical derivations can be conducted by thinning the BPP of UAVs according to the backhaul probability and then conducting the backhaul unaware association analysis. It is clearly seen that, as the number of UAVs increases, the coverage probability of backhaul aware instantaneous association exceeds that of backhaul aware transmission. The gap gets smaller with the increase of the UAVs height. Although the coverage probability is better for backhaul aware association, this scheme suffers from excessive handovers (i.e., ping pong effect) between UAV association and BS association. Such behavior highlights the validity of the proposed aware transmission scheme, especially for high UAVs altitudes.
\vspace{-0.8cm}
\subsection{Discussions and Design Insights}
	\vspace{-0.3cm}
In this section, we investigate the impacts of different system parameters to illustrate generic design guidelines and recommendations. Our results reveal that the height and number of assisting UAVs should be carefully adjusted to achieve optimal system performance. Furthermore, the characteristics of the terrestrial network must be taken into consideration while increasing the number of deployed UAVs. An interesting observation found in Fig.~\ref{fig:aware_ha_beta} is that an optimal height exists for the UAVs at which the coverage probability is maximized. This result is due to the applied probabilistic LOS/NLOS propagation model, implying that transmitting signals from higher altitude can benefit from increased possibility of LOS links. The optimal height is the point at which the UEs strike a balance between benefiting from the improved quality of received signals due to LOS and maintaining a good backhaul link between the UAVs and the terrestrial network. A trade-off can also be observed in Fig.~\ref{fig:coverage_aware_unaware100} regarding the number of deployed UAVs, which demonstrates that dense UAV deployment improves the UE performance if the UAVs are hovering in the low altitudes region and the backhaul links are guaranteed. For high heights, sparse UAV deployment is necessary to lower the interference.

Based on the interpretation of Fig.~\ref{fig:coverage_lambda}, we can notice that, for dense terrestrial BSs deployments, assisting the network with UAVs has a limited contribution on the coverage of UEs. The obtained results in Fig.~\ref{fig:misalignment_error} reveal the importance of considering the backhauling quality of the deployed UAVs as a key parameter in the design of the hybrid network. The advantages of assisting the terrestrial network with aerial BSs would not appear if these BSs are not supported with sufficient backhauling. In addition, adding UAVs with bad backhaul links would harm the UEs performance. Fig.~\ref{fig:coverage_fraction} offers the network operator the possibility of upgrading only a fraction of the existing BS network to provide backhauling for UAVs and optimized coverage for UEs; thereby reducing the needed infrastructure cost. Finally, supporting the UAVs with sufficient processing capabilities to be aware of the backhaul link quality would further improve the coverage.
\vspace{-0.6cm}
\section{Conclusion}\label{sec:conclusion}
\vspace{-0.2cm}
This paper uses stochastic geometry to assess the performance of mmWave backhauling for UAVs in a hybrid aerial-terrestrial cellular network considering key system parameters, e.g., UAVs height and number, and beamforming alignment error. We addressed both the backhaul unaware and the backhaul aware transmission scenarios. After characterizing the backhaul probability, the association and transmission probabilities and the distance distributions, we obtained an expression for the overall coverage probability and validated our results using Monte-Carlo simulations. Our results show that the quality of the backhaul link has a significant impact on the UE experience, and directly affects the deployment of UAVs for assisting terrestrial networks.
\vspace{-1.3cm}