A Comprehensive Analysis of 5G Heterogeneous Cellular Systems operating over κ-μ Shadowed Fading Channels

—Emerging cellular technologies such as those proposed for use in 5G communications will accommodate a wide range of usage scenarios with diverse link requirements. This will necessitate operation over a versatile set of wireless channels ranging from indoor to outdoor, from line-of-sight (LOS) to non-LOS, and from circularly symmetric scattering to environments which promote the clustering of scattered multipath waves. Unfortunately, many of the conventional fading models lack the ﬂexibility to account for such disparate signal propagation mechanisms. To bridge the gap between theory and practical channels, we consider κ – μ shadowed fading, which contains as special cases the majority of the linear fading models proposed in the open literature. In particular, we propose an analytic framework to evaluate the average of an arbitrary function of the signal-to-noise-plus-interference ratio (SINR) over κ – μ shadowed fading channels by using an orthogonal expansion with tools from stochastic geometry. Using the proposed method, we evaluate the spectral efﬁciency, moments of the SINR, and outage probability of a K -tier heterogeneous cellular network with K classes of base stations (BSs), differing in terms of the transmit power, BS density, shadowing, and fading characteristics. Building upon these results, we provide important new insights into the network performance of these emerging wireless applications while considering a diverse range of fading conditions and link qualities.


I. INTRODUCTION
T O MEET the ever-increasing demand for data on the move, all major telecommunications companies, as well as global standardization entities, are actively driving the research and development of the fifth generation (5G) of wireless communications.It is forecast that this new networking paradigm will provide 1000 fold gains in capacity over the next decade and data rates exceeding 10 Gigabit/s while achieving latencies of less than 1 millisecond [1], [2].To make this possible, 5G communications will utilize densely deployed small cells to achieve high spectral efficiency while harnessing all available spectrum resources, including opportunities offered by millimeter-wave frequencies.Key to the successful operation of 5G communications will be the unification of dissimilar networking technologies.This will create a diverse range of link requirements and the necessity for wireless devices to operate over a versatile set of channels ranging from indoor to outdoor, from line-of-sight (LOS) to non-LOS (NLOS), and from homogeneous diffuse scattering to those which promote the clustering of scattered multipath waves.
A range of tools developed within the framework of stochastic geometry have been used to capture the irregularity and heterogeneity of 5G wireless networks with considerable success.Specifically, stochastic geometry assumes that the locations of all wireless nodes are endowed with a spatial point process [3].Such an approach captures the topological randomness in the network geometry, allows the use of well-established mathematical tools, offers high analytical flexibility and achieves an accurate performance evaluation [4].A common assumption made within this scheme is that the nodes are distributed according to a Poisson point process (PPP).Using this supposition, the probability density function (PDF) of the aggregate interference and the outage probability were analyzed for cellular networks in [5] and [6], which were then generalized to the case of heterogeneous cellular networks (HetNets) in [7]- [11]. 1uch of the existing published work on stochastic geometry has focused on the Rayleigh distribution as the de facto small-scale fading model, owing to its simplicity and tractabil-ity.Several approaches have been proposed to derive the signal-to-noise-plus-interference ratio (SINR) distributions for general fading environments.For instance, in [15]- [19] the conversion method, which is based on displacement theorem, was used.This method treats the channel randomness as a perturbation in the location of the transmitter and transforms the original network with arbitrary fading into an equivalent network without fading.Although the conversion method can be applied to any fading distribution, it is more tractable for handling large-scale shadowing effects.Specifically, if one applies the conversion method to small-scale fading, the resulting equivalent model will have no fading, which limits the use of Laplace transform-based approach.An alternative approach to the analysis of general fading scenarios relies on the series representation method [20], [21].This approach expresses the interference functionals as an infinite series of higher order derivative terms given by the Laplace transform of the interference power.While the series representation method provides a tractable alternate for handling general fading, it often leads to situations where it is difficult to derive closed form expressions.Numerically evaluating a higher order derivative is also complex and prone to floating-point rounding errors [22].
Aside from the small-scale fading, random shadowing due to obstacles in the local environment or human body movements (in the case of user equipments) can impact link performance by causing fluctuations in the received signal.Shadowing affects the transmission performance, which will be especially pertinent in a dense network or millimeterwave links.Hence, the combined effect of small-scale and shadowed fading needs to be properly incorporated in 5G communications design.In this respect, composite channel models have been proposed in [23]- [27].In [23], the shadowed Nakagami fading distribution was first proposed by combining Nakagami-m multipath fading and lognormal distributed shadowing.Later, [24] introduced the generalized-K model by approximating the shadowing model in [23] using the gamma distribution to improve analytical tractability.Traditional composite channel models (referred to as multiplicative shadow fading models) assume that the shadowing affects the dominant components and the scattered waves equally, whereas, in practice, shadowing often only occurs on the dominant components, which gives rise to a different kind of composite model, often referred to as a LOS shadow fading model.To model shadowing in LOS channels, [25] proposed the Rician shadowed fading model by assuming a Rician distribution for the multipath fading and Nakagami-m distribution for the LOS shadowing.More recently, [26], [27] proposed κ-μ shadowed fading model by assuming κ-μ multipath fading with shadowing of the dominant component.
The κ-μ shadowed fading model is an attractive proposition, not just due to its excellent fit to the fading observed in a range of real-world applications (e.g.device-to-device [27], underwater acoustic [28], body-centric fading channels [29], etc.) but also its extreme versatility.More precisely, it is able to account for most of the popular fading distributions utilized in the literature.Motivated by the comprehensive nature of the κ-μ shadowed fading model, we use it along with a stochastic geometric framework to derive the downlink SINR distribution of a typical user in a K -tier HetNet with K classes of BSs, differing in terms of the transmit power, BS density, shadowing and fading characteristics.We evaluate the average of an arbitrary function of the SINR, which can be easily applied to other network metrics.For instance, it may be utilized to evaluate any performance measure that can be represented as a function of SINR, e.g., the spectral efficiency, outage probability, moments of the SINR, and error probability.
The main contributions of this paper may be summarized as follows.
1) The main difficulty with incorporating generalized fading models into stochastic geometry-based analyses is the lack of tractability in expressing the PDF of the interference.In general, it is more convenient to express the metrics of interest in terms of the Laplace transform of the interference.Nonetheless, this presents significant challenges when extending the analyses from Rayleigh fading to the more general fading models.We overcome this problem by deriving the Laplace transform of the interference over κ-μ shadowed channels in closed form, using which the distribution of interference is explicitly characterized.It is worth highlighting that this model encompasses as special cases, the majority of the fading models proposed in the literature, including Rayleigh, Rician, Nakagami-m, Nakagami-q, One-sided Gaussian, κ-μ, η-μ, and Rician shadowed distribution to name but a few.2) We use tools from stochastic geometry to evaluate the distribution of the SINR, coverage probability and average rate for κ-μ shadowed fading.We also propose a numerically efficient method to calculate the average of an arbitrary function of the SINR.3) We present numerical results which provide useful insights into the performance of cellular networks under different fading conditions.In particular, we observe the trade-off relation between the rate and average SINR based on the channel parameters, such as the intensity of dominant signal components, the number of scattering clusters, and shadowing effect.This information will be of paramount importance to those responsible for designing future 5G network infrastructure to ensure an adequate quality-of-service for all users.This paper is organized as follows.In Section II, the system model and assumptions are introduced.We then apply an orthogonal expansion to the κ-μ shadowed PDF in Section III, characterize the interference distribution in Section IV, and introduce a novel analytical framework in Section V.In Section VI, we present both numerical and simulated results to validate the analysis.Finally, Section VII concludes this paper.

