Mounia Hamidouche

Ejder Baştug

Jihong Park
email: jihong@es.aau.dk

Laura Cottatellucci
email: laura.cottatellucci@eurecom.fr

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

Downlink Performance of Dense Antenna Deployment: To Distribute or Concentrate?

Keywords: massive MIMO, small cells, energy efficiency, coverage probability, downlink, stochastic geometry, 5G

Massive multiple-input multiple-output (massive MIMO) and small cell densification are complementary key 5G enablers. Given a fixed number of the entire base-station antennas per unit area, this paper fairly compares (i) to deploy few base stations (BSs) and concentrate many antennas on each of them, i.e. massive MIMO, and (ii) to deploy more BSs equipped with few antennas, i.e. small cell densification. We observe that small cell densification always outperforms for both signal-tointerference ratio (SIR) coverage and energy efficiency (EE), when each BS serves multiple users via L number of subbands (multi-carrier transmission). Moreover, we also observe that larger L increases SIR coverage while decreasing EE, thus urging the necessity of optimal 5G network design. These two observations are based on our novel closed-form SIR coverage probability derivation using stochastic geometry, also validated via numerical simulations.

I. INTRODUCTION

Performance improvements such as having wider coverage, higher user data rate, higher energy efficiency and lower latency are under an intensive investigation for the design of the fifth generation (5G) and beyond cellular architectures. These improvements are fundamental to achieve the dramatic growth of connected devices and the tremendous amount of data in applications such as voice, videos, and games [START_REF]Cisco visual networking index: Global mobile data traffic forecast update 2016-2021[END_REF], as well as applications in wireless virtual-reality [START_REF] Baştug | Towards interconnected virtual reality: Opportunities, challenges and enablers[END_REF]. Novel innovative network technologies are used to meet the required performances. First, transmission with massive multiple-input multiple-output (MIMO) [START_REF] Hoydis | Massive MIMO in the UL/DL of cellular networks: How many antennas do we need?[END_REF] is considered as a candidate technology for 5G. The key feature of this technology is the use of a large number of antennas at the base station compared to the number of users. The more antennas the base stations are equipped with, the better the performance is in terms of data rate and energy consumption [START_REF] Wong | Performance enhancement of multiuser MIMO wireless communication systems[END_REF]. A second promising technology is small cell networks [START_REF] Hoydis | Green small-cell networks[END_REF], that consists of a dense number of small cell base stations in a given area. Due to the short distance between the base station (BS) and the user terminals, small cell networks have a low This research has been supported by the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for ComplexNetwork Engineering), and the U.S. National Science Foundation under Grant CCF-140922.

path loss thus yielding a low power consumption which can improve the energy efficiency (EE). For a given total number of BS antennas, interesting strategies in the deployment of these two technologies are either i) low density deployment of base stations with many antennas, i.e, massive MIMO or ii) a higher density deployment of BSs equipped with fewer antennas. Understanding which strategy is preferable is one of the goals of this work.

In fact, many works analyzed the impact of the massive amount of antennas on the EE. In particular, the work in [START_REF] Björnson | Optimal design of energy-efficient multi-user MIMO systems: Is massive MIMO the answer?[END_REF] solves the EE maximization problem for a multi-cell multiuser MIMO network and shows that small cells yield higher EE. In [START_REF] Björnson | Optimal design of energy-efficient multi-user MIMO systems: Is massive MIMO the answer?[END_REF], the authors give insights on how the number of antennas at the BS must be chosen in order to uniformly cover a given area and attain maximal EE. Altough many works study massive MIMO and small-cell densification, very few have focused on comparing their performance. Our goal is to analyze which one of the two technologies perform better in terms of the coverage probability and energy efficiency.

A comparison has been recently presented in [START_REF] Nguyen | Massive MIMO versus small-cell systems: spectral and energy efficiency comparison[END_REF]. The massive MIMO and small-cell systems were compared in terms of spectral and energy efficiency bounds. The authors observe via simulations that for the average spectral efficiency, small-cell densification is favourable in crowded areas with moderate to high user density and massive MIMO is preferable in scenarios with low user density. In contrast to the analysis in [START_REF] Nguyen | Massive MIMO versus small-cell systems: spectral and energy efficiency comparison[END_REF], we derive exact expressions of the coverage probability and EE by assuming other constraints on the model then we compare between massive MIMO and small cell networks in terms of these two metrics. One of the intersting constraint is to assume multi-carrier transmission in which the total bandwidth is divided into L ≥ 1 sub-bands. Then, instead of studying the downlink performance when each BS serve a single user at each time/frequency, we consider that each BS is scheduled to serve simultaneously multiple users on each sub-band. We also cancel the interference by using zero-forcing (ZF) processing. In addition, instead of introducing massive MIMO and smallcell systems separately, we examine the problem with a single system model by varying the number of BS antennas under the constraint of a fixed total number of BS antennas per unit area. In this work, we derive analytic expressions for the coverage probability and EE using a stochastic geometry approach [START_REF] Baccelli | Stochastic geometry and wireless networks: Volume ii applications[END_REF]. The key feature of this approach is that the base station positions are all independent which allows to use tools from stochastic geometry.

The rest of this paper is organized as follows. Section II details our system model of downlink transmission using linear processing ZF under perfect channel state information (CSI) at each base station. General expressions for coverage probability and EE are derived in Section III. In Section IV, numerical results are used to validate the theoretical analysis and make comparisons between massive MIMO and small cell densification for both coverage probability and EE metrics. Finally, the major conclusions and implications are given in Section V.

The following notation is used in this paper. The expectation operation with respect to a random variable and the absolute value are denoted by E{.} and |.|, respectively. We denote by I M the M ×M identity matrix, and we use CN (0, Σ) to denote a circularly symmetric complex Gaussian distribution with zero-mean and covariance matrix Σ. The Gamma function is denoted as Γ(.). The bold lower-case letters as h represent vectors, whereas the bold upper-case as H are matrices.

II. NETWORK MODEL

The cellular network consists of BSs independently distributed according to a homogeneous Poisson point process (PPP) Φ of intensity λ BS (measured in BSs/km 2 ), and is depicted in Figure 1. Each base station is equipped with an array of M antennas. We consider an independent collection of single antenna mobile users, located according to another independent stationary PPP Ψ with intensity λ UE . We assume that each user connects to its closest BS, namely each BS serves the users which are located within its Voronoi cell [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF]. In this section, assuming perfect CSI at each BS we study the signal model for downlink system. With this goal in mind, we first consider a typical user, which is connected to a tagged BS (BS 0 ). Since user locations are translation-invariant, we consider that the typical user is always located at the origin. This typical user's received signal y is then given as

y = r -α 0 h H 0 x 0 Desired signal + BSi∈Φ\{BS0} r -α i h H i x i Interference +n, (1)
where the stochastic vector h i ∈ C M denotes the small scale fading between the i-th base station to the typical user. It follows a complex generalized Gaussian distribution denoted as h i CN (0, I M ). The channel is considered to be noisy, with the Gaussian noise n of variance K P added to the received signal. The variable r i is the distance from the typical user to its closest base station BS i and α ≥ 2 is the path-loss exponent. We denote x i ∈ C M an arbitrary symbol transmitted from the i-th base station. In addition to the received signal model in [START_REF]Cisco visual networking index: Global mobile data traffic forecast update 2016-2021[END_REF], we define the signal-to-interference-plus-noise ratio (SINR) of the typical user as

SINR = r -α 0 h H 0 x 0 BSi∈Φ\{BS0} r -α i h H i x i + n . (2)
We suppose that the BSs must be deployed to match a given finite user density of λ UE UEs/km 2 , then each base station serves in average K = λUE λBS users. The total bandwidth W is divided into L ≥ 1 sub-bands. Therefore, K = K L ≤ M users are simultaneously served on each sub-band by each base station. We assume that the total number of antennas λ BS M is fixed and should be deployed in a given area. Based on this assumption, for simplicity we set λ BS M = λ UE , then the number of antennas in each base station is given by M = λUE λBS which is always equal or greater than K. We suppose that P = E[x H i x i ] is the average transmit power per base station which is given as Pmax λBS , where P max is the maximum power used when all antennas are concentrated on a single base station. In this scenario, to cancel out the interference while boosting the desired signal power, each BS applies ZF transmission to simultaneously serve K single antennas. Let s i,k CN (0, 1) be the message (the symbol) determined for the user k from the i-th base station. Then, the i-th BS multiplies the data symbol s i,k destined for the k-th user by w i,k . Therefore, the linear combination x i of the symbols transmitted by the i-th base station intended for the K users is

x i = K k=1 w i,k s i,k , (3)
where w i,k ∈ C M ×1 is ZF beamforming vector. Then, the received SINR at the typical user can now be expressed as

SINR = r -α 0 |h H 0 w 0,1 | 2 BSi∈Φ\{BS0} r -α i K k=1 |h i w i,k | 2 + K P = r -α 0 S I r + K P , (4)
where

I r = BSi∈Φ\{BS0} r -α i g i with g i = K k=1 |h i w i,k | 2
denotes the interference channel power and S = |h H 0 w 0,1 | 2 is the desired channel power. Now, we introduce the metrics we will investigate in the next sections.

Definition 1 (Coverage Probability). The coverage probability of a typical user is the the probability that the SINR received by the user is larger than a predefined threshold T such as P cov (T ) = P(SINR > T ).

(

) 5
Definition 2 (Energy Efficiency). The EE is defined as

