Ejder Baştug

Mehdi Bennis
email: bennis@ee.oulu.fi

Marios Kountouris
email: marios.kountouris@huawei.com

Mérouane Debbah
email: merouane.debbah@centralesupelec.fr

Edge Caching for Coverage and Capacity-aided Heterogeneous Networks

Keywords: edge caching, clustered point processes, heterogeneous cellular networks, stochastic geometry

A two-tier heterogeneous cellular network (HCN) with intra-tier and inter-tier dependence is studied. The macro cell deployment follows a Poisson point process (PPP) and two different clustered point processes are used to model the cacheenabled small cells. Under this model, we derive approximate expressions in terms of finite integrals for the average delivery rate considering inter-tier and intra-tier dependence. On top of the fact that cache size drastically improves the performance of small cells in terms of average delivery rate, we show that rate splitting of limited-backhaul induces non-linear performance variations, and therefore has to be adjusted for rate fairness among users of different tiers.

I. INTRODUCTION

Caching in heterogeneous cellular networks (HCN) significantly improves the system performance, and is of cardinal importance especially in limited-backhaul settings [START_REF] Baştug | Cacheenabled small cell networks: Modeling and tradeoffs[END_REF], [START_REF] Shanmugam | Femtocaching: Wireless content delivery through distributed caching helpers[END_REF]. However, existing stochastic geometry-based analyses of caching in heterogeneous networks [START_REF] Yang | Analysis on cache-enabled wireless heterogeneous networks[END_REF], [START_REF] Chen | Cooperative caching and transmission design in cluster-centric small cell networks[END_REF] ignore the impact of limited-backhaul and consider multiple tiers of mutually independent point processes. To remedy to this situation, this paper analyzes the benefits of edge caching in such wireless deployment conditions.

Consider a heterogeneous network consisting of mobile users, clustered cache-enabled small base stations (SBSs), macro base stations (MBSs) and central routers (CRs) for the backhaul. For SBSs, we shall consider two different topologies: 1) coverage-aided and 2) capacity-aided deployments. Coverage-aided deployment follows a Poisson hole process (PHP) while capacity-aided deployment is modeled using a Matérn cluster process (MCP). This models intertier and intra-tier dependence, respectively. Due to the multitier dependence and the non-tractability of the considered point processes, an exact calculation of the interference and key performance metrics seems unfeasible. Hence, we restrict ourselves to approximations in the form of finite integrals and validate their behavior via numerical results.

This hierarchical model exploiting random spanning trees hides the algorithmic details of such a complex network, yet This research has been supported by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering), the project BESTCOM, the Academy of Finland CARMA project and TEKES grant (2364/31/2014).

allows for tractable mathematical expressions to characterize the overall performance of such an heterogeneous network. We show that while the average delivery rate of small cell users (SUs) can be improved by adding more storage per SBS, the backhaul rate splitting among tiers induces non-linear performance variations and requires adjustments for the rate fairness.

II. SYSTEM MODEL

The system model is composed of MBSs, cache-enabled clustered SBSs, CRs and mobile users. We focus on two different topologies as in [START_REF] Deng | Heterogeneous cellular network models with dependence[END_REF].

A. Coverage-Aided Topology

We model a heterogeneous cellular network which consists of MBSs and cache-enabled SBSs. The MBSs are modeled by an independent homogeneous planar Poisson point process (PPP) of intensity λ mc , denoted by Φ mc = {y i } i∈N where y i denotes the location of the i-th MBS. Additionally, the potential SBSs are located according to another independent homogeneous PPP of intensity λ sc , denoted by Φ sc = {x i } i∈N where x i represents the location of the i-th SBS. We suppose that each MBS has an exclusion region which is made of a disk with radius R c centered at the position of MBS. Assuming that the aim of SBSs is to fill the coverage holes of MBSs to provide a better service to users, these SBSs are deployed outside of the exclusion regions of MBSs. Therefore, the deployed SBSs form clusters according to a Poisson hole process (Cox process) as follows [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF].

Definition 1 (Clustering process of coverage-aided SBSs). Let Φ mc be a homogeneous PPP of intensity λ mc for MBSs and Φ sc be an independent and homogeneous PPP of intensity λ sc for potential SBSs, with λ sc > λ mc . For each y ∈ Φ mc , remove all the points in

Φ sc ∩ B(y, R c ) (1)
where B(y, R c ) is the ball of radius R c centered at y. Then, the remaining points of Φ sc form clusters, known as the Poisson hole process Φ sc and represent the deployed SBSs. Moreover, this process has the intensity of On the other hand, CRs are distributed in the plane according to an independent homogeneous PPP of intensity λ cr , denoted by Φ cr = {u i } i∈N . These routers are in charge of providing broadband Internet connection to MBSs and SBSs via backhaul links. Mobile user terminals are also positioned in the whole plane according to an independent homogeneous PPP of intensity λ ut , denoted by Φ ut = {z i } i∈N , with λ ut > λ mc and λ ut > λ sc . The snapshots of these point processes are given in Fig. 1.

