Analysis of the Delay Distribution in Cellular Networks by Using Stochastic Geometry

In this paper, with the aid of the mathematical tool of stochastic geometry, we introduce analytical and computational frameworks for the distribution of three different definitions of delay, i.e., the time that it takes for a user to successfully receive a data packet, in large-scale cellular networks. We also provide an asymptotic analysis of one of the delay distributions, which can be regarded as the packet loss probability of a given network. To mitigate the inherent computational difficulties of the obtained analytical formulations in some cases, we propose efficient numerical approximations based on the numerical inversion method, the Riemann sum, and the Beta distribution. Finally, we demonstrate the accuracy of the obtained analytical formulations and the corresponding approximations against Monte Carlo simulation results, and unveil insights on the delay performance with respect to several design parameters, such as the decoding threshold, the transmit power, and the deployment density of the base stations. The proposed methods can facilitate the analysis and optimization of cellular networks subject to reliability constraints on the network packet delay that are not restricted to the local (average) delay, e.g., in the context of delay sensitive applications.


A. Motivation and Related Works
The fifth generation (5G) of wireless networks introduces ultra reliable low latency communication (URLLC) as one of its main use cases which promises to provide low latency and ultrahigh reliability for mission critical applications, such as the industrial Internet, the smart grid, remote surgery, and intelligent transportation systems.Motivated by this, delay and reliability emerge as important performance metrics to consider, in addition to traditional key performance indicators, such as the capacity and the spectral efficiency, in communication systems design.
According to the 3rd Generation Partnership Project (3GPP) [1], the end-to-end delay of 5G communications is required to be less than 1 ms, which is about 1/200 of the delay requirement for the fourth generation (4G), and the target reliability must be as high as 1 − 10 −2 to 1 − 10 −7 .
In cellular networks, the delay is tightly related to the transmission success probability (i.e., the coverage probability) which is typically determined by the signal-to-interference ratio (SIR) or the signal-to-noise-plus-interference ratio (SINR) at the user side.Recently, to take into account the strong interplay between the transmit power and the density of the base stations (BSs) for optimal cellular networks planning, a new definition of the coverage was introduced by the authors of [2], which depends on the SIR and the average signal-to-noise ratio (ASNR).Since the SIR and SINR are functions of the received signal as well as the interference from the BSs in the network, they strongly depend on the inter-distances among the BSs and the users.However, characterizing their exact expressions is challenging due to the complex nature of the deployment of BSs and users locations.Leveraging the theory of stochastic geometry, fundamental works such as [3] and [4] show that many performance metrics, such as the coverage and rate can be derived in a tractable way by modeling the spatial distribution of transmitters and receivers as homogeneous Poisson point processes (PPPs).Thanks to its simplicity and tractability, this modeling approach has attracted a lot of interest and, as a result, stochastic geometry has emerged as a major tool for modeling the spatial distribution of wireless networks.In 5G and beyond networks, the role of stochastic geometry is envisioned to be more prominent due to the increased heterogeneity and complexity of future wireless networks [5].
A fundamental measure for the delay in a large-scale wireless network is the local delay, which is defined as the average time, measured in the number of time slots, that is required for a successful packet transmission between a typical transmitter-receiver pair.With the aid of stochastic geometry, the authors of [6]- [8] characterize the local delay in a mobile ad-hoc network under different transmitter-receiver association criteria, such as bipolar, nearest-neighbor, nearesttransmitter, and nearest-receiver models.The local delay in heterogeneous cellular networks is derived in [9], [10].The authors of [9], in particular, compute the local delay and energy efficiency under a finite local delay constraint in a downlink heterogeneous cellular network with a discontinuous transmission (DTX) scheme, where a BS is randomly turned off in a given time slot to reduce the interference and energy consumption.In [10], on the other hand, the authors study the trade-off between the delay and reliability by taking into account the random BS activity.It is also worth mentioning that the local delay can be characterized as the −1-th moment generating function of the SIR/SINR meta distribution.For Poisson cellular networks and SIR-based coverage, it is given by the inverse of a hypergeometric function [11].The meta distribution of the SINR-based coverage criterion in cellular networks is derived in [12]- [14].
Despite the relative ease of analysis and tractability, the local delay has some inherent limitations.First of all, under some circumstances, a network may exhibit a phenomena called wireless contention phase transition [6], according to which the local delay is infinite due to the contributions from a few users located in "bad spots".Therefore, the evaluation of the performance based on the local delay alone renders networks with an infinite delay unreliable, despite the fact that these networks may provide good service in terms of delay for most of their users.In 5G and beyond networks, furthermore, mission-critical applications such as vehicular communications and remote surgery are envisioned to be among the most important use cases in which, rather than the average delay, more emphasis is put on the delay distributions and percentiles.Unfortunately, these fine-grained information cannot be provided by the traditional definition of local delay, since it provides information only about the average delay.

