Performance of Dense Wireless Networks in 5G and beyond Using Stochastic Geometry

Abstract

Device density in cellular networks is expected to increase considerably in the near future. Accordingly, the access point (AP) will be equipped with massive multiple-input multiple-output (mMIMO) antennas, using collimated millimeter-wave (mmW) and sub-THz communications, and increasing the bandwidth to accommodate the growing data rate demands. In this scenario, interference plays a critical role and, if not characterized and mitigated properly, might limit the performances of the network. In this context, this paper derives the statistical properties of the aggregated interference power for a cellular network equipping a mMIMO cylindrical array. The proposed statistical model considers the link blockage and other network parameters such as antenna configuration and device density. The findings show that the characteristic function (CF) of the aggregated interference power can be regarded as a weighted mixture of two alpha-stable distributions. Furthermore, by analyzing the service probability, it is found that there is an optimal configuration of the array depending on the AP height and device density. The proposed statistical model can be part of the design of dense networks providing valuable insights for optimal network deployment and resource management and scheduling.

1. Introduction

The fifth-generation (5G) cellular network has been recently deployed with unprecedented communication performance, i.e., a 10–100 times higher data rate, 1ms latency, and much higher area throughput [1,2]. The upcoming sixth-generation (6G) cellular network promises to further improve current performance by at least one order of magnitude [3]. To meet such requirements, it is necessary to operate on multiple frontiers, e.g., increase the bandwidth, cell density, and transceiver efficiency.

Current cellular networks operate at the sub-6 GHz band, which is heavily congested [3]. Recently, millimeter waves (mmW) and sub-THz frequencies (30–300 GHz) have enjoyed substantial interest due to the large unexploited spectrum [4,5,6]. However, propagation at these frequencies experiences higher path and penetration loss, making the links prone to blockage [7]. A solution to these challenges is to use a beam-type communication based on massive multiple-input multiple-output (mMIMO) systems and increase the cell density [8]. However, as the device density increases, interference emerges as one of the main challenges to be characterized and mitigated. The characterization of the aggregate interference power and coverage analysis in millimeter-wave (mmW) networks has been investigated in recent years but only in some specialized settings, under line-of-sight (LOS) and non-line-of-sight (NLOS) propagation. Among different models [9], $\alpha$-stable distribution has been widely accepted for modeling the aggregated interference power in a wireless network generated by PPP-distributed UEs [10,11]. The heavy-tailed behavior has been validated experimentally in the context of the Internet of Things (IoT) [12,13]. Generally, $\alpha$-stable distribution does not have a defined PDF model, except for specific cases. In [14], the exact PDF for some of the values of $\alpha$ has been derived in the context of wireless communications. The $\alpha$-stable distribution when the interferers transmit in multiple sub-bands is studied in [15]. In all of these studies, the considered scenario is a PPP on a 2D plane. The case of PPP in a 3D volume has also been modeled as an $\alpha$-stable distribution in [10]. However, in cellular networks, the distribution of the UEs is not a homogeneous PPP in a 3D volume, but rather a PPP distribution of UEs on a 2D plane, where the AP is placed at an arbitrary height h. When the height of the AP is considered, the heavy tail of the distribution becomes much shorter, similar to the interference power in cognitive radio networks, where the secondary UEs are not allowed to transmit when they are within the interference region [16,17].

The impact analysis of the user’s height was studied in Ref. [18]. The array’s height h plays an important role that needs to be analyzed, namely, in view of the 5G mmW and 6G use cases. Therefore, it is necessary to extend the Stochastic Geometry (SG) framework to the 3D framework of antennas and user equipments (UEs) (see e.g., [19]). SG provides a preferred framework in network modeling to perform coverage and rate performance analysis [17,20,21,22,23]. The impact of the antennas’ height in a 3D SG for ultra-dense networks proves an upper limit on network performance that depends on the path-loss model parameters [24]. Although the impact of the antenna and user equipment (UE) height difference has been studied [25], the existence of an optimum array height has not been deduced in dense networks. Here, we derive the analytical model of aggregated interference and show that the optimum array height depends on the path-loss model but also on the user’s density, array type, and size. Similarly, the impact of the height on low-altitude aerial platforms [26,27] and in unmanned aerial vehicles [28,29] proves that there are optimum altitudes maximizing the coverage probability according to some specific scenarios.

In the literature, it is common to assume that the fading follows Rayleigh distribution (or at least the interference link has a Rayleigh distribution) or Nakagami-m distribution. The coverage probability has a tractable form as a Laplace function of the aggregate interference power (see, e.g., [17,30,31,32]). However, for general fading, using the Laplace function of the aggregate interference is not possible anymore, and thus, the coverage and rate analysis cannot be expressed in a tractable way. To overcome this problem, in ref. [20], multiple techniques are introduced, including some methods to find the probability density function (PDF) of the aggregate interference calculated from the corresponding characteristic function (CF).

In previous works, interference distribution for single antennas and the uniform linear array (ULA) (e.g., in [30]) configuration was investigated most. The contribution of this paper is the usage of the 3D SG framework in the uplink for homogeneous Poisson point process (PPP) with a density $\lambda$ on a pseudo-3D geometry (sometimes referred to as 2.5D geometry [33,34]) where UEs lie in a plane (say ground) and the access point (AP) with $N_c\times N_v$ uniform cylindrical array (UcylA) (i.e., a set of $N_v$ half-wavelength rings of $N_c$ uniform circular array (UCA) antennas/each) has the height $h\ge 0$. Figure 1 illustrates the setup for the computation of the properties of the aggregate interference from the ensemble of $U{E}_1,U{E}_2,\dots U{E}_i,\dots$ when the array at the AP is pointing toward the user of interest $U{E}_0$. We show that by being able to characterize the aggregated interference, one can optimally adjust the configurations of the AP to increase the number of users served and, accordingly, the coverage probability.

We analytically prove in Section 4 that the aggregated interference power onto an array of antennas located at an arbitrary height $\left(h\right)$ has a CF is decomposable into a mixture of two stable distributions (a skewed alpha-stable and a Gaussian distribution). The closed-form CF for $h>0$ is another novelty to compute the PDF or cumulative density function (CDF) of the aggregated interference power, without cumbersome numerical integration. The scenario considered herein for users and interferers accounts for either LOS or NLOS propagation, and also the possibility of link blockage as is typical in mmW and sub-THz systems. The analytical model enables (Section 3.4) the evaluation of the trade-off in array height selection and the impact of the propagation scenarios. The UcylA with $h>0$ generalizes the previous works on ULA configuration, which was investigated mostly using coverage analysis (see, e.g., [30,35]). In array engineering, UcylA can be designed using multifaceted array in a far more practical way, and the results in [36] support the conclusion that any results for curved arrays apply to faceted ones.

Based on the analytical model, the main results can be summarized as follows: (i) in most scenarios with small interferers’ density, it is more beneficial to adopt a UCA rather than UcylA, while for large interferer density a UcylA would be preferred; (ii) there is an optimum AP height that depends on propagation and interferer density $\lambda$; (iii) at AP height $h=0$, the aggregated interference power is alpha-stable distributed, and for $h\to \infty$, the limit becomes Gaussian, but for any arbitrary height h, it is decomposable into two stable distributions; (iv) blockage probability impacts the service probability for small $\lambda$, but less if counting the average number of users served within a region; (v) a connection with multiple paths (LOS and NLOS) is more beneficial for small $\lambda$, as for large $\lambda$, the interference is too large.

This paper is organized as follows. We present the system model, including the signal and array gain models, in Section 2. The CF of the aggregate interference for UcylA in LOS links is described in Section 3.

Section 4 contains the statistical characterization of the interference power. In Section 5, the CF is extended considering NLOS paths, noise power, and blockage, and we conclude the paper in Section 6.

2. System Model

The scenario shown in Figure 1, where the AP has a cylindrical array with $N=N_c\times N_v$ antennas in total. The UEs are uniformly distributed following a homogeneous PPP with density $\lambda$ that denotes the mean number of active UEs per square meter. The spatial channel of the mmW and 6G sub-THz systems are purely directional (see, e.g., [37]), and the LOS (or NLOS) link is affected by the path-loss modeled in terms of UE-AP distance $d_{LOS}$ (or $d_{NLOS}$), and faded amplitude ${\beta }_{LOS}$ (or ${\beta }_{NLOS}$). The propagation attenuation model for LOS and NLOS is $\beta /d^b$ with amplitude path-loss $b\ge 1$. The array of antennas is uniformly cylindrical, the isotropic radiating antennas are arranged into a set of $N_v$ UCAs with $N_c$ antennas each, and antenna spacing is half the wavelength; that is, the two arrangements of antennas are such that the corresponding beamforming of UCAs reduce the interference angularly, and the vertical arrangement of the rings (acting as vertical ULAs), tilts the beam to improve the capability to reduce the near interferers when pointing toward far-away UEs. Note that here the UEs are considered on the ground, which means that the AP height is the height difference of UEs and AP (pseudo 3D or 2.5D geometry).