λ sc = λ sc exp(-λ mc πR 2 c ). (2

B. Capacity-Aided Topology

Suppose a two-tier heterogeneous cellular network consisting of MBSs and cache-enabled SBSs. The MBSs are distributed according to an independent homogeneous PPP of intensity λ mc , denoted by Φ mc = {y i } i∈N where y i denotes the position of the i-th MBS. The SBSs are placed in hot-spots to sustain the demand of highly concentrated users, according to an independent Matérn cluster process Φ sc = {x i } i∈N whose parent PPP Φ sc has intensity λ sc . The process is given as follows [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF].

Definition 2 (Clustering process of capacity-aided SBSs). Let Φ sc be a parent process modelled by a homogeneous PPP of intensity λ sc . Then, the clustering process of SBSs is given by

Φ sc = x ∈Φ sc N x (3)
where N x is a Poisson number of independent and identically distributed (i.i.d.) points with mean c, distributed uniformly in the ball B(x , R c ). Then, the process Φ sc is called Matérn cluster process Φ sc and has intensity of

λ sc = λ sc c. (4)
For the users, we suppose that mobile users (both macro cell users (MUs) and SUs) are distributed on the two-dimensional Euclidean plane, however SUs are highly concentrated in hotspot regions served by SBSs. From this motivation, we assume that all mobile users are distributed according to a Cox process Φ ut = {z i } i∈N [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF]. In particular, the population centers of radius R c are drawn from the parent PPP Φ sc , where SUs in these clusters are distributed uniformly at random and are covered by SBSs forming a Matern cluster process Φ sc . By doing so, these mobile users are (on average) covered by their SBSs deployed in these hot-spot areas.

Recalling c as the average number of SBSs per cluster, the density of active SUs per cluster is then

λ ut-s = c πR 2 c .
The MUs distributed in the rest of the network follow a PPP with density λ ut-m and are served by their own MBSs. For convenience, we consider that each MBS serves only one MU on average and the same holds for each SBS and its user. With this consideration in mind, the densities of MUs and SUs are equal to that of the MBSs and SBSs, respectively. Under this setting, the MUs and SUs form a Cox process with density λ ut = λ mc + λ sc , clustered in hot-spots and uniformly distributed elsewhere. On the other hand, we consider that CRs are modeled by an independent homogeneous PPP of intensity λ cr , denoted by Φ cr = {u i } i∈N . A snapshot of the topology is depicted in Fig. 2.

C. Signal Model, Connectivity and Backhaul

The downlink transmissions of MBSs and SBSs occur at the same frequency with reuse factor 1. The transmit power is P mc for each MBS and P sc for each SBS. All nodes (MBSs, SBSs and users) have single antenna. Having an MBS positioned at y and receiver at z (or simply call as transmitter y and user z), the channel coefficient is denoted by h y,z ∈ C. In case of SBS as a transmitter, the channel coefficient between transmitter x and user z is given by g x,z ∈ C. All the channel power coefficients are i.i.d. exponential random variables (Rayleigh fading) with

E[|h y,z | 2 ] = 1 and E[|g x,z | 2 ] = 1.
Supposing that the downlink rate of the typical user is a function of received signal-to-interference ratio (SIR), the target rate of signaling MBS and SBS are given by τ mc and τ sc respectively.

Each user is either associated to the nearest MBS or nearest SBS. The backhaul connection of each base station is provided from its nearest CR. Supposing that each CR has a sufficiently high-capacity broadband Internet connection, MBSs and SBSs are connected to their nearest CRs via errorfree wired backhaul links. In particular, each CR has a total capacity of C cr , which is an exponentially distributed random variable with mean µ.

D. Caching

We assume that the global content popularity distribution of users follow a power law defined as [START_REF] Newman | Power laws, Pareto distributions and Zipf's law[END_REF]

f pop (f, η) = (η -1) f -η , f ≥ 1, 0, f < 1, (5)
where f indicates a point in the support of the corresponding content and the parameter η models the steepness of the distribution. In this work, we assume that f pop is perfectly known at the SBSs. In f pop (f, η), the contents in the interval [1, F ) are cacheable and called as catalogue, whereas the remaining part [F, ∞] is called as non-cacheable contents (i.e., voice traffic, online gaming and sensor information). An interval [f, f + ∆f ) in the support of f pop (f, η) is considered as the probability of f -th content. Indeed, we assume that each content in the catalogue has a fixed length and called as chunk. Each chunk can belong to a part of cacheable video file, audio or picture and so on. In fact, storing/distributing fixed-length chunks is one of the key principle in content centric networks [START_REF] Jacobson | Networking named content[END_REF]. Therefore, even though we use the term "content" in the paper, the function f pop (f, η) is actually a chunk popularity distribution. Each SBS has a storage capacity of F sc contents/chunks, with 1 ≤ F sc ≤ F .

E. Hierarchical Model

The coverage and capacity-aided deployments of SBSs together with MBSs, users and CRs can be modeled as random stationary graphs (hierarchical spanning trees) [START_REF] Suryaprakash | Modeling backhaul deployment costs in heterogeneous radio access networks using spatial point processes[END_REF]. In particular, a random hierarchical tree whose root node is a CR located at u is given by Ψ = {(u, v u )}, where v u is a mark vector containing all random variables associated with the CR. In particular, the mark vector v u contains information of MBSs and SBSs (with their users) which are associated to the CR at u, such as:

• N mc : The number of MBSs connected to the CR located at u. The vector r mc ∈ C 1×Nmc is the relative position vector of those MBSs connected to the CR at u. Therein, each element r u,y represents the distance from CR at u to MBS at y. These positions are conditioned on N mc . • N mu : The number of users connected to an MBS located at y ∈ {r mc }. The vector r mu ∈ C 1×Nmu is the relative positions of those MUs which are conditioned on N mu . Each element r y,z distance from MBS at y to its user z. • N sc : The number of SBSs connected to the CR located at u. The vector r sc ∈ C 1×Nsc is the relative positions of those SBSs connected to the CR at u. Here, each element r u,x represents the distance from CR at u to SBS at x. These positions are conditioned on N sc . • N su : The number of users connected to an SBS located at x ∈ {r sc }. The vector r su ∈ C 1×Nsu is the relative positions of users connected to the SBS at x ∈ {r sc }.

They are conditioned on N su . Each element r x,z represents distance from SBS at x to its user at z.

III. PERFORMANCE METRICS

For the performance metrics of coverage and capacity-aided deployments, we first start by defining SIR. Definition 3. The SIR of a typical user connected to an MBS (namely typical MU) located at random position y is defined as

SIR mu P mc h y (y) I mm + I sm (6)
where (y) = y -α is the standard power-law path loss function (unless otherwise stated) with exponent α, I mm = yi∈Φmc\{y} P mc h yi (y i ) is the cumulative interference from other MBSs except the serving cell at y, and

I sm = xi∈Φsc P sc g xi (x i ) is the total interference from clustered
SBSs. Similarly, the SIR of a typical user connected to an SBS (namely typical SU) located at random position x is given by

SIR su P mc g x (x) I ss + I ms (7)
where The amount of backhaul rate allocated to typical users are defined by the following policy. Definition 4 (Backhaul Rate Splitting Policy). Following the hierarchical model, suppose that a typical user located at z ∈ {r mu } (typical MU) is connected to the MBS at y ∈ {r mc }, and this MBS is connected to the nearest CR at u. Then, the rate of backhaul link to the MBS at y is given as

I
R mu γC cr E N mc N mu , (8)
where γ ∈ [0, 1] is a fraction of capacity allocated to the MBSs. In case of SBS, in a similar vein, a typical user located at z ∈ {r su } is connected to the SBS at x ∈ {r sc } whose CR is at u. Then, the rate of backhaul link to the SBS is given as

R su (1 -γ)C cr E N sc N su . ( 9
)
The expectation above is taken over the topology generated by the point processes of users and base stations. We now define our main performance metric as follows.

Definition 5 (Delivery Rate). The delivery rate of a typical user connected to an MBS is defined as

R mu τ mc , if log(1 + SIR mu ) > τ mc and R mu > τ mc , 0, otherwise. (10)
Similarly, the delivery rate of a typical user connected to an SBS is defined as

R su      τ sc , if log(1 + SIR su ) > τ sc and R su > τ sc , τ sc , if log(1 + SIR su ) > τ sc and f z ∈ ∆ x , 0, otherwise, (11)
where f z represents the content requested by the typical SU and ∆ x is the cache of the SBS.

We are now ready to give the expressions for the average delivery rate of typical MUs and SUs.

IV. MAIN RESULTS

Theorem 1 (Average Delivery Rate of Typical MU). The average delivery rate of a typical user connected to the nearest MBS cell in coverage-aided deployment is approximated as

R(cov) mu ≈ τ mc B (cov) 1 B (cov) 2 (12)
where B , and F (x, y; z; w) is the hypergeometric function. For capacity-aided deployment, we have

R(cap) mu ≈ τ mc B (cap) 1 B (cap) 2 (13) where B
LI mm (s) = exp -sπλmcPmc(2/α) 1 -2/α r 2-α mc × F 1, 1 -2/α; 2 -2/α; -sPmcr -α mc (15) LI sm (s) = exp -λ sc (sPsc) 2/α π 2 (2/α) sin(π 2 α ) - πR 2 c Amc(s, Rc) (16)
Amc(s, Rc) = 1 πR 2 c 2π 0 rmccosϕ+ √ R 2 c -r 2 mc sin 2 ϕ 0 rdrdϕ 1 + s -1 P -1 sc r α (17) B (cov) 2 = 1 -exp - τmcλcr λmc + λ sc exp(-λmcπR 2 c ) µγλ 2 mc λut (18)
B (cap) 1 = Rc 0 e - (e τmc -1) Pmc r -α mc LI mm e τmc -1 Pmcr -α mc × LI sm e τmc -1 Pmcr -α mc k ν rmc ν k-1 e -(rmc/ν) k drmc (19) LI mm (s) = exp -sπλmcPmc(2/α) 1 -2/α r 2-α mc × F 1, 1 -2/α; 2 -2/α; -sPmcr -α mc (20) LI sm (s) = exp -λ sc R 2 1 -exp(-cν(s, y)) dy (21) ν(s, y) = R 2 f (x) 1 + (sPsc (x -y)) -1 dx (22) f (x) = 1 πR 2 c , if x < Rc, 0, otherwise. (23)
B (cap) 2 = 1 -exp - τmcλcr µγλut-m (24)
Theorem 2 (Average Delivery Rate of Typical SU). The average delivery rate of the typical user connected to the nearest SBS in coverage-aided deployment is approximated as

R(cov) su ≈ τ sc C (cov) 1 C (cov) 2 + τ sc C (cov) 1 C (cov) 3 (25) -τ sc C (cov) 1 C (cov) 2 C (cov) 3 (26)
where

C (cov) 1 , C (cov) 2 and C
(cov) 3

are given in (29), ( 33) and (34) respectively. Therein, Laplace transforms and other related function definitions are given below C

LI ss (s) = exp -λ sc R 2 1 -exp(-cν(s, x)) dx R 2 exp(-cν(s, x)) f (x)dx (36) ν(s, x) = R 2 f (y) 1 + (sPsc ˜ (y -x)) -1 dy (37) LI ms (s) = exp -λmc (sPmc) 2/α π 2 (2/α) sin(π 2 α ) (38)
C (cap) 2 = 1 -exp - τmcλcr µγλut-s (39)
C (cap) 3 = 1 -1 + Fsc 1-η (40)
V. NUMERICAL RESULTS AND DISCUSSION

In this section, we validate our expressions via Monte-Carlo simulations. First, the impact of backhaul rate splitting ratio γ on the average delivery rate is given in Figs. (3a) and (4a) for coverage and capacity-aided deployments respectively. Indeed, as seen from the figures, a dramatical increase in average delivery rate occurs which confirms our intuitions. The rate of decrement for the typical SU is relatively slow compared to the typical MU. In order to find a balance between average delivery rate of typical MUs and SUs, ensuring rate fairness, the plots show that one has to set γ carefully. In all of these cases, we observe that having caching capabilities at SBSs improves system performance in terms of average delivery rate. In other words, a heterogeneous network consists of MBSs and cache-enabled SBSs allows higher average delivery rates while ensuring fairness between users at different tiers.

Second, the impact of storage size F sc on the average delivery rate is given in Figs. (3b) and (4b) for coverage and capacity-aided deployments respectively. It is shown that increasing storage size of SBSs both in coverage and capacityaided deployments yields higher average delivery rates. The increment of storage size in coverage-aided deployment is more visible compared to capacity-aided deployment and allows typical SUs to achieve higher rates than typical MUs. Given the fact that caching monotonically improves the overall system performance of SUs, more storage does not seem necessary if one considers a linear cost for storage size, even though more storage is desirable for improving the average delivery rate.

Figure 1 :

 1 Figure 1: An illustration of the coverage-aided deployment.

Figure 2 :

 2 Figure 2: An illustration of the capacity-aided deployment.

  ss = xi∈Φsc\{x} P sc g xi (x i ) is the cumulative interference from other clustered SBSs except the serving cell at x, and I sm = yi∈Φsc P mc h yi (y) is the total interference from MBSs.

  given in (14) and (18) respectively. Therein, Laplace transforms and other related function definitions are given below B (cov) 1

Remark 1 .

 1 in (19) and (24) respectively.Proof. See Appendix C.1 in[START_REF] Baştug | Distributed caching methods in small cell networks[END_REF]. These expressions are cumbersome but numerically easy to compute. The terms B 1 and B 2 capture the downlink and backhaul behaviour respectively. e -(rmc/ν) k drmc (14)

Remark 2 .

 2 in (35), (39) and (40) respectively.Proof. See Appendix C.2 in[START_REF] Baştug | Distributed caching methods in small cell networks[END_REF]. The terms C 1 , C 2 and C 3 incorporate the downlink, backhaul and caching aspects respectively.

F 1 , 1 - 2 ( 2

 1122 e -(rsc/ν) k drsc (29) LI ss (s) = exp -sπλ sc Psc(2/α) 1 -2/α r 2-α sc × 2/α; 2 -2/α; -sPscr -α sc (30) LI ms (s) = exp -λmc (sPmc) 2/α π e -(rsc/ν) k drsc (35)

Figure 3 :Figure 4 :

 34 Figure 3: Evolution of average delivery rate in coverage-aided deployment. λcr = 1.0×10 -5 , λmc = 1.5×10 -5 , λ sc = 5.5×10 -5 , λut = 12.8×10 -5 unit/m 2 ; Pmc = 16, Psc = 3 Watt; τmc = τsc = 4 bits/s/Hz; α = 4 ; Rc = 80 meters; µ = 30 bits/s/Hz; γ = 0.6; f0 = 500 GByte; Fsc = 4 GByte; η = 1.45.

  )

	100
	50
	0
	-50
	-100
	-100	-50		0	50	100
	CR	MBS	SBS	exclusion boundary	mobile user


Senior Member, IEEE, Mehdi Bennis, Senior Member, IEEE Walid Saad

Fellow, IEEE Merouane Debbah

Unmanned Aerial Vehicle with Underlaid Device-to-Device Communications: Performance and Tradeoffs

In this paper, the deployment of an unmanned aerial vehicle (UAV) as a flying base station used to provide on the fly wireless communications to a given geographical area is analyzed. In particular, the co-existence between the UAV, that is transmitting data in the downlink, and an underlaid device-todevice (D2D) communication network is considered. For this model, a tractable analytical framework for the coverage and rate analysis is derived. Two scenarios are considered: a static UAV and a mobile UAV. In the first scenario, the average coverage probability and the average sum-rate for the users in the area are derived as a function of the UAV altitude and the number of D2D users. In the second scenario, using the disk covering problem, the minimum number of stop points that the UAV needs to visit in order to completely cover the area is computed. Simulation and analytical results show that, depending on the density of D2D users, optimal values for the UAV altitude exist for which the average sum-rate and the coverage probability are maximized. Moreover, our results also show that, by enabling the UAV to intelligently move over the target area, the overall communication rate and coverage probability can be significantly improved. Finally, in order to provide a full coverage for the area of interest, the tradeoff between the coverage and delay, in terms of the number of stop points, is discussed.

. Last but not least, in regions or countries in which building a complete cellular infrastructure is very expensive, deploying UAVs is highly beneficial as it removes the need for towers and cables. In order to reap the benefits of UAV deployments for communication purposes, one must address a number of technical challenges that include performance analysis, channel modeling, optimal deployment, and resource management, among others [5]-[15]. The most significant existing body of work on UAV communications focuses on air-to-ground channel modeling [5]-[8]. For instance, in [5] and [6], the probability of line of sight (LOS) for air-to-ground communication as a function of the elevation angle and average height of buildings in a dense urban area was derived. The air-to-ground path loss model has been further studied in [7] and [8]. As discussed in [8], due to path loss and shadowing, the characteristics of the air-to-ground channel are shown to depend on the height of the aerial base stations.

 considered the use of UAVs to compensate for the cell overload and outage in cellular networks.

, the authors investigated how to optimally move UAVs for

I. INTRODUCTION

The use of unmanned aerial vehicles (UAVs) as flying base stations that can boost the capacity and coverage of existing wireless networks has recently attracted significant attention [START_REF] Bucaille | Rapidly deployable network for tactical applications: Aerial base station with opportunistic links for unattended and temporary events absolute example[END_REF] and [START_REF] Zhan | Wireless relay communications with unmanned aerial vehicles: Performance and optimization[END_REF]. One key feature of a UAV that can potentially lead to the coverage and rate enhancement is having line-of-sight (LOS) connections towards the users. Moreover, owing to their agility and mobility, UAVs can be quickly and efficiently deployed to support cellular networks and M. Mozaffari and W. Saad are with Wireless@VT, Department of ECE, Virgina Tech, Blacksburg, VA, USA. Emails: {mmozaff,walids}@vt.deu. M. Bennis is with CWC -Centre for Wireless Communications, Oulu, Finland, Email: bennis@ee.oulu.fi. M. Debbah is with Mathematical and Algorithmic Sciences Lab, Huawei France R & D, Paris, France, Email:merouane.debbah@huawei.com. improving connectivity of ad hoc networks. However, [START_REF] Han | Optimization of MANET connectivity via smart deployment/movement of unmanned air vehicles[END_REF] only focused on an ad-hoc network and assumed that the UAV have complete information about the location of nodes. In [START_REF] Jiang | Optimization of UAV heading for the ground-to-air uplink[END_REF], considering static ground users, the optimal trajectory and heading of UAVs equipped with multiple antennas for ground to air uplink scenario was derived.

For scenarios in which there is limited or no infrastructure support, beyond the use of UAVs, there has been considerable recent works that study the use of direct device-to-device (D2D) communications between wireless users over the licensed spectrum [START_REF] Yaacoub | Energy-efficient device-to-device communications in LTE public safety networks[END_REF]. Such D2D communications has been shown improve coverage and capacity of existing wireless networks, such as cellular systems. In particular, in hotspot areas or public safety scenarios, D2D will allow users to communicate directly with one another without significant infrastructure. D2D communications are typically deployed using underlaid transmission links which reuse existing licensed spectrum resources [START_REF] Doppler | Device-to-device communication as an underlay to LTE-advanced networks[END_REF]. Therefore, deploying a UAV over a spectrum band that must be shared with an underlaid D2D network will introduce important interference management challenges. In the literature, there are some studies on the coexistence of the underlaid D2D and cellular communications with a single base station [START_REF] Lee | Power control for D2D underlaid cellular networks: Modeling, algorithms, and analysis[END_REF]. Furthermore, the authors in [START_REF] Shalmashi | Energy efficiency and sum rate tradeoffs for massive MIMO systems with underlaid device-to-device communications[END_REF] and [START_REF] Lin | The interplay between massive MIMO and underlaid D2D networking[END_REF] exploited the interplay between the massive MIMO and underlaid D2D communications in a single cell. However, none of theses prior works studied the coexistence of UAVs and underlaid D2D communications. In particular, a comprehensive analytical analysis to evaluate this coexistence in terms of different performance metrics, such as coverage and rate, is lacking in the current state-of-the-art [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF], [START_REF] Han | Optimization of MANET connectivity via smart deployment/movement of unmanned air vehicles[END_REF], [START_REF] Lee | Power control for D2D underlaid cellular networks: Modeling, algorithms, and analysis[END_REF]- [START_REF] Lin | The interplay between massive MIMO and underlaid D2D networking[END_REF].

Compared to the previous studies on the coexistence of D2D and cellular networks such as [START_REF] Shalmashi | Energy efficiency and sum rate tradeoffs for massive MIMO systems with underlaid device-to-device communications[END_REF] and [START_REF] Lin | The interplay between massive MIMO and underlaid D2D networking[END_REF], the presence of an aerial UAV base station along with D2D links introduces new challenges. First, the channel modeling between the UAV and ground users will no longer be a classical fading channel, instead, it will be based on probabilistic LOS and NLOS links [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF], [START_REF] Feng | Modelling the likelihood of line-of-sight for air-to-ground radio propagation in urban environments[END_REF], while the channel between a base station and the users will still follow a Rayleigh fading model. Second, unlike conventional, fixed base stations, the height of a UAVs is adjustable and this impacts the channel characteristics and the coverage performance. Third, the potential mobility of a UAV introduces new dimensions to the problem and the impact of such mobility on D2D and network performance must be analyzed. The prior studies on UAVs such as [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF]- [START_REF] Han | Optimization of MANET connectivity via smart deployment/movement of unmanned air vehicles[END_REF] have not addressed the third challenge. More specifically, the interplay between UAVs and D2D communications and the existing challenges and tradeoffs have not been investigated in these literature. To our best knowledge, this paper will provide the first comprehensive fundamental analysis on the performance of UAV communication in the presence of underlaid D2D links.

The main contribution of this paper is to analyze the coverage and rate performance of UAV-based wireless communication in the presence of underlaid D2D communication links.

In particular, we consider a network in which a single UAV must provide downlink transmission support to a number of users within a given area. In this area, a subset of the devices is also engaged in D2D transmissions that operate in an underlay fashion over the UAV's transmission.

We consider two types of users, namely downlink users (DUs) which receive data from the UAV, and D2D users which communicate directly with one another. Here, the UAV must communicate with the DUs while taking into account the potential interference stemming from the underlaid D2D transmissions. For this network, we analyze two key cases: static UAV and mobile UAV. Using tools from stochastic geometry, for both scenarios, we derive the average downlink coverage probabilities for DUs and D2D users and we analyze the impact of the UAV altitude and density of the D2D users on the overall performance. For the static case, we find the optimal values for the UAV altitude which leads to a maximum coverage probability for DUs.

In addition, considering both DUs and D2D users, an optimal altitude which maximizes the average sum-rate is computed. Our results demonstrate that the optimal UAV altitude decreases as the density of D2D users increases. The results show that a maximum average sum-rate can be achieved if the UAV altitude is appropriately adjusted based on the D2D users density.

Furthermore, for a given UAV altitude, we show that an optimal value for the number of D2D users that maximizes the average sum-rate exists.

For the mobile UAV case, we assume that the UAV can travel over the area while stopping at some given locations in order to serve the downlink users. Considering retransmissions at different time instances, we derive the overall coverage probability. Then, using the disk covering problem, we find a minimum number of stop points that the UAV needs to to completely cover the area. This can be interpreted as the fastest way to cover the whole area with a minimum required transmit power. In addition, we analyze the tradeoff between the number of stop points, which is considered as delay here, and the coverage probability for the downlink users. We show that, in order to enhance the coverage for DUs, the UAV should stop in more locations over the target area which can, in turn, lead an increased delay. For example, our results show that for a given density of D2D users, to increase the DU coverage probability from 0.4 to 0.7, the number of stop points should be increased from 5 to 23. Furthermore, the number of stop points is shown to significantly depend on the number of D2D users. For instance, if the average number of D2D users in the area increases from 50 to 100, in order to maintain the DUs' coverage requirement, the number of stop points should be increased from 20 to 55. Finally, we prove that the overall coverage probability for both DUs and D2D users can be improved by moving the UAV.

The rest of this paper is organized as follows. Section II presents the system model and describes the air-to-ground channel model. In Section III, coverage probabilities for DUs and D2D users are provided for a single static UAV. Section IV presents the performance evaluation for one mobile UAV which is used to provide full coverage for the target area. Section V presents the simulation results while Section VI draws some conclusions.

II. SYSTEM MODEL

Consider an area with a radius R c in which a number of users are spatiality distributed according to a Poisson point processes (PPP) [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF], and a UAV (at low altitude platform) is used to serve a subset of those users. In this network, the users are divided into two groups: downlink users located based on a PPP Φ A with density λ du (number of users per m 2 ) and D2D users whose distribution follows a PPP Φ B with density of λ d (number of D2D pairs per m 2 ). Note that, the average number of users in a given area is equal to the density of the users multiplied by the size of the area. Here, we focus on the downlink scenario for the UAV and we assume that the D2D users communicate in an underlay fashion. Furthermore, we assume that a D2D receiver connects to its corresponding D2D transmitter pair located at a fixed distance away from it in an isotropic direction [START_REF] Lee | Power control for D2D underlaid cellular networks: Modeling, algorithms, and analysis[END_REF]. Therefore, the received signals at the D2D receiver include the desired signal from the D2D transmitter pair and interference from the UAV and other D2D transmitters. A downlink user, on the other hand, receives the desired signal from the UAV but it also experiences interference from all the D2D transmitters. For tractability as discussed in [START_REF] Shalmashi | Energy efficiency and sum rate tradeoffs for massive MIMO systems with underlaid device-to-device communications[END_REF],

we also consider the interference from D2D transmitters located outside the area with the radius of R c . This assumption removes the concern stemming from the boundary effect in which users located at the cell boundary receives less interference than those who are closer to the center. However, we only evaluate the coverage and rate performance of users located inside the area.

The signal to interference plus noise ratio (SINR) expression for a D2D receiver is

γ d = P r,d I c d + I u + N , ( 1
)
where P r,d is the received signal power from the D2D transmitter, I c d is the total interference from other D2D users, I u is the interference from the UAV, and N is the noise power. Moreover, we have:

P r,d = P d d -α d 0 g 0 , ( 2
)
I c d = ∑ i̸ =0 P d d i -α d g i , (3)
I d = ∑ i P d d i -α d g i , (4)
where the index i = 0 is used for the selected D2D transmitter/receiver pair, g 0 and g i are,

respectively, the channel gains between a D2D receiver and its corresponding D2D transmitter, and the i th interfering D2D transmitters. For the D2D transmission, we assume a Rayleigh fading channel model [START_REF] Lee | Power control for D2D underlaid cellular networks: Modeling, algorithms, and analysis[END_REF], [START_REF] Lin | The interplay between massive MIMO and underlaid D2D networking[END_REF] and [START_REF] Afshang | Modeling and performance analysis of clustered device-to-device networks[END_REF]. P d is the D2D transmit power which is assumed to be fixed and equal for all the users, d i is the distance between a D2D receiver and the i th D2D transmitter, d 0 is the fixed distance between the D2D receiver and transmitter of the selected D2D pair, and α d is the path loss exponent between D2D users. Note that the received signal powers as well as the noise power are normalized by a path loss coefficient.

The SINR expression for a DU user that connects to the UAV is

γ u = P r,u I d + N , ( 5
)
where P r,u is the received signal power from the UAV and I d c is the total interference power from D2D transmitters.

A. Air-to-ground channel model

As discussed in [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF] and [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF], the ground receiver receives three groups of signals including LOS, strong reflected non-line-of-sight (NLOS) signals, and multiple reflected components which cause multipath fading. These groups can be considered separately with different probabilities of occurrence as shown in [START_REF] Holis | Elevation dependent shadowing model for mobile communications via high altitude platforms in built-up areas[END_REF] and [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF]. Typically, it is assumed that the received signal is categorized in only one of those groups [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF]. Each group has a specific probability of occurrence which is a function of environment, density and height of buildings, and elevation angle. Note that the probability of having the multipath fading is significantly lower than the LOS and NLOS groups [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF]. Therefore, the impact of small scale fading can be neglected in this case [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF].

One common approach to modeling air-to-ground propagation channel is to consider LOS and NLOS components along with their occurrence probabilities separately as shown in [START_REF] Hourani | Modeling air-to-ground path loss for low altitude platforms in urban environments[END_REF] and [START_REF] Holis | Elevation dependent shadowing model for mobile communications via high altitude platforms in built-up areas[END_REF].

Note that for NLOS connections due to the shadowing effect and the reflection of signals from obstacles, path loss is higher than in LOS. Hence, in addition to the free space propagation loss, different excessive path loss values are assigned to LOS and NLOS links. Depending on the LOS or NLOS connection between the user and UAV, the received signal power at the user location is given by [9]

P r,u =    P u |X u | -αu LOS connection, ηP u |X u | -αu NLOS connection, (6)
where P u is the UAV transmit power, |X u | is the distance between a generic user and the UAV, α u is the path loss exponent over the user-UAV link, and η is an additional attenuation factor due to the NLOS connection. Here, the probability of LOS connection depends on the environment, density and height of buildings, the location of the user and the UAV, and the elevation angle between the user and the UAV. The LOS probability can be expressed as follows [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF]:

P LOS = 1 1 + C exp(-B [θ -C]) , ( 7
)
where C and B are constant values which depend on the environment (rural, urban, dense urban, or others) and θ is the elevation angle. Clearly, θ = 180 π × sin -1

( h |Xu| ) , |X u | = √ h 2 + r 2 and
also, probability of NLOS is P NLOS = 1 -P LOS .

As observed from [START_REF] Feng | Path loss models for air-to-ground radio channels in urban environments[END_REF], the LOS probability increases as the elevation angle between the user and UAV increases.

Given this model, we will consider two scenarios: a static UAV and a mobile UAV. For each scenario, we will derive the coverage probabilities and average rate for DUs and D2D users.

Once those metrics are derived, considering the D2D users density, we obtain optimal values for the UAV altitude that maximize the coverage probability and average rate.

III. NETWORK WITH A STATIC UAV

In this section, we evaluate the coverage performance of the scenario in which one UAV located at the altitude of h in the center of the area to serve the downlink users in the presence of underlaid D2D communications. Clearly, in such a scenario, considering the uniform distribution of users over the area, placing the UAV in the center of the cell is an optimal deployment.

A. Coverage probability for D2D users

Consider a D2D receiver located at (r, φ), where r and φ are the radius and angle in a polar coordinate system assuming that the UAV is located at the center of the area of interest. Note that considering ( 6) and ( 7), the coverage probability for a user located at (r, φ) is also a function of the UAV altitude, h. In this case, the coverage probability can be derived as follows:

Theorem 1. For underlay D2D communication, the coverage probability for a D2D receiver connecting to the D2D transmitter located at a fixed distance away from it is given by:

P cov,d (r, φ, β) = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βd α d 0 N P d ) × ( P LOS (r) exp ( -βd α d 0 P u r -αu P d ) + P NLOS (r) exp ( -βd α d 0 ηP u r -αu P d )) . ( 8
)
Proof.