B. Contributions
Motivated by the limitations of the local delay and the need for more relevant metrics for characterizing the delay for the applications envisioned in 5G and beyond communications, we focus our attention of the distribution of the delay rather than on the average delay.More precisely, we introduce and compute analytical frameworks for three definitions of delay distribution.Each definition of delay distribution uniquely captures fine-grained delay information for the users.
To overcome some computational challenges of the obtained analytical formulations, we also propose efficient numerical approximations for computing the delay distributions in practice.To the best of our knowledge, this work offers the first comprehensive analytical formulation of the delay distribution in large-scale cellular networks by using stochastic geometry.
In summary, the main contributions of our paper are as follows: • We introduce a framework for analyzing a cellular network modeled as a homogeneous PPP where the active-inactive status of each BS depends on the existence of a user in the cell of the BS.We introduce different coverage criteria and compute the local delay for each criterion.In this paper, we consider three coverage criteria: 1) SIR-based coverage: Under this criterion, a transmission is successful if the SIR at the user is greater than a given decoding threshold.
2) SINR-based coverage: Under this criterion, a transmission is successful if the SINR at the user is greater than a given decoding threshold.
3) SIR+ASNR-based coverage: Under this criterion, a transmission is successful if the SIR at the user is greater than a given decoding threshold and the ASNR is greater than a given detection threshold.
• We consider three different delay distribution metrics and derive analytical expressions for each distribution and for the SIR-based, SINR-based, and SINR+ASNR-based definitions of coverage.In particular, the following delay distribution metrics are analyzed: 1) The first distribution is formulated as the expectation, over all network spatial realizations, of the complementary cumulative distribution function (CCDF) of the required number of time slots for successful transmission over all possible channel realizations.
2) The second distribution is formulated as the CCDF of the expectation of the required number of time slots for successful transmission conditioned on a fixed network realization.
3) The third distribution is formulated as the CCDF, of the CCDF conditioned on the network realization, of the required number of time slots for successful transmission.
• As for the first distribution, we analyze the asymptotic packet loss probability.We focus our attention on the event that a typical user never receives its packet regardless of the waiting time.We show that the packet loss probability is zero when the SIR-based or the SINR-based coverage is considered, and it is non-zero when the SIR+ASNR-based coverage is considered.
• To overcome some numerical difficulties that arise in the computation of the obtained delay distribution, we introduce efficient numerical approximations.In particular, we show that the second and the third delay distributions can be approximated by using the numerical inversion of the Laplace transform and (when the SIR-based or the SINR-based coverage are considered) the Beta distribution.We also show that the first delay distribution can be expressed in terms of an integration of the third distribution, and thus it can be efficiently approximated by using the Riemann sum.

C. Paper Organization
The rest of this paper is organized as follows.In Section II, we introduce the system model.
In Section III, we give preliminary results and definitions for the three delay distributions.In Expectation operator, probability measure Expectation over network realization, expectation over fading realization Probability over network realization, probability over fading L(λBS, λMT) = 1 − (1 + λMT/(3.5λBS))−3.5  Active probability of a base station P , S, I, W Transmit power, signal, interference, noise γ, θ decoding threshold, detection threshold ∆Φ, P Φ cov Conditional number of time slots, conditional coverage probability Section IV, we provide the analysis of the local delay under different coverage criteria and show that, under some conditions, it is finite under the SIR-based coverage, while it is infinite for the other coverage criteria.In Sections V and VI, we introduce analytical expressions of the delay distributions for different coverage criteria.In particular, in Section V, we provide the asymptotic analysis of the first distribution which corresponds to the packet loss probability.In Section VII, we propose efficient numerical approximations for the delay distributions based on techniques such as the Euler-sum method, the Beta distribution, and the Riemann sum, so as to overcome some inherent numerical difficulties in the evaluation of the delay distributions.
Finally, we provide numerical results in Section VIII to validate our findings, and we conclude the paper in Section IX.The main notation used in this paper is listed in Table I.

