SIR Coverage Analysis in Multi-Cell Downlink Systems With Spatially Correlated Queues

In this article, we analyze the signal-to-interference ratio (SIR) coverage probability of a multi-cell downlink system with random traffic arrivals. Based on a theoretical model, in which the base stations (BSs) and devices are deployed as independent Poisson point processes (PPPs), we first present the sufficient and necessary conditions of the $\varepsilon $ -stable region for the queueing system. Then, by taking into account the impact of spatially queueing interactions among BSs, we focus on the steady status of the queues and introduce a two-queue-length approximated model for the system. Specifically, in the studied new model, the BSs are separated into two sets in the steady state: the long-queue BS set and the short-queue BS set. The locations of the former are modeled by a Neyman-Scott process and those of the latter are modeled as a residual hole process. By applying the first-order statistic approximation, we further approximate the hole process by a homogeneous PPP with the same density. To model the deployment of the BSs precisely, the related parameters in the approximated model are fitted. Using tools from stochastic geometry and queueing theory, we derive the SIR coverage probability in the steady state. An iterative algorithm is proposed to calculate the active probability of the BSs. To reduce the computational complexity, Beta distribution is applied to approximate the probability density function of the service rate in each iteration of the algorithm. Finally, the effect of the fitting parameters and the accuracy of our analysis are presented via Monte Carlo simulations.


I. INTRODUCTION
The massive connectivity of mobile devices, the demand for high data rates, along with the scarce spectrum have imposed complicated interference situations to modern wireless systems. To stipulate the quality-of-service for the devices in the sophisticated wireless networks, coverage probability is no doubt an important metric. To provide a unified mathematical paradigm that characterizes the coverage probability, stochastic geometry has been widely used in various wireless systems, e.g., the device-to-device networks [1]- [4], the broadcast networks [5]- [8], and the heterogeneous networks [9]- [12], etc.. The associate editor coordinating the review of this manuscript and approving it for publication was Chien-Fu Cheng .
Nowadays, with the development of the multi-media and the evolution of mobile applications, a large amount of data from different sources with different traffic types are transmitted simultaneously in the same communication systems. In these systems, the impact of the data traffic on the network services is no longer negligible. On the one hand, the service requirements of different types of applications may vary greatly. For instance, for latency-sensitive applications, e.g., immersive virtual reality and automotive driving, the delay is required to be less than 2 ms [13], whereas the ones with a lower degree of sensitivity upon latency, e.g., video conference only requires the delay to be less than 100 ms [14]. On the other hand, the extremely increased network throughput requires us to explore the potential of the performance gains through the cross-layer designs. The conventional approaches of network performance analysis [5]- [7] heavily rely on the full buffer assumption, i.e., each base station (BS) is backlogged and keeps transmitting all the time. In comparison, interacting queues are considered in [15] with a simplified collision model, which allows one to quantify the network performance with a combination of the queueing theory and stochastic geometry. In the wireless networks, the spatial distance and channel gain related interference affects the transmission success probability of the packets (or the states of the interacting queues). Conversely, the buffer states in the previous time slot also affect the activation of the interfering BSs and the interference of the devices. In view of this, the interactions among the queues are dependent on both the spatial and temporal factors. Considering the random arrivals of the packets, a recent line of research has introduced a spatiotemporal model to analyze the coverage probability from a joint spatial-temporal aspect [16]- [18]. To simplify the analysis, all the previous works assume that the queueing evolution is independent and identically distributed (i.i.d.) among all the BSs, which is commonly known as mean field approximation.
However, in practical systems, multiple BSs sharing the same spectrum inevitably interfere with each other due to the broadcasting nature of wireless channels. As such, the queueing status of the BSs are correlated in both the time domain and space domain among all the BSs. In fact, even the traffic arrival pattern at each BS is the same, the spatial interaction of queues can lead to a very different queue length after a sufficiently long time. Specifically, for the BSs located in a crowded area, the severe mutual interference during the transmissions slows down the rate of the services and results in longer queue lengths. The long queue lengths increase the active probability of the BSs and further strengthen the mutual interference. On the contrary, for those BSs located far away from their neighbors, the queue lengths are always short, which further lighten the interference and improve the service rate [18]. An example of such a phenomenon is illustrated in Figure 1. As observed, when the queues are stable, the BSs with long queue lengths are gathered in three clusters, whereas the BSs with short queue lengths are scattered in other uncrowded space. Note that the effect of spatial interaction is stronger in the networks with denser spatial deployment.
This paper aims to study the SIR coverage probability for the multi-cell downlink system. Specifically, the spatiotemporal correlations amongst the queues are considered. The temporal correlations of the queues are handled by analyzing the dynamic interacting queues as in [18] whereas the spatially queueing interactions are revealed by introducing a new model of the BSs. Inspired by the observations of Figure 1, we adopt a new framework that separates the BSs into the longqueue BSs and short-queue BSs and considering the two sets of the BSs separately.

