---
abstract: 'How would a cellular network designed for maximal energy efficiency look like? To answer this fundamental question, tools from stochastic geometry are used in this paper to model future cellular networks and obtain a new lower bound on the average uplink spectral efficiency. This enables us to formulate a tractable uplink energy efficiency (EE) maximization problem and solve it analytically with respect to the density of base stations (BSs), the transmit power levels, the number of BS antennas and users per cell, and the pilot reuse factor. The closed-form expressions obtained from this general EE maximization framework provide valuable insights on the interplay between the optimization variables, hardware characteristics, and propagation environment. Small cells are proved to give high EE, but the EE improvement saturates quickly with the BS density. Interestingly, the maximal EE is achieved by also equipping the BSs with multiple antennas and operate in a “massive MIMO” fashion, where the array gain from coherent detection mitigates interference and the multiplexing of many users reduces the energy cost per user.'
author:
- |
    Emil Bj[ö]{}rnson, *Member, IEEE*, Luca Sanguinetti, *Senior Member, IEEE*,\
    and Marios Kountouris, *Senior Member, IEEE* [^1]
bibliography:
- 'IEEEabrv.bib'
- 'refs.bib'
title: 'Deploying Dense Networks for Maximal Energy Efficiency: Small Cells Meet Massive MIMO'
---

