Channel Hardening of IRS-Aided Multi-Antenna Systems: How Should IRSs Scale?

Ali Bereyhi, Member, IEEE,
Saba Asaad, Member, IEEE,
Chongjun Ouyang, Student Member, IEEE,
Ralf R. Müller, Senior Member, IEEE,
Rafael F. Schaefer, Senior Member, IEEE and
H. Vincent Poor, Life Fellow, IEEE

Saba Asaad, Ali Bereyhi and Ralf Müller are with the Friedrich-Alexander Universität; emails: {saba.asaad<EMAIL_ADDRESS>Chongjun Ouyang is with the Beijing University of Posts and Telecommunications; email<EMAIL_ADDRESS>H. Vincent Poor is with the Princeton University; email<EMAIL_ADDRESS>This work was supported by the German Research Foundation, Deutsche Forschungsgemeinschaft (DFG), under Grant No. MU 3735/7-1, and in part by the German Federal Ministry of Education and Research (BMBF) under the Grant 16KIS1242.

$\;$Unlike active array antennas, intelligent reflecting surfaces (IRSs) are efficiently implemented at large dimensions. This allows for traceable realizations of large-scale IRS-aided MIMO systems in which not necessarily the array antennas, but the passive IRSs are large. It is widely believed that large IRS-aided MIMO settings maintain the fundamental features of massive MIMO systems, and hence they are the implementationally feasible technology for establishing the performance of large-scale MIMO settings. This work gives a rigorous proof to this belief. We show that using a large passive IRS, the end-to-end MIMO channel between the transmitter and the receiver always hardens, even if the IRS elements are strongly correlated.

For the fading direct and reflection links between the transmitter and the receiver, our derivations demonstrate that as the number of IRS elements grows large, the capacity of end-to-end channel converges in distribution to a real-valued Gaussian random variable whose variance goes to zero. The order of this drop depends on how the physical dimensions of the IRS grow. We derive this order explicitly. Numerical experiments depict that the analytical asymptotic distribution almost perfectly matches the histogram of the capacity, even in practical scenarios.

As a sample application of the results, we use the asymptotic characterization to study the dimensional trade-off between the transmitter and the IRS. The result is intuitive: For a given target performance, the larger the IRS is, the less transmit antennas are required to achieve the target. For an arbitrary ergodic and outage performance, we characterize this trade-off analytically. Our investigations demonstrate that using a practical IRS size, the target performance can be achieved with significantly small end-to-end MIMO dimensions.

Reconfigurable intelligent surface, asymptotic channel hardening, large-system analysis, ergodic capacity, outage probability.

§ INTRODUCTION

Channel hardening is a property of massive mimo systems indicating that the random fading process in a mimo channel becomes a deterministic effective path-loss, as the  dimensions of the mimo channel grow large at least at one side; see for instance <cit.>. From the information-theoretic point of view, this is the key property of massive mimo systems which leads to their significant performance gains in various aspects. These gains are rather well-studied in the literature; e.g., <cit.>.

Despite the promising performance, massive mimo systems are still considered to be implementationally intractable. This follows from the fact that deploying large active antenna arrays imposes high hardware cost and energy consumption, as well as unrealistic design requirements <cit.>. The issue becomes more crucial in practice, since several recent studies show that the fundamental features of massive mimo systems[These features are mainly described by the channel hardening.] do not hold with practical dimensions in several settings, e.g.,  mimo systems with highly-correlated channels, scenarios with small scattering angular spread and cell-free networks; some detailed studies can be found in <cit.> and <cit.>.

irs-aided mimo communication has been recently proposed as an implementationally-tractable solution to realize large-dimensional mimo systems <cit.>. An irs is a two-dimensional surface, with a large number of reconfigurable passive reflective elements. These elements are capable of tuning the wireless propagation dynamically to achieve various design objectives: For instance, in indoor millimeter-wave applications, by installing irs out of the transmission site, e.g., on interior walls or room ceilings, the blockage issue of the millimeter-wave communication can be reduced significantly, due to the link established between the transmitter and the users through the irs <cit.>. As other examples of design targets which are achieved by employing irs in the system, we can name spatial interference suppression, accurate three-dimensional beamforming, and providing a more favorable propagation environment <cit.>.

As mentioned, unlike classical large-scale mimo technologies which require large active arrays of antennas, large irs-aided mimo systems, i.e., mimo settings with large reflecting surfaces, are implementationally feasible. This follows the fact that irs passively reflect the incoming signals and do not employ radio frequency chains. As a result, the hardware cost and energy consumption at large scale decrease notably <cit.>. Inspired by this appealing advantages, irs are employed to a develop practically-tractable alternative designs for several large-scale mimo technologies, such as massive mimo <cit.>, cognitive radio <cit.>, noma <cit.> and simultaneous wireless information and power transform systems <cit.>. The key motivation behind these designs comes from this widely-accepted belief: The key features of a large-scale mimo technology can be achieved via a small-dimension mimo setting assisted by large irs. Although under idealistic assumptions this belief seems to be intuitive, it is not easy to conclude its validity in practical irs-aided settings with highly-correlated irs elements and multiple propagation links between the transmitter and the receiver. This work aims to give a rigorous proof for this believe by considering the fundamental channel hardening property in a practical irs-aided multi-antenna setting.

§.§ Main Objective and Contributions

This work gives an answer to this intriguing question: How does the end-to-end channel in irs-aided multi-antenna systems harden, when only the irs size grows large? We answer this question by deriving the asymptotic distribution of the channel capacity expression in an irs-aided fading miso system[It is worth mentioning that the results readily extend to cases with multiple receive antennas. A miso setting is assumed mainly to keep the derivations tractable. This assumption however impacts neither the analytical approach nor the final conclusions of this study.] whose number of transmit antennas is fixed and rather small, when the number of irs elements grow asymptotically large. Our derivations demonstrate that for any irs covariance matrix, the channel capacity converges in distribution to a  real-valued Gaussian distribution whose mean increases unboundedly large, and whose variance tends to zero, as the irs size grows asymptotically large.

The analytical derivations of this study gives the following answer to our target question: Even with a strongly-correlated irs, the end-to-end channel hardens regardless of the transmit array size[In fact, it hardens asymptotically with only a single transmit antenna.]. Interestingly, this property is achieved without any need for frequent tuning of irs elements.

In addition to the main results, this study has several other contributions which are briefly highlighted below:

* We characterize the distribution of the channel capacity for a generic irs-aided fading miso setting in which both direct and reflecting links are available between the transmitter and receiver. We further take into account the impact of line-of-sight channels and consider an arbitrary scaling of the irs area and a general correlation among the reflecting elements. To the best of our knowledge, the characterization of channel hardening for such setting has been left open in the literature; see for instance the recent studies in <cit.> and <cit.>.

* For a given irs and transmit array architecture, we determine the exact distribution of the end-to-end snr for an arbitrary choice of phase-shifts at the irs. We give a universal upper-bound for the mean and variance of this distribution, and propose an effective choice of irs phase-shifts that guarantees the end-to-end channel hardening.

* We validate our analytical derivations through numerical experiments. Our investigations show that the asymptotic distribution fits almost perfectly the histogram of channel capacity for practical irs sizes[And even much smaller sizes of irs.]. We further conduct several numerical experiments to investigate channel hardening for classical and extreme irs correlation scenarios.

* To characterize the speed of channel hardening, we derive a lower-bound on the mean and a uniform upper-bound on the variance of the channel capacity. Invoking the bounds, we show that for any covariance matrix at the irs whose largest eigenvalue grows sub-linearly[The definition of sub-linear growth becomes clear in the forthcoming sections.] with the irs size, the limit of the variance, when the irs size grows large, is zero. Using extreme-case investigations, we demonstrate that this conclusion is further extended to most extreme scenarios with the largest eigenvalue growing linearly and the physical dimension of irs being fixed.

* As a sample application, we use the main results to investigate the dimensional trade-off between the transmitter and the irs. The result is intuitive: By increasing the irs size, the transmit array size required to achieve a given target performance reduces. We derive the trade-off curve analytically two performance metrics, namely the ergodic capacity and outage probability.

§.§ Related Work

Channel hardening in irs-aided mimo systems is studied in the literature for some restricted settings: In <cit.>, channel hardening is discussed in an isotropic scattering environment considering an irs-aided setting with a single-antenna transmitter and receiver, i.e., a siso setting. The analysis is extended to scenarios with multiple irs and uncorrelated Rician fading in <cit.>. The concept of channel hardening is further addressed in an alternative way in <cit.>, where the authors show that in an irs-aided siso system with generic Nakagami-$m$ fading channels, the end-to-end channel between the transmitter and the receiver, in the presence of a direct link, becomes nearly deterministic. The results of this study further demonstrate that increasing the number of irs elements, as well as reducing the fading severity, enhances the channel hardening property.

In addition to classical point-to-point settings, channel hardening has been further discussed in alternative scenarios. For instance, channel hardening is investigated in <cit.> for an irs-aided siso noma system, where the irs is assumed to be realized by the intelligent omni-surfaces technology, and the reflecting elements are considered to be correlated in general. The results show that for such settings, the average achievable rate converges asymptotically to that of uncorrelated channels, when the number of irs elements goes to infinity. Another example is a mimo setting with a fully stand-alone irs-based transmitter. For such settings, the authors of <cit.> show that by growing the size of the transmitter unit, the user channels become deterministic and mutually orthogonal.

Despite connections to the above literature, the most related lines of work to our study are those given in <cit.> and <cit.>: Starting with the study in <cit.>, the authors show that the conventional form of channel hardening is not valid for an irs-aided mimo setting. The result depicts that by increasing the number of transmit antennas, while keeping the number of irs elements fixed, the hardening property does not hold. This finding follows from the fact that the channel between the irs and transmitter is identical for all receivers. As a result, by growth of transmit array size, the randomness of the fading process grows large, and hence the effective end-to-end channel does not converge to a deterministic channel.  In their following work, i.e.,  <cit.>, the authors invoke the Lindeberg-Feller central theorem <cit.> to demonstrate that the channel hardening in irs-aided systems occurs when the number of irs elements grows large[The difference of <cit.> to this work is illustrated shortly after.].

The study in <cit.> derives the distribution of the input-output mutual information for an irs-aided mimo setting in which the communication is carried out only through the irs, i.e., there is no direct link between the transmitter and receiver. The analysis invokes random matrix theory to characterize the setting in an asymptotic regime, in which both the numbers of irs elements and transmit antennas grow large with a fixed-ratio. The authors further utilize the results to derive a closed-form expression for the outage probability and the optimal values of irs phase-shifts.

Although the studies in <cit.> extend earlier results on channel hardening to a wider scope of settings, they consider several simplifying assumptions which impact the validity of the final conclusions in practical scenarios. For instance, they consider the channel coefficients of reflecting elements to be uncorrelated and ignore the line-of-sight between the irs and other terminals in the network. These assumptions are rather unrealistic, as in many use-cases of irs-aided systems, the reflecting elements are assumed to be closely packed on the irs, and the end-nodes are located in the line-of-sight of the irs <cit.>. Moreover, <cit.> assumes the complete blockage of the direct path between the transmitter and the receiver which despite its validity in some use-cases, restricts the applicability of the results and makes the comparison between a large irs-aided mimo system and a massive mimo system unfair. The asymptotic regime considered in <cit.>  is further impractical, as we are often interested in irs-aided mimo settings with large irs dimensions, but rather small arrays at the transmitter and receiver sides. In this work, we deviate from these simplifying assumptions and characterize the asymptotic channel hardening principle for a generic irs-aided scenario.

§.§ Notation and Organization

Scalars, vectors and matrices are shown by non-bold, bold lower-case, and bold upper-case letters, respectively. The notation $\mathbf{H}^{\mathsf{H}}$ indicates the transposed conjugate of $\mathbf{H}$. An $N\times N$ identity matrix is denoted by $\mathbf{I}_N$ and $\boldsymbol{1}_N$ is an $N\times N$ matrix of all-ones. We use the notation $\dbc{\mathbf{H}}_{nm}$ to refer to the entry of $\mathbf{H}$ at the $n$-th row and $m$-th column. The function $\mathrm{Q}\brc{x}$ is the standard $\mathrm{Q}$-function, i.e.,

\begin{align}
\mathrm{Q} \brc{x} = \frac{1}{\sqrt{2\pi}}\int_{x}^{+\infty} \mathrm{e}^{-\tfrac{u^2}{2}} \mathrm{d} u.
\end{align}

The mathematical expectation is denoted by $\Ex{\cdot}$, and the notation $\mathcal{CN}\left( \eta,\sigma^2\right) $ represents the complex Gaussian distribution with mean $\eta$ and variance $\sigma^2$.

The rest of the manuscript is organized as follows: Section <ref> describes the system model and formulates the problem. The main analytical results along with several numerical investigations are then presented in Section <ref>. The dimensional trade-off between the transmitter and the irs is investigated in Section <ref>. Section <ref> provides the detailed derivations of the main results. Finally, Section <ref> concludes the manuscript.

§ PROBLEM FORMULATION

Consider a miso system in which a bs with an array-antenna of size $M$ transmits data to a single-antenna user. An irs with $N$ reflecting elements is further employed to modify the propagation environment between the bs and the user: Each element of the irs reflects its received signal after applying a phase-shift. The signal received by the user is then the superposition of two components: One that is received through the direct path between the bs and the user, and one that is reflected by the irs.

§.§ System Model

As suggested by the literature <cit.>, we assume that the bs and the irs are equipped with planar arrays. Namely, we assume the antenna-array at the bs and the array of reflecting elements on the irs are rectangular planar arrays with $M_\xx$ and $N_\xx$ horizontal elements and $M_\yy$ and $N_\yy$ vertical elements, respectively, such that $M=M_\xx M_\yy$ and $N=N_\xx N_\yy$. The antenna elements at the bs are assumed to be distanced with $\ell_\xx$ and $\ell_\yy$ on the horizontal and vertical axes, respectively. The horizontal and vertical distances at the irs are further denoted by $d_\xx$ and $d_\yy$, respectively.

To model the direct path[Note that the direct path is different from the line of sight.], we consider a classical scenario in which the los link is blocked via obstacles or mobility of the user. It is further assumed that $\ell_\xx$ and $\ell_\yy$ are set large enough, such that the spatial correlation can be ignored[This is a rational assumption, as we do not assume $M$ to be large.]. We therefore adopt the Rayleigh fading model in which the channel coefficient of the direct path between antenna element $m$ at the bs and the user is modeled as

\begin{align}
h_{\mathrm{d},m}=\sqrt{\alpha_\mathrm{d} A_M } \tilde{h}_{\mathrm{d},m}.
\end{align}

Here, $\alpha_\mathrm{d}$ is a distance-dependent path-loss, $A_M$ denotes the area of a single element on the planar antenna array at the bs, i.e.,  $A_M = \ell_\xx \ell_\yy$, and $\tilde{h}_{\mathrm{d},m}$ follows a circularly-symmetric complex Gaussian distribution with zero mean and unit variance, i.e., $\tilde{h}_{\mathrm{d},m}\sim \mathcal{CN}\left(0, 1\right)$. For sake of brevity, we define the direct channel vector as $\mathbf{h}_\mathrm{d}=[h_{\mathrm{d},1}, \cdots, h_{\mathrm{d},M}]^\trp$.

Due to the flexibility in deploying the irs, it is reasonable to assume that the link between the bs and the irs is dominated by a los component. We denote the channel spanning from the bs to the irs with ${\mathbf{T}}\in \mathbb{C}^{N \times M}$  where $t_{nm}=[{\mathbf{T}}]_{nm}$
denotes the channel coefficient between the $m$-th transmit antenna  and the $n$-th irs element. In particular, $t_{mn}$ is given by

\begin{align}
t_{nm}=\sqrt{\alpha_{\mathrm{s}} A_N }\bar{t}_{nm},
\end{align}

where $\alpha_{\mathrm{s}}$ represents the distance-dependent path-loss, $A_N$ is the area of a single reflecting element, i.e., $A_N = d_\xx d_\yy$, and $\bar{t}_{nm}$ denotes the los components.

Unlike the direct and the bs-to-irs links, the link between the user and the irs has often both los and nlos components. Moreover, the reflecting elements at the irs are often closely distanced, such that the given area of the irs is filled with a large number of elements. As a result, the spatial correlation in this channel cannot be ignored. We hence use the Rician fading model with spatial correlation to model this channel: Let $h_{\mathrm{r}, n}$ denote the coefficient of the channel between the $n$-th reflecting element and the user. This coefficient is modeled as

\begin{align} \label{eq:hr}
h_{\mathrm{r}, n}=\sqrt{\alpha_{\mathrm{ r}} A_N }\left(\sqrt{\frac{\kappa_\mathrm{r}}{\kappa_\mathrm{r}+1}}\bar{h}_{\mathrm{r}, n} +\sqrt{\frac{1}{\kappa_\mathrm{r}+1}}\tilde{h}_{\mathrm{r}, n}\right),
\end{align}

where $\alpha_{\mathrm{ r}}$ and $\kappa_\mathrm{r}$ are  the distance-dependent path-loss and the Rician factor[In general,  larger $\kappa_\mathrm{r}$ means that the channel is more deterministic.], respectively. The coefficient $\bar{h}_{\mathrm{r}, n}$ represents the los, and $\tilde{h}_{\mathrm{r}, n}$ models the small-scale fading process in the nlos link. To capture the spatial correlation among the irs elements, we consider a general covariance matrix $\mR\in \setC^{N\times N}$ for the vector $\tilde{\mathbf{h}}_{\mathrm{r}}=[\tilde{h}_{\mathrm{r}, 1}, \cdots, \tilde{h}_{\mathrm{r}, N}]^\trp$ which depends on $d_\xx$ and $d_\yy$, distribution of scatterers, and the radiation pattern. For the case with isotropic scattering in front of the irs, one can determine $\mR$ explicitly from <cit.>. Consequently, $\tilde{\mathbf{h}}_{\mathrm{r}}$ is modeled as a zero-mean complex Gaussian vector with covariance $\mR$, i.e., $\tilde{\mathbf{h}}_{\mathrm{r}} \sim \mathcal{CN}\left(\boldsymbol{0}, \mR\right)$. Note that due to power normalization, we assume that $[\mR]_{nn} =1$ for $n\in[N]$.

In general, the effective area of a single element on an array depends on the wave-length, and hence the given expressions for the array areas, i.e., $A_N = d_\xx d_\yy$ and $A_M=\ell_\xx \ell_\yy$, are not exact if the elements are distances longer than a half wave-length. Nevertheless, we assume the neighboring elements on the transmit array and the irs to be distanced less than a half wave-length, i.e., $d_\xx,d_\yy,\ell_\xx,\ell_\yy \leq \lambda/2$.

The los components are described by the array responses: Let $\lambda$ be the wavelength and define the following operators

\begin{align}
i_M\brc{m} &=\brc{m-1} \;  \mathrm{mod} \; M_\xx,  \\
i_N\brc{n} &= \brc{n-1}  \;  \mathrm{mod} \;  N_\xx,
\end{align}

where $x \;  \mathrm{mod} \;  L$ determines $x$ modulo $L$. Moreover, let

\begin{align}
j_M\brc{m} &= \left\lfloor \displaystyle \frac{m-1}{M_\xx} \right\rfloor , \\
j_N\brc{n} &= \left\lfloor \displaystyle \frac{n-1}{N_\xx} \right\rfloor.
\end{align}

We now define an exponent function at azimuth angle $\varphi$ and elevation angle $\theta$ for a given element $m$ on the bs array and element $n$ on the irs array, respectively, as follows:

\begin{align}
\Phi_m\left( \varphi,\theta\right) &= {i_M\brc{m} \ell_\xx \cos{\theta} \sin{\varphi}  +  j_M\brc{m} \ell_\yy \sin{\theta} },\\
\Pi_n\left( \varphi,\theta\right) &={i_N\brc{n} d_\xx \cos{\theta} \sin{\varphi} +  j_N\brc{n} d_\yy \sin{\theta} },
\end{align}

Consequently, the transmit and irs array responses are given at $\brc{\varphi,\theta}$ respectively by <cit.>

\begin{align}
\mathbf{a}_M\left( \varphi,\theta\right) &= \left[ \mathrm{e}^{   \tfrac{2\pi\jj }{\lambda} \Phi_1\left( \varphi,\theta\right)  }, \cdots, \mathrm{e}^{  \tfrac{2\pi\jj }{\lambda} \Phi_M\left( \varphi,\theta\right) } \right]^\trp,\\
\mathbf{a}_N\left( \varphi, \theta\right) &= \left[ \mathrm{e}^{ \tfrac{2\pi\jj }{\lambda}  \Pi_1\left( \varphi,\theta\right)  }, \cdots, \mathrm{e}^{  \tfrac{2\pi\jj }{\lambda}  \Pi_N\left( \varphi,\theta\right) } \right]^\trp.
\end{align}

Let $\bar{\mathbf{T}} \in \setC^{N\times M}$ be the matrix of los channel coefficients between the bs and the irs whose entry $\brc{n,m}$ is $\bar{t}_{nm}$. Define further the los channel vector between the irs and the user as $\mathbf{\bar{h}}_{\mathrm{r}}=[\bar{h}_{\mathrm{r}, 1}, \cdots, \bar{h}_{\mathrm{r}, N}]^\trp$. We can hence write

\begin{align}
\bar{\mathbf{T}} &= \mathbf{a}_N\left( \varphi_{\mathrm{r}1},\theta_{\mathrm{r}1}\right) \mathbf{a}_M\left( \varphi_{\mathrm{t}2} , \theta_{\mathrm{t}2} \right)^\mathsf{H}, \label{a1}\\
\mathbf{\bar{h}}_{\mathrm{r}} &= \mathbf{a}_N\left( \varphi_{\mathrm{t}1} , \theta_{\mathrm{t}1} \right), \label{a2}
\end{align}

where $\brc{\varphi_{\mathrm{r}1},\theta_{\mathrm{r}1}}$ is the aoa at the irs, the pair $\brc{\varphi_{\mathrm{t}1} , \theta_{\mathrm{t}1} }$ denotes the aod from the irs, and $\brc{\varphi_{\mathrm{t}2} , \theta_{\mathrm{t}2} }$ is the aod from the bs.

The signal received by the user is given by

\begin{align}\label{eq:rx}
y=\sum_{m=1}^{M} h_m x_m + z,
\end{align}

where $z$ is zero-mean and unit-variance complex Gaussian noise, and $x_m$ denotes the symbol sent by the $m$-th antenna element at the bs. The transmitted symbols satisfy

\begin{align}
\sum_{m=1}^{M} \Ex{\abs{x_m}^2} \leq \rho,
\end{align}

for some total transmit power $\rho$. The coefficient $h_m$ in (<ref>) defines for the end-to-end effective channel between the $m$-th bs antenna and the user which is given by

\begin{align}
h_m%&=h_{\mathrm{d}, m}+\sum_{n=1}^{N} \phi_n h_{\mathrm{r}, n} t_{nm}\\
&=h_{\mathrm{d}, m}+\sum_{n=1}^{N} \mathrm{e}^{-\jj\beta_n} h_{\mathrm{r}, n} t_{nm},
\end{align}

where $\beta_n$ denotes the phase-shift with element $n$ at the irs.

§.§ Asymptotic Scaling of IRS

We intend to characterize the irs-aided setting in a large-system limit in which the transmit-array size $M$ remains fixed and the number of reflecting elements at the irs grows large. In general, the growth in the number of reflecting elements can impact the physical dimension of the irs in various forms. To illustrate this point, let $A_{\rm IRS}$ denote the total area of the irs. Assuming the reflecting elements to be symmetrically inserted on the surface, we can write

\begin{align}
A_{\rm IRS} =  N A_N.
\end{align}

The scaling of $A_{\rm IRS}$ in terms of $N$ is illustrated at best by considering the two most extreme scenarios.

* One extreme case is to assume that the number of irs elements grows large and the distances between neighboring elements on the irs are kept fixed and large enough, e.g., $d_\xx = d_\yy = \lambda/2$. In this case, $A_N$ is fixed, and therefore the total area of the irs grows linearly in $N$. Depending on the wave-length, this scaling can lead to unrealistic physical dimensions of irs for large choices of $N$.

* Another extreme case is when we set the total area of the irs fixed, and reduce the distance between the neighboring elements by the growth of $N$. In this case. $A_{\rm IRS}$ does not scale in $N$, and $A_N$ shrinks reverse-linearly with $N$. In practice, however, the physical dimensions of a typical reflecting element cannot be set below a certain threshold. Hence, for a very large choice of $N$, assuming this extreme case can also be unrealistic.

In practice, with larger numbers of elements on the irs, the distance among neighboring elements is reduced. This leads to a smaller effective area for a reflecting element and therewith to a higher spatial correlation. The distance can however be reduced down to a certain level. We can hence conclude that the scaling of the irs area is in general neither linear nor fixed, but of some intermediate order; see some related discussions in <cit.>. To take this point into account, we consider a basic scaling model for the irs area. Namely, we assume that the area of each irs element scales with $N$ as

\begin{align}
A_N = \frac{A_0}{N^{q}}
\end{align}

for some constant $A_0$ and $0\leq q \leq 1$. Consequently, the total area of the irs scales with $N$ as

\begin{align}
A_{\rm IRS} =  A_0 N^{1-q}.
\end{align}

This basic scaling addresses the limiting scenarios in between the two extreme cases.

* $q=0$ is the idealistic case which corresponds to the first extreme case, i.e., the case in which the distances between neighboring elements are kept fixed.

* $q=1$ addresses the latter case in which the total area of the irs is fixed.

By setting $0< q < 1$, an intermediate scaling is considered in which the irs surface grows sub-linearly by $N$.

It is worth mentioning that for $q\neq 1$, this scaling model implies that by sending $N\rightarrow \infty$, the area of irs also grows asymptotically large. This observation can lead to this conclusion that with $q\neq 1$, the far-field derivation for the irs array response is no longer valid in the large-system limit. To avoid such inconsistency, let $d_{\min}$ denotes the minimum distance between two terminals in the network, i.e., the minimum of the bs-to-user, bs-to-irs and irs-to-user distances. We assume that $d_{\min}$ is bounded uniformly from below as

\begin{align}
d_{\min} > D_0N^{{\gamma}/{2}},
\end{align}

for some $D_0$ and $\gamma > 1-q$. By considering this assumption, we guarantee that the far-field derivations are valid through the asymptotic analyses. From the implementational point of view, this is a realistic assumption, as the distances are often in orders of tens or hundreds of meters while the irs dimensions are often a finite multiple of the wave-length.

§.§ Input-Output Mutual Information

Let $\mathbf{h}=[h_1, \cdots, h_M]^\trp$ denote the end-to-end channel vector. Considering Gaussian signaling, the mutual information between the input and output of the end-to-end channel, for a given realization of $\mathbf{h}$, is given by

\begin{align}
\mathcal{I}\left( \mathbf{h}, \rho\mathbf{Q}\right) =\log_2 \left(1+ \rho \mathbf{h}^\mathsf{H} \mathbf{Q}\mathbf{h}\right),
\end{align}

where $\mathbf{Q} \in \setC^{M \times M}$ is the transmit covariance matrix. From the power constraint, we know that $\mathbf{Q}$ satisfies $\tr{\mathbf{Q}} = 1$.

In this channel, the optimal covariance matrix which maximizes the mutual information is given by mrt <cit.>. We hence focus on the optimal case which determines the capacity of the end-to-end channel. To this end, we define the maximum mutual information[We show it by $\mathcal{C}$, as it determines the end-to-end capacity.] as

\begin{align}\label{eq:I}
\mathcal{C} =\log_2 \left(1+ \rho M \Gamma \right),
\end{align}

which determines the end-to-end channel capacity. Here, $\Gamma$ is given by

\begin{align}
\Gamma = \frac{ \norm{ \mathbf{h} }^2}{M}.
\end{align}

In the sequel, we refer to $\Gamma$ as the end-to-end snr gain per transmit antenna. Due to fading, $\Gamma$ and $\mathcal{C}$ are random. Our main goal is to characterize the statistics of these random variables analytically, when $N$ is asymptotically large.

§ ASYMPTOTIC CHANNEL HARDENING

Intuitively, channel hardening in massive mimo systems refers to the following phenomenon: A fading mimo channel behaves almost deterministically, as the number of antennas at one end grows large[In general, it is enough to have a large antenna array at one end.]. This fundamental property is characterized in the seminal work <cit.> in which a tight approximation for the distribution of the input-output mutual information is determined. The result shows that the variance of the mutual information shrinks rapidly while its mean grows large.

Invoking the literature, it can be straightforwardly shown that the end-to-end channel $\mathbf{h}$ hardens, as the number of transmit antennas $M$ grows unboundedly large. We are however interested in a different asymptotic regime; namely, a scenario with finite transmit antennas but unboundedly large number of irs elements, i.e., fixed $M$ and $N\rightarrow\infty$. Our main result shows that in this alternative asymptotic regime, the channel still hardens, under a very mild constraint on the correlation among the irs elements.

Let the area of the irs scale with $N$ as $A_{\rm IRS} =  A_0 N^{1-q}$ for some fixed $A_0$ and $0\leq q <1$, and the phase-shifts be set to

\begin{align}
\beta_n^\star= \frac{2\pi }{\lambda} \brc{\Pi_n\left( \varphi_{\mathrm{r}1} ,\theta_{\mathrm{r}1} \right)  + \Pi_n\left( \varphi_{\mathrm{t}1} ,\theta_{\mathrm{t}1} \right)}.  \label{eq:beta_n_star}
\end{align}

Assume the maximum eigenvalue of the irs covariance matrix $\mR$, denoted by $\lambda_{\max}$, grows with $N$ sub-linearly, i.e.,

\begin{align}
\lim_{N\rightarrow\infty} \frac{\lambda_{\max}}{N} = 0,
\end{align}

and satisfies

\begin{align}
\lim_{N\rightarrow\infty} \frac{1}{\lambda_{\max} A_{\rm IRS}} = 0. \label{eq:const_IRS}
\end{align}

Then, the maximum mutual information $\mathcal{C}$ is well approximated by a real Gaussian random variable whose mean and standard deviation are given by

\begin{align}
\mu_{\mathcal{C}} &=   \log_2 \left(1 +\rho M \mu \right) \\
\sigma_{\mathcal{C}} &= \frac{ \displaystyle \rho M  \log_2 \mathrm{e} }{1+ \rho M \mu }  \sqrt{ \omega \eta + \eta + \frac{M-1}{M} \alpha_{\mathrm{ d}} A_M }
\end{align}

respectively, for

\begin{align}
\mu &=  \alpha_{\mathrm{ d}} A_M+ \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 + \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} ,\\
\eta &=  \frac{\alpha_{\mathrm{ d}}A_M}{M} +  \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} ,\\
\omega &=  2\kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 + \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} ,
\end{align}

