I Introduction
In recent years, unmanned aerial vehicles (UAVs), commonly known as drones, have become a hotspot for wireless communications due to its unique attributes such as lowcost, mobility, and flexible reconfiguration[1]. In the meantime, in the process of standardization for 5G/B5G networks, UAVs are gradually being considered as an critical candidate to support diverse applications, such as reconnaissance, remote sensing, or working as temporal base stations[2]. The popularity of UAVs motivates the researchers to explore the opportunities for integrating UAVs into the existing wireless networks.
In the existing wireless networks, it is difficult to achieve universal connectivity due to the severe path loss and excessive intercell interference. One efficient approach to improve the coverage in currently deployed terrestrial cellular networks is to equip the UAVs as aerial base stations (ABS), augmented with the functionalities of terrestrial base stations (TBSs)[3]. As compared to terrestrial cellular networks, one distinct feature of UAV communications is that the AirtoGround (A2G) links are more likely to experience lineofsight (LoS) propagation which offers lower attenuation[4]. To further exploit the spectrum efficiency of the UAVaided networks, nonorthogonal multiple access (NOMA) has attracted much attention for its capability of serving multiple user equipments (UEs) at different qualityofservice (QoS) requirements in the same resource block[5, 6]. The key idea of NOMA is to employ a superposition coding (SC) at the transmitter and successive interference cancellation (SIC) at the receiver[7], which provides a good tradeoff between the throughput of the system and the UEs fairness. Therefore, by adopting NOMA techniques, the achievable spectrum efficiency of the networks can be improved.
On the other hand, using highfrequency band and densification will be two key capacityincreasing techniques for cellular networks, such as millimeter wave (mmWave) communications[8] and small cells. Deploying terrestrial mmWave small cells will offer high capacity when a connection is available. Motivated by the benefits of UAVaided networks, NOMA techniques, and mmWave transmissions, in this work we consider a UAVaided heterogeneous network (HetNets) where a mmWave terrestrial network coexists with a sub6GHz NOMA enabled aerial network. Note that the use of mmWave and microWave resources simultaneously is a feature of 5G/B5G networks, and their distinctive carrier frequencies avoid the intertier interference[9].
Ia Related Work and Motivation
Modeling and analyzing cellular networks with the aid of stochastic geometry has been widely adopted due to its accuracy and tractability. In the studies of UAVaided networks, the authors of [10] derived the coverage probability for a finite ABS network by modeling the locations of ABSs as a uniform binomial point process (BPP). The authors [11] analyzed the downlink coverage performance of UAVaided cellular networks when the UEs are clustered around the projections of ABSs on the ground. A framework was propoesed in [12] to analyze the behaviors of a ABSs network under a realistic A2G channel model which incorporates the LoS and nonlineofsight (NLoS) links. This work was further extended in [13] where the network comprises both ABSs and TBSs. Besides, instead of considering the average probabilistic path loss in most of works, the authors of [12, 13] considered more realistic LoS and NLoS transmissions, respectively. Multitier UAVaided networks were presented in [3, 14, 15]. Specifically, the authors of [3] and [14] proposed the multitier drone architecture based on the standard terrestrial path loss model. Furthermore, the multitier UAVaided networks based on the transmitteroriented or receiveroriented rules under a realistic A2G channel model were studied in [15]. The authors of [16] advocated a pair of strategies in UAVaided NOMA networks, i.e., the UAVcentric strategy for offloading actions and the usercentric strategy for providing emergency communications. In [17], a multipleinput multipleoutput (MIMO)NOMA enabled UAV network was proposed, where the outage probability and ergodic rate were evaluated in the downlink scenario.
Regarding the literature of stochastic geometry based HetNets systems, a system containing sub6GHz macrocells and mmWave small cells was analyzed in [18], where the macrocells provide universal coverage and the small cells provide high data rate when the mmWave LoS link is available. The authors of [19] studied the decoupled association in a sub6GHz and mmWave deployment from the resource allocation perspective. Building upon the above research contributions and the analytical tools of stochastic geometry, we propose an architecture of UAVaided HetNets where mmWave terrestrial networks coexist with a sub6GHz aerial networks, which has not been well studied in the literature. In contrast to the previously reported designs of UAVaided HetNets [13, 15], our proposed architecture poses three additional challenges: i) The NOMA techniques causes additional interference from the connected ABS to the served UE; ii) The channel ordering needs to be determined under the unique characteristics of A2G channels; iii) The UE association policy needs to be carefully designed under the existence of mmWave TBSs and sub6GHz ABSs.
IB Contributions and Organization
The primary contributions of this paper can be summarized as follows:

By taking advantage of unique attributes of UAVs and high transmission rate of mmWave, we propose a new model of HetNets where the sub6GHz NOMA enabled ABSs overlay the mmWave TBSs. We model the GroundtoGround (G2G) and A2G links incorporating the impact of LoS and NLoS path loss attenuations. We consider the LoS and NLoS transmissions separately, where two independent nonhomogeneous poisson point processes (PPPs) are formulated.

We develop a flexible association policy to address the coexistence of NOMA enabled ABSs and mmWave TBSs. Under this policy, we first derive the analytical expressions for the distance distributions given that the typical UE is associated with a TBS, a NLoS ABS or a LoS ABS.

We derive exact analytical expressions for the UAVaided HetNets in terms of coverage probability and spectrum efficiency. Additionally, the closedform coverage probability expressions are obtained in mmWave tier.

We provide the basic power allocation guidelines for the NOMA enabled networks, in which the targeted signaltointerferenceplusnoise ratio (SINR) threshold of the typical UE and the fixed UE both determine the coverage probability of a typical UE. We also provide insights on the HetNets design by numerical results, which demonstrate that our proposed NOMA enabled HetNets is capable of achieving superior performance compared with the conventional OMA enabled HetNets.
The rest of this paper is organized as follows. In Section II, the HetNets model and the association strategy are introduced. In Section III, new analytical expressions for distance distributions and association probabilities are derived. Then the coverage probability and spectrum efficiency of the network are investigated in Section IV and Section V, respectively. Our numerical results are demonstrated in Section VI, which is followed by the conclusions in Section VII.
Ii Network model
Iia Network Description
In this work, we present a two tier downlink UAVaided hybrid HetNets system shown in Fig. 1. In tier 1, the TBSs provide wireless connectivity to the ground UEs, where the spatial distribution of TBSs is modeled as an HPPP with density . In tier 2, the ABSs are deployed to enhance the coverage or boost the capacity. We assume that all the ABSs hover at a height and their horizontal locations form an HPPP with density . It is worth mentioning that the analysis in this network model is also applicable for the ABSs with different altitudes[10]. Specially, in tier 1 the TBSs are equipped with multiple antennas and the mmWave band is utilized to provide fast data rate in shortrange small cells, while in tier 2 the ABSs adopt sub6GHz and NOMA techniques in order to improve the coverage and freedom to serve multiple UEs. All the UEs are assumed to be equipped with a single antenna. Without loss of generality, the analysis is conducted based on a typical UE positioned at the origin. All the symbol notations are list in Table I.
IiB Channel Characteristics and Directional Beamforming in mmWave tier
Notation  Description 

