Stochastic Geometric Analysis of Energy-Efficient Dense Cellular Networks

Dense cellular networks (DenseNets) are fast becoming a reality with the rapid deployment of base stations (BSs) aimed at meeting the explosive data traffic demand. In legacy systems however this comes with the penalties of higher network interference and energy consumption. In order to support network densification in a sustainable manner, the system behavior should be made 'load-proportional' thus allowing certain portions of the network to activate on-demand. In this work, we develop an analytical framework using tools from stochastic geometry theory for the performance analysis of DenseNets where load-awareness is explicitly embedded in the design. The model leverages on a flexible cellular network architecture where there is a complete separation of the data and signaling communication functionalities. Using the proposed model, we identify the most energy- efficient deployment solution for meeting certain minimum service criteria and analyze the corresponding power savings through dynamic sleep modes. Based on state-of-the-art system parameters, a homogeneous pico deployment for the data plane with a separate layer of signaling macro-cells is revealed to be the most energy-efficient solution in future dense urban environments.


I. INTRODUCTION
Ultra-dense deployment of base stations (BSs), relay nodes, and distributed antennas is considered a de facto solution for realizing the significant performance improvements needed to accommodate the overwhelming future mobile traffic demand [2]. While legacy wireless communication systems are fast approaching the information-theoretic capacity limits, dense cellular networks (DenseNets) can push data rates further by shortening the transmitter-receiver distance and serving fewer users per cell [3]. The extremely populated topology of DenseNet, a.k.a. heterogeneous cellular network (HetNet), raises several technical challenges, including managing the interference and keeping the energy expenditure in check [4].
Understanding the interference behavior in DenseNets is challenging due to the rapid, irregular, and overlapping placement of nodes. In addition, in contrast to existing macro-cells where different parts of spectrum is allocated to neighboring cells, DenseNets employ an aggressive frequency reuse strategy where different nodes can access the same spectrum; thus highlighting the importance of interference management for facilitating efficient spectrum utilization [5]. On the other hand, legacy cellular networks and transmission technologies are designed and dimensioned to meet the coverage and capacity requirements in peak traffic conditions. This approach threatens the commercial viability of deploying many more network nodes which would substantially increase the total capital, operational, and environmental expenditure [6]. An extensive design overhaul towards a flexible cellular network architecture with load-proportional energy consumption behavior is hence required going forward into the future [7].
In recent years, several collaborative initiatives, such as the GreenTouch consortium [8], have focused on analyzing the achievable spectral and energy efficiency performance of network densification as well as other promising solutions, using mostly system-level simulations. This approach is inline with the traditional system planning where Monte-Carlo simulations are utilized for drawing conclusions on the cellular network performance. However, due to the inherent characteristics of DenseNets, the simulation-based investigations have become extremely resource-intensive. In order to reduce the underlying complexity associated with DenseNet planning, tractable and computationally-efficient analytical models are deemed necessary for depicting the fundamental bounds and trade-offs. Tools from applied probability theory, in particular stochastic geometry and point processes, are well-suited for characterizing the key performance metrics of DenseNets with random topologies, see [9], [10], and [11].
Despite insightful efforts, however, the common set of assumptions for studying cellular networks using stochastic geometry theory are benign, since rigorous analysis based on a direct signal-to-interference-plus-noise ratio (SINR) probability density function (pdf) approach is challenging [12], [13]. Most existing works thus resort to a limiting Rayleigh fading channel model which makes it possible to derive the SINR pdf [14], [15]. In [12], the authors utilize the non-direct moment-generating-function (MGF) methodology from [16] to characterize the average rate of always-full-buffer multi-tier cellular networks over arbitrary fading interference channels. However, neglecting network load by assuming that every deployed BS is always-transmitting leads to an unrealistic fullyloaded interference field which severely limits the achievable gains jointly in terms of throughput, deployment cost, and energy efficiency [17], [18], [19].
In [20], the authors provided a comprehensive stochastic geometry-based framework for the study of coverage probability and energy efficiency in HetNets using a full traffic load model with random BS switching. The authors in [17], however, showed that a significant gain in coverage probability performance can be achieved through conditional thinning of the interference field as a function of the users' density. Recently, in [19], a generalization of the interference-thinning approach for the modeling of load with multiple resource blocks was proposed. Moreover, the work in [21], with a focus on energy-efficient design, investigated the optimal deployment density in load-dependent two-tier cellular networks by numerically fitting the Poisson-Voronoi cell sizes using the Gamma function. A similar modeling approach based on the Gamma approximation of the random cell sizes was in addition adopted for the analysis of energy efficiency in HetNets in [22], and for the study of the energy-rate region in singletier cellular deployments in [23], respectively. Besides the lack of a formulation for the exact distribution of the cell sizes, the impact of SNR operating regions, artificial-bias, and non-Rayleigh fading on the achievable performance was not investigated in [21]. Furthermore, the dynamic switching of nodes (sleep mode mechanism) utilized in the existing works such as [20], [21], and [22] is limited due to the global coverage constraint associated with the inflexible traditional cellular network architecture in which the BSs must simultaneously provide coverage and capacity.
Recently, a major design overhaul in the cellular network architecture with a complete separation of the control and data infrastructures has been proposed through the Beyond Green Cellular Green Generation (BCG 2 ) project within the Green-Touch consortium [24], [25]. The control (signaling) network is responsible for providing continuous global coverage so that communication services can be requested by users at any given time and location. A sparse overlay of large signaling-only cells with extended range is considered to be the preferred solution in terms of network deployment and energy expenditure costs for this functionality. The data (capacity) network, on the other hand, is in charge of delivering communication services by intelligently activating resources on-demand based on an optimal selection of the access device which can meet the service requirements at minimum energy cost. The network layout for the capacity plane can in general be heterogeneous but pinpointing the optimal deployment solution is challenging and largely remains an open problem.
In [18], we incorporated the notion of inherent spatialcorrelations and load-proportionality in the design and analysis of multi-tier cellular networks. A computationally-efficient stochastic geometry-based framework for calculating the average rate of a typical user in the HetNet paradigm was accordingly provided. It was shown that the average rate performance of realistic correlated load-aware HetNets is significantly more optimistic that the state-of-the-art results which consider BSs to be either always or independently transmitting. In this paper, we extend the work in [18] to identify the most energyefficient combination of BS densities and the corresponding power savings through dynamic BS switching in multi-tier cellular networks. It should be noted that stochastic geometry Coverage Plane Capacity Plane Fig. 1: Green cellular architecture based on the separation of the control and data networks. theory cannot provide insight into the precise locations of the network nodes. The proposed methodology can however give us valuable information concerning the optimal type and number of BSs needed to satisfy certain user requirements. This approach can also be used to identify how many BSs can be switched off for the purpose of energy savings under the fluctuations in traffic volume.
Here, we focus on the green cellular network architecture in Fig. 1, where global coverage is provided by sparse large signaling-only cells and capacity is injected on-demand using dense data-only BSs. This allows for greater flexibility in utilizing sleep modes in the capacity plane which is no longer constrained by the global coverage constraint. Considering randomly-deployed cellular networks, we incorporate the notion of load-proportionality and correlated interfering sources by optimally and exclusively associating every user equipment (UE) to a data-only BS which provides the greatest reward under arbitrary shadowing characteristics. Closed-form expressions for the statistics of the received signal power and aggregate network interference over Nakagami-m fading channels are accordingly provided towards efficient computation of the average rate. An optimization problem for computing the optimal deployment solution that minimizes the total energy expenditure whilst satisfying a minimum rate requirement under a given network load is hence formulated and tackled using exhaustive search algorithms. For the special case of homogeneous interference-limited DenseNets, we provide new closed-form bounded solutions of the average rate and optimal BS density. Strategic sleep modes are then utilized for realizing power savings according to the temporal fluctuations in the traffic volume. The validity of the proposed analytical framework and its advantages in terms of preserving energy over the state-of-the-art fully-loaded and interference-thinningbased models are demonstrated via Monte-Carlo trials.
Several useful design guidelines are concluded from our findings. In general, we illustrate that the minimum deployment density required to satisfy the traffic requirements is significantly smaller in realistic load-proportional cellular networks compared to most existing results which typically assume independence among the BS activities. Furthermore, we show that on the contrary to the fully-loaded and interferencethinning approaches, our analytical model closely matches the actual optimal deployment density. The implications of these trends on the overall energy consumption and efficiency of the radio access network are accordingly depicted. In addition, the results confirm the promising potential of network densification towards effective offloading of traffic from largecells onto small-cells. The artificial-expansion of the smallcells coverage range, on the other hand, is shown to only further improve the performance when experiencing lower network loads. In addition, we demonstrate the limitations of large-cells in interference-dominant operating regions; a similar trend is also observed for small-cells in noise-dominant regions. Under anticipated traffic and rate requirements for dense urban cellular environments in the year 2020, the optimal deployment solution for the capacity plane is calculated to be a homogeneous pico network, capable of delivering peak power savings of near 15 kW/km 2 over a conventional macroonly cellular data network.
The remainder of the paper is organized as follows. The system model and mathematical preliminaries are described in Section II. In Section III, a computationally-efficient analytical framework for the evaluation of the average rate is provided. In Section IV, an optimization problem for identifying the optimal deployment solution is formulated and the corresponding power savings with dynamic switching of BSs is discussed. In Section V, theoretical and simulation studies are conducted towards unveiling network design pointers. Finally, the paper is concluded in Section VI.
Notation: E x {.} denotes the expectation operator with respect to random variable x; P(x) is the probability of event x; P x (.) represents the pdf of random variable x; M x (z) = E x {exp(−zx)} is the MGF of random variable x; |.| is the modulus operator; . corresponds to the Euclidean distance; Γ(x) =  with spatial densities λ (u) and λ We consider a co-channel deployment with universal frequency reuse where each operating BS equally allocates resources in terms of time or frequency slots to its associated UEs [26], [27]. This implies that there is no interference from transmissions associated with the same BS. For the sake of analytical tractability, we assume all tier-t BSs to have equal transmit power P tx t , artificial-biasing weight β t , and path-loss intensity α t . Let Y t,l,k , H t,l,k , χ t,l,k , and P rx t,l,k denote the Euclidean distance, fading power gain, shadowing power gain, and received signal power at the k-th UE from the l-th tier-t BS, respectively. In addition, the constant additive noise power is denoted by η. Note that the framework can be extended to the case where all nodes are equipped with multiple antennas, see, e.g., [28], [29].
In this work, we employ a cellular association and loadbalancing strategy where every active user exclusively connects to the closest BS of a certain tier which provides the strongest shadowed received signal power. This can be mathematically formulated as where ϕ t,l,k is a binary decision variable depicting whether or not the k-th UE is served by the l-th tier-t BS and constraints in (2) and (3) ensure that each UE is exclusively associated with exactly one BS. Accordingly, the optimal binary decision variables of all UEs are selected. The corresponding SINR of the k-th UE served by the l * -th tier-t * can hence be expressed as where and Based on the results from [30] and considering identical distribution across links, shadowing effects can be interpreted as random displacements in the original BSs locations using new transformed PPPs Φ iff E χ 2/αt t < ∞, ∀t ∈ T . As an application example, we consider Log-Normal shadowing with mean µ t (dB) and standard deviation σ t (dB), where t ∈ T . Note that small scale fading does not impact cell selection as it can be averaged or equalized using narrowband partitioning schemes such as orthogonal-frequency division multiplexing (OFDM). In this work, we consider independent unit-mean Nakagami-m fading for intended and interfering links. In this case, the pdf and MGF of the normalized channel power gain between the l-th tier-t BS and k-th UE are respectively expressed as [31] and where m t , t ∈ T , is the Nakagami-m fading parameter which can fit a wide-range of stochastic fading models. The ideal energy consumption behavior of a cellular network is load-proportional where the whole system power usage varies linearly according to the network load, i.e., from

