
\subsection{SINR of each link}
Knowing the distance $r_k(t)$ between the two communicating nodes, the SINR of each link in LoS/NLoS can be defined as
\begin{equation}
    \chi_{k,n|r_{k}(t)}(t) = \frac{P_k G_{k,M} L^{-1}_{k,n}(r_k(t))|h_{k,n|r_k(t)}|^2}{I_{k}(t) + \sigma^2}, \tag{9}
\end{equation}
with $P_{SD} = P_{SR} = P_T$, $P_{RD} = P_R$, $h_{k,n|r_k(t)}$ being the small-scale fading of the considered link, and
\begin{align}
   & L_{SD,n}(r) = L_{Gn}(1+r), \nonumber\\
   & L_{SR,n}(r) = L_{RD,n}(r) = L_{An}\left(\sqrt{r^2+H_R^2}\right). \nonumber
\end{align}
The aggregate interferences of the three links are given by
\begin{equation}
I_k(t) = \sum_{\bm{\phi}_j \in \bm{\Phi}_{k,I}(t)} P_k G_{k,m} L_{k,n_j}^{-1}(\rVert \bm{\phi}_j \rVert)|h_{k,n_j|\rVert\bm{\phi}_j\rVert}|^2, \tag{10}
\end{equation}
with $n_j = \{L,N\}$ depending on whether the link between $\bm{\phi}_j$ and the node of interest is in LoS or in NLoS, and $h_{k,n_j,j|\rVert\bm{\phi}_j\rVert}$ being the small-scale fading of this link.\smallskip

The interferers' PPs are defined as 
\begin{equation}
\bm{\Phi}_{k,I}(t) = \bm{\Phi}_k \backslash \{\tilde{\bm{\phi}}_{k}(t)\},\tag{11}
\end{equation}
with $\bm{\Phi}_{SD} = \bm{\Phi}_{SR} = \bm{\Phi}_T$, and $\bm{\Phi}_{RD} = \bm{\Phi}_R$. Even if the selected TBS can vary at different time instants, with instantaneous handovers, \cite{Banagar2019} has shown that the distribution of $r_{SR}(t)$ is time-invariant.\smallskip

Knowing $r_{SD}$ and $r_{SR}(t)$, the densities of the interferers' PPs for the direct link and the first hop of the relay link are given by
\begin{alignat}{10}
&\lambda_{SD,I|r_{SD}}(r,t) &&=\lambda_T \: u(r-r_{SD}), \tag{12}\\
&\lambda_{SR,I|r_{SR}(t)}(r,t) &&= \lambda_T \: u(r-r_{SR}(t)). \tag{13}
\end{alignat}

\begin{figure*}[b]
\rule{\linewidth}{0.5pt}
\begin{multline}
  \mathcal{A}_{SRD} = \exp{-\pi\lambda_T H_R^2}\left(1-\exp{-\pi\lambda_R\: v^2t^2}\right) \\+ 2\pi\lambda_R\:\exp{-\pi\lambda_T(H_R^2+v^2t^2)}\int_{vt}^\infty r\:\exp{-\pi(\lambda_R+\lambda_T)r^2 + 2\pi\lambda_T vt\:r}\mathrm{d}r\label{eqn:probsel} \tag{18}
\end{multline}
\end{figure*}
\vspace*{-0.0cm}

\subsection{PDF of the initial distances between the selected nodes}
Since the TBSs' and UAVs' HPPPs are independent, $r_{SD}$ and $r_{RD}(t)$ are independent. It is also the case for $r_{SR}(t)$ and $r_{RD}(t)$. However, $r_{SR}(t)$ is dependent of $r_{SD}$ since there is a non zero probability for the same TBS to be selected by the typical UE and the selected RN, as illustrated by Lemma 1.
\medskip