,  HPPP of TBSs with density , HPPP of ABSs with density 
,  PPP of LoS ABSs with density , PPP of NLoS ABSs with density 
,  Scurve parameters 
,  Probability of LoS/NLoS links under the horizontal distance 
, ,  mmWave LoS radius, ABS height, fixed UE distance 
, ,  Minimum distance between the typical UE and the NLoS ABSs, LoS ABSs and TBSs 
, ,  Path loss exponent between the typical UE and NLoS ABS, LoS ABS and TBS 
, ,  Additive loss exponent between the typical UE and NLoS ABS, LoS ABS and TBS 
, ,  Nakagami fading parameters between the typical UE and NLoS ABS, LoS ABS and TBS 
,  NOMA coefficients 
, ,  TBS transmit power, ABS tranmit power, noise power 
In the first tier of the networks, due to the deployment of mmWave, the transmitted signals suffer from attenuation due to the obstacles, and the blockage effect can not be neglected. Here we adopt a tractable equivalent LoS ball model [20] to characterize the blockage effect, which enables fast numerical computation and simplifies the analysis. We define a LoS radius , which represents the maximum distance between a UE and its potential mmWave TBS, and the LoS probability is one within and zero outside this radius. It has been shown that the LoS ball model can fit the real environment properly and provide enough analytical accuracy compared with other blockage models [20]. Regarding the NLoS paths, it has been pointed out in [21] that the impact of NLoS signals and NLoS interference can be ignored in mmWave networks. Hence, we will focus on the analysis where the typical UE is associated with a LoS TBS. Then the process for the case the TBSs located inside the LoS Ball can be expressed as . As a result, the path loss in the TBS tier can be expressed as
(1) 
where is the communication distance, is the unit step function, and are the additive loss and path loss exponents, respectively. We also characterize the small scale fading with the Nakagami
fading, where the channel gain follows the Gamma distribution with parameter
.In this work, multiple antennas are equipped at the TBSs to accomplish the directional beamforming, and we adopt a 3D sectorized model [22, 23]. The directivity gain is given by , where is the antenna 3dB beamwidth for the azimuth orientation in the horizontal direction and is the antenna 3dB beamwidth for the elevation angles in the TBS, with mainlobe gain and sidelobe gain . By adjusting the antenna direction toward the corresponding UE, the UE benefits from the high mainlobe gain . Moreover, since we evaluate the average directivity gain in our systems, the effect of misalignment is ignored in the rest of this paper. Considering the interfering transmission, the beamforming gain and its association probability can be expressed as follows
(2) 
where we assume that the azimuth angle
is uniformly distributed in the range of
, and the depression angle is uniformly distributed in the range for all interfering transmissions.IiC Channel Characteristics in the NOMA Enabled tier
In the second tier of the networks, the channel between the ABS and the UE is highly affected by the density and altitude of obstacles in the environment, then the A2G channels contain both LoS and NLoS links. In this paper, we adopt a measurement based on the probabilistic model for LoS/NLoS propagation [4], which is suitable for sub6GHz scenarios. The probability expressions of LoS/NLoS links are defined as and , respectively. The LoS link is shown as
(3) 
where and are referred to as the Scurve parameters which are related to the transmission environment, and denotes the horizontal distance between the ABS and the UE. Note that a higher ABS altitude results in a higher LoS probability due to fewer obstacles. Accordingly, the NLoS link probability is given by
Since each link between the ABS and the UE is either in a LoS or NLoS condition with probability and , the set of ABSs can be divided into two independent nonhomogeneous PPPs and with , which denote the LoS ABSs and the NLoS ABSs, respectively. The corresponding densities of and with respect to the horizontal distance from the typical UE are given by and , respectively. In the following, we use and to denote the horizontal locations of the LoS and NLoS ABSs, respectively.
We consider different channel parameters for LoS/NLoS links in the A2G channels. The additive loss and path loss exponents of the LoS link in the A2G channels are denoted by and , and accordingly we introduce and in the NLoS state. Therefore, the channel gains are denoted as and , which follow the Gamma distribution with parameter and , respectively.
IiD UE Association
In this UAVaided HetNets, a UE is allowed to access the mmWave tier or the NOMA enabled tier in order to provide the best coverage. The flexible UE association is based on the maximum longterm averaged received power at the UE of each tier. Intuitively, the typical UE will choose to connect to the BS which has the minimum distance to the UE for both of the tiers.
IiD1 Average Received Power in mmWave Tier
Denoting as the minimum horizontal distance between the typical UE and the TBSs. We assume that the transmit power of all the TBSs is . Thus, the average received power at the UE connected to the TBS is given by
(4) 
where is the directional beamforming gain.
IiD2 Average Received Power in the NOMA Enabled Tier
In the NOMA enabled tier, we adopt UE pairing to implement NOMA in order to reduce the complexity[24]. Compared with the UE association in the OMA scheme, NOMA allocates different power levels to multiple UEs by exploiting power sparsity. The locations of the UEs are also not predetermined due to the random spatial topology of the stochastic model. As such, we always assume that a near UE is chosen as the typical one first no matter it lies in a LoS/NLoS state[25], and we denote and as the minimum distance between the typical UE and the LoS/NLoS ABSs. Then the average received power at the UE connected to the LoS/NLoS ABS can be expressed as
(5) 
and
(6) 
respectively, where denotes the transmit power of all the ABSs, and denotes the power allocation factor for the near UE.
IiE SINR Analysis
Due to the fact that the TBS tier and ABS tier utilize distinctive carrier frequencies, the signals in these two tiers do not affect each other.
IiE1 mmWave Tier Transmission
The SINR at the typical UE when it is connected to a TBS at a distance can be expressed as
(7) 
where , and is the interference from the TBS tier. is the channel gain between the typical UE and the serving TBS, and refer to the channel gain and the distance between typical UE and the TBS (except for the serving BS), respectively. The value and probability of can be obtained through (2). denotes the noise power. Both and follow the distribution of .
IiE2 NOMA Enabled Tier Transmission
In the ABS tier, without loss of generality, we consider that each ABS is associated with one UE in the previous round of the UE association process. For simplicity, we follow the assumption in [25] where the distance between the associated UEs and the connected ABSs are the same, which are arbitrary values and denoted as . Since the path loss is more dominant compared with the small scale fading, we apply the SIC operation at the near UE side. However, it is not predetermined that the typical UE is the near UE or the far UE, we have the following near UE case and far UE case. We first assume that the typical UE is in a LoS state.
Near UE in a LoS state case
When the typical UE is the near UE in a LoS state, i.e., , the typical UE will first decode the information of the fixed UE to the same LoS ABS with the following SINR
(8) 
where denotes the power allocation factor for the far UE which satisfies and . denotes the interference from the LoS ABSs and denotes the interference from the NLoS ABSs. is the channel gain between the typical UE and the associated LoS ABS. and refer to the channel gain and the distance between the typical UE and LoS ABS (except for the serving ABS), respectively. and refer to the channel gain and the distance between the typical UE and NLoS ABS , respectively. and follow the distribution of . follows the distribution of .
If the information of the fixed UE can be decoded successfully, the typical UE will decode its own message with the following SINR
(9) 
For the fixed UE (far UE) served by the same ABS, the signal can be decoded by treating the message of the typical UE as interference, then the SINR for the fixed UE can be expressed as
(10) 
where refers to the channel gain between the fixed UE and the serving LoS ABS.
Far UE in a LoS state case
On the other hand, when the typical UE has a larger distance to the serving LoS ABS than the fixed UE, i.e., , the fixed UE will first decode the information of the typical UE with the following SINR
(11) 
Once the information of the typical UE can be decoded successfully, and by applying the SIC technique, the SINR to decode its own message at the fixed UE is given by
(12) 
For the typical UE that connects to the same LoS ABS, the SINR can be expressed as
(13) 
NLoS state case
When the typical UE is associated to a NLoS ABS, the typical UE can be a near UE or a far UE as well. The SINR analysis in NLoS state is similar to that of the LoS state case, and we omit the details here.
Iii Relevant distance distributions and Association Analysis
In this section, we focus on providing the distribution of the distances between the typical UE and the serving TBS, NLoS ABS and LoS ABS, respectively, in the HetNets system. Furthermore, we derive the expressions for the association probabilities. At last, the distance distributions are characterized given that the typical UE is associated with a TBS, a NLoS ABS or a LoS ABS.
Iiia Distance Distributions of the nearest BSs
Lemma 1.
Proof.
Using a similar method to Lemma 1 of [13], the above expressions can be obtained. ∎
IiiB Distance of the nearest interfering BSs
Due to the deployment of mmWave and the special channel characteristics of A2G channels, the distance of the nearest interfering BSs is not easy to observe. The following remarks show clear insights on the locations of the nearest interfering BSs, which will be useful to derive the main results of this paper.
IiiB1 The typical UE is associated with a NLoS ABS
When the typical UE is associated with a NLoS ABS, we have the following lemma and assumption, which help us to derive the minimum distance of the interfering LoS ABS.
Lemma 2.
The probability that the typical UE has at least one TBS in can be calculated by .
Assumption 1.
When there exists a TBS in , the typical UE never associates with a NLoS ABS.
If there exists a TBS in , and the typical UE is associated with a NLoS ABS. Then the average received power at a height connecting to a NLoS ABS should not be smaller than the average received power at a distance connecting to a TBS, resulting in the condition satisfied. With a huge path loss and additive loss in NLoS channels, this condition is not satisfied normally.
Remark 1.
Given that the typical UE is associated with a NLoS ABS at a distance , the minimum distance of the interfering LoS ABS is given by
(17) 
IiiB2 The typical UE is associated with a LoS ABS
When the typical UE is associated with a LoS ABS at a distance , depending on whether there exists a TBS in , we have two cases. We first denote , which is the distance between the typical UE and a LoS ABS when its average received power is the same as that the average received power from a NLoS ABS at a distance of , i.e., . We then denote as the maximum association distance between the typical UE and a LoS ABS when there exists a TBS in . As such, if there exists a TBS in , then only the LoS ABSs which lie in the range of can be associated to the typical UE.
Remark 2.
Given that the typical UE is associated with a LoS ABS at a distance , the minimum distances of the interfering TBS and NLoS ABS are given by
(18) 
If there exists a TBS in , and is satisfied, we have
(19) 
Otherwise when , we have
(20) 
If there does not exist a TBS in , can be expressed as
(21) 
Proof.
The proof is similar to the proof in Remark 1, therefore it is omitted here. ∎
IiiB3 The typical UE is associated with a TBS
When the typical UE is associated with a TBS, the condition should be satisfied. We first denote as the distance between a UE and a TBS when its average received power is the same as the average received power of a LoS ABS at a distance of , and similarly we denote .
Remark 3.
Given that the typical UE is associated with a TBS at a distance , the minimum distances of the interfering LoS ABS and NLoS ABS are shown below
Condition  