P cov,d (r, φ, β) = P [γ d ≥ β] = P [ P d d -α d 0 g I c d + I u + N ≥ β ] = P [ g ≥ βd α d 0 (I c d + I u + N ) P d ] (a) = E Iu,I c d [ exp( -βd α d 0 (I c d + I u + N ) P d ) ] (b) = E Iu [ exp( -βd α d 0 I u P d ) ] E I c d [ exp( -βd α d 0 I c d P d ) ] exp ( -βd α d 0 N P d ) , ( 9
)
where g is an exponential random variable with a mean value of one (i.e. g ∼ exp(1)), (a)

follows from the exponential distribution of g based on the Rayleigh fading assumption, and taking the expectation over I u and I c d (as random variables).

Step (b) comes from the fact that I u and I c d are independent because the interference stems from different sources which are spatially uncorrelated.

Here, E Iu and E I c d are given by:

E Iu [ exp( -βd α d 0 I u P d ) ] = P LOS (r) exp ( -βd α d 0 P u r -αu P d ) + P NLOS (r) exp ( -βd α d 0 ηP u r -αu P d ) , ( 10
)
E I c d [ exp( -βd α d 0 I c d P d ) ] = E d i ,g i [ ∏ i exp( -βd α d 0 P d P d d i -α d g i ) ] = E d i [ ∏ i E g i [ exp(-βd α d 0 d i -α d g i ) ] ] (a) = exp   -2πλ d ∞ ∫ 0 ( 1 -E g [ exp(-βd α d 0 r -α d g) ]) rdr   (b) = exp   -2πλ d ∞ ∫ 0   1 - ∞ ∫ 0 exp(-gβd α d 0 r -α d -g)dg   rdr   = exp   -2πλ d ∞ ∫ 0 βd α d 0 r 1-α d 1 + βd α d 0 r -α d dr   = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) ) , ( 11
)
where in (a) we used the probability generating functional (PGFL) of PPP. Note that for a point process Φ the PGFL is defined as