\textit{\textbf{Lemma 1:} The probability for the same TBS to be selected by the typical UE and the selected RN, conditioned on $r_{SD}$ and $r_{RD,0}$ is given by
\begin{align}
    &\Pr{\bm{\tilde{\phi}}_{SD},\bm{\tilde{\phi}}_{SR}(t)|r_{SD},r_{RD,0}}{\bm{\tilde{\phi}}_{SD} = \bm{\tilde{\phi}}_{SR}(t)}\label{eqn:probOfdependency} \tag{14} \\ 
    &= \exp{-\pi\lambda_T\left(r_{SD}^2+r_{RD}^2(t)\right)}I_0\left(2\pi\lambda_T\:r_{SD}r_{RD}(t)\right). \nonumber
\end{align}}
\textit{Proof:} see Appendix A. \hfill{\small$\blacksquare$}\medskip

Therefore, in order to simplify the mathematical expressions, Approximation 1 is introduced.\medskip

\textit{\textbf{Approximation 1:} We assume that $r_{SR}(t)$ is independent with $r_{SD}$.}
\medskip

Since \eqref{eqn:probOfdependency} decreases with $\lambda_T$, Approximation 1 becomes more accurate when the TBS' density increases because the dependence between $\tilde{\bm{\phi}}_{SD}$ and $\tilde{\bm{\phi}}_{SR}(t)$ decreases.\medskip

\textit{\textbf{Lemma 2:} With the nearest neighbour association rule, the PDFs of the initial distances between the couple of nodes of the first hop, second hop and direct link are given by
\begin{equation}
   f_{r_{k,0}}(r_{k,0}) = 2\pi\lambda_k\:r_{k,0}\: \exp{-\pi\lambda_k r_{k,0}^2},\label{eqn:PDFrSD} \tag{15}
\end{equation}
with $\lambda_{SD} = \lambda_{SR} = \lambda_T,\ \lambda_{RD} = \lambda_R$.}\smallskip

\textit{Proof:} see Appendix B. \hfill{\small$\blacksquare$}\medskip

With the nearest neighbour association rule, we define the regions of association $\Omega_{SD}(t)$ and $\Omega_{SRD}(t)$ as
\begin{alignat}{10}
    &\Omega_{SD}(t) &&:= \left\{(r_{SD},r_{RD,0})\:|\: r_{SD}^2 \leq r_{RD}^2(t) + H_R^2\right\}, \tag{16}\\
    &\Omega_{SRD}(t) &&:= \mathbb{R}_+^2\:\backslash\:\Omega_{SD}(t). \tag{17}
\end{alignat}
Therefore, at time $t\geq 0$, for a given couple of initial distances $\bm{x}=(r_{SD},r_{RD,0})$, the typical UE selects the direct or relay link respectively if $\bm{x}\in\Omega_{SD}(t)$ or $\bm{x}\in\Omega_{SRD}(t)$. This is illustrated in Figure \ref{fig:associationRegionsNNAMob}. The integration of the PDFs on these regions gives the probability of association with the direct and relay link. The result is presented in Lemma 3.\medskip