A. RELATED WORKS
The coverage probability is widely analyzed in the downlink multi-cell cellular systems [5]- [7]. In [5], the independent PPP placed BSs were considered and a tractable expression for the coverage probability was derived under a new framework for the multi-cell signal-to-interference-plus-noise ratio (SINR). The work in [6] proposed a generic uniform density plane-entity method to evaluate the inter-cell interference. Based on the results, the lower bounds of the coverage probability and the achievable throughput were respectively derived under three different BS distribution models. Reference [7] considered the Poisson-Poisson cluster processes model for the BS locations and presented a numerically computable coverage probability of the system. All the previous works assumed that the BSs are always in the full-buffer state.
To analyze the effect of the traffic, a variety of works consider the wireless systems with a spatiotemporal model. The work in [16] evaluated the delay performance of three different scheduling methods in the heterogeneous cellular networks. By leveraging properties from Geo/PH/1 queues, [17] analyzed the signal-to-interference ratio (SIR) performance under three different transmission strategies for uplink Internet of things (IoT) embedded cellular networks. Based on the SINR meta distribution, [18] derived the coverage probability for users to achieve different SINR levels in the downlink cellular networks. With the spatiotemporal model, [19] characterized the link quality of the nodes and further designed the uncoordinated multiple access strategies for wireless networks with massive transmitters.

B. CONTRIBUTION
In this paper, we study the SIR coverage probability for a multi-cell downlink system, in which the BSs and the devices are scattered according to independent Poisson point processes (PPPs). Similar to [18], we set a random arrival of the traffic and a dynamic interacting queue at each BS. The packet arrivals are assumed to follow independent Bernoulli processes. However, different from [18] and [19] that apply the mean field approximation, we develop a new model to characterize the spatial correlation of the queueing states. Specifically, we divide the BSs into two sets based on their VOLUME 8, 2020 steady queue states, i.e., the long-queue BS set and shortqueue BS set. As will be shown later, the locations of the long-queue BSs are approximated by a Neyman-Scott process (NSP), whereas that of the short-queue BSs form a residual hole process.
The technical contributions are two folds. First, we approximate the residual hole process with an independent PPP based on the first-order statistical approximation [20] and analyze the SIR coverage probability when the queues are in the steady state. Using tools from stochastic geometry and queueing theory, we derive a closed-form expression for the conditioned Laplace transform of the long-queue BSs interference. Regarding the conditioned Laplace transform of the short-queue BSs interference, an algorithm is proposed to update the active probability of the BSs iteratively. Second, to obtain a tractable expression for the coverage probability, Beta distribution is adopted to approximate the probability density function (PDF) of the service rate in each iteration of the algorithm. Simulation results are presented to validate the accuracy of the approximated model, discuss the effect of the fitting parameters and show the convergence of the algorithm. Besides, the simulations also verify that the analytical result by applying the proposed approach matches well with the real network deployment.
The remainder of the paper is organized as follows. Section II describes the system model of the multi-cell downlink system. Section III analyzes the stable region for the queues and presents the two-queue-length model. In Section IV, we derive the coverage probability on the condition of stable queues. Numerical results are shown in Section V and the conclusion is presented in Section VI. Figure 2, we consider the transmissions in a multi-cell downlink system which consists of multiple BSs and devices all equipped with a single antenna. The spatial locations of the BSs follow a homogeneous PPP b of density λ b . The devices are served by the nearest BS with a constant transmission power P. In this network, the spectrum is divided into orthogonal subchannels and reused in each cell for improving the spectrum efficiency. Similar to the setting in [21], we focus on a typical subchannel and assume that the density of the devices is large enough such that each BS has one associated device on the subchannel. 1 Suppose the typical device and its serving BS locate at the origin and position z ∈ b , respectively. At time slot t, the SIR of the typical device is expressed as

As illustrated in
where h z ∼ exp(1) is the small scale fading from the BS located at z to the origin that follows the exponential distribution with density one, · is the Euclidean norm, α is the pass loss exponent, and I is the normalized interference power 1 Following the idea in [18], our analysis can be easily extended to the case that multiple devices share one subchannel for each BS. received at the typical device, which will be analyzed later in Section IV.
In the media access control (MAC) layer, we use a discretetime queueing system to model the traffic dynamics. Specifically, at any time slot t, each BS has a potential of the packet arriving with probability (packet arrival rate) ξ ∈ [0, 1]. The incoming packet is stored in an infinite-sized buffer. Then, the packets in the buffer are scheduled by the first-in-first-out (FIFO) policy, i.e., if the buffer is non-empty, the head-of-line packet of the buffer is transmitted from the BS to the intended device. The packet is successfully transmitted when the SIR received at the intended device exceeds a decoding threshold θ, i.e., γ z,t ≥ θ. After a successful transmission, the target packet is removed from the buffer and the transmission for a new packet begins. Otherwise, the failed transmission will resume in the next time slot until successful. We define the conditional successful transmission probability of the serving BS at time slot t as Similar to [16], we consider the BSs with empty buffers entering into the inactive state to save energy and reduce the interference. Here, we define a binary indicator ζ z,t to describe the transmission state of a generic BS located at z ∈ b in time slot t as the following ζ z,t = 1, the BS located at z in time t is active, 0, otherwise.
In the next section, we will first analyze the conditions under which the queues across the network can be stable.
Then, under such conditions, we will consider the spatial correlation among the queues and present an approximated two-queue-length BS model.