and $\bar{\alpha}_N$ being defined as

\begin{align}
\bar{\alpha}_N=  \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} A_N^2  . \label{eq:bar_alpha}
\end{align}

The proof is given in Section <ref>.

Proposition <ref> considers three main constraints on the system setting; namely, it restricts the scaling order of the irs area, assumes sub-linearly growing $\lambda_{\max}$ and the constraint in (<ref>). In this respect, it is worth mentioning few remarks.

Proposition <ref> excludes the limiting case of $q=1$, for analytical rigor. Nevertheless, as we see in the forthcoming section, the analytical expressions are still rather accurate for $q=1$.

* Constraining the growth rate of the maximum eigenvalue $\lambda_{\max}$ to be sub-linear can be interpreted as follows: $\lambda_{\max}$ is uniformly bounded from above as $\lambda_{\max} \leq aN^u$ for some real scalars $a$ and $0\leq u <1$. This is in general not a strong constraint, as we can write

\begin{align}
\lambda_{\max} \leq \tr{\mR} = N.
\end{align}

As we illustrate later on, the exponent $u$ is mutually coupled with $q$. This is seen by considering the two extreme cases of $q$. For $q=0$, it is feasible to have $\mR=\mI_N$, and hence $\lambda_{\max} = 1$ meaning that $u=0$. The other extreme case is $q \rightarrow 1$ in which $\mR$ tends to a rank-one matrix. In this case $\lambda_{\max} = N$, and hence the uniform upper-bound is valid with[Similar to $q$, we exclude the linear case, i.e., $u=1$ for analytical rigor.] $a = u =1$. For less rank-deficient covariance matrices,  $u$ can be something in $\dbc{0,1}$. Similar to the constraint on the physical dimension of the irs, we show later that despite assuming $0\leq u <1$, Proposition <ref> is still valid for the extreme case of $u=1$.