III. AVERAGE RATE PERFORMANCE
In this section, we provide a framework for calculating the average communication rate achievable by an arbitrary user in the DenseNet paradigm. The Shannon channel capacity formula, i.e., log 2 (1 + SINR) b/s/Hz, is applicable here assuming capacity-achieving codes are used for the operating instantaneous SINR. Note that the model can be easily adjusted to capture other modulation/coding schemes by adding a SINR gap to the instantaneous rate formula, i.e., log 2 1 + SINR The average rate in b/s/Hz of an arbitrary UE k assumed to be located at the origin can be mathematically formulated by where ϕ t * ,l * ,k is used to denote the probability of UE k being tagged to the l * -th tier-t * BS, and E{log 2 (1 + γ t * ,l * ,k )} is the average rate of UE k conditioned on its association to the l * -th tier-t * BS, respectively.
In [16], Hamdi showed that the capacity evaluation of wireless communication systems can be greatly simplified by expressing the averages like E log 2 in terms of the MGFs of the independent random variables X and Y , i.e., Through extending this result to a stochastic geometry-based settings, the average rate expression in (10) can be expressed as in (11) where P Ŷ t * ,l * ,k (R) is the pdf of the random distance, and M P rx t * ,l * ,k |R and M I agg,k |R denote the conditional MGFs of the intended signal power and aggregate network interference, respectively. The pdf of transmitter-receiver distance and tier connection probability can be respectively calculated through the analytical expressions in (12) and (13) [18] where (a) follows from cases with equivalent path-loss exponent across all different tiers. The MGF of the intended signal power over Nakagami-m fading channels can be readily computed using the following expression Furthermore, we can derive a closed-form bounded expression for the aggregate network interference MGF -considering the inherent spatial-correlations in the activities of loadproportional BSs -as in the following lemma.
Lemma 1. The aggregate network interference MGF in spatially-correlated load-proportional heterogeneous DenseNets over Nakagami-m fading interference channels is given bỹ and In the case of Rayleigh fading interference channels, A t reduces to For the special case of α t = 4, ∀t ∈ T , the above can be further simplified to Proof: See Appendix A.
It should be highlighted that adopting the proposed generalized analytical framework allows for the efficient computation of average rate boundR(.) in spatially-correlated loadproportional multi-tier cellular networks over Nakagami-m fading channels through double integral operations. The com-mon direct pdf-based approach, on the other hand, involves manifold integral computations and is hence significantly more resource-intensive. For homogeneous DenseNet deployments, the average rate boundR(.) expression can be reduced to a single-integral format when considering interference-limited Rayleigh fading channels, as illustrated in the following lemma. Note that the tier index t is accordingly removed from the system parameters for single-tier deployment scenarios whenever the context is clear.
The average rate bound for interference-limited homogeneous DenseNets over Rayleigh fading channels can be expressed as where csc 2π The above integral can be further simplified for the special case of path-loss exponent being equal to four asR or alternativelyR Proof: See Appendix B.