2.1. Array Gain Model

Each AP equipped with the array of antennas is positioned at height h from the ground at $(0,0)$ planar coordinates as in Figure 1. The array gain for the UcylA in far-field $G(\varphi ,\theta )$ depends on the elevation angle ($\theta$) and azimuth ($\varphi$), which in turn depends on the number of antennas partitioning between $N_v$ and $N_c$. The beamforming for the cylindrical array is conveniently decomposed into the design of two compound arrays, and thus the array gain $G(\varphi ,\theta )=G_c\left(\varphi \right)G_v\left(\theta \right)$ is separable into the UCA gain $G_c\left(\varphi \right)$ and vertical ULA gain $G_v\left(\theta \right)$ [38,39]. The approximation holds true in UcylA when using separable weightings [40,41]. The beamforming used here is the conventional method, which is optimum for uniform interference, and the array gains for half-wavelength inter-element spacing either for UCA and ULA are [41,42]:

$$\begin{array}{ccc}\hfill G_c\left(\varphi \right)& =& J_0\left(\frac{N_c}{2}\sqrt{{(cos\varphi -cos{\varphi }_o)}^2+{(sin\varphi -sin{\varphi }_o)}^2}\right),\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill G_v\left(\theta \right)& =& \frac{sin[\pi (sin\theta -sin{\theta }_o)N_v/2]}{N_{v}sin[\pi (sin\theta -sin{\theta }_o)/2]},\hfill \end{array}$$

(2)

where (${\varphi }_o,{\theta }_o$) denotes the pointing azimuth and elevation pair to the intended UE, and for UCA the gain approximation (1) is obtained using the zero-order Bessel function of the first kind $J_0(.)$, which can be shown to be accurate for $N_c\ge 16$. Note that array gains are normalized for convenience to have $G({\varphi }_o,{\theta }_o)=1$. The beam widths along the two angles, $\Delta \varphi$ and $\Delta \theta$, are inversely proportional to $N_c$ and $N_v$, respectively. The elevation beam width is further distorted by the effective array aperture, which makes the beam width scale with the cosine of the tilt angle (stretching effect): $\Delta \theta /cos(\theta -{\theta }_o)$.

2.2. Signal Model and Service Probability

Let x be the transmitted signal; the signal received by the AP with beamforming pointing toward the UE of interest with angles ${\varphi }_o=0$ and radial distance $R_o$ is

$$y=\frac{{\beta }_o}{{\left(R_o^2+h^2\right)}^{\frac{b}{2}}}G({\varphi }_0,{\theta }_0)x+\iota +w,$$

(3)

where $w\sim CN(0,{\sigma }_w^2/N)$ is the additive Gaussian noise with power ${\sigma }_w^2/N$ after the array gain, and

$$\iota =\sum _{i=1}^{\infty }\frac{{\beta }_i}{{\left(R_i^2+h^2\right)}^{\frac{b}{2}}}G({\varphi }_i,{\theta }_i)x_i.$$

(4)

is the aggregated interference originated from PPP-distributed interfering UEs with density $\lambda$; all signals generated by all UEs are $x_i\sim CN(0,1)$. We assume that the aggregate interference power $I={\left|\iota \right|}^2$ is typically $\mathbb{E}\left[I\right]\gg {\sigma }_w^2/N$ as in the macro-cell network [22]. This assumption is relaxed in Section 5, since in high frequencies, the noise is not negligible in small cells, especially in the presence of blockage [43].

The outage analysis depends on the PDF of the aggregated interference (4) for the PPP distribution of UEs, each having polar coordinates $(R_i,{\theta }_i)$. The PDF of I without any array-gain for interference mitigation (here is a specific case with $G^2({\varphi }_i,{\theta }_i)=1$ for any i) and $h=0$ has been extensively investigated in the literature (see, e.g., [11,44,45]). The scope here is to evaluate the PDF of aggregated interference I for UcylA and $h\ge 0$. The $N_c\times N_v$ UcylA is the most general case as UCA is when $N_v=1$ and the single antenna is when $N_c=N_v=1$, as shown in

Table 1. Array configurations and reference to the characteristic function (CF) of aggregated interference $\mathsf{\Psi }\left(\omega \right)=\mathbb{E}\left[e^{j\omega I}\right]$.

	UcylA	UCA $\mathit{h}>0$	UCA $\mathit{h}=0$	Isotropic $\mathit{h}>0$	Isotropic $\mathit{h}=0$
$N_c$	>1	>1	>1	1	1
$N_v$	>1	1	1	1	1
h	≥0	≥0	0	≥0	0
	Section 3.1 CF: (15)	Section 3.2 CF: (19)	Section 3.2 CF: (23)	Section 3.2 CF: (25)	Ref. [11]

, with the corresponding references or sections for the analytic form of CF. Note that the density $\lambda$ refers to the number of active UEs per square meter coexisting on the same time-frequency and, depending on the specific radio resource allocation strategies, is likely to be meaningfully lower than the effective crowd density [46]. The signal-to-interference ratio (SIR) is defined as

$$\mathrm{SIR}=\frac{|{\beta }_0{|}^2{G}^2({\varphi }_0,{\theta }_0)}{I{\left(R_o^2+h^2\right)}^b}.$$

(5)

The average probability of successful connection experienced by the UE of interest in $({\varphi }_o,R_o)$ depends on a certain threshold T, on fading fluctuation $|{\beta }_o{|}^2$, and on the overall interference I. The service probability from the distribution of I is

$$P_s\left(R_o\right|R_{max})=F_I\left(\frac{|{\beta }_0{|}^2{G}^2({\varphi }_0,{\theta }_0)}{T{\left(R_o^2+h^2\right)}^b}\right),$$

(6)

where $F_I\left(x\right)=Pr(I<x)$ is the CDF. The service probability (6) is for the interference power I, which accounts for the randomness of the position of the interferers according to the PPP model within a certain radius $R_{max}$ and fading. The fluctuations of the interferers $|{\beta }_i|$ are embodied in the CF derivations (Section 3). The reference user in $(R_o,{\varphi }_o)$ is considered as deterministic for the computations of the (conditional) service probabilities, but whenever necessary for the unconditional probability, it can be assumed as PPP distributed in the same way as for the other interferers. Notice that the service probability $P_s\left(R_o\right|R_{max})$ depends on the CDF of the aggregated interference I to be evaluated for the PPP distribution of active interferers, the array-type and its height h, as evaluated in the next section based on CF analysis. For the unconditional service probability accounting for the fading of the user of interest, one should evaluate the expectation

$${\tilde{P}}_s\left(R_o\right|R_{max})={\mathbb{E}}_{|{\beta }_0{|}^2}\left[P_s\left(R_o\right)\right].$$

(7)

Differently from the interference analysis that aggregates multiple (and many) interfering contributions into the CF of I, (7) depends on the specific PDF of ${\beta }_0$ thus making the service (or outage) analysis distribution-dependent for the UE of interest; see e.g., [30,47]. Therefore, the fast-fading contribution ${\beta }_0$ is neutrally modeled here as ${\beta }_0=1\times e^{jϵ}$ with $ϵ\sim \mathcal{U}(0,2\pi )$, and further derivations for arbitrary fading ${\beta }_0$ are out of the scope here.

Remark 1.

The service probability (6) is defined by neglecting the noise power. It can be shown that, with the configurations used in this paper for transmitter power $P_{tx}>$ 20 dBm, the noise effect does not affect the results, while the standard UE power in the 3GPP standard (refer to ETSI TS 138 101-2 V15.5.0 2019-05) is higher than $P_{tx}=$ 20 dBm. Noise is more relevant when considering the propagation blockage (Section 5.4) The results regarding the service probability prior to Section 5.3 are all based on interference, as defined in (6).

3. CF of the Aggregated Interference in LOS

In this section, the CF of the aggregate interference is derived in the presence of LOS links first for the general case of $N_c\times N_v$ UcylA, and is then simplified for the UCA and single antenna. The focus of this paper is to analyze the behavior of the interference and coverage with respect to the height and arrangements of the antenna array. The mmW channel modeling is first LOS-only, and is then enriched with mixed LOS and NLOS links in Section 5.

3.1. Uniform Cylindrical Arrays

$N_c\times N_v$ UcylA is composed of $N_v$ uniformly spaced rings consisting of $N_c$ antennas each arranged in a cylinder shape (Figure 1). The interference power I originated from a coverage radius $R_{max}\to \infty$ represents the largest possible interference for a density $\lambda$ and thus it is the upper bound of the interference I when $R_{max}<\infty$, the service probability (6) depends on $R_{max}$, and it is lower bounded: $P_s\left(R_o\right|R_{max})\ge P_s\left(R_o\right|R_{max}\to \infty )$. For the computation of the CF of the aggregate interference, one should consider the entire gain pattern of the UcylA (1) and (2). Let $R_{max}\to \infty$, the aggregate interference power is