A. STABLE CONDITION
Since different devices suffer different interference levels, the number of the packets in the queue of each BS varies over different time slots. Based on Loyen's Theorem [22], the queue of an isolated system is stable if and only if the averaged services rate is larger than the average packet arrival rate over all the system operating time. Unfortunately, the strict stability that all the queues are finite-length is not achievable in the large scale Poisson network (except for the trivial case ξ = 0). For instance, when some randomly deployed devices locate at a poor SIR covered area, those devices will have an infinite queue length in the long run. Thus, given packet arrival rate ξ , successful transmission probability µ b z,t and system operating time T , we follow [23] and introduce the concept of ε-stability as follows Definition 1: For any 0 ≤ ε ≤ 1, the ε-stable region S ε and the critical arrival rate ξ c,ε of the static Poisson network are respectively defined as We remark that the system is ε-stable if and only if ξ ≤ ξ c,ε . In fact, it is non-trivial to obtain the exact expression of ξ c,ε . Instead, the sufficient and the necessary condition for ε-stability are given by the following proposition Proposition 1 [18]: The sufficient and necessary condition for the network stability is respectively given by ξ ≤ sup ξ ∈ R + : ξ ≤ sup ξ ∈ R + : where Im{·} is the imaginary part of the complex number, is the hypergeometry function [5].
The sufficient or necessary condition holds by letting all the BSs keep transmitting signals or regarding all the packets are successfully transmitted within the operating time. Comparing with the actual scenario, the interference suffered by the devices is more severe in the former case whereas the interference is more slight in the latter case. More details can be found in [18].  Table 1 presents the percentage of the ε-stable queues in the case ε = 0.01 by the numerical simulations. Table 2 shows more details of the simulation parameters. As seen, the queue becomes more unstable as the value of packet arrival rate ξ or the decoding threshold θ goes up. In this paper, we always focus on the cases that the queues are ε-stable.