EE = ASE AEC = Area Spectral Efficiency [bit/symbol/km 2 ] Average Energy Consumption [Joule/symbol/km 2 ] , (6)
where the area spectral efficiency (ASE) is expressed as in which λ BS represents BS density, K is the number of users that are served by the BSs and E[R] is the average data rate of users. Moreover, average energy consumption (AEC) is defined as similar to [START_REF] Zhang | Energy efficient downlink cooperative transmission with BS and antenna switching off[END_REF], that is

ASE = λ BS KE[R], (7)
AEC = P η + M P c + K 3 P pre + P 0 , (8)
where η denotes the power amplifier efficiency, P c is the circuit power per antenna, which indicates the energy consumption of the corresponding RF chains. The term K 3 P pre accounts for the energy consumption for precoding which is related to the number of users served simultaneously by each base station.

The term P 0 is the non-transmission power, which accounts for the energy consumption of baseband processing.

Remark 1. The goal in EE is to maximize the performance while minimizing the energy consumption which are two conflicting operations in 5G.

III. PERFORMANCE ANALYSIS

In this section we give the analytical expression of the coverage probability on a typical mobile user. Afterwards, we shall give the expression of the EE. We start by stating the following Lemma that shall be used to derive our results.

Lemma 1 ([9]

). The probability density function (PDF) of a typical user's association distance r 0 is

f r0 (r) = dF r0 (r) dr = e -λBSπr 2 2πλ BS r. (9)
Proof. Each user is connected to its closest base station, then all the interfering base stations are farther than a distance r.

Since the Poisson distribution helps in describing the chances of occurrence of a number of events in a given space, then the probability that no base station is closer than a distance r within an area πr 2 , is e -λBSπr 2 , expressed as

P(r 0 > r) = F r0 (r) = P[No BS within πr 2 ] = e -λBSπr 2 (10)
Therefore, the PDF results from the derivative of the cumulative distribution function F r (R).

A. Coverage Probability

Before deriving the expression of the coverage probability, we first derive the Laplace transform of both interference and desired signal. The desired channel power S is distributed as Γ(M -K + 1, 1) [START_REF] Dhillon | Downlink MIMO Hetnets: Modeling, ordering results and performance analysis[END_REF]. For the interfering signal, as w i,k is a unit-norm vector and independent of h i , then |h i w i,k | 2 is a squared-norm complex Gaussian, which is exponential distributed. For tractability we neglect the correlation between w i,k for different k, then the channel gain g i is the sum of K independent exponential distributed random variables which follows Γ(K, 1). Lemma 2 (Laplace Transform of Interference). The Laplace transform of interference L Ir (s)=E[e -sIr ] is

LI r (s) = exp -πλBSr 2 0 -1 + 2F1 K, -2 α ; -2 α + 1; -sr -α 0 , ( 11
)
where 2F1 is the Gauss-Hypergeometric function.

Proof. See Appendix A.

Lemma 3 (Laplace Transform of the Desired Signal). The Laplace transform of the desired signal

L S (s) = E[e -sS ] is L S (s) = 1 1 + s M -K+1 . ( 12
)
Proof. See Appendix B.

Combining the previous results given in Lemmas 1, 2 and 3 with the proof techniques proposed in [START_REF] Baccelli | Stochastic analysis of spatial and opportunistic ALOHA[END_REF], an expression for the coverage probability can be derived and it is given in the following theorem.

Theorem 1 (Coverage Probability in Downlink). The coverage probability at a typical mobile user in the general cellular network model described above is

BS density λBS 4

Path-loss exponent α 4

P cov (T ) = r0>0 ∞ -∞ L Ir 0 (i2πr α 0 T s) exp - i2πr α 0 T K P s × L S (-i2πs) -1 i2πs f r0 (r 0 )dsdr 0 , (13)
where f r0 (r 0 ) is the PDF of the distance between the typical user and the tagged base station (BS 0 ), L Ir 0 (.) is the Laplace transform of the interference and L S (.) is the Laplace transform of the desired signal (Lemmas 1, 2 and 3 respectively).

Proof. See Appendix C.

B. Energy Efficiency

To facilitate the analysis of the EE, we consider fixed modulation and coding schemes for each user by considering a fixed SINR threshold T as in [START_REF] Park | Coverage and rate of downlink sequence transmissions with reliability guarantees[END_REF], providing the average rate E[R] as a function of downlink coverage probability as

E[R] = log(1 + T )P cov (T ) = log(1 + T ) r0>0 ∞ -∞ L Ir 0 (i2πr α 0 T s)× exp - i2πr α 0 T K P s L S (-i2πs) -1 i2πs f r0 (r 0 )dsdr 0 , (15)
By plugging the average rate expression into [START_REF] Nguyen | Massive MIMO versus small-cell systems: spectral and energy efficiency comparison[END_REF], we obtain the ASE. Then, the expressions of EE can be readily obtained and its final expression is given in ( 14) on the top of this page.

IV. NUMERICAL RESULTS

In this section, we conduct Monte-Carlo simulations to validate the analytical expressions of coverage probability and EE of our multi-user MIMO system. The default parameter setting is given in Table I and shall be used unless otherwise stated.

Figure 2a illustrates the coverage probability expression provided in Theorem 1 as a function of target SINR for three different path-loss values: α ∈ {3, 4, 5}. The curves reveal that the coverage probability obtained by simulation behaves exactly as the analytical results which confirm the accuracy of our theoretical expressions. Also, increasing α increases the coverage probability because the interference power decreases faster as a function of α than the power signal [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF].

Figure 2b shows the coverage probability as a function of BS density. Several important observations can be made from the results of this figure. First, for any given L, the coverage probability can be greatly improved by increasing the BS density, meaning that distributed network densification is preferable over massive MIMO. Second, we note that the coverage probability increases as L increases, thus showing the importance of multi-carrier transmissions. However, it is shown that for large L, both massive MIMO and small cells provide the same coverage which is confirmed according to Figure 2c. This phenomenon occurs since the number of users K served by each BS decreases and become very small compared to the number of its own antennas M .

Figure 2d shows the impact of varying λ BS on the EE. The curves show that for any given L, the EE can also be improved by increasing the BS density, meaning that small cells is also preferable over massive MIMO. In contrast with the coverage probability, we can observe that in small cells scenario, decreasing the number of sub-bands L increases the EE but in massive MIMO scenario, large number of sub-bands provides the highest EE. Finally, we notice also that both massive MIMO and small cells provide the same performance gain when L becomes very large which means that the EE increase when the number of users served by the base station approach the number of its base station antennas M . All these observations show that the coverage and EE are conflicting such that improvements in one objective lead to degradation in the other objective for a fixed number of sub-bands L.

V. CONCLUSIONS Our proposed model was based on stochastic geometry, where the BS and user locations were distributed according to PPP. Using tools from stochastic geometry, we derived the coverage probability and EE expressions for downlink scenario. The coverage probability expression was validated via Monte-Carlo simulations. The comparison showed that for any given sub-band L, small cell densification is preferable over massive MIMO if both coverage probability and EE should be increased. However, increasing L improves the coverage probability but decreases the EE which shows that the two metrics are conflicting. An interesting future work is therefore to introduce a multi-objective optimization framework and find (possibly) optimal number of sub-bands and BS density that maximize the coverage probability and EE jointly for the uplink and downlink scenarios. Future work will also look at multiple antenna terminals which is still an open topic. interference is L Ir (s) = E[e -sIr ], where the average is taken over both the spatial PPP and the interference distribution is expressed as follows:

EE = λBSK r 0 >0 ∞ -∞ LI r 0 (i2πr α 0 T s) exp - i2πr α 0 T K P s × LS(-i2πs) -1 i2πs fr 0 (r0)dsdr0 log(1 + T ) P η + M Pc + K 3 Ppre + P0 (14)
L Ir (s) = E Φ,{gi}   exp   -s BSi∈Φ\{BS0} r -α i g i     , = E Φ,{gi}   BSi∈Φ\{BS0} exp(-sg i r -α i )   , (a) = E Φ   BSi∈Φ\{BS0} E gi exp(-sg i r -α i )   , (b) = E Φ   BSi∈Φ\{BS0} 1 Γ(K) ∞ 0 e -gi(sr -α i +1) g K-1 i dg i   , (d) = E Φ   BSi∈Φ\{BS0} 1 (1 + sr -α i ) K   , = exp -2πλ BS ∞ r0 1 - 1 (1 + sv -α ) K vdv (16)
where the step (a) follows from the i.i.d distribution of g i and further independence from the point process Φ. The step (b) follows from the PDF of g i Γ(K, 1) given as

f (g) = 1 Γ(K) g K-1 e -gi . (17)
Moreover, the step (c) follows from the computation of the integral by the means of integration by parts. The last step follows from the probability generating functional of the PPP with intensity λ [START_REF] Haenggi | Stochastic geometry for wireless networks[END_REF], which states that for some function f (x)

we have

E   BSi∈Φ\{BS0} f (x)   = exp -λ R 2 (1 -f (x)dx) .
(18) The inside integral can be evaluated by using the change of variables v -α → y and we obtain the result.

APPENDIX B PROOF OF LEMMA 3

The Laplace transform of the desired signal is L S (s) = E[e -sS ], where the average is taken over the desired signal distribution expressed as follows:

L S (s) = E S [e -sS ] = 1 Γ(M -K + 1) ∞ 0 S M -K e -S(1+s) dS, a = 1 1 + s M -K+1 , (19)
where the first step follows from the PDF of the desired signal

