Modeling D2D Underlaying Cellular Network for Hotspot Communications with Poisson Cluster and Hole Processes

This paper develops a new approach to the modeling and analysis of device-to-device (D2D) underlaying multi-tier cellular network for the dense hotspot communications, which consist of macro base stations (MBSs), pico BSs (PBSs), femto BSs (FBSs). A typicl user equipment (UE) can work either in D2D mode or cellular mode. Considering the dense hotspot communications, this work employs Poisson point process (PPP) to model the locations of MBSs and PBSs, and uses Poisson cluster process (PCP) to model the ones of UEs and FBSs. The locations of PBSs are also modeled as the centers of hotspots, referred to as the centers of PCPs. UEs and FBSs cluster around the common parent process PBSs. To guard the cluster-edge UEs, the clustered-UE classiﬁcation and modiﬁed fractional frequency reuse (FFR) are jointly used, by which both the UEs and FBSs are classiﬁed two sets, cluster-center UEs and cluster-edge UEs, cluster-center FBSs and cluster-edge FBSs, respectively. The total frequency band is divided into two orthogonal segments, one of which is shared by D2D devices, cluster-edge FBSs, and PBSs, and the other segment is shared by cluster-center FBSs and MBSs. For such clustered multi-tier network, by using the methods from PPP, PCP, and PHP, this paper presents a tractable approach for modeling and analyzing the performance of cellular and D2D networks and gives the statistical descriptions of the experienced interference at a typical receiver by using the approximated Poisson hole processes (PCP) theory. This yields the derivations of the coverage probabilities of both the D2D receivers and cellular destinations. In additon, during the analysis of cellular UEs, to derive the coverage probabilities, this paper specially constructs a UE association criterion as well as the derivations of both the association probabilities and the statistical descriptions of association distances for cluster-center and cluster-edge UEs. The simulations results exploit the effect of various network parameters on the network performance and give the insights in terms of the proposed schemes as well as the comparison between cluster-center and cluster-edge UEs.