PGFL = E [ ∏ x∈Φ f (x) ] . ( 12
)
For a PPP with intensity λ the PGFL is equal to exp

( -λ ∫ S [1 -f (x)]dx
) .

Also, the second step (b) is based on the exponential distribution of the channel gain (∼ exp(1)).

Finally, using ( 9), ( 10) and ( 11) Theorem I is proved.

From this theorem, we can make several key observations. First, given that the UAV is at the center of the target area, as r or equivalently the distance of a D2D user from the UAV increases, the D2D coverage probability in [START_REF] Holis | Elevation dependent shadowing model for mobile communications via high altitude platforms in built-up areas[END_REF] increases. This is because the interference power from the UAV is lower at higher distances and hence the D2D users located at the cell (target area) boundary have higher coverage probability than those which are closer to the center. Second, the D2D coverage probability in [START_REF] Holis | Elevation dependent shadowing model for mobile communications via high altitude platforms in built-up areas[END_REF] decreases when the UAV transmit power increases. To cope with this situation, the D2D users can increase their transmit power or reduce the fixed distance parameter (D). In addition, decreasing the D2D user density improves the coverage probability due to decreasing the interference.

Note that the result presented in Theorem I corresponds to the coverage probability for a D2D user located at (r, φ). To compute the average coverage probability in the cell, we consider a uniform distribution of users over the area with