$$I=\sum _{i=1}^{\infty }\frac{|{\beta }_i{|}^2}{{\left(R_i^2+h^2\right)}^b}{\tilde{G}}^2(R_i,{\varphi }_i),$$

(8)

where the served UE is in $({\theta }_o,{\varphi }_o)=(0,0)$ for analytical notation convenience, and the beamforming gains are reformulated in terms of azimuth (${\varphi }_i$) and elevation (${\theta }_i=arctan(h/R_i)$) angles

$$\tilde{G}(R_i,{\varphi }_i)=G_c^2\left({\varphi }_i\right)G_v^2(arctan(h/R_i))$$

(9)

The fluctuation power $|{\beta }_i{|}^2$ is independent of interfering users and identically distributed (iid). The randomly distributed interfering UEs in ${\varphi }_i\in [0,2\pi )$ can be partitioned into a set of K disjoint angular sectors ${\mathsf{\Phi }}_1,{\mathsf{\Phi }}_2,\dots$ such that ${\cup }_k{\mathsf{\Phi }}_k\equiv [0,2\pi ),$ where K is large enough so that the array gain $G_c\left({\mathsf{\Phi }}_k\right)$ in each sector can be considered as constant. The aggregate interference power (8) is

$$\begin{array}{cc}\hfill I& =\sum _{k=1}^K\sum _{{\varphi }_i\in {\mathsf{\Phi }}_k}\frac{|{\beta }_i{|}^2}{{\left(R_i^2+h^2\right)}^b}{\tilde{G}}^2(R_i,{\varphi }_i)\simeq \sum _k{I}_k,\hfill \end{array}$$

(10)

where

$$I_k=\sum _{i\in {\mathsf{\Phi }}_k}\frac{|{\beta }_i{|}^2}{{\left(R_i^2+h^2\right)}^b}{\tilde{G}}^2(R_i,{\mathsf{\Phi }}_k).$$

(11)

The asymptotic equality in (10) is due to the fact that every user i with azimuth ${\varphi }_i$ that is located within the k-th angular sector is assumed to have the array gain ${\tilde{G}}^2(R_i,{\mathsf{\Phi }}_k)$, as defined in (11). This approximation asymptotically converges to the true aggregated interference power for $k\to \infty$, while for k that is large enough (around twice the number of horizontal antennas), it is quite accurate. The CF for the interference within the kth angular sector follows on from Campbell’s theorem as in [11]

$${\mathsf{\Psi }}_{I_k}\left(\omega \right)=exp\left(-\frac{2\pi }{K}\lambda {\int }_0^{\infty }\left[1-{\mathsf{\Psi }}_{{\left|\beta \right|}^2}\left(\omega \frac{{\tilde{G}}^2(r,{\overline{\varphi }}_k)}{{\left(r^2+h^2\right)}^b}\right)\right]rdr\right),$$

(12)

where ${\mathsf{\Psi }}_{{\left|\beta \right|}^2}\left(\omega \right)=\mathbb{E}[e^{{j\omega \left|\beta \right|}^2}]$ is the CF of ${\left|\beta \right|}^2$. In order to solve the integral (12), it is more convenient to define $\omega {\left(r^2+h^2\right)}^{-b}=t$. By replacing the variable r with t and changing the boundaries of the integral accordingly, i.e., $t\to 0$ for $r\to \infty$ and $t\to \left|\omega \right|/h^{2b}$ for $r\to 0$, after some mathematical manipulations, the (12) becomes (for $\alpha =1/b$)

$${\mathsf{\Psi }}_{I_k}\left(\omega \right)=exp\left(-\frac{\pi }{K}\lambda \alpha {\left|\omega \right|}^{\alpha }{\int }_0^{\left|\omega \right|/h^{2b}}{\displaystyle \frac{1-{\mathbb{E}}_{{\left|\beta \right|}^2}[e^{{jt\left|\beta \right|}^2{\tilde{G}}^2(f\left(t\right),{\mathsf{\Phi }}_k)\mathrm{sign}\left(\omega \right)}]}{t^{\alpha +1}}}dt\right).$$

(13)

where $f\left(t\right)={\left({\left|\omega \right|}^{1/b}{t}^{-1/b}-h^2\right)}^{1/2}$ follows on from the conversion from variable r to t. (13) is based on radius $r\to \infty$; if the radius is limited, the lower bound of the integral of (13) would be $t\to \left|\omega \right|/{(h^2+R_{max}^2)}^b$ for $r\to R_{max}$. Since the array gain is a function of t, this expression can only be solved numerically. One way to make the CF tractable is by uniformly decomposing the elevation angle ${\theta }_v\in (0,\pi /2]$ into M angular sectors of ${\Delta }_v=\pi /2M$ width; each sector is centered in ${\tilde{\theta }}_m=\frac{\pi }{2}(\frac{2(M-m)+1}{2M})$; and the width ${\Delta }_v$ is small enough to let the array gain in every angular sector (2) be constant $G_v^2({\tilde{\theta }}_m)$. The array gain is constant on every annulus (ring)-shaped areas with unequal widths (non-uniform ring division for uniform elevations ${\tilde{\theta }}_m$). These rings are centered in the radial distance of the intersection of the mth bisector ${\tilde{\rho }}_m=h/tan({\theta }_m-\frac{\pi }{4M})$, so that the array gain in the kth wedge and mth ring is ${\tilde{G}}^2({\tilde{\rho }}_m,{\mathsf{\Phi }}_k)$. Within the kth wedge and mth ring, the array gain is constant and the interference originated $I_{k,m}$ has the CF

$${\mathsf{\Psi }}_{I_{k,m}}\left(\omega \right)=exp\left(-\frac{\pi }{K}\alpha \lambda {\left|\omega \right|}^{\alpha }{\int }_{{\tau }_{m+1}}^{{\tau }_m}\left[\frac{1-{\mathbb{E}}_{{\left|\beta \right|}^2}[e^{{jt\left|\beta \right|}^2{\tilde{G}}^2({\overline{\rho }}_m,{\mathsf{\Phi }}_k)\mathrm{sign}\left(\omega \right)}]}{t^{\alpha +1}}\right]dt\right),$$

(14)

where ${\tau }_m=\left|\omega \right|/{(h^2+{\rho }_m^2)}^b$. The CF of the aggregate interference statistically independent on all K wedges and M rings is ${\mathsf{\Psi }}_I\left(\omega \right)={\displaystyle \prod _{m=0}^{M-1}}{\displaystyle \prod _{k=0}^{K-1}}{\mathsf{\Psi }}_{I_{k,m}}\left(\omega \right)$. Given the solution of the integral ${\int }_0^{\left|\omega \right|/h^{2b}}\left[\frac{1-e^{j\mu t}}{t^{\alpha +1}}\right]dt=\underset{\epsilon \to 0}{lim}{\left.{(-j\mu )}^{\alpha }\mathsf{\Gamma }(-\alpha ,-j\mu t)-\frac{1}{\alpha t^{\alpha }}\right]}_{t=\epsilon }^{t=\left|\omega \right|/h^{2b}}$, for any constant and real $\mu$, and $0<\alpha <1$, the CF of the aggregate interference can be shown to reduce to

$${\mathsf{\Psi }}_I\left(\omega \right)=exp(-\lambda 2\pi R_{max}^2)exp\left(-\frac{\pi \lambda }{C_{\alpha }}{\left|\omega \right|}^{\alpha }(1-j\mathrm{sign}\left(\omega \right)tan\frac{\pi \alpha }{2})P_G(\omega ,\alpha ,h)\right),$$

(15)

with $C_{\alpha }^{-1}=\mathsf{\Gamma }(1-\alpha )cos(\pi \alpha /2)$ and

$$\begin{array}{cc}\hfill P_G(\omega ,\alpha ,h)=& \frac{{\overline{\beta }}^{2\alpha }}{K}\sum _{k=0}^{K-1}\sum _{m=0}^{M-1}{\tilde{G}}^{2\alpha }({\overline{\rho }}_m,{\mathsf{\Phi }}_k)\hfill \\ & \left[P(-\alpha ,-j\omega \xi ({\overline{\rho }}_m,{\mathsf{\Phi }}_k,h))-P(-\alpha ,-j\omega \xi ({\overline{\rho }}_{m+1},{\mathsf{\Phi }}_k,h))\right].\hfill \end{array}$$

(16)

We used a compact notation for different moments of fading ${\overline{\beta }}^c={\mathbb{E}[|\beta |}^c]$, and $P(x,z)={\int }_0^z{t}^{x-1}{e}^{-t}dt/\mathsf{\Gamma }\left(x\right)$ in (16) is the normalized incomplete Gamma function ratio and

$$\xi (r,\varphi ,h)=\frac{{\tilde{G}}^2(r,\varphi ){\overline{\beta }}^2}{{\left(h^2+r^2\right)}^b}.$$

(17)

The relationship (16) for $M\to \infty$ and $R_i\le R_{max}<\infty$, reduces after some calculus, to