f (S) = 1 Γ(M -K + 1) S M -K e -S dS. (20)
The step (a) follows from the computation of the integral by the means of integration by parts.

APPENDIX C PROOF OF THEOREM 1

The first part of the proof follows by conditioning on the nearest BS being at a distance r 0 from the typical user. Then, the probability of coverage is P cov (T, λ BS , α) = r0>0 e -πλBSr 2 0 P( r -α 0 S I r + K P > T ) 2πλ BS r 0 dr 0 .

(21) To evaluate P( r -α 0 S Ir+ K P > T ), we use the proof techniques proposed in [START_REF] Andrews | A tractable approach to coverage and rate in cellular networks[END_REF] and [START_REF] Baccelli | Stochastic analysis of spatial and opportunistic ALOHA[END_REF]. Some assumptions are required for this computation: A1) The desired signal S admits a square integrable density. A2) Either the interference I r or the noise admits a density which is square integrable. The interference and the noise are independent, then the second assumption imply that I r + K P admits a PDF f Ir+ K P (y) that is square integrable. Therefore, the coverage probability is expressed as follows: Using Fubini's theorem [START_REF] León | Fubini theorem for anticipating stochastic integrals in Hilbert space[END_REF], and moving the expectation inside Combining the last expression and (21) gives the result stated in Theorem 1.

P (SINR > T ) = P

Figure 1 :

 1 Figure 1: An illustration of the network model.

Figure 2 :

 2 Figure 2: Evolution of coverage probability with respect to the a) the target SINR. b) the BS density in BS/km 2 , and c) number of sub-bands L. The subfigure d) shows the evolution of EE with respect to the BS density in BS/km 2 .

  1 [0≤y≤S(T r -α 0 ) -1 ] dy . (22)where the step (a) above follows from the definitionP(a ≤ X ≤ b) = b a f (x)dx, (23)where f is the density function of the variable X. Using the Plancheral-Parseval theorem[START_REF] Brémaud | Mathematical principles of signal processing: Fourier and wavelet analysis[END_REF], we obtain

  LI r (-2iπT r α 0 s) exp -2iπT r α 0 K P s LS(e 2iπys ) -1 2iπs ds.

Table I :

 I DEFAULT SIMULATION PARAMETERS.

	System Parameter	Symbol	Value
	Power amplifier	η	0.318
	Circuit power per antenna	Pc	14.8 W
	Energy consumption for precoding	Ppre	1.74 W
	Nr. of Sub-bands	L	1
	Target SINR	T	1 dBm
	Non-transmission power	P0	65.8 W
	Maximum average power per BS	Pmax	40 dBm
	Users density	λUE	32 per km 2


Kenza Hamidouche
email: kenza.hamidouche@centralesupelec.fr

Saad Walid
email: walids@vt.edu

Merouane Debbah
email: merouane.debbah@huawei.com

Walid Saad

Mérouane Debbah

Multi-Games for LTE and WiFi Coexistence over Unlicensed Channels

published or not. The documents may come

Introduction

One promising approach to overcome the scarcity of the radio spectrum is to enable tomorrow's small cell networks (SCNs) to exploit simultaneously their licensed cellular bands along with the unlicensed, WiFi bands [START_REF] Liu | Small cell traffic balancing over licensed and unlicensed bands[END_REF][START_REF] Zhang | LTE-unlicensed: The future of spectrum aggregation for cellular networks[END_REF][START_REF] Zhang | Coexistence of wi-fi and heterogeneous small cell networks sharing unlicensed spectrum[END_REF][START_REF] Li | Modeling and analyzing the coexistence of wi-fi and lte in unlicensed spectrum[END_REF][START_REF] Hamidouche | A multi-game framework for harmonized lte-u and wifi coexistence over unlicensed bands[END_REF].

Although offloading part of the LTE traffic to the unlicensed bands can considerably increase the performance of cellular networks, the disparity of the medium access protocols that are used by the WiFi users (WUs) and dual mode SBSs, rise new challenges. In fact, the access to the spectrum in WiFi is based on carrier sense multiple access with collision avoidance (CSMA/CA) while LTE uses orthogonal frequency division multiple access (OFDMA). In many countries, there are no imposed regulations and the operators need to develop new scheduling mechanisms known as LTE-U to extend the LTE into the unlicensed channels.

In this regard, many works have proposed and analyzed new mechanisms for the deployment of LTE-U over unlicensed channels. The work in [START_REF] Chen | Echo state networks for self-organizing resource allocation in LTE-U with uplinkdownlink decoupling[END_REF] formulated the unlicensed spectrum allocation problem with uplink-downlink decoupling as a noncooperative game in which the SBSs are the players that select the unlicensed channels over which they serve their users. The authors in [START_REF] Gu | Exploiting student-project allocation matching for spectrum sharing in LTE-unlicensed[END_REF] formulated the unlicensed spectrum allocation problem as a student-project matching problem with externalities.The work in [START_REF] Kang | Mobile data offloading through a third-party wifi access point: An operator's perspective[END_REF] proposed a distributed traffic offloading scheme for LTE-U scenarios with a single base station. In [START_REF] Chen | Embedding lte-u within wi-fi bands for spectrum efficiency improvement[END_REF], the authors proposed hyper-access points (HAP) in which the functions of both an LTE SBS and a Wi-Fi access point (AP) are combined in the same HAP. In [START_REF] Galinina | 5G multi-RAT LTE-WiFi ultra-dense small cells: performance dynamics, architecture, and trends[END_REF], the authors analyze the performance dynamics of integrated LTE-U and WiFi when sharing the unlicensed bands.

Despite being interesting, in all of these existing works, the WUs are either assumed to transmit on fixed channels or the allocation of the unlicensed channels to the WUs is determined without accounting for the existing SBSs.

The main contribution of this paper is to propose a new unlicensed spectrum allocation approach that accounts for the impact of the SBSs' transmit power and the assigned channels to each of the WUs, on one another. We formulate a multigame [START_REF] Hamidouche | A multi-game framework for harmonized lte-u and wifi coexistence over unlicensed bands[END_REF] with two different types of games in which the WUs are considered as leaders and the dual-mode SBSs are the followers. On the followers side, the power allocation problem at the SBSs on the unlicensed channels is formulated as a noncooperative. On the other hand, we formulate the channels allocation problem as a one-sided matching game with externalities [START_REF] Abdulkadiroglu | House allocation with existing tenants[END_REF], in which the WUs predict the transmit power of the SBSs and autonomously select the channel over which they serve their traffic. For the formulated multi-game, we prove that there exists a unique Debreu equilibrium (DE) for the power allocation problem. Moreover, we propose a new distributed matching algorithm for the assignment of the WUs to the unlicensed channels and prove that the proposed algorithm converges to a unique stable matching that is in the core.

System Model

Consider a wireless network composed of a set N of N users served by a set S of S dual-mode small base stations (SBSs) and a set W of W WiFi users (WUs). The SBSs and WUs can transmit their content over a set C of C unlicensed channels. Here, each of the channels can be allocated to multiple WUs at a given time but only one WU can transmit based on the LBT scheduling scheme. These same bands when available, are allocated to the SBSs allowing them to aggregate LTE and LTE-Unlicensed (LTE-U) to improve the quality-of-service (QoS) of the served users.

The SBSs only serve the downlink traffic over unlicensed channels while in the WiFi network, both downlink and uplink transmissions can occur. The set of W WUs is denoted W. The achievable throughput by a WiFi user w in the network depends on the number of active WUs and SBSs that transmit over the same unlicensed channel c ∈ C, which we denote R w (W c , S c ), where W c is the number of active WUs on channel c and S c is the number of SBSs that transmit over channel c with a positive transmit power. The faction of time each WU w uses a channel c, also known as the channel busy time, is given by [START_REF] Lingzhi | Online channel selection and user association in high-density wifi networks[END_REF]:

t wc (W c , S c ) = l wc R w (W c , S c ) + γ wc , (1)
where l w is the size of transmitted file by the WU w and γ wc is the channel access overhead and the protocol overhead which depends on the type of the file. The WiFi network throughput R(W c , S c ) is given by:

R w (W c , S c ) = P s w L P s w T s + P c T c + P σ T σ , (2)
where

P s w = τ w ∏ W c +S c i =w (1 -τ i )
is the probability of a successful transmission which corresponds to the probability of having only one transmission over the channel and τ w is the probability of transmission by a given WU w. L is the average size of a packet file, T s is the average time needed to transmit a packet of size L,

P σ = ∏ W c i=1 (1 -τ i )
is the probability of the channel being idle, T σ is the duration of the idle period, P c is the probability of collision, and T c is the average time spent in the collision. The WUs detct collisions based on the power or interference they sense on the unlicensed channels. In fact, if the sensed interference from all other WUs and SBSs over a channel c exceeds a given threshold I th , the WUs consider the channel as busy and back-off. For the SBSs, the achievable physical throughput by a user j that is served by SBS i is given by:

R ic = ω c log 1 + p ic |h ic | 2 σ 2 + ∑ l∈S\i p lc |h lc | 2 + ∑ j∈W α jc p jc |h jc | 2 , (3)
where p ic is the transmission power of SBS/WU i through channel c, |h ic | is the channel gain of SBS/WU i over channel c, and ω c is the bandwidth of band c. α jc is a boolean that is related to the CSMA/CA transmission mode indicating that WU j is transmitting on channel c when α jc = 1, and α jc = 0 otherwise.

Game Formulation and Analysis