* Considering the scaling order of the irs area $A_{\rm IRS}$ and maximum eigenvalue $\lambda_{\max}$, the constraint in (<ref>) restricts $u$ and $q$ to satisfy $u\geq q$.  One may find this constraint intuitively valid following from the above discussions on the coupling of $u$ and $q$. We later on show that this is in fact the case for the covariance matrix derived for the standard Rayleigh model in <cit.>.

§.§ Numerical Investigations

Before starting with derivations, let us confirm the accuracy of Proposition <ref> for practical system dimensions through some numerical experiments. To this end, we consider a basic scenario in which the bs is equipped with a $2 \times 2$ planar array and the irs surface consists of $N=256$ reflecting elements which are aligned on a rectangular surface with $N_\xx = 8$ horizontal elements and $N_\yy = 32$ vertical elements. The elements at both bs and irs are distanced with $\ell_\xx = d_\xx = \ell_\yy = d_\yy = \lambda /2$.

For sake of brevity, we set $\alpha_{\mathrm{ d}} A_M = \alpha_{\mathrm{ r}} A_N= \alpha_{\mathrm{ s}} A_N=1$, and $\log \kappa_{\mathrm{ r} } = 0$ dB. For the covariance matrix $\mR$, we consider a Rayleigh fading model and invoke <cit.>. The entry $\brc{n,n'}$ of $\mR$ is hence given by

\begin{align}
\dbc{ \mR }_{n n'} = \mathrm{sinc} \brc{\frac{2}{\lambda} \sqrt{D_\xx^2 + D_\yy^2 } } \label{eq:RR}
\end{align}

with $D_\xx$ and $D_\yy$ being

\begin{align}
D_\xx &= d_\xx \dbc{i_N\brc{n} - i_N\brc{n'}},\\
D_\yy &= d_\yy \dbc{j_N\brc{n} - j_N\brc{n'}}.
\end{align}

The aoa and aod are further set to $\brc{\varphi_{\mathrm{r}1},\theta_{\mathrm{r}1}} = \brc{\pi/6,\pi/3 }$, $\brc{\varphi_{\mathrm{t}1},\theta_{\mathrm{t}1}} = \brc{\pi/8,2\pi/3 }$ and $\brc{\varphi_{\mathrm{t}2},\theta_{\mathrm{t}2}} = \brc{\pi/7,\pi/5 }$. The power constraint is set to $\rho=1$.

For this setting, we collect $10^5$ realizations of the channel and determine the input-output mutual information for every realization. The histogram of the collected samples is shown in Fig. <ref>. As the figure shows, the histogram closely matches the Gaussian distribution. We now plot the analytical distribution given by Proposition <ref> and compare it to the Gaussian distribution fitted to the histogram. The result is shown in Fig. <ref>. As the figure demonstrates, the analytical result of Proposition <ref> matches almost perfectly to the empirical density.

Histogram of the mutual information and the fitted Gaussian distribution for the covariance matrix $\mR$ given in (<ref>).

Probability density of the maximum mutual information and the Gaussian distribution fitted to the numerical simulations.

We now conduct a new experiment. We consider the same setting and let the irs size to vary as $N = N_\xx^2$ while changing from $N_\xx = 8$ to $N_\xx = 36$, gradually. Note that in this setting, the area of the irs grows linearly[Such scaling is not realistic for very large irs. We discuss more realistic scaling scenarios in the forthcoming parts of this section.] in $N$, i.e., $t = 0$, and hence the limit in (<ref>) is satisfied. For each choice of $N$, we collect $500$ samples, determine numerically the mean and variance of $\mathcal{C}$ and plot it against $N$ in Figs. <ref> and <ref>.

The numerical results in Figs. <ref> and <ref> are compared with the closed-form expressions given in Proposition <ref>. As the figures show, the analytical expressions closely track the numerical results. The figures further demonstrate the hardening of the channel in terms of the irs size. As $N$ grows large, the mean mutual information grow large, while the variance drops.

§.§ Channel Hardening Order

The numerical experiments in the previous section consider the case of $q=0$, i.e., when the irs area grows linearly in $N$. Although this can be considered feasible for a small or moderate number of reflecting elements, it is not a realistic scaling for asymptotically large irs. In fact, in such scenarios, the total area of the irs is limited, and hence, the area grows sub-linearly, i.e., $q > 0$. We address this point by investigating the speed of channel hardening in this section.

We are mainly interested to find out how fast the channel hardens with respect to $N$ for a given scaling of the irs area. In general, the speed depends on the correlation among the reflecting elements, and therefore the physical dimensions of the irs. This argument is analytically characterized in Proposition <ref>.

Numerical simulations and the analytic expression for the mean mutual information against $N$.

Numerical simulations and the analytic expression for the variance of the mutual information against $N$.

Let the sub-linearly growing maximum eigenvalue of $\mR$ be uniformly bounded from above as

\begin{align}
\lambda_{\max} \leq a N^u,
\end{align}

for some non-negative real $a$ and $0 \leq u < 1$. Then, there exists an integer $N_0$, such that for $N \geq N_0$ the mean and variance of the maximum mutual information is bounded uniformly as

\begin{align}
\mu_{\mathcal{C}} &\geq b + \brc{1-q}  \log_2 N\\
\sigma_{\mathcal{C}}^2 &\leq {\frac{c}{N^{1-u}}},
\end{align}

for some non-negative real scalars $b$ and $c$.

We start the proof by noting that

\begin{align}
\bar{\alpha}_N &=  \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} A_N^2\\
&= \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} \frac{A_{\rm IRS}^2}{N^2}
= \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} \frac{A_0^2}{N^{2q}}. \label{eq:aplpha_N}
\end{align}

From Proposition <ref>, since $\mu \geq \bar{\alpha}_N\kappa_{\mathrm{ r} }N^2$, we have

\begin{align}
\mu_{\mathcal{C}} %&=   \log_2 \left(1 +\rho M \eta + \rho M  \bar{\alpha}_N\kappa_{\mathrm{ r} }N^2 \right)\\
&\geq   \log_2 \left( \rho M \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 \right).
% &= \log_2 \left( \rho M  \kappa_{\mathrm{ r} } \right) + \log_2 \bar{\alpha}_N N^2.
\end{align}

Using  (<ref>), we can hence conclude that

\begin{align}
\mu_{\mathcal{C}} %&=   \log_2 \left(1 +\rho M \eta + \rho M  \bar{\alpha}_N\kappa_{\mathrm{ r} }N^2 \right)\\
%&\geq   \log_2 \left( \rho M \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 \right)\\
&= \log_2 \left( \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} \rho M  \kappa_{\mathrm{ r} } A_0^2  \right) + \log_2 N^{2-2q}.
\end{align}

Noting that the first term on the right hand side does not grow in $N$, we conclude that

\begin{align}
\mu_{\mathcal{C}} &\geq  b + \brc{1-q}\log_2 N,
\end{align}

by setting $b$ to

\begin{align}
b = \log_2 \left( \frac{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} }{1+\kappa_\mathrm{r}} \rho M  \kappa_{\mathrm{ r} } A_0^2  \right) .
\end{align}

To further bound the variance, we note that $\norm{\mathbf{\bar{h} }_{\mathrm{ r}}}^2 = N$. As a result, $\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}}/N$ determines the Rayleigh quotient of $\mR$ at $\mathbf{\bar{h} }_{\mathrm{ r}}$ which is bounded from above by the maximum eigenvalue of $\mR$ <cit.>. We hence have