f (r, φ) = r πR 2 c , 0 ≤ r ≤ R c , 0 ≤ φ ≤ 2π 1 ,
and we find the average over the area. Then, the average coverage probability for D2D users will be

Pcov,d (β) = E r, φ [P cov,d (r, φ, β)] = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βd α d 0 N P d ) Rc ∫ 0 E Iu [ exp( -βd α d 0 I u P d ) ] f (r, φ)drdφ = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βd α d 0 N P d ) Rc ∫ 0 E Iu [ exp( -βd α d 0 I u P d ) ] 2r R 2 c dr. (13)
From ( 13), we can see that the average coverage probability for D2D users increases as the size of the area , R c , increases. In fact, when the UAV serves a larger area, the average distance of D2D users from the UAV increases and on the average they receive lower interference from it. Next, we provide a special case for [START_REF] Rohde | Interference aware positioning of aerial relays for cell overload and outage compensation[END_REF] in which the UAV has a very high altitude or very small transmit power.

Corollary 1. For P u = 0 or h → ∞, the average coverage probability for the D2D users is simplified to

Pcov,d (β) = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βd α d 0 N P d ) , ( 14
)
Note that, the result in Corollary 1 corresponds to the coverage probability in overlay D2D communication in which there is no interference between the UAV and the D2D transmitters.

B. Coverage Probability for Downlink Users

In this section, we first derive an approximation for the downlink users' coverage probability.

Theorem 2. The average coverage probability for DUs in the cell is approximated as

Pcov,du (β) ≈ Rc ∫ 0 P LOS (r)A I ( P u X u -αu β -N ) 2r R 2 c dr + Rc ∫ 0 P NLOS (r)A I ( ηP u X u -αu β -N ) 2r R 2 c dr, (15)
where for T > 0,

A I (T ) = ( 1 -πλ d Γ(1+2/α d ) α d -2 ( T P d ) -2/α d ) exp ( -πλ d ( T KP d ) -2/α d Γ(1 + 2/α d )
) .

Also,

Γ(t) = ∞ ∫ 0 x t-1 e -x
dx is the gamma function [START_REF] Artin | The gamma function[END_REF].

Proof. The coverage probability for a cellular user located at (r, φ) is written as

P cov,du (r, φ, β) = P [γ u ≥ β] = P LOS (r)P [ P u r -αu I d + N ≥ β ] + P NLOS (r)P [ ηP u r -αu I d + N ≥ β ] = P LOS (r)P [ I d ≤ P u r -αu -βN β ] + P NLOS (r)P [ I d ≤ ηP u r -αu -βN β ] . (16)
Note that there is no closed-form expression for the cumulative distribution function (CDF) of the interference from D2D users [START_REF] Ganti | A stochastic geometry approach to the interference and outage characterization of large wireless networks[END_REF] and [START_REF] Weber | Transmission capacity of wireless ad hoc networks with outage constraints[END_REF]. Here, we provide lower and upper bounds for the CDF of interference. First, we divide the interfering D2D transmitters into two subsets:

   Φ 1 = {Φ B |P d d i -α d g i ≥ T }, Φ 2 = {Φ B |P d d i -α d g i ≤ T }, ( 17
)
where T is a threshold which is used to derive the CDF of the interference from D2D users. Now, considering the interference power from D2D users located in Φ 1 and Φ 2 as I d,Φ 1 and

I d,Φ 2 , we have P [I d ≤ T ] = P [I d,Φ 1 + I d,Φ 2 ≤ T ] ≤ P [I d,Φ 1 ≤ T ] = P [Φ 1 = 0] = E [ ∏ Φ B P(P d d i -α d g i < T ) ] = E [ ∏ Φ B P(g i < T d i α d P d ) ] (a) = P [ ∏ Φ B 1 -exp(- T d i α d P d ) ] (b) = exp   -λ d ∞ ∫ 0 exp(- T r α d P d )rdr   = exp ( -πλ d ( T P d ) -2/α d Γ(1 + 2/α d ) ) , ( 18
)
where (a) and (b) come from the Rayleigh fading assumption and PGFL of the PPP.

The upper bound is derived as follows:

P [I d ≤ T ] = 1 -P [I d ≥ T ] = 1 - ( P [I d ≥ T |I d,Φ 1 ≥ T ] P [I d,Φ 1 ≥ T ] + P [I d ≥ T |I d,Φ 1 ≤ T ] P [I d,Φ 1 ≤ T ] ) = 1 - ( P [I d,Φ 1 ≥ T ] + P [I d ≥ T |I d,Φ 1 ≤ T ] P [I d,Φ 1 ≤ T ] ) = 1 - ( 1 -P [Φ 1 = 0] + P [I d ≥ T |I d,Φ 1 ≤ T ] P [Φ 1 = 0] ) = P [Φ 1 = 0] ( 1 -P [I d ≥ T |Φ 1 = 0] ) . ( 19
)
Also,

P [I d ≥ T |Φ 1 = 0] (a) ≤ E [I d ≥ T |Φ 1 = 0] T = 1 T E [ ∑ Φ P d d i -α d g i 1(P d d i -α d g i ≤ T ) ] = 1 T E d i [ ∑ Φ P d d i -α d E g i [ g i 1(g i ≤ T d i α d P d ) ] ] = 1 T E d i      ∑ Φ P d d i -α d      T d i α d P d ∫ 0 ge -g dg           = 2πP d λ d T ∞ ∫ 0 r -α d      T r α d P d ∫ 0 ge -g dg      rdr = 2πλ d Γ(1 + 2/α d ) α d -2 ( T P d ) -2/α d . ( 20
)
where (a) is based on the Markov's inequality which is stated as follows: for any non-negative integrable random variable X and positive L, P (X ≥ L) ≤ E[X] L . Also, 1(.) is the indicator function which can only be equal to 1 or 0. Hence, the lower (L I ) and upper (U I ) bounds for the CDF of interference become