A. Network Model
We consider the downlink of a K -tier HetNet where randomly distributed small-cell BSs, such as pico or femtocell BSs, are overlaid on a network of macrocell BSs.The BSs of each tier may differ in terms of transmit power, spatial density and cell-selection bias.The locations of the k-th tier BSs are modeled by an independent, homogeneous PPP k with density λ k and the union of K point processes constitutes the K -tier HetNet = ∪ k∈K k where K = {1, 2, . . ., K }.The locations of the UEs are modeled by a homogeneous PPP (u) with density λ (u) that is independent of .Orthogonal multiple access is employed at each cell by allocating mutually orthogonal resource blocks to each user equipment (UE), implying no intra-cell interference within a cell.Without loss of generality, we assume that a typical UE is located at the origin and each BS has an infinitely backlogged queue.The received power at a typical UE from a k-th tier BS x k ∈ k2 is given by where a multiplicative channel model H x k = h x k χ x k with large-scale shadowing χ and small-scale fading h is utilized in the second equality, P k is the transmit power of the k-th tier BS, α is the path-loss exponent (α > 2) and τ is the path-loss intercept at a link-length x = 1.

B. Cell Association Policy
We assume a general cell association model where all BSs allow open access and each UE connects to the BS that provides the highest long-term biased received power (LRP) 3without small-scale fading as written below where B j is the bias connecting to the j -th tier BS (B j > 0) and a change of variable, i.e., y = χ − 1 α j x, is applied in the last equality.For a single tier network, (2) is equivalent to connecting with the closest BS.
Due to the displacement theorem [18, Lemma 1], the mapping between x and y converts a PPP j = {x} with density λ j into a new homogeneous PPP (e) j = {y} with density λ (e) j = λ j E[χ δ j ] where δ = 2 α .Thereby, the original network model with large-scale shadowing χ can be equivalently expressed as the network (e) = ∪ j ∈K (e) j without large-scale shadowing where the effect of large-scale shadowing is now incorporated through an appropriate scaling in the density λ j → λ (e) j .Given that the serving BS belongs to the k-th tier, the SINR at a typical UE can be formulated as follows where d = denotes equivalence in distribution, which follows from [18, Lemma 1], x * k represents the location of the associated k-tier BS, \{x * k } denotes the set of interfering BSs, Pj = P j P k represents the ratio between the transmit power of the interfering and serving BSs and N = N P k = N 0 W τ P k is determined by the noise power spectral density N 0 , bandwidth W , transmit power of the associated BS P k , and the reference path-loss τ at a unit distance.Similarly, we denote B j = B j B k as the ratio between the bias factors of the interfering and serving BSs and I j ∈K y∈ (e) j \{y * k } Pj h y y −α as the aggregate interference normalized by the transmit power of the serving BS.Since the cell association policy in (2) is independent of the small-scale fading distribution h, the probability that a typical UE connects to the k-th tier BS, denoted as P k , and the PDF of the link length y * k can be evaluated as below where δ = 2 α and (4) follows directly from [7, Lemma 1] and [18, Lemma 2].