IV. OPTIMAL DEPLOYMENT SOLUTION
From the mobile operators point of view, the commercial viability of network densification depends on the underlying capital and operational expenditure [32]. While the former cost may be covered by taking up a high volume of customers, with the rapid rise in the price of energy, and given that BSs are particularly power-hungry, energy efficiency has become an increasingly crucial factor for the success of DenseNets [33]. Generally, there are two main approaches to enhance the energy consumption of cellular networks: (1) improvement in hardware and (2) green system architecture design. Evidently, improving the power consumption of hardware is important. The potential gain, however, will be similar to the business-asusual case. Hence, the majority of improvement would have to come from green cellular network design.
Due to the recent advances in hardware technology, it has been made possible for wireless transceivers to consume varying power levels under different operational modes [34]. These include BS sleep, idle, transmit, and receive modes which can be accordingly adjusted based on the daily fluctuations in the traffic volume for the purpose of preserving energy. Defining a quantitative BS power model is however challenging given that one needs to take into consideration the particular components configurations. The following linear power model is however shown to be a reasonable approximation [35] where for tier-t BSs, ∆ (p) t is the reciprocal of the power amplifier drain efficiency, P (e) t is the circuit power, and P (s) t is the power in sleep mode. This linear power model accounts for the different specifications and architectures of long-termevolution (LTE) BSs including macro, micro, and pico types. The sleep mode power consumption (when there is nothing to transmit) is also included in this model to reflect upon a promising energy savings mechanism associated with future BSs. A set of power values for a default operating scenario can be found in [35,TABLE 2]. Note that the inclusion of backhauling constraints complicates the analysis here and is therefore left for future work. The reader is referred to [36] for a comprehensive study on modeling and tradeoffs of backhauling in HetNets using stochastic geometry.
In the remaining parts of this correspondence, we utilize the proposed stochastic geometric framework in order to pinpoint the most energy-efficient deployment solution and hence analyze the achievable gain in energy efficiency by incorporating dynamic sleep modes. The radio planning task under consideration is concerning a service provider interested in figuring out the most energy-efficient deployment solution for providing a minimum average rate to its users. Note that we assume global coverage is maintained through deploying, or utilizing the already in place, legacy macro-cells; the reader is referred to [32] for information on the operational characteristics and energy savings procedures in the separated coverage plane.
In order to compute the optimal combination of BS densities towards minimizing total energy expenditure, we formulate the following optimization problem min.
where R 0 is the minimum rate requirement and C t is the energy expenditure associated with tier-t BSs in active (transmission) mode. Note that R is a strictly increasing monotonic function in λ (a proof is provided in Appendix C), hence, the optimization problem under consideration has a unique solution. However, R is a highly complex function which involves for each tier two improper integrals and an infinite series sum. As a result, an exact closed-form solution cannot be obtained. The optimization problem under consideration should therefore be tackled numerically.
Exhaustive search algorithms are well-suited for tackling the problem considering the rate function derivative is not available analytically and its accurate evaluation (e.g., using finite differences) is resource-intensive. In the case of homo-geneous DenseNets, the dimension of freedom is reduced to one and the optimization task is equivalent to finding the root of a univariate function. Hence, Brent's algorithm [37] is a reasonable method of choice which at its best (worst) provides super-linear (linear) convergence to the solution. On the other hand, the non-linear constrained multidimensional optimization problem in the case of heterogeneous DenseNets can be tackled using heuristic downhill simplex method [38] with penalty function. The algorithm is typically solvable in exponential time [39]. The reader is referred to [40] for detailed descriptions of the above algorithms operation and code in C language.
In order to gain an analytical insight into the effect of different operational settings on the most energy-efficient deployment solution, we focus on the problem of finding the optimal BS density in homogeneous DenseNets. Specifically, we show that it is possible to develop lower-bound and upper-bound closed-form expressions of the average rate for interference-limited cases with the path-loss exponent being equal to four.