B. TWO-QUEUE-LENGTH BS APPROXIMATION
Unlike the full-buffer case, where the BSs always keep transmitting packets, the active states of the BSs in our work are changed with the dynamic queues. Intuitively, in a static Poisson network, a long-queue BS tends to have a long period of time operating in the active state, which incurs a high interference to the nearby BSs and consequently decrease their successful transmission probability. Accordingly, the queue lengths of those nearby BSs are to be prolonged. In view of this, instead of assuming the queues evolve independently from each other in space domain [18], we consider that the queues of the BSs are spatially correlated. As illustrated in Figure 1, the spatial correlation in the buffer states result in a cluster pattern of the long-queue BS locations.
Spatial factors aside, the traffic pattern also has a noticeable effect on the correlation of queueing states. The correlation of the queueing states can be easily verified by the following simulation instance. Let a t and b t respectively be the queue lengths of two BSs at time slot t, a = (a 1 , . . . , We define the correlation coefficient between a and b as follows: Suppose the BSs and the devices are deployed over a square area of 1000 m 2 with BSs density λ b = 0.1 BS/m 2 , decoding threshold θ = 0 dB and system operating time T = 1000. Figure 3 depicts the correlation coefficient of the queue lengths between the serving BS and its' five nearest interfering BSs. As shown, the correlation of all the interfering BSs increases dramatically as the traffic load of the system VOLUME 8, 2020 goes up. Specifically, in the randomly placed BSs scenario, the correlation coefficient increases from 0.2 to 0.7 as the packet arrival rate ξ changing from 0.5 packets/slot to 0.9 packets/slot. Besides, the correlation of the queues is more evident in the randomly placed BSs scenario than the hexagonally placed BSs scenario. Nearly 75% increment of the correlation between the serving BS and the interfering BS occurs when the BS locations changing from hexagon into random.
The above instance implies that the spatial correlation of the queues between the serving BS and the interfering BSs in its' neighbor cannot be ignored. Note that directly deriving the Laplace transforms in Section IV with spatial correlated queues is not accessible. Instead, we approximate the system with a two-queue-length model motivated by the observations above. Specifically, the original PPP placed BSs are separated into two sets based on their stable states of the queue lengths: the long-queue BSs set and the short-queue BSs set. Then two different point processes are respectively used to model the locations of the two BS sets.
Denote l and s as the locations of the long-queue and short-queue BSs, respectively. To be concrete, we have b = l ∪ s . We approximate the locations of the longqueue BSs by an NSP: a clustering process whose parent process is a PPP l c with intensity λ p , whereas the offspring point process l x for each parent x ∈ l c is the conditionally independent bivariate Gaussian process with zero mean and variance σ 2 I ∈ R 2 . The number of the offsprings is a Poisson random variable (r.v.) with meanm. On the other hand, the locations of the short-queue BSs satisfy a residual hole process with intensity λ s b = λ b −mλ p . Mathematically, the exact characterization of the hole process and its' dependency with NSP is not available. In order to obtain a tractable expression of the coverage probability, we assume that the NSP and the residual hole process are independent and apply the firstorder method to approximate the hole process [20]. That is, the locations of the short-queue BSs are approximated by a PPP with the same intensity λ s b . The accuracy of these approximations will be validated in Figure 4.

IV. COVERAGE PROBABILITY
In this section, we present the main technical part of this work, i.e., deriving the coverage probability of the employed network in the steady state. 2 To begin with, let us denote y x as the index of the BS located in cluster l x with x and y being the cluster center of l x and the relative position from cluster center x to BS y x , respectively. Suppose the typical device is served by a long-queue BS located at y 0x 0 , the normalized interference of the typical device at time slot t contains the following three parts: • Intra-cluster interference: the normalized interference from other long-queue BSs in the same cluster • Inter-cluster interference: the normalized interference from the long-queue BSs in other clusters • short-queue BS interference: the normalized interference from the short-queue BSs By plugging I = I l-intra + I l-inter + I l-short , (8), (9) and (10) into (1), the coverage probability attained at the typical device in the steady state is given by where γ y 0x 0 is the SIR of the device linked to BS y 0x 0 ∈ b in the steady state, r = x 0 + y 0 is the serving distance from the serving BS to the origin, v 0 = x 0 is the distance from cluster center x 0 to the origin, s = θr α , L z (s | r) is the conditioned Laplace transform of r.v. z at point s in the steady state, f R (r | v 0 ) is the conditioned serving distance distribution, and . On the other hand, if the serving BS is a short-queue BS y 0 , the normalized interference of the typical device contains two parts: • Interference from long-queue BSs: the normalized interference from the long-queue BSs • Interference from short-queue BSs: the normalized interference from other short-queue BSs By plugging the equation I = I s-long + I s-short , (12) and (13) into (1), when the typical device is served by short-queue BS y 0 , the coverage probability in the steady state is therefore given by where f R (r) = exp(−λ s b πr 2 )2πλ s b r [5]. In the next two subsections, we would like to derive the mathematical expressions of the conditioned serving distance distribution and the conditioned Laplace transform in (11) and (14) so as to obtain the coverage probability of the serving BS in the steady state.

A. DISTANCE DISTRIBUTION
Suppose W out and S x respectively be the set of the distance from the BSs that lies outside the serving distance r in cluster l x 0 to the origin and the distance from the long-queue BSs in cluster l x (x = x 0 ) to the origin. Let v = x , then we have the following Lemma: Lemma 1: w out ∈ W out , s ∈ S x are conditionally i.i.d. and Proof: The proof of Lemma 1 can be found in [24]. We remark that Lemma 1 gives the distance distribution for the intra-cluster interfering BSs and inter-cluster interfering BSs, which will be used to calculate the conditional Laplace transform L I l-intra (s | v 0 , r) and L I l-inter (s | r) in the next subsection. Now we focus on calculating the conditioned distribution of the serving distance f R (r | v 0 ). We first consider a circle area B(0, R) centered at the serving BS y 0x 0 with radius R. Let˜ l x 0 = l x 0 ∩ B(0, R) be the set of the long-queue BS in both cluster l x 0 and area B(0, R);S x 0 be the set of the distance from the long-queue BSs in˜ l x 0 to the origin. The number of BSs in˜ l x 0 equals to k = ˜ l x 0 . Define a set of the distance {s(i)} k i=1 such that the values of the k elements in setS x 0 are sorted in ascending order, i.e., s(1) ≤ . . . ≤ s(k). Since the typical device links with the nearest BS, we have the serving distance r = s (1). The conditional serving distance distribution can be derived by using the properties of order statistics and enlarging the radius of B(0, R) to infinity. The main results are presented in the following Lemma Lemma 2: The distribution of the serving distance conditioned on v 0 = x 0 is given by

B. LAPLACE TRANSFORM
Using the distance distribution derived in the previous subsection, we now calculate the conditional Laplace transform of the interference. Two cases that whether the serving BS is a long queue length BS or a short queue length BS are separately considered.
We first analyze the case that the serving BS is a long queue length BS. For the long-queue BSs, we have ζ y x 0 ,t = 1 in the large time horizon, i.e., t 1. Hence, lim t→∞ P(ζ y x 0 ,t = 1) is a constant that equals to one. In the sequel, we focus on the steady state, i.e., t → ∞, and drop the subscript t. Given that the serving BS belongs to the cluster center which is the closest to the origin among all other cluster centers, the conditional Laplace transform of intra-cluster interference in the steady state is given by the following Lemma Lemma 3: Suppose the serving BS is a long-queue BS. The conditional Laplace transform of the intra-cluster interference and the inter-cluster interference in the steady state are respectively given by where the expression of f W out (w out | v 0 , r) and f U (u | v) is obtained by (15) and (16), respectively. Proof: See Appendix B. Comparing with (18), the Laplace transform of the intercell interference dose not related to the order of the distance. Besides, it gets one more integration that relates to the distance of the cluster center. VOLUME 8, 2020 Unlike the long-queue BSs, the active state of a shortqueue BS varies across different time slots. We define q y = lim t→∞ P(ζ y,t = 1) as the active probability for BSs y ∈ s in the steady state and assume that ζ y for all BSs y ∈ s are independent. Suppose the serving BS is a longqueue BS, we now present the conditional Laplace transform of the short-queue BSs interference in the steady state.
Lemma 4: When the serving BS is a long-queue BS, the conditional Laplace transform of the short-queue BSs interference in the steady state is Proof: See Appendix C. A complete expression of (20) requires us to calculate the term q y . To do this, let us first denote µ b y to be the conditional SIR coverage probability (service rate) of the devices served by BS y ∈ s in the steady state, i.e., µ b y = P γ y ≥ θ | b . Furthermore, we denote f b y (u) as the PDF of µ b y . The expression of q y can be derived as follows.
Lemma 5: The active probability for BS y ∈ s in the steady state is given by Proof: Equation (21) is a standard result of Geo/G/1 queue [25].
Lemma 5 shows that the service rate of the short-queue BSs in the steady state is always larger than the packet arrival rate, i.e., µ b y > ξ . It is non-trivial to obtain the closed-form expression of f b y (u) (or q y ). Instead, we obtain the value of f b y (u) (or q y ) iteratively. Let us denote f b y,n (u) and µ b y,n as the expression of f b y (u) and the value of µ b y in the n-th iteration, respectively. Now we focus on calculating f b y,n . Define Y b y,n = ln(µ b y,n ) = ln(P(γ y,n ≥ θ) | b ), then we have the following Lemma Lemma 6: The l-th moment generation of Y b y,n is given by where q y,n is obtained by (21) with changing f b y (u) into f b y,n−1 (u).
Proof: See Appendix D. We remark that according to the relationship Y b y,n = ln(µ b y,n ), we have M Y b y,n (l) = E(e lY b y,n ) = E({µ b y,n } l ). Therefore, to obtain the mean and the variance of r.v. µ b y,n , it is necessary to calculate the moment generation function of Y b y,n , i.e., Suppose F b y,n (u) to be the cumulative distribution function (CDF) of µ b y,n . As seen shortly, the result in Lemma 6 will be used to derive F b y,n (u). A closed form expression of F b y,n (u) can be obtained by Euler's formula and Gil-Pelaez theorem [26] F b y,n (u) = P(Y b y,n < ln u) The value of f b y,n (u) is therefore obtained by Hence, given the active probability of the BSs q y,0 = ξ , the value of q y is obtained by calculating (22), (23) and (24), iteratively. More details can be found in Algorithm 1. n ← n + 1 7: until q y,n − q y,n−1 ≤ ; 8: Output q y = q y,n .
Due to the high computing complexity of (23), we approximate f b y,n (u) by the PDF of a Beta distribution with the same mean and variance of µ b y,n [18]. That is, in the n-th iteration, we approximate µ b y,n byμ b y,n ∼ Beta(a y,n , b y,n ), with parameters a y,n and b y,n satisfies (1)} 2 = a y,n b y,n (a y,n + b y,n ) 2 (a y,n + b y,n + 1) .
By solving (26) and (27), we obtain In summary, by applying Iterative Algorithm and plugging (15), (16) into (18), (19) and (20), we have derived the expression of the related conditional Laplace transforms in the case that the serving BS is a long queue length BS. Now, we consider the case that the serving BS is a shortqueue BS. Comparing (12) with (9), the interfering longqueue BSs is almost the same but adding a typical cluster x 0 ∈ l c of the BSs. Using Slivnyak's theorem [5], we have L I s-long (s | r) = L I l-inter (s | r).
Similarly, comparing (13) with (10), the interfering shortqueue BSs is almost the same but removing a point y 0 from PPP s , therefore, when the serving BS is a short-queue BS, the conditional Laplace transform of the short-queue BSs interference is Theorem 1: The coverage probability of the multi-cell downlink system with typical BS z ∈ b is given by The result follows the law of total probability for the typical BS.