\begin{align}
\bar{\alpha}_N	\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}}&\leq   \bar{\alpha}_N\lambda_{\max} N.
\end{align}

Let $\lambda_{\max}$ be bounded uniformly from above by $aN^u$ for some $0\leq u <1$. This concludes that $\mu$, $\eta$ and $\omega$ in Proposition <ref> are uniformly bounded from above respectively by

\begin{align}
\mu^{\rm U } \brc{N} &= \mathcal{O}\brc{ N^{2-2q} }\\
\eta^{\rm U } \brc{N} &= \mathcal{O}\brc{ N^{1+u-2q} }\\
\omega^{\rm U } \brc{N} &= \mathcal{O}\brc{ N^{2-2q} }.
\end{align}

We now start with bounding $\sigma_{\mathcal{C}}$. Since $1+ \rho M \mu >  \rho M \mu$, we can initially bound the variance as

\begin{align}
\sigma_{\mathcal{C}} &< \frac{ \displaystyle  \log_2 \mathrm{e} }{\mu }  \sqrt{ \omega \eta + \eta + \frac{M-1}{M} \alpha_{\mathrm{ d}} A_M } =\Xi_N \sqrt{\frac{ \eta  }{\mu }}
\end{align}

where $\Xi_N$ is given by

\begin{align}
\Xi_N=   \sqrt{ \frac{1}{\mu \eta}\brc{ \omega \eta + \eta + \frac{M-1}{M} \alpha_{\mathrm{ d}} A_M } } \log_2 \mathrm{e}
\end{align}

The term $\Xi_N$ is uniformly bounded from above by a constant. We can hence conclude that there exists an integer $N_1$ and a non-negative real $c_0$, such that for $N\geq N_1$, we have

\begin{align}
\sigma_{\mathcal{C}}
&\leq c_0  \sqrt{ \frac{ \eta }{ \mu  } } %\stackrel{\mathrm{(a)}}{\leq} c_1 \frac{\sqrt{\eta}}{N}
\end{align}

Noting that $\mu \geq \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2$, we can further bound the above upper-bound as

\begin{align}
\sigma_{\mathcal{C}}
&\leq c_1 \frac{\sqrt{\eta}}{N^{1-q}} \label{upper:1}
\end{align}

where $c_1$ being defined as

\begin{align}
c_1 = \sqrt{\frac{c_0 \brc{1+\kappa_{\mathrm{ r} }}}{\alpha_{\mathrm{ r}} \alpha_{\mathrm{ s}} \kappa_{\mathrm{ r} }}}.
\end{align}

Since $\eta$ is uniformly bounded with $\eta^{\rm U } \brc{N}$, we conclude that there exists an integer $N_2 \geq N_1$ and a non-negative real $c_2$, such that for $N \geq N_2$, we have

\begin{align}
\frac{ \eta }{  N^{2-2q}  } \leq \frac{c_2 }{N^{1-u}}.
\end{align}

The proof is finally completed by setting $c = c_1^2 c_2$ and $N_0$ to be $N_0 = N_2$.

Proposition <ref> formulates an intuitive behavior. The higher the correlation among the irs elements is, i.e., larger $u$, the slower the end-to-end channel hardens. It further indicates that the fastest hardening is achieved by $u=q=0$. This corresponds to reflecting arrays whose number of strongly correlated elements does not scale with $N$ while the area of the irs grows linearly in $N$. From implementation viewpoint, this means that the neighboring reflecting elements on the irs are well-distanced, and the physical dimensions of the irs grow with the number of reflecting elements, such that the distance between two neighboring elements remains constant.

We now validate Proposition <ref> through a numerical experiment. We consider the covariance matrix derived in <cit.> for a planar irs whose elements are spaced with $\lambda/2$ in both directions, i.e., $d_\xx = d_\yy = \lambda/2$. This means that for this array $q=0$. The covariance matrix in this case is specified via (<ref>). We consider a square array, i.e., $N_\xx = N_\yy$ and let $N$ grow gradually from $N= 64$ to $N=1296$. For this sequence of covariance matrices, $\lambda_{\max}$ is plotted against $N$ in Fig. <ref>. We further use the curve fitting toolbox of MATLAB, i.e.,  <cit.>, to fit the collected data to a curve of form $\lambda_{\max} = a N^u$. The result is shown with a dashed line in the figure for which $a=0.83$ and $u=0.25$.

$\lambda_{\max}$ against $N$ for the classic irs covariance matrix derived in <cit.>. The dashed line shows the fitted curve which reads $\lambda_{\max} \approx 0.83\sqrt[4]{N}$, i.e., $u = 0.25$.

Given the results in Fig. <ref>, Proposition <ref> suggests that the variance of the mutual information in this case is uniformly bounded from above by $cN^{-0.75}$ for some $c$. This is shown in Fig. <ref>, where we compare $\sigma_{\mathcal{C}}^2$ against the upper bound with $c=1.9$. The numerical results show consistency with Proposi-tion <ref>. Interestingly, the suggested upper bound gives a pessimistic approximation of the hardening speed. In fact the true variance, drops much faster than the order of the upper bound.

Comparing the channel hardening speed to the uniform upper bound suggested by Proposition <ref>.

§.§ Scaling of IRS Area and Covariance Matrix

As indicated earlier, $u$ and $q$ are mutually coupled. This is easily seen by sending $q\rightarrow 1$. In this case, the total area of the irs is constant. This means that in the limit $N\rightarrow \infty$, the covariance matrix converges to a rank-one matrix and hence

\begin{align}
\lambda_{\max} = \tr{\mR} = N,
\end{align}

or equivalently $u\rightarrow 1$; see also <cit.>. To understand this mutual coupling analytically and find out how realistic the condition (<ref>) is, i.e., $u\geq q$, we conduct a numerical experiment. For a given $q$, we repeat the previous experiment. This means that we consider square irs of size $N$, i.e., $N_\xx = N_\yy = \sqrt{N}$, and let $N$ grow gradually from $N= 64$ to $N=1296$. The distance between two neighboring elements on the irs in this case is

\begin{align}
d_\xx = d_\yy = \frac{\lambda}{2N^{q/2}}.
\end{align}

Invoking the least-squares method, we numerically find $u$ such that $\lambda_{\max} \leq aN^u$ for some real $a$. We then change $q$ from $q=0$ to $q=1$ and plot $u$ against it. The result is shown in Fig. <ref> along with the region $u \geq q$ which corresponds to the constraint given by (<ref>).

Eigenvalue order exponent $u$ against the area order exponent $q$. The gray area is the feasible region $u\geq q$ where (<ref>) is satisfied.

As the figure suggests, for the covariance matrix, derived for the Rayleigh fading model, $u\geq q$. This is intuitive, since in both extreme cases of $q =0$ and $q\rightarrow1$, this constraint is satisfied. The result further show a linear growth of $u$ in $q$. This observation further suggests that Propositions <ref> and <ref> are further valid approximations in the limits $q\rightarrow 1$ and $u\rightarrow 1$. We confirm this intuition through numerical analysis of this limiting case in the sequel.

§.§ The Extreme Case of Rank-One Covariance

As mentioned, the sub-linearity constraint in Propositions <ref> and <ref> is not a strong constraint. In fact, noting that $\lambda_{\max} \leq N$, one can conclude that the sub-linearity constraint is satisfied in all settings, but the extreme case with a finite-rank covariance matrix whose rank does not grow with $N$, and therefore $\lambda_{\max}$ scales linearly in $N$. From the mutual coupling between $u$ and $q$, in this case we further have $q=1$ which corresponds to an irs whose total area is kept fixed.

Inspired by our observation in Section <ref>, we conjecture that Propositions <ref> and <ref> give accurate characterizations, even in an extreme case with a linearly growing $\lambda_{\max}$. To confirm this conjecture, we run simulations for the same setting considered in Section <ref> while replacing the covariance matrix with $\mR = \boldsymbol{1}_N$, i.e., the matrix of all ones, and correspondingly setting $q=1$, i.e., assuming

\begin{align}
A_N = \frac{A_0}{N}.
\end{align}

for some real $A_0$. To achieve this scaling, we set the distance between two neighboring irs elements to be

\begin{align}
d_\xx = d_\yy = \frac{\lambda}{2\sqrt{N}}
\end{align}

which means that $A_0 = \lambda^2/4$. We further set $\alpha_{\mathrm{ r}} = \alpha_{\mathrm{ s}} = 1/A_0$ to make the comparison to the first experiment in Section <ref> fair. The remaining parameters are the same as those considered in Section <ref>.

The matrix $\mR = \boldsymbol{1}_N$ is a rank one matrix whose maximum eigenvalue reads $\lambda_{\max} = N$, and thus scales linearly with $N$. One can observe this case, as an extreme case of the covariance matrix in (<ref>), in which the distance between the most outer two elements of the irs is still significantly smaller than the wavelength[In other words, an irs whose physical dimensions are significantly smaller than the wavelength.].

Fig. <ref> shows the numerically-evaluated histogram of $\mathcal{C}$, as well as the properly scaled version of the fitted Gaussian probability density and the distribution suggested by Proposition <ref>. The histogram is evaluated for $10^5$ channel realizations. As the figure shows, the histogram in this case is rather loosely approximated by the distribution given in Proposition <ref>, and compared to the almost-perfect match in Figs. <ref> and <ref>, it shows some approximation error. This approximation error follows from the lower speed of channel hardening in this case. In fact, from Proposition <ref>, we know that by sending $q,u\rightarrow 1$, the mean and variance of $\mathcal{C}$ scale with $\mathcal{O}\brc{1}$ in $N$, and hence, the approximation in Proposition <ref> becomes inaccurate[This is the main reason that Proposition <ref> excludes the case $q=u=1$.].

Histogram and the fitted Gaussian density of $\mathcal{C}$ for the extreme case of $\mR = \boldsymbol{1}_N$.  The scaled version of the analytical density function given by Proposition <ref> is further plotted with a solid line.

Although the approximation by Proposition <ref> is inaccurate, the histogram in Fig. <ref> shows a rather small standard deviation for $\mathcal{C}$. This is an interesting observation confirming our earlier finding in Section <ref> which indicates that the upper bound given by Proposition <ref> is rather pessimistic.This finding can be analytically supported as follows: The sub-linearity constraint in Propositions <ref> and <ref> comes from utilizing the inequality  $	\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} \leq   \lambda_{\max} N$ to bound the variance $\sigma_{\mathcal{C}}^2$. Nevertheless, for rank-deficient covariance matrices, this bound is often asymptotically very loose. This means that despite the inaccuracy of Proposition <ref> for exactly linear scaling, i.e., $q=u=1$, the channel hardening still occurs in this case.

We now invoke the above heuristic conclusion to further modify Proposition <ref> for linearly scaling $A_{\rm IRS}$ and $\lambda_{\max}$. In fact, assuming that $\mathcal{C}$ is still Gaussian in this case, we can conclude that the approximation error in this case comes from the inaccuracy of the variance in  Proposition <ref> for cases with linear scaling. The following pseudo-proposition gives a modified approximation for the distribution of $\mathcal{C}$ which is also accurate for linearly scaling scenarios. Interestingly, by assuming sub-linear scaling, the approximation reduces to the limiting distribution given by Proposition <ref>.

For linearly scaling $A_{\rm IRS}$ and $\lambda_{\max}$, i.e., cases with $q=u=1$, the maximum mutual information $\mathcal{C}$ is well approximated by a real Gaussian random variable whose mean is $\mu_{\mathcal{C}}$ and whose variance is

\begin{align}
\hat{\sigma}^2_{\mathcal{C}} &=  \frac{\kappa_{\mathrm{ r} }}{\kappa_{\mathrm{ r} } + \vartheta }  \sigma^2_{\mathcal{C}} ,
\end{align}