Lemma 3. The lower-bound and upper-bound closed-form expressions of the average rate in interference-limited homogeneous
DenseNets over Rayleigh fading channels with path-loss exponent being equal to four are respectively given bỹ Proof: See Appendix D.
For the above case, we derive closed-form solutions for the minimum number of BSs per unit area λ (b) * needed to satisfy the service requirement based on the numerical roots of high-order polynomial functions. The results, capturing the worst-and best-case scenario of the deployment solution, are respectively presented in the following lemma.
Lemma 4. The bounded solutions of the optimal BS density in interference-limited homogeneous DenseNets over Rayleigh fading channels with path-loss exponent being equal to four are expressed as and Proof: See Appendix E.
It is important to note from the best and worst BS deployment density expressions, respectively derived in (29) and (30), that the optimal network density λ (b) * is directly proportional to the network load λ (u) . The tightness of the developed bounded expressions is analyzed and compared to theoretical (using exhaustive search algorithms) and Monte-Carlo simulation results in the next section. Note that the analytical tractability of the techniques used to derive these bounds quickly diminishes considering multi-tier deployments. For heterogeneous DenseNets with equivalent operational parameters across different tiers, i.e., equal energy cost, biasing weight, and transmission power, however, we respectively arrive at similar bounded optimal network density expressions with In order to estimate the optimal energy savings from loadproportional network behavior in a given cellular environment, we utilize the established theoretical framework to identify the number and type of BSs that maintain the rate requirement for the users' density over different times of the day. The power savings in Watts using dynamic sleep modes at a given time of the day can then be computed by utilizing the linear hardware model as follows where λ , t ∈ T , are the optimal tier-t network densities under full and partial (depending on the hour) network loads, respectively. It is important to note that although this framework cannot determine an optimal topology for a given area, it can provide useful information on how many BSs can be switched off with the temporal variations in the traffic volume.