To ensure a harmonious coexistence between the SBSs and the WUs, we formulate a multi-game which is a new game-theoretic framework that we introduced in [START_REF] Hamidouche | A multi-game framework for harmonized lte-u and wifi coexistence over unlicensed bands[END_REF].

A multi-game is a hierarchical game in which multiple interconnected sub-games are formulated. Different from classical Stackelberg games, the games can be of different types and the sets of players are not necessarily disjoint. In our context, the WUs are regarded as leaders in the formulated multi-game, that have the priority in using the unlicensed spectrum. On the other hand, the SBSs are considered as followers that adapt their traffic on unlicensed spectrum based on the leaders' traffic load.

Followers' Game

A set of SBSs compete over the unlicensed channels to transmit their content over the unlicensed bands with a predefined QoS requirement depending on the type of the served data. The goal of each SBS i is to allocate its transmission power on each channel to serve its requests of a total size l i during a limited fraction of time that is constrained by the required time for the WUs to transmit their content successfully. We define p min i = [p min i1 , ..., p min iC ] as the minimum required transmit power of an SBS i over all the channels to meet the delay constraint, given the power alloca-tion of other SBSs over the same unlicensed channels. The goal of each SBS is to minimize its transmit power over the unlicensed channels while accounting for the served WUs over the same channels. More formally, this problem can be defined as follows:

p min i = arg min p i ∈R C + C ∑ c=1 p ic (p -i ), (4a)
subject to

l ic w c log(1 + γ(p ic , p -ic )) ≤ T max -t wc , (4b)
∑ c p ic ≤ p max i , (4c)
p ic ≤ p max c , c ∈ C, (4d)
where p -i is the power allocation of all other SBSs except SBS i over all the channels c ∈ C, and l ic is the amount of traffic that SBS i wants to transmit over the unlicensed channel c. The power allocation p ic of SBS i over each channel c can be obtained by solving the following problem:

p ic = arg min {p ic } C ∑ c=1 p ic (p -i ), (5a)
s.t. R ic = w c log(1 + γ(p ic , p -ic )) ≥ l ic T max -t wc , ∑ c p ic ≤ p max i , p ic ≤ p max c ,
where p max c is the maximum transmit power that is allowed over the unlicensed channel c, and t wc is the fraction of time the WUs require to transmit their packets successfully over the selected unlicensed channel c. R ic is the mean achievable rate by the LTE users served by SBS i over the unlicensed channel c as defined in [START_REF] Zhang | Coexistence of wi-fi and heterogeneous small cell networks sharing unlicensed spectrum[END_REF] and can be rewritten as:

R ic = ω c log 1 + p ic |h ic | 2 N + ∑ l∈S\i p lc |h lc | 2 , (6)
where the noise plus interference term from the WUs

N = σ 2 + ∑ j∈W α jc p jc |h jc | 2 ,
can be seen as a constant, which is independent of the power allocation of the SBSs. This is because the WUs define choose their channels in an a the leaders' game and the SBSs respond to these actions by determining their power allocation. We consider the noise plus interference term from WUs as a given term N in (6) because the interference from the WUs transmitting over the same channel c as the SBS i is a constant, in the considered power allocation problem.

The power allocation problem (5a) is equivalent to the following problem:

p ic = arg min {p ic } ∑ c p ic (7a) s.t. w c log(1 + γ(p ic , p -ic )) = l ic T max -t wc , ∑ c p ic ≤ p max i , p ic ≤ p max c ,
where equality holds for the rate requirements in the first constraint in [START_REF] Gu | Exploiting student-project allocation matching for spectrum sharing in LTE-unlicensed[END_REF]. This problem can be formulated as a noncooperative game given by:

G(t w ) = {S, {P i } i∈S , {v i (p i , p -i )} i∈S }, (8)
where the set of SBSs S corresponds to the set of players, v i (p i , p -i ) is the cost function per SBS and given by:

v i (p i , p -i ) = ∑ c∈C p ic , (9)
and P i is the strategy set of SBS i given by:

P i (p -i ) = p ic ∈ R + : ∑ c∈C p ic ≤ p max i , p ic ≤ p max c , l i w log(1 + γ(p i , p -i )) ≤ T max -t w , ∀c ∈ C , (10)
and v i (p i , p -i ) is the cost function of SBS i that corresponds to the transmitted power as defined by the objective function [START_REF] Gu | Exploiting student-project allocation matching for spectrum sharing in LTE-unlicensed[END_REF]. Due to the dependence between the player's i strategy set P i (p -i ) on other players' strategies p -i , the solution of the game G is the Debreu equilibrium also known as the generalized Nash equilibrium [START_REF] Debreu | A social equilibrium existence theorem[END_REF], which can be defined as follows: Definition 1. A strategy profile p DE is a Debreu equilibrium (DE) of the game G(t w ) if, for all the SBSs i ∈ S, we have

p DE i ∈ P i (p -i ) with v i (p DE i , p DE -i ) ≤ v i (p i , p DE -i ), (11)
for all p i ∈ P(p DE -i ).

When all the SBSs except SBS i select their transmit power strategies p -i , the best choice of SBS i consists in a selected strategy from the best-response strategies set B(p -i ) which is given by:

B i (p -i ) = arg min p i ∈P i (p -i ) v i (p i , p -i ). (12)
Thus, the DEs can be derived by solving the fixed point problem [START_REF] Lingzhi | Online channel selection and user association in high-density wifi networks[END_REF] for all the S SBSs. The resulting set of DEs might be empty as it might contain many DEs. Therefore, to analyze the hierarchical game, we first need to establish the existence and uniqueness of a DE in the formulated generalized noncooperative game. The time constraint bound in this case:

w c log(1 + γ(p ic , p -ic )) = l ic T max -t wc . (13)
Thus, the minimum required transmission power of an SBS i to serve its requests without exceeding the time limit is given by:

p * ic = min {p max ic , p DE ic }, (14)
where, p DE ic = arg min

p ic ∈R + v i (p ic , p -i ). (15)
The following proposition provides the closed-form solution of the optimal power allocation and shows that this DE is unique.

Proposition 1. The DE power allocation for SBS i over channel c to its served LTE users is given by:

p DE ic = N |h ic | 2 • 1 -2 -α ic ∑ S j=1 2 -α jc -S + 1 . ( 16
)
Given the single SBS power constraint p max i and the power constraint on each unlicensed channel c, the rate requirement is only feasible if the following conditions are met. Proposition 2. The DE power allocation p ic is feasible under the power constraints if and only if,

∑ c N |h ic | 2 • 1 -2 -α ic ∑ S j=1 2 -α jc -S + 1 ≤ p max i , (17)
and, max

1≤ j≤S N(1 -2 -α jc ) p max c |h ic | 2 + S -1 < S ∑ j=1 2 -α jc < S. (18)
Given the existence and the uniqueness of the DE for the power allocation game on the unlicensed channels, the goal of the WUs is to choose the unlicensed channels over which they serve their traffic. Next, we formulate and solve this problem as the high-level, leaders game.

Leaders' Game

Given the power allocation at the SBSs on each channel, the WUs have to choose the channels over which they transmit their traffic. This WUs-channels allocation problem can be formulated as a house allocation problem with existing tenants [START_REF] Abdulkadiroglu | House allocation with existing tenants[END_REF].

The house allocation problem with existing tenants is one-sided matching that is composed of a set of houses and a number of agents that want to rent a house while some of the houses are already occupied and their tenants may choose to participate or not in the assignment process. Similarly, some agents may already have a house while others do not. In our context, the WUs are the agents and the unlicensed channels correspond to the houses. The unlicensed allocation game can be defined as a tuple P = (S, C ∪ {c 0 }, {κ wc } ∀w∈W,c∈C , {≺ w } w∈W , µ) with each element defined as follows:

• The set W of WUs represents the set of players.

• The set A = C ∪ {c 0 } of unlicensed channels is the set of actions that can be selected by each of the WUs, where c 0 corresponds to the case in which none of the unlicensed channels in chosen by the considered WU. • v wc is the utility of WU w when serving its traffic over unlicensed channel c. • ≺ w is the preference relation of the WUs. The preference relation ≺ w is transitive and complete. We use c ≺ w c to denote that WU w prefers to serve its content over channel c than serving it over channel c for c = c . • µ is the result of the actions selected by all the WUs.

The house allocation problem is a one-sided matching in which only the WUs have preferences over the channels while the unlicensed channels do not participate in the game by taking strategic decisions. This characteristic captures the fact that unlicensed spectrum is free and can be accessed by any user and thus, there does not exist any entity that can act on behalf of the unlicensed bands. The outcome of the house allocation problem with existing tenants can be defined as follows.

Definition 2. A matching between the WUs and unlicensed bands problem µ is a mapping from the set W ∪ C into the set W ∪ C such that for every w ∈ W and

c ∈ C: i) µ -1 (w) is contained in C and µ(c) is contained in W, ii) |µ -1 (w)| ≤ 1 for all w ∈ W, iii) |µ(c)| ≤ q c for all c ∈ C, iv) c ∈ µ -1 (w) if and only if w ∈ µ(c),
where q c is the maximum number of WUs that can be served over channel c.

The value of q c is not predefined at the channels and depends on the amount of LTE traffic that each WU decides to serve over that channel as well as the amount of traffic that each WU decides to serve over the unlicensed channels. Definition 2 states that a WU w can only select one unlicensed band µ -1 (w) while an unlicensed band c can serve multiple WUs µ(c), depending on its capacity and the WiFi traffic load. Before setting the assignment of the WUs to the unlicensed channels, each WU needs to specify its preferences over the unlicensed channels based on its utility function. The externalities in the formulated matching problem appear in the throughput of a given WU that depends on the assigned users to each channel.