\begin{figure*}[b]
\rule{\linewidth}{0.5pt}
\begin{alignat}{2}
&\mathcal{I}_{SD}^{(a)}(\beta) &&= \int_0^{H_R}\mathcal{P}_{SD|r_{SD}=r}(\beta)f_{r_{SD}}(r) \:\mathrm{d}r\label{eqn:I1}\tag{20}\\
&\mathcal{I}_{SD}^{(b)}(\beta,t) &&= \exp{-\pi\lambda_R \left(v^2t^2-H_R^2\right)} \int_{H_R}^\infty \mathcal{P}_{SD|r_{SD}=r}(\beta) f_{r_{SD}}(r)\: \exp{-\pi\lambda_R \left(r^2 + 2vt\sqrt{r^2-H_R^2}\right)}\mathrm{d}r\label{eqn:I2}\tag{21}\\
&\mathcal{I}_{SRD}^{(a)}(\beta,t) &&= \exp{-\pi\lambda_T H_R^2}\int_0^{\infty} \mathcal{P}_{SR|r_{SR}(t)=r}(\beta) f_{r_{SR}(t)}(r) \: \mathrm{d}r\int_{0}^{vt} \mathcal{P}_{RD|r_{RD,0}=r}(\beta,t) f_{r_{RD,0}}(r)\:\mathrm{d}r\label{eqn:I3}\tag{22}\\
\nonumber&\mathcal{I}_{SRD}^{(b)}(\beta,t) &&= \exp{-\pi\lambda_T \left(H_R^2+v^2t^2\right)}\int_0^{\infty} \mathcal{P}_{SR|r_{SR}(t)=r}(\beta) f_{r_{SR}(t)}(r) \: \mathrm{d}r\\
&&&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\cdot\int_{vt}^\infty \mathcal{P}_{RD|r_{RD,0}=r}(\beta,t)f_{r_{RD,0}}(r)\:\exp{-\pi\lambda_T (r^2-2vt\:r)}\mathrm{d}r\label{eqn:I4}\tag{23}
\end{alignat}
\rule{\linewidth}{0.5pt}
\begin{equation}
    \mathcal{P}_{k|n,r_{k}(t)}(\beta,t) = \frac{1}{2} + 
    \frac{1}{\pi}\int_0^\infty\frac{1}{\tau} 
    \mathrm{Im}\left[\mathcal{L}_{|h_{k,n|r_k(t)}|^2}\left(-j\tau P_kG_{k,M}L_{k,n}^{-1}(r_k(t))\right)
    \mathcal{L}_{I_{k}(t)|r_{k}(t)}(j\tau\beta) \: \exp{-j\tau\beta \sigma^2}\right]\mathrm{d}\tau\label{eqn:PSDRSD} \tag{25}
\end{equation}
\rule{\linewidth}{0.5pt}
\begin{equation}
    \mathcal{L}_{I_k(t)|r_k(t)}(s) = \exp{-2\pi \int_{0}^\infty \left(1 - \sum_{n=\{L,N\}} p_{k,n}(r) \mathcal{L}_{|{h}_{k,n|r}|^2}\left(s\:P_k G_{k,m}L_{k,n}^{-1}(r) \right)\right) \lambda_{k,I|r_k(t)}(r,t)\: r\mathrm{d}r} \label{eqn:LIT} \tag{27}
\end{equation}
\end{figure*}

\textit{\textbf{Lemma 3:} With the nearest neighbour association rule, the probability $\mathcal{A}_{SRD}$ that the typical UE selects the relay link at time $t\geq 0$ is given by \eqref{eqn:probsel}.}
\vspace*{0.1cm}

\textit{Proof:} see Appendix C. \hfill{\small$\blacksquare$}\medskip

\begin{figure}[H]
\begin{center}
\begin{adjustbox}{width = 0.8\linewidth}
\begin{tikzpicture}
    \node (RSD) at (1,4) {$r_{SD}$};
    \node (RRD) at (7,1) {$r_{RD,0}$};
    \node[anchor=north] (vt) at (2.5,1) {$vt$};
    \node[anchor=east] (HR)  at (1,2) {$H_R$};
    \draw[->,thick] (1,1) -- (RSD);
    \draw[->,thick] (1,1) -- (RRD);	
    \draw[thick,-] (HR) -- (2.5,2);
    \draw[-,black!20!white,dashed] (vt) -- (2.5,2);
    \draw[thick,-,dashed,black!50!white] (2.5,2) -- (2.5,4.2);
    \draw[thick,-,dashed,black!50!white] (2.5,2) -- (6.5,2);
    \node[fill=white] (SD)  at (5.25,2) {$\Omega_{SD}(t)$};
    \node[black!50!white] (ISDA) at (3,1.5) {$\mathcal{I}_{SD}^{(a)}$};
    \node[black!50!white] (ISDB) at (5.5,3) {$\mathcal{I}_{SD}^{(b)}$};
    \node[fill=white] (SRD) at (2.7,3.5) {$\Omega_{SRD}(t)$};
    \node[black!50!white] (ISRDA) at (1.8,2.8) {$\mathcal{I}_{SRD}^{(a)}$};
    \node[black!50!white] (ISRDB) at (3.8,4.2) {$\mathcal{I}_{SRD}^{(b)}$};
    \draw[thick,samples=100,domain=0:2.9,shift={(2.5,2)},
    postaction={decorate,decoration={text along path,text align={right},raise=1ex,text={{\scriptsize $r_{SD}^2=r_{RD}^2(t)+H_R^2$}{}}}}] plot ({\x}, {sqrt(\x^2+1)-1});