Introduction {#sec:intro}
============

The biggest challenges for next generation wireless communication systems (5G) are to support the ever-growing demands for higher data rates and to ensure a consistent quality of service (QoS) throughout the entire network [@Andrews2014a]. To meet these demands, the network area throughput (in bit/s/km$^2$) needs to increase by a [factor of 1000]{} over the next 10–15 years [@Qualcomm]. At the same time, the power consumption of the information and communication technology (ICT) industry and the corresponding energy-related pollution are becoming major societal and economical concerns [@Fehske2011a]. Credited sources foresee that, to meet such a $1000 \times$ higher data traffic without increasing the ICT footprint, new technologies that improve the overall energy efficiency (EE) by $1000 \times$ must be developed [@Greentouch2013]. Hence, higher network area throughput on the one hand and less power consumption on the other are seemingly contradictory 5G requirements [@Bjornson2014c]. There is a broad consensus that these goals can only be jointly achieved by a substantial network densification.

Two promising technologies towards network densification are small-cell networks [@Hoydis2011c; @Sanguinetti2015bJSAC] and massive MIMO systems [@Marzetta2010a; @Larsson2014a; @Yang2013b]. The first technology relies on an ultra-dense and irregular operator-deployment of low-cost and low-power base stations (BSs), with higher density where the user load is higher. Bringing the BSs and user equipments (UEs) closer to each other can increase the area throughput, while significantly reducing the radiated signal power. However, this may come at the price of a substantial increase of the circuit power consumption (per km$^2$) due to the larger amount of hardware and infrastructure [@Bjornson2013e]. In contrast, the massive MIMO technology aims at evolving the conventional BSs by using arrays with a hundred or more small dipole antennas. This allows for coherent multi-user MIMO transmission where tens of users can be multiplexed in both the uplink (UL) and the downlink (DL) of each cell. It is worth observing that, contrary to what the name “massive” suggests, massive MIMO arrays are rather compact; 160 dual-polarized antennas at 3.7 GHz fit into the form factor of a flat-screen television [@Vieira2014a]. In massive MIMO systems, the area throughput is improved by the multiplexing gain, while the array gain from coherent processing allows for major reductions in the emitted power. Similar to small-cell networks, however, the potential throughput gains from massive MIMO come from deploying more hardware (i.e., multiple antenna branches per BS), which in turn increases the circuit power consumption per BS.

In summary, both densification technologies can improve the area throughput and reduce the radiated power, but at the cost of deploying more hardware infrastructure. Hence, the overall EE of the network can only be improved if these benefits and costs are properly balanced. The main objective of this paper is to develop an analytical framework for designing dense cellular networks for maximal EE, as well as to provide guidelines for practical deployment when (among others) the spectral efficiency (SE), radiated power, BS density, number of antennas and UEs per BS, channel estimation, and the circuit power consumption are taken into account from the outset.

Related Works
-------------

The EE of cellular networks has been defined from different perspectives in the last decade [@Feng2013a; @Li2011a]. One of the most common definitions is a benefit-cost ratio, where the service quality per area unit (a.u.) is compared with the associated energy costs. In this paper, the following general definition is used: $$\begin{aligned}
&\mathrm{EE} = \\ &\frac{\textrm{Area spectral efficiency [bit/symbol/km}^2]}{\textrm{Transmit power + Circuit power per a.u.~[J/symbol/km}^2]}\end{aligned}$$ which yields to an EE metric measured in bit/J.[^2]

Most of the early works that analyzed the EE, as defined above, have focused on the single-cell case wherein the interference from other cells is neglected. Examples in this context are [@Miao2013a; @Hu2014a; @Bjornson2013e; @Ha2013a; @Yang2013a; @Mohammed2014a; @Bjornson2014b; @Bjornson2015a]. In particular, in [@Miao2013a] the authors consider the UL power allocation of multi-user MIMO systems and show that the EE is maximized when specific UEs are switched off. The UL was studied also in [@Hu2014a], where the EE was shown to be a concave function of the number of BS antennas, $M$, and the UE rates. The DL has been investigated in [@Bjornson2013e; @Ha2013a; @Yang2013a], whereof [@Bjornson2013e] and [@Ha2013a] show that the EE is a concave function of $M$, while a similar result is shown for the number of UEs, $K$, in [@Yang2013a]. Unfortunately, in these works, the system parameters were optimized by means of simulations that (although useful) do not provide a complete picture of how the EE is affected by $M$ and $K$. Recently, [@Mohammed2014a] derived the optimal $M$ and $K$ for a given UL sum rate, while [@Bjornson2015a] derived the optimal $M$, $K$ and UE rates to maximize the EE jointly for the UL and DL. These two works utilize a realistic power consumption model where the dissipation in circuits depends non-linearly on the above parameters. Simple closed-form expressions for the EE-maximizing parameter values and their scaling behaviors are derived for zero-forcing (ZF) processing with perfect channel state information (CSI) and verified by simulations for other processing schemes under imperfect CSI.

The EE analysis of multi-cellular networks is much more involved than in the single-cell case due to the complicated network topology and the arising inter-cell interference. The simplest approach is to rely on heavy Monte Carlo simulations. Attempts in this direction can be found in [@Kurniawan2012] and [@Liu2013]. Unfortunately, Monte Carlo simulated results are often anecdotal since one cannot separate fundamental properties from behaviors induced by parameter selection. Alternatively, simplified network topologies can be considered, such as the Wyner model [@Wyner1994a] or the symmetric grid-based deployment [@Bjornson2015a]. While attractive for its analytical simplicity, the Wyner model does not well capture the essential characteristics of real and practical networks [@Xu2011a]. Similarly, the symmetric grid-based models cannot capture the irregular structure of small-cell network deployments.

The need for developing tractable (yet reasonably accurate) models for future dense networks has increased the interest in random spatial models; more particularly, in using tools from stochastic geometry [@Haenggi-2009], wherein the BS locations form a realization of a spatial point process—typically a Poisson point process (PPP). A major advantage of this approach is the ability of providing tractable expressions for key performance metrics such as the coverage probability and the average SE of the network. A few prior works have also derived EE-related performance metrics and showed how these depend on the BS and UE densities; for example, [@Cao2012a] compares the deployment of two types of single-antenna BSs, while [@Li2014a] studies the EE when multi-antenna BSs serve one UE each. In [@Soh2013], the authors focus on the analysis and design of energy-efficient heterogeneous cellular networks (HetNets) through the deployment of small cells and sleeping strategies. The effect of cross-tier BS cooperation and clustering on the EE in the DL is studied in [@Nie2014]. In contrast to the UL analysis in our paper, the vast majority of work on the EE in multi-cell networks using stochastic geometry has focused on the DL. Furthermore, prior work usually aims at deriving closed-form SE and EE expressions, which often turn out to be involved and cumbersome to optimize, thereby providing little or no useful insights on the optimal system design.

Major Contributions
-------------------

This paper considers the UL of a multi-cell multi-user MIMO network in which the BSs are distributed according to a homogenous PPP of intensity $\lambda$. Each BS is equipped with $M$ antennas and communicates with $K$ single-antenna UEs uniformly distributed within its coverage area. Coherent detection based on maximum ratio combining (MRC) is used at the BSs, based on CSI acquired from UL pilot signaling with each pilot symbol being reused in $1/\beta$ of the cells. In addition, we assume that the UEs are prone to practical transceiver hardware impairments and apply UL power control based on statistical channel inversion (with proportionality coefficient $\rho$), to resolve the near-far effects that otherwise would prevent communication [@Bjornson2016a]. Tools from stochastic geometry and classical statistics are used to compute a new lower bound of the average SE. This expression is then used to define the EE metric by also using the power consumption model developed in [@Bjornson2014b], which not only accounts for the radiated power but also for the operating power consumption required by (among others) analog transceiver chains, digital processing, and backhaul infrastructure.

Within the above conditions, an EE maximization problem is formulated under the assumption that a given average SE target per UE must be met with equality, to guarantee good service quality. This problem is solved analytically with respect to the tuple $\theta = (\beta,\lambda, \rho,M,K)$ of optimization variables. The resulting closed-form expressions provide valuable design insights on the interplay between the different system parameters and the various components of the power consumption model. We show analytically and numerically that the radiated power is negligible in most practical cases and that massive MIMO configured BSs appear naturally since this technology can protect the desired signals from inter-cell interference and share the circuit power between many UEs.

Outline and Notation
--------------------

The remainder of this paper is organized as follows. Section \[sec:system\] introduces the system model for the network under investigation, including hardware impairments at the UEs and CSI estimation at the BSs. A new and tractable lower bound on the average SE is derived in Section \[sec:ASE\]. It is used in Section \[sec:problem\] to formulate the EE maximization problem and to obtain closed-form expressions for the EE-maximizing pilot reuse factor, BS density, and transmission power as well as the number of UEs and antennas per BS. The theoretical analysis is corroborated by numerical results in Section \[sec:simulations\] while the major conclusions and implications of this paper are drawn in Section \[sec:conclusion\].

The following notation is used throughout the paper. The notation $\mathbb{E}\{ \cdot \}$ indicates the expectation with respect to a random variable, whereas $\| \cdot \|$ and $| \cdot |$ stand for the Euclidean norm and absolute value, respectively. We let ${\mathbf{I}}_{M}$ denote the $M \times M$ identity matrix, whereas we use $\mathcal{CN}(\cdot,\cdot)$ to denote a multi-variate circularly-symmetric complex Gaussian distribution. We use $ \mathbb{C}$, $ \mathbb{Z}_+$, and $\mathbb{R}$ to denote the sets of all complex-valued numbers, all positive integers, and all real-valued numbers, respectively. The Gamma function is denoted by $\Gamma(\cdot)$.

System Model {#sec:system}
============

![Illustration of one realization of BS positions from the homogeneous PPP $\Phi_{\lambda}$ and of UEs uniformly distributed in the corresponding Poisson-Voronoi cells.[]{data-label="figureVoronoi"}](figureCells){width="0.95\columnwidth"}

-3mm

-3mm

We consider the UL of a cellular network that is designed to serve a dense heterogeneous distribution of UEs. This is modeled using the stochastic geometry framework adopted in [@Andrews2011a] in which the BSs are distributed spatially in $\mathbb{R}^2$ according to a homogeneous PPP $\Phi_{\lambda}$ of intensity $\lambda$ (measured in BSs per $\textrm{km}^2$). More precisely, this means that in any area of size $A$ (measured in $\textrm{km}^2$), the number of BSs is a Poisson distributed stochastic variable with mean value $\lambda A$. The BSs are uniformly and independently distributed over the area. Each of them is equipped with an array of $M$ antennas and serves $K$ single-antenna UEs, which are selected at random from a (potentially) very large set of UEs in the cell since we treat $K$ as an optimization variable. We assume that each UE connects to its closest BS such that the coverage area of a BS is its Poisson-Voronoi cell; see Fig. \[figureVoronoi\] for an illustration. The UEs are assumed to be uniformly distributed in the Poisson-Voronoi cell of their serving BS. Note that the geographic locations of UEs and BSs are correlated under this model, such that small cells serve more UEs per $\textrm{km}^2$ than larger cells. We interpret this as a judicious network deployment where the BSs are matched to a heterogenous distribution of UEs.[^3]

The translation invariance of PPPs allows us to perform statistical performance analysis for a *typical UE* (located for example in the origin), which is statistically representative for any other UE in the network [@Baccelli2008a]. Assume that this typical UE has the arbitrary index $k$ and is connected to a BS that is called the *typical BS* and is denoted as $\mathrm{BS}_0 \in \Phi_{\lambda}$. Then, the following basic properties hold (see for example [@Weber2010a]):

\[lemma:distance-distribution\] The distance $d_{00k}$ from the typical UE to its serving BS is Rayleigh distributed as $d_{00k} \sim \mathrm{Rayleigh} \big( \frac{1}{\sqrt{2\pi \lambda}} \big)$. The BSs of the other cells form a homogenous PPP $\Psi_{\lambda} = \Phi_{\lambda} \setminus \{ \mathrm{BS}_0 \}$ in $\mathbb{R}^2 \setminus \{ {\mathbf{x}} \in \mathbb{R}^2 : \| {\mathbf{x}} \| < d_{00k} \}$.

Using the notation introduced in Lemma \[lemma:distance-distribution\], each interfering BS, $\mathrm{BS}_j \in \Psi_{\lambda}$, is situated at the geographical position determined by its spatial point $\mathrm{BS}_j \in \mathbb{R}^2$. Since in most cases we are interested only in the index of the interfering BS, the shorter notation $j \in \Psi_{\lambda}$ is used henceforth for $\mathrm{BS}_j \in \Psi_{\lambda}$.

Channel Model and Power-Control Policy
--------------------------------------

We consider a UL transmission protocol where the time-frequency resources are divided into blocks of $T_c$ seconds and $W_c$ Hz. This leaves room for a total number of $S = T_c W_c$ transmission symbols for pilot signaling and data. The channel response ${\mathbf{h}}_{ljk} \in \mathbb{C}^M$ between $\mathrm{BS}_l$ and UE $k$ in cell $j$ is modeled as block-fading such that ${\mathbf{h}}_{ljk}$ has a constant stochastic realization within a block and takes independent realizations across blocks. Since the number of BS antennas is treated as an optimization variable in this work, we need a channel model that is reasonable for both small and large arrays. For this purpose, we make use of a Rayleigh fading model as it well-matches non-line-of-sight measurements in both cases [@Gao2015a]. Hence, we consider Rayleigh fading channels with ${\mathbf{h}}_{ljk} \sim \mathcal{CN}({\mathbf{0}}, \omega^{-1} d_{ljk}^{-\alpha} {\mathbf{I}}_M)$, where $d_{ljk}$ is the propagation distance (measured in km) and $\alpha> 2$ is the pathloss exponent. The parameter $\omega$ accounts for the pathloss at a reference distance of 1 km and can thus be used to model distance-independent propagation losses (such as wall penetration).

Power control is an essential mechanism in the UL of any multi-user MIMO system; the finite dynamic range of analog-to-digital converters at the BSs otherwise creates near-far blockage where weak signals from distant UEs drown in stronger signals from nearby UEs. Statistical channel inversion is a convenient power-control policy [@Bjornson2016a], where UE $i$ in cell $j$ uses the transmit power $$\label{p_{ji}}
p_{ji} = \rho  \omega d_{jji}^{\alpha}$$ with $\rho \geq 0$ being a power-control coefficient to be designed.[^4] The average transmit power at UE $i$ in cell $j$ will then be $$\label{p_{ji}_average}  
\mathbb{E}\{p_{ji}\} = \rho  \omega \mathbb{E}\{d_{jji}^{\alpha}\} = \rho  \omega \frac{ \Gamma({\alpha}/{2} +1 ) }{(\pi \lambda)^{\alpha/2}}$$ by using the distance distribution from Lemma \[lemma:distance-distribution\]. A key feature of the statistical power-control policy is that it provides the same average effective channel gain of $\mathbb{E}\left\{{p_{ji} \|{\mathbf{h}}_{jji}\|^2} \right\} = M {\rho}$ to all UEs (irrespective of their locations) without requiring any rapid feedback mechanism to provide UEs with instantaneous CSI. Note that this does not mean that the signal-to-interference-plus-noise ratio (SINR) is the same for all UEs, since the interference varies in the network.

Transceiver Hardware Impairments
--------------------------------

Any practical wireless transceiver is affected by hardware impairments; for example, unavoidable clock drifts in local oscillators, finite-precision digital-to-analog converters, amplifier non-linearities, non-ideal analog filters, etc. Such impairments can be partially mitigated by compensation algorithms [@wenk2010mimo], but not fully removed. Although negligible at low signal-to-noise ratios (SNRs), hardware impairments characterize the maximal achievable SE at high SNRs [@Bjornson2013c]. The vast majority of prior works on the design of energy-efficient networks assumes ideal transceivers, although hardware impairments might have a fundamental impact also on the EE. Motivated by this fact, the distortion noise arising from hardware impairments is included in the system under investigation. We concentrate on UE impairments, which is expected to be the dominating effect since small BS arrays can use high-grade hardware and massive arrays are resilient to hardware impairments [@Bjornson2014a].

Denote by $s_{ji} \in \mathbb{C}$ an arbitrary symbol transmitted by UE $i$ in cell $j$ and assume that it is normalized to unit power (i.e., $\mathbb{E}\{ | s_{ji}|^2 \}=1$). Similar to [@wenk2010mimo; @Zhang2012a; @Bjornson2013c], the hardware impairments are modeled as a reduction of the desired signal power by a factor of $1-\epsilon^2$, with $0 \leq \epsilon < 1$, and by adding Gaussian distortion noise that carries the removed power. More precisely, $s_{ji} $ is replaced with $\sqrt{1-\epsilon^2}s_{ji} + e_{ji}$, where $e_{ji} \sim \mathcal{CN}(0,\epsilon^2)$ is the additive UL distortion noise caused at UE $i$ in cell $j$. In these circumstances, the received signal ${\mathbf{y}}_0 \in \mathbb{C}^{M}$ at $\mathrm{BS}_0$ takes the form $$\label{eq:system-model-received}
\begin{split}
{\mathbf{y}}_0 &= \sum_{i=1}^{K} \sqrt{p_{0i}} (\sqrt{1-\epsilon^2}s_{0i} + e_{0i}) {\mathbf{h}}_{00i}  \\ &+ \sum_{j \in \Psi_{\lambda}}
\sum_{i=1}^{K} \sqrt{p_{ji}} (\sqrt{1-\epsilon^2}s_{ji} + e_{ji})  {\mathbf{h}}_{0ji}  + {\mathbf{n}}_{0}
\end{split}$$ where ${\mathbf{n}}_{0} \sim \mathcal{CN}(0,\sigma^2 {\mathbf{I}}_M)$ models receiver noise with variance $\sigma^2$. The parameter $\epsilon$ is referred to as the *level of impairments* and is tightly connected to the error vector magnitude (EVM), which is a common quality measure of transceivers.[^5] The 3GPP LTE standards specify EVM requirements in the range $[0.05,0.175]$, where higher-order modulations are supported if the EVM is below the lower part of this range [@Holma2011a]. Note that the ideal hardware case is obtained by simply setting $\epsilon = 0$.

Channel Acquisition
-------------------

As it is known, coherent processing of the received signal ${\mathbf{y}}_0$ requires knowledge of the UL channel vectors ${\mathbf{h}}_{00i}$ for $i=1,\ldots,K$. For this purpose, we assume that $B$ out of the $S$ symbols in each UL block are used for pilot transmission. This means that in the whole network there exists $B$ orthogonal pilot symbols that are shared among the cells. We assume that each BS picks $K \le B$ pilot symbols uniformly at random in each block, to avoid cumbersome pilot coordination. This means that on average $K/B$ of the cells reuse any given pilot symbol. We call $\beta = \frac{B}{K} \ge 1$ the pilot reuse factor such that $\beta K= B\le S$. In this setting, the pilot symbol, sent by the typical UE and received at the typical BS, is interfered with by the subset of cells in which there is another UE (still indexed by $k$ without loss of generality) reusing the same pilot symbol. This is modeled through a binary stochastic variable $\chi_{0kj} \in \{0,1\}$, where $\chi_{0kj}=1$ means that UE $k$ in cell $j$ uses the same pilot as the typical UE and thus causes pilot contamination. This occurs with probability $K/B$. Similarly, $\chi_{0kj}=0$ means that there is currently no pilot contamination from cell $j$, and it occurs with probability $1-K/B=1-1/\beta$. Therefore, the typical UE $k$ transmits a single pilot symbol (e.g., $s_{0k}=1$) and $\mathrm{BS}_0$ receives the following signal: $$\label{eq:received-pilot}
\begin{split}
{\mathbf{z}}_{0k} &= \sqrt{p_{0k}} (\sqrt{1-\epsilon^2} + e_{0k}) {\mathbf{h}}_{00k} \\ &+ \sum_{j \in \Psi_{\lambda} } \chi_{0kj} \sqrt{p_{jk}} (\sqrt{1-\epsilon^2} + e_{jk}) {\mathbf{h}}_{0jk} + {\bf{n}}_{0}.
\end{split}$$ We compute the minimum mean-squared error (MMSE) estimate of ${\mathbf{h}}_{00k}$ from the observation ${\mathbf{z}}_{0k}$.

$$\begin{aligned}
 \tag{8}  \label{eq:effective-SINR}
\frac{ M(1-\epsilon^2)^2   }{  \!\!\!\!\Bigg( K + {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } + \frac{\sigma^2}{  \rho}  \Bigg) \!\!\Bigg(  1   + \frac{1}{\beta} {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{\rho} \Bigg) 
+ M (1-\epsilon^2) \Bigg( \frac{1}{\beta} {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}  \!\!
\left( \frac{ d_{jjk}^{\alpha} }{ d_{0jk}^{\alpha}} \right)^2  + \epsilon^2 \Bigg)  }\end{aligned}$$

$$\begin{aligned}
 \tag{10} \label{eq:average-SINR}
\frac{ M (1-\epsilon^2)^2  }{ 
\Big( K + \frac{\sigma^2}{  \rho}  \Big) \! \Big(  1   + \frac{ 2}{\beta( \alpha-2)}  + \frac{\sigma^2}{\rho} \Big) \! + \!
 \frac{ 2K}{\alpha-2} \! \left(  1  + \frac{\sigma^2}{\rho} \right) \! + \!
\frac{K}{\beta} \! \left( \frac{4}{(\alpha-2)^2 } +
\frac{1 }{ \alpha-1} \right) \!
+  M (1-\epsilon^2) \Big( \frac{ 1}{\beta( \alpha-1)} \!+\! \epsilon^2 \Big)  }\end{aligned}$$

The MMSE estimate of the typical UE’s channel to its serving $\mathrm{BS}_0$ is $$\label{eq:MMSE-estimator}
  \hat{{\mathbf{h}}}_{00k} = 
\frac{  \sqrt{ \frac{1-\epsilon^2}{\rho \omega d_{00k}^{\alpha}} }}{ 1   +  \sum_{j \in \Psi_{\lambda} }    \chi_{0kj}    \frac{d_{jjk}^{\alpha}}{ d_{0jk}^{\alpha} } + \frac{\sigma^2}{ \rho} } {\mathbf{z}}_{0k}.$$ The estimation error $\Delta{\mathbf{h}}_{00k} = {\mathbf{h}}_{00k}  - \hat{{\mathbf{h}}}_{00k}$ is distributed as $\Delta{\mathbf{h}}_{00k}  \sim\mathcal{CN}({\mathbf{0}},{\mathbf{C}}_{00k})$, where the estimation error covariance matrix ${\mathbf{C}}_{00k} \in \mathbb{C}^{M \times M}$ is given by $$\label{eq:MMSE-error-cov}
{\mathbf{C}}_{00k} = \frac{1}{\omega d_{00k}^{\alpha}}  \left(1 - \frac{  1-\epsilon^2 }{ 1   +  \sum_{j \in \Psi_{\lambda} }    \chi_{0kj}   \frac{d_{jjk}^{\alpha}}{ d_{0jk}^{\alpha} } + \frac{\sigma^2}{\rho} }  \right) {\mathbf{I}}_M.$$

This follows from standard results on MMSE estimation; see [@Kay1993a Chapter 15.8].

Average Spectral Efficiency {#sec:ASE}
===========================

We assume that MRC processing is used at the BSs for data recovery, since MRC is computationally efficient and performs well in both small cells (with $K=1$) and massive MIMO (with $M \gg K$). In particular, the symbol transmitted by the typical UE is detected at $\mathrm{BS}_0$ by correlating the received signal in with the MMSE estimate $\hat{{\mathbf{h}}}_{00k}$; that is, $r_{0k} = \nu_{00k}\hat{{\mathbf{h}}}_{00k}^{\Htran}  {\mathbf{y}}_0$ is the received signal after MRC processing with $\nu_{00k} \in \mathbb{C}$ being a scaling factor (see Appendix A). Unfortunately, the ergodic capacity for a network (such as the one under investigation) in which only imperfect CSI is available at the receivers and the inter-cell interference is modeled as a shot-noise process is not known yet [@Weber2010a]. Hence, in this work the SE expression is obtained using classical achievable lower bounds on the ergodic capacity that are particularly common in the massive MIMO literature [@Marzetta2010a; @Yang2013b; @Bjornson2016a]. The following result holds for any given realization of $\Psi_{\lambda}$ and the UE locations.

\[lemma:SINR-ideal\] With MRC processing and any given realization of $\Psi_{\lambda}$ and the UE locations, a lower bound on the ergodic capacity \[bit/symbol/user\] of the typical UE is $$\label{eq:ergodic-capacity}
\Big( 1 - \frac{\beta K}{S} \Big) \log_2 \! \left( 1 + \mathrm{SINR}_{0k}  \right)$$ where the pre-log factor $\big( 1 - \frac{\beta K}{S} \big) $ accounts for pilot overhead and $\mathrm{SINR}_{0k}$, in at the top of the page, is the effective SINR.

The proof is given in Appendix A and proceeds as follows. First, the lower-bounding approach from [@Jose2011b] and [@Medard2000a] is used to obtain the mutual information in the presence of imperfect CSI. Then, the statistical properties of the MMSE estimator and power-control policy are exploited to compute the average with respect to the channels in the presence of pilot contamination.

The expression in Lemma \[lemma:SINR-ideal\] holds for any $\beta\geq 1$ such that $\beta K \leq S$, since the pilot signals need to be contained in a coherence block. Observe also that $\beta K$ does not need to be an integer since an arbitrary $\beta K$ can be achieved by switching (for appropriate fractions of time) between the closest smaller integer, $\lfloor \beta K \rfloor$, and the closest larger integer, $\lceil \beta K \rceil$.

The average SE per UE can be obtained by taking the expectation of with respect to the PPP $\Psi_{\lambda}$ and the UE locations. This would require heavy numerical evaluations of integrals. As a key contribution of this paper, in the following we provide a tractable and tight lower bound on the average SE.

\[prop:average-SE\] If MRC is employed, a lower bound on the UL average SE \[bit/symbol/user\] is $$\underline {\mathrm{SE}} = \Big( 1 - \frac{\beta K}{S} \Big)\! \log_2 \! \left( 1 + \underline {\mathrm{SINR}} \right)\label{eq:average-SE}$$ where $\underline {\mathrm{SINR}} $ is given in at the top of the page.

The proof is given in Appendix B and is based on Jensen’s inequality, which leads to tractable expressions since only moments of the terms in need to be computed.

The tightness of the SE bound in Proposition \[prop:average-SE\] is demonstrated later on by means of simulations (see Fig. \[figureBSdensity\]). The numerator of $\underline{\mathrm{SINR}} $ in scales with $M$ due to the array gain from coherent processing, but the hardware impairments cause the multiplicative loss $(1-\epsilon^2)^2$. The last term in the denominator scales with $M$ and also with the sum of $\frac{1}{\beta (\alpha-1)}$ and $\epsilon^{2}$. The former term is called coherent pilot contamination since it accounts for interference received from UEs that use the same pilot symbol as the typical UE, while the latter is due to the distortion noise emitted from the typical UE itself. Many of the interference terms in increase with $K$ since having more UEs lead to more transmit power per cell. Note that is independent of the BS density $\lambda$, due to the power-control policy in . The average transmit power per UE in is, however, proportional to $\lambda^{-\alpha/2}$ so that less power is needed (to sustain a fixed SNR per UE) as the network becomes denser.

Problem Statement and Energy Efficiency Optimization {#sec:problem}
====================================================

As described in Section \[sec:intro\], we concentrate on the UL EE defined as the benefit-cost ratio between the area spectral efficiency (ASE) \[bit/symbol/$\textrm{km}^2$\] and the area power consumption (APC) \[J/symbol/$\textrm{km}^2$\]. Using the novel and tight lower bound from Proposition \[prop:average-SE\], in this work the ASE is obtained as $$\mathrm{ASE} = \lambda K \, \underline{\mathrm{SE}}. \label{eq:def-ASE}$$ To specify the APC, we begin by observing that with the adopted power-control policy the average radiated power per UE is $$\frac{S-(\beta K -1)}{S} \mathbb{E} \{ p_{ji} \} 
= \left( 1- \frac{\beta K - 1}{S} \right)  \rho  \omega \frac{ \Gamma({\alpha}/{2} +1 ) }{(\pi \lambda)^{\alpha/2}}$$ where we have used and the fact that each user transmits one pilot symbol and $S-\beta K$ data symbols per block. Then, we observe that the APC must account not only for the radiated power, but also for the dissipation in analog hardware, digital signal processing, backhaul signaling, and other overhead costs (such as cooling and power supply losses). A detailed and generic model that takes all these factors into account was recently proposed in [@Bjornson2015a] and is such that (using the same notation as in [@Bjornson2015a] for simplicity) the APC is computed as $$\begin{aligned}
 \notag
\mathrm{APC}  &=  \lambda \bigg( \left( 1- \frac{\beta K - 1}{S} \right) \frac{\rho \omega}{\eta}  \frac{ \Gamma({\alpha}/{2} +1 ) }{ (\pi \lambda)^{\alpha/2} } K  \\ & +\! \mathcal{C}_0 \! +\! \mathcal{C}_1 K \!+\! \mathcal{D}_0 M \!+\! \mathcal{D}_1 M K \bigg)  \!+\! \mathcal{A} \cdot \mathrm{ASE} \label{eq:def-APC}\end{aligned}$$ where $\eta \in (0,1]$ is the linear power amplifier efficiency, $\mathcal{C}_0$ models the static power consumption at a BS, and $\mathcal{D}_0 M$ models the power consumption of the BS transceiver chains, which scales with the number of BS antennas. Moreover, $\mathcal{C}_1 K + \mathcal{D}_1 M K$ models the power consumed at the UEs and also by the signal processing tasks at the BS, except for the coding and decoding that are proportional to the number of bits (with proportionality constant $\mathcal{A}$). Notice that has an intuitive polynomial structure in $M$, $K$, and $\rho$. We refer to [@Bjornson2015a] for further modeling details. The forthcoming analysis holds for any positive values of the above parameters, but some examples are later given in Table \[table\_parameters\_hardware\].

The objective of this work is, for any given frame length ($S$), propagation parameters ($\alpha,\omega$), and hardware characteristics ($\eta,\epsilon, \mathcal{A},\mathcal{C}_0,\mathcal{C}_1,\mathcal{D}_0,\mathcal{D}_1$), to find the tuple of parameters $\theta = (\beta, \rho,\lambda,K,M)$ that solves the following constrained EE maximization problem: $$\label{eq:main-optimization-problem}
\begin{split}
{{\underset{{\theta \in \Theta}}{\mathrm{maximize}}}} &\quad \mathrm{EE}(\theta) = \frac{\mathrm{ASE}(\theta)}{\mathrm{APC}(\theta)}\\
\mathrm{subject} \,\, \mathrm{to}  &   \quad \underline{\mathrm{SINR}} = \gamma 
\end{split}$$ where $\Theta$ is the feasible parameter set defined as $$\Theta = \{\theta: \, \rho\geq 0, \lambda \geq 0, \beta \geq 1, (M,K) \in \mathbb{Z}_+, K\beta \leq S\}$$ with $K\beta \leq S$ being the upper limit on the pilot signaling overhead. The parameter $\gamma>0$ in is used to impose an average SE constraint of $\log_2(1+\gamma)$ \[bit/symbol/user\], where the average is computed with respect to both BS and UE locations. This constraint is needed to obtain a network with an acceptable throughput, since unconstrained EE maximization often leads operating points with very low SEs per UE; this is later on illustrated in Fig. \[figureBSdensity\] where the EE increases as $\gamma$ decreases.

In the remainder of this section, we analyze the EE maximization problem in to expose fundamental behaviors and to develop an algorithm for solving it.

Feasibility {#subsec:feasibility}
-----------

Due to the unavoidable inter-cell interference in cellular networks, the optimization problem is only feasible for some values of $\gamma$. This feasible range is obtained as follows:

The optimization problem is feasible if $$\label{eq:SINR-bound}
\gamma < \frac{S (\alpha -1) (1-\epsilon^2) }{ 1 + \epsilon^2 S (\alpha -1)}.$$

Observe that $\underline {\mathrm{SINR}}$ in is a monotonically increasing function of $M$. Therefore, an upper limit can be computed by letting $M \rightarrow \infty$. This yields $$\underline {\mathrm{SINR}}
\xrightarrow{M \to \infty}
 \frac{(1-\epsilon^2)^2}{(1-\epsilon^2) \big( \frac{1}{\beta( \alpha-1)} + \epsilon^2 \big)},$$ which is an increasing function of the optimization variable $\beta$. Under the constraint $\beta K \leq S$, the maximal value of $\beta$ is $S$ and this value is obtained for $K=1$. Putting all these facts together and using simple algebra lead to .

The above lemma shows that the maximal SINR level is limited only by the hardware impairments, through $\epsilon$, and by the severity of the pilot contamination, which is determined by the pathloss exponent $\alpha$ and the coherence block length $S$. These are the only limiting factors as $M \rightarrow \infty$. Assume for example that $\epsilon = 0.05$ [@Bjornson2013c; @Bjornson2014a] and consider the relatively conservative propagation parameters $\alpha=3$ and $S = 200$. For these numbers, the upper limit of the average SE per UE is $\log_2(1+199.5) \approx 7.65$, which is substantially higher than the SE of contemporary systems [@Holma2011a]. This means that the optimization problem is feasible in most cases of practical interest.

$$\begin{aligned}
  \label{eq:EEbeta}  \tag{23}
\mathrm{EE}(\beta^{\star}) = \frac{ K (1  \! -\! \frac{K}{S} \frac{B_1 \gamma}{M (1-\epsilon^2)^2  -B_2 \gamma} )  \log_2(1 \! + \! \gamma)}{ ( 1  \! + \!  \frac{1}{S} 
 \! -\!  \frac{K}{S} \frac{B_1 \gamma}{M (1-\epsilon^2)^2  -B_2 \gamma} ) \frac{K \rho \omega}{\eta}\frac{ \Gamma(\frac{\alpha}{2} +1 ) }{ (\pi \lambda)^{\alpha/2} }   \! + \!  \mathcal{C}_0  \! + \!  \mathcal{C}_1 K  \! + \!  \mathcal{D}_0 M  \! + \! \mathcal{D}_1 M K    \! + \!  \mathcal{A} K (1  \! -\! \frac{K}{S} \frac{B_1 \gamma}{M (1-\epsilon^2)^2  -B_2 \gamma} )  \log_2(1 \! + \! \gamma)  }\end{aligned}$$

$$\begin{aligned}
 \label{eq:EE-limit}  \tag{25}
\mathrm{EE}_{\infty} =  \frac{ K (1 - \frac{K}{S} \frac{\bar{B}_1 \gamma}{M (1-\epsilon^2)^2  -\bar{B}_2 \gamma} )  \log_2(1+\gamma)}{  \mathcal{C}_0 + \mathcal{C}_1 K + \mathcal{D}_0 M + \mathcal{D}_1 M K   + \mathcal{A} K (1 - \frac{K}{S} \frac{\bar{B}_1 \gamma}{M (1-\epsilon^2)^2  -\bar{B}_2 \gamma} )  \log_2(1+\gamma)  }\end{aligned}$$

Optimal Pilot Reuse Factor $\beta$
----------------------------------

We begin by deriving the optimal value of the pilot reuse factor $\beta$ when the other optimization variables are fixed.

\[th:optimal-beta\] Consider any set of $\{\rho, \lambda, M, K\}$ for which the problem is feasible. The SINR constraint in is satisfied by selecting $$\label{eq:beta-optimal}
\beta^\star = \frac{B_1 \gamma}{M(1-\epsilon^2)^2 - B_2 \gamma}$$ where $$\begin{aligned}
 \label{eq:B1}
B_1 &= \left( \frac{4K}{(\alpha-2)^2 } +
\frac{K +M (1-\epsilon^2) }{ \alpha-1} +  \frac{ 2 ( K + \frac{\sigma^2}{  \rho} )}{ \alpha-2} \right) \\
B_2 &=  \left( K + \frac{\sigma^2}{  \rho} +  \frac{ 2K}{\alpha-2} \right) \left(  1  + \frac{\sigma^2}{\rho} \right) + (1-\epsilon^2) \epsilon^2 M . \label{eq:B2}\end{aligned}$$

By gathering the terms that contain $\beta$ and using , the SINR constraint can be rewritten as $$\label{eq:gamma-equation}
\gamma = \frac{ M (1-\epsilon^2)^2  }{ B_1/\beta + B_2   }$$ with $B_1$ and $B_2$ given by and . We then obtain by solving for $\beta$.

The above theorem provides insights on how the EE-optimal pilot reuse factor $\beta^\star$ depends on the other system parameters. Firstly, recall that increasing $\beta$ translates into allocating a larger portion of each UL block for pilot transmission, so that each pilot symbol is on average only used in $1/\beta$ of the cells in the network. This leads to higher channel estimation accuracy and less coherent pilot contamination. Secondly, $\beta^\star$ is an increasing function of $B_1$ and also of $B_2$, since a larger $B_2$ makes the denominator smaller. Consequently, Theorem \[th:optimal-beta\] shows that to guarantee a certain average SINR, $\beta^\star$ must increase with $K$. This is intuitive since more UEs per cell means more inter-cell interference, which can be partially suppressed by increasing the estimation accuracy and reducing the pilot contamination; namely, using a larger $\beta$. Similarly, $\beta^\star$ decreases with $\rho$ since higher transmit powers reduce the detrimental impact of noise, leading to higher estimation accuracy and to a more interference-limited regime. Moreover, $\beta^\star$ is a decreasing function of $M$ since an improved array gain makes the system less sensitive to interference and estimation errors. Increasing the pathloss exponent $\alpha$ leads to a smaller $\beta^\star$ (since $B_1$ and $B_2$ are reduced), which is natural since inter-cell interference decays more quickly. Thirdly, the fact that $\beta \geq 1$ implies that we can only achieve values of $\gamma$ for which $ \frac{B_1 \gamma}{M(1-\epsilon^2)^2 - B_2 \gamma} \geq 1$, otherwise even $\beta=1$ would provide an SINR higher than $\gamma$.

Optimal BS Density and Radiated Power
-------------------------------------

Next, we optimize the BS density and radiated power. Plugging the optimal $\beta^\star$ from Theorem \[th:optimal-beta\] into , the EE maximization problem reduces to $$\begin{aligned}
 \label{eq:main-optimization-problem-modified}
{{\underset{{\rho, \lambda \geq 0, \,\,  M,K \in \mathbb{Z}_+}}{\mathrm{maximize}}}} &\quad \mathrm{EE}(\beta^{\star}) \\
\,\,\,\, \,\, \mathrm{subject} \,\, \mathrm{to} \,\,\, & \quad \frac{B_1 \gamma}{M(1-\epsilon^2)^2 - B_2 \gamma} \geq 1 \notag 
\\
& \quad  \frac{B_1 \gamma}{M (1-\epsilon^2)^2 -B_2 \gamma}  \leq \frac{S}{K} \notag\end{aligned}$$ with $\mathrm{EE}(\beta^{\star}) $ given in at the top of the page.

The optimal values for the BS density $\lambda$ and the power-control coefficient $\rho$ are given in the following theorem.

\[th:optimal-lambda\] Define $\rho = \lambda \tilde{\rho}$ for $\tilde{\rho}>0$ and consider any set of $\{\tilde{\rho},M,K\}$ for which the problem is feasible. Then, $\mathrm{EE}(\beta^{\star})$ is a monotonically increasing function of $\lambda$ and is maximized as $\lambda \to \infty$. The average transmit power per UE in will then go to zero.

The objective function $\mathrm{EE}(\beta^{\star})$ in is a monotonically increasing function of $\lambda$ since the transmit power term decreases with $\lambda$ proportionally to $\rho/\lambda^{\alpha/2} = \tilde{\rho}/\lambda^{\alpha/2-1}$. The EE is also an increasing function of the factor $(1  \! -\! \frac{K}{S} \frac{B_1 \gamma}{M (1-\epsilon^2)^2  -B_2 \gamma} )$ and this expression also increases with $\lambda$ since $B_1$ and $B_2$ are decreasing with $\rho = \lambda \tilde{\rho}$. Therefore, $\mathrm{EE}(\beta^{\star})$ is maximized as $\lambda \to \infty$.

This theorem proves that from an EE perspective it is preferable to have as high BS density as possible. This might be unexpected since smaller cells lead to more interfering UEs in the vicinity of each cell, but this issue is resolved by the assumed power-control policy that gradually reduces the transmit power as the BS density increases. The main consequence of letting $\lambda$ grow large is thus that the transmit power becomes negligible as compared to the circuit power in each cell.

Clearly, an infinitely high BS density is infeasible in practice. However, the numerical results in Section \[sec:simulations\] show that the asymptotic limit is almost achieved already at the modest density of $\lambda = 10$ BS/km$^2$. There are two factors that restrain the practical BS density: 1) the dimensionality of the BS equipment that limits the inter-BS distances; and 2) the UE density that limits the number of BSs that can serve $K$ UEs each. In ultra-dense networks, only a subset of all BSs are turned on at a given point in time; load balancing can be used to make sure that each active BS serves the most energy-efficient number of UEs, while the remaining ones are placed in sleep mode. This is further analyzed in Section \[subsec:given-UEdensity\].

$$\begin{aligned}
 \label{eq:c1}  \tag{31}
\bar{c}^\star=  \frac{a_1  + a_3  + \sqrt{  a_1  a_3  + a_1^2+\frac{ a_1 a_2 a_4 }{a_5} + \frac{a_0  a_3 a_4 }{a_5} - \frac{a_0 a_1  a_4 }{a_5}  - \frac{a_0^2  a_3 a_4  }{a_2 a_5} + \frac{ a_0 a_1 a_3 }{a_2} + \frac{a_0 a_3^2}{a_2} } }{ a_2-a_0 }\end{aligned}$$

Optimal Number of Antennas and UEs per BS {#subsec:opt-M-K}
-----------------------------------------

By using Theorem \[th:optimal-lambda\], the EE maximization problem in further reduces to $$\begin{aligned}
 \label{eq:main-optimization-problem-modified2}
{{\underset{{M,K \in \mathbb{Z}_+}}{\mathrm{maximize}}}} &\quad  \mathrm{EE}_{\infty} \\
\,\,\,\, \,\, \mathrm{subject} \,\, \mathrm{to} \,\,\, & \quad  \frac{\bar{B}_1 \gamma}{M (1-\epsilon^2)^2 -\bar{B}_2 \gamma}  \geq 1 \notag \\
& \quad  \frac{\bar{B}_1 \gamma}{M (1-\epsilon^2)^2 -\bar{B}_2 \gamma}  \leq \frac{S}{K} \notag\end{aligned}$$ where $\mathrm{EE}_{\infty}$ is given in at the top of the page and we have defined $$\begin{aligned}
 \label{eq:B1bar}
\bar{B}_1 &= K \left( \frac{4}{(\alpha-2)^2 } +
\frac{1}{ \alpha-1} +  \frac{ 2}{ \alpha-2} \right)  +
\frac{M (1-\epsilon^2) }{ \alpha-1} , \\
\bar{B}_2 &=  K\left(  1 +  \frac{ 2}{\alpha-2} \right) + (1-\epsilon^2) \epsilon^2 M . \label{eq:B2bar}\end{aligned}$$ To find the optimal values for $M$ and $K$, an integer-relaxed version of is first considered where $M$ and $K$ can be any positive scalars. The integer-valued solutions are then extracted from the relaxed problem. For analytic tractability, we replace $M$ with $\bar{c} = {M}/{K}$, which is the number of BS antennas per UE. For a given $\bar{c}$, the EE-maximizing value of $K$ is found as follows.

\[th:optimal-K\] Consider the optimization problem where $M$ and $K$ are relaxed to be real-valued variables. For any fixed $\bar{c} = {M}/{K}> 0$ such that the relaxed problem is feasible for some $K$, the EE is maximized by $$\begin{aligned}
 \label{eq:optimal-K}
K^{\star} = \frac{ \sqrt{\left(G\mathcal C_{0}\right)^{2}+\mathcal{C}_{0} \mathcal D_{1}\bar c + \mathcal{C}_{0} G \left(\mathcal C_{1} + \mathcal D_{0}\bar c \right)}   - G\mathcal C_{0} }{\mathcal D_{1}\bar c+ G\left(\mathcal C_{1} + \mathcal D_{0}\bar c \right)}\end{aligned}$$ where $$\begin{aligned}
G = \frac{1}{S} \frac{ \left( \frac{4\gamma}{(\alpha-2)^2 } +
\frac{\gamma}{ \alpha-1} +  \frac{ 2\gamma}{ \alpha-2}  
 \right)  + \frac{ \gamma (1-\epsilon^2) }{ \alpha-1}  \bar{c} }{ 
 (1-\epsilon^2) \left( 1- (1+\gamma) \epsilon^2 \right) \bar{c}
  - \left(  1 +  \frac{ 2}{\alpha-2} \right) \gamma }.\end{aligned}$$

If $\bar{c}$ is given, then we want to maximize (using the notation of the theorem and $\mathcal{R}=\log_2(1+\gamma)$) $$\begin{aligned}
 \label{eq:utility-K}
 \frac{ K (1 - K G ) \mathcal{R}}{  \mathcal{C}_0 + ( \mathcal{C}_1+ \mathcal{D}_0 \bar c ) K  + \mathcal{D}_1 \bar c K^2   + \mathcal{A} K (1 - K G ) \mathcal{R}  } .\end{aligned}$$ It is straightforward to show that is a quasi-concave function of $K$ (e.g., using the approach in the proof of Lemma 3 in [@Bjornson2015a]). The maximizing value $K^{\star}$ is thus obtained by taking the first derivative of and equating to zero. This leads to the expression in . The first constraint in is independent of $K$ and thus fulfilled if the problem was feasible. The solution also satisfies the second constraint, namely $K \leq {1}/{G}$, since the EE in is negative for $K > {1}/{G}$ and thus is not maximized in thus range (and the zero objective function for $K=0$ would then be higher).

Similarly, if $K$ is fixed, then the EE-maximizing value of $\bar{c}$ is obtained as follows.

\[th:optimal-c\] Consider the optimization problem where $\bar{c} = {M}/{K}$ and $K$ are relaxed to be real-valued variables. For any fixed $K> 0$ such that the relaxed problem is feasible, the EE is maximized by $\bar{c}^\star$ in , at the top of the page, if this solution satisfies the first inequality constraint in . Otherwise, the EE is maximized by $$\begin{aligned}
 \label{eq:c2} 
\bar{c}^\star= \frac{ \gamma \left( 1 +\frac{4}{(\alpha-2)^2 } + \frac{1}{ \alpha-1} +  \frac{ 4}{ \alpha-2}   \right)}{ (1-\epsilon^2) \left( 1- (1+\gamma) \epsilon^2 \right) - \frac{ \gamma (1-\epsilon^2) }{\alpha-1}  }.\end{aligned}$$ The following parameters are used in and : $$\begin{aligned}
a_0 &= \frac{ \gamma K (1-\epsilon^2) }{S(\alpha-1)} \\ 
a_1 &=  \frac{K}{S} \left( \frac{4\gamma}{(\alpha-2)^2 } + \frac{\gamma}{ \alpha-1} +  \frac{ 2\gamma}{ \alpha-2}   \right) \\
a_2 &=   (1-\epsilon^2) \left( 1- (1+\gamma) \epsilon^2 \right) \\
a_3 & = \left(  1 +  \frac{ 2}{\alpha-2} \right) \gamma \\
a_4 & =  \mathcal{C}_0 + \mathcal{C}_1 K \\
a_5 & = \mathcal{D}_0 K + \mathcal{D}_1 K^2.\end{aligned}$$

Using the notation introduced in the theorem, the objective function in reduces to $$\begin{aligned}
 \label{eq:utility-c}
 \frac{ K (1 -  \frac{  a_0  \bar{c} + a_1  }{ a_2 \bar{c} - a_3} )  \log_2(1+\gamma)}{  a_4 + a_5 \bar{c}      + \mathcal{A} K (1 -  \frac{  a_0  \bar{c} + a_1  }{ a_2 \bar{c} - a_3} )   \log_2(1+\gamma)  } \end{aligned}$$ which can be easily shown to be a quasi-concave function of $\bar{c}$. The maximizing value $\bar{c}^{\star}$ is thus obtained by taking the first derivative of and equating to zero. This yields the solution in . Note that $\bar{c}^\star$ needs to satisfy the first constraint in . If this is not the case with , then the EE is monotonically increasing for all feasible $\bar{c}$ and the largest value is obtained when there is equality in the first constraint, which occurs at the value in . Finally, we notice that $\bar{c}^\star$ satisfies the second constraint in automatically because the maximum cannot give a negative EE.

The above two theorems show how $K$ and $\bar{c}$ (and also $M=\bar{c}K$) are related at the EE-optimal points. It turns out that $K^{\star}$ decreases with $\bar{c}$, proportionally to $\sqrt{1/\bar{c}}$ when $\bar{c}$ grows large. Similarly, $M^{\star} = \bar{c}^{\star}\!K$ is found to increase with $K$. The intuition is that more BS antennas should be deployed to control the interference when more UEs are served. In addition to this, the results of Theorems \[th:optimal-K\] and \[th:optimal-c\] reveal how the power consumption coefficients impact $K^{\star}$ and $M^{\star}$, respectively. In particular, from it is found that $K^{\star}$ increases with the static power consumption $\mathcal{C}_0$, while it decreases with $\mathcal{C}_1$, $\mathcal{D}_0$, and $\mathcal{D}_1$ that are the terms of the APC that increase with $K$ and $M$. Similarly, $M^{\star}$ increases with $\mathcal{C}_0$ and $\mathcal{C}_1$, but decreases with $\mathcal{D}_0$ and $\mathcal{D}_1$. The simple intuition behind these scaling behaviors is that more hardware should be turned on (i.e., BS antennas and UEs) only if the increase in circuit power has a marginal effect on the total APC. Similarly, with a larger static consumption we can afford to turn on more BS antennas and UEs since the relative power cost is lower.

Using Theorems \[th:optimal-K\] and \[th:optimal-c\], we can devise an alternating optimization algorithm to solve the integer-relaxed EE maximization problem:

-2mm

[|c|c|c||c|c|c|]{} **System Parameter & **Symbol & **Value & **Hardware Parameter & **Symbol & **Value\
************

Coherence block length & $S$ & $400$ & Coding/decoding/backhaul & $\mathcal{A}$ & $1.15 \cdot 10^{-9} \, \mathrm{[J/bit]}$\
Pathloss exponent & $\alpha$ & $3.76$ & Static power consumption & $\mathcal{C}_0$ & $10 \, \mathrm{W} \cdot \tau \,\, \mathrm{[J/symbol]}$\
Propagation loss at 1 $\textrm{km}$ & $\omega$ & $130 \, \mathrm{dB}$ & Circuit power per active user & $\mathcal{C}_1$ & $0.1 \, \mathrm{W} \cdot \tau \,\, \mathrm{[J/symbol]}$\
Power amplifier efficiency & $\eta$ & $0.39$ & Circuit power per BS antenna & $\mathcal{D}_0$ & $0.2 \, \mathrm{W} \cdot \tau \,\, \mathrm{[J/symbol]}$\
Level of hardware impairments & $\epsilon$ & 0.05 & Signal processing coefficient & $\mathcal{D}_1$ & $1.56 \cdot 10^{-10} \, \mathrm{[J/symbol]}$\
Symbol time & $\tau$ & $\frac{1}{2 \cdot 10^7} \,\, \mathrm{[s/symbol]}$ & Noise variance & $\sigma^2$ & $10^{-20} \, \mathrm{[J/symbol]}$\

1.  Find a feasible starting point $(M,K)$ to ;

2.  Optimize $K$ for a fixed $M$ using Theorem \[th:optimal-K\];

3.  Optimize $M$ for a fixed $K$ using Theorem \[th:optimal-c\];

4.  Repeat 2)–3) until convergence is achieved.

This algorithm converges since the EE has a finite upper bound and the EE increases monotonically in each iteration. In fact, it converges to the global optimum of the relaxed problem.

The integer-relaxed version of the EE maximization problem is quasi-concave. Hence, the alternating algorithm converges to the global optimal solution $(K^{\star \star},M^{\star \star}) \in \mathbb{R}^2$ to the integer-relaxed version of the problem.

The alternating optimization algorithm is guaranteed to converge to a local maximum according to [@Grippo2000a Proposition 4]. Suppose that there are multiple local maxima; for example, $(M_1,K_1)$ and $(M_2,K_2)$. By setting $M = \frac{M_2-M_1}{K_2-K_1} K + \frac{K_2 M_1- K_1 M_2}{K_2-K_1}$ and varying $K$, we can search on the line between these local maxima. For this choice of $M$, it is straightforward to show that the EE is a quasi-concave function of $K$. Hence, only one of the points can be a true local maximum, which is a contradiction to the existence of multiple local maxima; thus, the algorithm converges to the only local/global optimum of the problem.

The real-valued solution $(K^{\star \star},M^{\star \star})$ obtained from the alternating optimization algorithm is a good starting point for finding the integer-valued global optimum to . In particular, the quasi-concavity implies that the integer-solution is contained in a convex level set around $(K^{\star \star},M^{\star \star})$. In many cases, it is one of the four integer points obtained by respectively taking the floor and ceiling of $K^{\star \star}$ and $M^{\star \star}$. In general, the global optimum is obtained by searching through all integers in the vicinity of the real-valued solution, keeping in mind that the EE is quasi-concave in all directions.

To summarize, the original problem formulation in has been solved through the following steps: $i$) by selecting $\beta$ to satisfy the SINR constraint; $ii$) by letting $\lambda \to \infty$ (which also makes $\rho \to 0$); $iii$) by devising an alternating optimization algorithm that provides the real-valued $M$ and $K$ that maximize the EE; $iv$) finally, by searching through the integer points in the vicinity of the real-valued solution, and capitalizing on the quasi-concavity. In the process of solving , Theorems \[th:optimal-beta\]–\[th:optimal-c\] have also exposed the fundamental interplay among the optimization variables and how they depend on the hardware characteristics and propagation parameters.