C. Channel Model
Due to the wide range of use cases provisioned for 5G communications, conventional cellular channel models which typically consider only a single source of shadowing (e.g.large-scale shadowing) are unlikely to be sufficient.In reality, it is probable that cellular applications will encounter multiple independent types of shadowing which may or may not occur concurrently.For example in the downlink scenario, the signal transmitted from the BS to the UE will undergo two key types of shadowing, the first of which is large-scale shadowing, denoted here by χ, which is induced due to large terrestrial objects e.g.buildings or hills, which can cause a random fluctuation in the total signal power.In cellular networks, the BSs are usually placed at elevated locations and are typically free from surrounding clutter.However, UEs are most often operated at lower levels and the LOS signal path is often obscured by local obstacles including the user's body itself.Therefore we consider a second type of shadowing which affects (i.e.randomly fluctuates) the dominant signal component.In this contribution, this LOS shadowed smallscale fading is denoted as h and is modeled as a κ-μ shadowed random variable [26], [27].Together, these two independent random processes create an extremely versatile channel model, H = hχ, which can incorporate a wide range of shadowing and fading scenarios.
1) Large-Scale Shadowing: The analysis presented in this paper is valid for any finite distribution of the large-scale shadowing χ and we summarize the three most commonly used large-scale shadowing distributions, namely the lognormal, gamma, and inverse-Gaussian distributions [30], with the , 0 = 10 ln (10) . ( (b) Gamma Shadowing: χ ∼ Gamma(k g , θ g ) with shape parameter k g and scale parameter θ g , (c) Inverse Gaussian Shadowing: χ ∼ I G(μ ig , λ ig ) with mean μ ig and shape parameter λ ig , where K n (z) is a modified Bessel function of the second kind.
2) Small-Scale Fading and LOS Shadowing: The κ-μ shadowed distribution is a very flexible model which contains as special cases the majority of the linear fading models proposed in the open literature, including Rayleigh, Rice (Nakagami-n), Nakagami-m, Hoyt (Nakagami-q), One-Sided Gaussian, κ-μ, η-μ and Rician shadowed to name a few [31] (See Table I).Because of this generality, the κ-μ shadowed fading model can be used to account for small-scale fading which originates due to LOS or non-LOS conditions, multipath clustering with circularly symmetric or elliptical scattering, and power imbalance between the in-phase and quadrature signal components.
The channel coefficient h of a κ-μ shadowed fading channel can be expressed in terms of the in-phase and quadrature components of the fading signal as follows where μ is the number of the multipath clusters, 4 ξ is a Nakagami-m distributed random variable with E ξ 2 = 1, X i and Y i are mutually independent Gaussian random variables with p i and q i are real numbers and d 2 = μ i=1 p 2 i + q 2 i is the power of the dominant components.
In the following, we summarize the key statistics of the κ-μ shadowed fading model which will be used in the network performance analysis conducted here.
Lemma 1: The PDF, fractional moment, and Laplace transform of h for a κ-μ shadowed channel are respectively given by h μ (1+κ)m , κ, μ, m and l are positive real-valued constants, (t) is the gamma function and 1 F 1 a b x is the confluent hypergeometric function.Proof: A detailed derivation of the PDF and Laplace transform expression is provided in [26].The fractional moment of h can be derived as follows where the PDF f h (x) from ( 10) is substituted in the first equality, a change of variables, i.e., t ← θ 2 −θ 1 θ 1 θ 2 x, and (48) is applied to the last equality.
Physically, κ = d 2 2μσ 2 represents the ratio between the total power of the dominant components and the total power of the scattered waves, μ denotes the real-valued extension of the number of multipath clusters, and m indicates the amount of shadowed perturbation in the dominant component as illustrated in Fig. 1 (a).Since the Laplace transform of the Nakagami-m distribution converges to lim m→∞ L h (s) = lim m→∞ (1 + s h/m) −m = e −s h , the dominant component becomes increasingly deterministic as m → ∞.Hence, a κ-μ shadowed fading channel with m → ∞ has a constant dominant power and is therefore equivalent to a κ-μ faded channel.
3) Combined Large-Scale Shadowing, Small-Scale Fading and LOS Shadowing: Since the κ-μ shadowed fading model includes small-scale fading and LOS shadowed fading as special cases, the proposed channel model H = hχ can be used to represent four different classes of fading environments as illustrated in Fig. 1 (b); namely 1) small-scale fading only if χ is constant, 2) small-scale fading with LOS shadowed fading only if h is either Rician shadowed or κ-μ shadowed and χ is constant, 3) traditional composite fading/shadowing if h is the result of small-scale fading only with randomly distributed χ, and 4) double shadowed fading conditions if h is the result of small-scale and LOS shadowed fading and χ is a random variable.
Remark 1 (Multiple Antenna Systems): If the BS is equipped with N c antennas and communicating with a single-antenna UE with zero-force beamforming, the corresponding channel between the BS and a UE can be represented as a summation of N c i.i.d.κ-μ shadowed random variables, which is a κ-μ shadowed random variable with fading parameters (κ, N c μ, N c m) [26], [31].Thereby, the theoretical analysis presented in this paper can be directly applied to multiple-antenna diversity systems.