V. PERFORMANCE EVALUATION
The aim of this section is to evaluate the average rate, optimal deployment density, and power savings of DenseNets, considering different combinations of large-cell macro and small-cell micro and pico BSs. We aim to quantify the impact of network-wide decisions which helps unveil important design pointers for optimal network management. In regards to the power model, we use the practical hardware values captured in [35]. To analyze the accuracy of the established theoretical model, we perform load-dependent Monte-Carlo simulations (see Appendix F).

A. Framework Validation and Impact of System Parameters
The performances of a mixed micro/pico system under different SNR and load levels using the interference-thinning model, proposed green framework, and Monte-Carlo simulations are shown in Fig. 4. A key point to observe is that the former approach, while being an improvement over the long-standing fully-loaded model, produces pessimistic performance values, particularly under light and moderate network loads. Furthermore, our proposed analytical model correctly provides a tight lower-bound fit of the actual performance curve. The gap between different evaluation tools is negligible under both heavy traffic, due to the full-loaded interference field, and low SNRs, due to noise power dominance over interference. To further illustrate the shortcomings of the uncorrelated interferers assumption, we plot the transmission probability of the micro and pico tiers in a heterogeneous DenseNet as a function of network load in Fig. 5.
The implication of the above trends on network cost is depicted in Fig. 6, where we calculate the optimal network density of an interference-limited macro-only system as a function of minimum rate demand using the interferencethinning model, the proposed green framework, Monte-Carlo trials, and the bounded approximations derived in eqs. (29) and (30). The applicability of our proposed bounded framework in capturing realistic scenarios is further confirmed as it provides a tight fit to the Monte-Carlo trials. The interference-thinning model, on the other hand, requires a much larger BS density to meet a particular rate requirement. E.g., from Fig. 6, to satisfy R 0 = 2.4 b/s/Hz, the interference-thinning model requires 10.7%, 82.0%, and 159.0% larger network density over the closed-form upper-bound solution, proposed green framework (with exhaustive search algorithms), and Monte-Carlo trials, respectively. Henceforth, the analysis is carried out using the proposed framework as we have extensively shown its advantages over the state-of-the-art models.
Next, we depict the impact of network densification using small-cells with different biasing values on average rate performance in Fig. 7. Firstly, we observe that performance improves nearly linearly by adding unbiased pico BSs, further confirming the promising potential of small-cells in offloading traffic from large-cells in congested areas. However, deploying small-cells with low artificial-bias deteriorates performance. The reason lies on the added intra-tier interference experienced by pico BSs without a significant reduction in the inter-tier interference from the micro BSs. We can thus infer that the performance gain from artificial expansion of small-cells range is negative or otherwise negligible under relative moderate and heavy traffic loads.