Introduction
In the past few years, due to the proliferation of high-speed multi-media applications and high-end devices such as smartphone, tablets, wearable devices, laptops, machine-to-machine communication devices, etc., the mobile data traffic has been growing exponentially to satisfy custom's huge and different link rate demand in different network services. Moreover, in future Internet-enabled high-end devices, it is expected to generate more mobile data traffic. All of these are strikingly transforming the traditional cellular network from the coverage-derived deployment of macrocells to a more user-centric capacity-derived deployment of several types of base stations (BSs), collectively called small cells. This resulting network architecture consists of one or more small cell tiers overlaid on the traditional macrocellular tier and is referred to as a multi-tier heterogeneous network (HetNet) [1]- [2]. Macro BSs (MBSs) having large coverage are usually used to control user equipment (UE) scheduling, resource allocation, and high-mobile UEs, whereas small BSs (SBSs) serve as the fundamental element for traffic offloading from MBSs, therefore act as the main capacity-driver for indoor and outdoor low-mobile UEs. As a result, SBSs are usually deployed in hotpots to serve tens of active users by offloading their traffic from MBSs. In small cells, the transmission power of pico BSs (PBSs) is typically up to 35dBm. However, because femtocells are user-deployed indoor BSs, the transmission power of femto BSs (FBSs) typically is less than 20dBm and the coverage range of FBSs is in the order of tens of meters [3]. However, in such densely deployed HetNet, the increase of SBSs puts a heavy burden on backhaul links which connect BSs with the core network, even causes network congestion. Existing works focus on offloading the data-hungry or computing intensive tasks to the cloud for execution. This mehtod will yield a long latency when mobile devices are connected to the cloud relayed by BS, which does not satisfy the requirements of delay sensitive tasks.
In this regard, device-to-device (D2D) based cellular traffic offloading architecture is an efficient and essential technique in the fifth generation (5G) and beyond 5G (5G/B5G) network [4]. In D2D offloading communications, mobile users download the cellular traffic contents from their nearby mobile devices that hold the same contents using D2D communications. D2D offloading communication allows proximity UEs to establish a direct communication link with each other by bypassing the BS so that the conventional two-hop cellular link is replaced by a direct D2D link. This short range communication has a number of advantages, e.g., high data rate, high energy and spectral efficiency, increased network capacity, low delay, improvement of the performance in cell edges [5]. D2D communication can be categorized into out-band and in-band D2Ds based on the used spectrum [6]. In the case of out-band mode, D2D transmitter (D2D-Tx) exploits the unlicensed spectrum. In the second case, D2D uses the licensed spectrum to communicate. In practice, D2D offloading communication should fully take into account the requirements and (or) constraints of both operators and users. From the operator standpoint, the more offloaded traffic, the better for cellular networks to avoid the network congestion. From the perspective of users, the users always expect their quality-of-service (QoS) requirements to be satisfied, which is also concerned by operators, since this directly yields profits to the operators. In D2D underlaring cellular network, the various resulting co-channel interference is a serious challenge and affects users' QoS. In particular, the UEs locating in cell edge area usually receive weaker signal strength from its serving MBS and stronger interference from co-channel MBSs. Moreover, the cell edge users (CEUs) suffer from the massive inter-tier interference from the deployment of small cells (such as pico and femto cells). As a result, in such HetNets, these inter-tier and intra-tier interferences yield that CEUs experience reduced coverage as well as the degradation of transmission rate. Correspondingly, to satisfy the QoS, the more bandwidth should be allocated to D2D users or cellular users. This also results in low spectrum efficiency. In fact, without proper interference management and effective resource allocation, directly deploying small cells and D2D network may worsen network performance.
To mitigate this problem and improve the performance of network, different techniques of spectrum management have been proposed in literature. The first method is fractional frequency reuse (FFR). In FFR, the total available bandwidth is partitioned into two orthogonal sets. Users in cell-center are allocated a common sub-band of frequencies while cell-edge users bandwidth is partitioned across cells based on a reuse factor. In total, this strict FFR scheme requires a total of sub-bands. UEs in cell-center do not share any spectrum with cell-edge UEs. Obviously, this strict FFR can reduce interference for both cell-center and cell-edge UEs and has been investigated widely. Specially, the works [7] and [8] focused on single-tier cellular network, and the works [9] and [10] extended to multi-tier networks. Different from strict FFR, the soft frequency reuse (SFR) allows a full frequency band reused in each cell, which is a variation of FFR. In general, SFR employs the same cell-edge bandwidth partitioning strategy as strict FFR, but the cell-center UEs are allowed to share sub-bands with cell-edge UEs locating in other cells. Because cell-center UEs share the bandwidth with neighboring cells, a power control is required . That is to say, the cell-center UEs transmit at lower power levels than the cell-edge UEs [11]. While SFR is more bandwidth efficient than strict FFR, it results in more interference to both cell-center and cell-edge UEs [12][13][14][15].
Another key technology of spectrum management to mitigate the inter-tier interference is shared spectrum allocation (SSA) [16]. As stated in [16], the SSA method is composed of orthogonal sharing and co-channel sharing methods. Orthogonal sharing divides the licensed spectrum into two disjoint fractions, i.e. one for macro cell layer and the other for small cell tier, thus eliminating the interference between layers (cross-layer interference). However, if we consider a mature macro cell network, there may be an insufficient number of sub-bands available for small cells. Orthogonal sharing can be inefficient in terms of spectrum reuse. Considering the scarcity of radio resource and the easy of deployment, co-channel sharing would be more efficient and preferable for operators due to increased spectral efficiency through spectral frequency reuse. However, co-channel sharing between layers is technology challenging since the capacity is limited by cross-layer interference. As a compromise between cross-layer interference and efficient spectrum reuse, partial co-channel sharing between layers has been considered. Partial co-channel sharing allows the macro cell layer to occupy the entire spectrum, while the small cell layers only share portion of the spectrum. In partial co-channel sharing, if the small cells access most of the macro cell spectrum, both unnecessary co-channel interference to the macro cell and poor spectrum efficiency may result. Contrarily, a cosing only a few spectrum could result in excessive small cell interference in the sharing spectrum. Therefore, the works [17][18][19] considered such SSA scheme.
In addition, for the performance evaluation, the researchers mostly relied on highly approximated models like Wyner model and grid-based hexagonal model. However, existing works showed that due to irregularity and randomness exhibiting by the BS locations under real deployments in multi-tier HetNets, stochastic geometry and point process spatial models have been identified as elegant tools for accurate modeling and tractable analysis of the systemlevel performance of heterogeneous cellular networks [20,21]. In stochastic geometry-based approach, the spatial locations of the large-scale wireless network are abstracted to convenient point processes so that the system-level energy-spectrum efficiency and coverage performance can be evaluated. A strong motivation for viewing the BS arrangement as a homogeneous point process can be drawn from the study of the works [22][23][24], which suggested that significant insights can be admired by bounding the downlink (DL) performance between the idea hexagonal grid mode and homogeneous Poisson point processes (PPP) based mode. Specially, Elsawy et al provided a comprehensive tutorial on stochastic geometry-based analysis for cellular networks [25]. Therefore, inspired by the above considerations, the works [9,14,15,[26][27][28][29] evaluated the performance of FFR and SFR based cellular networks by using stochastic geometry and point process based model. However, these works only modeled the locations of BSs and UEs in each tier as the general independent PPPs. The assumption of independence between BSs and UEs may not be quite accurate in some cases/UE-content based 5G HetNets since there exists correlation between UEs and BSs, especially SBS [19]. In such architecture, one can envision the small cells being deployed to serve cluster of UEs by offloading from MBSs. Such models have been used by some standardization bodies, such as third Generation Partnership Project (3GPP). 3GPP has considered the clustered configurations in which locations of the users and small-cell BSs are coupled, in addition to the uniformly distributed UEs [19]. The authors in [30] have proved that the model in which UEs are distributed according to Poisson cluster process (PCP) around PPP distributed SBSs closely resemblances the 3GPP configuration of single SBS per user hotspot in a HetNets. Therefore, user-centric deployment of small cells is becoming an important part of future capacity-derived wireless architectures [31]. Though using PCP exploits the correlation between UEs and SBSs, it could be an unrealistic choice for modeling the scenario in which locations of active nodes exhibit spatial separation due to interference management. In certain case, to manage interference, the UEs are categorized into cell-center and cell-edge users. This leads to the locations of cell-center or cell-edge user partial separation. Obviously, for such scenarios where locations of active nodes exhibit spatial separation, it is an unrealistic choice to model such scenario by using PPP. In this case, a more appropriate model is Poisson Hole process. As stated in [31][32][33][34][35][36], PHP is formed by carrying out holes from a baseline PPP, where the centers of the holes are assumed to follow an independent PPP. The remaining points of the baseline PPP are said to form a PHP.
The distances of a typical D2D-Rx to its serving cluster-center and cluster-edge D2D-Txs The association distances of a cluster-center UE with FBS, MBS,and PBS The association distances of a cluster-edge UE with FBS and MBS SIR CPs of cluster-center and cluster-edge D2D-Rxs, total CP of a D2D-Rx The PDF and CDF of random variable X Therefore, motivated by the above considerations and the huge capacity of D2D offloading, by jointly employing PPP, PCP, and PHP, this paper presents a tractable framework for modeling and analyzing the hotspot communication scenario of a three-tier D2D underlaying heterogeneous network which consists of MBSs, PBSs, and FBSs. A typical UE can work either in cellular mode, called cellular UE, or in D2D mode, called D2D UE. The D2D communications are responsible for offloading cellular network. The locations of MBSs and PBSs are modeled as independent PPPs.The PCP and PHP are used to model the locations of FBSs and UEs. Specially, in this paper, the Thomas cluster process (TCP) is used. On one hand, the PPP distributed PBSs are assumed to be the common center of clustered FBSs and UEs so that the locations of FBSs and UEs are modeled as PCPs of parent process consisting PBSs. On the other hand, with a given coverage radius R 2 of PBSs, the UEs scattered around each PBS are classified as cluster-center UEs and cluster-edge UEs. The classification yields that the cluster-center UEs and cluster-edge UEs are modeled by PCP and PHP. To restrict interference, the SSA-based interference management is provided. In SSA, the total available bandwidth is portioned into two orthogonal sub-bands. The cluster-center sub-band is shared by PBSs and all D2D users, but the cluster-edge sub-band is shared by MBSs. As a result, the interference experienced at D2D receivers (D2D-Rxs) is restricted and the offloading capacity of D2D users is improved. At the same time, to enhance the spectrum efficiency, inspired by the essential low-power of FBSs, the two sub-bands are shared by cluster-center and cluster-edge FBSs, respectively. For such network architecture with D2D offloading, by resorting to the method from PCP and PHP, the interference distributions at each type of receiver are obtained. Finally, the coverage performance is calculated. The presented simulations validate the derivations, but also show a clear insight in terms of the effect of the key network parameters. It is found that the presented method improves D2D users, but also cellular users. For a quick reference, the frequently used notations are summarized in Table 1. Unless otherwise specified, we also use b (0, a) denoting a circle of radius a and B (a1, a2) denoting an annulus of radii a1 and a2.
2 Network mode and assumptions Fig.1 shows a three-tier heterogeneous cellular network with underlaying D2D communications, which consists of MBSs, PBSs, and FBSs. The whole network operates over sub-6GHz frequency band with total available bandwidth W . MBSs are equipped with massive MIMO of N M antennas and PBSs are equipped with traditional MIMO of N P antennas. All FBSs and UEs are equipped with single antenna having low emplement complexity. MBSs and PBSs employ the simple linear zero-forcing beamforming (ZFBF) to communicate simultaneously with S M and S P UEs, respectively. One UE can work in either D2D communication mode or cellular network communication mode. Specially, a D2D link is established only when the D2D-Rx falls inside the coverage range of D2D-Tx and requested content can be found in its counterpart. For simplicity, throughout this work it is assumed that a UE contains the requested content of its counterpart with probability p and a D2D UE is a D2D-Tx with probability q. The transmit powers of MBSs, PBSs, FBSs, and D2D-Txs are P M , P P , P F , and P D , respectively. The coverage radii of MBSs, PBSs, FBSs, and D2D-Txs are R 1 , R 2 , R 3 , and R 4 . Before going to more technical details, it is assumed that a typical UE is randomly selected. According to the Slivnyak-Moche's celebrated theorem, the fact that a point process on a Borel space is Poisson if and only if its reduced Palm distributions agree almost everywhere with the original distribution theorem. Therefore, it is reasonable that the analysis is conducted on the typical UE located at origin.

Network model and distribution of distances
The spatial locations of MBSs and PBSs are modeled as independent Poisson point processes (PPPs) Φ M and Φ P with density λ M and λ P on the Euclidean plane, respectively. By letting x M denote the distance from a typical UE to the nearest MBS, x M ∈ Φ M , the PDF of the distance x M is given by f where . is the Euclidean norm. Further, with a more realistic scenario, this work considers that pico cells (PBSs) are expected to be deployed in crowed areas or hotspots to patch coverage dead zones. Without loss of generality, it is assumed that the densely deployed PBSs are the centers of hotspots. The locations of FBSs are modeled by TCP Θ F = ∪ x∈Φ P N x of spatial density λ F with parent process Φ P , where the offspring N x generated for a given parent x ∈ Φ P [37][38][39]. In each cluster, the offspring N x are scattered around the parent point x ∈ Φ P according to a symmetric  [37]. At the same time, to improve the performance of UEs locating inside cell edge, we employ a PHP to model the points locating inside the edge of cluster [35,36]. Specially, with a given coverage radius R 2 of PBSs in the parent process Φ P , the set of the points of Θ F falling inside coverage region of PBSs is modeled by where b(y, D) is a circle (ball) of radius D centered at y. Then, the set of the remaining points of FBSs is modeled by a PHP , called cluster-edge FBSs. Straightforwardly, the locations of UEs follow arbitrary independent point process Θ U of spatial density λ U . To capture the natural coupling property, UEs are independently scattered around parent point process Φ P according to TCP having the points' average number c D and the variance σ 2 D of distance distribution. The points set of Θ U falling inside the coverage region of PBSs is modeled by , called cluster-center UEs, and remaining points of UEs is modeled by a PHP For FBSs, a FBS is chosen randomly in the representative cluster, of which cluster center located at x P0 ∈ Φ P . Then, the distance from a typical UE to the FBS in the representative cluster is written as x R F . In general, the statistical description of the distance x R F greatly depends on the common representative cluster center x P0 ∈ Φ P . That being said, if we condition on the location of the common representative cluster center x P0 ∈ Φ P , the distance x R F can be characterized by Rician distribution [40], i.e., f x R , where I 0 (.) is the modified Bessel function of the first kind with order zero and the variance σ 2 F of distances. However, the results in [41] and [42] indicated that the correlation among the distances in the same cluster is sufficiently weak and can be ignored. More specially, since the distances from the representative cluster center x P0 ∈ Φ P to the typical FBS and UE are independent and non-identically distributed Gaussian random variables with variances σ 2 F and σ 2 D , respectively, the resulting summarized distance x R F of the typical UE to FBS is also Gaussian with variance σ 2 F +σ 2 D . Therefore, we can use Rayleigh fading of variance σ 2 F +σ 2 D to approximate the distribution of x R F . That is to say, the PDF of . The CDF is given by . However, in this work we consider the scenario where a typical UE could be associated with the nearest FBS in the same or representative cluster and denote the corresponding distance as x F . In this case, the CDF of the distance x F from the typical UE to the nearest FBS is given by F . By taking the derivative of F x F (x), the PDF of the nearest distance x F is given by It is further assumed that all wireless channels in this D2D underlaying network are independent and are subject to both path loss and small-scale fading. All small-scale fading channel coefficients are modeled as zero-mean complex Gaussian random variables. Thus, by assuming that the typical UE is associated with one arbitrary access point at y n , then we have that y n = y n is the distance of the typical UE to the access point. The corresponding path loss is modeled as β y n −α , where α is the path loss exponent with a typical value 2 ∼ 6 and β is a frequency dependent constant value. For the channel with the single-antenna transmitter, the small-scale fading channel gain is expressed as h y n and follows Rayleigh fading with unit power, i.e., h y n ∼ exp(1), and for the channel from the MIMO transmitter, the equivalent small-scale fading channel gain is expressed as g y n and obeys Gamma distribution, i.e., g y n ∼ Γ(., .). Because of the severe interference from other transmitters, an interference-limited HetNet is considered so that the additive white noise is ignored.

Separated spectrum allocation of classified UEs
As depicted in Fig.1, it is easy to see that a typical UE can work in D2D mode or cellular mode. Under cellular mode, a typical UE can also be associated with MBS, PBS, or FBS. Moreover, with a given association, we have that a typical UE that locates inside cluster is referred to as cluster-center user, and a typical UE that locates outside cluster is referred to as cluster-edge UE. All of these cases result in the complicated types of UEs. Firstly, for the set of D2D transmitters, we have the cluster-center D2D-Tx Φ CDT and cluster-edge D2D-Tx Φ EDT processes. For the set of D2D-Rxs, we have cluster-center D2D-Rx Φ CDR and cluster-edge D2D-Rx Φ EDR . Secondly, we have cluster-center FBS process Φ CFB and cluster-edge FBS processes Φ EFT . Thirdly, for a typical UE that is associated with MBS and FBS, we have the cluster-center macro (femto) UE process Φ CM(F)U and cluster-edge macro (femto) UE process Φ EM(F)U , respectively. Because of the fact that when the typical UE is associated with PBS, the typical PUE always locates inside the coverage of cluster. We have the cluster-center pico UE process Φ CPU . With the given UE classification, we propose a SSA strategy that can preserve from additional interference at UEs. In SSA strategy, the total available total bandwidth W is partioned into two orthogonal sub-bands W 1 and W 2 , respectively. By definding a bandwidth portion factor ξ , it is achieved that W 1 = ξW and W 2 = (1 − ξ )W . With this sub-band partion, we propose the spectrum access scheme as shown in Fig.2. We have that the sub-band W 1 is allocated to all PBSs, cluster-edge FBSs, and all D2D users.The sub-band W 2 is allocated to MBSs and cluster-center FBSs. Because the orthogonal spectrum sharing is used, the interference between the macro cluster-edge UEs and pico cluster-center UEs is eliminated effectively. Besides this orthogonal spectrum sharing, it is easy to find that due to the inverse spectrum allocation between PBS (MBS) and FBSs, the considered co-channel sharing scheme is able to restrict inter-tier interference. Moreover, the cluster-center and cluster-edge FBSs adopt the orthogonal so that the cotier interference is restricted. Finally, because of D2D-Txs having low transmission power, the proposed SSA scheme allows all D2D-Txs to share the sub-band W 1 . That is to say, the D2D-Txs adopt partial co-channel sharing.

Coverage performance of D2D network
We first perform the detailed analyses with respect to the performance of D2D network, in which a typical UE is associated with another UE that has the required content of the typical UE with D2D mode probability p and falls inside the coverage of its counterpart. To this end, this section first provides the laplace transforms (LTs) of interferences experienced at a typical D2D-Rx by considering the proposed cluster-based UE classification and spectrum sharing allocation. Then, we would derive the SIR coverage probability. Note that, considering the propagation loss of signals, only the dominant interferences that practically contribute the interference at a typical UE are considered, i.e., out-ofrange transmitters to a typical UE are excluded.

Cluster-center D2D-Rxs
Without loss of generality, let the typical D2D-Rx associated with the D2D-Tx locate at x C D . Observing Fig.2, it is achieved that all D2D-Txs share the sub-band W 1 with PBSs and cluster-edge FBSs.This observation yields that the received interference I CDU at the typical cluster-center D2D-Rx is formulated as where I CDT CDU is the interference from cluster-center D2D-Txs, I EDT CDU is the one from cluster-edge D2D-Txs, I P CDU is the one from PBSs, and I EF CDU is the one from cluster-edge FBSs. In I CDT CDU , considering the fact the typical UE receiver is associated with a cluster-center UE transmitter, we have that the possible cluster centers would be located in the region x P ∈ Φ P ∩ b (0, R 2 + R 4 ) and the interference is from all cluster-center D2D-Txs locating at B where R 2 and R 4 represent the coverage radii of PBSs and D2D-Txs, respectively, and B x P stands for the set of D2D-Txs inside the cluster centered at x P . Hence, it is easy to mathematically achieve the formulation of I CDT CDU as In (2), the interference from cluster-edge D2D-Txs is formulated as where Φ EDT represents the set of cluster-edge D2D-Txs. The interference from PBSs is written as where g x P is the small-scale fading interfering channel power gain, which follows gamma distribution. The interference from cluster-edge FBSs is formulated as where Φ EFT is the set of cluster-edge FBSs. Comparing (4) and (6), it is not difficult to find the interferences I EF CDU and I EDT CDU having the completely similar forms. Then, combining the spectrum allocation scheme given by Fig.2, we can characterize the aggregate interference I CDU at the typical cluster-center D2D-Rx and have Lemma 1.

Lemma 1:
The LT of the interference I CDU experienced at the typical cluster-center D2D-Rx is given by where (s), and L I EF CDU (s) are given by Proof: See Appendix A

Cluster-edge D2D-Rxs
Fig .2 shows that cluster-edge D2D-Rxs also share sub-band W 1 with PBSs and cluster-edge FBSs. Then, the interference at the typical cluster-edge D2D-Rx locating at x E D is written as In (12), I EDT EDU is the interference from cluster-edge D2D-Txs except the one connected with the typical D2D-Rx lo- is the interference from clustercenter D2D-Txs at the typical UE whose access point D2D-Tx falls cluster-edge area. As stated previously, the typical D2D-Rx can fall inside cluster. This yields that in I CDT EDU , all possible interfering cluster centers would locate in B (R 2 − R 4 , R 2 + R 4 ). Then, by taking the cluster center x P ∈B, we further define B x p denoting the set of interfering D2D-Txs in the cluster centered at x p . As a result, I CDT EDU is formulated as I P EDU is the aggregate interference from PBSs. Because the classification of UEs is based on the locations of transmitters, we have that all possible interfering PBSs fall inside the region . Therefore, the interference I P EDU is written as , the interference from cluster-edge FBSs, is written as I EF EDU = ∑ y n ∈Φ EFT ∩b(0,R 3 ) P F h y n β y n −α . Lemma 2 is achieved, which gives the LT of interference Lemma 2: The LT of the aggregate interference I EDU experienced at the typical cluster-edge D2D-Rx is characterized by where (s), and L I P EDU (s) are given by Proof: See Appendix B.

Coverage probability of D2D network
With Lemma 1 and Lemma 2, we can calculate the coverage probability of D2D network transmission. As illustrated in Fig.1 and Fig.2, under D2D communication mode, when the typical cluster-center (or cluster-edge) D2D-Rx locating at x

C(E) D
is considered, the received SIR by the typical cluster-center (or cluster-edge) D2D-Rx is expressed as is the corresponding access distance between the typical D2D-Rx and its serving D2D-Tx and is . With the fact that the smallscale fading gain between D2D-Tx and D2D-Rx satisfies h (1), we have the coverage probability of the typical cluster-center (or cluster-edge) D2D-Rx given by Therefore, we achieve Theorem 1 that gives the achievable coverage probability of a typical D2D-Rx.
Theorem 1 Under D2D mode, where the typical cluster-center (cluster-edge) D2D-Rx is associated with a D2D-Tx, the achievable SIR coverage probability is where L I CDU (s) and L I EDU (s) for the cluster-center and cluster-edge D2D-Rxs are given by (7) and (14). The total coverage probability of the typical UE is C Tot where Λ In D and Λ Out D are the probabilities that a typical D2D-Rx locates inside cluster and outside cluster, respectively.

Coverage performance of cellular network
This section performs the analysis for the system under the case where the typical UE works in cellular mode with probability 1 − p. In this case, a typical cellular UE can also be associated with MBS, PBS, or FBS, which greatly depends on the strength of the received signals from cellular BSs. As commonly used, we adopt the long-term average received signal power model for UE association. In addition, a typical UE could also fall inside the cluster center or cluster edge. Without loss of generality, we define x C(E) z as the association distance of the typical UE to its serving BS z, z ∈ {F, P, M}.Then, using the method given by [43] and [44], we have Lemma 3 and Lemma 4.
Lemma 3: Given that the typical cluster-center UE is associated with a BS in the tier z ∈ {F, M, P}, the PDFs of the association distances x C F , x C M , and x C P are respectively given by where we define variables B zk = B z /B k , B z is the identical bias factor for tier z, P zk = P z /P k , and , and Λ C P are the association probabilities of the typical cluster-center cellular UE with FBS, MBS, and PBS, and given by Lemma 4: Given that the typical cluster-edge UE is associated with a BS in the tier z ∈ {F, M}, the PDFs of the association distances x E F and x E M are respectively given by where Λ E F and Λ E M are the association probabilities of the typical cluster-edge cellular UE with FBS and MBS, and given by

Association with FBS
We first consider the case where the typical UE is associated with a FBS from Θ F . In Θ F , the points falling inside the coverage region of PBSs is modeled by TCP and remaining points From Fig.2, because the cluster-center FBSs share sub-band W 2 with MBSs, the interference at a typical cluster-center FUE from Φ CFU is formulated as I CFU = I CFT CFU + I MT CFU . I CFT CFU is the interference from the cluster-center FBSs at the typical cluster-center FUE and I MT CFU is the one from MBSs. Similar to I CDT CDU , we have that all possible interfering cluster centers (PBSs) would locate in b (0, R 2 + R 3 ). By assuming the cluster center PBS locations x P ∈ b (0, R 2 + R 3 ), it is achieved that the typical cluster-center FUE receives interferences from all FBSs locating inside y c ∈ B x P ∩ B x C F , R 3 \x C F except the serving FBS locating at x C F for this typical FUE. Therefore, the cluster-center interference I CFT CFU is formulated as It is easy to see the LT of I CDT CDU and the one of I CFT CFU are similar, which is given by (8). The interference I MT CFU at the typical FUE from MBSs is formulated as , g x m is the small-scale fading interfering channel power gain, which follows gamma distribution. The LT of the interference I MT CFU is provided in Appendix C. Having Lemma 5. Lemma 5: The LT of the aggregate interference I CFU = I CFT CFU + I MT CFU experienced at the typical cluster-center FUE is characterized by where L I CFT CFU (s) is given by x C F and the association distance x C F is characterized by the PDF (21).
Proof: See Appendix C.
From Fig.2, we see that when FBSs locate outside clusters, these cluster-edge FBSs share sub-band W 1 with all PBSs and all D2D-Txs. Therefore, the interference at the typical cluster-edge FUE is written as I EFU = I EFT EFU + I CDT EFU + I EDT EFU . I EFT EFU is the interference from FBSs locating outside clusters except the serving FBS locating at x E F , which is written as denotes the set of interfering cluster-edge FBSs that fall inside the annulus of radii x E F and R 3 except the serving cluster-edge FBS locating at x E F . I CDT EFU is the interference from cluster-center D2D-Txs. Similar to I CDT EDU formulated by (13), in this case, to formulate the interference I CDT EFU , the locations of all possible interfering cluster centers are required firstly, which would be located in the annulusB (R 2 , R 2 + R 4 ). Then, by taking x P ∈B (R 2 , R 2 + R 4 ) as the set of all D2D-Txs in the cluster centered at x P and combining the association criterion for cluster-edge FUEs, it is easy to formulate the interference Note that, here we consider the fact that under cellular communications, a typical cellular UE experiences the interference from all D2D-Txs. For the interference I EDT EFU from the cluster-edge D2D-Txs at the typical FUE, I EDT EFU is formulated as I EDT EFU = ∑ y n ∈Φ EFT ∩b(0,R 4 ) P D h y n y n −α .
Therefore, we have Lemma 6.
F is the association distance from the typical cluster-edge FUE to its serving FBS and is characterized by the PDF (27).
Proof: See Appendix D.

Association with PBSs
Different from the UEs associated with FBSs and D2D-Txs, under the considered network model, the typical PUE associated with the PBS locating at x C P , is always located inside the coverage of cluster, i.e., cluster-center UE. From  Fig.2, we see that over the sub-band W 1 , the interference at the typical PUE is formulated as In (39), I P CPU is the interference at the typical PUE from other PBSs. With the association distance x C P between the typical PUE and its serving PBS, it is easy to achieve that all interfering PBSs fall inside the annulusB x C P , R 2 . The interference I P CPU is given by I P CPU = ∑ x P ∈Φ P ∩B( x C P ,R 2) \x C P P P g x p β x P −α , where the small scale fading interfering channel power gain g x P follows gamma distribution [63]-[65], i.e., g x P ∼ Γ (S P , 1) with the PDF f g x P (x) = 1 Γ(S P ) x S P −1 e −x . I EFT CPU is the interference at the typical PUE from cluster-edge FBSs, which is formulated as I EFT CPU = ∑ y n ∈Φ EFT ∩B(D F P (x C P ),R3) P F h y n y n −α . Note that, with the used association condition B P we have the distance from the nearest interfering FBS is D F P x C P = B FP P FP G FP 1 / α x C P . I CDT CPU is the interference from cluster-center D2D-Txs. Similar to I CDT CDU , it is easy to have that the interference I CDT CPU is written as I EDT CPU is the interference from cluster-edge D2D-Txs, which is given by I EDT CPU = ∑ y n ∈Φ EDT ∩b(0,R 4 ) β P D h y n y n −α . Then, we have Lemma 7.
where D F P x C P = B FP P FP G FP Proof: See Appendix E.

Association with MBSs
Since MBSs and cluster-center FBSs share sub-band W 2 , we have the interference at the typical cluster-center MUE is written as I CMU = I M CMU + I CFT CMU , where I M CUE is the interference from MBSs except the serving MBS to this typical cluster-center MUE. With the access distance x C M between the typical UE and its serving MBS at x C M , we can formulate the interference I M CMU as In I CMU , I CFT CMU is the interference from cluster-center FBSs at the typical intra-cluster MUE and is formulated as Hence, the LT of the interference I CFT CMU is calculated as  (46) and (47) Hence, the LT of the interference I CFT EMU is calculated as

Coverage probabilities of cellular UEs
Considering the fact that a cellular UE could be associated with one of the three types of transmitters, MBS, PBS, and FBS, and also could fall inside the cluster center or cluster edge, we use the following general form to model the received SIR by a typical cellular UE.
where P B,r x V B , B ∈ {F, P, M}, V ∈ {C, E}, denotes the received instantaneous signal power by the typical cellular UE from its serving BS locating at x V B with the association distance x V B , I V BU is the total experienced interference at the typical UE. Similar to (19), this yields that the general form of the SIR coverage probability of the typical cellular UE is formulated as With the general form of coverage probability, we first focus on the FBS downlink transmission, i.e., B = F,V = C(E). From Fig.1 and Fig.2, when the typical cluster-center (or cluster-edge) FUE is considered, the received instantaneous signal power by the typical FUE is given by P F,r x The small-scale fading channel gain follows h X C(E) F ∼ exp (1), the interference I C(E)FU are given by Lemma 5 and Lemma 6, respectively.
Using the similar method as (20) for the D2D network, it is easy to have the conditional SIR coverage probability C However, because of the fact that when the typical UE is associated with PBS, the typical PUE is always located inside the coverage of cluster, in (49) we have B = P and V = C. In this case, the received SIR by the typical PUE is denoted by SIR C P , the received instantaneous signal power is given by P P,r x C P = (P P /S P ) g X C P β x C P −α . With the fact that PBSs are equipped with N P antennas and simultaneously serve S P UEs, we have the small-scale fading channel gain g X C P following gamma distribution with parameter N P − S P + 1, i.e., g X C P ∼ Γ (N P − S P + 1, 1). The association distance is x C P and the experienced interference I PUE is given by Lemma 7. As a result, the DL received SIR coverage probability C C P x C P by the typical PUE can be achieved, of which proof is presented in Appendix F. Theorem 2 Under cellular model, a typical UE could be associated with one of the three types of transmitters, MBS, PBS, and FBS, and also could fall inside the cluster center or cluster edge. When the typical cluster-center (clusteredge) UE is associated with FBS, the SIR coverage probabilities C C(E) F is given by where L I C(E)FU (s) are given by Lemma 5 and Lemma 6, respectively, the PDFs of the association distances x C(E) F are given by (21) and (27), respectively. When the typical cluster-center UE is associated with PBS, the SIR coverage probabilities C C P is given by where k n 1 , n 2 , n 3 , n 4 = k! n 1 !n 2 !n 3 !n 4 ! is the well-known multinomial coefficient, ∑ n 1 +n 2 +n 3 +n 4 =k (...) denotes overall ktuple of nonnegative integers (n 1 , n 2 , n 3 , n 4 ) that satisfy the constant condition n 1 + n 2 + n 3 + n 4 = k. The LTs of the interferences I P CPU , I EFT CPU , I CDT CPU , and I EDT CPU are given by Lemma 7. The PDFs of the association distance x C P is given by (23). When the typical UE is associated with the MBS located at x   (22) and (28), respectively.
Proof: See Appendix F.

Simulations and numerical results
Based on the above derivations and analyses, the simulated and numerical results are presented in terms of the coverage performance of D2D and cellular networks under various network parameters. Unless otherwise specified, during the analysis, a 3-tier D2D underlaying heterogeneous cellular network is considered. It is assumed that a typical UE works in D2D mode with probability p = 0.5 and a D2D user is D2D-Tx with probability q = 0.5. Over the whole network, the path loss exponents of all links are assumed to identical and taken as α = 2. The locations of MBSs and PBSs are modeled as independent PPPs with densities λ M = (1 ∼ 10)MBSs/(π × Km 2 ) and λ P = 10λ M , respectively. The points in PCPs Θ F and Θ U are scattered around the common parent points with the variances σ 2 F = 10 ∼ 1000 and σ 2 D = 10 ∼ 1000. The average numbers of active FBSs and UEs arec D = 40 andc F = 10.The total available bandwidth is B = 20MHz. For a quick view, Table 2 lists the detailed description of parameters and corresponding values.

Average numbers of UEs and FBSs in a clusterc
Probability that a typical UE is D2D destination p = 0.5 Probability that that a D2D user is a D2D-Tx q = 0.5

Coverage performance of D2D communications
By taking the above parameters, Fig.3 first presents the comparison analysis of the coverage probability of D2D network. Fig.3 (a) provides the coverage probabilities C C D and C E D of cluster-center and cluster-edge UEs versus the coverage radius R 2 of PBSs. It is easy to see that when the typical D2D-Rx locates in the coverage region of PBS (cluster-center D2D user), the SIR coverage probability C C D of this cluster-center D2D-Rx decreases with R 2 under the case where R 2 is small. Then, with the continuous increase of R 2 , the coverage probability C C D increases slowly. Unlike the coverage probability C C D of cluster-center D2D user, we see that when the typical D2D-Rx does not locate inside (a) P P =35dBm (b) P P =30dBm Fig. 3: Coverage probability of D2D network the coverage region of PBS (cluster-edge UE), the SIR coverage probability C E D of this D2D-Tx transmission always decreases over the entire range of R 2 so that we have C C D > C E D when the coverage radius R 2 is large enough. Note that, we see that it is possible to achieve C C D < C E D when R 2 is small, e.g., R 2 = 400m. However, when the coverage radius R 2 is small enough, the result C C D > C E D is always achieved. For these observations, we have the following explanations. On one hand, when the coverage radius R 2 of PBS is small, the increase of R 2 results in the increase of the number of PBSs and FBSs in the coverage range of D2D-Txs. The increased PBSs and FBSs impose the increased interference on the typical D2D-Rx so that the coverage performance C C D decreases. Besides PBSs and FBSs, in this case the number of D2D-Txs also increases. On the other hand, when the coverage radius R 2 is large, the number of randomly distributed PBSs is saturated and the FBSs and D2D-Txs impose great and dominant impact on the performance C C D . In this case, the increase of R 2 indicates the increase of the average distance from FBSs and D2D-Txs to the typical D2D-Rx because both FBSs and UEs are modeled as TCPs and the numbers of FBSs and D2D-Txs in a typical cluster are constant. The increase of average distances yields the increase of path loss and the decrease of interference powers so that the coverage probability C C D increases. For the coverage probability C E D of the cluster-edge D2D-Rx, we have that when the typical D2D-Rx falls outside the coverage of a typical cluster, the interference from the cluster-center transmitters always increases with R 2 so that the coverage probability C E D decreases. This observation is further displayed in Fig.3 (b). Moreover, comparing Fig.3 (a) and (b), we see that when the transmission power P P is small, the SIR coverage probability C C D of the cluster-center D2D-Rx is greater than the one C E D of the cluster-edge D2D-Rx over the entire region of R 2 due to the decreased interference from PBSs. Finally, it is also achieved that the transmission power P P of PBSs imposes greater impact on C C D than C E D due to the main interference considered. As a result, we have the insight that, to guard the cluster-edge D2D-Rxs, the small R 2 is desired. On the contrary, to guard the cluster-center D2D-Rxs, we should increase R 2 .
Based on the observations in Fig.3, in Fig.4 we exploit the joint impact of both the coverage radii R 2 and R 4 of PBSs and D2D-Txs. From Fig.4, it is not difficult to see that as the coverage radius R 4 increases, the coverage probability C C D of the typical cluster-center D2D-Rx increases quickly, while the one C E D of the typical cluster-edge D2D-Rx decreases. The reason is that the increase of the coverage radius R 4 of the D2D-Txs indicates more D2D-Txs far from the cluster-center D2D-Rx so that the experienced interference at the typical D2D-Rx decreases. Contrarily, when the D2D-Rx falls outside the coverage range of PBSs, we have that the increase of R 4 indicates that the typical D2D-Rx suffers from more strong interference. However, from Fig.4 (c) we have the total SIR coverage probability C Tot D of a typical D2D-Rx increases with R 4 when R 2 is large, but decreases when R 2 is small. As a result, we have the total coverage performance of D2D network is dominated by cluster-center transmission, especially for large R 2 . When the gain of the total coverage performance is considered, we see that the achievable gain of the cluster-center D2D-Rx increases with the coverage radius R 2 and R 4 when R 2 is large. However, for the coverage probability C E D of the typical cluster-edge D2D-Rx, the effect of R 2 on the gain of coverage performance is negligible. This observation indicates that it is required to carefully select R 4 when it is small. We have that, to guard the cluster-center D2D-Rx, the large R 4 is preferable, to guard the cluster-edge D2D-Rx, the small R 4 should be chosen.
To achieve comprehensive insights of the effect of coverage radii R 2 and R 4 , a 3D plot is provided in Fig.5 with respect to coverage probability versus both R 2 and R 4 . As a result, a clearer observation is achieved in terms of the changing direction of the coverage probabilities C C D and C E D of the cluster-center and cluster-edge D2D-Rx, which can not be achieved in both Fig.3 and Fig.4. We see that when R 4 is small relatively, the coverage probability C C D decreases with R 2 . However, when R 4 is large relatively, the coverage probability C C D decreases with R 2 only when R 2 is small, and then C C D increases with the continuously increased R 2 . At the same time, Fig.5 (b) shows that over the entire ranges of both R 2 and R 4 , the coverage probability C E D of the typical cluster-edge D2D-Rx always decreases with both R 2 and R 4 . The observations from Fig.5 can be interpreted with the similar reason as in Fig.3 and Fig.4. Fig.6 gives the coverage probabilities versus the transmit power P D of D2D-Txs. As expected, we see that when the transmit power P D is small, both the coverage probabilities C C D and C E D decreases with P D due to the increased interference from D2D-Txs. However, we also see that when the transmission power P D is large, both the coverage probabilities C C D and C E D increase. The reason is that in this case, the signal from D2D-Txs is greatly stronger than the ones from MBSs, PBSs, and FBSs. Moreover, there exists a value of P * D , at which both the curves of C C D and C E D cross. When the transmission power P D is smaller than P * D , the coverage probability C E D is greater than the one C C D ; otherwise C E D is smaller than C C D . For these observations, we have the following explanations. Because of the main interference considered, the typical D2D-Rx that is associated with cluster-center D2D-Tx suffers from more severe interference than the one that is associated with cluster-edge D2D-Tx. However, when the transmission power P D is large, the system reduces to interference-free one due to the relatively small interference from FBSs and PBSs. In this case, the path loss of cluster-center D2D links is smaller than the one of cluster-edge links. Hence, we have the coverage probability C C D is greater than one C E D . For these observations, Fig.6 (b) and (c) further validate. Comparing Fig.6 (a)∼(c), we have that both the coverage probabilities C C D and C E D increase as P P decreases due to the decreased interference. Note that, in this numerical analysis the range of P D is 5 ∼ 24dBm. In practice, the transmission power P D is usually small. With above investigations, Fig.7 gives a comprehensive comparison by providing the 3D plot of (a) P P = 35dBm (b) P P = 30dBm (c) P P = 25dBm Fig. 6: Effect of D2D-Txs' transmission power P D coverage probabilities versus both P P and P D . In this figure, the similar results as in Fig.6 can be achieved. Hence, the detailed interpretation for this figure is ignored.

Coverage performance of cellular communications
In this subsection, we present the coverage analysis of cellular network. Fig.8 investigates the coverage probabilities of FUEs. Specially, Fig.8 (a) exploits the impact of the key parameters of femto cell, i.e., σ 2 F andc F . We see that the coverage probability of cluster edge FUE increases with the variance σ 2 F of FBSs' distribution. The reason is that the increase of σ 2 F indicates that more FBSs are far away from its cluster center so that more FBSs locate inside (a) Cluster-center D2D-Rx (b) Cluster-edge D2D-Rx Fig. 7: Comparison of coverage probabilities versus both P P and P D the cluster-edge region and the typical cluster-edge FUE experiences the decreased interference. We also see that the achieved performance gain increases with the average numberc F . Fig.8 (b) simultaneously investigates the coverage probabilities of both cluster-center and cluster-edge FUEs. The relationship between the FUE's coverage probability and the density λ M of MBSs is exploited. Because the MBSs and cluster-center FBSs share the sub-band W 2 , the coverage probability of the typical cluster-center FUE decreases with the density λ M due to huge cross-tier and cochannel interference from MBSs. Straightforwardly, we have that the density λ M of MBSs has less impact on the cluster-edge FUE. However, in practical implementation, the density of MBSs is very small, the interference from which is limited. Hence, the impact of the interference from MBSs on the coverage probabilities of FUEs is small. This validates the proposed SSA scheme. At the same time, Fig.8 (b) indicates that the achievable performance gain is related to the coverage radius of FBSs. We see that with the used system setting, the coverage probabilities' gap decreases with the density λ M due to the cluster-center FUEs sharing the sub-band W 2 with MBSs. Fig.9 investigates the coverage probabilities of PUE and MUE. From Fig.9 (a) that presents the coverage probability of the typical cluster-center PUE, we see that the transmit power P P imposes very limited impact on the coverage probability C C P of the typical PUE. As achieved previously, the coverage radius R 2 greatly affects the coverage probability C C P . The decrease of R 2 causes the increase of the coverage probability C C P because of the fact that the decrease of the coverage radius R 2 indicates the decrease of the interference PBSs and the decrease of interferences from other transmitters. At the same time, the interfering D2D users and cluster-edge FUEs also decrease. Considering the fact that the typical MUE can fall either inside cluster center or cluster edge, Fig.9 (b) simultaneously investigates the coverage probabilities of both the cluster-center and cluster-edge MUEs. We first observe the coverage probability of cluster-edge MUE. We have that due to the severe intra-tier interference from MBSs, it is a direct result that the coverage probability C E M decreases with the transmit power P M of MBSs. However, after observing the coverage C C M , the completely different results can be obtained. We have that the coverage probability C C M increases with the decrease of the radius R 2 . Moreover, the impact of the coverage R 2 on C C M is distinctive due to the decrease for the interference from other transmitters.

Conclusions
With the proliferation of high-speed multi-media applications and high-dense internet-of-things devices, the 5G/B5G networks have being faced the emergence of dense hotspot communications. A key consequence of such hotspots is the heterogeneity of BSs deployment and the user-BS coupling. With these motivations, by exploiting the poential of D2D communication, this paper considered a realistic scenario for clustered heterogeneous network, which consists of macro BSs, pico BSs, and femto BSs with underlaying D2D links. Under this clustered communications, the UEs are scattered around the center of hotspots. This works modeled the locations of MBSs and PBSs as PPPs and modeled the ones of UEs as PCP. To further efficiently offloading and balancing load, the low-power deployed FBSs are also distributed in hotspot and modeled as a PCP. The two PCPs of UEs and FBSs have the same parent process consisting of the set of PBSs. To restrict the severe interference from inter-tier and intra-tier transmitters, a clustered-UE classification was developed based on the coverage radius of PBSs. As a sequence, the UEs and FBSs outsider the coverage of PBSs were modeled two independent Poisson hole processes. With the clustered-based device classification, the paper developed a modified fractional frequency reuse scheme. The spectrum allocation guards both PBSs and FBSs by allocating different sub-bands. A UE can work either in D2D mode or cellular mode. Then, by using the methods from PPP, PCP, and PHP, this paper presented a tractable approach for modeling and analyzing the performance of cellular and D2D networks and gave the statistical descriptions of the experienced interference at a typical receiver by using the approximated PHP theory. Secondly, with an effective performance metric, this paper gave the coverage probabilities of both the D2D receiver and cellular destinations.
Appendix A: Proof of Lemma 1 With (3), the LT of the interference I CDT CDU is written as where (a) is obtained from the independent assumption of all interfering links, (b) is from the Rayleigh fading assumption h y c ∼ exp (1), (c) is obtained by applying probability generating function (PGFL) for a single cluster in TCP, i.e., [45]. Note that the average number of D2D-Txs per cluster is pqc D , where p is the probability that a UE has the required content by D2D-Rx and q is the probability that a D2D user is a D2D transmitter. Then, by applying the distribution of offspring points (UEs) and Cartesian to polar coordinates conversion, we have the L I CDT CDU (s) is further calculated by With the definition of I EDT CDU by (4) and the assumption of independence of all transmitters, the LT of the interference I EDT CDU from cluster-edge D2D-Txs is calculated as follows.
where (a) follows from the fact that only the dominant interference is considered, (b) is obtained by applying probability generating function (PGFL) of Φ U ; E R 2 = y∈Φ P b (y, R 2 ) is the coverage region of holes of radius R 2 and E C R 2 is the complement of E R 2 , (c) is from Cartesian to Polar coordinates conversion. In (56), A I EDT CDU is further calculated as where (d) is obtained by applying (3.194.1) in [46], 2 F 1 ( . , . ; . ; .) is Gauss hypergeometric function defined by (9.14.2) in [46]. The term B I EDT CDU is written as where (e) follows by applying the cosine-law: r 2 + x P 2 − 2r x P cos θ = R 2 2 , the interference at the typical UE from PHP with single hole is considered, (f) obtained by applying PGFL of PPP Φ P . Observing (58), we have that it is challenging to gain a tractable and simple form. To obtain a tractable but yet accurate form, we consider the following approximation. Note that, due to path loss, the effect of holes that are closest to the typical point will be much more significant compared to the holes that are further away. Therefore, to derive an easy-to-use form of (58), we only consider one hole, the one that is closest to the typical point, and ignore other holes. Denoting the location of the closest hole by x P 1 , the interference field in this case is Φ P ∩ b (y, R 2 ). Since the typical UE (D2D-Rx) always lies outside the hole, the closest point of Φ P is at least a distance R 2 from it. Using this fact and the law of PPP, it is achieved that the PDF of the distance v 1 = x P 1 can be shown to be With this consideration and the assumption that only the main interference is considered in this work, from (e) in (58), we have that (58) can be formulated by using the following tractable form.
Hence, substituting (57) and (59) into (56), the LT L I EDT CDU (s) is given by In (5), because we assume that each PBS is equipped with N P antennas and simultaneously transmits S P data streams with equal power assignment, with linear zero-forcing beamforming, we have that the channel gain following g X P ∼ Γ (S P , 1), whose PDF is given by f g X P (x) = 1 Γ(S P ) x S P −1 e −x . Then, the LT of the interference I P CDU is calculated by where (a) is obtained by using the result (3.351.1) in [46]. Then, by further appling (3.194.2) in [46], we have The LT of I EF CDU can be easily achieved by replacing the variables P D , R 4 and pλ U in (60) with P F , R 3 and λ F . Therefore, we have the following derivation.
Appendix B: Proof of Lemma 2 With the similar line to (12), the LT of the interference I EDT EDU is calculated as follows.
By using the independence, it is obtained that As defined in (56), E R 2 = y∈Φ P b (y, R 2 ) is the coverage region of holes of radius R 2 and E C R 2 is the complement of E R 2 . In (65), using Gauss hypergeometric function 2 F 1 ( . , . ; . ; .) defined by (9.14.2) in [37], the first part is given by (57). By taking the second part as B I EDT EDU , we have where (a) is achieved by using the similar line as (58). With the main interference only considered, we exploit the effect of the closest hole located at x P 1 . Since the typical D2D-Rx always lies outside the holes with reduced radius R 2 − R 4 , the closest point of Φ P is at least a distance R 2 − R 4 from it. Using this fact and the theory from stochastic geometry [45], the PDF of the distance v 1 = x P can be shown to be Therefore, (66) is given by Substituting (68) into (65), L I EDT EDU (S) is achieved, for clarity, and is given by In (12), I CDT EDU is the interference from cluster-center D2D-Txs at the typical UE whose access point D2D-Tx falls cluster-edge area. As stated previously, because the classification of UEs is based on the locations of transmitters, the typical D2D-Rx can fall inside cluster. This yields that in I CDT EDU , all possible interfering cell centers would be located in B (R 2 − R 4 , R 2 + R 4 ). Then, by taking x P ∈B, we further define B x p denoting the set of interfering D2D-Txs in the cluster centered at x p . As a result, I CDT ED is formulated as Hence, the LT of the interference I CDT EDU is calculated as follows.
Comparing (70)  In (12), I P EDU is the aggregate interference from PBSs. Because the classification of UEs is based on the locations of transmitters, we have that all possible interfering PBSs fall inside the region x P ∈ Φ P ∩ B (R 2 − R 4 , R 2 ). Therefore, the LT of interference I P EDU is calculated by Appendix C: Proof of Lemma 5 The LT of the interference I MT CFU is firstly provided. Under cellular mode, when the typical cluster-center UE is associated with a FBS, according to the used UE association criterion, we have This yields that the closest distance of interfering MBSs is D M F (x F ) = B MF P MF G MF 1 / α x C F . Note that, the association distance is x C F . Therefore, the interference I MT CFU at the typical FUE from MBS is formulated as where g x M is the small-scale fading interfering channel power gain, which follows gamma distribution. Because we assume that each MBS is equipped with N M antennas and communicates simultaneously with S M UEs, it is achieved that the PDF of g x M is given by Therefore, the LT of the interference I MT CFU is calculated as follows.
where (a) follows from the independent assumption of points in Φ M and Θ F , (b) is obtained by using (3.351.3) in [46], (c) is obtained by applying probability generating function (PGFL) for a single cluster in TCP, i.e., [45]. By using the binomial theorem, (72) is further calculated by Appendix D: Proof of Lemma 6 The LT of I EFT EFU from the cluster-edged FBSs is calculated by where x E F is the association distance from the typical cluster-edge FUE to its serving FBS. In (74), by using the result (3.194.1) in [46], the first part is written as By ignoring the effect of holes overlap and using the fact E R 2 dy n 1+(sP F β ) −1 y n α = ∑ x P ∈Φ P B(x P ,R 2 )c F dy n 1+(sP F β ) −1 y n 2 [32,33], the second part of (74) is given by For a give hole process B(x P , R 2 ), we have x P − R 2 y n x P + R 2 [32,33]. Then, using the cosin-law r 2 + v 1 2 − 2v 1 r cos(θ (r)) = R 2 2 and r = y n , the above expression is further given by At the same time, the effect of holes that are close to the cluster-edge FUE will be much more significant compared to the holes that are farther away. Therefore, to achieve an easy-to-use derivation, we only consider one hole; the one that is closest to the cluster-edge FUE and ignore the other holes. Since the typical cluster-edge FUE always lies outside the holes with the related radius R 2 , the closest point of Φ P is at least a distance R 2 from it. Using these fact, by assuming the closest point's distance as v 1 = x P , x P ∈ Φ P , we have the PDF of v 1 given by Then, deconditioning the expression (77) with respect to the distribution With the similar line as (70), the LT of the interference I CDT EFU is written as With (74) and (78), the LT of the interference I EDT EFU is written as where λ D2D−T x = pqc D λ P . Specially, by taking the first part of (80) as A I EDT EFU , we have Considering the typical cluster-edge FUE always lies outside the holes, the closest point of Φ P is at least a distance R 2 from it. Hence, taking the second part of (80) asB I EDT EFU and it is given by Substituting (82) and (81)  The interference from PBSs is written as I P CPU = ∑ x P ∈Φ P ∩B( x C P ,R 2) P p S P g x p β x P −α . With the similar argument as (62), by using g x P ∼ Γ (S P , 1) with the PDF f g x P (x) = 1 Γ(S P ) x S P −1 e −x , it is easy to immediately write the LT of I P PUE as To derive a tractable form for L I P CPU (s), we rewrite the integral in the expression of (83) as Using (3.191.6) in [46], the first part in (84) is calculated as 1 α x C P 2α+αS P (s (P P /S P )) S P (S P + 2/α) × 2 F 1 S P , S P + 2 α ; S P + 2 α + 1; − x C P α s (P P /S P ) − R 2α+αS P 2 (s (P P /S P )) S P (S P + 2/α) × 2 F 1 S P , S P + 2 α ; S P + 2 α + 1; −(R 2 ) α s (P P /S P ) Similar to (74) and (78), the LT of the interference I EFT CPU = ∑ y n ∈Φ EF ∩B(D F P ( x C P ),R3) P F h y n y n −α is expressed as = E Φ P exp −λ PcF (R 2 ∩B(D F P ( x C P ),R3))∩E C R 2 1 − 1 1 + sP F h y n β y n −α dy n = E Φ P exp −λ PcF (R 2 ∩B(D F P ( x C P ),R3))∩E C R 2 1 1 + sP F h y n β y n α dy n = E Φ P exp −λ PcF B (D F P ( x C P ),R3) 1 1 + sP F h y n β y n α dy n + λ PcF E C R 2 1 1 + sP F h y n β y n α dy n 2πλ PcF rdr  (s) given by (8), the LT of I CDT CPU = ∑ x P ∈Φ P ∩b(0,R 2 +R 4 ) ∑ y c ∈B X P ∩b(0,R 4 ) P D h y C β x P + y C −α is written as In (91), the first part is given by R 4 0 rdr 1+(sβ P D ) −1 r α = 1 α R 2 4 2 F 1 1, 1 − 2 α ; 1 − 2 α + 1; −(sβ P D ) −1 R α 4 . With the consideration that the cluster-edge D2D-Tx is determined by its associated UE, the nearest point is at least a distance R 2 − R 4 from it. The second part is expressed as E Φ P exp pqc D λ P E R 2 dy n 1 + (sβ P D ) −1 y n 2 = E Φ P exp ∑ x P ∈Φ P E R 2 pqc D λ P dy n 1 + (sβ P D ) −1 y n Obviously, L I EDT CPU (s) = L I EDT CD (s).

Appendix F: Proof of Theorem 2
Here, we only present the proof of the coverage probability C C P x C P for the association with PBS. When the typical UE is associated with the PBS locating at x C P with association distance x C P , the small-scale fading channel gain g X C P follows gamma distribution with parameter N P − S P + 1, i.e., g X C P ∼ Γ (N P − S P + 1, 1) = Γ (G P , 1), the PDF of g X C P is f g X C P (x) = 1 Γ(G P ) x G P −1 e −x , the coverage probability C C P x C P of the PBS DL transmission is written as where k n 1 , n 2 , n 3 , n 4 = k! n 1 !n 2 !n 3 !n 4 ! is the well-known multinomial coefficient, ∑ n 1 +n 2 +n 3 +n 4 =k (...) denotes overall ktuples of nonnegative integers (n 1 , n 2 , n 3 , n 4 ) that satisfy the constant condition n 1 + n 2 + n 3 + n 4 = k. Observing (94), it is found that the expression (94) contains the common term E I n e −sI . To further calculate the coverage probability  Spectrum allocation of D2D underlaying three-tier heterogeneous network Figure 3 Coverage probability of D2D network Coverage probability of D2D networks Coverage probabilities of PUE and MUE