$$P_G(\omega ,\alpha ,h)=\frac{2b}{{(-j\omega )}^{\alpha }\mathsf{\Gamma }(-\alpha )}{\int }_0^{R_{max}}\left({\int }_0^{2\pi }{e}^{j\omega \xi (r,\varphi ,h)}\frac{d\varphi }{2\pi }\right)rdr.$$

(18)

This relation (18) completes the CF of aggregated interference (15). The PDF and the CDF of the interference I are obtained by numeric inversion. For $h\to 0$, the UcylA degenerates into UCA on the ground, and the distribution of I is alpha-stable for $R_{max}\to \infty$ (Section 3.2). As can be observed, the general amplitude fluctuations ${\beta }_i$ of the interferers are embodied in the derivation of the CF (whose values depend on the specific fading model chosen); however, for the numerical calculations throughout the paper, the fading power ${\overline{\beta }}^2$ is neglected because it is averaged out in the SG over the summation of the interference for all UEs. This assumption holds as long as the average fading is $\mathbb{E}\left[\right|\beta \left|\right]=1$ [48] (the fading power for the user of interest $|{\beta }_0{|}^2$ does not average out, and the results depend on the choice of fading model).

Remark 2.

Although the CF of aggregated interference (15) for UcylA depends on (18), its numerical computation has some trade-offs. The granularity of numerical integration basically depends on the main beamwidth: the azimuth $K=N_c$ is a safe choice with good accuracy, while for the elevation angle, $M=N_v/$2 to $M=N_v/$4 is acceptable ($M=N_v/$2 is a safe choice for large λ, e.g., $\lambda >$ 0.1). Thus, the complexity reduces with respect to massive integration, and in fact, this method would be a good way to reduce the computation complexity. Alternatively, whenever one uses a beam gain approximation model, such as the flat-top model [31], the summation reduces straightforwardly to two terms, and this would be quite fast in terms of computation for network analyses.

Remark 3.

The array gain (2) of vertical ULA is critical for the analytical tractability of the CF derivation, and the use of $N_v>$ 1 could be questionable for a height that is too small. The array gain (2) holds true when the array aperture is compact compared to the array height h to have a plane wave-front. In practice, for mmW communications at a frequency of approximately 30 GHz (or larger), the wavelength is approximately 1cm (or smaller), and for $N_v=$ 10 antennas, the array aperture for half-wavelength antenna spacing is 5 cm (or less). In the scenario in Figure 1 the height should be above the people heights, and for $h>2$, m the approximation that the array aperture is compact (2 m ≫ 5 cm) holds true.

3.2. Specific Cases

UCA: It is a special case of UcylA for $N_v=1$ (i.e., $G_v\left(\theta \right)=1$). The statistical distribution of aggregated interference I for UCA can be adapted by considering $R_{max}\to \infty$, although it can be extended to $R_{max}<\infty$. In UCA, there is no radial mitigation of the interference, but it is only along the azimuth. After simplifying the relation (15) and resolving the singularity (Appendix A), it yields:

$${\mathsf{\Psi }}_I\left(\omega \right)=exp\left(-\frac{\pi \lambda }{C_{\alpha }}{\left|\omega \right|}^{\alpha }(1-j\mathrm{sign}\left(\omega \right)tan\frac{\pi \alpha }{2})P_{G_c}(\omega ,\alpha ,h)+\pi \lambda h^2\right),$$

(19)

where

$$P_{G_c}(\omega ,\alpha ,h)=\frac{{\overline{\beta }}^{2\alpha }}{K}\sum _{k=1}^K{G}_c^{2\alpha }\left({\mathsf{\Phi }}_k\right)P\left(-\alpha ,-j\omega \frac{G_c^2\left({\mathsf{\Phi }}_k\right){\overline{\beta }}^2}{h^{2b}}\right).$$

(20)

The limit for $K\to \infty$ angular sectors is

$$P_{G_c}(\omega ,\alpha ,h)=\frac{{\overline{\beta }}^{2\alpha }}{2\pi }\underset{0}{\overset{2\pi }{\int }}{G}_c^{2\alpha }\left(\varphi \right)P\left(-\alpha ,-j\omega \frac{G_c^2\left(\varphi \right){\overline{\beta }}^2}{h^{2b}}\right)d\varphi ,$$

(21)

and this integral can be evaluated numerically. Considering $h\to 0$ as a special case, $P_{G_c}(\omega ,\alpha ,h)\to 1$, so the frequency dependence of $P_{G_c}(\omega ,\alpha ,h)$ vanishes, and the aggregated interference in the kth angular sector is skewed alpha-stable

$$G_c^2\left({\mathsf{\Phi }}_k\right)\sum _{i\in {\mathsf{\Phi }}_k}\frac{|{\beta }_i{|}^2}{R_i^{2b}}\sim \mathcal{S}\left(\alpha =\frac{1}{b},{\gamma }_k\right),$$

(22)

where ${\gamma }_k=\frac{\Delta {\varphi }_k}{2}\lambda G_c^{2\alpha }\left({\mathsf{\Phi }}_k\right)\frac{{\overline{\beta }}^{2\alpha }}{C_{\alpha }}$ and $C_{\alpha }^{-1}=\Gamma (1-\alpha )cos(\pi \alpha /2)$. $\mathcal{S}(\alpha ,\gamma )$ denotes the skewed stable distribution with characteristic exponent $\alpha \in (0,2]$, unitary skewness, and scale parameter (or dispersion) $\gamma \ge 0$ with a characteristic function