B. Optimal Deployment Solution
Here, we investigate the optimal deployment solution under traffic and rate requirements anticipated for the year 2020. Based on the study of mature markets within GreenTouch [41] the peak traffic volume in the busy hour of a dense urban region in 2020 is 702 Mbits/sec/km 2 . Further, a discrete daily  traffic profile ranging from 20% to 140% of the average load in dense urban environments is utilized [41], see Fig. 8. Using the above parameter values, the active UEs density at 100% load level is calculated to be 84.87 UEs/km 2 . The recent Ofcom Market Report states that the average rate requirement for users on fourth-generation (4G) services in 2020 is expected to reach 2.4 b/s/Hz. Moreover, typical propagation values in dense urban environments are selected with Nakagamim fading with m = 2 (i.e., two-antenna transmit diversity Rayleigh), Log-Normal shadowing with µ t = 0 and variance σ 2 t = 6, and path-loss exponent α t = 4. The optimal BS densities, which meet users demand under different relative load levels in high and low SNR operating regions, of two different single-tier micro and pico DenseNets are depicted in Fig. 9. The optimal network density is shown to vary significantly in different hours, e.g., for the noisedominant case of the micro-only DenseNet, the optimal BS density at 20% load corresponding to 4-6 am (morning time) is 90.14% less than the value at 140% load experienced between 4-10 pm (night time). Furthermore, because of the negligible impact of transmit power in the interference-dominant cases, the optimal deployment density in the two different micro-only and pico-only systems are almost the same. This trend highlights the major disadvantage of deploying energyhungry large-cells in dense interference-dominant scenarios of the future. On the other hand, in noise-dominant regions, a greater number of pico BSs is required to meet the service requirements, e.g., at 20% load, λ (b) * of the pico-only system is 2.94 times greater than the optimal BS density of the singletier micro deployment.
We now turn to the more challenging problem of optimal deployment solution in multi-tier DenseNets. Considering all combinations of macro, micro, and pico BSs, we employ an exhaustive search algorithm to compute the most energyefficient deployment solution, e.g., the ratio of the energy cost in transmission mode of a micro BS over a pico BS is Cm Cp = 9.8085 and a macro BS over a pico BS is C M Cp = 30.7514. Our findings interestingly reveal that a small-cell deployment with only pico BSs with the values previously provided in Fig. 9 is the optimal solution for minimizing total energy expenditure. For comparison purposes, we identify the second most energyefficient deployment solution as a heterogeneous sparse micro and dense pico network with the values provided  b/s/Hz over the optimal homogeneous pico DenseNet. Note that by further reducing the noise power down to zero watts, we record no significant changes to the optimal type and number of BSs obtained.
Note that the practical feasibility of the above solutions is justified given the already in place legacy cellular networks are utilized for satisfying the global coverage constraint. However, a "practically optimal" solution could arguably accommodate the existence of legacy macro-cells in both data and control planes. In Fig. 10, we depict the optimal pico BSs density needed to satisfy the rate requirement in mixture legacy cellular networks overlaid with different amounts of macrocells. We observe that the use of high-power macro-cells in the data plane for dense urban environments has a detrimental impact considering a higher amount of small-cells is needed to meet the average rate constraint. For example, under a relative maximum load of 140%, the optimal deployment density of small-cells is around 3.43% lower in a homogeneous pico system over a mix macro/pico DenseNet with λ (b) M = 0.1 BSs/km 2 . It can therefore be inferred that in such environments the use of existing legacy cellular networks must be directed solely towards providing global coverage.