Numerical Examples {#sec:simulations}
==================

Numerical results are used in this section to validate the average SE expression provided in Proposition \[prop:average-SE\] and the fundamental results established in Section \[sec:problem\]. The results are obtained for the parameter setting summarized in Table \[table\_parameters\_hardware\], where we assume for instance a coherence block length of $S=400$ (e.g., $T_c=4$ ms and $W_c = 100$ kHz) and a bandwidth of 20 MHz (giving the symbol time $1/(2\cdot 10^7)$ s). The hardware parameters are inspired by a variety of prior works; for example, [@Bjornson2015a] and references therein. The simulations were performed using Matlab and the code is available for download at <https://github.com/emilbjornson/maximal-EE>, which enables reproducibility as well as simple testing of other parameter values.

Optimizing the Energy Efficiency
--------------------------------

We begin by recalling that Theorem \[th:optimal-lambda\] proves that the EE maximization problem in is solved when the BS density is infinitely large; that is, $\lambda \to \infty$. We now illustrate how large $\lambda$ needs to be for applying this asymptotic result in practice. To this end, Fig. \[figureBSdensity\] shows the EE as a function of $\lambda$, and the other optimization variables are optimized numerically for each given value of $\lambda$. Three different SINR constraints are considered in Fig. \[figureBSdensity\]: $\gamma \in \{  1, \, 3, \, 7 \}$ which corresponds to the average SEs $\log_2(1+\gamma) \in \{  1, \, 2, \, 3 \}$. In all three cases, the EE is computed using both the lower bound on the average SE in Proposition \[prop:average-SE\] and an upper bound obtained by averaging over the instantaneous SE derived in Lemma \[lemma:SINR-ideal\] using Monte Carlo simulations.[^6]

Several important observations can be made from the results presented in Fig. \[figureBSdensity\]. Firstly, there is only a small gap (that reduces as $\gamma$ takes larger values) between the lower and upper bounds, and the curves behave exactly the same for any value of $\gamma$. This validates the accuracy of the SE expression provided in Proposition \[prop:average-SE\]. Secondly, the EE can be greatly improved by increasing the BS density, meaning that small cells are a promising solution for maximal EE deployment. However, the gain from increasing the BS density saturates in the interval from $\lambda = 10$ to $\lambda = 100$ BS/km$^2$, which roughly corresponds to an average inter-BS distance of 100–315 meters. This is not small compared to contemporary urban deployments. In other words, EE maximization based on letting $\lambda \rightarrow \infty$ is expected to give representative results in most practical dense deployments. Thirdly, we note that the EE decreases as $\gamma$ increases, which is why it is important to specify a target SINR; otherwise the EE maximizing operating point might be very spectrally inefficient and therefore useless from a practical user service perspective.

![Energy efficiency (in Mbit/J) as a function of the BS density in $\mathrm{BS}/\mathrm{km}^2$, for different SINR constraints. Results are obtained by both computing upper bounds from Monte-Carlo simulations and using the closed-form lower bound from Proposition\[prop:average-SE\].[]{data-label="figureBSdensity"}](figureBSdensity){width=".95\columnwidth"}

-3mm

-3mm

We proceed by further studying the impact of varying $M$ and $K$. To this end, we let $\lambda \rightarrow \infty$ as prescribed by Theorem \[th:optimal-lambda\] and fix the SINR constraint at $\gamma = 3$ or, equivalently, the average SE per data symbol at $\log_2(1+\gamma) = 2$. Fig. \[figure3d\] shows the EE lower bound as a function of $M$ and $K$ when $\beta$ is optimized according to . The global EE maximum gives an EE of $10.156$ Mbit/J and is achieved by $(M^\star,K^\star) = (91,10)$ using the pilot reuse factor $\beta^\star = 7.08$. Interestingly, this is a configuration that falls within the class of massive MIMO setups [@Marzetta2010a; @Larsson2014a; @Vieira2014a]. The intuition behind this result is that the strong inter-cell interference in dense deployments is efficiently mitigated by MRC when the BSs are equipped with many antennas and also by using a substantial pilot reuse factor (to protect against channel estimation errors and pilot contamination).

The results obtained with the alternating optimization algorithm from Section \[subsec:opt-M-K\] are also shown in Fig. \[figure3d\]. The initialization point was set to $(M,K) = (20,1)$. As seen, the algorithm converges after three iterations to the real-valued solution $(M^{\star \star},K^{\star \star}) = (91.6,10.1)$ with an EE of $10.157$ Mbit/J. This real-valued operating point gives only a 0.009% higher EE than the integer-valued solution, thus showing that the EE performance is quite flat around the global optimum.

To study the global optimum even further, the pie diagram in Fig. \[figurePieDiagram\] shows the relative size of each term in the APC of (the transmit power is not shown since it takes negligible values). The dominating terms are the static power consumption $\mathcal{C}_0$ and the consumption $\mathcal{D}_0 M$ of the BS transceiver chains. These seem to be the main factors to improve in order to make the hardware more energy-efficient. Observe that the above results might be different for other values of $\mathcal{C}_0$, $\mathcal{C}_1$, $\mathcal{D}_0$, and $\mathcal{D}_1$.

![Energy efficiency (in Mbit/J) for $\gamma = 3$. The global optimum is star-marked, while the convergence of the alternating algorithm from Section \[subsec:opt-M-K\] is indicated with circles.[]{data-label="figure3d"}](figure3d){width=".95\columnwidth"}

-3mm

![Relative size of each term in the APC of at the global optimum from Fig. \[figure3d\].[]{data-label="figurePieDiagram"}](figurePieDiagram){width=".7\columnwidth"}

-3mm

Impact of Transceiver Hardware Impairments
------------------------------------------

Next, we consider the saturation regime where $\lambda \to \infty$ and exemplify the impact of transceiver hardware impairments on the EE. Fig. \[figureImpairments\] shows the EE as a function of $\epsilon$, the level of hardware impairments. As expected, the EE decreases with $\epsilon$ since the desired signal power in decays as $(1-\epsilon^2)^2$. The loss is marginal for $\gamma = 1$, but it can be relatively large when $\gamma$ increases. This is in line with the previous results in [@Bjornson2013c], wherein it is shown that hardware impairments greatly affect the channel capacity in the high SNR regime (i.e., for large $\gamma$) while their impact is negligible in the low SNR regime. The results of Fig. \[figureImpairments\] indicate that for the investigated network the EE loss due to hardware impairments is negligible for $\epsilon \leq 0.1$ when $\gamma \in \{1, \, 3 \}$. Since these values correspond to the operating points that give the highest EE (see Fig. \[figureBSdensity\]), we may conclude that modest levels of hardware impairments have a negligible impact on networks designed for high EE.

![Energy efficiency (in Mbit/J) as a function of the level of transceiver hardware impairments, $\epsilon$, for different SINR constraints.[]{data-label="figureImpairments"}](figureImpairments){width=".95\columnwidth"}

-3mm

EE Maximization for a Given UE Density {#subsec:given-UEdensity}
--------------------------------------

So far, the BSs have been deployed to match an unlimited heterogeneous UE distribution. In particular, we have assumed that each BS serves $K$ UEs such that a BS density of $\lambda$ BS/km$^2$ corresponds to an average of $K \lambda$ UEs per km$^2$. Hence, as the BS density grows large we also let the average UE density grow large. To validate the plausibility of this model, suppose that we instead deploy the BSs to match a fixed average UE density of $\mu$ UE/km$^2$. Mathematically, this amounts to solving with the additional constraint $$\label{eq:UE-density-constraint}
\mu = K \lambda.$$ We will study how such an extra constraint affects the results, taking into account that future average UE densities from $\mu=10^2$ UE/km$^2$ (in rural areas) to $\mu=10^5$ UE/km$^2$ (in shopping malls) have been predicted in the METIS project [@METIS_D11_short].

![Energy efficiency (in Mbit/J) vs. the UE density $\mu$. The EE is optimized according to with the extra constraint $\mu = K \lambda$, or only with respect to $(\lambda,\beta,\rho)$ for given $M$ and $K$.[]{data-label="figureUEdensity_EE"}](figureUEdensity_EE){width=".95\columnwidth"}

-3mm

![Optimized BS density (in $\mathrm{BS/km}^2$) vs. the UE density (in $\mathrm{UE/km}^2$). The system is optimized in the same way as in Fig. \[figureUEdensity\_EE\].[]{data-label="figureUEdensity_BS"}](figureUEdensity_BS){width=".95\columnwidth"}

-3mm

Fig. \[figureUEdensity\_EE\] shows the EE has a function of the UE density $\mu $ for the average SINR level $\gamma = 3$, while Fig. \[figureUEdensity\_BS\] shows the corresponding BS density. The design parameters $M$, $K$, $\beta$, $\lambda$, and $\rho$ are optimized numerically according to , with the additional constraint . Two reference cases are also considered: single-user single-input multiple-output (SIMO) transmission with $(M,K) = (10,1)$; and massive multi-user MIMO transmission with $(M,K) = (91,10)$, which was shown in Fig. \[figure3d\] to be optimal as $\lambda \to \infty$. Clearly, only $\beta$, $\lambda$, and $\rho$ are optimized in the two reference cases.

Several important observations can be made from the results presented in Figs. \[figureUEdensity\_EE\] and \[figureUEdensity\_BS\]. Firstly, the EE level becomes independent of the UE density when $\mu$ is large enough; the saturation occurs for $\mu \geq 100$ in the optimal and the fixed massive MIMO cases, while it occurs for $\mu \geq 2$ in the single-user SIMO case. The future UE density range $\mu \in [10^2, 10^5]$ predicted in [@METIS_D11_short] is fully covered in these saturation regimes. In all practical relevant cases we can, according to Fig. \[figureUEdensity\_BS\], simply scale the BS density linearly with the UE density, while using the same number of BS antennas and UEs per BS. Similarly, one can turn on and off BSs to account for variations in the UE load. Secondly, the fixed massive MIMO configuration $(M,K) = (91,10)$ achieves the maximal EE in the saturation regime, while it is inefficient as compared to single-user MISO for low UE densities (e.g., for $\mu \leq 8$). In contrast, single-user SIMO transmission performs reasonably well at low UE densities, but saturates earlier and at an EE level that is $3.14 \times$ lower than the maximal EE. More importantly, Fig. \[figureUEdensity\_BS\] shows that the single-user case requires a $10 \times$ higher BS density in the saturation regime, which might greatly reduce the deployment cost.

In summary, the transmission power appears to be negligible, as compared to the circuit power, for networks that can handle future UE densities. By adding massive multi-user MIMO capability to the BSs, the EE can be increased by a few hundred percentages, and the even more substantial benefit is that it can reduce the BS density with an order of magnitude; this is likely to be a key property to achieve a cost and energy efficient network densification.

Conclusion {#sec:conclusion}
==========

Network densification is the key for achieving high EE in cellular networks, and can be realized by either having many BSs or many antennas per BS. In order to find the optimal densified network configuration, we formulated an UL EE maximization problem under the assumption of a stochastic BS deployment based on Poisson point processes. By deriving a new lower bound on the average SE in the network and using a state-of-the-art power consumption model, the EE expression became tractable and was maximized analytically with respect to the density of BSs, the transmit power levels, the number of BS antennas and users per cell, and the pilot reuse factor. The closed-form expressions provided general guidelines on the optimal operating regimes and exposed the fundamental interplay among the optimization variables, hardware characteristics, and propagation environment. The analysis focused on the UL, but similar results are expected in the DL due to the UL/DL duality concept.

The analysis shows that reducing the cell size is undoubtable the way towards high EE, but the positive effect of increasing the BS density saturates when the circuit power dominates over the transmission power. A further leap in EE can typically be achieved by adding extra BS antennas to multiplex several UEs per cell; the numerical examples resulted in an EE maximum with 91 BS antennas and 10 UEs, which is a massive MIMO setup. The EE gains come from suppressing intra-cell interference and by sharing the circuit power costs among multiple UEs. We stress that a hundred antennas can be deployed also at a small-size BS, since we consider simple dipoles and possibly higher frequencies than in contemporary networks. Moreover, the analysis shows that a large pilot reuse factor can protect against inter-cell interference and can be tailored to guarantee a certain average SE.

While this paper focused on the ratio between spectral efficiency and energy costs, it is straightforward to also include the deployment cost, site renting, and other economical factors; the annual depreciation expense can be turned into an energy-equivalent form (by multiplying with the energy cost in J/Euro) and incorporated in the coefficients $\mathcal{C}_0$, $\mathcal{C}_1$, and $\mathcal{D}_0$ in . This is likely to mainly increase the constant term $\mathcal{C}_0$, thus pushing the optimal operating point towards having fewer BSs, with more antennas and UEs per cell.

$$\begin{aligned}
 \notag
&\mathbb{E} \{ \mathrm{SINR}_{0k}^{-1} \} \\ \notag
& = \frac{ \mathbb{E} \left\{ \left( K + {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } + \frac{\sigma^2}{  \rho}  \right) \left(  1   + \frac{1}{\beta} {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{\rho} \right) \right\}
+ M (1-\epsilon^2) \left(  \mathbb{E} \left\{ \frac{1}{\beta} {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}  
\left( \frac{ d_{jjk}^{\alpha} }{ d_{0jk}^{\alpha}} \right)^2 \right\} + \epsilon^2 \right) }{ M (1-\epsilon^2)^2  } \\ \notag
&=   \frac{1}{ M (1-\epsilon^2)^2} \left(
\bigg( K + \frac{\sigma^2}{  \rho}  \right) \left(  1  + \frac{\sigma^2}{\rho} \right) 
+\mathbb{E} \left\{ {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } \right\} \left(  1  + \frac{\sigma^2}{\rho} \right) +   \left( K  + \frac{\sigma^2}{  \rho}  \right) \mathbb{E} \bigg\{ \frac{1}{\beta} {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} }  \bigg\} \\
&+ 
 \mathbb{E} \bigg\{ \frac{1}{\beta}  \bigg(  {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} }  \bigg) \bigg(   {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} }  \bigg) \bigg\}
+ M (1-\epsilon^2) \mathbb{E} \bigg\{ \frac{1}{\beta} {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}  
\left( \frac{ d_{jjk}^{\alpha} }{ d_{0jk}^{\alpha}} \right)^2 \bigg\} + M (1-\epsilon^2) \epsilon^2 \bigg).  \label{eq:inverse-SINR} \tag{51}\end{aligned}$$

Appendix A: Proof of Lemma \[lemma:SINR-ideal\] {#app:lemma-SINRs .unnumbered}
===============================================

Let ${\mathbf{v}}_{00k} = \nu_{00k} \hat{{\mathbf{h}}}_{00k}$ denote the MRC detector. Then, the received signal is given by $$\begin{aligned}
 \notag
& {\mathbf{v}}_{00k}^{\Htran} {\mathbf{y}}_0 = \sum_{i=1}^{K} \sqrt{p_{0i}} (\sqrt{1-\epsilon^2}s_{0i} + e_{0i}) {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00i}  \\ & + \sum_{j \in \Psi_{\lambda}} \notag
\sum_{i=1}^{K} \sqrt{p_{ji}} (\sqrt{1-\epsilon^2}s_{ji} + e_{ji})  {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{0ji}  + {\mathbf{v}}_{00k}^{\Htran} {\mathbf{n}}_{0} \\ \notag &= \sqrt{p_{0k}} \sqrt{1-\epsilon^2}s_{0k} \mathbb{E} \{ {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \} + 
\sqrt{p_{0k}} e_{0k} {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \\ \notag &+ \sqrt{p_{0k}} \sqrt{1-\epsilon^2}s_{0k} \left( {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} - \mathbb{E} \{ {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \} \right)   &  \\ \notag &+ 
 \sum_{i \neq k} \sqrt{p_{0i}} (\sqrt{1-\epsilon^2}s_{0i} + e_{0i}) {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00i}  + \notag \\ &+ \sum_{j \in \Psi_{\lambda}}
\sum_{i=1}^{K} \sqrt{p_{ji}} (\sqrt{1-\epsilon^2}s_{ji} + e_{ji})  {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{0ji}  + {\mathbf{v}}_{00k}^{\Htran} {\mathbf{n}}_{0}. \label{eq:processed-signal}\end{aligned}$$ The ergodic capacity of this channel is defined as the maximal mutual information, which unfortunately cannot be computed exactly under imperfect CSI. However, we can compute a rigorous lower bound using the well-established bounding technique from [@Jose2011b] and [@Medard2000a]. More precisely, the first term in is treated as the only desired signal while the other terms are uncorrelated with the first one and treated as worst-case Gaussian noise in the detection (a worst-case assumption). This leads to the mutual information $\log_2(1+\mathrm{SINR}_{0k})$ where $\mathrm{SINR}_{0k}$ is the ratio between the desired power and the “noise” power: $$\begin{aligned}
 \notag
&\mathrm{SINR}_{0k} = p_{0k} (1-\epsilon^2) |\mathbb{E} \{ {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \}|^2 \\ & \notag \Big/ \Big(
p_{0k} (1-\epsilon^2) \left( \mathbb{E} \{ | {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} |^2 \} - |\mathbb{E} \{ {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \}|^2 \right)   \\ &+
p_{0k} \epsilon^2 \mathbb{E} \{ |{\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} |^2 \} +
 \sum_{i \neq k} p_{0i} \mathbb{E} \{ |{\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00i} |^2 \} \notag
\\  &+ \sum_{j \in \Psi_{\lambda}}
\sum_{i=1}^{K} p_{ji} \mathbb{E} \{ | {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{0ji} |^2 \}
+ \mathbb{E} \{ |{\mathbf{v}}_{00k}^{\Htran} {\mathbf{n}}_{0}|^2 \}  \Big) \label{eq:SINR-derivation}\end{aligned}$$ where we have used the fact that $\mathbb{E}\{ | \sqrt{1-\epsilon^2}s_{ji} + e_{ji}  |^2  \} = 1-\epsilon^2 + \epsilon^2 = 1$. To proceed further with the computation of the expectations in , we define $$\tau_{0k} = 1   + \sum_{l \in \Psi_{\lambda} }  \chi_{0kl}  \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{ \rho}$$ and notice that the MMSE estimate of ${\mathbf{h}}_{0jk}$ (for $j=0$ or $j  \in \Psi_{\lambda}$ such that $\chi_{0kj} = 1$) takes the form $$\label{eq:MMSE-estimator-impairments}
 \begin{split}
  \hat{{\mathbf{h}}}_{0jk} &= 
\frac{  1}{ \tau_{0k} d_{0jk}^{\alpha} } \sqrt{ \frac{(1-\epsilon^2) d_{jjk}^{\alpha}}{\rho \omega } }{\mathbf{y}}_{0k}^{\textrm{pilot}} \\ &\sim \mathcal{CN}\left( {\mathbf{0}},  \frac{ (1-\epsilon^2)  d_{jjk}^{\alpha}   }{ \tau_{0k} \omega d_{0jk}^{2\alpha}} {\mathbf{I}}_M \right)
\end{split}$$ whereas the independent estimation error is distributed as $$\Delta {\mathbf{h}}_{0jk} \sim \mathcal{CN}\left( {\mathbf{0}}, \frac{1}{\omega d_{0jk}^{\alpha}} \left( 1 -  \frac{ (1-\epsilon^2) d_{jjk}^{\alpha}    }{ \tau_{0k} d_{0jk}^{\alpha}} \right) {\mathbf{I}}_M \right).$$ For normalization purposes, the scaling factor $\nu_{00k}$ is selected as $$\label{eq:nu-expression}
\nu_{00k} = \tau_{0k} \sqrt{\frac{  \omega d_{00k}^{\alpha}  }  { (1-\epsilon^2) \rho M  }}.$$ Therefore, we have that $$\label{eq:express1}
\begin{split}
& p_{0k} |\mathbb{E} \{ {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{00k} \}|^2  \\ &=\rho \omega d_{00k}^{\alpha} 
\left| \sqrt{\frac{  \omega d_{00k}^{\alpha}  }  { (1-\epsilon^2) \rho M  }} \mathbb{E} \{ \tau_{0k}  \hat{{\mathbf{h}}}_{00k}^{\Htran} \hat{{\mathbf{h}}}_{00k} \} \right|^2 \\ & =  (1-\epsilon^2)   M,
\end{split}$$ which follows by using and the orthogonality principle of MMSE estimators as well as the adopted power-control policy in . Since the noise is independent of the channels, we also have $$\label{eq:express2}
\begin{split}
 \mathbb{E} \{ | {\mathbf{v}}_{00k}^{\Htran} {\mathbf{n}}_{0}|^2 \} &= \sigma^2  \mathbb{E} \{ \nu_{00k}^2 \|\hat{{\mathbf{h}}}_{00k}\|^2 \} =  
  \frac{  \sigma^2 }{ \rho } \mathbb{E} \{ \tau_{0k} \}  \\ &=   \frac{  \sigma^2 }{ \rho }  \left( 1   + \frac{1}{\beta} \sum_{l \in \Psi_{\lambda} }  \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{ \rho} \right).
  \end{split}$$ Next, we consider the terms in that are affected by pilot contamination. For convenience, we let $\chi_{0k0}=1$ since the typical UE uses the same pilot as itself. Then, for all users $i$ in cells $j$ that are not causing pilot contamination (i.e., for $i \neq k$ or for $i=k$ and $\chi_{0kj}\ne1$), we have $$\label{eq:no-pilot-cont}
\begin{split}
& p_{ji} \mathbb{E} \{ | {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{0ji} |^2 \} =
\frac{\rho \omega d_{jji}^{\alpha} }{ \omega d_{0ji}^{\alpha} } \mathbb{E} \{ \nu_{00k}^2 \| \hat{{\mathbf{h}}}_{00k}  \|^2 \} \\ &= 
\frac{d_{jji}^{\alpha} }{d_{0ji}^{\alpha}}
 \mathbb{E} \{ \tau_{0k} \} = \frac{d_{jji}^{\alpha} }{d_{0ji}^{\alpha}} \left( 1   + \frac{1}{\beta} \sum_{l \in \Psi_{\lambda} }  \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{ \rho} \right)
 \end{split}$$ where we have taken into account that $\hat{{\mathbf{h}}}_{00k}$ is independent of ${\mathbf{h}}_{0ji}$ and the last equality follows from averaging over the stochastic pilot selection (with $ \mathbb{E} \{ \chi_{0kl} \} = 1/\beta$ by design).

For all users $i$ in cells $j$ such that $i=k$ and $\chi_{0kj}=1$, the channel estimate $\hat{{\mathbf{h}}}_{00k}$ is correlated with ${\mathbf{h}}_{0ji}$ due to pilot contamination and thus we obtain $$\begin{aligned}
 \label{eq:with-pilot-cont}
& p_{jk} \mathbb{E} \{ | {\mathbf{v}}_{00k}^{\Htran} {\mathbf{h}}_{0jk} |^2 \} \\ &= \rho \omega d_{jjk}^{\alpha}  \Bigg( \frac{d_{jjk}^{\alpha}  d_{00k}^{\alpha} }{d_{0jk}^{2\alpha} } \mathbb{E} \{ \nu_{00k}^2 | \hat{{\mathbf{h}}}_{00k}^{\Htran} \hat{{\mathbf{h}}}_{00k} |^2 \} \notag \\ & \qquad + \mathbb{E} \{ \nu_{00k}^2 | \hat{{\mathbf{h}}}_{00k}^{\Htran} \Delta {\mathbf{h}}_{0jk} |^2 \}  \Bigg) \notag\\
&=  \omega d_{jjk}^{\alpha}  \Bigg[ \frac{d_{jjk}^{\alpha}  d_{00k}^{\alpha} }{d_{0jk}^{2\alpha} }
 \frac{ (1-\epsilon^2 )  }{   \omega d_{00k}^{\alpha} }  (1+M) \notag \\ & \qquad + 
\frac{1}{\omega d_{0jk}^{\alpha}} \left(  \mathbb{E} \{ \tau_{0k} \} -  \frac{ (1-\epsilon^2) d_{jjk}^{\alpha}    }{ d_{0jk}^{\alpha}} \right) \Bigg] \notag \\
 &=  \left( \frac{d_{jjk}^{\alpha} }{d_{0jk}^{\alpha} } \right)^2 \! \! (1-\epsilon^2) M
  +  \frac{d_{jjk}^{\alpha} }{d_{0jk}^{\alpha} }  \! \left( 1   + \frac{1}{\beta} \sum_{l \in \Psi_{\lambda} }  \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } + \frac{\sigma^2}{ \rho} \right)\end{aligned}$$ where the first equality follows from the expansion ${\mathbf{h}}_{0jk} = \hat{{\mathbf{h}}}_{0jk} + \Delta {\mathbf{h}}_{0jk} $ and the fact that $\hat{{\mathbf{h}}}_{0jk} = \sqrt{\frac{d_{jjk}^{\alpha}  d_{00k}^{\alpha} }{d_{0jk}^{2\alpha} }
} \hat{{\mathbf{h}}}_{00k}$. The second equality follows from using [@Bjornson2015b Lemma 2] and , whereas the third equality simplifies the expression taking into account that $ \mathbb{E} \{ \chi_{0kl} \} = 1/\beta$. Since the second term in is equal to , it follows that it is only the first term in that comes from pilot contamination and its probability to appear is $ \mathbb{E} \{ \chi_{0kl} \} = 1/\beta$. Plugging – into the result in the lemma follows.

Appendix B: Proof of Proposition \[prop:average-SE\] {#app:prop:average-SE .unnumbered}
====================================================

The goal is to compute a tractable lower bound on the average SE $( 1 - \frac{B}{S} ) \mathbb{E}\{ \log_2 ( 1 + \mathrm{SINR}_{0k}  ) \}$, where the expectation is taken with respect to the PPP and UE locations. To this end, the Jensen’s inequality is applied as $\mathbb{E}\{ \log_2(1+\frac{1}{\mathrm{SINR}_{0k}^{-1} }) \} \geq  \log_2(1+\frac{1}{\mathbb{E}\{ \mathrm{SINR}_{0k}^{-1} \}})$ to move the expectation inside the logarithm. The expectation of the inverse SINR is then expanded as in at the top of the page.

To proceed further, we focus on the BSs in a circular area of finite radius $r$ and consider wrap around in the radial domain to keep the translation invariance. Note that by letting $r \rightarrow \infty$ a PPP in the whole of $\mathbb{R}^2$ is obtained, and there is no bias since the intra-cell distances are assumed to be distributed exactly as in $\mathbb{R}^2$. The reason for doing this is that the average number of interfering BSs in the limited circular area is finite and given by $\pi (r^2 - \mathbb{E}\{ d_{00k}^2 \} ) = \pi (r^2 - \frac{1}{\pi \lambda} )$, where we have used the fact that $d_{00k} \sim \mathrm{Rayleigh} \big( \frac{1}{\sqrt{2\pi \lambda}} \big)$ (see Lemma \[lemma:distance-distribution\]) to obtain $$\label{eq:dook_distance}
\mathbb{E}\{ d_{00k}^\nu \} =  \frac{\Gamma(\nu/2+1)}{(\pi \lambda)^{\nu/2}}$$ for any $\nu > - 2$. Two of the expectations in contain summations over the PPP $\Psi_{\lambda}$, but with different powers of the terms. To compute both, we let $\kappa=1$ or $\kappa=2$ and notice that $$\begin{aligned}
 \notag
&\mathbb{E} \left\{ {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}  
\left( \frac{ d_{jjk}^{\alpha} }{ d_{0jk}^{\alpha}} \right)^{\kappa} \right\}  \stackrel{(a)}{=} \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right)  \mathbb{E} \left\{ \left( \frac{ d_{jjk}^{\alpha} }{ d_{0jk}^{\alpha}} \right)^{\kappa} \right\}  \\ & \notag \stackrel{(b)}{=} \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right)   \mathbb{E}\left\{  \mathbb{E}\left\{ 
d_{0jk}^{-\kappa \alpha} 
  \Big| d_{jjk} \right\} d_{jjk}^{\kappa \alpha} \right\} \\ \notag
&\stackrel{(c)}{=} \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right)  \mathbb{E}\left\{ 
d_{jjk}^{\kappa \alpha} \int_{d_{jjk}}^{r} x^{-\kappa \alpha} \frac{2x}{r^2-d_{jjk}^2} dx  \right\} \\ &
\stackrel{(d)}{=} \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right) \mathbb{E}\left\{ d_{jjk}^{\kappa \alpha} \frac{2}{r^2-d_{jjk}^2} 
\frac{d_{jjk}^{2-\kappa \alpha} - r^{2-\kappa \alpha}}{\kappa \alpha-2} \right\} \notag \\  
&  
\xrightarrow{(e): \, r \to \infty} \frac{2 \pi \lambda}{\kappa \alpha-2 } \mathbb{E}\left\{ 
d_{jjk}^{2} \right\} = \frac{2 \pi \lambda}{\kappa \alpha-2 } \frac{\Gamma(2)}{\pi \lambda} = \frac{2 }{\kappa \alpha-2} \label{eq:derivation-moment}\end{aligned}$$ where $(a)$ follows from computing the average number of BSs in the circular area and by letting $j$ be an arbitrary BS index in $\Psi_{\lambda}$ and $k$ be an arbitrary UE index (since each term in the sum has the same marginal distribution). Next, $(b)$ divides the expectation into one conditional expectation where $d_{jjk}$ is given and one outer expectation with respect to $d_{jjk}$. Step $(c)$ and $(d)$ compute the inner expectation by utilizing that the BSs of other cells are uniformly distributed in the circle of radius $r$ at distances larger than $d_{jjk}$. The limit $r \rightarrow \infty$ in $(e)$ is taken both inside and outside the expectation, which is allowed since the dominated convergence theorem is satisfied.[^7] The final expression is obtained by computing $\mathbb{E}\{ 
d_{jjk}^{2} \}$ as in and exploiting the typicality, which implies that $d_{jjk}$ and $d_{00k}$ have the same marginal distribution. Using similar techniques as in yields $$\mathbb{E} \left\{ {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } \right\} = K \frac{ 2}{\alpha-2}.$$ The only remaining term in is $$\label{eq:second-moment}
\begin{split}
& \mathbb{E} \left\{ \left(  {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} }   \right) \left(  {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right) \right\} \\ &= {\overset{K}{\underset{i=1}{\sum}}} 
\mathbb{E} \left\{ {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}
{{\underset{{l \in \Psi_{\lambda} \setminus \{ j \} }}{\sum}}} \frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right\} 
\\ &+ {\overset{K}{\underset{i=1}{\sum}}} \mathbb{E} \left\{ {{\underset{{l \in \Psi_{\lambda} }}{\sum}}} \frac{d_{lli}^{\alpha}}{ d_{0li}^{\alpha} } \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } 
\right\} 
\end{split}$$ where the equality follows from computing the expectations separately over sets that contain different BS indices. The first term in is computed as\
$$\label{1000}
\begin{split}
&{\overset{K}{\underset{i=1}{\sum}}} 
\mathbb{E} \left\{ {{\underset{{j \in \Psi_{\lambda} }}{\sum}}}
{{\underset{{l \in \Psi_{\lambda} \setminus \{ j \} }}{\sum}}} \frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} } \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right\} \\ & = K 
\left( \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right) \right) \left( \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right) -1 \right)  \\ & \quad \times
\mathbb{E}\left\{ d_{jji}^{\kappa \alpha} \frac{2}{r^2-d_{jji}^2} 
\frac{d_{jji}^{2-\kappa \alpha} - r^{2-\kappa \alpha}}{\kappa \alpha-2} \right\} \\ & \quad \times
\mathbb{E}\left\{ d_{llk}^{\kappa \alpha} \frac{2}{r^2-d_{llk}^2} 
\frac{d_{llk}^{2-\kappa \alpha} - r^{2-\kappa \alpha}}{\kappa \alpha-2} \right\}  \\ 
& \xrightarrow{r \to \infty}  
 K \pi^2 \lambda^2 \! \left( \frac{2}{\alpha-2}\right)^2 \mathbb{E}\left\{ d_{jji}^{2}  \right\} 
\mathbb{E}\left\{ d_{llk}^{2} \right\} \\& = K \pi^2 \lambda^2 \! \left( \frac{2}{\alpha-2}\right)^2 \! \left(\frac{\Gamma(2)}{\pi \lambda}\right)^2 \!= K \! \left( \frac{2}{\alpha-2}\right)^2
\end{split}$$ by first computing the expectation with respect to the number of terms in the PPP for a finite radius $r$ and exploiting that $j \neq l$. By letting $r \rightarrow \infty$ we arrived at an expression that is basically twice the expression previously given in . The second term in is bounded as $$\begin{aligned}
 \notag
 & {\overset{K}{\underset{i=1}{\sum}}} \mathbb{E} \left\{ {{\underset{{l \in \Psi_{\lambda} }}{\sum}}} \frac{d_{lli}^{\alpha}}{ d_{0li}^{\alpha} } \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right\} \! \stackrel{(a)}{=} \! K \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right)  \notag
 \mathbb{E} \left\{ \frac{d_{lli}^{\alpha}}{ d_{0li}^{\alpha} } \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right\} \\ &\notag \stackrel{(b)}{\leq} K \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right)  \sqrt{
 \mathbb{E} \left\{ \frac{d_{lli}^{2\alpha}}{ d_{0li}^{2\alpha} } \right\}   \mathbb{E} \left\{ \frac{d_{llk}^{2\alpha}}{ d_{0lk}^{2\alpha} } \right\} } \\ & \stackrel{(c)}{=}  K \lambda \pi \left(r^2 - \frac{1}{\pi \lambda} \right) \mathbb{E} \left\{ \frac{d_{llk}^{2\alpha}}{ d_{0lk}^{2\alpha} } \right\} \notag \\ & \xrightarrow{(d): \, r \to \infty} 
K \frac{2 }{2 \alpha-2 } = K \frac{1 }{ \alpha-1 } \label{1001}\end{aligned}$$ where $(a)$ computes the expectation with respect to the number of terms in the PPP for a finite $r$, $(b)$ follows from Hölder’s inequality, and $(c)$ exploits that the two expectations are equal. At this point, the expression is basically the same as the second expression in , thus $(d)$ follows from the same steps as those taken in . Plugging and into we eventually obtain $$\label{eq:second-moment2}
\begin{split}
& \mathbb{E} \left\{ \left(  {{\underset{{j \in \Psi_{\lambda}}}{\sum}}} {\overset{K}{\underset{i=1}{\sum}}} 
\frac{d_{jji}^{\alpha}}{ d_{0ji}^{\alpha} }   \right) \left(  {{\underset{{l \in \Psi_{\lambda} }}{\sum}}}   \frac{d_{llk}^{\alpha}}{ d_{0lk}^{\alpha} } \right) \right\} \\ & \leq  K \left( \frac{4}{(\alpha-2)^2 } +
\frac{1 }{ \alpha-1} \right).
\end{split}$$ We have now computed all the expectations in , exactly or as upper bounds. Using all these expressions eventually lead to the achievable lower bound in .

[^1]: E. Björnson is with the Department of Electrical Engineering (ISY), Linköping University, Linköping, Sweden (emil.bjornson@liu.se). L. Sanguinetti is with the University of Pisa, Dipartimento di Ingegneria dell’Informazione, Pisa, Italy (luca.sanguinetti@iet.unipi.it) and is also with the Large Systems and Networks Group (LANEAS), CentraleSupélec, Université Paris-Saclay, 3 rue Joliot-Curie, 91192 Gif-sur-Yvette, France. M. Kountouris is with the Mathematical and Algorithmic Sciences Lab, France Research Center, Huawei Technologies Co. Ltd., France (marios.kountouris@huawei.com). E. Björnson is supported by ELLIIT and an Ingvar Carlsson Award from the Swedish Foundation for Strategic Research. L. Sanguinetti is funded by the People Programme (Marie Curie Actions) FP7 PIEF-GA-2012-330731 Dense4Green and is also supported by the ERC Starting Grant 305123 MORE. Preliminary versions [@Bjornson2015icc; @Bjornson2015spawc] of this paper were presented at the IEEE International Conference on Communications (ICC) 2015 and at the IEEE Workshop on Signal Processing Advances in Wireless Communications (SPAWC) 2015.

[^2]: Some prior works have considered erroneous EE metrics in bit/J/Hz, by forgetting to scale both the numerator and the denominator with the symbol rate (or the bandwidth) that was explicitly (or implicitly) assumed in order to compute the noise power and the dimensionless SNR value. These works are not mentioned in this section for obvious reasons.

[^3]: Alternatively, the UEs could have been distributed according to an independent homogeneous PPP [@Li2014a], leading to independent BS and UE locations. This is less sensible since smaller cells would on average have fewer UEs than larger cells, which contradicts the main principle of densifying networks mainly at places with high user loads.

[^4]: We show later that the transmit power is low at EE-optimized operating points. To simplify the notation, we have therefore not specified any maximal transmit power constraints at the users.

[^5]: The EVM is usually defined as the ratio between average distortion magnitude and signal magnitude, which becomes $\epsilon/\sqrt{1-\epsilon^2} \approx \epsilon$ with our notation.

[^6]: This is an upper bound since we only consider an average of 1000 closest interfering BSs, while the exact result requires an infinite number of interferers in $\mathbb{R}^2$.

[^7]: The dominated convergence theorem can be applied since $|d_{jjk}^{\alpha} /d_{0jk}^{\alpha}| \leq 1$.