L I (T ) = ( 1 - 2πλ d Γ(1 + 2/α d ) α d -2 ( T P d ) -2/α d ) exp ( -πλ d ( T P d ) -2/α d Γ(1 + 2/α d ) ) , ( 21
)
U I (T ) = exp ( -πλ d ( T P d ) -2/α d Γ(1 + 2/α d ) ) . ( 22
)
Thus, we have

L I (T ) ≤ P{I d ≤ T } ≤ U I (T ).
Here, for simplicity, we approximate P{I d ≤ T } with the average of its lower and upper bounds:

P{I d ≤ T } ≈ L I (T ) + U I (T ) 2 = A I (T ). (23)
Finally, using ( 15) and ( 23), the average coverage probability for the DUs is found as per Theorem 2.

From Theorem 2, we can first see that, for T >> P d , given that e -x ≈ 1 -x when x → 0,

we have U I (T ) = L I (T ) ≈ 1 -πλ d ( T P d ) -2/α d Γ(1 + 2/α d ).
This means that the approximation in ( 23) becomes tighter for lower transmit power of D2D users. Moreover, from ( 21) and ( 22), when λ d → ∞, the number of D2D users tends to infinity and U I = L I = 0. Consequently, the downlink users experience an infinite interference from the D2D users which results in Pcov,du = 0.

As per Theorem 2, increasing R c decreases the average coverage probability for the downlink users. However, higher R c results in a higher D2D average coverage probability. Moreover, the average coverage probability for downlink users decreases as the density of the D2D users increases. In this case, to improve the DUs coverage performance, one must increase P u or reduce R c . Next, we derive the DU coverage probability in the absence of the D2D users.

Proposition 1. Assuming there is no interference from D2D users, we have P d = 0, and, then, the average coverage probability for the downlink users can be expressed by

Pcov,du (β) = ∫ min[( Pu βN ) 1/αu ,Rc] 0 P LOS (r) 2r R 2 c dr + ∫ min[( ηPu βN ) 1/αu ,Rc] 0 P NLOS (r) 2r R 2 c dr. ( 24
)
Proof. For a DU located at (r, φ), the coverage probability in absence of D2D users becomes

P cov,du (r, φ, β) = P [γ u ≥ β] = P LOS (r)P [γ u ≥ β|LOS] + P NLOS (r)P [γ u ≥ β|NLOS] = P LOS (r)1 [ r ≤ ( P u βN ) 1/αu ] + P NLOS (r)1 [ r ≤ ( ηP u βN ) 1/αu ] , ( 25
)
The average coverage probability is computed by taking the average of P cov,du (r, φ, β) over the cell with the radius R c .

P cov,du (r, φ, β) = E r, φ [P cov,du (r, φ, β)] = ∫ min[( Pu βN ) 1/αu ,Rc] 0 P LOS (r) 2r R 2 c dr + ∫ min[( ηPu βN ) 1/αu ,Rc] 0 P NLOS (r) 2r R 2 c dr. ( 26
)

C. Average sum-rate

Now, we investigate the average achievable rates for the DUs and D2D users which can be expressed as in [START_REF] Shalmashi | Energy efficiency and sum rate tradeoffs for massive MIMO systems with underlaid device-to-device communications[END_REF]:

Cdu = W log 2 (1 + β) Pcov,du (β), (27)
Cd = W log 2 (1 + β) Pcov,d (β), (28)
where W is the transmission bandwidth. Considering the whole DUs and D2D users in the cell, the average sum-rate, Csum , can be derived as a function of the coverage probabilities and the number of users as follows:

Csum = R c 2 πλ du Cdu + R c 2 πλ d Cd . ( 29
)
Assuming µ = λ du λ d , we have

Csum = λ d R c 2 π [ µ Pcov,du (β) + Pcov,d (β) ] W log 2 (1 + β), ( 30
)
where R c 2 πλ d and R c 2 πλ du are the number of DUs and D2D users in the target area respectively.

From (30), observe that, on the one hand, Csum is directly proportional to λ d , but on the other hand, it depends on the coverage probabilities of DUs and D2D users which both are decreasing functions of D2D user density. Therefore, in general increasing λ d does not necessarily enhance the rate. Note that, considering (13), ( 16) and (30), for both λ d → 0 and λ d → ∞ cases the average sum-rate tend to zero. Hence, there is an optimum value for λ d that maximizes Csum .

According to (28), Csum is a function of the coverage probability and a logarithmic function of the threshold (β). The former is a decreasing function of β whereas the latter is an increasing function of β. In other words, although increasing the threshold is desirable for the rate due to increasing the logarithmic function, it also reduces the coverage probability. Therefore, in order to achieve a maximum rate, a proper value for the threshold must be derived.

IV. NETWORK WITH A MOBILE UAV

Now, we assume that the UAV can move around the area of radius R c in order to provide coverage for all the downlink users in the target area. In particular, we consider a UAV that moves over the target area and only transmits at a given geographical location (area) which we hereinafter refer to as "stop points". Each stop point represents a location over which the UAV stops and serves the present downlink users. Here, our first goal is to minimize the number of stop points (denoted by M ) and determine their optimal location. Note that, as the UAV moves, it can have a different channel to a user at different time instances. The objective of the UAV is to cover the entire area and ensure that the coverage requirements for all DUs are satisfied with a minimum UAV transmit power and minimum number of stop points. In other words, we find the minimum number and location stop points for the UAV to completely cover the area. We model this problem by exploiting the so-called disk covering problem [START_REF] Kershner | The number of circles covering a set[END_REF]. In the disk covering problem, given a unit disk, the objective is to find the smallest radius required for M equal smaller disks to completely cover the unit disk. In the dual form of the problem, for a given radius of small disks, the minimum number of disks required to cover the unit disk is found.

In Figure 2, we provide an illustrative example to show the mapping between the mobile UAV communication problem and the disk covering problem. In this figure, the center of small disks can be considered as the location of stop points and the radius of the disk is the coverage radius of the UAV. Using the disk covering problem analysis, in Table I, we present, for different number of stop points, the minimum required coverage radius of a UAV for completely covering the target area [START_REF] Kershner | The number of circles covering a set[END_REF], [START_REF] Tóth | Thinnest covering of a circle by eight, nine, or ten congruent circles[END_REF]. Thereby, using the dual disk covering problem, for a given maximum coverage radius of a UAV, we can find the minimum number of stop points for covering the entire area. The detailed steps for finding the minimum number of stop points are provided next.

First, we compute the coverage radius of the UAV based on the minimum requirement for the DU coverage probability. The coverage radius is defined as the maximum radius within which the coverage probability for all DUs (located inside the coverage range) is greater than a specified threshold, ϵ. In this case, the UAV satisfies the coverage requirement of each DU which is inside its coverage range. The maximum coverage radius for the UAV at an altitude h transmitting with a power P u will be given by:

R m = max{R|P cov,A (β, R) ≥ ε, P u , h} = P -1 cov,A (β, ε), ( 31
)
where ε is the threshold for the average coverage probability in the cell (area covered by the UAV). Note that a user is considered to be in coverage if it is in the coverage range of the UAV.

The minimum required number of stop points for the full coverage is

{ L = min{M }, P cov,du (r, φ, β) ≥ ε, ( 32
)
where M represents the number of stop points, the second condition guarantees that the area is completely covered by the UAV, and L is the minimum value for the number of stop points if the following condition holds:

R min,L ≤ R m ≤ R min,L-1 → min{M } = L. ( 33
)
By using Table I, we see that, R min,L-1 and R min,L are, respectively, the minimum radius required to cover the entire target area with L -1 and L disks. After finding the minimum M , we can reduce the UAV transmission power such that the coverage radius decreases to the minimum required radius (R min,L ). In this way, the UAV transmit power is minimized. Thus we have

P u,min = argmin Pu {P -1 cov,du (β, ε) = R min,L |h}, ( 34
)
where P u,min is the minimum UAV transmit power. Thereby, the minimum number of stop points leads to a full coverage at a minimum time with a minimum required transmit power.

In summary, the proposed UAV deployment method that leads to the complete coverage with a minimum time and transmission power proceeds as follows. First, depending on the parameters of the problem such as density of users and threshold, we compute the maximum coverage radius of a UAV at the optimal altitude that can serve the DUs. Second, considering the size of target area, using the disk covering problem, we find the minimum required number of transmission points along with the coverage radius at each point. Third, we reduce the transmission power of UAV such that its maximum coverage radius becomes equal to the required coverage radius found in the previous step. Using the proposed method, the target area can be completely covered by the UAV with a minimum required transmit power and minimum number of stop points.

Coverage radius of the UAV

) M = 1, 2 R c M = 3 √ 3 2 R c M = 4 √ 2 2 R c M = 5 0.61R c M = 6 0.556R c M = 7 0.5R c M = 8 0.437R c M = 9 0.422R c M = 10 0.398R c M = 11 0.38R c M = 12 0.361R c
Next, we derive the overall coverage probability for a typical D2D user in the M time instances for the mobile UAV and the static UAV cases. In other words, we consider the network in M time instances in which the UAV and D2D users have M retransmissions, and compare the overall achievable coverage performance for the D2D users in the mobile UAV and static UAV scenarios.