V. NUMERICAL RESULTS
In this section, we provide Monte Carlo simulation results to verify the accuracy of our analysis and compare the performance of the proposed method with that proposed in [18], which serves as a benchmark. Simulation parameters are based on the 3GPP standard [27]. In the simulation, we consider a square study area with 1000 m 2 . Following the definition of dense networks [28], we deploy the BSs and the devices in the study area as independent PPPs with density λ b = 0.1 BS/m 2 . We consider both large-scale fading and small-scale fading. The path loss exponent is set as α = 3.8. The small scale fading of all links are i.i.d. and follow complex Gaussian distribution with zero mean and unit variance. In the MAC layer, the packets arriving at each node follow the independent Bernoulli process with parameter ξ . The maximum system operating time is set as T = 1000 time slots. The iteration tolerance in Algorithm 1 is set by = 10 −5 . For the sake of clarity, the parameter settings are listed in Table 2.
By iteratively running the system to the maximum operating time, the long-queue and short-queue BSs are extracted by adopting the K-means method [29]. To characterize the pattern of the spatial point process, we introduce Ripley's K function [30]. Given a randomly chosen event and a distance t, Ripley's K function K (t) is defined as the expected number of the events (except the chosen event) within the radius t circle centered at the chosen event divided by the number of the events in the unit area. Mathematically, the Ripley's K functions of NSP and PPP are respectively given as follows:  Figure 8.