C. Power Savings via Sleep Modes
Finally, we analyze and compare the total power consumption and energy savings gain of different deployment scenarios in the dense urban environment under consideration. with and without BS sleep modes. The results confirm the optimality of the homogeneous pico DenseNet in terms of total energy efficiency, e.g., at peak traffic level, the mixed micro/pico deployment consumes more than 4.11% power than the optimal solution. Both of these solutions, however, are considerably more energy-efficient compared to any other combinations of BS densities such as the mix macro/pico system. E.g., the optimal pico DenseNet is capable of realizing peak power savings of near 15 kW/km 2 over a stand-alone macro-cell deployment; additional 20 W/km 2 and 40 W/km 2 improvements compared to the heterogeneous micro/pico and macro/pico DenseNets, respectively. It can also be seen from Fig. 11 that due to the large difference in the power usage of idle and sleep states, a significant reduction in total energy consumption can be achieved by activating BSs on-demand. Specifically, as shown in Fig. 12, by powering down BSs, relative to operating under full-load, peak energy efficiency gains of 35.36%, 36.57%, and 33.0% at night time, and average daily gains of 11.78%, 12.19%, and 10.97% for the pico-only, mixed micro/pico, and mixed macro/pico systems are respectively recorded.

VI. CONCLUSIONS
We have provided a comprehensive theoretical framework for performance evaluation and optimization of dense cellular networks. By incorporating the notion of load-proportionality in a flexible separated data/control plane cellular network architecture, we identify the most energy-efficient deployment solution for satisfying a minimum rate requirement under a given traffic level. The validity of the proposed green framework and its advantages over state-of-the-art fully-loaded and interference-thinning models in terms of pinpointing the optimal deployment density required for meeting users' demands was confirmed through extensive Monte-Carlo trials. Under a dense urban environment in the year 2020, the optimal deployment solution for the capacity plane was found to be a populated homogeneous pico network capable of realizing power savings of up to near 15 kW/km 2 compared to a traditional stand-alone macro cellular network.

APPENDIX A AGGREGATE NETWORK INTERFERENCE STATISTICS
The MGF of the aggregate network interference, considering a disk of radius D around the reference user and then taking the limit as D → ∞, can be derived as in (A.1) where (a) follows from applying Jensen's inequality which leads to is from the independence of the tiers of BSs with N (b) t being the total number of potentially interfering tier-t sources and j being an arbitrary tier-t source, (c) is from using a Binomial distribution with parameters (κ t , ρ t ) to characterize N (b) t , using the uniformly-distributed locations of the interfering sources being the distance of closest tier-t source, and (d) can be derived by taking the limits as D → +∞, κ t → +∞, ρ t → 0, and utilizing the Poisson limit theorem with t,s , finally, (e) is obtained by taking the average with respect to the Gamma-distributed fading power gain of the arbitrary interferer using