III. LAGUERRE POLYNOMIAL SERIES EXPANSION OF THE κ -μ SHADOWED DISTRIBUTION
As we can see from (10), the κ-μ shadowed distribution includes the hypergeometric function which often leads to a computationally complex performance evaluation.Due to the inherent mathematical intractability, limited work has been conducted which considers κ-μ shadowed fading in the context of stochastic geometry.Most notably, in [32], the author approximated a κ-μ shadowed random variable using a gamma distributed random variable based on second-order moment matching, but the accuracy of this approximation can not be guaranteed for all fading parameters.In [33], the authors analyzed a cellular network over κ-μ shadowed fading where they represented the confluent hypergeometric function by its Although the series representation converges locally, it is valid only for integer-valued parameters a and b, the radius of convergence diverges over different combinations of parameters, and it is computationally complex to evaluate.As illustrated in Fig. 2, there are noticeable discrepancies between the approximation methods proposed in [32] and [33] and the exact PDF for several cases, limiting their application. 5o overcome this problem, we adopt the generalized Laguerre polynomial expansion proposed in [34] and [35] that is analogous to the Fourier series: As a Fourier series can represent any PDF in terms of harmonic bases, we use a generalized Laguerre polynomial as an orthogonal base and simplify the PDF and CDF of the κ-μ shadowed fading model as given below.
Lemma 2: The PDF and CDF of the channel coefficient h for the κ-μ shadowed fading model can be expressed in series expression form as follows Fig. 2. Numerical evaluation of the κ-μ shadowed fading distribution using the PDFs given in ( 10), ( 12), [32], and [33] where κ, μ and m are positive real-valued parameters, L is the generalized Laguerre polynomial of degree n and order μ − 1 at x, 0 ≤ x < ∞, γ (μ, x) is the lower incomplete gamma function, the coefficients C n , c i,n , and b i,n are calculated as written below and E h j can be obtained from (11) by substituting nonnegative integer index j to the index l.Proof: See Appendix II.Remark 2: If μ and m are positive integers, then by using [36,Th. 1], the expression in (12) can be simplified to a single summation with finite terms as follows.
where A 1 j , A 2 j , B j are given in [36, eq. (6)].(12) and (15) imply that κ-μ shadowed fading is the result of a linear combination of gamma distributed random variables, which follows a gamma mixture distribution.To represent κ-μ shadowed fading as a gamma mixture model, a double summation with infinite terms is required for real valued μ and m, whereas for integer valued μ and m, only a single summation with finite terms is necessary.