A. Cellular Network Modeling
We consider a downlink cellular network where the BSs are modeled as points of a homogeneous PPP, denoted by Φ BS , with density λ BS while the mobile terminals (MTs) are modeled as another homogeneous PPP, denoted by Φ MT , with density λ MT .We assume that Φ BS and Φ MT are independent of each other.Thus, the whole network can be modeled as a point process Φ which is the joint set of BSs and MTs, i.e., Φ = Φ BS ∪ Φ MT .Due to the properties of homogeneous PPPs, Φ forms another homogeneous PPP with the intensity λ = λ BS + λ MT .Each BS transmits with constant power denoted by P .Leveraging Slivnyak's theorem [15,Th. 1.4.5],we focus our analysis on the performance as seen by a typical user located at the origin, denoted by MT 0 , whose serving BS is denoted by BS 0 .Throughout the paper, the subscripts 0, i and n identify the intended link, a generic interfering link, and a generic BS-to-MT link.It is worth mentioning that, although the network is modeled as a random variable, it is static.This means that, for a given network realization, the locations of the MTs and BSs do not change over time.

B. Transmission Scheme
We consider that time is divided into discrete slots with equal duration and each transmission attempt occupies one time slot.At the beginning of a time slot, a typical user, MT 0 , requests a packet to its serving BS, BS 0 .Depending on the medium access control (MAC) policy, the BS 0 may or may not handle the request in the same time slot where the request is made.In this paper, we assume that BS 0 will attempt a transmission of the packet requested by MT 0 immediately.If the attempt fails due to an insufficient received signal quality, a re-transmission is scheduled in the next time slot until the user gets the requested packet.

C. Channel Modeling
For each BS-to-MT link, we model the wireless channel with a distance-dependent pathloss and fast-fading that takes into account the small-scale fading along the transmission path.
Additionally, we consider that all BS-to-MT links are mutually independent and identically distributed (i.i.d.).We assume that the fading is quasi-static, that is, the fading coefficient is constant over a period of a time slot and it varies from one time slot to another according to the i.i.d.Rayleigh fading model.The power fading coefficient between MT 0 and its serving BS, BS 0 is denoted by h 0 , while the coefficient between MT 0 and an interfering base station, The path-loss from a BS and an MT is modeled as where K > 0 is the path-loss constant, r > 0 is the distance between the BS and the MT, and α > 2 is the path-loss exponent.

D. Cell Association Criterion
The cell association criterion is based on the highest average received power without considering the impact of fading.Thus, the association depends only on the deterministic path loss within a fixed network.Let BS n ∈ Φ BS denote a generic BS of the network.Thus, given a typical user, MT 0 , its serving BS, BS 0 , is identified as follows: Let us denote the distance from MT 0 to BS 0 by r 0 .Due to the randomness of the network, r 0 is a random variable whose probability density function (PDF) is given as follows: where indicates that the only random variable in evaluating the probability is the network geometry Φ.In non-fully loaded scenarios, there exist some BSs that are not associated with any MT in the network.In such cases, the BSs are said to be inactive.On the other hand, if a BS is associated with at least one MT in its cell, it is said to be active.Given λ BS and λ MT , the probability of a randomly selected BS being active, denoted by L(λ BS , λ MT ), is [16]: Due to the static assumption of the network, for a given network realization Φ, the active/inactive statuses of all BSs remain the same all the time.Therefore, the time-dependence of the interference is a function of only the power fading coefficients h i .Let Φ I be the set of all interfering BSs, i.e. all active BSs other than BS 0 .Since Φ I is constructed from a homogeneous PPP with intensity λ BS , the active probability of a generic BS is L(λ BS , λ MT ), and a single point at the origin is removed, it forms an inhomogeneous point process [17] with density λ * (r) defined as:

III. PRELIMINARIES
In this section, we introduce necessary concepts for further analysis, such as the conditional coverage probability (CCP), the local delay, as well as the definitions of three delay distributions.