The goal of each WU w is to serve its content within the time duration T max . The utility of a WU w when transmitting over channel c is given by:

v wc (c, p c , µ(c)) = T max -tic (p c ) -t wc (µ(c), S c ), (19)
where p c = [p 1c , ..., p Sc ] is the transmit power of all the SBSs over channel c, and t wc (µ(c), S c ) is given in [START_REF] Liu | Small cell traffic balancing over licensed and unlicensed bands[END_REF]. Assuming all the SBSs transmit at the same time, tic is the maximal duration during which the SBSs transmit and is given by: tic (p c ) = arg max

t ic {t ic = l i w c log(1 + γ c (p c )) , ∀i ∈ S c }, (20)
where S c is the set of SBSs that decide to transmit over the unlicensed channel c and t ic is the fraction of time during which an SBS i uses channel c. From (19), we can see that the utility of a WU w not only depends on the set of WUs µ(c) that are assigned to channel c, but also on the interference generated by other SBSs transmitting over the same channel. Based on the defined utility function, each WU can define its preference relation ≺ w over the set of channels C, such that for any two channels c, c ∈ C with c = c , and two matchings µ, µ ∈ W × C, s ∈ µ(c), w ∈ µ (c ):

(µ, c) ≺ w (µ , c ) ⇔ v wc (a -w , µ, S c ) < v wc (a -w , µ , S c ), (21)
To solve the matching problem, we are interested in finding a desirable matching outcome that is Pareto optimal and in the core that can be defined as follows.

Definition 3. A matching µ is Pareto optimal if there does not exist another matching µ under which, at least one of the WUs can improve it utility while none of the WUs will degrade their utility under matching µ compared to matching µ. Definition 4. A matching µ is in the core of the one-sided matching (W, C, ≺ , {u wc } ∀w∈W,c∈C∪c 0 ), if there is no coalition of WUs, W ⊆ W, and a matching µ such that: i) µ -1 (W ) ∈ {c 0 } ∀i∈W for all s ∈ W , ii) µ -1 (w) w µ -1 (w) for all w ∈ W , iii) µ -1 (w) ≺ w µ -1 (w) for some w ∈ W .

In the formulated problem, the channels can serve multiple WUs while in the original house allocation problem only one agent can be assigned to a given house. Thus, the existing algorithms do not account for the externalities that appear in the utility function (19) of the WUs as it depends on the WUs that are assigned to each of the channels. Moreover, the existing top trading cycle algorithm [START_REF] Abdulkadiroglu | House allocation with existing tenants[END_REF] for solving such problems is centralized and cannot be applied for the assignment of WUs to the unlicensed channels in which the algorithms should be decentralized due to the large number of access points. Next, we propose a new decentralized algorithm to solve the formulated one-sided matching game with externalities.

Proposed Algorithm

In the proposed algorithm, we solve the conflict between the WUs that want to access the same channel, by ranking the WUs randomly. Following the defined order, if a WU is not associated to its most preferred unlicensed channel that allows it to strictly increase its current utility, the WU w 0 sends a request to one of the WUs that is assigned to its most preferred channel called w 1 . Upon receiving a request, if WU w 1 is not matched to its most preferred unlicensed channel, it sends a request to a WU w 2 that is assigned to its most preferred channel. It also includes a list that contains the WU w 1 from which it received a request and WU w 0 that has sent a request to WU w 1 , thus, including all the WUs that lead to the initial requesting WU. At the end of the requests process, each WU checks the existence of a cycle in the received list from its requesting WUs. Once all the WUs emit their requests, the WUs that have initiated the requests procedure, check if a cycle is detected in the list they receive. If so, the WU with the highest order accepts the requests of its preceding WU in the list and transmit the list by removing itself from it. Each WU in the cycle does the same thing until the last WU of the list. Then, all the WUs that belong to a cycle leave the matching game and the WUs update their preferences list based on the remaining WUs. The assignment process is then repeated among the remaining WUs and their associated channels.

Theorem 1. The WUs-channels assignment that results from the proposed algorithm is in the core.

Theorem 2. The outcome of the proposed matching algorithm is the unique matching in the core.

Theorem 3. The proposed algorithm is Pareto efficient.

Conclusions

In this paper, we have addressed the problem of coexistence over unlicensed bands. In particular, we have formulated a multi-game in which the WUs are considered as leaders that have priority when accessing the unlicensed bands, and the SBSs as followers. In the leaders game, the WUs select the unlicensed channel using which they transmit their content while anticipating the SBSs' possible reactions. In the followers game, the SBSs follow the leaders decisions and respond to it by determining their transmit power over each of the channels.


Hyesung Kim
email: hskim@ramo.yonsei.ac.kr

Jihong Park
email: jihong@es.aau.dk

Mehdi Bennis
email: bennis@ee.oulu.fi

Seong-Lyun Kim
email: slkim@ramo.yonsei.ac.kr

Mérouane Debbah
email: merouane.debbah@huawei.com

Ultra-Dense Edge Caching under Spatio-Temporal Demand and Network Dynamics

Keywords: Edge caching, mean-field game, stochastic geometry, spatio-temporal dynamics, ultra-dense networks, 5G

This paper investigates a cellular edge caching design under an extremely large number of small base stations (SBSs) and users. In this ultra-dense edge caching network (UDCN), SBS-user distances shrink, and each user can request a cached content from multiple SBSs. Unfortunately, the complexity of existing caching controls' mechanisms increases with the number of SBSs, making them inapplicable for solving the fundamental caching problem: How to maximize local caching gain while minimizing the replicated content caching? Furthermore, spatial dynamics of interference is no longer negligible in UDCNs due to the surge in interference. In addition, the caching control should consider temporal dynamics of user demands. To overcome such difficulties, we propose a novel caching algorithm weaving together notions of mean-field game theory and stochastic geometry. These enable our caching algorithm to become independent of the number of SBSs and users, while incorporating spatial interference dynamics as well as temporal dynamics of content popularity and storage constraints. Numerical evaluation validates the fact that the proposed algorithm reduces not only the long run average cost by at least 24% but also the number of replicated content by 56% compared to a popularity-based algorithm.

I. INTRODUCTION

Upcoming 5G systems are expected to become ultra-dense networks (UDNs) due to relentless user demand growth [START_REF] Kamel | Ultra-Dense Networks: A Survey[END_REF], [START_REF] Samarakoon | Ultra Dense Small Cell Networks: Turning Density Into Energy Efficiency[END_REF]. Limited backhaul capacity is however unable to cope with a huge number of small base stations (SBSs) in a UDN. In this respect, edge caching is a promising UDN enabler that alleviates backhaul congestion by storing popular content during off-peak hours at the network edges, i.e. SBSs [START_REF] Bas ¸tuǧ | Living on the Edge: The Role of Proactive Caching in 5G Wireless Networks[END_REF]- [START_REF] Tamoor-Ul-Hassan | Caching in Wireless Small Cell Networks: A Storage-Bandwidth Tradeoff[END_REF], yielding to an ultra-dense edge caching network (UDCN). This paper aims at proposing a UDCN caching control algorithm by solving the following three fundamental issues in UDCNs.

1. Caching control complexity. Edge caching network performance relies on how much amount of users request the content cached in SBSs, i.e. cache hits. If user demand is fully correlated, caching the most popular files would be an optimal control policy. However, user demand in reality is heterogeneous, and therefore leads to the fundamental tradeoff of edge caching: maximizing individual user's cache hits while minimizing replicated content caching at SBSs. Solving this trade-off in a centralized manner is too complicated to be implemented in practice, so preceding works suggest semidistributed caching control algorithms where each SBS only exploits the local caching decisions of neighboring SBSs [START_REF] Wang | Dynamic Femtocaching for Mobile Users[END_REF]. Unfortunately, the number of neighboring SBSs in a UDCN is much larger than that of a traditional network, and thus existing algorithms are computationally infeasible under ultradense SBS environment. To make these worse, the following two spatio-temporal aspects of UDCNs further increase the caching control complexity.

2. Temporal dynamics of user demand. Content popularity among users evolves over time, so demand prediction error exists when caching files at SBSs. Minimizing this prediction error straightforwardly brings about huge complexity increase, but nevertheless UDCN caching control cannot overlook this misprediction impact. The reason is that numerous users amplify individual misprediction loss in a UDCN. A large number of SBSs also lead to significant redundant caching due to misprediction, increasing congestion at both backhauls and cache storages.

3. Spatial dynamics of interference. In the downlink, interference amount and variance increase linearly with the number of SBSs [START_REF] Haenggi | Stochastic Geometry for Wireless Networks[END_REF]. This significant impact of interference should thus be carefully considered in UDCN caching control design. Many preceding works apply somewhat simplified interference modelings such as the average interference [START_REF] Liu | Will Caching at Base Station Improve Energy Efficiency of Downlink Transmission?[END_REF] with fixed SBS locations. Incorporating the interference distribution that specifies random locations of SBSs however incurs additional complexity in UDCN caching control.

On the basis of these issues, our objective is to propose a UDCN caching control algorithm that incorporates spatiotemporal demand, cache storage and interference dynamics with low computational complexity while minimizing the long run average (LRA) cost. The cost is defined such that it increases with (i) backhaul use, (ii) cache storage use, and (iii) content replication, while decreasing with (iv) spectral efficiency.

