{'fraction_non_alphanumeric': 0.04367241700114298, 'fraction_numerical': 0.012944897796443476, 'mean_word_length': 4.20359612724758, 'pattern_counts': {'":': 0, '<': 22, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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.

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 Hypergeometric function, Gamma function, indicator function ϕZ (·), MY (·) Characteristic function of Z, moment-generating function of Y ℑ{·}, ℜ{·} Operators returning the imaginary part, real part of the argument 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 homoge- 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 P Φ [·] 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: SINR-based coverage, 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.
Under the SIR-based coverage criterion, the local delay has the following closed-form expression: Proof. According to (17), under the SIR-based criterion, we have G(r 0 ) = 1. Thus, (19) becomes: The integral is finite when Otherwise, it is infinite since the exponential term of the integrand function grows without bound. Computing the integral completes the proof.

Corollary 3. Under the SINR-based or the SIR+ASNR-based criterion, the local delay is infinite.
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 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.
Proof. F sir 1 (τ ) and F sinr 1 (τ ) follow by inserting G(r 0 ) in (17)  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 θ )

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., is fixed. Therefore, given a pair BS 0 and MT 0 , the MT 0 will never receive its packet (i.e., it experiences a packet loss), if r 0 > r * where r * = P KW θ 1/α (see Fig. 1 for illustration). In contrast, the condition for which the user is served according to the SINR-based coverage criterion is P |h 0 |/Kr α 0 i∈Φ I P |h i |/Kr α i +W > γ. Due to the random fading coefficients |h 0 | and |h i |, there is a non-zero probability for this condition to be fulfilled, which results in a zero probability of packet loss for a long enough waiting time.

A. Analytical Formulation of the F 2 and F 3 Distributions
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 . 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 P Φ cov −1 . Let (16), also, we have F 3 (x, τ ) = P Φ 1 − P Φ cov τ ≥ x = P Z ≤ log 1 − x 1/τ = F Z 1 − x 1/τ . 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 whereF Y (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: (41) Proof. Define Y = − log(P Φ cov ). From (15), the F 2 distribution can be written as: 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: 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 P Φ cov , i.e., a = µb/(1 − µ) and b = µ(1 − µ) 2 /ν − (1 − µ). [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 6. As remarked in
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 F 1 distribution, i.e., 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. 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 λ 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).     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. 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,

A. Local Delay
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.     By virtue of the linearity of the expectation operator, we can rewrite the local delay in (12) as follows:

D. F 3 Distribution
The proof follows from the definition of probability mass function (PMF) of the discrete random When the SINR-based coverage criterion is considered, the CCP is defined as follows: where (a) comes from the definition of SINR and P Φ cov from the assumption that the random variable h 0 is exponentially distributed with unit mean, (b) is obtained by substituting the total interference I and by taking into account that the term exp(−γW ℓ(r 0 )/P ) is independent of I, (c) follows from the fact that the fading coefficients h i are mutually independent, and (d) follows because h i is an exponentially distributed random variable with unit mean.
The CCP for the SIR-based coverage criterion can be derived by setting W = 0. When the SIR+ASNR-based coverage criterion is considered, the CCP is: where (e) holds because W is independent of h 0 and I. The proof follows using similar steps as for the SINR-based coverage and by noting that P [W ≤ P/(Kr 0 α θ)|Φ] = ½ r 0 ≤ P 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: where (a) is obtained by substituting the CCP given in Lemma 1 and (b) follows by applying the PGFL theorem and using the notable integral The proof follows from deconditioning the expectation on r 0 .
APPENDIX E

PROOF OF COROLLARY 4
As for F sir+asnr 1 (τ ), we first substitute G(r 0 ) = ½ r 0 ≤ P KW θ 1/α into (25) and obtain: where (a) follows from two steps: (i) splitting the sum for k = 0 and for k ≥ 1 as well as using the fact that ½ r 0 ≤ P KW θ 1/α 0 = 1, and (ii) by taking into account that, for k ≥ 1, 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.
2) Case 2: {a i } converges to M = 1: Assume that {a i } converges to M = 1, that is, log(a i ) = −∞, 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) = −∞.
Without loss of generality, let us assume that a 1 > K for some 0 < K < 1. Due to the fact that {a i } is an ordered sequence, we have a i > K for all i and, hence 1/a i < 1/a 1 = 1/K.
Multiplying both sides with a i − 1 < 0, we have a i −1 a i > a i −1 K . According to [26,Eq. (1)], for any a i > 0, we obtain a i −1 a i ≤ log (a i ) ≤ a i − 1. Then, we obtain a i −1 K ≤ log (a i ) ≤ a i − 1. To prove (59), we need to prove two implications: (1) log(a i ) = −∞. First, let us assume that ∞ i=1 log(a i ) = −∞ holds. From the inequality a i −1 log(a i ) = −∞. Based on these results, we conclude that (59) (and hence (58)) holds for when {a i } converges to M = 1. This completes the proof of (58).
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: We assume that no two interferers have exactly the same distance from the typical user, i.e., r 1 < r 2 < . . . < r N . Denote a n = 1 + γ r 0 rn α −1 . Thus, {a n } is an ordered sequence where 0 < a 1 < a 2 < . . . < 1. Clearly, P Φ cov > 0 if N < ∞. Assume N = ∞. According to (58), for P Φ cov to be non-zero, it must hold that ∞ n=0 (1−a n ) < ∞. Denote b n = γ r 0 rn α / 1 + γ r 0 rn α .