TABLE 4.
Fitting parameters for the cases in Figure 9.
However, for the events with unknown point processes, the closed-form expression of Ripley's K functions is unavailable. To address this issue, based on the definition in [30], the estimated Ripley's K function is given bŷ whereλ = N A is the estimated density of the target point process with N equals to the observed number of the points in the study area A; 1 + (·) is the indicator function (same as the one in (15)); d i,j is the distance between point i and point j; and w i,j is the weight that related to point i and point j. Suppose that there is a circle centered at point i and passing through point j, w i,j equals to the portion of the circumference of that circle that falls in the study area. In our system, suppose N l and N s is the number of the longqueue BSs and the short-queue BSs in the study area, respectively. Replacing parameters (N , ) in (35) with (N l , l b ) and (N s , s b ), we obtain the estimated Ripley's K function of the long-queue BSs and the short-queue BSs i.e.,K l (t) andK s (t), respectively.
In Monte Carlo simulations, the coverage probability of serving BS z is estimated bŷ where n tran is the transmission number at the transmitter and n recv is the number of the packets that are successfully received at the intended receiver. The result is obtained by averaging over 1000 independent Monte Carlo simulations. Parameters λ p ,m, σ 2 of NSP for the long-queue BSs and λ s b of PPP for the short-queue BSs in Algorithm 1 control the portion of the long-queue and the short-queue BSs as well as the patterns of the two point processes. Similar to the idea in [31], the coverage probability of (32) is obtained by carefully choosing the aforementioned parameters, i.e., λ p ,m, σ 2 . Table 3 and Table 4 give the fitting results of the parameters in different cases. VOLUME 8, 2020 In addition, a comparative result in [18], named as ''IID queue approximation'', is also presented with the assumption of temporal related but spatial independent queues.
To present the spectral efficiency (SE) that the system can provide for the devices, we follow [32] and define the averaged potential user SE by We remark that the potential SIR of the successfully transmitted device is regarded as constant θ. Therefore, (37) is not the accurate expression but the lower bound of SE for the target device.