To this end, we tackle the aforementioned triple difficulty in reverse order by using stochastic geometry (SG) and meanfield game (MFG) theory. SG first plays a key role to specify the interference dynamics at a randomly selected user, i.e. a typical user. In UDNs, the interference normalized by BS densification becomes deterministic [START_REF] Park | Spatio-Temporal Network Dynamics Framework for Energy-Efficient Ultra-Dense Cellular Networks[END_REF], [START_REF] Park | User-Centric Mobility Management in Ultra-Dense Cellular Networks under Spatio-Temporal Dynamics[END_REF], which allows us to simplify the spatial dynamics of UDCNs. Based on this result, MFG captures temporal dynamics of demand and its corresponding caching behavior by using two partial differential equations (PDEs), Hamilton-Jacobi-Bellman equation (HJB) and Fokker-Planck-Kolmogorov equation (FPK). In a UDCN, remarkably, no matter how many SBSs and users exist, only a single pair of HJB-FPK can describe the entire caching control and the resultant LRA cost dynamics [START_REF] Guéant | Mean field games and applications[END_REF], achieving fixed computational complexity. Solving these HJB and FPK equations provides the desired optimal UDCN caching control.

Related works. A recent study on caching [START_REF] Perabathini | Caching at the Edge: a Green Perspective for 5G Networks[END_REF] considers interference dynamics using SG for the spatial domain. Caching policy adaptive to time-varying user request [START_REF] Gregori | Wireless Content Caching for Small Cell and D2D Networks[END_REF] is proposed to enhance the local caching gain. The authors of [START_REF] Hamidouche | Mean-Field Games for Distributed Caching in Ultra-Dense Small Cell Networks[END_REF] considered temporal dynamics of SBSs' states, replicated contents, and interference on caching control in UDN environment. They, however, neglected temporal evolution of content popularity and spatial dynamics of interference determined by user locations. None of the works on caching jointly considers the spatio-temproal dynamics of user demand and network (interference and SBSs' backhaul and storage capacities) in ultra-dense scenarios.

Under these spatio-temporal dynamics, it is hard to improve the caching gain while minimizing the replicated content caching [START_REF] Chun | Selfish Caching in Distributed Systems: A Game-Theoretic Analysis[END_REF]. The replication problem becomes more severe in UDCN, since the increased number of candidate SBSs serving the same user induces that only one of them transmits the cached data. Furthermore, the amount of replicated contents varies according to the aforementioned dynamics. A caching strategy dealing with the dynamically evolving replication remains an open problem in UDNs.

Contributions. The contributions of this paper are summarized as below.

• A novel caching control algorithm for UDCNs is proposed, which minimizes LRA cost while incorporating spatio-temporal demand and interference dynamics with fixed computational complexity, independent of the numbers of SBSs and users (see Proposition 1). Specifically, the algorithm does not need to get full knowledge of other SBSs' caching strategies or storage states. • The proposed algorithm's gains in LRA cost and replicated contents are numerically validated (Figs. 5 and6). • The unique mean-field equilibrium (MFE) [START_REF] Guéant | Mean field games and applications[END_REF] is verified by numerical evaluation (Fig. 4) as well as by analysis (Proposition 1).

II. NETWORK MODEL AND SPATIO-TEMPORAL DYNAMICS

A. Ultra-dense network

We consider a UDN with SBS density λ b and user density λ u . Their locations follow independent homogeneous Poisson point processes (PPPs), respectively. We set λ b λ u , according to the definition of UDN [START_REF] Park | Tractable Resource Management with Uplink Decoupled Millimeter-Wave Overlay in Ultra-Dense Cellular Networks[END_REF]. User i receives signals from SBSs within a reception ball b(y i , R) centered at y i with radius R, which represents the average distance that determines a region where the received signal power is larger than noise floor, as depicted in Fig. 1. When R goes to infinity, the reception ball model is identical to a conventional PPP network.

Transmitted signals from SBSs experience path-loss attenuation. The attenuation from the k-th SBS coordinates z k to the i-th user coordinates

y i is l k,i = min(1, ||z k -y i || -α ),
where α is the path-loss exponent. The transmitted signals experience independent and identically distributed fading with the coefficient g k,i (t). We assume that the coefficient is not temporally correlated. The channel gain

h k,i from SBS k to user i is |h k,i (t)| 2 = l ki |g k,i (t)| 2 . The received signal power is given as S(t) = P |h k,i (t)| 2
, where P denotes transmit power of an SBS.

The SBS directionally transmits signal by using N a number of antennas. Its beam pattern at a receiver follows a sectored uniform linear array model [START_REF] Venugopal | Interference in Finite-Sized Highly Dense Millimeter Wave Networks[END_REF] where the main lobe gain is N a with beam width θ Na = 2π/ √ N a assuming that side lobes are neglected and the beam center points at the receiver.

B. Edge caching model

Let us assume that there is a set N consisting of N SBSs within the reception ball with radius R. Each SBS k ∈ N has data storage unit size C k,j assigned for content j. The storage units allow SBSs to download contents a priori from a server connected with non-ideal backhaul link as depicted in Fig. 1. We assume that there is a set M consisting of M contents, and the server has all the contents files in the set M. Users request content j from the contents set M with probability x j , and the size of the content j is denoted by L j . The goal of SBS k is to determine a vector of caching probabilities p k (t) = {p k,1 (t), ..., p k,j (t), ...p k,M (t)}, where p k,j (t) ∈ [0, 1] is a probability that SBS k downloads content j through its own backhaul link at time t. The control variable p k,j (t) can also be interpreted as a fraction of the file in case that each content file is encoded using a maximum distance separable dateless code [START_REF] Bioglio | Optimizing MDS Codes for Caching at the Edge[END_REF].

According to the definition of UDN, multiple SBSs can serve a user [START_REF] Park | Tractable Resource Management with Uplink Decoupled Millimeter-Wave Overlay in Ultra-Dense Cellular Networks[END_REF] unlike traditional unidentified networks as show in Fig. 1. We assume that a reference user is associated to one of the SBSs storing the requested content within the threshold distance R. If there are multiple SBSs having the same content within the region, the serving SBS is randomly chosen. If no SBS has cached the requested content, one of the SBSs within R fulfills the user request by downloading it from the server through the backhaul.

C. Spatio-temporal dynamics of user demand and network

Evolution law of time-varying content popularity. The probability that users request a content, representing popularity, varies over time in reality. We use Ornstein-Uhlenbeck process [START_REF] Mériaux | Stochastic Differential Games and Energy-Efficient Power Control[END_REF] to model this temporal fluctuation as a stochastic differential equation (SDE):

dx j (t) = (u j -a j )dt + ηdW j (t), (1)
where x j (t) is the probability that users request content j at time t, u j and a j respectively denote popularity increment and decrement for content j, η is a positive constant, and W j (t) is a Wiener process.

Temporal dynamics of cache storage size. The remaining cache storage capacity varies according to the instantaneous caching policy. Let us assume that SBS k discards content files at a rate of µ k,j in order to make a space for downloading other data. Considering the discarding rate, we model the evolution law of storage unit as follows

dQ k,j (t) = (µ k,j -L j p k,j (t))dt, (2)
where Q k,j (t) denotes the remaining storage size dedicated to content j of SBS k at time t, and L j is data size of content j. L j p k,j (t) represents instant data size of content j downloaded by SBS k at time t. Spatial dynamics of interference. In UDN, an SBS having no associated user within its coverage becomes dormant, not transmitting any signal. Locations of users determine activation of SBSs, characterizing spatial distribution of interference. Consider a randomly selected typical user. We assume that active SBSs have always data to transmit regardless of their own caching policy. It indicates that the aggregate interference is imposed by active SBSs with active probability p a . Assuming that p a is homogeneous over SBSs yields [START_REF] Yu | Downlink Capacity and Base Station Density in Cellular Networks[END_REF]. It provides that density of arbitrary interfering SBSs is equal to p a λ b . Then, active SBSs' coordinates, consisting a set Φ R (p a λ b ), determine dynamics of interference I f (t).

p a ≈ 1 -[1 + λ u /(3.5λ b )] -3.5
I f (t) = |Φ R (paλ b )| k P |h k,i (t)| 2 , (3)
The signal-to-interference-plus-noise (SINR) with N a number of transmit antennas is

SINR(t) = NaP |h(t)| 2 / σ 2 + θN a 2π NaI f (t) . (4)
In the following section, we present the spatially averaged version of I f (t).

III. PROBLEM FORMULATION

The goal of each SBS k is to determine its own caching probability p * k,j (t) for content j in order to minimize a long run average (LRA) cost. The LRA cost is determined by other SBSs' caching policies, states of content request probability, wireless channel, backhaul capacity, and remaining storage size. Caching strategies of other SBSs also vary according to the spatio-temporal dynamics.

A. Interactions: Replication and Interference

There are inherent interactions among SBSs with respect to their own caching strategies. These interactions depend on replicated contents and interference, which are major bottleneck for optimizing distributed caching. Our purpose is to estimate these interactions in a distributed fashion without full knowledge of other SBSs' states or actions.

Replication. As shown in Fig. 1, there may be replicated contents downloaded by multiple SBSs within R. The overlapping contents increases redundant cost due to inefficient resource utilization [START_REF] Chun | Selfish Caching in Distributed Systems: A Game-Theoretic Analysis[END_REF]. It is worth noticing that the number of replicated contents is determined by other SBSs' caching strategies. We define the replication function I r k,j (p -k,j (t)) as the expected value of the overlapping content size per unit storage size C k,j .