A. Conditional Coverage Probability
Consider the typical user, MT 0 , located at the origin.The signal from the serving BS, BS 0 , and the aggregated signals from the interfering BSs, ∪ BS i = Φ I , are denoted by S and I, respectively, and they are given as follows: Using these notations, the SIR, SINR, and ASNR at MT 0 are given, respectively, as follows: where W > 0 is a deterministic real number which represents the variance of the Gaussian noise power density.Note that SIR, SINR, and ASNR are random variables.In particular, the randomness of SIR and SINR comes from the network geometry Φ and the fading coefficients h 0 and h i , while the randomness of ASNR comes only from Φ.
In general, the CCP is defined as the probability of the event that MT 0 is successfully served by BS 0 in a given time slot, with the (spatial) realization of Φ being fixed.The specific expression of the CCP depends on the criterion used to define a successful transmission (i.e., the coverage probability).Throughout this paper, we consider three different coverage criteria: (i) SIR-based, (ii) SINR-based, and (iii) SIR+ASNR-based.To this end, let P Φ cov denote the CCP for a given realization of Φ and let P h [•] and E h [•] be the probability and expectation operator that treat the fading coefficients h 0 and h i as the only random variables (while the other variables are fixed).
Using these notations, the definition of the CCP for each criterion is given as follows.

1) SIR-based coverage:
Under this criterion, the MT 0 is said to be covered by BS 0 if, given Φ, the SIR is greater than a given decoding threshold γ.Thus, the CCP is formulated as: 2) SINR-based coverage: Under this criterion, the MT 0 is said to be covered by BS 0 if, given Φ, the SNIR is greater than a given decoding threshold γ.Thus, the CCP is formulated as: If W = 0 in (7), the SINR is equal to the SIR, and thus (9) simplifies to (8).
3) SIR+ASNR-based coverage: Under this criterion, the MT 0 is said to be covered by BS 0 if, given Φ, the SIR is greater than a given decoding threshold γ and the ASNR is greater than a given detection threshold θ.Thus, the CCP is formulated as: If W = 0 in (7), the ASNR becomes infinite.Thus, the second condition in (10) always holds and the definition of the CCP coincides with (8).

B. Local Delay
Given Φ, let ∆ Φ be the number of time slots required for a successful transmission, i.e., until MT 0 is covered by BS 0 and can successfully decode the data packet intended to it.Due to the time-independence of the channels between MT 0 and all the BSs in the network, the success probability for any time slot is equal to P Φ cov .In this case, the packet transmission from BS 0 to MT 0 can be regarded as a Bernoulli trial with success probability P Φ cov .Therefore, ∆ Φ is a geometrically distributed random variable whose probability mass function (PMF) and expected value are given as follows: The local delay, denoted by D, is defined as the expected number of time slots required for successful transmission.Thus, it can be obtained by spatially averaging ∆ Φ , i.e.,

C. Delay Distribution
As discussed in Section I, the local delay may have some limitations for characterizing the delay reliability of cellular networks, e.g., it provides information only on the average delay.To overcome these limitations, we propose three delay distribution metrics, where each metric captures some unique statistical information.The mathematical definitions of the delay distributions and their relations to relevant performance metrics in wireless communications are elaborated as follows.
1) F 1 distribution: Given τ ∈ Z + , the F 1 distribution is defined as follows: In the context of wireless communications, F 1 (τ ) is the probability that a typical user experiences a delay greater than the predetermined time deadline τ .Therefore, it can be regarded as some kind of delay violation probability.Equivalently, it can also be regarded as the expected fraction of users in the network that experience a delay greater than τ time slots.It is worth mentioning that the local delay can be expressed in terms of the F 1 distribution.
Corollary 1.The local delay and the F 1 distribution are related as follows: Proof.See Appendix A.
2) F 2 distribution: Given T ∈ R + , the F 2 distribution is defined as follows: In the context of wireless communications, F 2 (T ) is the probability that the expected delay, conditioned on a network realization and averaged over the fading channels, is greater than the threshold T .Similar to the F 1 distribution, it can be seen as a kind of delay violation probability.
3) F 3 distribution: Given τ ∈ Z + and x ∈ [0, 1], the F 3 distribution is defined as follows: To understand the physical interpretation of F 3 (x, τ ) in wireless communications, consider a reliability metric R(p, τ ) which is a measure of the system capability to transmit a packet within a deadline τ and with success probability greater than p [18].