where $\mu_{\mathcal{C}}$ and $\sigma^2_{\mathcal{C}}$ are given in Proposition <ref>, and

\begin{align}
\vartheta = \frac{\alpha_{\mathrm{ d}} A_M}{2 \bar{\alpha}_N N^2} \label{eq:vartheta}
\end{align}

The proof is given in Section <ref>.

For the numerical experiment of Fig. <ref>, we further plot the approximation by Pseudo-Proposition <ref> in Fig. <ref>. As the figure demonstrates, the proposed approximated density tracks the empirically-evaluated histogram more closely[The interested reader can check that for $u,q<1$, we have $\hat{\sigma}_{\mathcal{C}}^2 \rightarrow \sigma_{\mathcal{C}}^2$, as we set $N\rightarrow \infty$.] compared with the limiting result given by Proposition <ref>.

Histogram and the fitted Gaussian density of $\mathcal{C}$ for the extreme case of $\mR = \boldsymbol{1}_N$.  The scaled version of the density functions given by Pseudo-Proposition <ref> and Proposition <ref> are plotted with solid and dotted lines, respectively.

§ IRS-TRANSMITTER DIMENSIONAL TRADE-OFF

The analytical results of this study enable us to address various design challenges in irs-aided systems. In this section, we focus on a particular application: We employ the results to investigate the dimensional trade-off between the bs and the irs. More precisely, we try to answer this fundamental question: How does the transmit array dimension, i.e., $M$, change, when we employ an irs to enhance the communication link?

§.§ Dimensional Trade-Off for Ergodic Capacity

We start our investigations by considering cases in which the nlos links experience a fast fading process. This means that the channel varies from a short block of symbol transmit intervals to another. In this case, the performance is best described by the ergodic capacity, defined as

\begin{align}
\bar{\mathcal{C}} = \Ex{\mathcal{C}}.
\end{align}

For a target ergodic capacity, the bs needs to be equipped with a certain number of transmit antennas, i.e., $M$. Intuitively, this number of required antennas is expected to reduce, as we increase the irs dimension. This draw a dimensional trade-off between $N$ and $M$: With a larger irs, the minimum number of required bs antennas to achieve a given target performance decreases. Proposition <ref> enables us to quantitatively formulate this dimensional trade-off. To this end, we first use Proposition <ref> and derive the ergodic capacity for a given $N$ and $M$ in a closed form as follows:

\begin{align}
\hspace*{-1mm}	\bar{\mathcal{C}} =\log_2\brc{1+\rho M \dbc{\alpha_{\mathrm{ d}} A_M + \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 + \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}}}}. \hspace*{-1mm}
\end{align}

It is hence readily concluded that for a given $N$ and the target ergodic capacity $\bar{\mathcal{C}}$, $M$ needs to satisfy

\begin{align}
M \geq \frac{2^{\bar{\mathcal{C}} } - 1}{
\rho \brc{\alpha_{\mathrm{ d}} A_M + \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 + \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}}}
}. \label{eq:Min_M_Erg}
\end{align}

The inequality in (<ref>) formulates the trade-off between $N$ and $M$: By increasing the physical dimensions of the irs, a smaller array is required at the bs. It further specifies the speed of this drop, i.e., $M$ drops proportional to $N^{2q-2}$.

We now conduct some experiments to investigate this trade-off numerically. To this end, a setting is considered in which an irs is distanced from the bs with $D_{\rm s} = 25$ m.  The receiver is further located at a random point at which its distances from the bs and the irs are $D_{\rm d} = 20$ m and $D_{\rm d} = 15$ m, respectively. We further set the aoa and aod to be the same as those considered in Section <ref>, and assume that $\rho=1$. The large-scale fading coefficients are generated as

\begin{align}
\alpha_i = \frac{\alpha_{\rm ref}}{D_i^{\varepsilon_i}}
\end{align}

for $i\in\set{\mathrm{ s},\mathrm{ d},\mathrm{ r}}$, where $\varepsilon_i$ is the path-loss exponent of link $i$ and $\alpha_{\rm ref}$ is the reference path-loss. For numerical simulations, we set $\log \alpha_{\rm ref} = 10$ dB, $\varepsilon_{\mathrm{ d}} = 3.5$ and $\varepsilon_{\mathrm{s}} = \varepsilon_{\mathrm{r}} = 2.3$.

The irs is considered to be a square planar array, i.e., $N_\xx = N_\yy$, whose elements are distanced with a half wavelength. This means that $q=0$. The covariance matrix of the nlos link is further set to be the Rayleigh covariance matrix, i.e., (<ref>). We assume that the phase-shifts of the reflecting elements are set according to Proposition <ref>; hence, the configuration of the transmit array is not required to be taken into account.

To investigate the dimensional trade-off between $M$ and $N$, we let the ergodic capacity to be set to a given target value $\bar{\mathcal{C}}$. The number of reflecting elements is then increased gradually starting from $N=64$ and ending at $N=1296$. For each $N$, we determine the real-valued lower bound on the minimum required $M$, denoted by $M_{\rm Erg}$, from (<ref>), i.e.,

\begin{align}
M_{\rm Erg} = \frac{2^{\bar{\mathcal{C}} } - 1}{
\rho \brc{\alpha_{\mathrm{ d}} A_M + \kappa_{\mathrm{ r} } \bar{\alpha}_N N^2 + \bar{\alpha}_N \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} }
\end{align}

Fig. <ref> plots $M_{\rm Erg}$ against the irs size $N$, for three various choices of $\bar{\mathcal{C}}$. The minimum required bs array size, denoted by $M_{\min}$, is then determined by quantizing $M_{\rm Erg}$ to the next larger integer.  As the figure shows, for rather large choices of $N$, the target ergodic capacity with only one antenna at the transmitter. Our further numerical investigations confirm the consistency of the result with simulations.

IRS-transmitter dimensional trade-off for target ergodic capacity $\bar{\mathcal{C}}$: The minimum required number of transmit antennas $M_{\min}$ is determined by quantizing $M_{\rm Erg}$ to the next larger integer.

§.§ Dimensional Trade-Off for Outage Probability

With a slow fading process, the channel remains fixed within a long sequence of symbol intervals. The appropriate metric for performance evaluation in this case is therefore the outage probability which is defined for a given target rate $\mathcal{R}$ as <cit.>

\begin{align}
\mathcal{P}_\mathrm{out} \brc{ \mathcal{R}} =\Pr \left\lbrace 	\mathcal{C} < \mathcal{R}\right\rbrace .
\end{align}

The dimensional trade-off can be also studied for this metric using the analytical results of Section <ref>. From Proposition <ref>, we can determine the outage probability as

\begin{align}
{\mathcal{P}}_\mathrm{out} \brc{ \mathcal{R}}= \mathrm{Q}\left( \frac{\mu_{\mathcal{C}} -\mathcal{R}}{\sigma_{\mathcal{C}} }\right),
\end{align}

with $\mathrm{Q}\brc{\cdot}$ being the standard $\mathrm{Q}$-function. To guarantee achieving ${\mathcal{P}}_\mathrm{out} \brc{ \mathcal{R}} \leq p_{\rm out}$ for a given target rate $\mathcal{R}$, we need to have

\begin{align}
\frac{\mu_{\mathcal{C}} -\mathcal{R}}{\sigma_{\mathcal{C}} } \geq  \mathrm{Q}^{-1} \brc{p_{\rm out}},
\end{align}

where $\mathrm{Q}^{-1}\brc{\cdot}$ is the inverse of the $\mathrm{Q}$-function with respect to composition, i.e., $\mathrm{Q}^{-1}\brc{\mathrm{Q}^{-1} \brc{p_{\rm out}}} = p_{\rm out}$.

For a given reflecting array and rate $\mathcal{R}$, it is readily shown, after few lines of calculations, that the target outage is achieved if $M\geq M_{\rm Out}$, where $M_{\rm Out}$ a positive real-valued solution to the following fixed-point equation in $\xx$:

\begin{align}
\mathcal{R} + \frac{C}{1 + \rho \mu \xx} \sqrt{B \xx^2 + \omega\alpha_{\mathrm{ d}} A_M \xx} = \log_2\brc{ 1+ \rho \mu \xx}.
\end{align}

Here, $\mu$ and $\omega$ are defined for the given irs as in Proposition <ref>. Moreover, $C$ and $B$ are given by

\begin{align}
C &= \rho\mathrm{Q}^{-1} \brc{p_{\rm out}} \log_2 \mathrm{e},\\
B &= \brc{\omega+1} \mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} + \alpha_{\mathrm{ d}} A_M.
\end{align}

Fig. <ref> demonstrates the irs-bs dimensional trade-off with respect to the outage probability for various target outages. Here, the setting is set to be the same as the one considered in Fig. <ref>. The target rate is moreover set to $\mathcal{R} = 3$. For $50\%$ outage, the trade-off curve recovers the one given in Fig. <ref> for $\bar{\mathcal{C}} = 3$. This follows from the symmetry of the Gaussian density which leads to this property that the median coincides with the mean. We further observe that by reducing the outage probability, the trade-off figure shifts upward. This is intuitive, as for lower-outages, we require a better end-to-end link. These upward shifts are however not significant and rather negligible for large choices of $N$. This is a direct result of the channel hardening: With large $N$, the capacity term $\mathcal{C}$ is almost deterministic, and hence the outage probability tends to a step function in $\mathcal{R}$.

IRS-transmitter dimensional trade-off for various target outages: Here, the target rate is set to $\mathcal{R} = 3$.

§ DERIVATION OF THE MAIN RESULT

This section proves the asymptotic result given in Proposition <ref>. The derivation follows three main steps:

* First, we derive the distribution of the end-to-end snr gain $\Gamma$ for an arbitrary $M$ and $N$.

* We then bound the mean and variance of the snr gain in terms of the spectrum of the covariance matrix $\mR$.

* Finally, we send $N\rightarrow\infty$ and show that the maximum mutual information converges in distribution to a Gaussian random variable, if the maximum eigenvalue of $\mR$ grows sub-linearly in $N$.

§.§ Distribution of End-to-End SNR

Theorem <ref> determines the statistics of the end-to-end snr gain $\Gamma$ for an arbitrary choice of system dimensions.

Let the vector of phase-shifts at the irs be set to $\bbeta = [\beta_1,\ldots,\beta_N]^\trp$. Define $E\brc{\cdot}$ as

\begin{align}
E\brc{\bbeta} &=    \bar{\alpha}_N\;
\mathbf{\bar{g}}_{\mathrm{r}}^\trp \mPhi \brc{\bbeta} {\mR} \mPhi^\her \brc{\bbeta}	\mathbf{\bar{g}}_{\mathrm{r}}^*,  \label{eq:E}
\end{align}

for matrix $\mPhi \brc{\bbeta} = \mathrm{diag} \set{ \mathrm{e}^{-\jj \beta_1},\ldots,\mathrm{e}^{-\jj \beta_N} }$, $\bar{\alpha}_N$ given in (<ref>) and $\mathbf{\bar{g}}_{\mathrm{r}} = \mathbf{a}_N\left( \varphi_{\mathrm{r}1} , \theta_{\mathrm{r}1} \right)$. Let $F_\zeta\brc{\bbeta}$ and $\Lambda\brc{\bbeta}$ be

\begin{align}
F_\zeta\brc{\bbeta} &=\bar{\alpha}_N\kappa_{\mathrm{ r} }
\abs{ \mathbf{\bar{g}}_{\mathrm{r}}^\trp  \mPhi \brc{\bbeta} 	\mathbf{\bar{h} }_{\mathrm{r} } }^2 + \zeta E\brc{\bbeta} \label{eq:F} \\
\Lambda\brc{\bbeta} &= \frac{\alpha_{\mathrm{ d}} A_M }{M}  + E\brc{\bbeta}. \label{eq:Lambda}
\end{align}

Then, the end-to-end snr gain $\Gamma$ is distributed with a generalized chi-square distribution of order $2M$ whose density is given by

\begin{align}
f_\Gamma\brc{\gamma}= \frac{\tau_M}{ \Lambda\brc{\bbeta} } \int_{-\infty}^{+\infty} f_1\brc{\frac{\gamma - v}{ \displaystyle  \Lambda\brc{\bbeta} } } f_0 \brc{\tau_M v } \mathrm{d} v,
\end{align}

for $\tau_M ={M}/{\alpha_{\mathrm{ d}} A_M }$ and the density functions

\begin{align}
f_0\brc{v}&=\frac{v^{M-2} }{\brc{M-2}!} \mathrm{e}^{-v}, \label{eq:f_0}\\
f_1\brc{v}&= \mathrm{e}^{-v - \lambda\brc{\bbeta}  } I_0 \brc{2 \sqrt{ \lambda\brc{\bbeta} v } } , \label{eq:f_1}
\end{align}

with $\lambda\brc{\bbeta} = {F_0\brc{\bbeta}} / {\Lambda\brc{\bbeta}}$, and denoting the modified Bessel function of the first kind and order zero.

Moreover, the mean and variance of $\Gamma$ are given by