I r k,j (p -k,j (t)) = 1 C k,j N r(j) |N | i =k p i,j (t), (5)
where p -k,j (t) is a vector of caching probabilities of all the other SBSs except SBS k, and N r(j) denotes the number of contents whose request probability is asymptotically equal to that of content j. Specifically, N r(j) is the cardinality of a set {m|m ∈ M such that |x m -x j | ≤ } where is sufficiently small. Interference. Interactions among SBSs through interference can be bottlenecks for optimizing distributed caching. To incorporate this spatial interaction, we adopt an analysis of [START_REF] Park | Spatio-Temporal Network Dynamics Framework for Energy-Efficient Ultra-Dense Cellular Networks[END_REF] on interference in UDN. It provides Îf (t), the interference normalized by SBS density and the number of antennas as follows:

Îf (t) = (λ u πR) 2 N -1 2 a λ -α 2 b 1+ 1 -R 2-α α -2 P E g [|g(t)| 2 ]. (6)
It gives us an average downlink spectral efficiency (SE) R k (t) in UDN as follows:

R k (t) = E S,I f [log(1+SINR(t))] ≈ ES log   1+ S k (t) σ 2 Naλ α/2 b +E I f[ Îf (t)]   , (7)
where σ 2 is the noise power. Eq. ( 7) allows us to investigate the effect of interference on the upper bound of an average SE.

B. Cost Functions

An instantaneous cost function J k,j (t) defines the LRA cost. It is affected by backhaul capacity, remaining storage size, spectral efficiency, and overlapping contents among SBSs. SBS k cannot download more than B k,j (t), the allocated backhaul capacity for downloading content j at time t.

If p k,j (t) < B k,j (t)
Lj , the backhaul cost φ k,j is model as φ k,j (p k,j (t)) = -log(B k,j (t) -L j p k,j (t)). Otherwise, the cost function φ k,j (p k,j (t)) goes to infinity, preventing the download rate of content j from exceeding the available backhaul rate. This form of cost function is widely used as in [START_REF] Hamidouche | Mean-Field Games for Distributed Caching in Ultra-Dense Small Cell Networks[END_REF].

As cached content files occupy the storage, it causes processing latency [START_REF] Cao | Fundamental Storage-Latency Tradeoff in Cache-Aided MIMO Interference Networks[END_REF] or delay to search requested files by users. We consider this storage cost as follows:

ψ k,j (Q k,j (t)) = γ(C k,j -Q k,j (t))/C k,j . (8)
Then, the global instantaneous cost is

J k,j (p k,j (t),p -k,j (t)) = φ k,j (p k,j (t))(1+I r k,j (p -k,j (t))) R k (t,I f (t))x j (t) + ψ k,j (Q k,j (t)). (9)
With this cost function, we define an average total cost over the long run, called long run average (LRA) cost.

C. Stochastic Differential Game for Edge Caching

The state of SBS k and content j at time t is defined as s k,j (t) = {x j (t), R k (t), Q k,j (t)}, ∀k ∈ N , ∀j ∈ M. The stochastic differential game (SDG) for edge caching is defined by (N , S k,j , A k,j , J k,j ) where S k,j is the state space of SBS k and content j, A k is the set of all caching controls {p k,j (t), 0 ≤ t ≤ T } admissible for the state dynamics, and J k,j is the LRA cost over a time window [0, T ] defined as follows: [START_REF] Park | User-Centric Mobility Management in Ultra-Dense Cellular Networks under Spatio-Temporal Dynamics[END_REF] where κ : S k,j → R, s k → κ(s k ) is the cost of having a remaining storage size at the end of the time window [0, T ]. Thus, the SDG is formulated as follows:

J k,j = E T t J k,j (p k,j (t), p-k,j(t)) dt + κ(s k (T )) ,
(P1) v k,j (t) = inf p k,j (t) J k,j (t). (11)
subject to dx j (t) = (u j -a j )dt + ηdW j (t),

dQ k,j (t) = (µ k,j -L j p k,j (t))dt. ( (12)
) 13
The existence of a Nash equilibrium of problem P1 is guaranteed if there exists a joint solution of the following coupled HJB equations for all k and j [START_REF]Stochastic Differential Equations[END_REF]:

0 = ∂tv k,j (t)+ inf p k,j (t) J k,j (p k,j (t), p-k,j(t))+ η 2 2 ∂ 2 xx v k,j (t) +(µ k -Ljp k,j (t)) (A) ∂Q k v k,j (t) + (uj -aj) (B) ∂xv k,j (t) . (14)
If the smoothness of the drift functions (A) and (B) in the dynamic equation ( 14) and the cost function [START_REF] Park | Spatio-Temporal Network Dynamics Framework for Energy-Efficient Ultra-Dense Cellular Networks[END_REF], we can assure that a unique solution of the equation ( 14) exists [START_REF]Stochastic Differential Equations[END_REF]. Unfortunately, it is complex to solve the coupled N × M HJB equations. We thus release complexity of this system in the next section.

D. Mean-field Game for Caching

It is hard to react to each individual SBS's caching strategy in solving the SDG P1. Mean-field game (MFG) theory enables us to transform these multiple interactions into a single interaction, called MF interaction, via MF approximation. According to [START_REF] Mériaux | Stochastic Differential Games and Energy-Efficient Power Control[END_REF], this approximation holds under the following conditions: a large number of players, excahngeability of players under the caching control and finite MF interaction.

Remark that the first condition corresponds to the definition of UDNs. Players in the SDG are said to be exchangeable or indistinguishable under the control p k,j (t) and the states of players and contents, if the player's control is invariant by their indices and decided by only their individual states. In other words, permuting players' indices cannot change their control policies. Under this exchageability, we can focus on a generic SBS by dropping its index k.

Respective MF interactions ( 5) and ( 6) should asymptotically converge to a finite value under the above conditions. MF replication (5) goes to zero when the number of contents per SBS is extremely large, i.e. M N . Such a condition implies that the cardinality of the set consisting of asymptotically equal content popularity goes to infinity. In other words, N r(j) goes to infinity yielding that I r k,j (p -k,j (t)) becomes zero. In terms of interference, MF interference converges as the ratio of SBS density to user density goes to infinity, i.e. N a λ α b /(λ u R) 4 → ∞ [START_REF] Park | Spatio-Temporal Network Dynamics Framework for Energy-Efficient Ultra-Dense Cellular Networks[END_REF]. Such a condition corresponds to the notion of UDN [START_REF] Park | Tractable Resource Management with Uplink Decoupled Millimeter-Wave Overlay in Ultra-Dense Cellular Networks[END_REF] or massive MIMO (N a → ∞). That is, the conditions enabling MF approximation inherently hold under ultra-densified caching networks.

To approximate interactions from other SBSs, we need a state distribution of SBSs and contents at time t, called MF distribution m t (x, Q). It is defined by a counting measure M

(N ×M ) t (x, Q) = 1 N M M j=1 N k=1 δ {xj (t),Q k (t)} .
We assume that the empirical distribution M (N ×M ) t (x, Q) converges to m t (x, Q), which is the density of contents and SBSs in state (x, Q). Note that we omit the SE R(t) from the density measure to only consider temporally correlated state without loss of generality.

The MF distribution m t (x, Q) is a solution of the following Fokker-Planck-Kolmogorov (FPK) equation.

0 = ∂ t m t (x, Q) -(u j -a j )∂ x m t (x, Q) + (µ k -L j p j (t))∂ Q m t (x, Q) + η 2 2 ∂ 2 xx m t (x, Q). ( 15
)
Let us denote the solution of the FPK equation ( 15) as m * t (x j , Q). Exchangeability and existence of the MF distribution allow us to approximate the interaction I r k,j (p -k,j (t)) as a function of m * t (x j , Q) as follows:

I r j (t, m * t (xj, Q)) = Q x m * t (xj, Q) C k,j N r(j) pj(t, x, Q)dxdQ, (16)
Remark that we can estimate the interaction from replication (16) without observing other SBSs' caching strategies. Thus, it is not necessary to have full knowledge of the states or the caching control policies of other SBSs. It provides that an SBS only need to solve a pair of equations, the FPK equation ( 15) and the following modified HJB one from ( 14):