\end{tikzpicture}
\end{adjustbox}
\end{center}
\caption{\small Regions of association with the nearest neighbour association rule for the second mobility scheme at time $t\geq 0$. The typical UE selects the direct or the relay link if the couple $(r_{SD},r_{RD,0})$ is respectively in $\Omega_{SD}(t)$ or $\Omega_{SRD}(t)$.}
\label{fig:associationRegionsNNAMob}
\end{figure}

\subsection{Coverage Probability}
Using the definition of the SINRs of each link, Lemmas 2 and Approximation 1, Proposition 1 gives the mathematical expression of the CP of the typical UE.
\medskip 

\textit{\textbf{Proposition 1:} With the second mobility scheme, considering the nearest neighbour association rule and selection between direct and relay link, the approximate CP of the typical UE at time $t\geq 0$ for a given SINR threshold $\beta$ is given by 
\begin{equation}
\mathcal{P}(\beta,t) = \mathcal{I}_{SD}^{(a)}(\beta) + \mathcal{I}_{SD}^{(b)}(\beta,t) + \mathcal{I}_{SRD}^{(a)}(\beta,t) + \mathcal{I}_{SRD}^{(b)}(\beta,t), \tag{19}
\end{equation}
where $\mathcal{I}_{SD}^{(a)}$, $\mathcal{I}_{SD}^{(b)}$, $\mathcal{I}_{SRD}^{(a)}$ and $\mathcal{I}_{SRD}^{(b)}$ are respectively given by \eqref{eqn:I1}-\eqref{eqn:I4}. The subscripts $SD$ and $SRD$ correspond to whether the considered couples $(r_{SD},r_{RD,0})$ are included in $\Omega_{SD}(t)$ or $\Omega_{SRD}(t)$, as shown in Figure \ref{fig:associationRegionsNNAMob}. In these expressions, $\mathcal{P}_{k|r_{k}(t)}$ are respectively the CP of each link conditioned on the distance between the two nodes of interest at time $t$.}\smallskip

\textit{Proof:} see Appendix D. \hfill{\small$\blacksquare$}\medskip

Moreover, the CPs of every link taken independently are obtained by removing the conditioning of $\mathcal{P}_{k|r_{k}(t)}$. Using the law of total probability, we have
\begin{equation}
    \mathcal{P}_{k|r_{k}(t)}(\beta,t) = \sum_{n=\{L,N\}}p_{k,n}(r_k(t))\:\mathcal{P}_{k|n,r_k(t)}(\beta,t). \tag{24}
\end{equation}
These expressions are functions of the CPs of each hop if these links are in LoS or NLoS, and respectively conditioned on the distances between the nodes of interest. The Gil-Pelaez theorem enables to compute these quantities based on their characteristic functions \cite{GilPelaez}. It enables to reduce the mathematical complexity of the framework (no multiple numerical derivations of the Laplace Transforms (LT) of the aggregate interferences for each link) using numerical integrations. Based on this theorem, Proposition 2 gives the expressions of $\mathcal{P}_{k|n,r_{k}(t)}$.
\medskip

\textit{\textbf{Proposition 2:} The CPs of each link in LoS or NLoS, conditioned on the distance between the two communicating nodes are given by \eqref{eqn:PSDRSD}. In these expressions, $\mathcal{L}_{|h_{k,n|r_k(t)}|^2}$ are the LTs of the small-scale fading channel power gain of each link, while $\mathcal{L}_{I_{k}(t)|r_{k}(t)}$ are the LTs of the aggregate interferences respectively conditioned on $r_{SD}$, $r_{SR}(t)$ and $r_{RD}(t)$.}\smallskip

\textit{Proof:} see Appendix E. \hfill{\small$\blacksquare$}\medskip