Assume that the relative location of the i th stop point with respect to the D2D user is (r i , h i )

where r i is the distance between the projection of the UAV on the ground and D2D user and h i is the UAV altitude. Clearly, the distance between the user and UAV is

|X u,i | = √ h 2 i + r i 2 .
As proved in Theorem 1, the coverage probability at the i th time instance or i th stop point is

P i cov,d (β) = exp ( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βD α d N KP d ) × E i Iu [ exp( -βd α d 0 I u KP d ) ] , ( 35
)
where

E i Iu [ exp( -βd α d 0 Iu KP d ) ] = P LOS,i (r i ) exp ( -βD α d Pu|X u,i | -αu P d ) + P NLOS,i (r i ) exp ( -βd α d 0 ηPu|X u,i | -αu P d
) ,

and

P LOS,i = 1 1+C exp ( -B [ 180 π ×sin -1 ( h |X u,i | ) -C ]) .
The overall coverage probability for a D2D user after M retransmissions assuming the UAV location is different in different retransmission times, is

P O,m cov,d (β) = 1 - M ∏ i=1 ( 1
-P i cov,d (β) ) . ( 36
)
Next, we derive the overall coverage probability for D2D users when the UAV is static. Similar to the dynamic UAV case, we consider M number of retransmissions at different time instances.

Theorem 3. The overall D2D coverage probability in M retransmissions considering the static UAV case is given by

P O,s cov,d (β) = P 2 × [ 1 -(1 -P 1,i ) M ] , ( 37
)
where P 1,i = exp

( -2π 2 λ d β 2/α d d 2 0 α d sin(2π/α d ) - βd α d 0 N KP d
) and

P 2 = E Iu [ exp( -βd α d 0 Iu KP d )
] .

Proof. For M retransmissions, when the UAV is static, we have to break the D2D coverage probability at each time instance in two components: the first part corresponds to the D2D users contribution and the second component shows the contribution of the UAV. Since the UAV is static, the second component is the same for all time instances but the second part is different due the Rayleigh fading channel. Assuming that the Rayleigh fading channels at different transmission time instances between D2D pairs are uncorrelated,

P i cov,d (β) = P 1,i × P 2 . ( 38
)
Then we have

P O,s cov,d (β|I u ) = P [ γ d,i ≥ β|I u , at least for one of i ∈ {1, ..., M } ] = 1 - ( P [γ d,i < β|I u ] ) M = 1 - ( 1 -P 1,i ) M . ( 39
)
Finally,

P O,s cov,d (β) = P O,s cov,d (β|I u ) × P 2 = P 2 × [ 1 -(1 -P 1,i ) M ] . ( 40
)
From Theorem 3, we can see that, when M → ∞, P O,s cov,d (β) → P 2 which is less than one. However, P O,m cov,d (β) → 1. In other words, in the static UAV case the average coverage probability never tends to one while in the mobile UAV case it can reach one for high values of M . In fact, a very high D2D coverage probability (close to one) for all the users is not achievable in the static UAV case. More specifically, D2D users in the coverage radius of the UAV are more susceptible to a constant high interference from the UAV. By changing the location of the UAV, interference generated by the UAV on the D2D users does not remain high constantly. This is due to the fact that the distance between a D2D user and the UAV changes over time. Thereby, a D2D transmitter which has a higher distance from the UAV, has a higher chance of successful transmission accordingly. Now, using the coverage probability expressions for DUs and D2D users, the average rates for both types of users considering M retransmissions are given by:

Cd (β) = 1 M Rc ∫ 0 2π ∫ 0 M ∑ i=1 C i d (r, φ, β) r πR 2 c drdφ, (41)
Cdu (β) = 1 M Rc ∫ 0 2π ∫ 0 M ∑ i=1 C i du (r, φ, β) r πR 2 c drdφ, (42)
where

C i d (r, φ, β) = P i cov,d (r, φ, β)×W log 2 (1+β) and C i du (r, φ, β) = P i cov,du (r, φ, β)×W log 2 (1+ β) .
Interestingly, increasing M has a different impact on the average rate of DUs and D2D users.

For higher values of M , a downlink user should wait for a longer time until the UAV becomes close to it and provides the required coverage. That is, having higher number of stop points for serving the downlink users results in a higher delay and hence the average rate of DUs decreases.

On the other hand, changing the number of stop points does not considerably change Cd (β). This is due to the fact that D2D users are not served by the UAV and increasing the number of stop points does not cause any delay for D2D communications. However, as will be discussed in the next section, the number of stop points improves the average overall coverage probability and reduces outage area where D2D transmissions are not successful.

V. SIMULATION RESULTS AND ANALYSIS

A. The static UAV scenario

First, we compare our analytical results of the coverage probabilities using numerical simulations. Table II lists parameters used in the simulation and statistical analysis. These parameters are set based on typical values such as in [START_REF] Hourani | Optimal LAP altitude for maximum coverage[END_REF] and [START_REF] Shalmashi | Energy efficiency and sum rate tradeoffs for massive MIMO systems with underlaid device-to-device communications[END_REF]. Here, we will analyze the impact of the various parameters such as the UAV altitude, D2D density, and SINR threshold on the performance evaluation metrics.

In Figures 3 and4, we show, respectively, the D2D coverage probability and approximation of DU coverage probability for different SINR detection threshold values. From these figures, we can clearly see that, the analytical and simulation results for D2D match perfectly and the analytical approximation for DU coverage probability and simulation results are very close.

Figures 3 and4 show that, by increasing the threshold, the coverage probability for D2D users and DUs will decrease. Section III, we can see that the rate depends on the coverage probability, which is a decreasing function of the threshold, β, and an increasing logarithmic function of it. Clearly, for high values of β, the received SINR cannot exceed the threshold and, thus, the coverage probabilities tend to zero. On the other hand, according to ( 27) and (28), as β increases, log 2 (1 + β) increases accordingly. However, since the coverage probability exponentially decreases but log 2 (1 + β) increases logarithmically, the average rate tends to zero for the high values of β. Furthermore, for β → 0, since log 2 (1 + β) tends to zero and the coverage probabilities approach one, the rate becomes zero. Hence, an optimum value for the SINR threshold for which the rate is maximized can exist. As can be seen from Figure 5, for the given parameters in Table I, the maximum rate is achieved for β = 4 and 8 for λ d = 10 -4 and 0.5 × 10 -4 , respectively.

Figure 6 shows the impact of D2D density on the sum-rate. In this figure, we can see that a low D2D density yields low interference. However, naturally, decreasing the number of D2D users in an area will also decrease the sum-rate. For high D2D density, high interference reduces the coverage probability and consequently the data rate for each user. However, since the sumrate is directly proportional to the number of D2D users, increasing the D2D density can also improve the sum-rate. According to the Figure 6, as the density of downlink users increases, the optimal λ d that maximizes the sum-rate decreases. This is due to the fact that, as λ du increases, the contribution of DUs in the average sum-rate increases and hence increasing the rate of each DU enhances the average sum-rate. To increase the rate of a DU, the number of D2D users as the interference source for DUs should be reduced. As a result, the optimal λ d decreases as as λ du increases. For instance as shown in the figure, by increasing λ du from 10 -4 to 4 × 10 -4 , the optimal λ d decreases from 0.9 × 10 -4 to 0.3 × 10 -4 .

It is important to note that the value of the fixed distance, d 0 , between the D2D pair significantly impacts the rate performance. Figure 7 shows the Csum as a function of the density of D2D users and d 0 . From this figure, we can see that, the rate increases as the fixed distance between a D2D receiver and its corresponding transmitter decreases. Moreover, the optimal D2D density which leads to a maximum Csum , increases by decreasing d 0 . In fact, for lower values of D we can have more D2D users in the network. For instance, by reducing d 0 from 8 m to 5 m, the optimum average number of D2D users increases by a factor of 3.