IV. DISTRIBUTION OF THE AGGREGATE INTERFERENCE
In this section, we calculate the Laplace transform of the aggregate interference for the κ-μ shadowed fading channel and characterize the distribution of the interference.As shall be seen in Section V, the Laplace transform of the aggregate interference is a crucial measure for evaluating network performance in stochastic geometry based analyses.
Lemma 3: Given that a typical UE is associated to the BS y * k located at y * k = r (or equivalently expressed as x * k using x = χ 1 α j y), the Laplace transform of the aggregate interference over a multiplicative channel with κ-μ shadowed fading and large-scale shadowing is calculated as where the subindex j ∈ K = {1, 2, . . ., K } represents the parameters for the j -th tier, Pj = is given by ( 5)- (7), W j (z) denotes the following expression is the Appell Hypergeometric function which is defined in (53), Appendix I [37].
Proof: See Appendix III.Some comments on the numerical computation of the Appell's function F 2 (•) are presented in Appendix IV.
By using a change of variable, i.e., sz −1 = r α , the Laplace transform of the interference can be expressed as which indicates that the aggregate interference follows a Stable distribution as described below.Note that the exclusion zone in the interference field is considered in (16) based on the condition y * k = r .Lemma 4: The aggregate interference over a multiplicative channel of κ-μ shadowed fading and large-scale shadowing follows a Stable distribution [3] with four parameters; namely, stability δ, (17).The fractional moment of the aggregate interference is given by for 0 < l < 2 α .Any moment with order above l > 2 α is undefined, i.e., becomes infinity.
Using (54), the Appell's function reduces to a Gauss hypergeometric function if one of the parameters is zero.Hence the expressions in ( 16) and ( 17) can be simplified as below.

V. THEORETICAL ANALYSIS OF THE PERFORMANCE MEASURES
In this section, we propose a novel method to compute E g (γ ) for an arbitrary function of the SINR g(γ ) using stochastic geometry.The original idea was proposed by Hamdi in [38] for Nakagami-m fading, and later in [39] for κ-μ and η-μ fading, which we further extend to κ-μ shadowed fading in this paper.By using the proposed method, one can evaluate any performance measures that are represented as a function of SINR (or SIR).For instance, the spectral efficiency, outage probability, moments of the SINR, and error probability can be expressed as averages of g(x) = log(1 + x), g(x) = I(x ≤ x 0 ), g(x) = x n , and g(x) = Q (x), respectively.

A. General Case and Main Result
Theorem 1: For a K -tier HetNet with κ-μ shadowed fading, E [g (SINR)] is given by where P k is derived in ( 4), SINR k represents the SINR when a typical UE is associated to the k-th tier BS y * k , C n is defined in (14), and ξ i represents the following integral the distribution f y * k (r ) is given by ( 4) and L I (s) is derived in (16).g μ+i (z) is defined as where we used the general Leibniz rule in the last equality.Proof: See Appendix VI.Theorem 1 is the most general result in this paper that evaluates arbitrary performance measures for a K -tier Het-Net, considering noise, interference, per-tier BS density, and different fading and shadowing characteristics across each tier.The analytic functions g(z) and g μ+i (z) for various performance measures are summarized in Table II. 6We also note that E [g (SINR k )] in (24) is computationally efficient; the computational complexity of ( 24) is same as a single summation expression since ξ i is independent of the index n.
Remark 3: If μ and m are positive integers, (15) can be utilized to achieve an expression analogous to Theorem 1, in terms of a single summation with finite terms as described where θ z dz and the coefficients A 1 j , A 2 j , B j are derived in [36, eq. (6)].The proof of ( 27) is omitted since it is analogous to Theorem 1.

B. Special Cases
Theorem 1 and E r e −r α N z L I (r α z) can be further simplified for a number of special conditions, such as a noise-limited scenario, an interference-limited scenario, or identical fading and shadowing parameters for all tiers, which are described below.
Lemma 6 (Fixed Path-Loss): If α = 4, the term E r e −r 4 N z L I r 4 z in (25) can be simplified as where denotes the following expression and W j (z) is defined in (17).If all tiers have identical fading parameters (κ, μ, m), then W (z) = W j (z) for any j ∈ K and ( 28) can be further simplified as the term E r e −r α N z L I (r α z) in ( 25) can be simplified as If all tiers have identical fading characteristics, then (30) reduces to a succinct form as