IV. ANALYTICAL FORMULATION OF THE LOCAL DELAY
In this section, we derive analytical expressions for the local delay for each coverage definition.
Then, we determine the condition under which the local delay is infinite.First, we present the following result on the formulation of the CCP.
Lemma 1.Let G(r 0 ) be a coverage criterion-dependent function defined as: For any coverage criterion with G(r 0 ) defined in (17), the CCP is formulated as follows: Proof.See Appendix B Using Lemma 1, we can formulate the local delay as follows.
Theorem 1.For any coverage criterion, the local delay is given by: Proof.See Appendix C.
The local delay in (19) is given as a simple integral.If we consider the SIR-based coverage, we can obtain a closed-form expression under some conditions.
Proof.According to (17), under the SIR-based criterion, we have G(r 0 ) = 1.Thus, (19) becomes: The integral is finite when Proof.According to (17), under the SINR-based criterion, we have G(r 0 ) = e − γW Kr α 0 P .Thus, the local delay in (19) reduces to: The integral expression is infinite due to the assumption α > 2. When the SIR+ASNR-based coverage is considered, we have G = ½ r 0 ≤ P KW θ 1/α and the local delay in (19) becomes: The last integral on the right hand side of ( 23) is infinite since ½ r 0 ≤ P KW θ 1/α becomes zero and hence D is infinite.This completes the proof.
Remark 1.Let γ * be the SIR critical threshold defined as From Corollary 2, under the SIR-based criterion, the local delay is finite and positive if SIR < γ * , and infinite otherwise.On the other hand, according to Corollary 3, the local delay is always infinite under the SINR-based or the SIR+ASNR-based criterion.

V. ANALYSIS OF THE F 1 DISTRIBUTION
In this section, we derive the general formulation of the F 1 distribution, as well as its analytical expression for each considered coverage criterion.We also provide the asymptotic analysis of the distribution which corresponds to the packet loss probability.

A. Analytical Formulation of the F 1 Distribution
The general formulation of the F 1 distribution is given in the following theorem.
Theorem 2. For any coverage criterion, the F 1 distribution is formulated as: Proof.See Appendix D.
Based on Theorem 2, the analytical expressions of the F 1 distribution for the three considered coverage criteria are given by the following corollary.

B. Asymptotic Analysis of the F 1 Distribution
In this section, we analyze the asymptotic behavior of the F 1 distribution as τ → ∞.In other words, we are interested in computing the following limit: P e = lim τ →∞ F 1 (τ ).Recall that F 1 (τ ) can be regarded as the delay violation probability within the deadline τ .Thus, in wireless communications, P e corresponds to the packet loss probability, i.e., the probability that a user never receives its packet within an infinite waiting time.We first introduce the following results.
Lemma 2. Under the SIR-based coverage criterion, we have P Φ cov > 0.
Proof.See Appendix F.
Proposition 1.For any coverage criterion, the packet loss probability is: Proof.See Appendix G.
Proposition 1 implies that the expected packet loss probability is equivalent to the probability of the CCP being zero, which is equivalent to the probability that G(r 0 ) is zero.
Remark 3. From Corollaries 5 and 6, we can evince the following: • Under the SIR-based or SINR-based criterion, as τ → ∞, all the users are guaranteed to be served eventually, even though the local delay in these cases can be infinite.
• On the other hand, since e −π λ BS( P KW θ ) 2/α > 0 if θ > 0, the packet loss probability under the SIR+ASNR-based criterion is always non-zero.This implies that there is a certain portion of the users who, on average, will not get their packet no matter how long they wait.
The non-zero probability of packet loss under the SIR+ASNR-based coverage criterion can be explained as follows.A necessary condition for a typical user MT 0 to be served by its tagged BS, BS 0 , in a given time slot, is that the ASNR is greater than the detection threshold, i.e., To obtain analytical expressions for the F 2 and F 3 distributions, we first state the following preliminary result.Lemma 3. Define the random variable Z = log(P Φ cov ).For any coverage criterion, the characteristic function of Z, denoted by ϕ Z (t), is given as follows: Proof.By definition, the characteristic function of Z = log(P Φ cov ) is: By substituting the conditional coverage probability P Φ cov in ( 18) into (33), we obtain: where (a) follows from separating the conditioning on Φ into Φ I and r 0 .The proof follows by applying the probability generating functional (PGFL) theorem to the point process Φ I [19] and by computing the expectation over r 0 .
Remark By substituting G(r 0 ) in ( 17) into (32) for each of the coverage criteria, we obtain: Theorem 3. The F 2 and F 3 distributions can be formulated as follows: Proof.Recall that ∆ Φ is a geometrically distributed random variable whose mean is 16), also, we have . The proof follows from the inversion theorem [20] and substituting z = − log T and z = 1 − x 1/τ for the F 2 and F 3 distributions, respectively.