\begin{align}
\mu_\Gamma \brc{\bbeta}  &= \alpha_{\mathrm{ d}} A_M  + F_1\brc{\bbeta} ,% \\
%\sigma_\Gamma^2 \brc{\bbeta}  &=  \frac{\alpha_{\mathrm{ d}} A_M }{M} \dbc{ 2 F_1\brc{\bbeta}  + \alpha_{\mathrm{ d}} A_M }   +  E\brc{\bbeta}\dbc{F_0\brc{\bbeta}  + F_1\brc{\bbeta}} ,
\end{align}

and (<ref>), given on the top of the next page, respectively.

\begin{align}
\sigma_\Gamma^2 \brc{\bbeta}  &=  \frac{\alpha_{\mathrm{ d}} A_M }{M} \dbc{ 2 F_1\brc{\bbeta}  + \alpha_{\mathrm{ d}} A_M }   +  E\brc{\bbeta}\dbc{F_0\brc{\bbeta}  + F_1\brc{\bbeta}} . \label{eq:sig_Gamma}
\end{align}

Let $\bbeta = [\beta_1,\ldots,\beta_N]^\trp$. For a given index $m\in [M]$, define $g_m\brc{\bbeta}$ be defined as

\begin{align}
g_m \brc{\bbeta} = \sum_{n=1}^{N} \mathrm{e}^{-\jj\beta_n} h_{\mathrm{r}, n} t_{nm}.
\end{align}

Here, $\beta_n$ and  $ t_{nm}$ are deterministic scalars, and $h_{\mathrm{r}, n}$ are jointly Gaussian for $n\in [N]$. We can hence conclude that $g_m\brc{\bbeta}$ is distributed Gaussian with mean and variance

\begin{align}
\mu_m \brc{\bbeta} & =\frac{ \mu_0 \brc{\bbeta}}{\sqrt{M}}
\mathrm{e}^{ - \jj \Phi_m\left( \varphi_{\mathrm{t}2} ,\theta_{\mathrm{t}2} \right) }, \label{eq:mu_m}\\
\tilde{\sigma}^2 \brc{\bbeta} &=  \bar{\alpha}_N\;
\mathbf{\bar{g}}_{\mathrm{r}}^\trp \mPhi \brc{\bbeta} \mR \mPhi^\her \brc{\bbeta}	\mathbf{\bar{g}}_{\mathrm{r}}^*,
\end{align}

respectively, where $\mathbf{\bar{g}}_{\mathrm{r}} = \mathbf{a}_N\left( \varphi_{\mathrm{r}1} , \theta_{\mathrm{r}1} \right)$, $\bar{\alpha}_N$ is defined in (<ref>),

\begin{align}
\mu_0 \brc{\bbeta} & =\sqrt{ M \kappa_\mathrm{r} \bar{\alpha}_N }
\mathbf{\bar{g}}_{\mathrm{r}}^\trp \mPhi \brc{\bbeta}	\mathbf{\bar{h}}_{\mathrm{r}},
\end{align}

$\mPhi \brc{\bbeta} = \mathrm{diag} \set{ \mathrm{e}^{-\jj \beta_1},\ldots,\mathrm{e}^{-\jj \beta_N} }$.

As $h_m$ for $m\in\dbc{M}$ have common term $g_m\brc{\bbeta}$, the entries of $\bh$ are correlated. We hence determine the covariance of $\bh$. To this end, we first note that

\begin{align}
\bmu\brc{\bbeta} = \dbc{\mu_1 \brc{\bbeta}, \ldots , \mu_M \brc{\bbeta}}^\trp = \mu_0 \brc{\bbeta}  \bv^*
\end{align}

where $\bv = \mathbf{a}_M\left( \varphi_{\mathrm{t}2} , \theta_{\mathrm{t}2} \right)/\sqrt{M}$. As the result, the covariance matrix is given by

\begin{align}
\mC &= \Ex{ \brc{\bh - \mu_0 \brc{\bbeta} \bv^*} \brc{\bh^\her - \mu_0^* \brc{\bbeta} \bv^\trp}},\\
&= \alpha_{\mathrm{ d}} A_M  \mI_M + M \tilde{\sigma}^2 \brc{\bbeta} \bv^* \bv^\trp = \mC_0^2,
\end{align}

where $\mC_0 = \sqrt{\alpha_{\mathrm{ d}} A_M } \; \mI_M + \varpi \brc{\bbeta} \bv^* \bv^\trp$ with $\varpi$ being

\begin{align}
\varpi \brc{\bbeta} =    \sqrt{\alpha_{\mathrm{ d}} A_M + M {\tilde{\sigma}^2 \brc{\bbeta} } } - \sqrt{\alpha_{\mathrm{ d}} A_M }.
\end{align}

We now represent the end-to-end channel $\bh$ as $\bh = \mC_0 \bh_0$, where $\bh_0 \sim \mathcal{CN}\left( \mC_0^{-1} \bmu\brc{\bbeta}, \mI_M\right)$. Consequently, we have

\begin{align}
\Gamma = \frac{1}{M} \bh_0^\her \mC \bh_0. \label{eq:Gam}%\brc{\mC_0 \bh_0 +\bmu\brc{\bbeta}}^\her \brc{\mC_0 \bh_0 +\bmu\brc{\bbeta}}\\
%\frac{\sigma^2\brc{\bbeta} }{M} V_{M}\brc{\bbeta}
%	\mathcal{I}_{\max}\left( \mathbf{h} \right) =\log_2 \left(1+ {\rho\sigma^2\brc{\bbeta} } V_{M}\brc{\bbeta}  \right),
\end{align}

To determine this weighted norm, we use the eigenvalue decomposition of $\mC$ to write $\mC = \mV^* \mSigma \mV^\trp$, where $\mV\in\setC^{M\times M}$ is a unitary matrix whose last column is $\bv$ and $\mSigma $ is an $M\times M$ diagonal matrix defined as

\begin{align}
\mSigma = \mathrm{diag} \set{ \alpha_{\mathrm{ d}} A_M , \ldots,\alpha_{\mathrm{ d}} A_M ,\alpha_{\mathrm{ d}} A_M +M  \tilde{\sigma}^2 \brc{\bbeta} }.
\end{align}

By replacing into (<ref>), we have

\begin{align}
\Gamma = \frac{\alpha_{\mathrm{ d}} A_M }{M}  \sum_{m=1}^{M-1} \abs{r_m}^2 + \frac{\alpha_{\mathrm{ d}}A_M +M \tilde{\sigma}^2 \brc{\bbeta} }{M}   \abs{r_M}^2 \label{eq:Gamma_1}%\brc{\mC_0 \bh_0 +\bmu\brc{\bbeta}}^\her \brc{\mC_0 \bh_0 +\bmu\brc{\bbeta}}\\
%\frac{\sigma^2\brc{\bbeta} }{M} V_{M}\brc{\bbeta}
%	\mathcal{I}_{\max}\left( \mathbf{h} \right) =\log_2 \left(1+ {\rho\sigma^2\brc{\bbeta} } V_{M}\brc{\bbeta}  \right),
\end{align}

where $r_m$ is the $m$-th entry of $\mathbf{r} = \mV^\trp \bh_0$. Noting that $\mV^\trp$ is unitary, we conclude that $\mathbf{r}\sim \mathcal{CN}\left(   \bmu_{\mathrm{r}}\brc{\bbeta}, \mI_M\right)$, where

\begin{align}
\bmu_{\mathrm{r}}\brc{\bbeta} = \mV^\trp\mC_0^{-1} \bmu \brc{\bbeta}.
\end{align}

Since $\mC = \mC_0^{2}$, the matrix $\mC_0^{-1}$ is decomposed as

\begin{align}
\mC_0^{-1} = \mV^* \sqrt{\mSigma}^{\; -1} \mV^\trp,
\end{align}

where $\sqrt{\mSigma}$ is a diagonal matrix whose first $M-1$ diagonal entries is $\sqrt{\alpha_{\mathrm{ d}} A_M }$ and the last one is $\sqrt{\alpha_{\mathrm{ d}} A_M }+\varpi \brc{\bbeta}$.

Hence, we can write

\begin{align}
\bmu_{\mathrm{r}}\brc{\bbeta} &= \mV^\trp \mV^* \sqrt{\mSigma}^{\; -1} \mV^\trp \bmu \brc{\bbeta}\\
&= \mu_0 \brc{\bbeta}  \sqrt{\mSigma}^{\; -1} \mV^\trp \bv^*.
\end{align}

Noting that $\bv$ is the last column of the unitary matrix $\mV$, we conclude that $\mV^\trp \bv^*$ is a standard base vector with the first $M-1$ being zero and the last entry being $1$. This concludes that $\mu_{\mathrm{r},m}\brc{\bbeta} = 0$ for $m\in \dbc{M-1}$ and

\begin{align}
\mu_{\mathrm{r},M}\brc{\bbeta} &= \frac{\mu_0 \brc{\bbeta} }{\sqrt{\alpha_{\mathrm{ d}}A_M }+\varpi \brc{\bbeta} }\\
&= \frac{\mu_0 \brc{\bbeta} }{\sqrt{\alpha_{\mathrm{ d}} A_M + M {\tilde{\sigma}^2 \brc{\bbeta} } }  } \label{eq:mu_r_M}
\end{align}

We now define independent random variables

\begin{align}
V_0 &=   \sum_{m=1}^{M-1} \abs{r_m}^2, \label{eq:v0}\\
V_1 &=    \abs{r_M}^2. \label{eq:v1}
\end{align}

By definition, $V_0$ is distributed chi-square with $2M-2$ degrees for freedom whose mean and variance are $M-1$ and whose density function is given by (<ref>).

The random variable $V_1$ is further distributed non-central chi-square with two degrees of freedom and non-centrality parameter

\begin{align}
\lambda\brc{\bbeta} = \abs{\mu_{\mathrm{r},M}\brc{\bbeta}}^2.
\end{align}

Consequently, the mean and variance of $V_1$ are $1+\lambda\brc{\bbeta}$ and $1+ 2 \lambda\brc{\bbeta}$, respectively, and its probability density is given by (<ref>) in Theorem <ref>.

The snr gain $\Gamma$ is determined in terms of $V_0$ and $V_1$ as

\begin{align}
\Gamma = \frac{\alpha_{\mathrm{ d}} A_M }{M}  V_0 + \frac{\alpha_{\mathrm{ d}} A_M +M \tilde{\sigma}^2 \brc{\bbeta} }{M}   V_1. \label{eq:Gamma_sum}
\end{align}

Noting that $V_0$ and $V_1$ are independent, the proof of Theorem <ref> is concluded after defining the functions $E\brc{\bbeta} = \tilde{\sigma}^2\brc{\bbeta}$ and $F_\zeta\brc{\bbeta} = \abs{\mu_0 \brc{\bbeta}}^2/M + \zeta E\brc{\bbeta}$.

§.§ Maximum Average End-to-End SNR

Using irs in the system, we are often interested in tuning the irs phase-shifts, such that the average end-to-end snr is maximized, i.e., $\mu_\Gamma \brc{\bbeta}$ is maximized over $\bbeta$. Considering Theorem <ref>, the maximum mean snr is not necessarily found in a closed-form. Nevertheless, closed-form lower and upper bounds can be derived in terms of the spectrum of $\mR$.

Let the phase-shifts at the  irs be set to the vector $ \bbeta^\star = [\beta^\star_1,\ldots,\beta^\star_N]^\trp$ with $\beta_n^\star$ being specified in (<ref>).

Then, the average end-to-end snr gain is bounded as

\begin{align}
\mu_\Gamma \brc{\bbeta^\star} &\leq  \xi_{\max} \brc{N} +  \kappa_{\mathrm{r}} \bar{\alpha}_N N^2 , \\
\mu_\Gamma \brc{\bbeta^\star} &\geq \xi_{\min} \brc{N} +  \kappa_{\mathrm{r}} \bar{\alpha}_N N^2
\end{align}

where $\xi_i \brc{N}$ for $i\in\set{\min , \max}$ is defined as

\begin{align}
\xi_i \brc{N} = {\alpha_{\mathrm{ d}} A_M+ \bar{\alpha}_N\lambda_i N } \label{eq:xi_i}
\end{align}

with $\bar{\alpha}_N$ being given in (<ref>), and $\lambda_{\min}$ and $\lambda_{\max}$ denoting the minimum and maximum eigenvalue of $\mR$, respectively. The variance of $\Gamma$ is further bounded as

\begin{align}
P_{\min}\brc{N} \leq \sigma_\Gamma^2 \brc{\bbeta^\star} &\leq P_{\max}\brc{N},
\end{align}

where $P_{i}\brc{N}$ is given for $i\in\set{\min,\max}$ in (<ref>) at the top of the next page.

\begin{align}
P_{i}\brc{N} =  \frac{ \alpha_{\mathrm{ d}} A_M }{M} +  \frac{2\alpha_{\mathrm{ d}} A_M \lambda_{i} }{M} \bar{\alpha}_N N
+ \brc{ \frac{2\alpha_{\mathrm{ d}} A_M \kappa_{\rm r} }{M} \bar{\alpha}_N +\bar{\alpha}_N^2 \lambda_{i}^2 } N^2
+ 2 \lambda_{i} \kappa_{\mathrm{ r} } \bar{\alpha}_N^2 N^3 \label{eq:upp_sig}
\end{align}