Lemma 8 (Noise-Limited Scenario With Fixed Path-Loss): If I
N and α = 4, then E r e −r 4 N z L I r 4 z can be simplified as Proof: The proof of Lemmas 6, 7 and 8 are given in Appendix VII.
Remark 4: If all tiers have identical fading characteristics and are interference-limited, the performance measure E [g (SINR)] can be expressed by using Lemma 5 as follows where C n , g μ+i (z) and W (z) are independent of the PPP density λ j .(33) provides an important insight into the system performance of a PPP-distributed cellular network with κ-μ shadowed fading and arbitrary large-scale shadowing.Specifically, any performance measure of a PPP-distributed HetNet that can be represented as a function of SINR, is independent of the BS transmit power P k , BS density λ k , and the number of tiers K .This invariance property was originally introduced in [6], [7], and [9] for Rayleigh fading.(33) generalizes this argument by proving that the invariance property holds for any linear small-scale fading and finite large-scale shadowing distribution.
Next, we apply Theorem 1, Lemmas 5, 6 and 7 to evaluate various performance measures.

C. Performance Measure 1: Spectral Efficiency
Spectral efficiency is defined in [7] as where P k is the tier association probability to the k-th tier evaluated by ( 3) and SINR k is the received SINR from the k-th tier BS.E [ln (1 + SINR k )] can be evaluated by using Theorem 1 with g(z) = ln(1 + z) and g μ+i (z) as follows [38] Given identical channel characteristics across each tier, the spectral efficiency reduces to where C n is given by ( 14), W (z) is derived in (17), is defined in (28) and K (z) denotes Remark 5: Theorem 1 can be applied to Rayleigh fading by letting μ = 1, C 0 = 1 and C n = 0 for n > 0 in (24).Then (35) can be further simplified as follows by using a change of variable, i.e., x = e t − 1.As expected, the above expression is identical to [7, eq. ( 27)].

D. Performance Measure 2: Moments of the SINR
Higher order moments of the SINR are crucial performance measures which have an important role in the characterization of network performance.E SINR l can be evaluated by using Theorem 1 with g(z) = z l and g μ+i (z) as For the case when we have identical channel characteristics across each tier, the fractional moment is simplified to where (x) n = (x+n) (x) is the Pochhammer symbol, the index l is a positive real-valued constant and K (z) is defined in (36).

E. Performance Measure 3: Outage Probability
The outage probability is defined in [6] as (40) for a predefined threshold T .Theoretically, one can use Theorem 1 to calculate (40) by approximating the step function with a smooth sigmoid function, i.e., g(z , where controls the sharpness.However, even with a smooth function, g μ+i (z) behaves like an impulse signal for a large derivation order μ + i [40].Hence, most numerical software will present a precision overflow while evaluating (25).
Instead of using Theorem 1, it appears more convenient to use the Gil-Pelaez's inversion based approach [41] as follows.
Step 1) P (SINR k < T ) in ( 40) can be represented in terms of the interference distribution as follows where the expectation in (41) averages over the link length r and the channel coefficient h.
Step 2) The CDF of the interference can be derived using the Gil-Pelaez's inversion as follows where I m(z) represents the imaginary part of a complex number z.