A. SIMILARITY OF THE APPROXIMATED MODEL
To illustrate the rationality of our approximated model, we first present the similarity between the point processes in the real deployed simulation and the approximated standard point processes. Figure 4 depicts the Ripley's K functions of the standard NSP and PPP as well as the estimated Ripley's K functions of the long-queue BSs and the short-queue BSs with θ = 10 dB and ξ = 0.05 packets/slot. In this figure, the dashed lines are the corresponding values ofK l (t) andK s (t), whereas the solid lines are the theoretical values K NSP (t) and K PPP (t) in (33) and (34), respectively. Solid line K NSP (t) (resp. K PPP ) is enveloped by two dotted lines, which are 0.025 and 0.975 quantiles of K NSP (t) (resp. K PPP ) estimated over 1000 simulations, respectively. As seen in Figure 4, the Ripley's K functions for NSP and PPP are apparently different especially when the distance becomes large. Comparing the dashed lines with their corresponded solid lines, there is a close match between the simulated locations of the long-queue BSs (respectively, short-queue BSs) and the NSP (respectively, PPP), which validates our approximation approach. Figure 5 shows the CDFs of the conditional SIR coverage probability F y,n obtained by the simulations and the CDFs of the approximated Beta distribution over different decoding thresholds. From Figure 5, we observe that given a value of CDF, the conditional SIR coverage probability decreases as the decoding threshold θ increases. the conditional SIR coverage probability decreases as the decoding threshold θ increases. For instance, when θ = 5 dB, 50% percentage of the devices obtain 0.9 coverage probability, whereas when θ = 15 dB, more than half of the devices cannot be successfully served anymore. The reason comes from the fact that higher decoding threshold increases the transmission failure probability of the packets. Comparing the two distributions, the CDFs of the Beta distribution match well with the one calculated by the simulation results under different decoding orders, which confirms the approximation method of the conditional SIR coverage probability. Figure 6 examines the effect of the density of parent point process λ p and the portion of long-queue BSsm λ p λ b on the value of coverage probability in (32). As seen, with the increase of fitting parameters λ p andm λ p λ b , the coverage probability keeps decreasing. This is due to the fact that the more long-queue BSs, the more interference suffered by the system, which results in a lower coverage probability. By comparing the two fitting parameters, we find that the portion of the long-queue BSs has a more significant impact on the coverage probability. For instance, fixing λ p = 0.005, the coverage probability of the serving BS gets more than 62.50% degradation (from 0.40 to 0.15) as the expected number of the BSs in each cluster improves 5 times (m λ p λ b changing from 0.1 to 0.5). By contrast, with the fixed number of the long-queue BSs (i.e., mλ p λ b = 0.5), the performance reduction is only 28.57% (from 0.21 to 0.15) when λ p gets 5 times increment from 0.001 to 0.005. Figure 7 depicts the effect of the density of parent point process λ p and the variance of the offspring point process σ 2 in NSP. Similar to the results in Figure 6, the coverage probability decreases with the increase of λ p and σ 2 . However, different from Figure 6, both of the two fitting parameters do not impact the performance of the coverage probability  apparently. Specifically, among all simulated values of σ 2 , the largest coverage probability reduction occurs to be 7.04% when λ p = 0.002. Similarly, with λ p increasing from 0.001 to 0.005, the maximum decrement of the coverage probability is 18.04% with σ 2 = 2. Figure 6 and Figure 7 demonstrate that the coverage probability is more sensitive to the portion of the long-queue BSs rather than λ p and σ 2 . This implies that to obtain a better approximated performance via Algorithm 1, the portion of the long-queue BSs plays an important role and should be selected carefully. Figure 8 depicts the coverage probability versus decoding threshold with ξ = 0.05 packets/slot. We compare the results of Simulation, IID queue approximation and Theorem 1. Specifically, the fitting parameters for Theorem 1 are given in Table 3. As observed from the figure, the coverage probabilities of the three methods decrease as decoding threshold θ increases. Moreover, Theorem 1 approaches Simulation well under all the considered decoding orders. The gap of the coverage probability between the two methods is always less than 0.05. The performance of IID queue approximation is better than Theorem 1 when θ ≤ 4 dB. This is due to the fact that when decoding threshold θ is small, the packets can be successfully transmitted with a relatively high probability and the buffers of the BSs are more likely to be empty. Therefore, the BSs with empty buffers shut down the transmissions and those active BSs are nearly spatially independent. However, Figure 8 shows that there is a noticeable gap between IID queue approximation method and Simulation when the decoding threshold θ becomes larger than 6 dB. Specifically, with θ = 12 dB, the coverage probability of IID queue approximation method is 1.8 times larger than Simulation. The gap between the two methods mainly attributes to the spatial correlation among the queueing states of the BSs located in geographical proximity, which leads to a prolonged transmission duration and raises the interference across the entire network. Overall, by incorporating an underlying clustered structure to the interfering point process, our method successfully addresses the issue caused by the spatial queueing interaction and gives an accurate characterization to the probability of coverage. Figure 9 compares the performance of the three methods with different packet arrival rates. To guarantee a stationary queue for the system, we set θ = 5 dB and packet arrival rate ξ changing from 0.1 packets/slot to 0.3 packets/slot based on Table 1. The values of the fitting parameters are given in Table 4. We observe that the coverage probability of all three methods decrease with the increase of the packet arrival rate. Specifically, when ξ increases from 0.1 packets/slot to 0.3 packets/slot, the coverage probabilities of Simulation, IID queue approximation and Theorem 1 decrease from 0.69 to 0.39, from 0.76 to 0.56 and from 0.69 to 0.36, respectively. Moreover, Figure 9 shows that the coverage probability of IID queue approximation is always larger than Simulation. Furthermore, the gap of the two curves continuous increasing from 10.17% to 45.47% as the packet arrival rate increases from 0.1 packets/slot to 0.3 packets/slot. The reason is that when ξ goes up, the increased traffic load results in the phenomenon that the interfering BSs are more and more close to the serving BS. Therefore, the impact of the spatial correlation between the BSs becomes more and more apparent and the high interferences among the BSs reduce the performance of the coverage probability. Among the three methods, Theorem 1 approaches the results of Simulation much better than IID queue approximation. The performance gap by applying Theorem 1 is always less than 5.63%, which manifests the validation of our algorithm. Table 3 and Table 4 show the value of the fitting parameters for the cases in Figure 8 and Figure 9, respectively. The basic idea of such settings corresponding to the intuition that the denser of the BSs, the more portion of the long-queue BSs mλ p λ b in the system, and hence results in a larger value of λ p and σ 2 . For example in Table 4, with ξ improving from 0.1 packets/slot to 0.3 packets/slot, the value of λ p ,m λ p λ b and σ 2 increases from 1e−3 to 5e−3, from 0.02 to 0.15 and from 1 to 3, respectively. Besides, based on the observations in Figure 6 and Figure 7, the coverage probability is more sensitive to the value ofm λ p λ b than λ p and σ 2 . For instance, in Table 3, when θ improves from 0 dB to 16 dB, the value ofm λ p λ b enlarges more than ten times relative to the original value 0.02. Figure 10 displays the convergence of active probability of the BSs in Algorithm 1 over different decoding thresholds, where the packet arrival rate is fixed as ξ = 0.05 packets/slot. From the figure, one can observe that the sequences of active probability {q y,n } n of all the three cases monotonically increase and converge fast. Specifically, the value of q y,n gradually increases to 0.18, 0.16 and 0.11 when θ equals to 15 dB, 10 dB and 5 dB, respectively and sequence {q y,n } n of all the three cases converges within five iterations. This phenomenon indicates that Algorithm 1 runs with a low computing burden. Comparing among the three curves, the larger decoding threshold θ is, the larger convergence value of sequence {q y,n } n can be obtained. This is due to the fact that given a large decoding threshold, the probability of failing packet transmission is relatively high, which prolongs the queue lengths and the active probability of the BSs. Figure 11 presents the averaged potential user SE obtained by (37) over different decoding thresholds. Three packet arrival rates ξ = 0.01, ξ = 0.03 and ξ = 0.05 packets/slot are considered. Different from the curves of coverage probability in Figure 8, the tendency of the averaged potential user SE changes with different values of ξ . For instance, when ξ = 0.01 packets/slot, the averaged potential user SE monotonically increases with the ascending of θ. While in case ξ = 0.05 packets/slot, the curve first increases then decreases when θ changing from 0 dB to 16 dB. This is due to the reason that potential SE is determined by the multiply of potential rate log 2 (1 + θ) and P(γ > θ). Although P(γ > θ) is the lowest when θ = 16 dB, the large value of θ would compensate the loss of coverage probability and results in a high averaged potential user SE. Even though, we would like to mention that this kind of high averaged potential user SE is obtained with the cost of the un-fairness among the devices. That is, most of the devices are failed to transmit the signals, while only a small portion of the devices are successfully transmitted with a high user SE.