We start the proof by stating the following lemma:

For the function $F_0\brc{\bbeta}$ defined in (<ref>), we have

\begin{align}
\max_{\beta_1,\ldots,\beta_N} F_0 \brc{\bbeta}  = \kappa_{\mathrm{r}} \bar{\alpha}_N N^2,
\end{align}

being obtained by setting $\beta_n = \beta^\star_n$ given in (<ref>) for $n\in \dbc{N}$.

The proof is given in the Appendix <ref>.

By setting $\beta_n = \beta^\star_n$ in Theorem <ref>, the mean and variance of the snr gain reduce to

\begin{align}
\mu_\Gamma \brc{\bbeta^\star} &=  \alpha_{\mathrm{ d}} A_M+  E\brc{\bbeta^\star} + \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2 ,
%\sigma_\Gamma^2 \brc{\bbeta^\star}  &=  \frac{\alpha_{\mathrm{ d}}}{M} \dbc{ 2\kappa_{\mathrm{r}}  \bar{\alpha}_N N^2+ 2 E\brc{\bbeta^\star} + \alpha_{\mathrm{ d}} }  +  %E\brc{\bbeta^\star}\dbc{ E\brc{\bbeta^\star} + 2 \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2} ,
%\sigma_\Gamma^2 \brc{\bbeta^\star} &= \frac{E\brc{\bbeta^\star}}{M} \brc{  E\brc{\bbeta^\star}+ 2  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2   }.
\end{align}

and the expression given in (<ref>) at the top of the next page, respectively.

\begin{align}
%	\mu_\Gamma \brc{\bbeta^\star} &=    E\brc{\bbeta^\star} + \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2 , \\
\sigma_\Gamma^2 \brc{\bbeta^\star}  &=  \frac{\alpha_{\mathrm{ d}} A_M }{M} \dbc{ 2\kappa_{\mathrm{r}}  \bar{\alpha}_N N^2+ 2 E\brc{\bbeta^\star} + \alpha_{\mathrm{ d}} A_M }  +  E\brc{\bbeta^\star}\dbc{ E\brc{\bbeta^\star} + 2 \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2} , \label{eq:sig_star}
\end{align}

\begin{align}
\sigma_{\Gamma,i}^2   &=  \frac{\alpha_{\mathrm{ d}} A_M }{M} \brc{ 2\kappa_{\mathrm{r}}  \bar{\alpha}_N N^2+ 2 \maE_{i} - \alpha_{\mathrm{ d}} A_M }  +  \brc{\maE_{i} - \alpha_{\mathrm{ d}} A_M } \brc{ \maE_{i} + 2 \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2 - \alpha_{\mathrm{ d}} A_M } , \label{eq:sig_max}
%\sigma_\Gamma^2 \brc{\bbeta^\star} &= \frac{E\brc{\bbeta^\star}}{M} \brc{  E\brc{\bbeta^\star}+ 2  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2   }.
\end{align}

The mean can hence be bounded as

\begin{align}
\maE_{\min} + \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2 \leq	\mu_\Gamma \brc{\bbeta^\star} \leq \maE_{\max} + \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2, \label{eq:Mu_B}
\end{align}

and the variance reads

\begin{align}
\sigma_{\Gamma,\min}^2 \leq \sigma_\Gamma^2 \brc{\bbeta^\star} &\leq  \sigma_{\Gamma,\max}^2, \label{eq:SigB} %\frac{E_{\max} }{M} \brc{  E_{\max}+ 2  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2   }. \label{eq:SigB}
\end{align}

with $\sigma_{\Gamma,i}^2$ being defined in (<ref>) at the top of the next page for $i\in \set{\min,\max}$. Here,

\begin{align}
\maE_{\max} &= \alpha_{\mathrm{ d}} A_M+  \max_{\bbeta\in\setR^N} E\brc{\bbeta} \\
\maE_{\min} &= \alpha_{\mathrm{ d}} A_M+  \min_{\bbeta\in\setR^N} E\brc{\bbeta}
\end{align}

whose values are given by the following lemma:

For the function $E\brc{\bbeta}$ defined in (<ref>), we have

\begin{align}
\xi_{\min} \brc{N} \leq \alpha_{\mathrm{ d}} A_M + E \brc{\bbeta}  \leq \xi_{\max} \brc{N}
\end{align}

for any $\bbeta\in\setR^N$, where the function $\xi_i\brc{N}$ is defined in (<ref>) for $i\in \set{\min,\max}$.

The proof is given in the Appendix <ref>.

By substituting the results of Lemma <ref> into the bounds in (<ref>) and (<ref>), the proof is concluded.

From the proof, it is straightforward to conclude that the upper-bounds given in Lemma <ref> are in general valid for any choice of irs phase-shifts. This means that, for any $\bbeta\in\setR^N$, we have

\begin{align}
\mu_\Gamma \brc{\bbeta} &\leq  \xi_{\max} \brc{N} +  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2 , \\
\sigma_\Gamma^2 \brc{\bbeta} &\leq P_{\max}\brc{N}.
\end{align}

Lemma <ref> explains the logic behind the choice of irs phase-shifts in Proposition <ref>. It is worth noting that the proposed phase-shifts in (<ref>) do not necessarily maximize the average end-to-end snr. Nevertheless, they are still good enough to guarantee the end-to-end channel hardening.

§.§ Proof of Proposition <ref>

The proof of Proposition <ref> follows readily from Theorems <ref> and <ref>: Using Theorem <ref>, we represent the end-to-end snr gain as $\Gamma = \mu_\Gamma \brc{\bbeta^\star}+\sigma_\Gamma \brc{\bbeta^\star} \tilde{\Gamma}$, where $\tilde{\Gamma}$ is a zero-mean and unit-variance generalized chi-square random variable denoting the centralized and normalized version of $\Gamma$. We now replace it in (<ref>) and use the Taylor series of the logarithm to write

\begin{align}
\mathcal{C} &=\log_2 \left(1+  \rho M \mu_\Gamma\brc{\bbeta^\star}+ \rho M \sigma_\Gamma\brc{\bbeta^\star} \tilde{\Gamma} \right),\\
&= \log_2 \left(1 +  \rho M \mu_\Gamma\brc{\bbeta^\star} \right)  + \log_2 \left(1+  \varsigma_N \tilde{\Gamma} \right)\\
&= \log_2 \left(1 +\rho M \mu_\Gamma\brc{\bbeta^\star}  \right)  +  \varsigma_N \log_2\mathrm{e}  \; \tilde{\Gamma} + \epsilon_N, \label{eq:taylor}
\end{align}

where $\varsigma_N$ is defined as[We use subscript $N$ to indicate its dependency on dimension $N$.]

\begin{align}
\varsigma_N = \frac{ \rho M \sigma_\Gamma \brc{\bbeta^\star}   }{1 + \rho M \mu_\Gamma \brc{\bbeta^\star}  },
\end{align}

and $ \epsilon_N$ is a polynomial in $\varsigma_N\tilde{\Gamma}$ satisfying

\begin{align}
\epsilon_N = \mathcal{O}\left(  {\varsigma_N^2 \tilde{\Gamma}^2} \right).
\end{align}

Invoking Theorem <ref>, we can bound $\varsigma_N$ from above as

\begin{align}
\varsigma_N &\leq \frac{ \rho M \sqrt{ P_{\max} \brc{N} }  }{1 + \rho M \brc{ \xi_{\min} \brc{N} +  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2} }\\
&\stackrel{\dagger}{\leq} \frac{ \rho M \sqrt{ P_{\max} \brc{N} }  }{1 + \rho M \brc{ \alpha_{\mathrm{ d}} A_M  +  \kappa_{\mathrm{r}}  \bar{\alpha}_N N^2} }
\end{align}

where $\dagger$ follows from $\xi_{\min} \brc{N}\geq \alpha_{\mathrm{ d}} A_M $.

We now focus on scaling with $N$. Since $\lambda_{\max}$ grows sub-linearly in $N$, we can conclude that there exist real scalars $\lambda_0$ and  $0\leq u <1$, such that $\lambda_{\max}$ is uniformly bounded by

\begin{align}
\lambda^{\rm U} \brc{N} &= \lambda_0 N^{u},
\end{align}

From the constraint $A_{\rm IRS} = A_0 N^{1-q}$, we can further conclude that

\begin{align}
A_N &= \frac{A_0}{N^q},
\end{align}

for some real $A_0$ and $0\leq q <1$, where the constraint in (<ref>) guarantees that $u\geq q$. As a result, $P_{\max} \brc{N}$ is uniformly bounded from above by

\begin{align}
P^{\rm U}\brc{N} &= b N^{3+u - 4q}.
\end{align}

These upper-bounds lead to this conclusion that $\varsigma_N$ is bounded as $\varsigma_N \leq \varsigma_N^{\rm U}$ for some $\varsigma_N^{\rm U}$ satisfying

\begin{align}
\varsigma_N^{\rm U} = \mathcal{O}\brc{N^{\tfrac{u-1}{2}}}.
\end{align}

We now note that $\tilde{\Gamma}$ is zero-mean and unit-variance, and $0\leq u < 1$. Thus, $\epsilon_N \rightarrow 0$ as $N$ grows large, and hence $\mathcal{C}$ is well approximated by the first two terms in (<ref>).

We now use the following lemma:

Let $F_0\brc{\bbeta}/\Lambda\brc{\bbeta}$ and $\sigma^2_\Gamma\brc{\bbeta}$ grow large. Then, we have

\begin{align}\label{eq:V}
\dfrac{\Gamma-\mu_\Gamma \brc{\bbeta}}{\sigma_\Gamma\brc{\bbeta} }\xrightarrow{\mathrm{ d}} \mathcal{N}\left( 0, 1\right).
\end{align}

The proof is given in the Appendix <ref>.

We first note that $\sigma^2_\Gamma\brc{\bbeta^\star}$ grows large with $N$. Moreover, $F_0\brc{\bbeta^\star} = \bar{\alpha}_N\kappa_{\mathrm{ r} } N^2$ and $\Lambda\brc{\bbeta^\star}$ is uniformly bounded from above with $ \mathcal{O}\brc{N^{1+u - 2q}}$ for some $0\leq u <1$. This concludes that $F_0\brc{\bbeta^\star}/\Lambda\brc{\bbeta^\star}$ grows large as $N$ increases. We hence use Lemma <ref> and conclude that $\tilde{\Gamma}$ converges in distribution to a zero-mean and unit-variance Gaussian random variable, as $N$ grows asymptotically large.

Finally by noting that $\bar{\mathbf{g}}_{\mathrm{ r}}^\trp \mPhi\brc{\bbeta^\star} = \mathbf{\bar{h}}_{\mathrm{r}}^\her$, we can write

\begin{align}
E\brc{\bbeta^\star}  =  \bar{\alpha}_N\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}}.
\end{align}

Substituting into (<ref>), Proposition <ref> is concluded after few lines of standard calculations.

§.§ Derivation of Pseudo-Proposition <ref>

From our numerical validations in Section <ref> we heuristically conclude that $\mathcal{C}$ is still well-approximated by a Gaussian random variable in the linearly scaling settings. Nevertheless, the variance of $\mathcal{C}$ in this case is not tightly approximated by Proposition <ref>. From the proof of Proposition <ref>, the inaccuracy comes from the fact that the necessary condition for Lemma <ref> is not satisfied and hence, $\tilde{\Gamma}$ does not converge in distribution to a zero-mean and unit-variance Gaussian random variable[This is why Proposition <ref> excludes the strictly-linear case.].

The exact distribution of $\tilde{\Gamma}$ in this case can be derived explicitly from Theorem <ref>. We however follow our heuristic and approximate it by a zero-mean Gaussian random variable of a smaller variance. Following the proof of Lemma <ref>, we show in Appendix <ref> that $\tilde{\Gamma}$ can be approximated in this case as

\begin{align}
\tilde{\Gamma} \approx \tilde{\Gamma}_\infty + \hat{\epsilon} \label{eq:approx_Gam}
\end{align}

where $\tilde{\Gamma}_\infty\sim \mathcal{N}\brc{0, {\sigma}_{\infty}^2 }$ with

\begin{align}
{\sigma}_{\infty}^2 = \frac{\kappa_{\mathrm{ r} }}{\kappa_{\mathrm{ r} } + \vartheta } \label{eq:sig_1_final}
\end{align}

and $\vartheta$ being given in (<ref>). The random variable $\hat{\epsilon}$ is further a centralized[That means $\hat{\epsilon}$ is zero-mean.] chi-square residual term. By replacing (<ref>) in the Taylor expansion (<ref>) and ignoring the residual term as well as higher-order terms, the approximation is derived.

§ CONCLUSIONS

Even with highly correlated reflecting elements, a large irs hardens the end-to-end channel of a mimo system with small transmit and receive array antennas. This conclusion follows from the main result of this study which shows that the mutual information between the input and output signals in an irs-aided system is almost perfectly approximated by a Gaussian random variable, whose mean increases with the irs size and whose variance vanishes to zero as the irs size grows large.