\subsection{LTs of the small-scale fading channel power gains}
\textit{\textbf{Lemma 4:} Under Rician fading, the LTs of the small-scale fading channel power gains of each link are given by
\begin{equation}
\mathcal{L}_{|h_{k,n|r}|^2}(s) = \frac{K_{k,n}(r)+1}{K_{k,n}(r)+1+s}\:\exp{-\frac{K_{k,n}(r)\:s}{K_{k,n}(r)+1+s}}.
\label{eqn:LTRicianFading} \tag{26}
\end{equation} 
} \smallskip
\textit{Proof:} see Appendix F. \hfill{\small$\blacksquare$}\medskip

\subsection{LTs of the aggregate interferences}
\textit{\textbf{Proposition 3:} The LTs of the aggregate interferences of each link at time $t \geq 0$, respectively conditioned on $r_{SD}$, $r_{SR}(t)$ and $r_{RD}(t)$ are given by \eqref{eqn:LIT}.}\smallskip

\textit{Proof:} see Appendix G. \hfill{\small$\blacksquare$}\medskip
\subsection*{Contributions}
The contributions of this paper are summarised here below:
\begin{enumerate}
    \item We consider a relay network composed of TBSs and UAV RNs in an infinite area, with positions modelled using two independent HPPPs. The UAV RNs are either hovering or are allowed to move following two mobility schemes presented in \cite{Banagar2019}. The association policy is the following: the typical UE is either served via a direct link (from a TBS), or via a relay link (TBS - UAV RN - UE). The selection of the link is performed using a nearest neighbour association rule. 
    \item We derive an approximate expression for the Coverage Probability (CP) of the typical UE, and the CP of each link taken independently. We resort to the Gil-Pelaez theorem \cite{GilPelaez} for an increased tractability.
    \item Different path loss and small-scale fading models are used depending on the considered link, in order to make the scenario more realistic. Every link is either in LoS or NLoS, based on a probability model developed for Air-to-Ground (A2G) and Ground-to-Ground (G2G) cases. The parameters of the fading and path loss of each link depend on its LoS/NLoS and A2G/G2G nature. Regarding the small-scale fading, Rician fading is considered, and the introduction of distance-dependent K factors is proposed.
\end{enumerate}
\smallskip

This paper is organised as follows: the system model is detailed in Section \ref{section:system_model}. Then, the coverage probability is derived in Section \ref{section:coverage_probability}. Finally, numerical results are presented in Section \ref{section:numerical_results}.


\subsection{Network topology}
We consider a typical UE located at $(0,0)$ in $\mathbb{R}^2$. The locations of the UAV RNs and TBSs are modelled as two independent HPPPs $\bm{\Phi}_R$ and $\bm{\Phi}_T$ of densities $\lambda_{R}$ and $\lambda_{T}$. The UAVs are all flying at a fixed altitude $H_R$, while the TBSs and the typical UE are at ground level. The transmit powers of the UAVs and TBSs are respectively given by $P_R$ and $P_T$. The network is assumed to be highly congested, such that all RNs and TBSs interfere with each other while transmitting.

\subsection{Mobility Schemes}
The UAVs are allowed to move following two mobility schemes in $\mathbb{R}^2$ presented in \cite{Banagar2019}, both at a fixed altitude.\smallskip

In the first mobility scheme, all the UAVs move independently at a constant speed $v$ in random directions in $\mathbb{R}^2$. It has been shown in \cite{Banagar2019} that the RNs' HPPP remains an HPPP with the same density at every time $t \geq 0$ with this scheme. Therefore, assuming an instantaneous handover between different RNs, the performance at the typical UE is not affected by the movements of the UAV RNs.\smallskip

In the second mobility scheme, all the UAVs move independently at a constant speed $v$ in random directions in $\mathbb{R}^2$, except the serving UAV, which moves towards the typical UE at position $(0,0)$ and then hovers above it. Knowing the initial distance between the selected RN and the typical UE, the density of the RNs interferers' Point Process (PP) is given by \eqref{eqn:densityMob} in that scenario \cite{Banagar2019}. With this mobility scheme, the distance between the typical UE and the selected RN at time $t\geq 0$ is given by 
\begin{equation}
    r_{RD}(t) = (r_{RD,0} - vt)\:u(r_{RD,0} - vt), \tag{2}