VII. APPROXIMATIONS FOR THE DELAY DISTRIBUTIONS
The analytical formulations for the distributions of the delay derived in the previous section may not be easy to be computed in some cases.In particular, the F 1 distribution in Theorem 2 is difficult to compute for large values of τ due to the binomial coefficients.The difficulty of computing the F 2 and F 3 distributions in Theorem 3 lies in numerically computing, with high precision, the integral expressions in (37) and (38).Motivated by these considerations, in this section, we develop accurate and efficient numerical approximations for computing the F 1 , F 2 , and F 3 distributions.

A. Approximation for the F 2 and F 3 Distributions via the Inverse Laplace Transform
To facilitate the numerical computation of the F 2 and F 3 distributions, we utilize the numerical inversion of the Laplace transform based on the Euler-sum method.This method is based on [21] and it was used in [22] to compute the outage probability of diversity systems over generalized fading channels.More recently, it was also used in [23] and [24].Using this technique, the CDF of a random variable Y can be approximated as where FY (s) is the Laplace transform of F Y (y), E(A, Q, N) is the approximation error that is determined by the three parameters A, N, and Q, and β n is defined as β n = 2 if n = 0 and β n = 1 if n > 0. In this paper, we use the same parameters as in [22], namely A = 10 log 10 ≈ 23.03, which guarantees a discretization error of the order of 10 −10 , as well as N = 21, and Q = 15 to ensure that the resulting truncation error is less than 10 −10 .
Proposition 2. The F 2 and F 3 distributions can be approximated as follows: Proof.Define Y = − log(P Φ cov ).From ( 15), the F 2 distribution can be written as: Since F Y (0) = P[− log(P Φ cov ) < 0] = P P Φ cov > 1 = 0, the Laplace transform of F Y (y), i.e., FY (y), can be formulated as FY (s) = ϕ Z (−is) s .The final approximation follows by substituting y = log(T ) and F (s) = ϕ Z (−is)/s into (39) and plugging the resulting expression of F Y (y) into (42).The proof for the F 3 distribution is obtained similarly by setting y = − log 1 − x 1/τ .Remark 5.The approximations in Proposition 2 have a singularity at T = 1 when computing the F 2 distribution and at x ∈ {0, 1} when computing the F 3 distribution.This, however, poses no problems, since F 2 (T = 1) = 1, F 3 (x = 0, τ ) = 1, and Therefore, there is no need to apply Proposition 2 for these specific values.

B. Approximation of the F 2 and F 3 Distributions by Using the Beta Distribution
The approximations of the F 2 and F 3 distributions provided in Proposition 2 require the computation of the characteristic function ϕ Z (•).If the SINR-based coverage criterion is considered, ϕ Z (•) is available only in integral form, as shown in ( 35).An alternative approach to approximate the F 2 and F 3 distributions is to utilize the beta distribution, which is shown in [11] to offer accurate estimates for distributions whose support lies in the entire range [0, 1].Proposition 3. Let µ and ν be be the mean and variance of P Φ cov , respectively, given as follows: As for the SINR-based coverage, the F 2 and F 3 distributions can be approximated as: where Proof.From ( 15) and ( 16), the F 2 and F 3 distributions can be written as F 2 (T ) = P P Φ cov ≤ 1/T and F 3 (x, τ ) = P P Φ cov ≤ 1 − x 1/τ , respectively.Under the SINR-based coverage criterion, P Φ cov lies in the range [0, 1].Thus, F 2 (T ) can be approximated with a Beta distribution whose shape parameters a and b are obtained by matching the mean µ and variance ν of Remark 6.As remarked in [25], the Beta approximation can be applied only if the original distribution falls within the class of Beta distributions, e.g., the domain of the random variable is the entire interval [0, 1].As shown in [24], under the assumption of SIR+ASNR-based coverage, on the other hand, the conditional coverage probability does not necessarily fall within this class.Therefore, the Beta approximation may not be applicable.
Remark 7. The Beta approximation in Proposition 3 relies on matching the shape of the Beta distribution to that of the original F 2 and F 3 distributions.This is different from the Euler-sum inversion method in Proposition 2, which is an efficient numerical computation of the original distribution with a known truncation error.Therefore, when both approaches are applicable, one is encouraged to use the Euler-sum inversion method for a better accuracy.