VI. CONCLUSION
We conducted a statistical study on the coverage probability of the multi-cell downlink system with the dynamic traffic pattern and spatially interacting queues. In terms of the queueing model, the randomness in the spatial and temporal domain is generally homogeneous. We first gave the sufficient and necessary conditions of the ε-stable region. Then, taking into account the spatial correlation among the queues, we extracted a two-queue-length non-homogeneous spatial structure to model the original system. By respectively approximating the locations of long-queue BSs and shortqueue BSs as NSP and PPP, we calculated the coverage probability of the system via an iterated algorithm. To further reduce the computational complexity, Beta distribution was introduced to approximate the PDF of the service rate in each iteration of the algorithm. Simulation results revealed that the proposed algorithm converges fast and the analytical result is sensitive to the portion of long-queue BSs. Numerical results also demonstrated the accuracy of the analysis and revealed that the proposed scheme offers significant performance gain compared to a benchmark result. In the future, one promising direction is to extend the current work to other dense scenarios, e.g., multi-user case, ultra-dense network, dense optical communication system, etc..

APPENDIXES APPENDIX A PROOF OF LEMMA 2
The conditional serving distance distribution of the typical device is given by The conditional Laplace transform of the intra-cluster interference at the typical device is where (a) follows from the law of total probability over a infinite space; B c (0, r) is the complementary area of B(0, r), f W (w | ·) is the conditional PDF of w = x 0 +y; (b) follows from the Poisson distribution of the BS number k and the averaged number of the BSsm in cluster l x 0 (one serving BS andm−1 interfering BSs); (c) follows from equality e x = ∞ k=0 x k k! , B c (0,r) f W (w | v 0 , r) dy = 1, and from changing variable x 0 + y into w out . Define˜ l x to be the set of the BSs that in both set l x and circle area B(0, R), the conditional Laplace transform of intercluster interference at the typical device is where u = x + y and v = x ; the terms of (a) and (b) in the first product of (39) follow a same idea of the induction (a)-(c) in (38); equation (c) is obtained by using the probability generating functional (PGFL) of PPP for parent set {x ∈ l c \x 0 }.

APPENDIX C PROOF OF LEMMA 4
When the serving BS is a long queue length BS, the conditional Laplace transform of the intra-cluster interference at the typical device is where v = y ; (a) follows the exponential distribution of h y ; (b) follows the law of total probability of ζ y ; (c) follows the PGFL of PPP for short-queue BSs set {y ∈ s }. He is currently the Cheng Tsang Man Chair Professor with the Singapore University of Technology and Design (SUTD). He is also the Head of ISTD Pillar, a Sector Lead of the SUTD AI Program, and the Deputy Director of the SUTD-ZJU IDEA. He has been actively involved in organizing and chairing sessions. His current research interests include wireless communications and networking, network intelligence, the Internet-of-Things, URLLC, and big data processing.
Dr. Quek has served as a member of the Technical Program Committee as well as symposium chairs in a number of international conferences. He is an He is currently a Professor with Beijing Jiaotong University, where he is also a Chief Scientist with the State Key Laboratory of Rail Traffic Control and Safety. His interests include wireless communications for railways, control theory and techniques for railways, and global system for mobile communications for railways systems. He received the Mao YiSheng Scientific Award of China, the Zhan TianYou Railway Honorary Award of China, and the Top ten Science/Technology Achievements Award of Chinese Universities. VOLUME 8, 2020