APPENDIX B SIMPLIFIED AVERAGE RATE EXPRESSION
Utilizing (14) and (15) in (11), considering η = 0 and m = 1, the bounded average rate expression in single-tier scenarios can be written in a double integral form as in (B.1). By employing a change of variables with z = u α ϕP tx and hence converting from Cartesian to polar coordinates with R = r sin(t) and u = r cos(t), (B.1) reduces to (B.2) By using the following integral identity (where x > 0) we arrive at the simplified average rate expression in (20). For the special case of α = 4, with variable substitution s = tan 2 (t) and some basic algebraic manipulations, (20) can be further simplified to obtain (21).

APPENDIX C MONOTONICITY ANALYSIS OF THE RATE FUNCTION
Without loss of generality, consider a homogeneous DenseNet with Rayleigh fading for the intended and interference channels and path-loss exponent being equal to four. Recall that the average rate in nat/s/Hz of an arbitrary user in this case can be expressed by To investigate the behavior of the rate function with respect to the deployment density, we differentiate using basic substitution the inside of the above integral with respect to λ (b) as where (C.2) follows given λ (u) only takes on positive values and s − arctan(s) + π 2 > 0 holds for positive values of s. This proves that the average rate is a strictly monotone function for BS deployment density. By performing a partial fraction decomposition of the above expressions we respectively obtaiñ R(.) = log 2 (e) 1 1 + ϕ π π 4 (ϕ + 1) − 1 × +∞ 0 2s − πϕ + 2 π 2 (s + 1) + s 2 + ϕ(π(ϕ + 1) − 2(s + 1)) 1 + ϕs 2 ds. t \{l * } E H t,l,k , Ŷ t,l,k e −zϕ t P tx t H t,l,k Ŷ t,l,k t E H t,j,k , Ŷ t,j,k e −zϕ t P tx t H t,j,k Ŷ t,j,k ρ t E H t,j,k , Ŷ t,j,k e −zϕ t P tx t H t,j,k Ŷ t,j,k From (D.6), (D.7), and using the rules of integration by parts, we arrive at the closed-form bounded expressions of the average rate in (27) and (28).

APPENDIX E OPTIMAL DEPLOYMENT DENSITY
The closed-form bounded expressions of the average rate in (27) and (28) are complex highly non-linear functions of the ratio of the BS over the UE spatial densities. As a result, it is not possible to directly derive an expression for the optimal deployment solution. By taking the limits as λ (u) → 0 (sparse traffic) and λ (u) → +∞ (full traffic), it can be readily deduced that 0 ≤ ϕ = 1 − λ (u) λ (b) ≤ 1. Hence, we can apply the lower-bound approximation ln(1 + ϕ) ≤ ϕ based on Taylor series expansion in (27) and (28), in order to rearrange R λ (u) , λ (b) , P tx , η, α, m, µ, σ = R 0 and develop upperbound and lower-bound closed-form approximations for the optimal network density λ (u) * based on the real positive real roots of high-order polynomial functions in (29) and (30).
APPENDIX F MONTE-CARLO SIMULATIONS 1) Set the UEs density, and for each tier, select transmit power, BSs density, path-loss exponent, Nakagami-m fading, and Log-Normal shadowing mean and variance. 2) Define a region of sufficiently large area around reference UE situated at the origin. 3) Generate the statistical numbers of tiers of BSs and UEs.

4) Deploy Uniformly-distributed heterogeneous BSs and
UEs around the specified area. 5) Generate Nakagami-m fading and Log-Normal shadowing gains for all links. 6) Optimally and exclusively associate every UE to a BS which provides the strongest received shadowed power. 7) Search through all BSs and if a BS is associated with one or more UEs it is active; otherwise is not transmitting. 8) Compute the aggregate network interference experienced by the reference UE using the sum of received signal powers from only the interfering BSs. 9) Calculate the reference UE SINR and average rate. 10) Repeat steps (3) to (9) for a sufficiently large number of times and take the average.