Figure 8 shows the coverage probability for DUs and D2D users as a function of the UAV altitude. From the DUs' perspective, the UAV should be at an optimal altitude such that it can provide a maximum coverage. In fact, the UAV should not position itself at very low altitudes, due to high shadowing and a low probability of LOS connections towards the DUs. On the other hand, at very high altitudes, LOS links exist with a high probability but the large distance between UAV and DUs results in a high the path loss. As shown in Figure 8, for h = 500 m the DU coverage probability is maximized. Note that from a D2D user perspective, the UAV creates interference on the D2D receiver. Therefore, D2D users prefer the UAV to be at an altitude for which it provides a minimum coverage radius. As seen in Figure 8, for h → ∞, the D2D users and 30 m. Note that the optimal h that maximizes the sum-rate depends on the density of DU and D2D users. From Figure 9, considering d 0 = 20 m as an example, we can see that for h > 1300 m, the average sum-rate starts increasing. This stems from the fact that the DU coverage probability tends to zero and, thus, only D2D users impact Csum . Hence, as the UAV moves up in altitude, the interference on D2D users decreases and Cd increases. Moreover, for 300 m < h < 1300 m, Figure 9 shows that the coverage probability and, consequently, the average rate for the downlink users decrease as the altitude increases. However, increasing the UAV altitude reduces the interference on the D2D users and improves the average rate for D2D users. In addition, in this range of h, since DUs have more contributions on Csum than the D2D users, Csum is a decreasing function of altitude.

B. The mobile UAV scenario

Next, we study the mobile UAV scenario. In this case, we can satisfy the coverage requirement for all the DUs. In fact, the UAV moves over the target area and attempts to serve the DUs at the stop points to guarantee that all the DUs will be in its coverage radius. Figure 10 shows the coverage radius of the mobile UAV when it is located at the optimal altitude as the D2D density varies. As expected, the coverage radius decreases as the D2D density increases. For instance, for ε = 0.6, when λ d increases from 10 -5 to 10 -4 , the coverage radius decreases from 1600 m to 300 m. Moreover, by reducing the minimum coverage requirement of DUs, the UAV can cover a larger area. For instance, reducing ϵ from 0.6 to 0.4 increases the UAV coverage radius from 290 m to 380 m for λ d = 10 -4 . Note that, since the main goal of the UAV is to provide coverage for the entire target area, to compensate for the low coverage radius, we should increase the number of stop points for serving the DUs and consequently a longer time is required for the full coverage.

In Figure 11, we show the minimum number of stop points as a function of the D2D user density. In this figure, we can see that, as expected, the number of stop points must increase when the density of D2D users increases. In fact, to overcome the higher interference caused by increasing the number of D2D users, the UAV will need more stop points to satisfy the DUs' coverage constraints. For instance, when λ d increases from 0.2 × 10 -4 to 0.8 × 10 -4 , the number of stop points must be increased from 3 to 8. Note that, when computing the minimum number of stop points for each λ d , we considered optimal values for the UAV altitude such that it can provide a maximum coverage for the DUs. Therefore, the UAVs altitude changes according to the D2D density. Moreover, as seen from Figure 11, the minimum number of stop points remains constant for a range of λ d . This is due to the fact that the number of stop points is an integer and hence, for different values of λ d , the integer value will be the same. However, although the minimum number of stop points for two different D2D densities are the same, the UAV can transmit with lower power in the case of lower D2D density.

In Figure 12, we show the minimum number of stop points as a function of the UAV altitude for λ d = 10 -4 . Figure 12 shows that, for some values of h which correspond to the optimal UAV altitude, the minimum number of stop points is minimized. For example, the range of optimal h for ϵ = 0.4 and ϵ = 0.6 is, respectively, 400 m < h < 500 m and 300 m < h < 350 m. As expected, the minimum number of stop points is lower for the lower value of ϵ.

Next, we compare the D2D coverage performance in the static and mobile UAV scenarios.

For a fair comparison, we consider the same number of retransmissions for both scenarios. In other words, the number of stop points is equivalent to the number of retransmissions.

Figure 13 shows the tradeoff between the downlink coverage probability and the delay which is considered to be proportional to the number of stop points. In Figure 13, we can see that, in order to guarantee a higher coverage probability for DUs, the UAV should stop at more locations.

As observed in this Figure , for λ d = 10 -4 , to increase the DU coverage probability from 0.4 to 0.7, the number of stop points should increase from 5 to 23. For a higher number of stop points, the UAV is closer to the DUs and, thus, it has a higher chance of LOS. However, on the average, a DU should wait for a longer time to be covered by the UAV that reaches its vicinity. In addition, as the density of D2D users increases, the number of stop points (delay) increases especially when a high coverage probability for DUs must be satisfied. For instance, if λ d increases from 0.5 × 10 -4 to 10 -4 , or equivalently from 50 to 100 for the given area, the number of stop points should increase from 4 to 9 to satisfy a 0.5 DU coverage probability, and from 20 to 55 for a 0.8 coverage requirement.

Figure 14 shows the overall coverage probability for a D2D user located at the center of the target area. As the number of retransmissions (stop points) increases, the overall coverage probability also increases for both static and mobile UAV cases. However, the coverage probability enhancement in the mobile UAV case is significantly higher than the static case. For example, for 5 retransmissions, as compared to the static UAV, we observe a 21% improvement in the overall D2D coverage probability by moving the UAV. Note that, a D2D user, prefers to be outside the coverage range of the UAV to experience a low interference from it. For the static UAV case, the coverage probability for a D2D user located within the coverage range of the UAV is low due to the high interference stemmed from the UAV. On the other hand, if the UAV moves, the interference on the D2D user decreases in the next time instances.

In Figure 15, we present the overall D2D coverage probability for the static and mobile UAV cases. We consider four stop points for the mobile UAV case and four retransmissions for the static UAV case. Figure 15 shows that, the variation of coverage probability at different locations for the static case is significantly higher than the mobile UAV case. The minimum coverage probability is 0.002 and 0.48 in the static and mobile UAV cases, respectively. From Figure 15, we can see that, the mean and standard deviation of coverage probability are 0.51 and 0.27 for the static case, and 0.59 and 0.06 for the mobile UAV case. More importantly, Figure 15a shows that, in the static case, the coverage probability at 41% of the locations is below 0.5 whereas, as we can see in Figure 15b, this value for the mobile UAV case is 16%. Hence, as compared to the static case, the mobile UAV provides a higher average overall coverage probability for the D2D users and more fairness in terms of coverage for the D2D users in different locations.

VI. CONCLUSIONS

In this paper, we have studied the performance of a UAV that acts as a flying base station in an area in which users are engaged in D2D communication. We have considered two types of users: in the network: the downlink users served by the UAV and D2D users that communicate directly with one another. For both types, we have derived tractable expressions for the coverage probabilities as the main performance evaluation metrics. The results have shown that a maximum average sum-rate can be achieved if the UAV altitude is appropriately adjusted based on the D2D users density. Furthermore, as compared to the static UAV case, moving the UAV enhances the overall coverage performance of both DUs and D2D users. In the mobile UAV scenario, using the disk covering problem, the entire target area (cell) can be completely covered by the UAV in a shortest time with a minimum required transmit power. Finally, we have analyzed the tradeoff

Figure 1 :

 1 Figure 1: Network model including a UAV, downlink users and D2D.

Figure 2 :

 2 Figure 2: Five disks covering problem.

Figure 5

 5 Figure 5 illustrates the average sum-rate (Gbps) versus the threshold for 1 MHz transmission bandwidth, λ du = 10 -4 , h = 500 m, and two different values of λ d . By inspecting (30) in

Figure 3 :

 3 Figure 3: D2D coverage probability vs. SINR threshold

Figure 4 :

 4 Figure 4: DU coverage probability vs. SINR threshold.

Figure 5 :

 5 Figure 5: Average sum-rate vs. SINR threshold.

4 λ

 4 du =0.5*10 -4

Figure 6 :Figure 7 :

 67 Figure 6: Average sum-rate vs. D2D density (number of D2D pairs per m 2 ).

Figure 9 shows

 9 Figure 9 shows Csum versus the UAV altitude for different values of the fixed distance, d 0 , the fixed distance between a D2D transmitter/receiver pair. The optimum values for the height which lead to a maximum Csum are around 300 m, 350 m, and 400 m for d 0 = 20 m, 25 m

Figure 8 :Figure 9 :

 89 Figure 8: Coverage probability vs. UAV altitude.

Figure 10 :

 10 Figure 10: Maximum UAV coverage radius vs. D2D density (number of D2D pairs per m 2 ).

Figure 11 :Figure 12 :

 1112 Figure 11: Number of stop points vs. D2D density.

4 Figure 13 :

 413 Figure 13: Minimum number of stop points vs. coverage probability (coverage-delay tradeoff)

Figure 14 :

 14 Figure 14: Overall D2D coverage probability vs. number of retransmissions.

Figure 15 :

 15 Figure 15: Overall D2D coverage probability vs. location of a D2D user.

Table I :

 I Number and radii of disks in the covering problem.

Number of stop points Minimum required coverage radius (R min

Table II :

 II Simulation parameters.

	Description	Parameter	Value
	UAV transmit power	P u	5 W
	D2D transmit power	P d	100 mW
	Path loss coefficient	K	-30 dB
	Path loss exponent for UAV-user link	α d	2
	Path loss exponent for D2D link	α u	3
	Noise power	N	-120 dBm
	Bandwidth	W	1 MHz
	D2D pair fixed distance	d 0	20 m
	Excessive attenuation factor for NLOS	η	20 dB
	Parameters for dense urban environment	B, C	0.136, 11.95

Note that the number of users has a Poisson distribution but their location follows the uniform distribution over the area.