The results of this study indicate that the channel hardening property is rather generic in irs-aided systems and is easily achieved by small transmit and receive array antennas. More interestingly, the phase-shifts of reflecting elements are not required to be updated frequently and should only be appropriately matched with the large-scale statistics of the wireless channels, i.e., aoa andaod. These findings reveal this interesting fact: Even with irs of moderate size whose physical dimensions are rather significantly smaller than the wavelength, we can still harden the end-to-end mimo channel between the transmitter and the receiver.

Earlier investigations in the literature have demonstrated the promising performance of massive mimo systems in various respects. Nevertheless, the complexity of implementing such systems left this technology challenging for commercial uses.  Along with earlier results in the literature, the result of this study indicate that using irs, the key features of the massive mimo technology can be achieved with rather small end-to-end dimensions.

Although this work studies a basic setting, the analytical results can still be employed to investigate various fundamental properties of irs-aided mimo systems. An example is given in Section <ref>, where we demonstrate the dimensional trade-off between the irs and the bs: The larger the irs is, the smaller the bs needs to be in order to achieve a target performance. Similar investigations and further extensions of the results to more advanced settings, and under more realistic assumptions, are interesting study directions. They are however out of the scope of this particular paper and are left as potential directions for future work.

§ PROOF OF LEMMA <REF>

Considering the definition of $F_0\brc{\bbeta}$, one can write

\begin{align}
F_0\brc{\bbeta} &=\bar{\alpha}_N\kappa_{\mathrm{ r} }
\abs{
\sum_{n=1}^{N} 	\mathrm{e}^{
\jj \brc{\Pi_n\left( \varphi_{\mathrm{r}1} ,\theta_{\mathrm{r}1} \right)  + \Pi_n\left( \varphi_{\mathrm{t}1} ,\theta_{\mathrm{t}1} \right) - \beta_n}
\end{align}

The expression on the right hand side includes the amplitude of a sum whose summands are complex numbers on the unit circle. As the result, the amplitude of the sum is maximized by setting all the summands in-phase. This is obtained by setting $\beta_n = \beta_n^\star$ for $\beta_n^\star$ given in (<ref>). In this case, the sum adds to $N$ and $F\brc{\bbeta^\star}$ is given by the expression given in Lemma <ref>.

§ PROOF OF LEMMA <REF>

Let $\btau\brc{\bbeta}  = \mPhi^\her\brc{\bbeta} \mathbf{g}_{\mathrm{r}}^*$. Hence, $E\brc{\bbeta}$ can be written as

\begin{align}
E\brc{\bbeta} &= \bar{\alpha}_N\;
\btau^\her\brc{\bbeta}  {\mR} \btau \brc{\bbeta} \label{eq:e_neu}
\end{align}

Since every entry of $\btau\brc{\bbeta}$ lies on the unit circle, we can write $\norm{\btau\brc{\bbeta}}^2 = N$ for any $\bbeta\in\setR^N$. Consequently, the function

\begin{align}
Q\brc{\bbeta} = \frac{1}{N} \btau^\her\brc{\bbeta} \mR \btau\brc{\bbeta}
\end{align}

determines the Rayleigh quotient of $\mR$ at $\btau\brc{\bbeta}$. We now use the fact that the Rayleigh quotient is bounded in terms of the eigenvalue of $\mR$, as follows <cit.>

\begin{align}
\lambda_{\min} \leq Q\brc{\bbeta} \leq \lambda_{\max},
\end{align}

where $\lambda_{\min}$ and $\lambda_{\max}$ are the minimum and maximum eigenvalue of $\mR$, respectively. Substituting into (<ref>), Lemma <ref> is concluded.

§ PROOF OF LEMMA <REF>

We start the proof by considering (<ref>). Using the definition of  $\Lambda\brc{\bbeta}$ in Theorem <ref>, we have

\begin{align}
\Gamma = \frac{\alpha_{\mathrm{ d}} A_M }{M}  \sum_{m=1}^{M-1} \abs{r_m}^2 + \Lambda\brc{\bbeta}  \abs{r_M}^2.
\end{align}

Here, $r_m$ for $m\in\dbc{M-1}$ are independent complex Gaussian random variables with zero mean and unit variance, and $r_M$ is a complex unit-variance Gaussian random variable independent of $r_m$ for $m\in\dbc{M-1}$. The mean of $r_M$, i.e., $\mu_{\mathrm{r},M}\brc{\bbeta}$, is given in (<ref>). Using the definitions given in Theorem <ref>, it can be shown that

\begin{align}
\abs{\mu_{\mathrm{r},M}\brc{\bbeta}} = \sqrt{\frac{F_0 \brc{\bbeta} }{\Lambda\brc{\bbeta}} }.
\end{align}

Let $r_M = \tilde{r}_M + \mu_{\mathrm{r},M}\brc{\bbeta}$ where $ \tilde{r}_M$ is the centralized form of $r_M$. We can hence write

\begin{align}
\Gamma = T_0  + 2 \Lambda\brc{\bbeta}  \Re\set{\tilde{r}_M \mu_{\mathrm{r},M}^*\brc{\bbeta} } + \abs{\mu_{\mathrm{r},M}\brc{\bbeta}}^2 \label{eq:Gamma_T}
\end{align}

where $T_0$ is defined as

\begin{align}
T_0 = \frac{\alpha_{\mathrm{ d}} A_M }{M}  \sum_{m=1}^{M-1} \abs{r_m}^2 + \Lambda\brc{\bbeta}  \abs{\tilde{r}_M}^2.
\end{align}

We now consider the normalized snr gain, i.e.,

\begin{align}
\hat{\Gamma} &= \frac{\Gamma}{\sigma_\Gamma\brc{\bbeta}}
\end{align}

By replacing $\Gamma$ with the expression in (<ref>), we have

\begin{align}
\hat{\Gamma} &= \Gamma_0 + \Gamma_1 + c \label{eq:Decomp}
\end{align}

where $\Gamma_0$ and $\Gamma_1$ are random expressions defined as

\begin{align}
\Gamma_0 &= \frac{T_0}{\sigma_\Gamma\brc{\bbeta}}  \\
\Gamma_1 &= 2 \frac{\Lambda\brc{\bbeta}}{\sigma_\Gamma\brc{\bbeta}}  \Re\set{\tilde{r}_M \mu_{\mathrm{r},M}^*\brc{\bbeta} }
\end{align}

and $c$ is a deterministic constant given by

\begin{align}
c = \frac{\abs{\mu_{\mathrm{r},M}\brc{\bbeta}}^2}{\sigma_\Gamma\brc{\bbeta}}.
\end{align}

Despite its complicated form of $\sigma^2_\Gamma\brc{\bbeta}$, we can use the fact that $F_0\brc{\bbeta} \leq F_1 \brc{\bbeta}$ in Theorem <ref>, and bound $\sigma^2_\Gamma\brc{\bbeta}$ as

\begin{align}
2 \Lambda\brc{\bbeta} F_0\brc{\bbeta} \leq \sigma^2_\Gamma\brc{\bbeta} \leq 2 \Lambda\brc{\bbeta} F_1\brc{\bbeta} + \alpha_{\mathrm{ d}}^2 A_M^2 .
\end{align}

Noting that $\alpha_{\mathrm{ d}}^2 A_M^2$ is fixed, we can conclude that growth of $\sigma^2_\Gamma\brc{\bbeta}$ guarantees that $ \Lambda\brc{\bbeta} F_0\brc{\bbeta}$ grows large.

We use the upper bound to write

\begin{align}
\Gamma_0   &\leq  \frac{T_0}{2\sqrt{\Lambda \brc{\bbeta} F_0\brc{\bbeta} }}\\
&= \frac{\displaystyle \alpha_{\mathrm{ d}} A_M \sum_{m=1}^{M-1} \abs{r_m}^2 }{2M \sqrt{\Lambda \brc{\bbeta} F_0\brc{\bbeta}}}   + \sqrt{ \frac{ \Lambda\brc{\bbeta} }{ F_0\brc{\bbeta} } } \abs{\tilde{r}_M}^2 .
\end{align}

Now, let $\sigma^2_\Gamma\brc{\bbeta}$  and $ F_0\brc{\bbeta} /  \Lambda\brc{\bbeta}$ grow large. We can then conclude that the upper bound in this case converges to zero in the mse. Since $\Gamma_0 \geq 0$, we can conclude that $\Gamma_0$ converges to zero.

The expression $\Gamma_1$ is further a real-valued zero-mean Gaussian random variable whose variance is given by

\begin{align}
\sigma_1^2 &=  \frac{2\Lambda^2\brc{\bbeta}}{\sigma_\Gamma^2\brc{\bbeta}}  \abs { \mu_{\mathrm{r},M}\brc{\bbeta} }^2. \label{sig_1}
\end{align}

Using the lower and upper bounds on the variance of $\Gamma$, we can bound $\sigma_1^2$ as

\begin{align}
\frac{F_0\brc{\bbeta}}{F_1\brc{\bbeta}} \leq \sigma_1^2 \leq 1.
\end{align}

Noting that

\begin{align}
F_1\brc{\bbeta} = F_0\brc{\bbeta} + \Lambda \brc{\bbeta} - \frac{\alpha_{\mathrm{ d}} A_M }{M}
\end{align}

we can further write

\begin{align}
\lim_{\tfrac{F_0\brc{\bbeta}}{\Lambda\brc{\bbeta}} \rightarrow \infty }	\frac{F_0\brc{\bbeta}}{F_1\brc{\bbeta}} = 1
\end{align}

which means that $\sigma_1^2 = 1$.  As a result, by growth of $\sigma^2_\Gamma\brc{\bbeta}$ and $F_0\brc{\bbeta} /  \Lambda\brc{\bbeta}$, the normalized snr gain $\hat{\Gamma}$ converges in distribution to a unit-variance real Gaussian random variable.

Finally, by noting that $\brc{\Gamma - \mu_\Gamma} / \sigma_\Gamma$ is the centralized form of $\hat{\Gamma}$, we can conclude that

\begin{align}
\dfrac{\Gamma-\mu_\Gamma \brc{\bbeta}}{\sigma_\Gamma\brc{\bbeta} }\xrightarrow{\mathrm{ d}} \mathcal{N}\left( 0, 1\right).
\end{align}

This concludes the proof.

§ APPROXIMATED SNR FOR LINEAR SCALING

For a linearly scaling irs area and $\lambda_{\max}$, the conditions in Lemma <ref> are not satisfied, and hence the given limit is not valid. Considering the proof in Appendix <ref>, this follows from the fact that due to linear scaling $\Gamma_0$ does not converge to zero, and accordingly the limit of $\sigma_1^2$ is not one. Considering (<ref>) in Appendix <ref>, we can write

\begin{align}
\tilde{\Gamma} &= \hat{\Gamma} - \frac{\mu_\Gamma\brc{\bbeta} }{\sigma_\Gamma\brc{\bbeta}} = \Gamma_1 + \hat{\epsilon},
\end{align}

where $\hat{\epsilon}$ is a centralized chi-square random variable[$\hat{\epsilon}$ is zero-mean, since $\tilde{\Gamma}$ is zero-mean]. The random variable $\Gamma_1$ is further zero-mean Gaussian whose variance is $\sigma_1^2$ given in (<ref>).

Considering linear scaling, i.e., $q=u = 1$, $\mR$ is a rank-one matrix and hence we can write $\mR = \mathbf{r} \mathbf{r}^\her$ for some $\mathbf{r}\in\setC^N$. As a result,

\begin{align}
\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} = \abs{\mathbf{h}^\her \mathbf{r} }^2.
\end{align}

Noting that entries of $\mathbf{h}^\her$ lie equidistantly on the unit-circle, one can approximately write $\mathbf{\bar{h}}_{\mathrm{r}}^\her \mR \mathbf{\bar{h}}_{\mathrm{r}} \approx 0$ for a typical[For instance, for the covariance matrix $\mR = \boldsymbol{1}_N$, $\mathbf{r}$ is a vector of all-ones, and hence this approximation is almost exact.] $\mathbf{r}$.

We further note that $\bar{\alpha}_N N^2$ converges to a constant with $q=1$. Considering this fact and using the above approximation, the variance $\sigma_1^2$ is approximated with $\sigma_{\infty}^2$  given in (<ref>).

mimo[MIMO]multiple-input multiple-output
miso[MISO]multiple-input single-output
siso[SISO]single-input single-output
csi[CSI]channel state information
awgn[AWGN]additive white Gaussian noise
iid[i.i.d.]independent and identically distributed
ut[UT]user terminal
bs[BS]base station
snr[SNR]signal-to-noise ratio
rf[RF]radio frequency
los[LoS]line of sight
nlos[NLoS]non-line of sight
irs[IRS]intelligent reflecting surface
ula[ULA]uniform linear array
mrt[MRT]maximum ratio transmission
noma[NOMA]non-orthogonal multiple access
mse[MSE]mean squared error