$$\mathbb{E}\left[e^{j\omega x}\right]=\left\{\begin{array}{cc}exp[-\gamma |\omega {|}^{\alpha }(1-j\mathrm{sign}\left(\omega \right)tan\frac{\pi \alpha }{2})]\hfill & \alpha \ne 1\hfill \\ {exp[-\gamma |\omega |}^{\alpha }(1+j\frac{2}{\pi }\mathrm{sign}\left(\omega \right)ln\left|\omega \right|\left)\right]\hfill & \alpha =1\hfill \end{array}\right..$$

The overall interference reduces to the sum of skewed stable random terms (straightforwardly from ([49], eq. (1.8)))

$$I=\sum _k{I}_k\sim \mathcal{S}\left(\alpha =\frac{1}{b},{\gamma }_c\right),$$

(23)

where the total dispersion for UCA becomes

$${\gamma }_c=\pi \lambda \frac{{\overline{\beta }}^{2\alpha }}{C_{\alpha }}·\frac{{\int }_0^{2\pi }{G}_c^{2\alpha }\left(\varphi \right)d\varphi }{2\pi },$$

(24)

assuming sectors $\Delta {\varphi }_k\to 0$. Therefore, the aggregated interference for UCA and $h=0$ is skewed alpha-stable. However, when increasing the height h, the distribution deviates from alpha-stable, as detailed later. Comparing this result with [11], one notices an additional term that depends on the UCA array gain $G_c^2\left(\varphi \right)$, which mitigates the mean level of interference in skewed stable distribution. As before, the fading powers of the interferers are embodied in the derivations, but they will average out in summation of the interference power over all of the UEs (the fading amplitude of each signal from each UE depends on the chosen fading model), while the fading power from the user of interest remains effective in the calculation of the service probability.

Point antenna: a single point antenna can be considered as a special case of a UCA, where $N_c=1$ leads to isotropic gain. The radiation gain $G_c\left({\mathsf{\Phi }}_k\right)=1$ for any azimuth, and the choice of $K=1$ can be chosen, as it is not necessary to decompose the area into different cones. Placing the antenna at height $h>0$, the CF of the aggregated interference power follows by simplifying (19) as

$${\mathsf{\Psi }}_I\left(\omega \right)=exp\left(-\frac{\pi \lambda }{C_{\alpha }}{\left|\omega \right|}^{\alpha }{\overline{\beta }}^{2\alpha }P(-\alpha ,-j\omega \frac{{\overline{\beta }}^2}{h^{2b}})(1-j\mathrm{sign}\left(\omega \right)tan\frac{\pi \alpha }{2})+\pi \lambda h^2\right).$$

(25)

This CF (25) generalizes the CF for $h=0$ in [11]. The term $P(-\alpha ,-j\omega \frac{{\overline{\beta }}^2}{h^{2b}})$, which depends on $\omega$, increases from the initial value for $h=0$.

3.3. Analysis of AP Height

To gain an insight into the effect of the height of the antenna, it is useful to evaluate the mean aggregated interference that, for simplicity, is for UCA. The mean $\mathbb{E}\left[I\right]$ follows from the CF properties (check Appendix B for derivation):

$$\mathbb{E}\left[I\right]=\frac{\pi \lambda {\alpha }^2}{(1-\alpha )}{h}^{2(1-b)}\overline{{\beta }^2}\overline{G_c^2}.$$

(26)

where $\overline{G_c^2}={\frac{1}{2\pi }}_0^{2\pi }{G}_c^2\left(\varphi \right)d\varphi$ is the mean power gain. The first and second derivatives of the mean interference power with respect to the AP height are

$$\begin{array}{cc}\hfill \frac{d\mathbb{E}\left[I\right]}{dh}& =-2\pi \lambda \alpha \overline{{\beta }^2}\overline{G_c^2}{h}^{1-2b}<0,\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \frac{d^2\mathbb{E}\left[I\right]}{d{h}^2}& =-2\pi \lambda (\alpha -2)\overline{{\beta }^2}\overline{G_c^2}{h}^{-2b}>0.\hfill \end{array}$$

(28)

The aggregated interference power decreases vs. h but increases with UE density $\lambda$, with the second derivative always being positive. On the other hand, the received signal power $P_{rx}$ from a target UE located at distance $R_0$ is

$$P_{rx}\left(h\right)=\frac{|{\beta }_0{|}^2}{{(R_0^2+h^2)}^b}.$$

(29)

The first and second derivatives with respect to the AP height are

$$\begin{array}{cc}\hfill \frac{d{P}_{rx}\left(h\right)}{dh}& =\frac{-2b|{\beta }_0{|}^{2}h}{{(R_0^2+h^2)}^{b+1}}<0,\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill \frac{d^2{P}_{rx}\left(h\right)}{d{h}^2}& =\frac{2b|{\beta }_0{|}^2\left((2b+1)h^2-R_o^2\right)}{{(R_0^2+h^2)}^{b+2}}.\hfill \end{array}$$

(31)

The first derivative is always negative. The second derivative is negative for $R_0 > h\sqrt{(2b+1)}$ and positive for $R_0 < h\sqrt{(2b+1)}$. Comparing the first and second derivatives of the mean aggregated interference power and received target power, given a specific set of parameters, such as the target user distance $R_0$, the signal-to-interference ratio (SIR) would have a concave behavior vs. the AP height h. This is numerically evaluated in the next section.

3.4. Numerical Validation on Aggregated Interference

The CF derived for each of the arrays is validated here by numerically computing the PDF and CDF using inverse Fourier methods tailored for statistical distributions to be accurate on the tails of the distributions [50,51]. The service probability (6) is the comparison metric adopted here for the validation of the CFs in the previous sections, by considering a LOS system with ${|{\beta }_o|}^2=1$ and a threshold of $T=1$ (or 0 dB). The amplitudes are normalized as $|{\beta }_o|=\mathbb{E}\left[\right|\beta \left|\right]=1$, since in the analysis of the SIR, the transmit power and the path-loss at the reference distance are not important as they are equal for the interferers and user of interest.

Remark 4.

Since the CF (15) is based on $R_{max}\to \infty$, it is impossible to validate it numerically. Instead, a maximum numerical radius $R_{max}^{\left(num\right)}$ is used. To match the simulation results with analytical ones, we limit the integral in (12) to $R_{max}=R_{max}^{\left(num\right)}$. The choice of a proper $R_{max}$ for a numerical example to explore the behavior of an infinite area is discussed in [48,52].

Figure 2 demonstrates the comparative analysis of UcylA and UCA, showing the $P_s\left(R_o\right)$ vs. radius $R_o$ using ${\mathsf{\Psi }}_I\left(\omega \right)$ for UcylA in (15) and for UCA in (19) where the total number of antennas is preserved in all cases $(N_c\times N_v=128)$ and the SIR threshold ${\left[T\right]}_{dB}=0$. The transmitting interferers are numerically generated as random PPP with a maximum radius $R_{max}^{\left(num\right)}$, as specified below for every Monte Carlo iteration, and are affected by the array gain $G(\varphi ,\theta )=G_c\left(\varphi \right)G_v\left(\theta \right)$ (see (1) and (2)) keeping the radial position fixed at $R_o$ for the UE of interest as aggregated interference is isotropic vs. azimuth. The analytical curves (shown in solid lines) are for $R_{max}=300$ m (unless otherwise mentioned), which are validated by simulation curves (shown in dashed lines with markers). Here, for the purpose of illustration, we have also shown the analytical curve (in yellow color) with $R_{max}\to \infty$, for the case of UCA with $\lambda ={10}^{-3}$, that lower bounds a simulation curve with $R_{max}^{\left(num\right)}=2000$ m. Although, in case of UCA, the analytical curves are a lower bound of $P_s\left(R_o\right)$ when compared to the numerical simulation with $R_{max}^{\left(num\right)}\ll R_{max}$, this margin is not considerable in the case of UcylA, since this array automatically cancels out the effect of far-away UEs. Service probability $P_s\left(R_o\right)$ decreases for increasing served UE position $R_o$ as interference from $\lambda$-density interferers dominates. For smaller density (here $\lambda ={10}^{-3}{\mathrm{m}}^{-2}$), the service probability $P_s\left(R_o\right)>0.5$, up to $70$ m for UCA and around 50 m for UcylA. It can be observed that in the given scenario with given parameters, the UCA seems to surpass the UcylA from the service probability point of view.

Increasing the path-loss exponent ($2b=3.6$ in Figure 3) affects the service probability as aggregated interference is more attenuated for far-away interferers, and it is more effective for denser users (i.e., for larger $\lambda$, the increase in the path-loss is more beneficial for $P_s\left(R_o\right)$, while it is detrimental for small $\lambda$). The rest of the figures in the paper are based on analytical formulas, and simulation results are not reported for the sake of the readability of the figures.

As an illustrative example of the optimum height based on relations (26) to (31), Figure 4 shows the received signal power and the mean aggregated interference power, where the transmit power of all the UEs is ${\left[P_{tx}\right]}_{dBm}=20$, and the normalized array gain for the user of interest located at distance $R_0$ is maximum (i.e., $G_c\left({\varphi }_0\right)=1$). The amplitude ${\beta }_o$ also considers the path-loss at distance 1 m: $|{\beta }_o|=\sqrt{P_{tx}}4\pi F_c/c$ for $c=3\times {10}^8$ m/s and $F_c$ is the carrier frequency (here $F_c=28$ GHz). The mean interference power follows (26). It can be observed that for large thresholds ${\left[T\right]}_{dB}=\{4,5\}$, the average interference power is larger than the signal power, and the system is in outage for every AP height. For smaller thresholds, the AP serves the target UE for a range of AP heights. For example, at ${\left[T\right]}_{dB}=0$, the serving range of AP height is approximately $1\mathrm{m}<h<33$ m. At around $h=1$ m, the difference of target UE signal power and mean aggregated interference is zero. For the AP height range $1\mathrm{m}<h<9\mathrm{m}$, this difference increases, and for the AP height range of $9\mathrm{m}\le h<33\mathrm{m}$, the difference decreases, while for $33\mathrm{m}<h$, the user is in outage. Finally, for a better demonstration of the concave behavior of SIR vs. AP height, Figure 5 is plotted, for two groups of UE density $\lambda =\{{10}^{-4},{10}^{-3}\}{\mathrm{m}}^{-2}$ and for the target UE distance of $R_0=\{20,30,40\}$ m. Three SIR threshold ${\left[T\right]}_{dB}=\{5,10,15\}$ guidelines are plotted as horizontal lines. It can be observed that for $\lambda ={10}^{-3}$, and a pre-defined threshold of ${\left[T\right]}_{dB}=5$, a target user located at $R_0=40$ m is always in the outage, while the target users located at $R_0=\{20,30\}$ m would experience a concave behavior vs. the AP height h.

Once the analytical model has been validated, one might investigate the $N_c$ vs. $N_v$ arrangement of UcylA for a given total number of antennas $N_c{N}_v$ (e.g., for the same complexity of the radio frequency circuitry). The cylinder arrangement of the UcylA can be tall ($N_c<N_v$), fat ($N_c>N_v>1$) or just a ring ($N_v=1$), and the optimum array geometry for service probability depends on different parameters, such as the SIR threshold (T), path-loss exponent (b), the antenna (or users) height h, and the directivity of every antenna element (not considered here). The metric used herein is the ratio of the served users in a certain area ($2\pi \lambda {\int }_0^{\overline{R}}{P}_s\left(r\right)rdr$) to the total average users ($\pi {\overline{R}}^2\lambda$), assumed here to be PPP distributed, which is also referred as average service probability

$${\overline{P}}_s=\frac{2{\int }_0^{\overline{R}}{P}_s\left(r\right)rdr}{{\overline{R}}^2}.$$

(32)

In the following examples, we maintain the total number of antennas $N_c\times N_v=256$, changing the ratio of UCA and vertical antennas. Figure 6 illustrates the number of users served vs. the antenna ratio ${log}_2(N_c/N_v)$, varying the path-loss exponent $2b=\{2,2.4,2.8,3.2,3.6\}$ for a medium UE density $\lambda ={10}^{-3}{\mathrm{m}}^{-2}$ and an extremely large UE density $\lambda =5\times {10}^{-2}{\mathrm{m}}^{-2}$. It can be seen that, in the case of a smaller UE density, it is preferable to use a UCA. In the case of a very large UE density, a UcylA would be preferred. This is because for a smaller UE density, there are not many interferers that are far from the AP, while their power is largely attenuated w.r.t interferers located closer to the AP. Thus, UCA outperforms by creating a narrower main lobe toward the interested UE, nullifying most of the other UEs located close to the AP but in different azimuths. However, when $\lambda$ is large, there are too many far-away interferers. Although the interference power of each one of them is highly attenuated, the summation of many far-away interferers is still too high. In this case, it is preferable to have a vertical beamforming component, which can filter out the effect of far-away UEs by creating a ring-shaped gain pattern around the AP. Furthermore, based on Figure 6, it can be seen that the increase in the path-loss exponent has an inverse effect on the two sets of curves. For the case of a small UE density, by increasing the path-loss exponent, the number of UEs served decreases, while in the case of a very dense network ($\lambda =5\times {10}^{-2}{\mathrm{m}}^{-2}$), by increasing the path-loss exponent, the number of UEs served increases. This is because, in the case of dense networks, higher path-loss is beneficial by highly attenuating the interference coming from the large crowd of interferers that are located far away, while in the case of a small UE density, high path-loss attenuates to a larger extent than the signal coming from the target user. It can be noted that for $2b=2$, the curve monotonously decreases, and it is the minimum for a UCA, while the other curves have a concave shape. This is due to the fact that $2b=2$ is a singularity point when the area where the interferers is infinite (i.e., $R_{max}\to \infty$). This can be inferred by setting $\alpha =1$ in relation (26), or from the literature on aggregated interference, e.g., [17,48]. This means that for $\alpha =1$, the average aggregated interference power diverges, and the AP does not serve any UE. In the example of Figure 6, the radius is $R_{max}=400$ m, and thus, the service probability is not zero.

Figure 7 illustrates the average number of the users served vs. the antenna ratio ${log}_2(N_c/N_v)$, varying the SIR threshold ${\left[T\right]}_{dB}=\{0,5,10\}$ for small and large path-loss exponents $2b=\{2,3.6\}$. It can be seen that in this scenario, by increasing the threshold T, UCA performs worse. Note that in this setup, the UEs are very dense ($\lambda =5\times {10}^{-2}{\mathrm{m}}^{-2}$). It can be shown that for smaller UE densities, UCA outperforms UcylA. Figure 8 shows a similar example of the average number of users served vs. the antenna ratio ${log}_2(N_c/N_v)$, fixing the threshold ${\left[T\right]}_{dB}=5$ and varying the AP height $h=\{5,10,15,20\}$ m. Clearly, by increasing the AP height, the usage of UcylA becomes more advantageous in terms of ${\overline{P}}_s$, and it is more meaningful to use vertical beamforming, while for smaller heights, it is preferable to use a larger circular array. However, it can be observed that, in this scenario, where UEs are dense, the service probability, in general decreases by increasing the height. Thus, it is preferable to use an AP with a smaller height.

Figure 9 shows the effect of UE density $\lambda$. It is seen that for areas that are not very dense, a UCA is marginally more preferable. However for more dense areas, a UcylA would be more preferable depending on different parameters. It can be shown that changing every one of the above-mentioned parameters affects the shape of the ${\overline{P}}_s$ curves. Furthermore, the constraint of ground UEs is also responsible for the results. UcylAs are favored when UEs have arbitrary heights (usually indoor UEs are at arbitrary heights, while when targeting high frquencies, penetration loss is too high).

A pragmatic conclusion from these evaluations is that for most real-life scenarios where UE densities are not too high, it is beneficial to invest in UCA arrangement rather than UcylA. There are two reasons supporting this result: (i) the most powerful interferers are in the vicinity of the AP, while the minimum separable angle by the ULA is $\Delta {\theta }_{min}\approx \frac{\lambda }{N_v{d}_{v}cos\left({\theta }_0\right)}$, where $d_v$ is the vertical inter-element spacing, and ${\theta }_0$ is the target elevation; (ii) although the resolution was not angle dependent, decomposing the whole elevation plane into small portions of the same width, provides emphasis in the same as near and to far-away UEs, while the most powerful UEs are closer ones. In the rest of the numerical examples in the paper, we focus on the usage of UCA, since the target threshold used is set to ${\left[T\right]}_{dB}=0$ and the interferer density used is not extremely high, which justifies the usage of a UCA over UcylA.

Remark 5.

Note that although the leftmost parts of the curves (i.e., corresponding to a ULA or a very tall UcylA) are shown in the figures, in practice, they are not feasible to deploy. Therefore, one might consider the range $N_c>N_v/$8 as practical solution.

4. Statistical Characterization of the Aggregated Interference Power

The aggregate interference I in Section 3.1 for UcylA and arbitrary height $h>0$ is complex to be computed analytically in a closed-form. Herein, we propose a method for the approximation of the aggregate interference CF. We show that the aggregated interference power for an array of antennas located at an arbitrary height can be approximated using a weighted mixture of two stable distributions, and we detail herein the equivalent CF.

In order to gain a deeper insight into the distribution at an arbitrary height, one can start from the Taylor series of the argument of the CF (${\mathsf{\Psi }}_I\left(\omega \right)$). For example, the series for a UCA with arbitrary height and $|{\beta }_i|=1$, follows from the CF (19) that, with some simplifications, yields:

$${\mathsf{\Psi }}_I\left(\omega \right)=exp\left(\pi \alpha \lambda \Xi (\alpha ,\omega )\right),$$

(33)

where the Taylor series is

$$\Xi (\alpha ,\omega )=\frac{1}{K}\sum _{k=1}^K\sum _{i=1}^{\infty }\frac{{\left(\frac{1}{h^{2b}}\right)}^{i-\alpha }}{(i-\alpha )i!}{\left(j\omega {\overline{\beta }}^2{G}_c^2\left({\mathsf{\Phi }}_k\right)\right)}^i.$$

(34)

Let us define $G_{avg}^{\left(z\right)}=\frac{1}{K}{\sum }_{k=1}^K{G}_c^{2z}\left({\mathsf{\Phi }}_k\right)$, which for $h=0$ relation (34) simplifies to

$$\Xi (\alpha ,\omega )={(-j\omega )}^{\alpha }\mathsf{\Gamma }(-\alpha ){\overline{\beta }}^{2\alpha }{G}_{avg}^{\left(\alpha \right)}.$$

(35)

For very large heights $h\gg 0$ in (33), the terms with a higher index i are negligible, and it can approximated with only the first two terms, which would make the CF (${\mathsf{\Psi }}_I\left(\omega \right)$) have a Gaussian distribution (the CF of a Gaussian distribution is $exp\left(j\mu \omega -\frac{{\sigma }^2}{2}{\omega }^2\right)$ with $\mu$ being the shift parameter and $\sigma$ is the standard deviation). In order to evaluate the appropriate CF for any $h>0$, one isolates the behavior vs. $j\omega$ from $\Xi (\alpha ,\omega )$ as:

$$\Xi (\alpha ,\omega )={\Xi }^{\prime }(\alpha ,\omega )+\frac{{\left(\frac{1}{h^{2b}}\right)}^{1-\alpha }}{1-\alpha }j\omega {\overline{\beta }}^2{G}_{avg}^{\left(1\right)},$$

(36)

where the second term (corresponding to $i=1$ in (33)) is a shift or location parameters. The behavior of the real part of $\Xi {(\alpha ,\omega )}^{\prime }$ reveals the corresponding exponent of $\omega$, for every defined $\omega$. The reason behind separating $\Xi (\alpha ,\omega )$ into two parts is that we need to omit the shift, to be able to visualize the exponent of the $\omega$ within the distribution. For large h, ${\Xi }^{\prime }(\alpha ,\omega )$ vs. $\omega$ behaves as ${\omega }^2$, and for very small h, behaves as ${\omega }^{\alpha }$, while for medium heights, it has two different slopes based on $\omega$. The transition $\omega$ where the behavior changes is $\overline{\omega }$, which depends on the height, the path-loss exponent, UE density and array gain. Having gained insight into the behavior vs. $\omega$ and $\overline{\omega }$, which is explained further on in the text, $\Xi (\alpha ,\omega )$ can be approximated as:

$$\Xi (\alpha ,\omega )=W_1\left(\omega \right)\times a\left(\omega \right)+W_2\left(\omega \right)\times \left(1-a\left(\omega \right)\right),$$

(37)

where $a\left(\omega \right)$ is a Heaviside step function, i.e., $a\left(\omega \right)=1$ for $\omega <\overline{\omega }$ and 0 otherwise (or some function with smoother transition), which acts as a switch between two cases with different behavior:

$$\begin{array}{ccc}\hfill W_1\left(\omega \right)& \approx & \frac{{\left(\frac{1}{h^{2b}}\right)}^{1-\alpha }}{1-\alpha }j\omega {\overline{\beta }}^2{G}_{avg}^{\left(1\right)}-\frac{{\left(\frac{1}{h^{2b}}\right)}^{2-\alpha }}{2(2-\alpha )}{\omega }^2{\overline{\beta }}^4{G}_{avg}^{\left(2\right)},\hfill \end{array}$$

(38)

$$\begin{array}{ccc}\hfill W_2\left(\omega \right)& =& {(-j\omega )}^{\alpha }\mathsf{\Gamma }(-\alpha ){\overline{\beta }}^{2\alpha }{G}_{avg}^{\left(\alpha \right)},\hfill \end{array}$$

(39)

where $W_2$ coincides with the skewed-stable distribution (23). Relation (37) means that the CF is can be decomposed as

$${\mathsf{\Psi }}_I\left(\omega \right)=exp\left(\pi \alpha \lambda W_1\left(\omega \right)a\left(\omega \right)\right).exp\left(\pi \alpha \lambda W_2\left(\omega \right)(1-a\left(\omega \right))\right).$$

(40)

Figure 10 demonstrates the behavior of the real part of $\Xi (\alpha ,\omega )$ vs. $\omega$ for a single isotropic antenna and set of AP heights $h=\{0,2,5,20,100\}$ m. It can be seen that by increasing the AP height, the breaking frequency $\overline{\omega }$ increases, while for extremely large AP heights, it tends to infinite, which is Gaussian distribution behavior. Figure 11 is the same analysis, comparing an isotropic antenna with a UCA consisting of $N_c=16$ isotropic antennas for the set of AP heights $h=\{0,5\}$ m. It is observed that the slopes of the curves are maintained, while the breaking frequency is increased when using a UCA.

One can approximate the CF for a UCA in every arbitrary height, given the knowledge on the transition frequency $\overline{\omega }$. Knowing the transition point $\overline{\omega }$, based on height h, one may characterize different $\omega$ with different statistical distributions. The numerical formulation of $\overline{\omega }$ vs. antenna height and other parameters can be investigated further. A rule of thumb for $\overline{\omega }$ can be achieved as follows. Let $S_i\left(\omega \right)$ denote the i-th term of the series (34) as:

$$S_i\left(\omega \right)=\frac{{\left(\frac{1}{h^{2b}}\right)}^{i-\alpha }}{(i-\alpha )i!}{\left(j\omega {\overline{\beta }}^2{G}_{avg}^{\left(i\right)}\right)}^i.$$

(41)

It is empirically observed that $\overline{\omega }$ can be approximately achieved by imposing the condition $S_2\left(\omega \right)=S_3\left(\omega \right)/{\overline{\beta }}^2$, and by solving for $\omega$, the $\overline{\omega }$ can be achieved. Please note that this formula holds for normalized array gain. The rationale behind these conditions lies in the fact that for an $\alpha$-stable distribution and for height $h=0$ m, all the terms for $i>1$ tend to infinite, while for Gaussian distribution, only the first three terms exist. In the mixture case, for $\overline{\omega }<\omega$ behavior resembles $\alpha$-stable. In the case of Figure 10, the calculated $\overline{\omega }$ for the AP heights of $h=\{2,5,10,100\}$ m are $\overline{\omega }=\{32,355,13070,858160\}$ respectively, and in Figure 11, the $\overline{\omega }$ for AP height $h=5$ m for two cases of $N_c=\{1,16\}$ are $\overline{\omega }=\{355,5688\}$, respectively. Based on the figures, it can be noted that these approximations are close to the real breaking points of the curves.

Figure 12 and Figure 13 (corresponding to AP height $h=2$ m and $h=30$ m, respectively) validate the proposed approximate model, with respect to true CF of a UCA described in (19) and the $\alpha$-stable distribution model [11] that coincides with (19) for $h\to 0$. It can be seen that the stable distribution might be appropriate only if the AP height and UE density are not very large. By increasing the UE density $\lambda$, for large AP height h, the approximate mixture model (40) outperforms the stable model. The difference in the median between the approximate model (dashed curves) and the true model (solid curves) in the dB scale, for $h=2$ m, is in the order of $0.1$ to $0.3$ dB, and for $h=30$ m, it is in the order of $1.1$ to $1.3$ dB, while the curves corresponding to stable distribution are dashed lines with a marker. From the computation complexity point of view, depending on the parameters, the approximate model (40) is computed very fast compared to (33). The computational complexity of Monte-Carlo simulations, on the other hand, depend on the UE density $\lambda$, while the true and approximate analytical models do not depend on the density $\lambda$.

In practical systems, the height is known, but other parameters, such as the density of active users $\lambda$, are not known, and some inaccuracies with respect to the ideal model might occur. We believe that the knowledge of a reasonable approximation of the distribution of the aggregated interference enables the measurement of the approximating alpha-stable distribution during multiple idle times of the communication intervals using any unsupervised learning method [53], which is a practical on-the-fly method.

5. Outage Analysis in the Presence of NLOS, Noise and Blockage

In the previous sections, we derived the CF of the aggregated interference power with different antenna array configurations, and characterized the distribution of the aggregated interference power. Moving toward the modeling of practical mmW and sub-THz systems, in this section, it is shown how the NLOS propagation, noise power, and blockage can be integrated into the model.

5.1. NLOS Paths

Previous studies and measurement campaigns [54] have shown that the NLOS clusters of rays are present in mmW communications where they form sparse multipath faded channels. These paths can increase the amount of interference, but at the same time would lead to a more useful signal being received in the case of the coherent reception of a signal. On the other hand, diversity is an efficient way to compensate for the blockage and to increase the reliability of communication systems.

5.2. NLOS Paths

Let L be the total number of paths that the signal arrives from the user of interest to the same AP; the $L\times 1$ set of signals $\mathbf{y}={[y_1,y_2,\dots ,y_L]}^T$ after the multi-beam beamforming to each of the paths from the user of interest transmitting x is

$$\mathbf{y}=\mathbf{h}x+\iota +\mathbf{w},$$

(42)

where ${\left[\mathbf{h}\right]}_{\ell }=h_{\ell }={\beta }_{o,\ell }/D_{\ell }^b$ for distance $D_{\ell }={(R_{o,\ell }^2+h^2)}^{1/2}$ corresponding to the ℓth paths of arrival in case $\ell =1$ it is the direct LOS link and $R_{o,1}=R_o$ is the geometric distance between the user of interest, while the distances for NLOS links are modeled later but $R_{o,\ell }\ge R_o$. The ensemble of the aggregated interference amplitudes from the PPP-distributed interferers is $\iota ={[{\iota }_1,{\iota }_2,\dots ,{\iota }_L]}^T$, which are independent and identically distributed (iid) random variables obtained from a set of L beamforming toward the distinct angles ${\varphi }_{o,1},{\varphi }_{o,2},\dots ,{\varphi }_{o,L}$ for LOS (${\varphi }_{o,1}$) and NLOS (${\varphi }_{o,2},\dots ,{\varphi }_{o,L}$) of the user of interest, so that adapting (4) to this case with multipath for interference is:

$$\iota =\sum _{i=1}^{\infty }\sum _{\ell =1}^L\frac{{\beta }_{i,\ell }}{{\left(R_{i,\ell }^2+h^2\right)}^{\frac{b}{2}}}{G}_c\left({\varphi }_{i,\ell }\right)x_{i,\ell }.$$

(43)

This assumption is justified by the interfering ray-paths on every beamforming that has different attenuations and phase shifts, and is thus independent. $\mathbf{w}={[w_1,w_2,\dots ,w_L]}^T$ is the collection of noise amplitudes. The LOS/NLOS links are shown in Figure 14 and using the Weyl model (similar to Saleh-Valenzuela [55] adapted for mmW [56,57]), where the indirect NLOS paths from the transmitter are reflected from a secondary point, which is uniformly distributed around a circle with a radius d around each of the transmitter’s location. Usually, at high frequencies, there are not many NLOS paths, which are typically $L=2-3$ [54,58,59]. The receiver for L paths, possibly with different (and likely delay-resolvable for mmW and sub-THz system with large bandwidth) delays, is expected to combine to maximize the service probability. If using the Maximal Ratio Combining (MRC) combiner for the L paths related to the user of interest affected by the interference powers $I_{\ell }={\left|{\iota }_{\ell }\right|}^2$, one obtains the following service probability analysis (see Appendix C for derivation and specific MRC notation):

$$P_{MRC}={\mathbb{P}}_I\left(\sum _{\ell =1}^L{\alpha }_{\ell }{I}_{\ell }<\frac{1}{T}{\left(\sum _{\ell =1}^L\frac{|{\beta }_{\ell }{|}^2}{{\overline{I}}_{\ell }^2{D}_{\ell }^{2b}}\right)}^2\right).$$

(44)

where $|{\beta }_{\ell }{|}^2$ is the fluctuation for the signal from the user of interest from the ℓth path, and the distribution of total interference $I_{MRC}={\sum }_{\ell =1}^L{\alpha }_{\ell }{I}_{\ell }$ follows on from the CF for iid interferers over the L-beamformers:

$${\mathsf{\Psi }}_{I_{MRC}}\left(\omega \right)=\prod _{\ell =1}^L{\mathsf{\Psi }}_I\left(a_{\ell }\omega \right).$$

(45)

The NLOS distances $D_{\ell }={(R_{o,\ell }^2+h^2)}^{1/2}$, and $R_{o,\ell }$ for $\ell >1$ are affected by the random angular position $\eta$ of the NLOS reflections around the radius d, as depicted in Figure 14:

$$R_{o,\ell }=d+\sqrt{R_{o,1}^2+d^2-2d{R}_{o,1}cos\left(\eta \right)}\forall \ell >1.$$

(46)

Pragmatically, the NLOS distance for the served user is dependent on the specific multipath model assumed here, and a convenient way to incorporate the NLOS attenuation into the additional path $R_{o,\ell }-R_{o,1}$ is to approximate this term is by considering the mean distance for NLOS ${\overline{R}}_{o,\ell }=E_{\eta }\left[R_{o,\ell }\right]$. The distance ${\overline{D}}_{\ell }={({\overline{R}}_{o,\ell }^2+h^2)}^{1/2}$ is for the NLOS $\ell >1$, and thus the relationship (44) is somewhat simplified by constant distances. Recalling that LOS/NLOS models hold for interferers’ $i_{\ell }$, the power is augmented by the NLOS components, and the service probability is reduced accordingly compared to $L=1$ (LOS-only). The CF for $L>1$ is derived in Appendix D, accounting for LOS and NLOS.

Figure 15 and Figure 16 show the average number of users served $M_s$ within a $\Delta \times \Delta$ square-shaped area (from (32) it is $M_s=\lambda {\Delta }^2{\overline{P}}_s$) when the AP with UCA is located at the center. The multipath increases the interference and the number of NLOS paths from 1 to 3 (or $L=1,2,3,4$). In Figure 15, the performance degrades as the density of the interferers is quite high (here $\lambda ={10}^{-2}{\mathrm{m}}^{-2}$); in contrast, in Figure 16, the aggregated interference is lower due to the lower density ($\lambda ={10}^{-3}{\mathrm{m}}^{-2}$), and thus, there is a clear benefit arising from the multipath that vanishes for a large height (here $h>25$ m). The performance of MRC shown in Figure 15 and Figure 16 is compared to the selection combining (SC), which select the path with the largest SIR ([60], Ch. 7). The derivation of SC is straightforward (not shown here). As expected, the MRC outperforms the SC, but one might observe that degradation becomes more severe for large multipath. In the remainder of the paper, we consider the multipath condition with both LOS and NLOS links.

5.3. Noise

Although the paper analyzes the aggregated interference power, in a 6G network with high path-loss it is inevitable to consider also the effect of the noise power. The SINR is defined as:

$$\mathrm{SINR}=\frac{|{\beta }_0{|}^2{G}^2({\varphi }_0,{\theta }_0)}{(I+{\sigma }_w^2/N){\left(R_o^2+h^2\right)}^b},$$

(47)

where $I={\left|\iota \right|}^2$ and $\iota$ is defined in (44). For the service probability in the presence of noise, it is easy to prove that relation (44) must be modified as

$${\mathbb{P}}_{service}={\mathbb{P}}_I\left(I_{tot}+N_{tot}<\frac{1}{T}{\left(\sum _{\ell =1}^L\frac{|{\beta }_{\ell }{|}^2}{{\overline{I}}_{\ell }{D}_{\ell }^{2b}}\right)}^2\right),$$

(48)

where $I_{tot}={\displaystyle \sum _{\ell =1}^L}{\alpha }_{\ell }{I}_{\ell }$ and $N_{tot}={\sigma }_n^2{\displaystyle \sum _{\ell =1}^L}{\alpha }_{\ell }$ are the total interference power and noise power respectively, after MRC at the AP with ${\alpha }_{\ell }={\left|{\beta }_{\ell }\right|}^2/{({\overline{I}}_{\ell }+\frac{{\sigma }_w^2}{N})}^2{D}_{\ell }^{2b}$.

5.4. Blockage

In mmW and 6G sub-THz systems, the waves are prone to blockage due to static and dynamic blockage. Static blockage [61,62,63] is caused by structures such as buildings and trees; self-blockage [63,64,65] is caused by the body holding the UE; and dynamic blockage [63] is caused by moving objects, humans or vehicles. The behavior of blockages and their impact on the coverage probability and system performance are different. Some of the blockage effects can be modeled in closed form, while some others need numerical methods. Here, we limit the analysis of the paper to the numerical evaluation of the coverage probability in the presence of blockage.

It is convenient to define a dummy binary variable ${\mu }_{\ell }$ for the ℓ-th path, which is ${\mu }_{\ell }=0$ when the link is blocked and ${\mu }_{\ell }=1$ when the link is free of any blockage, such that $Prob({\mu }_{\ell }=0)=P_B$ is invariant on every path, and $P_B$ is the probability of blockage. The blockage per path can be incorporated into (48) as

$${\mathbb{P}}_{service}={\mathbb{P}}_I\left({\tilde{I}}_{tot}+{\tilde{N}}_{tot}<\frac{1}{T}{\left(\sum _{\ell =1}^L\frac{|{\beta }_{\ell }{|}^2{\mu }_{\ell }}{{\overline{I}}_{\ell }{D}_{\ell }^{2b}}\right)}^2\right),$$

(49)

where ${\tilde{I}}_{tot}={\displaystyle \sum _{\ell =1}^L}{\tilde{\alpha }}_{\ell }{I}_{\ell }{\mu }_{\ell }$ and ${\tilde{N}}_{tot}={\sigma }_n^2{\displaystyle \sum _{\ell =1}^L}{\tilde{\alpha }}_{\ell }$ are the total interference and noise power computed respectively, after MRC at the AP, with ${\tilde{\alpha }}_{\ell }={\mu }_{\ell }{\alpha }_{\ell }$ being a scaling parameter for a given path, and it is obviously switched off by ${\mu }_{\ell }$ when the path is blocked. The distribution of the aggregated interference power I, can be achieved by considering that blockage modifies the average density of interferers to $\lambda (1-P_B)$.

The goal is to assess the effect of the impact of blockage $P_B$ on service probability. At each snapshot, one link can be either blocked or available. In order to calculate the service probability, taking into account the blockage probability, one must numerically evaluate (49), for $Prob({\mu }_{\ell }=0)=P_B$ for $\forall \ell$.

Figure 17 is the average service probability and Figure 18 is the number of UEs served; both figures are versus blockage probability $P_B$ for small ($\lambda =5\times {10}^{-3}{\mathrm{m}}^{-2}$) and large ($\lambda =5\times {10}^{-2}{\mathrm{m}}^{-2}$) UE density, and a varying number of paths ($L=1,4$), and path loss ($2b=2,2.8$). In Figure 17, for a small $\lambda$ and LOS path ($L=1$), the blockage probability reduces the service probability, while for $L=4$, the service probability is more robust even for large $P_B$, and this is due to the diversity of the multipath against the interference. On the other hand, when $\lambda =5\times {10}^{-2}{\mathrm{m}}^{-2}$, a multipath channel with $L=4$ severely degrades the performance. Numerical analysis shows the impact of the blockage on the service probability in Figure 17, and the average number of user served in Figure 18. The total number of users within a $\Delta \times \Delta$ square-shaped area is $M_s=\lambda {\Delta }^2{\overline{P}}_s$. This means that for small $P_B$ any increase of one decade of UE density results in a smaller variation in the number of users served. Recall that the number of UEs served in Figure 18 refers to the UEs allocated in the same spectrum, and thus by knowing the average number of UEs assignable on the same spectrum region, one can pre-design the largest number of users that a resource scheduler can expect to assign (not covered in this paper).

6. Conclusions

In this paper, we derived tractable expressions for the characteristic function of the aggregate interference power for the homogeneous distribution of UEs using the SG framework for $N_c\times N_v$ UcylA placed at an arbitrary height. We proved that for $h>0$, the distribution of the aggregated interference could be analytically approximated by a decomposable mixture of two distributions: skewed alpha-stable and Gaussian. The numerical analysis validates the results derived for the array configurations, including the case when $N_v=1$ for UCA. The analysis of the average service probability vs. UcylA height shows that there are different trade-offs to exploit. Furthermore, the appropriate array geometry depends on different environment and propagation parameters. The impact of multipath was evaluated analytically, thus showing the trade-off due to the increased aggregated interference and the diversity of the UE of interest. Blockage allows a realistic analysis of mmW and 6G systems, and the blockage analysis has showed that there are several design insights to exploit.

Future work could consider the extension to distributed antenna systems (DAS) with different types of multi-AP coordination. The availability of the aggregated interference distributions in analytic form opens the possibility to explore multi-AP cooperation, which is otherwise, possible via massive simulations.