\end{equation}
with $r_{RD,0}=r_{RD}(0)$, and $u$ being the Heaviside step function.

\subsection{Association Rule}
For sake of mathematical tractability, the nearest neighbour criterion is used as association rule to select the serving nodes: 
\begin{equation}
    \tilde{\bm{\phi}} = \argmin_{\bm{\phi}_i \in \bm{\Phi}} \left\{\left\rVert \bm{\phi}_i - \bm{x} \right\rVert \right\}. \tag{3}
\end{equation}
This rule is also applied to make the selection between the best direct link (TBS to UE) and the best relay link (TBS to UE through a RN). Since the selected links can vary at different time instants with mobile UAVs, we assume instantaneous handovers between different UAVs or between direct and relay links.\smallskip

In summary, the typical UE selects the direct link if the distance with its closest TBS $\tilde{\bm{\phi}}_{SD}$ is smaller than the distance with its closest RN $\tilde{\bm{\phi}}_{RD}(t)$. Otherwise the relay link is selected. In that case, the selected RN is served by its closest TBS $\tilde{\bm{\phi}}_{SR}(t)$. This scenario is illustrated in Figure \ref{fig:scenario}.

\begin{figure}[H]
    \centering
    \includegraphics[width=0.7\linewidth]{Images/scenario.PNG}
    \caption{\small Illustration of the scenario: the typical UE is at $(0,0)$ and selects either the direct or the relay link depending on the distances between the nodes. The RNs and TBSs are respectively represented by blue and red dots. The selected nodes are circled in green.}
    \label{fig:scenario}
\end{figure}
\vspace*{-0.3cm}

\subsection{Relay Transmission Protocol}
In this paper, the UAVs adopt the Decode-and-Forward (DF) protocol, executed over two time instants. At the first time instant, all the TBSs are transmitting to UEs or RNs. Thus the TBSs interfere with each other at that moment. At the second time instant, all the RNs re-transmit the signal to the UEs. Therefore, the RNs are interfering with each other at that instant. The Signal-to-Interference-plus-Noise-Ratio (SINR) of the relay link is given by \cite{Laneman2004} 
\begin{equation}
    \chi_{SRD}(t) = \min\left(\chi_{SR}(t),\chi_{RD}(t)\right),
    \label{eq:SINRDF} \tag{4}
\end{equation}
where $\chi_{SR}(t)$ is the SINR between the TBS and the RN, and $\chi_{RD}(t)$ is the SINR between the RN and the typical UE. 