VI. NUMERICAL RESULTS
In this section, we present numerical evaluations of the theoretical results and compare them with Monte-Carlo simulations.All of the numerical results presented in this paper were obtained by using the Julia language which provides fast computation times and a straight-forward syntax that is similar to Matlab [42].The analytical results are plotted as lined curves without markers, whereas the simulation results are represented by markers without a line.Both Figs. 3 and 4 show that the numerical results accurately match the simulation results in every scenario.In our analysis, we considered a two-tier HetNet with BS intensity λ 1 = 2λ 2 , transmit power P 2 = P 1 − 20 dBm, path-loss exponent α = 4 and lognormal distributed χ with (μ l , σ l ) = (0, 4).We assumed identical fading and shadowing parameters for both tiers.
Fig. 3 compares the spectral efficiency across a wide range of channel parameters within an interference-limited environment.We observed the following patterns in the spectral efficiency: 1) Strong dominant components or weak scattered components (large κ) achieve a higher rate.2) Rich scattering with a large number of multipath clusters (large μ) achieves a higher rate.3) Strong LOS shadowing (small m) achieves a higher rate if μ ≤ 1 and vice versa.Large κ and μ parameters indicate strong LOS link and rich multipath components, respectively, both collectively contribute toward increasing the average rate, as illustrated in Figs 3 (a)-(b).On the contrary, the average rate does not show a monotonic behavior for changes in the m parameter.A small value of m implies strong random fluctuation of the dominant component, which decreases not only the received signal power but also the interference level.is weak (small κ) or the number of multipath clusters is small (small μ), LOS shadowing subsides as m increases, which increases the interference power, thus deteriorating the received SINR as well as the average rate.Figs 4 (a)-(b) compare the spectral efficiency versus the macro BS intensity λ 1 .As mentioned in Remark 4, the spectral efficiency becomes invariant for a large BS intensity λ 1 .In a dense network with a large BS intensity, the aggregate interference becomes significantly larger than the noise power, achieving an interference-limited condition.Additionally we observe that the BS intensity required to reach the rate asymptote is inversely proportional to the operating SNR level.For a high SNR regime, the average rate reaches the asymptote around λ 1 = 10 −2 , whereas in a low SNR regime, a large number of BSs (λ 1 ≥ 10 −1 ) are required to obtain sufficiently larger interference power than the noise.In a sparse network with small BS intensity λ 1 , different fading parameters, such as m, do not affect the spectral efficiency significantly.However, in a dense network with a large BS intensity λ 1 , different fading parameters have a notable effect on the rate, as illustrated in Figs 4 (a)-(b).This pattern is due to (z) in (29), which is proportional to λ 0 (1 + W (z)) and W (z) is a function of the channel parameters.If λ 0 (or λ j ) is large, any difference in the channel parameters will be emphasized, affecting the spectral efficiency via (z).However, if λ 0 is small, any difference in the channel parameters will be unnoticed due to λ 0 (1 + W (z)), thus achieving nearly identical network performance.Figs 4 (c)-(d) compare the average SINR across a range of channel parameters assuming an interference-limited environment.We observe that large κ and μ parameters jointly achieve a higher rate similar to Fig. 3.We also notice that for large κ (κ ≥ 15) and weak shadowing (large m), a higher SINR level is achieved.

VII. CONCLUSION
In this paper, we have considered a cellular network in which the signal fluctuation is the result of large-scale and LOS shadowing to encapsulate the diverse range of channel conditions that can occur in 5G communications.We applied a Laguerre polynomial series expansion to represent the κ-μ shadowed fading distribution as a simplified series expression.Based on the series expressions, we then proposed a novel stochastic geometric method to evaluate the average of an arbitrary function of the SINR over κ-μ shadowed fading channels.The proposed method is numerically efficient, can be easily applied to other network models, and can evaluate any performance measure that can be represented as a function of SINR.Using the proposed method, we have evaluated the spectral efficiency, moments of the SINR, bit error probability and outage probability of a K -tier HetNet with K classes of BSs, differing in terms of the transmit power, BS density, shadowing characteristics and small-scale fading.Furthermore, we provided numerical results and investigated the performance over a range of channel parameters and observed that a dominant LOS component (large κ), rich scattering environment (large μ) and weak shadowing condition (large m) collectively lead to high spectral efficiency.Finally, it is worth remarking that the analytical framework proposed in this paper can be applied to practical use cases of 5G communications, where Rayleigh fading fails to fully capture the diverse nature of the underlying channel.