Proof.
The proof is similar to the proof in Remark 1 and Remark 2, therefore it is omitted here. ∎
IiiC Association Probability
We first study the probability that typical UE is associated with a NLoS ABS.
Lemma 3.
The probability that a typical UE connects to a NLoS ABS can be calculated as
(22) 
Proof.
Using a similar method to Lemma 2 of [13], and considering there does not exist a TBS in , the above expressions can be obtained. ∎
We then study the probability that the typical UE is associated with a LoS ABS.
Lemma 4.
The probability that the typical UE has at least one LoS ABS when can be calculated by .
Lemma 5.
The probability that a typical UE connects to a LoS ABS can be calculated as
(23) 
where and are the association probability related to the cases when there does not exist a TBS and there exists a TBS in , respectively. and are given by
(24) 
and
(25) 
respectively, where , if , or
(26) 
and
(27) 
Proof.
See Appendix A. ∎
Finally, we study the probability that the typical UE is associated with a TBS.
Proposition 1.
The probability that a typical UE connects to a TBS can be calculated as
(28) 
IiiD Conditional distance distribution of the serving BS
Now, denoting , and as the minimum distance between the typical UE and the serving BS given that it is associated with the NLoS ABS, LoS ABS and TBS, respectively. The following lemmas characterize the distributions of , and .
Lemma 6.
The probability density function (PDF) of is given by
(29) 
Proof.
Denoting as the event that the typical UE is associated with a NLoS ABS. Then we have
(30) 
where can be derived as
(31) 
Then taking the first order derivative, we can obtain the expressions in (29). ∎
Given that the typical UE is associated with a LoS ABS, depending on whether there exists a TBS in , and utilizing the results in Remark 2, we have the following lemma.
Lemma 7.
The PDF of is given by or , where is the PDF of when there does not exist a TBS in , and is the PDF of when there exists a TBS in . and are given by
(32) 
When the condition is satisfied,
(33) 
When the condition is satisfied,
(34) 
Proof.
The proof is similar to that in Lemma 6, therefore it is omitted here. ∎
Given that the typical UE is associated with a TBS, and we consider a typical condition when . Utilizing the results in Remark 3, we have the following lemma.
Lemma 8.
When the condition is satisfied, the PDF of is given by
(35) 
Iv Coverage Probability Analysis
The coverage probability is defined as that a typical UE can successfully transmit signals with a targeted SINR. To begin with, we derive the Laplace transform of the interference.
Iva Laplace Transform of Interference
Since the ABSs and the TBSs utilize different frequencies, the intertier interference can be avoided. We first consider the case that the interference signals are from the TBSs.
Lemma 9.
The Laplace transform of interference from TBSs to a typical UE is given by
(36) 
where , , and is the incomplete Beta function.
Proof.
See Appendix B. ∎
Then we consider the scenario where the interference signals are from ABSs.
Lemma 10.
The Laplace transform of interference from ABSs to a typical UE is given by
(37) 
where , , , and are given as follows.
Condition  

Proof.
The interference from ABSs contains both the interference from LoS ABSs and interference from NLoS ABSs , then the Laplace transform of equals to . Following the same steps in Lemma 9, we can obtain the above expressions. ∎
IvB Coverage Probability in mmWave Tier
Given that the typical UE is associated with a TBS, the coverage probability is defined as
(38) 
where is the targeted SINR of the typical UE.
Lemma 11.
The exact and approximated expression of can be expressed as
(39) 
where , and