\subsection{Beamforming model}
The TBSs and UAVs are equipped with antenna arrays to generate directional beams and perform transmit and receive beamforming. This leads to an SINR increase for the same transmitted power, enabling to compensate for a high path loss. As proposed in \cite{Belbase2018,Biswas2016,Ma2019}, a simple sectored antenna model is introduced:
\begin{equation}
    G(\psi,\theta) = \left\{
    \begin{matrix}
    G_{M} &\mathrm{if} |\psi| \leq \frac{\psi_{bw}}{2}\ \mathrm{and}\ |\theta| \leq \frac{\theta_{bw}}{2}\\
    G_{m} &\mathrm{otherwise}
    \end{matrix}
    \right., \tag{5}
\end{equation}
where $G_{M}$ is the maximum beamforming gain obtained when the elevation angle $\theta$ and the azimuth angle $\psi$ are respectively within the half elevation and azimuth power beamwidth $\theta_{bw}$ and $\psi_{bw}$, and $G_{m}$ is the minimum beamforming gain obtained in the other cases. In practice, it will be assumed that the beam of the target receiver is perfectly aligned with the beam of the transmitter and misaligned with the beams of the interferers, leading thereby to the array gains listed in Table \ref{tab:beamforming_gains}. The maximum and minimum beamforming gains of the TBSs and RNs are respectively denoted by $G_{TM}$, $G_{Tm}$, $G_{RM}$ and $G_{Rm}$.
\begin{table}[H]
\begin{center}
\begin{tabular}{|l||l|l|}
    \hline
         & \multicolumn{1}{c|}{Target} & \multicolumn{1}{c|}{Interference} \\\hline\hline
            {TBS-UE} & $G_{SD,M}=G_{TM}$ & $G_{SD,m}=G_{Tm}$ \\\hline
        {TBS-RN} & $G_{SR,M}=G_{TM}G_{RM}$ & $G_{SR,m}=G_{Tm}G_{Rm}$ \\\hline
        {RN-UE} & $G_{RD,M} = G_{RM}$ & $G_{RD,m} = G_{Rm}$ \\\hline
    \end{tabular}
    \caption{\small Beamforming gains for the different links.}
    \label{tab:beamforming_gains}
\end{center}
\end{table}
\vspace*{-0.6cm}
\subsection{Propagation model}
The channel model takes into account path loss and small-scale fading. The parameters of the different models are different for A2G and G2G links, in LoS or NLoS. We also introduce a model for the LoS probability.\smallskip

\paragraph{Path loss}
We consider two linear path loss models for A2G and G2G links \cite{Khawaja2019}:
\begin{equation}
    L_{ln}(r)\ [dB] = \beta_{ln} + 10\:\alpha_{ln} \: \log_{10}\left(r\right), \tag{6}
\end{equation}
where $r$ is the distance between the two communicating nodes, $l = \{A,G\}$ for A2G or G2G links, and $n = \{L,N\}$ for links in LoS or NLoS. The parameters $\beta_{ln}$ are the floating intercepts while the parameters $\alpha_{ln}$ are the path loss exponents. In order to avoid convergence issues with the numerical integration of the expressions for G2G links, as in \cite{Khoshkholgh2019}, the distance $r$ is replaced by $1+r$ in the path loss model.\smallskip

\paragraph{Small-scale fading}
We assume that the propagation is affected by Rician fading. Knowing the ground distance $r$ (in two dimensions) between the two communicating nodes, it is modelled by a random variable $h_{ln|r}$, with normalised power and a distance-dependent K factor $K_{ln}(r)$, different for A2G/G2G links in LoS/NLoS. We suppose that the small-scale fading of different links are independent, and that the random variables at different time instants for the same link are independent.\smallskip

\paragraph{LoS probability (G2G links)} 
In this paper, the 3rd Generation Partnership Project (3GPP) $d_1/d_2$ model is selected, with a LoS probability given by \cite{Sun2015}
\begin{equation}
p_{GL}(r) = \min\left(\frac{d_1}{r},1\right) \left(1 - \exp{-\frac{r}{d_2}}\right) + \exp{-\frac{r}{d_2}}, \tag{7}
\end{equation}
where $r$ is the ground distance. The parameters $d_1$ and $d_2$ are respectively called the near-field and far-field critical distances.\smallskip

\paragraph{LoS probability (A2G links)} A commonly used LoS probability function for A2G links is given by \cite{Al-Hourani2014}
\begin{equation}
    p_{AL}(r) = \frac{1}{1 + a\  \exp{-b\left(\frac{180}{\pi}\arctan\left(\frac{H_R}{r}\right)-a\right)}},
    \label{eqn:LOSprobA2G} \tag{8}
\end{equation}
where $r$ is the ground distance, and $a$ and $b$ are two fitting parameters, directly linked to three other parameters with a physical meaning \cite{Al-Hourani2014}: the ratio of built-up land area to the total land area, the mean number of buildings per unit area, and a scale parameter that describes the buildings' heights distribution according to a Rayleigh Probability Density Function (PDF).\smallskip

\paragraph{Noise} We suppose that every link is affected by a white Gaussian noise of constant power $\sigma^2$.
\section{Introduction}
\input{Introduction.tex}

\section{System Model}
\input{System_Model.tex}

\section{Coverage Probability}
\label{section:coverage_probability}
\input{Coverage_Probability.tex}

\section{Numerical Results}
\input{Numerical_Results.tex}

\section{Conclusion}
\label{section:conclusion}
\input{Conclusion.tex}

\bibliographystyle{ieeetr}