APPENDIX I
In this appendix, we summarize the operational equalities of the special functions, which are used in this paper. 7First, the generalized Laguerre polynomial of degree n and order β has the following functional identities The following properties of the hypergeometric function hold for real constants a, b and c (50) The lower incomplete gamma function γ (s, x) = x 0 t s−1 e −t dt has the following series representation and functional identity for arbitrary positive real constant s x s+n e −x (s The binomial coefficient can be defined for real constants x, y using the gamma function as Appell's function F 2 (•) is defined via the Pochhammer symbol (x) n = (x+n) (x) as follows (53) 7 Most of the expressions in Appendix I were introduced in [43], except for (53) and (54), which were proved in [37].
Appell's function can be reduced to the hypergeometric function using the following properties The following integration holds under the following constraints d > 0 and Gauss-Laguerre quadratures can be used to evaluate the following integral for a given analytic function g(x) as where x n and w n are the n-th abscissa and weight of the N-th order Laguerre polynomial.

APPENDIX II
In this appendix, we provide a proof of Lemma 2. The PDF of h for κ-μ shadowed fading in (10) can be represented in the orthogonal series expansion form as where we applied the Laguerre polynomial series expansion in [34, eq. ( 9)] and the coefficient C n is evaluated by substituting (10) as follows [34, eq. ( 8)] The integral I 1 can be simplified by using the series representation of L μ−1 n (x) in (45) as follows where we used (10) to express the integral as the PDF of the κ-μ shadowed fading in the last equality.Then, by substituting (59) into (58), the coefficient C n in ( 14) can be derived after algebraic manipulation.The series expansion form in (57) can be further simplified by using (45) and (52) as follows which achieves (12).
The CDF of h can be evaluated as follows where we used (57) in the second equality, utilized (46) in the third equality, applied a change of variable, i.e., n and (45) in the last equality.The coefficient b i,n can be simplified by using (52) as then the CDF in ( 13) can be subsequently obtained.This completes the proof.

APPENDIX III
In this appendix, we provide a proof of Lemma 3. Due to (2), all interfering BSs within the j -th tier are located farther The Laplace transform of the interference from the j -th tier is derived in (64), as shown at the bottom of this page, where we represented the distance to the serving BS as y * k = r in the second equality, applied a change of variable, i.e., s Pj B j hl −α = t, in the third equality, then used integration by parts.The first part of the expectation term in (64) is evaluated by (65), as shown at the bottom of this page, where we used the PDF of κ-μ shadowed fading with a change of variable, i.e., h θ 1 = t in the first equality, applied (51) to the second equality, utilized the integration (55) in the last equality [37], A = 1−θ 1 /θ 2 1+θ 1 sr −α and B = θ 1 sr −α 1+θ 1 sr −α .The second part of the expectation term in (64) follows directly by using the Laplace transform of the κ-μ shadowed channel coefficient (10).By denoting sr −α = z, ( 16) and ( 17) can be obtained.This completes the proof.

APPENDIX IV
All numerical results provided in this paper were obtained by using the Julia language [42] where we used the appellf2 function in the SymPy package [44] to evaluate the Appell's function in (17).However, if one needs to use MATLAB, where a native Appell's function library does not exist yet, W j (z) can be accurately approximated by using the Gauss-Laguerre Quadrature as (66), as shown at the bottom of this page, where f x n and w n are the n-th abscissa and weight of the Nth order Laguerre polynomial.Since the approximation error converges rapidly to zero [43], (66) provides a numerically accurate and efficient approximation to W j (z).

APPENDIX VI
In this appendix, we provide a proof of Theorem 1.The average of an arbitrary function of the SINR where we denoted b = r α (I + N), applied integration by parts, defined g k (z) in ( 26), and g k (0) = 0, for k < μ + i − 1 g(0), for k = μ + i − 1. (73) Then, the average of an arbitrary function of the SINR is given by (70), as shown at the top of this page, where we used n i=0 (−1) i n i = 0 in the second equality.This completes the proof.

APPENDIX VII
In this appendix, we provide proof of Lemmas 6, 7 and 8.By substituting (4) and ( 16) in (25), the expectation term E r e −r α N z L I (r α z) can be evaluated as (71), as shown at the top of the previous page, where we used a change of variable, i.e., t = πr 2 in the last equality.If α = 4, then (28) can be achieved by applying (50).Given an interference-limited condition, (71) reduces to (75) (32) readily follows by substituting W j (z) → 0 and α = 4 in (28).This completes the proof.
Based on Figs 3 (c)-(d), it can be observed that the average rate is an increasing function of m if and only if μ > 1 and κ > 16, where the actual threshold values are determined by the network configuration.Otherwise, the rate is a decreasing function of m.If either the dominant LOS component
k where y * k denotes the associated k-th tier BS and Pj = P j P k is the transmit power ratio between the interfering and serving BS B j P j y −α < B k P k y * k −α for any y ∈

h x 0 x 0 g
−αI +N is evaluated in (69), as shown at the top of this page, where(12) is used in the second equality, a change of variable, i.e., xr −α I +N = z, is utilized in the third equality, and (14) is employed in the last equality.(69) can be evaluated as follows μ+i (z)e −bz dz,

TABLE I SPECIAL
CASES OF THE κ -μ SHADOWED FADING MODEL fractional moment for each distributions summarized below, where l is a positive real number.(a)Lognormal Shadowing: χ ∼ L N (μ ln , σ 2 ln ) with mean μ ln and standard deviation σ ln ,

TABLE II THEORETICAL
FRAMEWORK FOR EVALUATING VARIOUS SYSTEM MEASURES