0 = ∂tvj(t)+ inf p j (t) Jj(pj(t), Ij(t, m * t (x, Q))+ η 2 2 ∂ 2 xx vj(t) +(µ -Ljpj(t))∂Qvj(t) + (uj -aj)∂xvj(t) . (17)
We can solve this equation in backward and find the solution of the stochastic optimization problem P1. It is worth mentioning that the number of PDEs to solve for one content is reduced to two from N . The respective processes of solving P1 in ways of SDG and MFG are depicted in Fig. 2.

Proposition 1. The optimal caching probability is given by:

p * j (t) = 1 L j B j (t) - 1 + I r j (t, m * t (x j , Q)) R(t, I f (t))x j (t)∂Qv * j + , (18)
where m * t (x, Q) and v * j (t) are the unique solutions of (15) and (17), respectively. Proof: Appendix.

The optimal caching probability p * j (t) is calculated by water-filling fashion, where the allocated backhaul capacity B j (t) determines the water level. Noting that SE R(t) increases with the number of antenna N a and SBS density λ b , SBSs cache more contents from server when wireless capacity is improved. The partial derivative ∂Qv j implies the effect of the current decision on the final LRA cost.

The respective solutions m * t (x, Q) and v * j (t) of ( 15) and ( 17) define the mean field equilibrium (MFE) which corresponds to the notion of Nash equilibrium in the MFG framework [START_REF] Guéant | Mean field games and applications[END_REF]. Note that existence and uniqueness of the optimal caching control is guaranteed. It verifies that the optimal caching algorithm converges to a unique MFE, when the initial conditions m 0 , x j (0), and Q(0) are given.

The optimal caching algorithm in Proposition 1 becomes more tractable when each SBS affords high storage capacity enough to ignore cost for storage usage. Thus, we can obtain the following corollary of Proposition 1.

Corollary 1. When the storage capacity is sufficiently high (C k,j → ∞), the optimal caching control is given by

p * j (t) = 1 L j B j (t) - 1 + I r j (t, m * t (x j , Q)) R(t, I f (t))x j (t)γ + , (19)
where γ is a constant.

Proof: When the storage capacity goes to infinity, the instantaneous cost function for storage usage ( 8) is given as follows: lim Cj →∞ ψ j (Q j (t)) = γ. Hence, the partial derivative of LRA cost with respect to the remaining storage size ∂Qv j also has a constant value of γ.

It is worth noticing that the optimal caching control becomes independent from the LRA cost dynamics with respect to the storage state.

IV. NUMERICAL RESULTS

In this section, we present numerical results for cases in which content request probability remains static and temporally varies according to its evolution law [START_REF] Kamel | Ultra-Dense Networks: A Survey[END_REF]. Let us assume full knowledge of contents request probability and their evolution laws. The initial distribution of the SBSs m 0 is given as a normal distribution N (0.7, 0.05 2 ). We assume that the storage size Q belongs to a set [0, 1]. Assuming Rayleigh fading with mean one, we set the parameters as follows: µ = 0.1, γ = 0.01, λ u = 0.001, λ b = 0.03, R = 10/ √ π, B = 1, N r(j) = 20, Q(0) = 0.7, η = 0.1, α = 4.

A. Static content request probablity

We numerically analyze the algorithm performance when the content request probability is static. In this case, the problem becomes simplified, because the caching control strategies does not depend on the evolution law of the content popularity. Specifically, in HJB [START_REF] Venugopal | Interference in Finite-Sized Highly Dense Millimeter Wave Networks[END_REF] and FPK [START_REF] Chun | Selfish Caching in Distributed Systems: A Game-Theoretic Analysis[END_REF] equations, the derivative terms with respect to content request probability x become zero.

Fig. 3 shows the evolution of the optimal caching probability p * (t) for a content with respect to storage state. We observe that the value of p * (t) is lower than the content request probability at all time slots and decreases when available storage space Q becomes small. The reason is that SBSs restrain replication of the content and save storage capacity for downloading more popular contents. Fig. 4 represents evolutions of MF distribution m * t (Q) for a content with different popularity. Note that m * t (Q) indicates the density of SBSs whose unoccupied storage size is equal to Q at time t. We observe that the unoccupied storage space of SBSs does not diverge from each other. It validates that the proposed algorithm achieves the MFE. In this equilibrium, the amount of downloaded content file decreases when the content popularity x becomes low. This tendency corresponds to the trajectory of the optimal caching probability in Fig. 3. Almost every SBS has downloaded the content over time, but not spent its whole storage. The remaining storage saturates even though the content popularity is equal to 0.9. It implies that SBSs reduce the downloading amount of the popular content in consideration of the expected replication amount, which increases with the popularity.

B. Dynamic content request probablity

We investigate the performance of the proposed MF caching algorithm under temporal dynamics of the content popularity [START_REF] Kamel | Ultra-Dense Networks: A Survey[END_REF]. Note that the temporal dynamics constrains the optimal caching strategy to follow the optimal trajectory [v * j (t), m * t (x, Q)], a set of solutions to the coupled equations ( 17) and [START_REF] Chun | Selfish Caching in Distributed Systems: A Game-Theoretic Analysis[END_REF]. Fig. 5 shows the LRA cost of the proposed MF caching algorithm compared to results of a popularity based algorithm and random caching one. The blue dot in the figure indicates the minimal LRA cost obtained by an exhaustive search. A slight disagreement between the proposed and exhaustive algorithm by full search appears from two reasons; one is that each segment size of time t and storage capacity Q is not infinitesimal, leading to inaccurate partial derivative values, and the other is the finite number of SBSs. Nevertheless, the proposed caching strategy reduces about 36% of the LRA cost as compared to an algorithm that determines caching probability in proportional to the request probability. This performance gain results from allowing the system to avoid the expected replication of contents and interference varying according to the spatio-temporal dynamics.

Fig. 6 shows the amount of replicated contents per storage usage as a function of initial content probability x(0). The proposed MF caching algorithm reduces caching content replication averagely 56% compared to popularity based caching.

However, MF caching algorithm yields the higher amount of replication than random caching does when the content request probability becomes high. The reason is that the random policy downloads contents regardless of their own popularity, so the amount of replication remains steady. On the other hand, MF caching increases the downloaded volume of popular contents.

V. CONCLUDING REMARKS

In this paper, we proposed an edge caching algorithm for ultra-dense edge caching networks (UDCNs), taking into account the spatio-temporal user demand and network dynamics. Leveraging the frameworks of MFG theory and SG, our algorithm enables SBSs to distributively determine their own caching strategy without full knowledge of other SBSs' caching strategies or storage state. The proposed MF caching algorithm reduces not only the long run average cost but also the replicated content amount compared to a caching control that is merely based on content popularity. Future work can extend to an optimal caching strategy to cope with imperfect knowledge of content popularity dynamics.

VI. APPENDIX: PROOF OF PROPOSITION 1

The optimal control strategy is the argument of the infimum term of the HJB equation [START_REF] Hamidouche | Mean-Field Games for Distributed Caching in Ultra-Dense Small Cell Networks[END_REF]. inf pj (t) J j (p j (t), I r j (t, m * t (x, Q))+

η 2 2 ∂ 2 xx v j (t) +(µ -L j p j (t))∂ Q v j (t) + (u j -a j )∂ x v j (t) (20)
Solving [START_REF] Yu | Downlink Capacity and Base Station Density in Cellular Networks[END_REF] for all time t is immediate by convexity of the objective function, of which the critical point is the unique optimal control p * j (t) calculated as described in [START_REF] Bioglio | Optimizing MDS Codes for Caching at the Edge[END_REF]. Note that p * j (t) is a function of the solutions of the equations ( 15) and [START_REF] Venugopal | Interference in Finite-Sized Highly Dense Millimeter Wave Networks[END_REF]. The expression of p * j (t) [START_REF] Bioglio | Optimizing MDS Codes for Caching at the Edge[END_REF] provides the final versions of the HJB and FPK equations composing MFG as follows: From these equations, we can find the values of v * j (t) and m * t (x, Q). Note that the smoothness of the drift functions and in the dynamic equation and the cost function [START_REF] Park | Spatio-Temporal Network Dynamics Framework for Energy-Efficient Ultra-Dense Cellular Networks[END_REF] assures the uniqueness of the solution [START_REF]Stochastic Differential Equations[END_REF].

0 = ∂tvj(t

Fig. 1 :

 1 Fig. 1: An illustration of ultra-dense edge caching networks (UD-CNs) with spatio-temproal user demand and network dynamics.

Fig. 2 :

 2 Fig. 2: Ultra-dense edge caching flow charts with and without utilizing MFG and SG.

Fig. 3 :

 3 Fig. 3: Evolution of the optimal caching control p * (t, Q) at the MF equilibrium when xj = 0.4. The initial MF distribution m0 is given as N (0.7, 0.05 2 ).

Fig. 4 :

 4 Fig. 4: Evolution of the MF distribution m * t (Q) that indicates the density of SBSs having the remaining storage space Q at time t with respect to different content popularity x, B = 1, N r(j) = 20, µ = 0.1, λu = 0.001, λ b = 0.03, R = 10/ √ π.

Fig. 5 :

 5 Fig. 5: Long run average costs of different caching strategies (Q(0) = 0.7, x(0) = 0.3, η = 0.1, a = 0.15, u = 0.1).

Fig. 6 :

 6 Fig. 6: The amount of replicated contents per storage usage (Q(0) = 0.7, η = 0.1, a = 0.15, u = 0.1).

  ) -log Bj(t) -Bj(t) -

					1+I r j (t,m * t (x,Q))	+
					R(t,I f (t))x j (t)∂Qv j
				R(t, I f (t))xj(t)
	×(1+I r j (t, m * t (x, Q))) +	α(Cj -Qj(t)) Cj	+ (uj -aj)∂xvj(t)
	+ µ-Bj(t)-	1+I r j (t, m * t (x,Q)) R(t, I f (t))xj(t)∂Qvj	+		∂Qvj(t)+	η 2 2	∂ 2 xx vj(t),
	0 = ∂tmt(x, Q) -(uj -aj)∂xmt(x, Q) +	η 2 2	∂ 2 xx mt(x, Q)
	+ µ k -Bj(t)-	1+I r j (t, mt(x,Q)) R(t, I f (t))xj(t)∂Qv * j		+	∂Qmt(x, Q).

ACKNOWLEDGEMENT

This research was supported by the Korean government through the National Research Foundation (NRF) of Korea Grant (NRF-2014R1A2A1A11053234), the Academy of Finland CARMA project, the FOGGY project and the ERC Starting Grant 305123 MORE (Advanced Mathematical Tools for Complex Network Engineering).