C. Approximation of the F 1 Distribution via the F 3 Distribution
The F 1 and F 3 distributions are related to each other through the following integral relation: where (a) follows from the definition of the and writing the outer expectation operator as the integral of the CDF.Let p = {p i } be an ordered sequences of real numbers between 0 and 1, i.e., 0 = p 0 < p 1 < . . .< p n−1 < p n = 1.Also, let ǫ = {ǫ i } be a sequence that satisfies p i−1 < ǫ i ≤ p i , ∀i ∈ {1, 2, . . ., n}.The F 1 distribution can be expressed as a limit of a sum of the F 3 distribution, i.e., Remark 8.By using (47) with a sufficiently large value of n, we can compute the F 1 distribution efficiently.By using this technique, we avoid working with large values of the binomial coefficients that appear in Theorem 2, which usually leads to numerical issues in evaluating the distribution.

VIII. NUMERICAL RESULTS
In this section, we show numerical results to validate the proposed analytical frameworks for computing the local delay and the delay distributions, as well as to substantiate the obtained findings on the packet loss probability.Unless otherwise stated, the simulation setup is summarized in Table II.The simulation results presented in this section are obtained through Monte  Carlo simulations, by generating 5, 000 realizations of the BS and MT, as well as simulating 5, 000 packet transmissions for each point process (spatial) realization.λ BS .The theoretical results are computed by using (21).From the figure, we can see that the local delay increases with the BS density and the decoding threshold.This reveals that, given a spatial realization of BSs and MTs, a typical MT is more likely to successfully delivers its data packets when the number of BSs increases due to shorter distances, on average, despite the MT is likely to experience more interference due to higher number of interfering BSs.We observe, in addition, that the local delay tends to infinity near the critical threshold given in (24).the figure, we can see that the analytical formulations agree with the simulation results.We also observe that the delay performance is better for smaller values of the threshold and higher values of the BS density.Fig. 4 shows the F 1 distribution as τ tends to infinity, i.e., it illustrates the packet loss probability.As evident from the figure, the probability goes to zero as τ goes to infinity when the SIR and SINR-based coverage criteria are considered, while it goes to a non-zero value when the SIR+ASNR-based coverage criterion is utilized.We also observe that the delay performance improves when the transmit power P increases.From the figures, we can observe that the delay performance is better for lower values of the threshold and for higher values of the BS density.We also observe that, in a fully-loaded scenario (i.e., when λ BS ≪ λ MT ), the impact of the noise (for the SINR-based and SIR+ASNR-based coverage criterion) is not apparent.This is because, in this case, the system is nearly interference-limited, i.e., the interference is a more dominant factor as compared to the noise.On the other hand, the impact of the noise is clearly seen in lightly-loaded scenario (i.e., when λ BS ≫ λ MT ), since several BSs may be turned off.

D. F 3 Distribution
Fig. 6 validates the computational framework for the F 3 distribution as a function of the number of time slots τ ∈ {5, 10} and the decoding threshold γ D when the BS density is λ BS = 0.1 λ MT .
The approximated results are computed by using the inverse Laplace transform method (for the SIR-based and SIR+ASNR-based coverage criteria) and the Beta approximation (for the SINRbased coverage criterion).From the figures, we can see that the analytical formulations agree     method, the Riemann sum, and the Beta distribution have been employed to circumvent some inherent numerical difficulties in computing the obtained analytical formulations for some case studies, as well as to render the computation of the proposed analytical methods more efficient.
Also, we have analyzed the asymptotic behavior for one of the proposed delay distributions, which provides information on the network packet loss probability.Finally, through extensive simulations, we have shown that the obtained analytical formulations and approximations are accurate compared to Monte Carlo simulation results, and have studied the impact on the delay performance of several design parameters, such as the decoding threshold, the transmit power, and the network deployment density.The proposed approaches can facilitate the analysis and optimization of cellular networks subject to reliability constraints on the network packet delay that are not restricted to the local (average) delay, e.g., for delay sensitive applications.

APPENDIX C PROOF OF THEOREM 1
By substituting ( 18) into (12), and separating the conditioning on Φ into r 0 and Φ I , we have: where (a) follows from the fact that the term (G (r 0 )) −1 is independent of Φ I .According to the probability generating functional (PGFL) theorem of a PPP [19], for a general function f (x) and a PPP Ψ over R 2 , we have: where λ(x) is the density function of the PPP.By expressing the whole domain R 2 in polar coordinate, we can write dx = rdrdφ where φ is the azimuth angle of a particular point with respect to the origin.Using this notation, we obtain: where (b) is obtained by substituting λ(r i , φ) with ( 5) .The proof follows by substituting (54) into (52) and by then calculating the last expectation over r 0 with the aid of the PDF in (3).

APPENDIX D PROOF OF THEOREM 2
Since ∆ Φ is a geometrically distributed random variable, its conditional CCDF is: From (13) and by separating the conditioning on Φ into r 0 and Φ I , we obtain: We start the proof by first proving the following intermediate result: if {a i } is an infinite sequence such that 0 < a 1 < a 2 < . . .< 1, then the following holds: First, we note that, due to the equivalence of implication, (58) is equivalent to ∞ i=1 a i ≤ 0 ⇐⇒ ∞ i=0 (1 − a i ) ≥ ∞.Also, since a i > 0 for all i, the latter implication is equivalent to: To prove (58), it is sufficient to prove (59).We consider two cases.
1) Case 1: {a i } converges to 0 < M < 1: Assume that {a i } converges to M where 0 < M < 1, that is, lim i→∞ a i = M < 1.Since a i < M for all i, we have ∞ i=1 a i < ∞ i=1 M (a) = 0 where (a) follows from the fact that 0 < M < 1.Since lim i→∞ a i = M < 1, we have lim i→∞ (1 − a i ) > 0.
2) Case 2: {a i } converges to M = 1: Assume that {a i } converges to M = 1, that is, where (b) follows from taking the logarithm of both sides of the equation and (c) follows from the fact that the logarithm of a product is equal to the sum of the logarithm as well as log(0) = −∞.
Let us rearrange the indices of the interferers i such that they are ordered from the closest to the furthest from the typical user.Then, under the SIR-based coverage criterion, the CCP in (8) can be written as follows: (60)

Corollary 3 .
is infinite since the exponential term of the integrand function grows without bound.Computing the integral completes the proof.Under the SINR-based or the SIR+ASNR-based criterion, the local delay is infinite.

Remark 2 .
From Corollary 4, we can see that, under the SIR-based or SIR+ASNR-based criterion, F 1 (τ ) is given in a closed-form expression.Under the SINR-based criterion, on the other hand, F 1 (τ ) is not formulated, in general, in a closed-form expression.Since F (k, α, γ) ≥ 1 and (26), (27), and (28) are expressed as a finite sum, the F 1 distribution is always finite for all coverage criteria, even though the corresponding local delay may be infinite.

Fig. 1 :
Fig.1: Illustration on the non-zero packet loss probability under the SIR+ASNR-based coverage criterion.Each dot and square represents a BS and an MT, respectively.Both the red and green MTs are associated to the black BS.The dashed circle represents the points with distance r * = P KW θ 1/α from the green BS.The red MT will never receive its packet due to failure in fulfilling the coverage condition.

Fig. 2 :
Fig. 2: Local delay as a function of the decoding threshold (SIR-based coverage criterion).

Fig. 2
Fig. 2 validates the analytical framework for the local delay that corresponds to the SIRbased coverage criterion, as a function of the threshold γ D for different values of the BS density

Fig. 4 :C. F 2 Fig. 5
Fig. 4: F 1 distribution for different coverage criteria as a function of number of time slots, decoding threshold, and BS density.

TABLE I :
Main symbols and functions used in the paper

TABLE II :
Simulation setup.