# Two-Timescale Design for Reconfigurable Intelligent Surface-Aided Massive
MIMO Systems with Imperfect CSI

Kangda Zhi, Cunhua Pan, Hong Ren, Kezhi Wang, Maged Elkashlan, Marco Di Renzo,
, Robert Schober, , H. Vincent Poor, , Jiangzhou Wang, and Lajos Hanzo
(Corresponding author: Cunhua Pan). Part of this work will be presented in the
IEEE SPAWC, 2021 [1]. K. Zhi, C. Pan, M Elkashlan are with the School of
Electronic Engineering and Computer Science at Queen Mary University of
London, London E1 4NS, U.K. (e-mail: k.zhi, c.pan,
maged.elkashlan@qmul.ac.uk). H. Ren is with the National Mobile Communications
Research Laboratory, Southeast University, Nanjing 210096, China. (e-mail:
hren@seu.edu.cn). K. Wang is with Department of Computer and Information
Sciences, Northumbria University, UK. (e-mail: kezhi.wang@northumbria.ac.uk).
M. Di Renzo is with Université Paris-Saclay, CNRS and CentraleSupélec,
Laboratoire des Signaux et Systèmes, Gif-sur-Yvette, France. (e-mail:
marco.direnzo@centralesupelec.fr). R. Schober is with the Institute for
Digital Communications, Friedrich-Alexander-University Erlangen-Nürnberg
(FAU), Germany (e-mail: robert.schober@fau.de). H. V. Poor is with the
Department of Electrical and Computer Engineering, Princeton University,
Princeton, NJ 08544 USA (e-mail: poor@princeton.edu). Jiangzhou Wang is with
the School of Engineering and Digital Arts, University of Kent, UK. (e-mail:
J.Z.Wang@kent.ac.uk). Lajos Hanzo is with the School of Electronics and
Computer Science, University of Southampton, Southampton, SO17 1BJ, U.K.
(e-mail: lh@ecs.soton.ac.uk.).

###### Abstract

This paper investigates the two-timescale transmission design for
reconfigurable intelligent surface (RIS)-aided massive multiple-input
multiple-output (MIMO) systems, where the beamforming at the base station (BS)
is adapted to the rapidly-changing instantaneous channel state information
(CSI), while the passive beamforming at the RIS is adapted to the slowly-
changing statistical CSI. Specifically, we first consider the ideal model with
spatial-independent Rician fading channels, which leads to tractable
expressions and admits analytical insights including the power scaling laws
and the impact of various system parameters. Then, we extend to the practical
model with spatial-correlated Rician fading channels and consider the
electromagnetic interference (EMI) received by the RIS, which is inevitably
produced by any uncontrollable sources in the environment. Under both models,
we first perform the linear minimum mean square error (LMMSE) estimator to
estimate the aggregated channel from the users to the BS, then apply the low-
complexity maximal ratio combining (MRC) detector, and finally derive the
lower bound of the achievable rate in a closed form. Besides, a novel gradient
ascent-based algorithm is proposed for solving the minimal user rate
maximization problem. Numerical results validate that the spatial-independent
model without EMI can achieve a high accuracy approximation for the practical
model when the element spacing of RIS is large and when the EMI is light. We
also find that with spatial correlation, the RIS could possess a greater
capacity to tailor the wireless environment. Furthermore, it is shown that
integrating RIS into massive MIMO systems brings significant gains when the
RIS is deployed close to some cell-edge users, while the gains of using RIS
for users distributed over a large area are relatively small.

###### Index Terms:

Reconfigurable intelligent surface (RIS), massive MIMO, two-timescale design,
channel estimation, spatial correlation, electromagnetic interference (EMI).

##  I Introduction

As an emerging candidate for next-generation communication systems,
reconfigurable intelligent surfaces (RISs), also termed intelligent reflecting
surfaces (IRSs), have attracted significant interest from both academia and
industry[2, 3]. An RIS is a reconfigurable engineered surface that does not
require active radio frequency (RF) chains, power amplifiers, and digital
signal processing units, and is usually made of a large number of low cost and
passive scattering elements that are coupled with simple low power electronic
circuits. By intelligently tuning the phase shifts of the impinging wave with
the aid of a controller, an RIS can constructively strengthen the desired
signal power or deconstructively weaken the interference signal, which yields
an appealing passive beamforming gain.

Compared with existing multi-antenna systems[4, 5, 6, 7, 8, 9], it has been
demonstrated that RIS-aided systems have the potential to achieve better
performance in terms of cost and energy consumption[10, 11, 12, 13, 14, 15,
16, 17]. Recently, RISs have been considered for integration into various
communication scenarios, such as terahertz, sub-terahertz, and millimeter-wave
systems[18, 19], simultaneous wireless information and power transfer
(SWIPT)[20], unmanned aerial vehicle (UAV) communications[21], cell-free
systems[22], physical-layer security[23, 24, 25], mobile edge computing
(MEC)[26, 27, 28], device-to-device (D2D) communications[29, 30] and so on.
Furthermore, the effectiveness of RIS-aided systems in the presence of
practical imperfections has been demonstrated in [31, 32, 33, 34].
Specifically, relying on imperfect instantaneous channel state information
(CSI), the robust transmission design of RISs was studied in [31, 32]. The
authors in [33] studied the RIS beamforming design by considering transceiver
hardware impairments. With the consideration of RF impairments and phase
noises, the authors in [34] successfully conducted theoretical analysis on the
fundamental tradeoffs between spectral and energy efficiency of the RIS
communication network. In addition, a valuable experimental investigation of
the RIS-assisted channel was carried out in [35].

While several benefits of RISs have been demonstrated in the above-mentioned
contributions, most of them considered the design of the passive beamforming
at the RIS based on the knowledge of instantaneous CSI in each channel
coherence interval. In practice, these instantaneous CSI-based schemes face
two challenges. The first one is the overhead for the acquisition of the
instantaneous CSI. Due to the passive nature of the reflecting elements that
constitute the RIS, many authors proposed to estimate the cascaded user-RIS-BS
channels instead of the separated user-RIS and RIS-BS channels[36, 37]. The
pilot overhead of these channel estimation schemes is still proportional to
the number of RIS elements. However, an RIS generally consists of a large
number of reflecting elements to reap its potential[38], which incurs a
prohibitively high pilot overhead. Secondly, in each channel coherence time
interval, the BS needs to calculate the optimal beamforming coefficients for
the RIS, and send them back to the RIS controller via dedicated feedback
links. For instantaneous CSI-based scheme, therefore, the beamforming
calculation and information feedback need to be executed frequently in each
channel coherence interval, which results in a high computational complexity,
feedback overhead, and energy consumption.

To address these two practical challenges, recently, Han _et al._[39] firstly
proposed a novel two-timescale based RIS design scheme, which facilitates the
deployment and operation of RIS-aided systems. This promising two-timescale
scheme was further analyzed in recent research works [40, 41, 42, 43, 44, 45,
46, 47, 48, 49, 50, 51]. In the two-timescale scheme, the BS beamforming is
designed based on the instantaneous aggregated CSI, which includes the direct
and RIS-reflected links. The dimension of this aggregated channel is the same
as for conventional RIS-free systems, which is not related to the number of
RIS reflecting elements. Hence, in the two-timescale scheme, the length of the
pilot sequences only needs to be larger than the number of users, which
significantly reduces the channel estimation overhead. More importantly, the
two-timescale scheme aims to design the RISs only based on long-term
statistical CSI, such as the locations and the angles of the users with
respect to the BS and the RIS, which vary much slower than the instantaneous
CSI in the sub-6 GHz bands. The phase shifts of the reflecting elements of the
RIS need to be redesigned only when the large-scale channel information
changes. Compared with instantaneous CSI-based designs that need to update the
phase shifts of the reflecting elements in each channel coherence interval,
therefore, RIS-aided designs based on statistical CSI can significantly reduce
the computational complexity, feedback overhead and energy consumption.

In addition, massive MIMO technology has been identified as the cornerstone of
the fifth generation (5G) and future communication systems[52, 53]. Massive
MIMO exploits tens or hundreds of BS antennas to serve multiple users
simultaneously. Due to the complexity of wireless propagation environments,
e.g., the presence of large blocking objects, however, the signal power
received at the end-users may be still too weak, and it may be insufficient to
support emerging applications that entail high date rate requirements, such as
virtual reality (VR) or augment reality (AR). Inspired by the capability of
RISs to customize the wireless propagation environment, a natural idea is to
integrate them into massive MIMO systems. By constructing alternative
transmission paths, it is envisioned that RIS-aided massive MIMO systems can
achieve significant performance gains, especially when the direct links
between the BS and the users are blocked by obstacles. In RIS-aided massive
MIMO systems, the transmission scheme needs to be carefully designed, and the
channel estimation overhead needs to be taken into account considering the
large channel dimension. The application of instantaneous CSI-assisted
schemes, in particular, may lead to a prohibitive complexity and overhead.
Instead, due to the reduced channel estimation and feedback overhead, the two-
timescale scheme is deemed more suitable for RIS-aided massive MIMO systems.

Even though RIS-aided massive MIMO systems have been investigated in some
recent works [54, 55, 49, 50], there still exist three key open issues.
Firstly, it is crucial to provide comprehensive performance analysis for RIS-
aided massive MIMO systems based on the two-timescale design and the imperfect
CSI. In the presence of channel estimation error, the impact of key system
parameters, the achievable rate scaling order, and the possible power scaling
laws are completely unknown. To answer these questions, it is necessary to
derive explicit information-theory expressions and obtain analytical insights.
Secondly, since spatial diversity is a key factor for massive MIMO systems and
considering the heights of the BS and the RIS, it is essential to adopt the
Rician channel model including line-of-sight (LoS) and non-LoS (NLoS)
components so that the impact of the degree of scatterers can be modeled and
analyzed. This also enables us to provide guidelines with respect to the
deployment of the RIS. Thirdly, there are some unique and realistic
characteristics that must be considered in RIS-aided systems, including the
spatial correlation at the RIS and the electromagnetic interference (EMI)
captured by the RIS. These two factors might bring significant difference to
conclusions while it has not been examined in RIS-aided massive MIMO systems
with two-timescale design and imperfect CSI. To be specific, due to the planar
structure of the RIS, spatial correlation inevitably exists among the diagonal
elements of the RIS [56], and therefore it is more practical to adopt the
correlated Rician channel model. Meanwhile, since the RIS is commonly
comprised of a large number of elements, it could easily capture the EMI from
its surrounding wireless environment, which is generated by any uncontrollable
external sources (e.g., signals from adjacent cells and natural background
radiation) [57, 58]. Therefore, the EMI re-radiated by a large RIS towards the
BS might deteriorate the channel estimation quality and reduce the end-to-end
SINR, especially when the RIS is large and when the useful signal power is
weak. The above-mentioned three questions motivate our work.

In this paper, we focus on the uplink (UL) two-timescale transmission of an
RIS-aided massive MIMO system that is subject to imperfect aggregated CSI. The
Rician channel model is adopted to model the impact of spatial diversity. To
gain some basic insights, we first utilize the ideal model with spatial-
independent Rician fading which admits tractable achievable rate expressions,
and it enables us to provide comprehensive theoretical analysis on the impact
of critical system parameters and power scaling laws. Then, to characterize
the exact performance in practical systems, the ideal model is extended to the
practical one with spatial-correlated Rician fading and EMI. In this context,
we further analyze the impact of spatial correlation and EMI on the achievable
rate and on the power scaling laws. Finally, a gradient ascent method is
proposed to solve the minimal user rate maximization problem which only relies
on the statistical CSI. The specific contributions of this paper are
summarized as follows.

  * •

To begin with, we utilize the spatial-independent Rician fading model for the
sake of obtaining useful insights. The aggregated channel is estimated by
relying on linear minimum mean square error (LMMSE) method and then the
performance of mean square error (MSE) and normalized MSE (NMSE) is analyzed.
Next, applying MRC detectors, we derive the closed-form expressions for the
use-and-then-forget (UatF) bound of the achievable rate. The derived results
hold for arbitrary numbers of BS antennas and RIS elements. Then, we analyze
the impact of crucial system parameters, asymptotic behaviors of the rate, and
the power scaling laws. We also gain useful insights on the special case of a
single user.

  * •

To gain a general understanding of practical RIS-aided systems, we further
extend our model to include the presence of spatial correlation at the RIS and
the EMI captured by the RIS. We obtain the LMMSE channel estimates and derive
the UatF bound of the achievable rate in a closed form. We reveal that spatial
correlation makes RIS possess a higher capacity to customize the wireless
environment and prove that severe EMI might break the power scaling laws.

  * •

For both the spatial-independent and spatial-correlated models, we propose a
novel gradient ascent method to solve the minimal user rate maximization
problem. We first apply a log-sum-exp approximation to obtain a smooth
objective function. Then, we derive the real gradient vectors with respect to
real angle vectors. The performance loss in the projection is avoided since
the objective function is periodic with angles and the unit modulus constraint
holds for all angles. Besides, in the special case of a single user, we obtain
the closed-form solutions which substantiate our power scaling law analysis in
the previous section.

  * •

Numerical results validate our analytical conclusions and demonstrate that the
insights based on ideal models achieve high-accuracy approximations for those
in practical models when the spatial correlation and the EMI are light.
Besides, our results imply that the RIS with large numbers of elements could
bring more gains in the presence of spatial correlation. We also find that
when the EMI is severe, the RIS-aided systems could be inferior to
conventional massive MIMO systems. Furthermore, the gain of integrating RIS
into massive MIMO systems is confirmed especially when the RIS is deployed
close to some cell-edge users.

The remainder of this paper is organized as follows. The performance analysis
based on spatial-independent channels without EMI is carried out in Section
II, III, and IV. Specifically, Section II proposes the system model. Section
III derives the LMMSE channel estimates and evaluates the channel estimation
quality. Section IV derives closed-form lower bound of the achievable rate and
analyze the corresponding power scaling laws. The extension to spatial-
correlated channels in the presence of EMI is studied in Section V. Section VI
proposes the gradient ascent-based algorithm for solving the minimal user rate
maximization problem. Section VII provides extensive numerical results and
Section VIII concludes the paper.

_Notations_ : Vectors and matrices are denoted by boldface lower case and
upper case letters, respectively. The transpose, conjugate, conjugate
transpose, and inverse of matrix $\bf X$ are denoted by ${\bf X}^{T}$, ${\bf
X}^{*}$, ${\bf X}^{H}$ and ${\bf X}^{-1}$, respectively. $\left[{\bf
X}\right]_{m,n}$ denotes the $(m,n)$th entry of matrix $\bf X$. The real,
imaginary, trace, expectation, and covariance operators are denoted by ${\rm
Re}\left\\{\cdot\right\\}$, ${\rm Im}\left\\{\cdot\right\\}$, ${\rm
Tr}\left\\{\cdot\right\\}$, ${\mathbb{E}}\left\\{\cdot\right\\}$, and ${\rm
Cov}\left\\{\cdot\right\\}$, respectively. The $l_{2}$ norm of a vector and
the absolute value of a complex number are denoted by $\left\|\cdot\right\|$
and $\left|\cdot\right|$, respectively. $\mathbb{C}^{M\times N}$ denotes the
space of $M\times N$ complex matrices. ${\bf I}_{M}$ and $\bf 0$ denote the
$M\times M$ identity matrix and all-zero matrix with appropriate dimension,
respectively. The operator $\bmod$ returns the remainder after division, and
$\lfloor x\rfloor$ denotes the nearest integer smaller than $x$. ${\bf
x}\sim\mathcal{CN}\left({\bar{\bf x}},{\bf C}\right)$ is a complex Gaussian
distributed vector with mean $\bar{\bf x}$ and covariance matrix ${\bf C}$.
$\mathcal{O}$ denotes the standard big-O notation. Besides, for ease of
reference, the main symbols used in this work are listed in Table I.

TABLE I: List of Main Symbols Symbol | Definition | Symbol | Definition
---|---|---|---
$M$/$N$/$K$ | Number of BS antennas/RIS elements/users | $p$ | Transmit power for each user
$\theta_{n}$ | Phase shifts of the $n$-th element | $\bm{\theta}$ | Angle vector equal to $\left[{\theta_{1}},{\theta_{2}},\ldots,{\theta_{N}}\right]^{T}$
$\bm{c}$ | Vector equal to $e^{j\bm{\theta}}$ | $\bf\Phi$ | RIS phase shifts matrix, $\bm{\Phi}=\operatorname{diag}(\mathbf{c})$
$\sigma^{2}$/$\sigma_{e}^{2}$, $\rho$ | Power of thermal noise/EMI, $\rho=\frac{\sigma_{e}^{2}}{\sigma^{2}}$ |  $\mathbf{x}$/$\mathbf{n}$/$\bm{v}$ | Signal/noise/EMI vector
$d_{ris}$/$d_{bs}$ | Element spacing of RIS/BS | $\lambda$ | Wavelength
$\tau$/$\tau_{c}$ | Lengths of pilot signal/coherence interval |  $\mathbf{s}_{k}$, $\mathbf{S}$ | User $k$’s pilot sequence, $\mathbf{S}=\left[\mathbf{s}_{1},\mathbf{s}_{2},\ldots,\mathbf{s}_{K}\right]$
$\mathbf{N}$/$\mathbf{V}$ | Noise/EMI collecting in $\tau$ time slots | $\gamma_{k}$ | pathloss of user $k$’s direct link
$\alpha_{k}$ | Pathloss of user $k$-RIS link | $\beta$ | pathloss of RIS-BS link
$\delta$ | Rician factor of RIS-BS link | $\varepsilon_{k}$ | Rician factor of user $k$-RIS link
$\mathbf{d}_{k}$, $\tilde{\mathbf{d}}_{k}$ | User $k$-BS direct link, $\mathbf{d}_{k}=\sqrt{\gamma_{k}}\tilde{\mathbf{d}}_{k}$ |  $\mathbf{h}_{k}$, $\overline{\mathbf{h}}_{k}$, $\tilde{\mathbf{h}}_{k}$ | User $k$-RIS link, comprised of $\overline{\mathbf{h}}_{k}$ and $\tilde{\mathbf{h}}_{k}$
$\mathbf{H}_{2}$/$\mathbf{H}_{c,2}$ | RIS-BS link without/with correlation |  $\tilde{\mathbf{H}}_{2}$/$\tilde{\mathbf{H}}_{c,2}$ | NLoS part of $\mathbf{H}_{2}$/$\mathbf{H}_{c,2}$
$\mathbf{q}_{k}$/$\mathbf{q}_{c,k}$ | Aggregated link without/with correlation |  $\mathbf{Q}$/$\mathbf{Q}_{c}$ | Matrix with the $k$-th column of $\mathbf{q}_{k}$/$\mathbf{q}_{c,k}$
$\hat{\mathbf{q}}_{k}$/$\hat{\mathbf{q}}_{c,k}$ | Channel estimate of $\mathbf{q}_{k}$/$\mathbf{q}_{c,k}$ |  $\hat{\mathbf{Q}}$/$\hat{\mathbf{Q}}_{c}$ | Matrix with the $k$-th column of $\hat{\mathbf{q}}_{k}$/$\hat{\mathbf{q}}_{c,k}$
$\mathbf{q}_{k}^{1}$-$\mathbf{q}_{k}^{4}$, $\underline{\mathbf{q}}_{k}$ | Notations defined in (II-A) |  $\hat{\mathbf{q}}_{k}^{1}$-$\hat{\mathbf{q}}_{k}^{4}$, $\underline{\hat{\mathbf{q}}}_{k}$ | Notations defined in (17)
$\mathbf{y}$/$\mathbf{y}_{c}$ | BS received signal without/with correlation |  $\mathbf{r}$/$\mathbf{r}_{c}$ | Decoded symbols from $\mathbf{y}$/$\mathbf{y}_{c}$
$\mathbf{Y}_{P}$/$\mathbf{Y}_{c,P}$ | Received pilot signals at the BS |  $\mathbf{y}_{p}^{k}$/$\mathbf{y}_{c,p}^{k}$ | Observation vector without/with correlation
$\mathbf{a}_{M}\\!(.)$/$\mathbf{a}_{N}\\!(.)$ | Array response vector for BS/RIS | $E_{u}$ | A constant used in power scaling laws
$c_{k}$, $\widehat{c}_{k}$ |  $c_{k}=\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}$, $\widehat{c}_{k}=\frac{\beta\alpha_{k}}{\delta+1}$ |  $a_{k1}$ \- $a_{k4}$ | Notations defined in Lemma 1 and Theorem 1
$e_{k1}$ \- $e_{k3}$ | Notations defined in Lemma 2 | $f_{k}(\bm{\Phi})$ | Scalar equal to $\mathbf{a}_{N}^{H}\bm{\Phi}\overline{\mathbf{h}}_{k}$
$\mathbf{A}_{k}$, $\mathbf{B}_{k}$ | Matrices defined in Theorem 1 | $\bm{\Upsilon}_{k}$ | Matrix defined in Theorem 3
$\underline{R}_{k}$/$\underline{R}_{c,k}$ | Use $k$’s rate without/with correlation |  $f(\bm{\theta})$/$f_{c}(\bm{\theta})$ | Approximated minimal user rate
$\mathbf{R}_{ris}$, $\\!\\!\mathbf{R}_{emi}$ | Spatial correlation matrices |  $\bm{{f}}_{d}(.)$, $\mathbf{z}_{k}(.)$ | Function defined in Lemma 4, 5
$f_{c,1}(\mathbf{\Phi})$, $f_{c,k,2}(\mathbf{\Phi})$ \- $f_{c,k,7}(\mathbf{\Phi})$, $f_{c,ki,8}(\mathbf{\Phi})$ \- $f_{c,ki,9}(\mathbf{\Phi})$ | Scalar functions defined in (95)
$\bm{f}_{c,1}^{\prime}(\bm{\theta})$, $\bm{f}_{c,k,2}^{\prime}(\bm{\theta})$ \- $\bm{f}_{c,k,7}^{\prime}(\bm{\theta})$, $\bm{f}_{c,ki,8}^{\prime}(\bm{\theta})$ \- $\bm{f}_{c,ki,9}^{\prime}(\bm{\theta})$ | Gradient vectors defined in Lemma 6
$E_{k}^{\text{signal }}$, $I_{ki}$, $E_{k}^{\text{leak }}$, $E_{k}^{\text{noise }}$ | Signal, interference, leakage, and noise in Theorem 2
$E_{c,k}^{\text{signal }}$, $I_{c,ki}$, $E_{c,k}^{\text{leak }}$, $E_{c,k}^{\text{emi }}$, $E_{c,k}^{\text{noise }}$ | Signal, interference, leakage, EMI and noise in Theorem 4

##  II System Model

To begin with, we study the RIS-aided massive MIMO systems by neglecting the
spatial correlation and the EMI. This is an ideal model but it could
effectively simplify the expressions, facilitate the analysis, and help us
gain some preliminary but useful insights. Although this ideal model will not
be valid in practical scenarios, our simulation results demonstrate that it
can be regarded as a high-accuracy approximation for the practical model when
the RIS element spacing is large (half-wavelength) and when the EMI is light.
Therefore, the conclusions and insights based on this ideal model remain
meaningful for understanding the properties and benefits of integrating RIS
into massive MIMO systems.

Figure 1: An RIS-aided massive MIMO system.

As illustrated in Fig. 1, we consider the UL transmission of an RIS-aided
massive MIMO system, where an RIS is deployed in the proximity of $K$ users to
assist their UL transmissions to the BS. For convenience, we denote the set of
users as $\mathcal{K}=\left\\{1,2,\ldots,K\right\\}$. The BS is equipped with
$M$ active antennas, the RIS comprises $N$ passive reflecting elements, and
the $K$ users are equipped with a single transmit antenna. The channels from
user $k$, $k\in\mathcal{K}$ to the BS, from user $k$, $k\in\mathcal{K}$ to the
RIS, and from the RIS to the BS are denoted by
$\mathbf{d}_{k}\in\mathbb{C}^{M\times 1}$,
$\mathbf{h}_{k}\in\mathbb{C}^{N\times 1}$, and
$\mathbf{H}_{2}\in\mathbb{C}^{M\times N}$, respectively. Additionally, we
define
$\mathbf{D}=\left[\mathbf{d}_{1},\mathbf{d}_{2},\ldots,\mathbf{d}_{K}\right]$
and
$\mathbf{H}_{1}=\left[\mathbf{h}_{1},\mathbf{h}_{2},\ldots,\mathbf{h}_{K}\right]$.

The RIS shapes the propagation environment by phase-shifting the impinging
signals. Its phase shift matrix is denoted by
$\mathbf{\Phi}=\mathrm{diag}\left\\{e^{j\theta_{1}},e^{j\theta_{2}},\ldots,e^{j\theta_{N}}\right\\}$,
where $\theta_{n}\in[0,2\pi)$ represents the phase shift of the $n$th
reflecting element. Based on these definitions, the cascaded user $k$-RIS-BS
channel can be written as $\mathbf{g}_{k}=\mathbf{H}_{2}{\bf\Phi}{\bf h}_{k}$,
and the cascaded channels of the $K$ users are collected in the matrix
$\mathbf{G}=[\mathbf{g}_{1},\mathbf{g}_{2},\ldots,\mathbf{g}_{K}]=\mathbf{H}_{2}{\bf\Phi}{\bf
H}_{1}\in\mathbb{C}^{M\times K}$.

The $K$ users transmit their data in the same UL time-frequency resource. For
ease of exposition, let $\mathbf{Q}=\mathbf{G}+\mathbf{D}=[{\bf q}_{1},{\bf
q}_{2},\ldots,{\bf q}_{K}]\in\mathbb{C}^{M\times K}$ denote the aggregated
instantaneous channel matrix from the users to the BS. Thereby, the signal
vector received at the BS is given by

$\displaystyle\mathbf{y}=\sqrt{p}{\bf
Q}\mathbf{x}+\mathbf{n}=\sqrt{p}\sum\nolimits_{k=1}^{K}{\bf
q}_{k}x_{k}+\mathbf{n},$ (1)

where $p$ is the average transmit power of each user,
$\mathbf{x}=[x_{1},x_{2},\ldots,x_{K}]^{T}$ are the transmit symbols of the
$K$ users, and $\mathbf{n}\sim\mathcal{CN}\left({\bf
0},\sigma^{2}\mathbf{I}_{M}\right)$ denotes the noise vector.

The BS applies a low-complexity MRC receiver to detect the transmitted
symbols. Before designing the MRC matrix, the channel $\bf Q$ has to be
estimated at the BS. A standard LMMSE estimator is employed to obtain the
estimated channel $\hat{\bf{Q}}$, as explained in the next section111Given the
LMMSE channel estimator, this work focuses on the study of new properties and
possible benefits after introducing the RIS into massive MIMO systems. It is
meaningful to investigate other channel estimators (such as least-squares and
element-wise MMSE[53, 59, 60]) and evaluate the trade-off between estimation
quality and implementation complexity. The comparison between different
channel estimation strategies will be left to our future work.. Relying on the
channel estimate, the BS performs MRC by multiplying the received signal $\bf
y$ with $\hat{\bf Q}^{H}$, as follows

$\displaystyle\mathbf{r}=\hat{\bf{Q}}^{H}\mathbf{y}=\sqrt{p}\hat{\bf{Q}}^{H}\mathbf{Q}\mathbf{x}+\hat{\bf{Q}}^{H}\mathbf{n}.$
(2)

Then, the $k$th element of the vector $\bf r$ can be expressed as

$\displaystyle
r_{k}=\sqrt{p}\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}x_{k}+\sqrt{p}\sum\nolimits_{i=1,i\neq
k}^{K}\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{i}x_{i}+\hat{\mathbf{q}}_{k}^{H}\mathbf{n},\;\;k\in\mathcal{K},$
(3)

where $\hat{\mathbf{q}}_{k}$ is the $k$th column of $\hat{\mathbf{Q}}$.

###  II-A Channel Model

Since the users may be located far away from the BS and a large number of
environmental blocking objects (i.e., blockages such as trees, vehicles,
buildings) may exist in the area of interest, the LoS path between the users
and the BS could be blocked. As in [39, 47, 48], we adopt the Rayleigh fading
model to describe the NLoS channel between the user $k$ and the BS, as follows

$\displaystyle\mathbf{d}_{k}=\sqrt{\gamma_{k}}\tilde{\mathbf{d}}_{k},\;\;k\in\mathcal{K},$
(4)

where $\gamma_{k}$ denotes the distance-dependent path-loss, and
$\tilde{\mathbf{d}}_{k}$ denotes the fast fading NLoS channel. The entries of
$\tilde{\mathbf{d}}_{k}$ are independent and identically distributed (i.i.d.)
complex Gaussian random variables, i.e.,
$\tilde{\mathbf{d}}_{k}\sim\mathcal{CN}\left({\bf 0},\mathbf{I}_{M}\right)$.

Considering that the RIS is often installed on the facades of high-rise
buildings and it could be placed near the users, the channels between the
users and the RIS have a high LoS probability. In addition, the RIS and the BS
are usually deployed at some heights above the ground, which implies that LoS
paths are likely to exist between the RIS and the BS. Therefore, as in [39,
47, 48, 49, 50], we adopt the Rician fading model for the user-RIS and RIS-BS
channels, as follows

$\displaystyle\mathbf{h}_{k}=\sqrt{\frac{\alpha_{k}}{\varepsilon_{k}+1}}\left(\sqrt{\varepsilon_{k}}\,\overline{\mathbf{h}}_{k}+\tilde{\mathbf{h}}_{k}\right),k\in\mathcal{K},$
(5)
$\displaystyle\mathbf{H}_{2}=\sqrt{\frac{\beta}{\delta+1}}\left(\sqrt{\delta}\,\overline{\mathbf{H}}_{2}+\tilde{\mathbf{H}}_{2}\right),$
(6)

where $\alpha_{k}$ and $\beta$ represent the path-loss coefficients,
$\varepsilon_{k}$ and $\delta$ are the Rician factors that account for the
ratio of the LoS power to the NLoS power of the corresponding propagation
paths. Furthermore, $\overline{\mathbf{h}}_{k}$ and
$\overline{\mathbf{H}}_{2}$ denote the LoS components, whereas
$\tilde{\mathbf{h}}_{k}$ and $\tilde{\mathbf{H}}_{2}$ represent the NLoS
components. For the NLoS paths, the components of $\tilde{\mathbf{h}}_{k}$ and
$\tilde{\mathbf{H}}_{2}$ are i.i.d. complex Gaussian random variables with
zero mean and unit variance. For the LoS paths, uniform linear array (ULA) and
uniform squared planar array (USPA) models are adopted for the BS and the RIS,
respectively. Hence, $\overline{\mathbf{h}}_{k}$ and
$\overline{\mathbf{H}}_{2}$ are, respectively, modelled as follows

$\displaystyle{\overline{\bf h}_{k}}={\bf
a}_{N}\left({\varphi_{kr}^{a},\varphi_{kr}^{e}}\right),k\in\mathcal{K},$ (7)
$\displaystyle{\overline{\bf H}_{2}}={\bf
a}_{M}\left({\phi_{r}^{a},\phi_{r}^{e}}\right){\bf
a}_{N}^{H}\left({\varphi_{t}^{a},\varphi_{t}^{e}}\right),$ (8)

where $\varphi_{kr}^{a}$ ($\varphi_{kr}^{e}$) is the azimuth (elevation) angle
of arrival (AoA) of the incident signal at the RIS from the user $k$,
$\varphi_{t}^{a}$ ($\varphi_{t}^{e}$) is the azimuth (elevation) angle of
departure (AoD) reflected by the RIS towards the BS, and $\phi_{r}^{a}$
($\phi_{r}^{e}$) is the azimuth (elevation) AoA of the signal received at the
BS from the RIS, respectively. Furthermore, ${{\bf
a}_{X}}\left({\vartheta^{a},\vartheta^{e}}\right)\in\mathbb{C}^{X\times 1}$
denotes the array response vector, whose $x$-th entry is

$\displaystyle\begin{aligned} &\left[{{\bf
a}_{M}}\left({\vartheta^{a},\vartheta^{e}}\right)\right]_{x}=\exp\left\\{j2\pi\frac{d_{bs}}{\lambda}(x-1)\sin\vartheta^{e}\sin\vartheta^{a}\right\\},\\\
&\left[{{\bf
a}_{N}}\left({\vartheta^{a},\vartheta^{e}}\right)\right]_{x}=\exp\left\\{j2\pi\frac{d_{ris}}{\lambda}\left({\lfloor\left({x}-1\right)/\sqrt{N}\rfloor\sin\vartheta^{e}\sin\vartheta^{a}+\left(\left({x}-1\right)\bmod\sqrt{N}\right)\cos\vartheta^{e}}\right)\right\\},\end{aligned}$
(9)

where $d_{bs}$, $d_{ris}$, and $\lambda$ denote the BS antenna spacing, RIS
element spacing, and the wavelength, respectively.

To simplify the notation, in the sequel, we denote ${\bf
a}_{N}\left({\varphi_{kr}^{a},\varphi_{kr}^{e}}\right)$, ${\bf
a}_{M}\left({\phi_{r}^{a},\phi_{r}^{e}}\right)$, and ${\bf
a}_{N}\left({\varphi_{t}^{a},\varphi_{t}^{e}}\right)$ simply by
$\overline{\mathbf{h}}_{k}$, ${\bf a}_{M}$, and ${\bf a}_{N}$, respectively.
Then, the aggregated channel from the user $k$ to the BS can be expressed as

$\displaystyle\mathbf{q}_{k}$
$\displaystyle=\mathbf{g}_{k}+\mathbf{d}_{k}=\mathbf{H}_{2}{\bf\Phi}{\bf
h}_{k}+\mathbf{d}_{k}$
$\displaystyle={\underbrace{\sqrt{c_{k}\delta\varepsilon_{k}}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}}_{\mathbf{q}_{k}^{1}}+\underbrace{\sqrt{c_{k}\delta}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\tilde{\mathbf{h}}_{k}}_{\mathbf{q}_{k}^{2}}+\underbrace{\sqrt{c_{k}\varepsilon_{k}}\,\tilde{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}}_{\mathbf{q}_{k}^{3}}+\underbrace{\sqrt{c_{k}}\,\tilde{\mathbf{H}}_{2}\mathbf{\Phi}\tilde{\mathbf{h}}_{k}}_{\mathbf{q}_{k}^{4}}}+\sqrt{\gamma_{k}}\,\tilde{\mathbf{d}}_{k}$
$\displaystyle\triangleq{\underline{\mathbf{q}}}_{k}+{\mathbf{d}}_{k},$ (10)

where
$c_{k}\triangleq\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}$,
and ${\underline{\mathbf{q}}}_{k}=\sum_{\omega=1}^{4}{\bf q}^{\omega}_{k}$.
Note that $\underline{\mathbf{q}}_{k}$ and ${\mathbf{d}}_{k}$ are mutually
independent.

##  III Channel Estimation

In this section, we use the LMMSE method to obtain the estimated aggregated
instantaneous channel $\hat{\bf Q}$. Specifically, the BS estimates the
aggregated channel matrix $\bf Q$ based on some predefined pilot signals. Let
$\tau_{c}$ and $\tau$ denote the length of the channel coherence interval and
the number of time slots used for channel estimation, respectively, where
$\tau$ is no smaller than $K$, i.e., $\tau\geq K$. In each channel coherence
interval, the $K$ users simultaneously transmit mutually orthogonal pilot
sequences to the BS. The pilot sequence of the user $k$ is denoted by ${\bf
s}_{k}\in\mathbb{C}^{\tau\times 1}$. By defining ${\bf S}=[{\bf s}_{1},{\bf
s}_{2},\ldots,{\bf s}_{K}]$, we have ${\bf S}^{H}{\bf S}={\bf I}_{K}$. Then,
the $M\times\tau$ pilot signals received at the BS can be written as

$\displaystyle\mathbf{Y}_{p}=\sqrt{\tau p}\mathbf{Q}{\bf S}^{H}+\mathbf{N},$
(11)

where $\tau p$ is the transmit pilot power, and $\bf N$ denotes the
$M\times\tau$ noise matrix whose entries are i.i.d. complex Gaussian random
variables with zero mean and variance $\sigma^{2}$. Multiplying (11) by
$\frac{{\bf s}_{k}}{\sqrt{\tau p}}$ and exploiting the orthogonality of the
pilot signals, the BS obtains the following observation vector for user $k$

$\displaystyle\mathbf{y}_{p}^{k}=\frac{1}{\sqrt{\tau
p}}\mathbf{Y}_{p}\mathbf{s}_{k}=\mathbf{q}_{k}+\frac{1}{\sqrt{\tau
p}}\mathbf{Ns}_{k}.$ (12)

The optimal estimate of the $k$-th user’s channel based on the observation
vector $\mathbf{y}_{p}^{k}$ can be determined based on the MMSE criterion,
which has been widely utilized in conventional massive MIMO systems [61, 62,
63]. In RIS-aided massive MIMO systems where Rician fading is considered for
all RIS-related channels, however, it is challenging to obtain the MMSE
estimator. This is because the cascaded user-RIS-BS channel $\bf G$ in RIS-
aided systems is not Gaussian distributed, but double Gaussian distributed
instead[64]. To obtain closed-form channel estimates, as is needed for an
insightful analysis, we adopt the sub-optimal but tractable LMMSE estimator.
This is because the LMMSE estimator only requires the knowledge of first and
second order statistics, and therefore does not need to know the exact channel
distributions. In the following lemma, we present the required statistics for
the channel vector $\mathbf{q}_{k}$ and the observation vector
$\mathbf{y}_{p}^{k}$.

######  Lemma 1.

For $k\in\mathcal{K}$, the mean vectors and covariance matrices needed for the
LMMSE estimator are given by222These long-term statistical properties (such as
geometric pathloss, Rician factors, and angles) can remain unchanged for a
long time in the sub-6 GHz bands. Therefore, they can be estimated readily
given a system setup and then be assumed to be perfectly known for brevity[53,
54, 40]. The investigation of imperfect statistical CSI will be left to our
future work.

$\displaystyle\mathbb{E}\left\\{\mathbf{q}_{k}\right\\}=\mathbb{E}\left\\{\mathbf{y}^{k}_{p}\right\\}=\sqrt{c_{k}\delta\varepsilon_{k}}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k},$
(13)
$\displaystyle\operatorname{Cov}\left\\{\mathbf{q}_{k},\mathbf{y}_{p}^{k}\right\\}=\operatorname{Cov}\left\\{\mathbf{y}_{p}^{k},\mathbf{q}_{k}\right\\}=\operatorname{Cov}\left\\{\mathbf{q}_{k},\mathbf{q}_{k}\right\\}=a_{k1}\mathbf{a}_{M}\mathbf{a}_{M}^{H}+a_{k2}\mathbf{I}_{M},$
(14)
$\displaystyle\operatorname{Cov}\left\\{\mathbf{y}_{p}^{k},\mathbf{y}_{p}^{k}\right\\}=\operatorname{Cov}\left\\{\mathbf{q}_{k},\mathbf{q}_{k}\right\\}+\frac{\sigma^{2}}{\tau
p}\mathbf{I}_{M}=a_{k1}\mathbf{a}_{M}\mathbf{a}_{M}^{H}+\left(a_{k2}+\frac{\sigma^{2}}{\tau
p}\right)\mathbf{I}_{M},$ (15)

where $a_{k1}=Nc_{k}\delta$ and
$a_{k2}=Nc_{k}\left(\varepsilon_{k}+1\right)+\gamma_{k}$ are two auxiliary
variables.

Proof: See Appendix B. $\blacksquare$

######  Theorem 1.

Using the observation vector ${\mathbf{y}}^{k}_{p}$, the LMMSE estimate
$\hat{\mathbf{q}}_{k}$ of the channel vector ${\mathbf{q}}_{k}$ is given by
333Note that $\hat{\mathbf{q}}_{k}^{1}$ and ${\mathbf{q}}_{k}^{1}$ are
actually the same due to the unbiased estimation. However, we define two
symbols in order to simplify the form of expressions and make the derivations
easier to understand (see (223) and (373) for example).

$\displaystyle\hat{\mathbf{q}}_{k}=\mathbf{A}_{k}\mathbf{y}_{p}^{k}+\mathbf{B}_{k}$
(16)
$\displaystyle\quad=\underbrace{\underbrace{\sqrt{c_{k}\delta\varepsilon_{k}}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}}_{\hat{\mathbf{q}}_{k}^{1}}+\underbrace{\left(Ma_{k3}+a_{k4}\right)\sqrt{c_{k}\delta}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\tilde{\mathbf{h}}_{k}}_{\hat{\mathbf{q}}_{k}^{2}}+\underbrace{\sqrt{c_{k}\varepsilon_{k}}\mathbf{A}_{k}\tilde{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}}_{\hat{\mathbf{q}}_{k}^{3}}+\underbrace{\sqrt{c_{k}}\mathbf{A}_{k}\tilde{\mathbf{H}}_{2}\mathbf{\Phi}\tilde{\mathbf{h}}_{k}}_{\hat{\mathbf{q}}_{k}^{4}}}_{\hat{\underline{\mathbf{q}}}_{k}}$
$\displaystyle\quad\quad+\sqrt{\gamma_{k}}\mathbf{A}_{k}\tilde{\mathbf{d}}_{k}+\frac{1}{\sqrt{\tau
p}}\mathbf{A}_{k}\mathbf{Ns}_{k},$ (17)

where

$\displaystyle\mathbf{A}_{k}=\mathbf{A}^{H}_{k}=a_{k3}\mathbf{a}_{M}\mathbf{a}_{M}^{H}+a_{k4}\mathbf{I}_{M},$
(18)
$\displaystyle\mathbf{B}_{k}=\left(\mathbf{I}_{M}-\mathbf{A}_{k}\right)\sqrt{c_{k}\delta\varepsilon_{k}}\,\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k},$
(19) $\displaystyle a_{k3}=\frac{a_{k1}\frac{\sigma^{2}}{\tau
p}}{\left(a_{k2}+\frac{\sigma^{2}}{\tau
p}\right)\left(a_{k2}+\frac{\sigma^{2}}{\tau p}+Ma_{k1}\right)},$ (20)
$\displaystyle a_{k4}=\frac{a_{k2}}{a_{k2}+\frac{\sigma^{2}}{\tau p}},$ (21)

and the NMSE of the estimate of $\mathbf{q}_{k}$ is

$\displaystyle\operatorname{NMSE}_{k}=\frac{\operatorname{Tr}\left\\{\operatorname{Cov}\left\\{{\bf
q}_{k}-\hat{\bf q}_{k},{\bf q}_{k}-\hat{\bf
q}_{k}\right\\}\right\\}}{\operatorname{Tr}\left\\{\operatorname{Cov}\left\\{\mathbf{q}_{k},\mathbf{q}_{k}\right\\}\right\\}}=\frac{\frac{\sigma^{2}}{\tau
p}\left(Ma_{k1}a_{k2}+a_{k2}^{2}+\left(a_{k1}+a_{k2}\right)\frac{\sigma^{2}}{\tau
p}\right)}{\left(a_{k2}+\frac{\sigma^{2}}{\tau
p}\right)\left(a_{k2}+\frac{\sigma^{2}}{\tau
p}+Ma_{k1}\right)\left(a_{k1}+a_{k2}\right)}.$ (22)

Proof: See Appendix C. $\blacksquare$

As is evident from Theorem 1, we only estimate the aggregated channel matrix
${\bf Q}\in\mathbb{C}^{M\times K}$ including the reflected and direct
channels, which has the same dimension as the user-BS channel matrix in
conventional massive MIMO systems. Therefore, we only require that the length
of the pilot sequences is no smaller than the number of users, i.e., $\tau\geq
K$. Compared to estimation methods that estimate the $MN$ individual channels
in RIS-aided communications[36, 37], the proposed method has lower overhead
and computational complexity.

######  Remark 1.

When $c_{k}=0,\forall k$, i.e., the RIS-assisted channels are absent, we have
$a_{k1}=0$, $a_{k2}=\gamma_{k}$, $a_{k3}=0$,
$a_{k4}=\frac{\gamma_{k}}{\gamma_{k}+\frac{\sigma^{2}}{\tau p}}$ and ${\bf
B}_{k}={\bf 0}$. In this case, the estimate in (16) reduces to
$\hat{\mathbf{q}}_{k}=\frac{\gamma_{k}}{\gamma_{k}+\frac{\sigma^{2}}{\tau
p}}\mathbf{y}_{p}^{k}$ and the MSE matrix in (178) reduces to ${\rm\bf
MSE}_{k}=\frac{\gamma_{k}\frac{\sigma^{2}}{\tau
p}}{\gamma_{k}+\frac{\sigma^{2}}{\tau p}}{\bf I}_{M}$, which, as expected, is
the same as the MSE in conventional massive MIMO systems [62]. If the RIS
channels only have the LoS components, i.e.,
$\delta,\varepsilon_{k}\to\infty,\forall k$, we can also obtain $a_{k1}\to 0$
and $a_{k2}\to\gamma_{k}$. In this case, the MSE matrix in (178) is again the
same as that in conventional massive MIMO systems. This is because the LoS
channels are deterministic and known, and, thus, they do not bring additional
estimation errors.

######  Corollary 1.

In the low pilot power-to-noise ratio regime, high pilot power-to-noise ratio
regime, and large $N$ regime, the asymptotic NMSE is, respectively, given by

$\displaystyle\lim\nolimits_{\frac{\sigma^{2}}{\tau
p}\rightarrow\infty}\mathrm{NMSE}_{k}\rightarrow 1,$ (23)
$\displaystyle\lim\nolimits_{\frac{\sigma^{2}}{\tau p}\rightarrow
0}\mathrm{NMSE}_{k}\rightarrow 0,$ (24)
$\displaystyle\lim\nolimits_{N\rightarrow\infty}\mathrm{NMSE}_{k}\rightarrow
0.$ (25)

Besides, consider that the power $p$ is scaled proportionally to $p=E_{u}/N$,
where $E_{u}$ denotes a constant. As $N\to\infty$, we have

$\displaystyle\lim\nolimits_{p=\frac{E_{u}}{N},N\rightarrow\infty}\operatorname{NMSE}_{k}<1.$
(26)

Proof: When $\frac{\sigma^{2}}{\tau p}\rightarrow\infty$ or
$N\rightarrow\infty$, by selecting the dominant terms in (22), which scale
with ${(\frac{\sigma^{2}}{\tau p})}^{2}$ or $N^{3}$, we arrive at (23) and
(25), respectively. Substituting $\frac{\sigma^{2}}{\tau p}=0$ into (22), its
numerator reduces to zero, which leads to (24). Replacing the power $p$ in
(22) with $p=E_{u}/N$, as $N\to\infty$, we can readily find that all the
dominant terms in the numerator are present in the denominator as well, which
results in (26). We omit the specific limit of (26) since it is a complex
expression but is simple to compute. $\blacksquare$

It is worth noting that NMSE values between $0$ (i.e., perfect estimation) and
$1$ (i.e., using the mean value of the variable as the estimate) quantify the
relative estimation error [53]. In conventional massive MIMO systems, a common
method for reducing the NMSE is to increase the length of the pilot sequence
$\tau$. In RIS-aided massive MIMO systems, Corollary 1 indicates that
increasing the number of RIS elements $N$ can play a similar role as
increasing $\tau$. Therefore, increasing the number of RIS elements not only
helps improve the system capacity, but also helps reduce the NMSE.
Additionally, (26) reveals that an RIS equipped with a large number of
reflecting elements $N$ can help the NMSE converge to a limit lower than one,
even for low pilot powers.

To better understand the impact of increasing $N$ for channel estimation, we
present the following asymptotic results.

######  Corollary 2.

When $\tau\to\infty$, we have $\hat{\mathbf{q}}_{k}\to{\mathbf{q}}_{k}$, which
implies that ${\bf e}_{k}\to\bf 0$ and therefore the channel estimation is
perfect. When $N\to\infty$, by contrast, we have

$\displaystyle\hat{\mathbf{q}}_{k}\to{\mathbf{q}}_{k}+\frac{1}{\sqrt{\tau
p}}\mathbf{Ns}_{k},$ (27) $\displaystyle{\bf e}_{k}={\bf q}_{k}-\hat{\bf
q}_{k}\to\frac{-1}{\sqrt{\tau p}}\mathbf{Ns}_{k},$ (28) $\displaystyle{\rm\bf
MSE}_{k}=\mathbb{E}\left\\{{\bf e}_{k}{\bf
e}_{k}^{H}\right\\}\to\frac{\sigma^{2}}{\tau p}{\bf I}_{M}.$ (29)

Proof: When $\tau\to\infty$ or $N\to\infty$, based on Theorem 1, we have
$a_{k3}\to 0$, $a_{k4}\to 1$, and $Ma_{k3}+a_{k4}\to 1$, which yields ${\bf
A}_{k}\to{\bf I}_{M}$. If $\tau\to\infty$, we further get $\frac{1}{\sqrt{\tau
p}}\to 0$, which completes the proof. $\blacksquare$

Although the NMSE converges to zero as $N\to\infty$ (see (25)), Corollary 2
shows that, in contrast to increasing $\tau$, the MSE of the LMMSE estimator
converges to a non-zero constant as $N\to\infty$. If we estimate the channel
$\mathbf{q}_{k}$ based on the least-squares (LS) estimator[53, (3.35)], it is
interesting to note that we will obtain the same results as in (27) and (29).
In general, the LS estimator, which does not exploit any prior channel
statistics, has worse estimation performance (higher MSE) than the LMMSE
estimator[37, 59, 53]. Therefore, Corollary 2 indicates that the MSE
performance of the LMMSE estimation converges towards an upper bound, which is
the MSE performance of the LS estimation, as $N\to\infty$. This result will be
validated in Section VII.

######  Corollary 3.

When the RIS-BS channel reduces to the purely Rayleigh channel (i.e.,
$\delta=0$), the estimated channel vector, MSE, and NMSE, respectively, reduce
to

$\displaystyle{\hat{\bf
q}}_{k}=a_{k4}\mathbf{I}_{M}\mathbf{y}^{k}_{p}=\frac{N\beta\alpha_{k}+\gamma_{k}}{N\beta\alpha_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}}\mathbf{y}^{k}_{p},$ (30)
$\displaystyle{\rm\bf{MSE}}_{k}=\frac{\left(N\beta\alpha_{k}+\gamma_{k}\right)\frac{\sigma^{2}}{\tau
p}}{N\beta\alpha_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau p}}\mathbf{I}_{M},$
(31) $\displaystyle\operatorname{NMSE}_{k}=\frac{\frac{\sigma^{2}}{\tau
p}}{N\beta\alpha_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau p}}.$ (32)

Proof: When $\delta=0$, we have $a_{k1}=0$,
$a_{k2}=N\beta\alpha_{k}+\gamma_{k}$, $a_{k3}=0$,
$a_{k4}=\frac{N\beta\alpha_{k}+\gamma_{k}}{N\beta\alpha_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}}$ and ${\bf B}_{k}={\bf 0}$. The proof follows by inserting these results
in Theorem 1 and (178) . $\blacksquare$

Corollary 3 corresponds to a scenario where a large number of scatterers exist
nearby the RIS and the BS, and the LoS path between the RIS and the BS is
negligible. Therefore, the RIS-BS channel is dominated by the NLoS paths. In
this case, both the MSE and NMSE have simple analytical expressions, which
help us better understand the conclusions drawn in Corollary 1 and Corollary
2. It is apparent that the MSE (represented by the trace of ${\rm\bf MSE}_{k}$
in (31)) and the NMSE (represented by ${\rm NMSE}_{k}$ in (32)) are decreasing
functions of the pilot power $\tau p$. As a function of $N$, on the other
hand, the MSE is an increasing function, while the NMSE is a decreasing
function. When $N\to\infty$, we have ${\rm\bf
MSE}_{k}\to\frac{\sigma^{2}}{\tau p}{\bf I}_{M}$ but ${\rm NMSE}_{k}\to{0}$.
Note that we can obtain the MSE and NMSE for conventional massive MIMO systems
by setting $N=0$ in (31) and (32). Therefore, the obtained result implies that
the MSE of RIS-aided massive MIMO systems is worse than the MSE of massive
MIMO systems without RISs, while the NMSE of RIS-aided massive MIMO systems is
better than the NMSE of massive MIMO systems without RISs. The reason is that
an RIS introduces $N$ additional paths to the system, but the pilot length
$\tau$ does not increase correspondingly, which increases the estimation
error. However, the presence of an RIS results in better channel gains, which
help decrease the normalized error.

Furthermore, if we reduce the power as $p=E_{u}/N$, as $N\to\infty$, the NMSE
in (32) converges to a limit less than one, as follows

$\displaystyle\lim\nolimits_{\delta=0,\,p=\frac{E_{u}}{N},\,N\rightarrow\infty}\operatorname{NMSE}_{k}\to\frac{{\sigma^{2}}}{\tau
E_{u}\beta\alpha_{k}+{\sigma^{2}}}<1.$ (33)

##  IV Analysis of the Achievable Rate

Based on the channel estimates provided in Theorem 1, closed-form expressions
of the lower bound of the achievable rate are derived and analyzed in this
section444To avoid verbose expressions, the term “lower bound of the
achievable rate” will be abbreviated as “achievable rate” in the sequel of
this paper.. In Section VI, the obtained analytical expressions are utilized
for optimizing the phase shifts of the RIS based on statistical CSI.

###  IV-A Derivation of the Rate

As in [65, 66, 59, 67], we utilize the so called UatF bound, which is a
tractable lower bound, to characterize the ergodic rate of RIS-aided massive
MIMO systems. First, we rewrite $r_{k}$ in (3) as

$\displaystyle
r_{k}=\underbrace{\sqrt{p}\,\mathbb{E}\left\\{\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}\right\\}x_{k}}_{\text{Desired
signal
}}+\underbrace{\sqrt{p}\left(\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}-\mathbb{E}\left\\{\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}\right\\}\right)x_{k}}_{\text{Signal
leakage }}+\underbrace{\sqrt{p}\sum\nolimits_{i=1,i\neq
k}^{K}\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{i}x_{i}}_{\text{Multi-user
interference }}+\underbrace{\hat{\mathbf{q}}_{k}^{H}\mathbf{n}}_{\text{Noise
}}.$ (34)

Then, we formulate the lower bound of the $k$-th user’s ergodic rate as
$\underline{R}_{k}=\frac{\tau_{c}-\tau}{\tau_{c}}\log_{2}\left(1+\mathrm{SINR}_{k}\right)$,
where the pre-log factor $\frac{\tau_{c}-\tau}{\tau_{c}}$ represents the rate
loss that originates from the pilot overhead, and the SINR is expressed as

$\displaystyle\mathrm{SINR}_{k}=\frac{p\left|\mathbb{E}\left\\{\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}\right\\}\right|^{2}}{p\left(\mathbb{E}\left\\{\left|\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}\right|^{2}\right\\}-\left|\mathbb{E}\left\\{\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{k}\right\\}\right|^{2}\right)+p\sum\limits_{i=1,i\neq
k}^{K}\mathbb{E}\left\\{\left|\hat{\mathbf{q}}_{k}^{H}\mathbf{q}_{i}\right|^{2}\right\\}+\sigma^{2}\mathbb{E}\left\\{\left\|\hat{\mathbf{q}}_{k}\right\|^{2}\right\\}}.$
(35)

To simplify the expression of $\underline{R}_{k}$, we define three auxiliary
variables $e_{k1}$, $e_{k2}$, and $e_{k3}$. These variables capture the
performance degradation due to the imperfect knowledge of the CSI.

######  Lemma 2.

For $k\in\mathcal{K}$, we have ${\rm Tr}\left\\{{\bf A}_{k}\right\\}=Me_{k1}$,
${\bf A}_{k}\overline{\bf H}_{2}=e_{k2}\overline{\bf H}_{2}$ and ${\rm
Tr}\left\\{{\bf A}_{k}{\bf A}_{k}\right\\}=Me_{k3}$, where

$\displaystyle e_{k1}\triangleq a_{k3}+a_{k4},$ (36) $\displaystyle
e_{k2}\triangleq Ma_{k3}+a_{k4},$ (37) $\displaystyle e_{k3}\triangleq
Ma_{k3}^{2}+2a_{k3}a_{k4}+a_{k4}^{2}.$ (38)

Furthermore, $e_{k1},e_{k2}$ and $e_{k3}$ are bounded in $\left[0,1\right]$.
When $\tau p\to\infty$ or $N\to\infty$, we have $e_{k1},e_{k2},e_{k3}\to 1$.
When $\tau p\to 0$, by contrast, we have $e_{k1},e_{k2},e_{k3}\to 0$.

Proof: See Appendix D. $\blacksquare$

In the following theorem, we derive a closed-form expression for the
achievable rate.

######  Theorem 2.

A lower bound for the ergodic rate of the $k$-th user is given by555The phase
shift matrix $\bf\Phi$ is assumed to be fixed when deriving the achievable
rate. After obtaining the achievable rate expressions, we would design
$\bf\Phi$ so that the derived rate is optimized.

$\displaystyle\begin{aligned}
&\underline{R}_{k}=\tau^{o}\log_{2}\left(1+\mathrm{SINR}_{k}\right),\\\
&\mathrm{SINR}_{k}=\frac{pE_{k}^{\rm
signal}\left(\mathbf{\Phi}\right)}{pE_{k}^{\rm
leak}\left(\mathbf{\Phi}\right)+p\sum\limits_{i=1,i\neq
k}^{K}I_{ki}\left(\mathbf{\Phi}\right)+\sigma^{2}E_{k}^{\rm{noise}}\left(\mathbf{\Phi}\right)},\end{aligned}$
(39)

where $\tau^{o}=\frac{\tau_{c}-\tau}{\tau_{c}}$, $E_{k}^{\rm
signal}\left(\mathbf{\Phi}\right)=\left\\{E_{k}^{\rm{noise}}\left(\mathbf{\Phi}\right)\right\\}^{2}$,

$\displaystyle
E_{k}^{\rm{noise}}\left(\mathbf{\Phi}\right)=M\left\\{\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}\delta\varepsilon_{k}+Nc_{k}\delta
e_{k2}+\left(Nc_{k}\left(\varepsilon_{k}+1\right)+\gamma_{k}\right)e_{k1}\right\\},\qquad\quad\quad\quad\quad$
(40) $\displaystyle\begin{array}[]{l}E_{k}^{\rm
leak}\left(\mathbf{\Phi}\right)=M\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}^{2}\delta\varepsilon_{k}\left\\{N\left(M\delta+\varepsilon_{k}+1\right)\left(e_{k2}^{2}+1\right)+2\left(Me_{k1}+e_{k2}\right)\left(e_{k2}+1\right)\right\\}\\\
\qquad\qquad\;\;\quad+M\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}\delta\varepsilon_{k}\left\\{\gamma_{k}+\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k2}^{2}\right\\}\\\
\qquad\qquad\;\;\quad+M^{2}N^{2}c_{k}^{2}\delta^{2}e_{k2}^{2}\\\
\qquad\qquad\;\;\quad+MN^{2}c_{k}^{2}\left\\{2\delta\left(\varepsilon_{k}+1\right)e_{k2}^{2}+\left(\varepsilon_{k}+1\right)^{2}e_{k3}\right\\}\\\
\qquad\qquad\;\;\quad+M^{2}Nc_{k}^{2}\left\\{\left(2\varepsilon_{k}+1\right)e_{k1}^{2}+2\delta
e_{k1}e_{k2}\right\\}\\\
\qquad\qquad\;\;\quad+MNc_{k}\left\\{c_{k}\left(2\delta
e_{k2}^{2}+\left(2\varepsilon_{k}+1\right)e_{k3}\right)+\left(2\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)\left(\delta
e_{k2}^{2}+\left(\varepsilon_{k}+1\right)e_{k3}\right)\right\\}\\\
\qquad\qquad\;\;\quad+M\gamma_{k}\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k3},\end{array}$ (48)

and

$\displaystyle\begin{array}[]{l}I_{ki}\left(\mathbf{\Phi}\right)=M^{2}\left|f_{k}(\mathbf{\Phi})\right|^{2}\left|f_{i}(\mathbf{\Phi})\right|^{2}c_{k}c_{i}\delta^{2}\varepsilon_{k}\varepsilon_{i}\\\
\qquad\quad+M\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}\delta\varepsilon_{k}\left\\{c_{i}\left(MN\delta+N\varepsilon_{i}+N+2Me_{k1}\right)+\gamma_{i}\right\\}\\\
\qquad\quad+M\left|f_{i}(\mathbf{\Phi})\right|^{2}c_{i}\delta\varepsilon_{i}\left\\{c_{k}e_{k2}\left(MN\delta
e_{k2}+N\varepsilon_{k}e_{k2}+Ne_{k2}+2Me_{k1}\right)+\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k2}^{2}\right\\}\\\
\qquad\quad+M^{2}N^{2}c_{k}c_{i}\delta^{2}e_{k2}^{2}\\\
\qquad\quad+MN^{2}c_{k}c_{i}\left\\{\delta\left(\varepsilon_{k}+\varepsilon_{i}+2\right)e_{k2}^{2}+\left(\varepsilon_{k}+1\right)\left(\varepsilon_{i}+1\right)e_{k3}\right\\}\\\
\qquad\quad+M^{2}Nc_{k}c_{i}e_{k1}\left\\{\left(\varepsilon_{k}+\varepsilon_{i}+1\right)e_{k1}+2\delta
e_{k2}\right\\}\\\
\qquad\quad+M^{2}c_{k}c_{i}\varepsilon_{k}\varepsilon_{i}e_{k1}\left(\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}e_{k1}+2\delta\operatorname{Re}\left\\{f_{k}^{H}(\mathbf{\Phi})f_{i}(\mathbf{\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\}\right)\\\
\qquad\quad+MN\left\\{\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)c_{i}\left(\delta
e_{k2}^{2}+\left(\varepsilon_{i}+1\right)e_{k3}\right)+\gamma_{i}c_{k}\left(\delta
e_{k2}^{2}+\left(\varepsilon_{k}+1\right)e_{k3}\right)\right\\}\\\
\qquad\quad+M\gamma_{i}\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k3},\end{array}$ (58)

with

$\displaystyle f_{k}({\bf\Phi})\triangleq$
$\displaystyle\mathbf{a}_{N}^{H}{\bf\Phi}\overline{\mathbf{h}}_{k}=\sum\nolimits_{n=1}^{N}e^{j\left(\zeta_{n}^{k}+\theta_{n}\right)},$
(59) $\displaystyle\zeta_{n}^{k}=$ $\displaystyle
2\pi\frac{d}{\lambda}\left(\lfloor(n-1)/\sqrt{N}\rfloor\left(\sin\varphi_{kr}^{e}\sin\varphi_{kr}^{a}-\sin\varphi_{t}^{e}\sin\varphi_{t}^{a}\right)\right.$
$\displaystyle\left.+((n-1)\bmod\sqrt{N})\left(\cos\varphi_{kr}^{e}-\cos\varphi_{t}^{e}\right)\right).$
(60)

Proof: See Appendix E. $\blacksquare$

The closed-form expression in Theorem 2 does not involve the calculation of
inverse matrices and the numerical computation of integrals. In contrast to
time-consuming Monte Carlo simulations, the evaluation of the rate based on
Theorem 2 has a low computational complexity even if $M$ and $N$ are large
numbers, as usually is in RIS-aided massive MIMO systems. More importantly,
Theorem 2 only relies on statistical CSI, i.e., the distance-based path loss
factors, the AoA and AoD, and the environment-related Rician factors, which
change slowly in the sub-6 GHz frequency band. By using the analytical
expression of the rate in (39) as an objective function for system design, we
are able to optimize the phase shifts of the RIS only based on long-term
statistical CSI.

By comparing the formulation in Theorem 2 with that given in [49, Theorem 1],
we can see that the impact of imperfect CSI is completely characterized by the
parameters $e_{k1},e_{k2},e_{k3}$ and $\frac{\sigma^{2}}{\tau p}$. In the
perfect CSI scenario, we can assume $\tau=\infty$, which leads to
$e_{k1}=e_{k2}=e_{k3}=1$ and $\frac{\sigma^{2}}{\tau p}=0$. Based on Theorem
2, we can analyze the rate performance of RIS-aided massive MIMO systems for
arbitrary system parameters. Even though these expressions are long, it would
be clear when focusing on the coefficients of $M$, $N$, and the variables of
$f_{k}(\mathbf{\Phi})$, $\forall k$. These are the key parameters since they
decide the order of magnitudes of signal, interference, signal leakage, and
the noise with respect to the number of BS antennas and RIS elements. For
example, since the interference term $I_{ki}$ scales on the order of
$\mathcal{O}(M^{2})$, it is easy to know that the RIS-aided massive MIMO
systems suffer from stronger multi-user interference than conventional massive
MIMO systems. In the following, we will provide comprehensive analysis of the
RIS-aided massive MIMO systems, including the asymptotic behavior of the rate
for large values of $M$ and $N$, the power scaling laws, and the impact of the
Rician factors. To this end, we begin with a useful lemma.

######  Lemma 3.

  * •

If $N=1$, for arbitrary $\bf\Phi$, we have
$\left|f_{k}(\mathbf{\Phi})\right|=1$ in (59).

  * •

If $N>1$, by optimizing $\bf\Phi$, the range of values
$0\leq\left|f_{k}(\mathbf{\Phi})\right|\leq N$ is achievable in (59).

  * •

If we configure the phase shifts of the RIS to achieve
$\left|f_{k}(\mathbf{\Phi})\right|=N$, unless the user $i,i\neq k$, has the
same azimuth and elevation AoA as the user $k$, the function
$\left|f_{i}(\mathbf{\Phi})\right|$ in (59) is bounded when $N\to\infty$.

  * •

Unless the user $i,i\neq k$, has the same azimuth and elevation AoA as the
user $k$, the term
$\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}$ is
bounded when $N\to\infty$.

Proof: See Appendix F. $\blacksquare$

###  IV-B Multi-user Case

In this section, we consider the general multi-user scenario, i.e., $K>1$.
Since any two users are unlikely to be in the same location, we assume that
the azimuth and elevation AoA of any two users are different, i.e.,
$\left({\varphi_{kr}^{a},\varphi_{kr}^{e}}\right)\neq\left({\varphi_{ir}^{a},\varphi_{ir}^{e}}\right)$.
To begin with, we investigate the asymptotic behavior of the rate in (39) for
large values of $M$ and $N$.

######  Remark 2.

From Theorem 2, we observe that, as a function of $M$, $E_{k}^{\rm
signal}\left(\mathbf{\Phi}\right)$, $E_{k}^{\rm
leak}\left(\mathbf{\Phi}\right)$ and $I_{ki}\left(\mathbf{\Phi}\right)$ behave
asymptotically as $\mathcal{O}\left(M^{2}\right)$. Therefore, the rate
$\underline{R}_{k}$ converges to a finite limit when $M\to\infty$. If, on the
other hand, we align the phase shifts of the RIS for maximizing the intended
signal for the user $k$, i.e., we set $\left|f_{k}(\mathbf{\Phi})\right|=N$,
then we have $\underline{R}_{k}\to\infty$ for user $k$, and
$\underline{R}_{i}\to 0$ for the other users $i\neq k$ as $N\to\infty$, based
on Lemma 3. In a multi-user scenario, this implies that it is necessary to
enforce some fairness among the users when designing the phase shifts of the
RIS.

Next, we study the power scaling laws of RIS-aided massive MIMO systems with
different Rician factors. Note that the Rician factor characterizes the fading
severity of the environment and the richness of scatterers in the environment.
The smaller the Rician factor is, the larger the number of scatterers in the
environment is. If the Rician factor is zero, we retrieve the Rayleigh fading
channel as a special case in which only the NLoS components exist. If the
Rician factor tends to infinity, the channel is deterministic and it is
characterized only by the LoS component. It is worth mentioning that, under
the assumption of imperfect CSI, decreasing the transmit power $p$ results in
a reduction of the power used for both the data and pilot signals.

We will analyze several scenarios for the RIS-BS and user-RIS channels. For
ease of exposition, we summarize the obtained power scaling laws as a function
of $M$ and $N$ in Table II.

TABLE II: Power scaling laws in the multi-user case. | $\left(\text{RIS-BS channel, user-RIS channels}\right)$
---|---
| $\left(\text{Rician, Rician}\right)$ | $\left(\text{Rician, Rayleigh}\right)$ | $\left(\text{Rayleigh, Rician}\right)$ | $\left(\text{Rayleigh, Rayleigh}\right)$
Imperfect CSI | $M$ | ${1}/{M}$ | ${1}/{M}$ | ${1}/{\sqrt{M}}$ | ${1}/{\sqrt{M}}$
$N$ | $\diagdown$ | $1/N$
Perfect CSI | $M$ | $1/M$
$N$ | $\diagdown$ | $1/N$

Specifically, the following notations are used. “Imperfect CSI” and “Perfect
CSI” are referred to the power scaling laws obtained for imperfect and perfect
CSI, respectively. By setting $e_{k1}=e_{k2}=e_{k3}=1$ and
$\frac{\sigma^{2}}{\tau p}=0$, which are obtained from $\tau\to\infty$, the
imperfect CSI setup reduces to the perfect CSI setup. The notation
“$\left(\text{Rician, Rician}\right)$” means that the RIS-BS channel and all
the user-RIS channels are Rician distributed, i.e., $\delta>0$ and
$\varepsilon_{k}>0,\forall k$. Similarly, the notation “$\left(\text{Rician,
Rayleigh}\right)$” means that the RIS-BS channel is Rician distributed and all
the user-RIS channels are Rayleigh distributed, i.e., $\delta>0$ and
$\varepsilon_{k}=0,\forall k$. The notations “$1/M$”, “$1/\sqrt{M}$” and
“$1/N$” imply that the rate tends to a non-zero value if the transmit power
scales down proportionally to $1/M$, $1/\sqrt{M}$ and $1/N$, respectively. We
mention, for completeness, that the readers interested in the power scaling
laws as a function of $M$ in conventional massive MIMO systems without RISs
may refer to [61] and [62]. Besides, we stress that the rate does not depend
on the RIS phase shift matrix $\bf\Phi$ if $\delta=0$ or
$\varepsilon_{k}=0,\forall k$, which will be proved in Corollary 5. In the
following, we mainly focus on the proof for the imperfect CSI case, since the
perfect CSI setup can be obtained in a similar manner, by setting
$e_{k1}=e_{k2}=e_{k3}=1$ and $\frac{\sigma^{2}}{\tau p}=0$.

######  Corollary 4.

(“$1/M$” for “(Rician, Rician)” and “(Rician, Rayleigh)”) Assume that the
transmit power $p$ is scaled as $p={E_{u}}/{M}$. For $M\to\infty$, the rate of
the user $k$, $k\in\mathcal{K}$, is lower bounded by

$\displaystyle\underline{R}_{k}\rightarrow\tau^{o}\log_{2}\left(1+\frac{E_{u}c_{k}^{2}\delta^{2}\left(\left|f_{k}(\mathbf{\Phi})\right|^{2}\varepsilon_{k}+Ne_{k2}\right)^{2}}{E_{u}E_{k}^{\rm
leak}\left(\mathbf{\Phi}\right)+E_{u}\sum\limits_{i=1,i\neq
k}^{K}I_{ki}\left(\mathbf{\Phi}\right)+\sigma^{2}c_{k}\delta\left(\left|f_{k}(\mathbf{\Phi})\right|^{2}\varepsilon_{k}+Ne_{k2}\right)}\right),$
(61)

where

$\displaystyle E_{k}^{\rm{leak}}\left(\mathbf{\Phi}\right)=$ $\displaystyle
N\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}^{2}\delta^{2}\varepsilon_{k}\left(e_{k2}^{2}+1\right)+\frac{\sigma^{2}}{\tau
E_{u}}\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}\delta\varepsilon_{k}e_{k2}^{2}$
$\displaystyle+N^{2}c_{k}^{2}\delta^{2}e_{k2}^{2}+\frac{\sigma^{2}}{\tau
E_{u}}Nc_{k}\delta e_{k2}^{2},$ (62) $\displaystyle
I_{ki}\left(\mathbf{\Phi}\right)=$
$\displaystyle\left|f_{k}(\mathbf{\Phi})\right|^{2}\left|f_{i}(\mathbf{\Phi})\right|^{2}c_{k}c_{i}\delta^{2}\varepsilon_{k}\varepsilon_{i}+N\left|f_{k}(\mathbf{\Phi})\right|^{2}c_{k}c_{i}\delta^{2}\varepsilon_{k}$
$\displaystyle+\left|f_{i}(\mathbf{\Phi})\right|^{2}c_{i}\delta\varepsilon_{i}e_{k2}^{2}\left(Nc_{k}\delta+\frac{\sigma^{2}}{\tau
E_{u}}\right)+N^{2}c_{k}c_{i}\delta^{2}e_{k2}^{2}+N\frac{\sigma^{2}}{\tau
E_{u}}c_{i}\delta e_{k2}^{2},$ (63) $\displaystyle e_{k2}=$
$\displaystyle\frac{Nc_{k}\delta}{\frac{\sigma^{2}}{\tau
E_{u}}+Nc_{k}\delta}.$ (64)

Proof: If $p=E_{u}/M$ and $M\to\infty$, we have $e_{k1}\to 0$, $e_{k3}\to 0$,
and $e_{k2}$ tends to (64). The proof is completed by substituting $p=E_{u}/M$
into Theorem 2 and retaining the non-zero terms whose asymptotic behavior is
$\mathcal{O}\left(M\right)$. $\blacksquare$

For a massive number of antennas, Corollary 4 shows that the rate of all the
users tends to a non-zero value when the transmit power scales as $p=E_{u}/M$.
From (61), we evince that the rate $\underline{R}_{k}$ is still non-zero if
$\varepsilon_{k}=0,\forall k$, i.e., all the user-RIS channels are Rayleigh
distributed. This proves the power scaling law “$1/M$” for the
“$\left(\text{Rician, Rayleigh}\right)$” setup in Table II. However, the rate
$\underline{R}_{k}$ in (61) reduces to zero if $c_{k}=0$ or $\delta=0$, i.e.,
the RIS-aided channels are absent or the RIS-BS channel is Rayleigh
distributed. This indicates that the power scaling law “$1/M$” does not hold
for these two case studies. Specifically, the considered system degenerates to
an RIS-free massive MIMO system with Rayleigh fading if $c_{k}=0,\forall k$.
In this case, it has been proven that the rate can maintain a non-zero value
when the power scales as $p=E_{u}/\sqrt{M}$ [62, (37)]. As for the power
scaling law for $\delta=0$, we first provide the rate expression when
$\delta=0$.

######  Corollary 5.

If the RIS-BS channel is Rayleigh distributed ($\delta=0$), the rate of user
$k$, $k\in\mathcal{K}$, is lower bounded by

$\displaystyle\underline{R}_{k}^{(\rm{NL_{1}})}=\tau^{o}\log_{2}\left(1+\frac{pE_{k}^{\mathrm{signal}}}{pE_{k}^{\mathrm{leak}}+p\sum_{i=1,i\neq
k}^{K}I_{ki}+\sigma^{2}E_{k}^{\mathrm{noise}}}\right),$ (65)

where

$\displaystyle E_{k}^{\rm{signal}}=$ $\displaystyle
M\left(Nc_{k}\left(\varepsilon_{k}+1\right)+\gamma_{k}\right)^{2}e_{k1},$ (66)
$\displaystyle E_{k}^{\rm{noise}}=$ $\displaystyle
Nc_{k}\left(\varepsilon_{k}+1\right)+\gamma_{k},$ (67) $\displaystyle
E_{k}^{\rm{leak}}=$ $\displaystyle
N^{2}c_{k}^{2}\left(\varepsilon_{k}+1\right)^{2}e_{k1}+MNc_{k}^{2}\left(2\varepsilon_{k}+1\right)e_{k1}$
$\displaystyle+Nc_{k}\left\\{c_{k}\left(2\varepsilon_{k}+1\right)+\left(2\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)\left(\varepsilon_{k}+1\right)\right\\}e_{k1}+\gamma_{k}\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k1},$ (68) $\displaystyle I_{ki}=$ $\displaystyle
N^{2}c_{k}c_{i}\left(\varepsilon_{k}+1\right)\left(\varepsilon_{i}+1\right)e_{k1}+MNc_{k}c_{i}\left(\varepsilon_{k}+\varepsilon_{i}+1\right)e_{k1}+Mc_{k}c_{i}\varepsilon_{k}\varepsilon_{i}\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}e_{k1}$
$\displaystyle+N\left\\{\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)c_{i}\left(\varepsilon_{i}+1\right)+\gamma_{i}c_{k}\left(\varepsilon_{k}+1\right)\right\\}e_{k1}+\gamma_{i}\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)e_{k1},$ (69)

and

$\displaystyle
e_{k1}=\frac{N\beta\alpha_{k}+\gamma_{k}}{N\beta\alpha_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}}.$ (70)

Proof: When $\delta=0$, we have $a_{k1}=0$,
$a_{k2}=N\beta\alpha_{k}+\gamma_{k}$, and $a_{k3}=0$. Then, we obtain
$e_{k3}=e_{k1}^{2}$, where $e_{k1}=a_{k4}$ is given in (70). Substituting
$\delta=0$ into Theorem 2 and using $e_{k3}=e_{k1}^{2}$, the proof follows
with the aid of some algebraic simplifications. $\blacksquare$

It is observed that the rate expression in Corollary 5 does not depend on
$\bf\Phi$. Therefore, with a purely NLoS RIS-BS channel, any RIS phase shift
matrix will result in the same ergodic rate. This is because the RIS phase
shift matrix $\bf\Phi$ is a unitary matrix and the entries of the NLoS channel
$\tilde{\bf H}_{2}$ are Gaussian distributed. Therefore, $\tilde{\bf
H}_{2}{\bf\Phi}$ has the same statistical properties as $\tilde{\bf H}_{2}$.
Likewise, there is no need to design the RIS phase shifts if all the user-RIS
links are purely NLoS. This conclusion is apparent from (39) by setting
$\varepsilon_{k}=0,\forall k$.

By analyzing the dominant terms of (65) when $M,N\to\infty$, we evince that
the rate increases without bound for all the users. This implies that fairness
among the users is guaranteed in this special case. As $N\to\infty$,
specifically, the dominant terms in (65) scale asymptotically as
$\mathcal{O}\left(N^{2}\right)$, and the rate converges to

$\displaystyle\underline{R}_{k}^{(\rm{NL_{1}})}$
$\displaystyle\to\tau^{o}\log_{2}\left(1+\frac{M\alpha_{k}}{\sum_{i=1}^{K}\alpha_{i}}\right),\text{
as }N\to\infty,$ (71)
$\displaystyle=\tau^{o}\log_{2}\left(1+{M}/{K}\right),\text{ if
}\alpha_{1}=\ldots=\alpha_{K}.$ (72)

From (71), we evince that the SINR,
$\frac{M\alpha_{k}}{\sum_{i=1}^{K}\alpha_{i}}$, does not depend on the pilot
power $\tau p$ and it increases linearly with $M$. Therefore, promising rate
performance can be achieved if $\delta=0$ and $N\to\infty$.

With the aid of Corollary 5, we investigate, in the following corollaries, the
power scaling laws as a function of $M$ and $N$ when $\delta=0$.

######  Corollary 6.

(“$1/\sqrt{M}$” for “(Rayleigh, Rician)” and “(Rayleigh, Rayleigh)”) If the
RIS-BS channel is Rayleigh distributed ($\delta=0$), and the power is scaled
as $p=E_{u}/\sqrt{M}$ with $M\to\infty$, the rate of user $k$,
$k\in\mathcal{K}$ tends to
$\underline{R}_{k}^{(\mathrm{NL_{1}})}\rightarrow\tau^{o}\log_{2}\left(1+\mathrm{SINR}_{k}\right)$,
where the effective SINR is given by

$\displaystyle\operatorname{SINR}_{k}=\frac{\tau
E_{u}^{2}\left(Nc_{k}\left(\varepsilon_{k}+1\right)+\gamma_{k}\right)^{2}}{\tau
E_{u}^{2}Nc_{k}^{2}\left(2\varepsilon_{k}+1\right)+\sum_{i=1,i\neq k}^{K}\tau
E_{u}^{2}c_{k}c_{i}\left\\{N\left(\varepsilon_{k}+\varepsilon_{i}+1\right)+\varepsilon_{k}\varepsilon_{i}\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}\right\\}+\sigma^{4}}.$
(73)

Proof: First, we substitute $p=E_{u}/\sqrt{M}$ into Corollary 5 and ignore the
terms that tend to zero as $M\to\infty$. Then, we divide the numerator and
denominator of the SINR by $\frac{N\beta\alpha_{k}+\gamma_{k}}{\sigma^{2}}$.
This yields (73) and the proof is completed. $\blacksquare$

From (73), we evince that the numerator of the SINR scales with
$\mathcal{O}\left(N^{2}\right)$, but the denominator of the SINR only scales
with $\mathcal{O}\left(N\right)$. Therefore, Corollary 6 indicates that the
rate scales logarithmically with $N$ if $p=E_{u}/\sqrt{M}$ and $M\to\infty$,
which is a promising result for RIS-aided massive MIMO systems. Besides, it is
worth noting that (73) will reduce to the same expression as in [62, Eq. (37)]
when $c_{k}=0,\forall k$.

######  Corollary 7.

(“$1/N$” for “(Rayleigh, Rician)” and “(Rayleigh, Rayleigh)”) If the RIS-BS
channel is Rayleigh distributed ($\delta=0$) and the power is scaled as
$p=E_{u}/{N}$ with $N\to\infty$, the rate of user $k$, $k\in\mathcal{K}$, is
lower bounded by

$\displaystyle\underline{R}_{k}^{(\mathrm{NL_{1}})}\rightarrow\tau^{o}\log_{2}\left(1+\frac{E_{u}M\beta\alpha_{k}}{\sum_{i=1}^{K}\left(E_{u}\beta\alpha_{i}+\frac{\alpha_{i}}{\alpha_{k}}\frac{\sigma^{2}}{\tau}\right)+\sigma^{2}\left(1+\frac{\sigma^{2}}{\tau
E_{u}\beta\alpha_{k}}\right)}\right).$ (74)

Proof: First, we substitute $p=E_{u}/N$ into Corollary 5. When $N\to\infty$,
we have
$e_{k1}\rightarrow\frac{\beta\alpha_{k}}{\beta\alpha_{k}+\frac{\sigma^{2}}{\tau
E_{u}}}$. Then, we remove the non-dominant terms that do not scale as
$\mathcal{O}\left({N}\right)$. By noting that
$c_{k}\left(\varepsilon_{k}+1\right)=\beta\alpha_{k},\forall k$, and dividing
the numerator and denominator of the SINR by $\beta\alpha_{k}$, we obtain
(74). This completes the proof. $\blacksquare$

Corollary 7 sheds some interesting insights. Firstly, we note that Corollary 6
has unveiled that the transmit power $p$ can only be reduced proportionally to
$1/\sqrt{M}$, while maintaining a non-zero rate, when $\delta=0$. Corollary 7,
on the other hand, proves that the transmit power can be reduced
proportionally to $1/N$, while maintaining a non-zero rate, when $\delta=0$.
This reveals the positive role of deploying RISs in massive MIMO systems.
Secondly, the obtained power scaling law does not depend on the Rician factors
of user-RIS links, i.e., $\varepsilon_{k},\forall k$. This implies that the
rate in (74) is the same for LoS-only and NLoS-only user-RIS channels.
Thirdly, in (74), the desired signal term in (74) scales as $\mathcal{O}(M)$
and the interference term scales as $\mathcal{O}(1)$. As a result, the rate
scales logarithmically with the number of BS antennas. When the number of
antennas is large, the power of interference will be relatively small compared
with the power of desired signal, and then a promising rate can be achieved
with the setup in Corollary 7. Therefore, the rich-scattering environment
between the RIS and the BS ($\delta=0$) is beneficial in RIS-aided massive
MIMO systems, since it can provide sufficient spatial multiplexing gains and
help mitigate the multi-user interference. Finally, (74) unveils that, if the
users are all located at the same distance from the RIS, i.e.,
$\alpha_{1}=\ldots=\alpha_{K}$, they all achieve the same rate. Therefore,
fairness can be guaranteed in this special case.

Corollary 7 sheds light on the achievable rate when the RIS-BS channel is
Rayleigh distributed ($\delta=0$). In the next corollary, we analyze the
opposite scenario in which the user-RIS channels are Rayleigh distributed
($\varepsilon_{k}=0,\forall k$).

######  Corollary 8.

(“$1/N$” for “(Rician, Rayleigh)”) Assume $\delta>0$. If the user-RIS channels
are Rayleigh distributed ($\varepsilon_{k}=0,\forall k$) and the power is
scaled as $p=E_{u}/{N}$ with $N\to\infty$, the rate of user $k$,
$k\in\mathcal{K}$, is lower bounded by

$\displaystyle\underline{R}_{k}^{({\rm
NL_{2}})}\rightarrow\tau^{o}\log_{2}\left(1+\frac{E_{u}Mc_{k}^{2}\left(\delta
e_{k2}+e_{k1}\right)^{2}}{E_{u}\left(E_{k}^{\mathrm{leak}}+\sum_{i=1,i\neq
k}^{K}I_{ki}\right)+\sigma^{2}c_{k}\left(\delta e_{k2}+e_{k1}\right)}\right),$
(75)

with

$\displaystyle E_{k}^{\mathrm{leak}}+\sum_{i=1,i\neq
k}^{K}I_{ki}=\sum_{i=1}^{K}c_{i}\left\\{Mc_{k}\delta^{2}e_{k2}^{2}+c_{k}\left(2\delta
e_{k2}^{2}+e_{k3}\right)+\frac{\sigma^{2}}{\tau E_{u}}\left(\delta
e_{k2}^{2}+e_{k3}\right)\right\\},$ (76) $\displaystyle
a_{k3}=\frac{c_{k}\delta\frac{\sigma^{2}}{\tau
E_{u}}}{\left(c_{k}+\frac{\sigma^{2}}{\tau
E_{u}}\right)\left(c_{k}+\frac{\sigma^{2}}{\tau E_{u}}+Mc_{k}\delta\right)},$
(77) $\displaystyle a_{k4}=\frac{c_{k}}{c_{k}+\frac{\sigma^{2}}{\tau E_{u}}}.$
(78)

Proof: It follows from Theorem 2 by setting $\varepsilon_{k}=0,\forall k$ and
$p=E_{u}/N$, and by keeping only the dominant terms for $N\to\infty$.
$\blacksquare$

Corollary 8 characterizes the achievable rate when the user-RIS channels are
rich scattering. The obtained performance trends are very different from those
unveiled in Corollary 7 (i.e., the RIS-BS channel is rich scattering). In
contrast to Corollary 7, in particular, both the desired signal and the
interference in (75) scale as $\mathcal{O}\left(M\right)$. As a result, if the
user-RIS channels are Rayleigh distributed, the rate in (75) is still bounded
from above even if the number of BS antennas is very large. Besides, it is not
hard to prove that the rate (75) will reduce to the same expression as (74) if
we further set $\delta=0$. This result confirms our conclusion in Corollary 7
that the scaling law is not related with the Rician factor $\varepsilon_{k}$
if $\delta=0$.

From Corollary 7 and Corollary 8, we can find that small $\delta$ is
beneficial for the power scaling laws. This is because a small $\delta$
corresponds to a high-rank RIS-BS channel which provides sufficient spatial
diversity for the multi-user communication. It is known that due to the
product pathloss effect, it is better to deploy the RIS either close to BS or
close to the users so that the product pathloss effect is mitigated[68, 69].
Our finding further reveals that for each kind of deployment strategy, it is
also crucial to guarantee the enough spatial diversity between the BS and the
RIS. When the RIS is deployed close to the users, $\delta$ could be small
since the Rician factor commonly decreases with the communication
distance[70]. Therefore, placing the RIS close to the user is still a good
choice since it guarantees a high rank for the RIS-BS channel. If the RIS is
deployed near the BS, $\delta$ could be large and the RIS-BS channel could be
rank-deficient. In this context, some other operations are needed to improve
the rank such as introducing some deterministic scatterers between the BS and
the RIS or placing the RIS very close to the BS so that the spherical wave
model is valid[46].

###  IV-C Single-user Case

In this subsection, we analyze the power scaling laws in the special case with
only one user, i.e., $K=1$. Without loss of generality, the user is referred
to as user $k$. Since no other user exists, the rate can be obtained from
Theorem 2 by ignoring the multi-user interference term, i.e., by setting
$I_{ki}\left(\mathbf{\Phi}\right)=0$. For analytical tractability, we further
assume that the number of RIS elements is large. In this scenario (single-user
and large $N$), it can be proved that the optimal phase shift design that
maximizes the rate corresponds to the condition
$\left|f_{k}(\mathbf{\Phi})\right|=N$. This statement is formally proved in
the next section (Theorem 5).

Therefore, by setting $I_{ki}=0$ and $\left|f_{k}(\mathbf{\Phi})\right|=N$ in
Theorem 2, we obtain that the power of the desired signal scales as
$\mathcal{O}\left(M^{2}N^{4}\right)$, the power of the signal leakage scales
as $\mathcal{O}\left(M^{2}N^{3}\right)$, and the power of the noise term
scales as $\mathcal{O}\left(MN^{2}\right)$. Therefore, the rate is bounded for
$M\to\infty$, but it can grow without bound for $N\to\infty$. For ease of
exposition, similar to the multi-user case, we summarize the obtained power
scaling laws in Table III. In the following, we report the proofs only for
some (those that lead to insightful design guidelines) system setups that are
summarized in Table III. The proof of each case study can, in fact, be
obtained by using analytical steps similar to the multi-user case. Finally, we
mention that the power scaling laws in the single-user case with perfect CSI
can be derived readily based on [39, Eq. (17)].

TABLE III: Power scaling laws in the single-user case. | $\left(\text{RIS-BS channel, user-RIS channel}\right)$
---|---
| $\left(\text{Rician, Rician}\right)$ | $\left(\text{Rician, Rayleigh}\right)$ | $\left(\text{Rayleigh, Rician}\right)$ | $\left(\text{Rayleigh, Rayleigh}\right)$
Imperfect CSI | $M$ | ${1}/{M}$ | ${1}/{M}$ | ${1}/{\sqrt{M}}$ | ${1}/{\sqrt{M}}$
$N$ | $1/N^{2}$ | $1/N$
Perfect CSI | $M$ | $1/M$
$N$ | $1/N^{2}$ | $1/N$

######  Corollary 9.

Consider a single-user system with $\left|f_{k}(\mathbf{\Phi})\right|=N$. If
the transmit power is scaled as $p={E_{u}}/{\left(MN^{2}\right)}$ with
$M,N\to\infty$, the rate is lower bounded by

$\displaystyle\underline{R}_{k}\to\tau^{o}\log_{2}\left(1+\frac{E_{u}}{\sigma^{2}}\frac{\beta\alpha_{k}\delta\varepsilon_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right).$
(79)

If the transmit power is scaled as $p={E_{u}}/{N^{2}}$ with $N\to\infty$, the
rate is lower bounded by

$\displaystyle\underline{R}_{k}\rightarrow\tau^{o}\log_{2}\left(1+\frac{E_{u}}{\sigma^{2}}Mc_{k}\delta\varepsilon_{k}\right).$
(80)

Proof: Let us set $p={E_{u}}/{\left(MN^{2}\right)}$,
$\left|f_{k}(\mathbf{\Phi})\right|=N$ and $I_{ki}=0$ in Theorem 2. The rate in
(79) follows by using the fact $e_{k1},e_{k2},e_{k3}\to 0$ and by retaining
the dominant terms that scale as $\mathcal{O}\left(MN^{2}\right)$ for
$M,N\to\infty$. Next, let us set $p={E_{u}}/{N^{2}}$,
$\left|f_{k}(\mathbf{\Phi})\right|=N$ and $I_{ki}=0$ in Theorem 2. The rate in
(80) follows by retaining the dominant terms that scale as
$\mathcal{O}\left(N^{2}\right)$ for $N\to\infty$.

$\blacksquare$

The SNRs in (79) and (80) do not depend on $\tau$, and except for a pre-log
scaling factor, the same SNR as for perfect CSI-based systems can be derived
readily from [39, Eq. (17)]. We evince, therefore, that $\tau=K=1$ is the
optimal pilot length based on (79) and (80). Therefore, the overhead for
channel estimation is relatively low. Furthermore, the rates in (79) and (80)
are increasing functions of the Rician factors $\delta$ and $\varepsilon_{k}$,
which unveils that LoS-dominated environments are favorable for RIS-aided
single-user systems. If both $\delta\to\infty$ and $\varepsilon_{k}\to\infty$,
(79) and (80) are maximized. On the contrary, if $\delta=0$ or
$\varepsilon_{k}=0$, we observe that (79) and (80) tend to zero. This implies
that the power scaling law $1/N^{2}$ does not hold anymore. In these two
cases, the transmit power can be scaled only proportionally to $1/N$ for
maintaining a non-zero rate when $N\to\infty$. Mathematically, the
corresponding power scaling laws can be proved from Corollary 7 and Corollary
8 by setting the multi-user interference to zero. As an example, the case
study for $\delta=0$ is analyzed in the following corollary.

######  Corollary 10.

Consider a single-user system with $\delta=0$. If the transmit power is scaled
as $p={E_{u}}/{N}$ with $N\to\infty$, the rate is lower bounded by

$\displaystyle\underline{R}_{k}^{(\mathrm{NL_{1}})}\rightarrow\tau^{o}\log_{2}\left(1+\frac{E_{u}M\beta\alpha_{k}}{E_{u}\beta\alpha_{k}+\frac{\sigma^{2}}{\tau}+\sigma^{2}\left(1+\frac{\sigma^{2}}{\tau
E_{u}\beta\alpha_{k}}\right)}\right).$ (81)

As $\tau$ increases, the denominator of the SNR of (81) decreases. Therefore,
the SNR of (81) is an increasing function of $\tau$. Therefore, $\tau=1$ is
not guaranteed to be optimal in a rich-scattering environment ($\delta=0$),
and a relatively large number of pilots signals may be needed. Thus, Corollary
10 again emphasizes that LoS environments are favorable for RIS-aided single-
user systems.

##  V Extension to Correlated Channels with EMI

In this section, we begin to extend the ideal model to the practical one which
takes into the consideration of spatial correlation at the RIS and the EMI
captured by the RIS. This section has two purposes. Firstly, we want to
analyze the impact of spatial correlation and EMI on the achievable rate. This
could help us gain a more general understanding of RIS-aided massive MIMO
systems. Secondly, we want to utilize the practical results in this section to
examine whether the conclusions drawn based on the ideal model have reference
values for practical systems. Note that only the correlation at the RIS is
considered, while the correlation at the BS can be neglected due to the
assumption of ULA with half-wavelength antenna spacing. However, this
assumption never holds for the RIS since it is a planar and its diagonal
elements are always correlated [56].

###  V-A Channel Model with Spatial Correlation

Recall that our previous analysis indicates that it is a good choice to place
the RIS close to the users. In this context, LoS components could dominate the
user-RIS channels, and therefore the Rician factor $\varepsilon_{k}$ could be
relatively large. For the sake of brevity, this section only focuses on the
scenario where user-RIS channels are fully LoS (i.e.,
$\varepsilon_{k}\to\infty$, $\forall k$).666 Many contributions have revealed
that the change of the rate is marginal when the Rician factor is larger than
10[39, 40]. Thus, the considered scenario actually serves as a tractable
approximation of scenarios with relatively large $\varepsilon_{k}$. The
extension to the more general scenarios with arbitrary Rician fading factors
$\varepsilon_{k},\forall k$ can be done using methods proposed in Section IV
and this section. The extension will be left to our future work. In the
following, we will present the generalized system model after introducing
spatial correlation and EMI. Note that a letter $c$ will be added in the
subscript to indicate the existence of spatial correlation.

Considering spatial correlation and EMI, the received signal at the BS is
updated as

$\displaystyle\mathbf{y}_{c}=\sqrt{p}\mathbf{Q}_{c}\mathbf{x}+\mathbf{H}_{c,2}\mathbf{\Phi}\bm{v}+\mathbf{n},$
(82)

where
$\bm{v}\sim\mathcal{CN}\left(\mathbf{0},\sigma_{e}^{2}\mathbf{R}_{emi}\right)$
denotes the EMI received at the RIS with spatial correlation matrix
$\mathbf{R}_{emi}$. EMI is reflected by the RIS and then passes through RIS-BS
channel $\mathbf{H}_{c,2}$ leading to the term
$\mathbf{H}_{c,2}\mathbf{\Phi}\bm{v}$ in (82).
$\mathbf{Q}_{c}=[\mathbf{q}_{c,1},\mathbf{q}_{c,2},\ldots,\mathbf{q}_{c,K}]\in\mathbb{C}^{M\times
K}$ denotes the correlated aggregated channel from $K$ users to the BS, where
$\mathbf{q}_{c,k}=\mathbf{H}_{c,2}\bm{\Phi}\mathbf{h}_{k}+\mathbf{d}_{k}$ is
the aggregated channel of user $k$. User $k$-RIS channel $\mathbf{h}_{k}$ and
RIS-BS channel $\mathbf{H}_{c,2}$ are respectively given by

$\displaystyle\mathbf{h}_{k}=\sqrt{\alpha_{k}}\overline{\mathbf{h}}_{k},$ (83)
$\displaystyle\mathbf{H}_{c,2}=\sqrt{\frac{\beta}{\delta+1}}\left(\sqrt{\delta}\overline{\mathbf{H}}_{2}+\tilde{\mathbf{H}}_{c,2}\right),$
(84)

where $\tilde{\mathbf{H}}_{c,2}=\tilde{\mathbf{H}}_{2}\mathbf{R}_{ris}^{1/2}$
and $\mathbf{R}_{ris}$ denotes the spatial correlation matrix of NLoS channel
components at the RIS. Assuming isotropic scattering environment for $\bm{v}$
and $\tilde{\mathbf{H}}_{c,2}$, the spatial correlation matrices
$\mathbf{R}_{emi}$ and $\mathbf{R}_{ris}$ at the RIS can be characterized as
$\mathbf{R}_{emi}=\mathbf{R}_{ris}=\mathbf{R}$ with[56, 57]

$\displaystyle\left[\mathbf{R}\right]_{a,b}=\operatorname{sinc}\left(\frac{2\left\|\mathbf{u}_{a}-\mathbf{u}_{b}\right\|}{\lambda}\right),1\leq
a,b\leq N,$ (85)

where $\left\|\mathbf{u}_{a}-\mathbf{u}_{b}\right\|$ denotes the distance
between the $a$-th and $b$-th elements of the RIS, which depends on the RIS
element spacing $d_{ris}$. Since $\operatorname{sinc}(\cdot)$ is even, it is
readily to find that $\mathbf{R}=\mathbf{R}^{H}$. For brevity, we define
$\widehat{c}_{k}=\frac{\alpha_{k}\beta}{\delta+1}$. Then, based on (83) and
(84), the spatial-correlated aggregated channel of user $k$ can now be
expressed as

$\displaystyle\begin{aligned}
&\mathbf{q}_{c,k}=\mathbf{H}_{c,2}\bm{\Phi}\mathbf{h}_{k}+\mathbf{d}_{k}=\sqrt{\widehat{c}_{k}\delta}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{k}+\sqrt{\widehat{c}_{k}}\tilde{\mathbf{H}}_{c,2}\bm{\Phi}\overline{\mathbf{h}}_{k}+\sqrt{\gamma_{k}}\tilde{\mathbf{d}}_{k}.\end{aligned}$
(86)

###  V-B Channel Estimation

Next, we derive the LMMSE channel estimate $\hat{\mathbf{q}}_{c,k}$ for user
$k$’s aggregated channel. In the channel estimation phase, the BS receives the
$M\times\tau$ pilot signal as follows

$\displaystyle\mathbf{Y}_{c,P}=\sqrt{\tau
p}\mathbf{Q}_{c}\mathbf{S}^{H}+\mathbf{H}_{c,2}\mathbf{\Phi}\mathbf{V}+\mathbf{N},$
(87)

where
$\mathbf{V}=\mathbf{R}_{emi}^{1/2}\tilde{\mathbf{V}}\in\mathbb{C}^{N\times\tau}$
and each element of $\tilde{\mathbf{V}}\in\mathbb{C}^{N\times\tau}$ is
independently distributed as $\mathcal{CN}\left({0},\sigma_{e}^{2}\right)$.
After correlating $\mathbf{Y}_{c,P}$ with $\mathbf{s}_{k}$, the observation
vector for user $k$’s channel ${\mathbf{q}}_{c,k}$ is given by

$\displaystyle\begin{aligned} \mathbf{y}_{c,p}^{k}=\frac{1}{\sqrt{\tau
p}}\mathbf{Y}_{c,P}\mathbf{s}_{k}=\mathbf{q}_{c,k}+\frac{\left(\mathbf{H}_{c,2}\mathbf{\Phi}\mathbf{V}+\mathbf{N}\right)\mathbf{s}_{k}}{\sqrt{\tau
p}}.\end{aligned}$ (88)

######  Theorem 3.

Based on $\mathbf{y}_{c,p}^{k}$, the LMMSE channel estimate for
$\mathbf{q}_{c,k}$ is given by

$\displaystyle\hat{\mathbf{q}}_{c,k}=\sqrt{\widehat{c}_{k}\delta}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}+\sqrt{\widehat{c}_{k}}\mathbf{\Upsilon}_{k}\tilde{\mathbf{H}}_{c,2}\mathbf{\Phi}\overline{\mathbf{h}}_{k}+\sqrt{\gamma_{k}}\mathbf{\Upsilon}_{k}\tilde{\mathbf{d}}_{k}+\frac{\mathbf{\Upsilon}_{k}\mathbf{H}_{c,2}\bm{\Phi}\mathbf{V}\mathbf{s}_{k}}{\sqrt{\tau
p}}+\frac{\mathbf{\Upsilon}_{k}\mathbf{N}\mathbf{s}_{k}}{\sqrt{\tau p}},$ (89)

where

$\displaystyle\mathbf{\Upsilon}_{k}=\mathbf{\Upsilon}_{k}^{H}=\left(\widehat{c}_{k}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\mathbf{\Phi}\overline{\mathbf{h}}_{k}+\gamma_{k}\right)\times$
$\displaystyle\left\\{\left(\widehat{c}_{k}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}+\frac{\sigma_{e}^{2}\beta\operatorname{Tr}\left\\{\mathbf{R}_{emi}\mathbf{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\right\\}}{\tau
p(\delta+1)}\right)\mathbf{I}_{M}+\frac{\sigma_{e}^{2}\beta\delta\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}}{\tau
p(\delta+1)}\right\\}^{-1}.$ (90)

Proof: See Appendix G. $\blacksquare$

Besides, applying [71, Eq. (12.21)], the MSE matrix can be obtained as

$\displaystyle\textbf{MSE}_{c,k}=\left(\widehat{c}_{k}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\mathbf{\Phi}\overline{\mathbf{h}}_{k}+\gamma_{k}\right)\left(\mathbf{I}_{M}-\mathbf{\Upsilon}_{k}\right).$
(91)

(91) embodies the impacts of spatial correlation and EMI on channel
estimation. On one hand, MSE could be degraded by EMI power $\sigma_{e}^{2}$
through the term of $\mathbf{\Upsilon}_{k}$. On the other hand, it can be seen
that unitary matrices $\mathbf{\Phi}^{H}$ and $\mathbf{\Phi}$ are not canceled
by each other due to the existence of spatial correlation matrices
$\mathbf{R}_{ris}$ and $\mathbf{R}_{emi}$. By contrast, if the spatial
correlation is neglected ($\mathbf{R}_{ris}=\mathbf{R}_{emi}=\mathbf{I}_{N}$),
the MSE matrix in (91) will no longer depend on $\bf\Phi$, and therefore we
cannot exploit any benefits from optimizing the RIS to improve the quality of
channel estimation. Therefore, compared with the spatial independence case,
RIS with spatial correlation has a greater ability to improve the quality of
channel estimation.

###  V-C Achievable Rate

Based on the estimated channel $\hat{\mathbf{q}}_{c,k}$, we are able to design
the MRC detector and then derive the corresponding UatF bound of the
achievable rate in the presence of spatial correlation and EMI. Multiplying
MRC decoder matrix
$\hat{\mathbf{Q}}_{c}^{H}=\left[\hat{\mathbf{q}}_{c,1},\ldots,\hat{\mathbf{q}}_{c,K}\right]^{H}$
to the received signal $\mathbf{y}_{c}$ in (82), the decoded symbols at the BS
are given by

$\displaystyle\mathbf{r}_{c}=\hat{\mathbf{Q}}_{c}^{H}\mathbf{y}_{c}=\sqrt{p}\hat{\mathbf{Q}}_{c}^{H}\mathbf{Q}_{c}\mathbf{x}+\hat{\mathbf{Q}}_{c}^{H}\mathbf{H}_{c,2}\mathbf{\Phi}\bm{v}+\hat{\mathbf{Q}}_{c}^{H}\mathbf{n}.$
(92)

Then, the $k$-th entry of $\mathbf{r}_{c}$ can be expressed as follows

$\displaystyle\begin{aligned}
r_{c,k}&=\sqrt{p}\mathbb{E}\left\\{\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}\right\\}x_{k}+\sqrt{p}\left(\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}-\mathbb{E}\left\\{\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}\right\\}\right)x_{k}\\\
&+\sqrt{p}\sum_{i=1,i\neq
k}^{K}\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,i}x_{i}+\hat{\mathbf{q}}_{c,k}^{H}\mathbf{H}_{c,2}\mathbf{\Phi}\bm{v}+\hat{\mathbf{q}}_{c,k}^{H}\mathbf{n}.\end{aligned}$
(93)

Accordingly, the SINR of user $k$ can be written as

$\displaystyle\operatorname{SINR}_{c,k}=\frac{pE_{c,k}^{\mathrm{signal}}}{pE_{c,k}^{\mathrm{leak}}+p\sum_{i=1,i\neq
k}^{K}I_{c,ki}+\sigma_{e}^{2}E_{c,k}^{\mathrm{emi}}+\sigma^{2}E_{c,k}^{\mathrm{noise}}},$
(94)

with desired signal of
$E_{c,k}^{\mathrm{signal}}=\left|\mathbb{E}\left\\{\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}\right\\}\right|^{2}$,
signal leakage of
$E_{c,k}^{\mathrm{leak}}=\mathbb{E}\left\\{\left|\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}\right|^{2}\right\\}-\left|\mathbb{E}\left\\{\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,k}\right\\}\right|^{2}$,
interference of
$I_{c,ki}=\mathbb{E}\left\\{\left|\hat{\mathbf{q}}_{c,k}^{H}\mathbf{q}_{c,i}\right|^{2}\right\\}$,
EMI of
$E_{c,k}^{\mathrm{emi}}=\mathbb{E}\left\\{\hat{\mathbf{q}}_{c,k}^{H}\mathbf{H}_{c,2}\mathbf{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\mathbf{H}_{c,2}^{H}\hat{\mathbf{q}}_{c,k}\right\\}$,
and noise of
$E_{c,k}^{\mathrm{noise}}=\mathbb{E}\left\\{\left\|\hat{\mathbf{q}}_{c,k}\right\|^{2}\right\\}$.

In the following, we will derive the closed-form UatF bound of the achievable
rate. To this end, for the sake of brevity, for $1\leq k,i\leq K$, some terms
are denoted as follows

$\displaystyle\begin{aligned}
&f_{c,1}(\bm{\Phi})=\operatorname{Tr}\left\\{\mathbf{R}_{ris}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\right\\},\quad
f_{c,k,2}(\bm{\Phi})=\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k},\\\
&f_{c,k,3}(\bm{\Phi})=\operatorname{Tr}\left\\{\mathbf{\Upsilon}_{k}^{2}\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right\\},\quad
f_{c,k,4}(\bm{\Phi})=\operatorname{Tr}\left\\{\bm{\Upsilon}_{k}^{2}\right\\},\\\
&f_{c,k,5}(\bm{\Phi})=\left|\operatorname{Tr}\left\\{\bm{\Upsilon}_{k}\right\\}\right|^{2},\quad
f_{c,k,6}(\bm{\Phi})=\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k},\\\
&f_{c,k,7}(\bm{\Phi})=\left|f_{k}(\bm{\Phi})\right|^{2},\quad
f_{c,ki,8}(\bm{\Phi})=\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\mathbf{\Upsilon}_{k}^{2}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\overline{\mathbf{h}}_{i},\\\
&f_{c,ki,9}(\bm{\Phi})=\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}.\end{aligned}$
(95)

Using notations in (95), we next present the closed-form expression of the
achiveable rate.

######  Theorem 4.

With spatial correlation and EMI, the UatF bound for user $k$’s achievable
rate is given by

$\displaystyle\underline{R}_{c,k}=\tau^{o}\log_{2}\left(1+\operatorname{SINR}_{c,k}\right),$
(96)
$\displaystyle\operatorname{SINR}_{c,k}=\frac{pE_{c,k}^{\mathrm{signal}}}{pE_{c,k}^{\mathrm{leak}}+p\sum_{i=1,i\neq
k}^{K}I_{c,ki}+\sigma_{e}^{2}E_{c,k}^{\mathrm{emi}}+\sigma^{2}E_{c,k}^{\mathrm{noise}}},$
(97)

where the signal term is
$E_{c,k}^{\mathrm{signal}}=\left(E_{c,k}^{\mathrm{noise}}\right)^{2}$ and the
noise term is

$\displaystyle
E_{c,k}^{\mathrm{noise}}=M\widehat{c}_{k}\delta\left|f_{k}(\bm{\Phi})\right|^{2}+\widehat{c}_{k}\operatorname{Tr}\left\\{\mathbf{\Upsilon}_{k}\right\\}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}+\gamma_{k}\operatorname{Tr}\left\\{\mathbf{\Upsilon}_{k}\right\\}.$
(98)

The EMI term is denoted by
$E_{c,k}^{\mathrm{emi}}=\frac{\beta}{\delta+1}\sum_{\omega=1}^{8}E_{c,k}^{\omega,\mathrm{emi}}$
where

$\displaystyle\begin{aligned}
E_{c,k}^{1,\mathrm{emi}}&=M^{2}\widehat{c}_{k}\delta^{2}f_{c,k,7}(\bm{\Phi})\mathbf{a}_{N}^{H}\bm{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\mathbf{a}_{N},\\\
E_{c,k}^{2,\mathrm{emi}}&=\left(\widehat{c}_{k}\delta
f_{c,k,2}(\bm{\Phi})+\frac{2\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,1}(\bm{\Phi})+\delta\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)\right)f_{c,k,3}(\bm{\Phi}),\\\
E_{c,k}^{3,\mathrm{emi}}&=\left(M\widehat{c}_{k}\delta
f_{c,k,7}(\bm{\Phi})+\left(\frac{\sigma^{2}}{\tau
p}+\gamma_{k}+\widehat{c}_{k}f_{c,k,2}(\bm{\Phi})+\frac{\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,1}(\bm{\Phi})\right)f_{c,k,4}(\bm{\Phi})\right)f_{c,1}(\bm{\Phi}),\\\
E_{c,k}^{4,\mathrm{emi}}&=\frac{\beta\delta^{2}\sigma_{e}^{2}}{\tau
p(\delta+1)}\operatorname{Tr}\left\\{\left(\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}{\bf\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\right)^{2}\right\\},\\\
E_{c,k}^{5,\mathrm{emi}}&=2\widehat{c}_{k}\delta\operatorname{Tr}\left\\{{\bf\Upsilon}_{k}\right\\}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}_{emi}\mathbf{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}\right\\},\\\
E_{c,k}^{6,\mathrm{emi}}&=\frac{2\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}\operatorname{Tr}\left\\{{\bf\Upsilon}_{k}\right\\}\operatorname{Tr}\left\\{\mathbf{R}_{emi}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}{\bf\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\right\\},\\\
E_{c,k}^{7,\mathrm{emi}}&=\widehat{c}_{k}f_{c,k,5}(\bm{\Phi})f_{c,k,6}(\bm{\Phi}),\\\
E_{c,k}^{8,\mathrm{emi}}&=\frac{\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,k,5}(\bm{\Phi})\operatorname{Tr}\left\\{\left(\mathbf{R}_{ris}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\right)^{2}\right\\}.\end{aligned}$
(99)

The interference term is $I_{c,ki}=\sum_{\omega=1}^{8}I_{c,ki}^{\omega}$ where

$\displaystyle\begin{aligned}
&I_{c,ki}^{1}=\gamma_{i}E_{c,k}^{\mathrm{noise}}+M^{2}\widehat{c}_{k}\widehat{c}_{i}\delta^{2}f_{c,k,7}(\mathbf{\Phi})f_{c,i,7}(\mathbf{\Phi}),\\\
&I_{c,ki}^{2}=\bigg{\\{}M\widehat{c}_{k}\widehat{c}_{i}\delta
f_{c,k,7}(\bm{\Phi})+\left(\widehat{c}_{i}\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)+\frac{\widehat{c}_{i}\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,1}(\bm{\Phi})\right)f_{c,k,4}(\bm{\Phi})\\\
&\qquad\quad+\frac{\widehat{c}_{i}\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,k,3}(\bm{\Phi})\bigg{\\}}f_{c,i,2}(\bm{\Phi}),\\\
&I_{c,ki}^{3}=\left\\{\widehat{c}_{k}\widehat{c}_{i}\delta
f_{c,ki,8}(\bm{\Phi})+\widehat{c}_{k}\widehat{c}_{i}f_{c,k,4}(\bm{\Phi})f_{c,i,2}(\mathbf{\Phi})\right\\}f_{c,k,2}(\bm{\Phi}),\\\
&I_{c,ki}^{4}=\left\\{\frac{\widehat{c}_{i}\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,1}(\bm{\Phi})+\widehat{c}_{i}\delta\left(\gamma_{k}+\frac{\sigma^{2}}{\tau
p}\right)\right\\}f_{c,ki,8}(\bm{\Phi}),\\\
&I_{c,ki}^{5}=\left\\{\widehat{c}_{k}\widehat{c}_{i}\left|\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{i}\right|^{2}+\frac{\widehat{c}_{i}\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,i,6}(\bm{\Phi})\right\\}f_{c,k,5}(\mathbf{\Phi}),\\\
&I_{c,ki}^{6}=2\widehat{c}_{k}\widehat{c}_{i}\delta\operatorname{Tr}\left\\{{\bf\Upsilon}_{k}\right\\}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}\right\\},\\\
&I_{c,ki}^{7}=\frac{\widehat{c}_{i}\beta\delta^{2}\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,ki,9}(\bm{\Phi}),\\\
&I_{c,ki}^{8}=\frac{2\widehat{c}_{i}\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}\operatorname{Tr}\left\\{{\bf\Upsilon}_{k}\right\\}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\mathbf{R}_{ris}\mathbf{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}{\bf\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\right\\}.\end{aligned}$
(100)

The signal leakage term is
$E_{c,k}^{\mathrm{leak}}=\sum_{\omega=1}^{8}E_{c,k}^{\omega,\mathrm{leak}}$,
where

$\displaystyle\begin{aligned}
E_{c,k}^{1,\mathrm{leak}}&=M\widehat{c}_{k}\delta\gamma_{k}f_{c,k,7}(\bm{\Phi}),\\\
E_{c,k}^{2,\mathrm{leak}}&=\left\\{M\widehat{c}_{k}^{2}\delta
f_{c,k,7}(\bm{\Phi})+\widehat{c}_{k}^{2}\delta
f_{c,kk,8}(\bm{\Phi})+\left(\widehat{c}_{k}^{2}f_{c,k,2}(\bm{\Phi})+2\widehat{c}_{k}\gamma_{k}+\frac{\widehat{c}_{k}\sigma^{2}}{\tau
p}\right)f_{c,k,4}(\bm{\Phi})\right\\}f_{c,k,2}(\bm{\Phi}),\\\
E_{c,k}^{3,\mathrm{leak}}&=\left\\{\widehat{c}_{k}\delta\gamma_{k}+\frac{\widehat{c}_{k}\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,1}(\bm{\Phi})+\frac{\widehat{c}_{k}\delta\sigma^{2}}{\tau
p}\right\\}f_{c,kk,8}(\bm{\Phi}),\\\
E_{c,k}^{4,\mathrm{leak}}&=\left\\{\gamma_{k}^{2}+\frac{\gamma_{k}\sigma^{2}}{\tau
p}+\frac{\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}\left(\gamma_{k}+\widehat{c}_{k}f_{c,k,2}(\bm{\Phi})\right)f_{c,1}(\bm{\Phi})\right\\}f_{c,k,4}(\bm{\Phi}),\\\
E_{c,k}^{5,\mathrm{leak}}&=\frac{\widehat{c}_{k}\beta\delta^{2}\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,kk,9}(\bm{\Phi}),\\\
E_{c,k}^{6,\mathrm{leak}}&=\frac{2\widehat{c}_{k}\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}\operatorname{Tr}\left\\{\mathbf{\Upsilon}_{k}^{H}\right\\}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}{\bf\Upsilon}_{k}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\mathbf{R}_{ris}\mathbf{\Phi}\overline{\mathbf{h}}_{k}\right\\},\\\
E_{c,k}^{7,\mathrm{leak}}&=\frac{\beta\delta\sigma_{e}^{2}}{\tau
p(\delta+1)}\left\\{\gamma_{k}+\widehat{c}_{k}f_{c,k,2}(\bm{\Phi})\right\\}f_{c,k,3}(\bm{\Phi}),\\\
E_{c,k}^{8,\mathrm{leak}}&=\frac{\widehat{c}_{k}\beta\sigma_{e}^{2}}{\tau
p(\delta+1)}f_{c,k,5}(\bm{\Phi})f_{c,k,6}(\bm{\Phi}).\end{aligned}$ (101)

Proof: See Appendix H. $\blacksquare$

Compared the rate expression $\underline{R}_{c,k}$ in Theorem 4 with
$\underline{R}_{k}$ in Theorem 2, we are able to find some insights with
respect to spatial correlation and EMI. The impact of spatial correlation on
the achievable rate is discussed in the following remark.

######  Remark 3.

After considering the spatial correlation at the RIS, the RIS would possess a
greater ability to tailor the wireless channels and improve the achievable
rate. This is because compared with rate expressions with spatial-independent
channels, the rate with spatial correlation possesses more terms involving the
phase shifts variable $\bf\Phi$. These terms can be exploited in the
optimization and therefore the RIS’s capability of modifying the propagation
environment is significantly enhanced. To be specific, take the term
$\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}$
for example. In the spatial-independent case, this term is fixed (equal to
$N$) without any possibility to be adjusted. However, it can be smartly
adjusted in the presence of spatial correlation. For the sake of exposition,
let us assume the most severe spatial-correlated case where
$\mathbf{R}_{ris}=\mathbf{1}_{N\times N}$ so that
$\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}=\left|\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{1}_{N\times
1}\right|^{2}$. Following the proof in Lemma 3, we have
$0\leq\left|\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{1}_{N\times
1}\right|^{2}\leq N^{2}$, which demonstrates the stronger and more flexible
adjustment ability of the RIS in the presence of spatial correlation.

Next, we would like to discuss the impact of EMI on the power scaling laws.
Due to the complex expressions in (99) and the undecided variable $\bf\Phi$,
the general conclusion cannot be given. However, by investigating some special
cases, we can present the following corollary based on the proof by
contradiction.

######  Corollary 11.

Power scaling laws shown in Table II are not guaranteed to hold in the
presence of EMI.

Proof: On one hand, it can be seen that both desired signal
$E_{c,k}^{\mathrm{signal}}$ and one of the EMI entries
$E_{c,k}^{1,\mathrm{emi}}$ in (99) scale on the order of
$\mathcal{O}\left(M^{2}\right)$. This means that if we reduce the power
proportionally to $p=1/M$, the SINR in (97) will tend to zero when
$M\to\infty$. On the other hand, to facilitate the examination of power
scaling laws with $N$, let us focus on the fully NLoS spatial-independent case
where $\delta=0$ and $\mathbf{R}_{ris}=\mathbf{R}_{emi}=\mathbf{I}_{N}$. In
this context, we can simplify $\mathbf{\Upsilon}_{k}$ as

$\displaystyle\operatorname{Tr}\\{\mathbf{\Upsilon}_{k}\\}=\frac{M(N\widehat{c}_{k}+\gamma_{k})}{N\widehat{c}_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}+\frac{N\sigma_{e}^{2}\beta}{\tau p}}.$ (102)

Then, we have
$E_{c,k}^{\mathrm{signal}}=\left(E_{c,k}^{\mathrm{noise}}\right)^{2}$ where

$\displaystyle
E_{c,k}^{\mathrm{noise}}=\frac{M(N\widehat{c}_{k}+\gamma_{k})^{2}}{N\widehat{c}_{k}+\gamma_{k}+\frac{\sigma^{2}}{\tau
p}+\frac{N\sigma_{e}^{2}\beta}{\tau p}}.$ (103)

From (103), it can be seen that if the power is scaled proportionally to
$p=1/N$, when $N\to\infty$, we have $E_{c,k}^{\mathrm{signal}}\to(\frac{\tau
M\widehat{c}^{2}_{k}}{\sigma_{e}^{2}\beta})^{2}$ which means that the
$pE_{c,k}^{\mathrm{signal}}=E_{c,k}^{\mathrm{signal}}/N\to 0$ and then the
SINR would tend to zero. This special case demonstrates that the power scaling
laws with respect to $N$ are not guaranteed to hold after considering EMI.
$\blacksquare$

The reason behind Corollary 11 is that when the transmit power $p$ is reduced
proportionally to $1/M$ or $1/N$, as $M$ or $N$ increases, the useful signal
power received by the RIS becomes weaker and weaker while the power of EMI
received by the RIS is unchanged. Thus, taking the strength of the useful
signal as a reference, the EMI becomes stronger and stronger. Finally, as
$M,N\to\infty$, the useful power is extremely weak and then the EMI power
dominates the received signal at the RIS. Note that even though the useful
signal is boosted with larger $M$ or $N$, the EMI is boosted as well. This
leads to a dominant EMI term in the received SINR and therefore beaks the
power scaling laws.

Nevertheless, the importance of power scaling laws does not lie in the
performance limits in the asymptotic region of $M,N\to\infty$. In practice,
neither the number of BS antennas nor the number of RIS elements can be
infinite. The meaning of investigating power scaling laws is that it answers
the question of whether we can help users save their transmit power by
increasing $M$ or $N$ while not sacrificing the rate performance too much.
Therefore, we are more interested in the performance of power scaling laws
when $M$ or $N$ is large but finite.

When $M$ or $N$ is large but finite, we can expect that the power scaling laws
still hold if the EMI is not severe. Specifically, the users in RIS-aided
systems share the same RIS-BS channel, and therefore there exists strong
multi-user interference when using MRC. Note that the existence of strong
interference has been revealed in Remark 2. Even though the EMI boosted by the
RIS could be much stronger than the thermal noise, it could still be much
smaller than the multi-user interference which is boosted by the RIS as well.
When $M$ or $N$ is not very large, the transmit power, which has been reduced
proportionally to $1/M$ or $1/N$, is small (not extremely small), but it could
still result in strong multi-user interference and this interference makes the
impact of EMI negligible. Therefore, in this region, the power scaling laws
could still hold. The exact impacts of EMI on the achievable rate and power
scaling laws will be examined in the simulations.

##  VI Design of the RIS Phase Shifts

In this section, we optimize the phase shifts of the RIS to maximize the
achievable rate derived in Theorem 2 and Theorem 4. Since the derived ergodic
rate depends only on statistical CSI, we do not need to update the phase
shifts design of the RIS until the long-term CSI varies, which results in a
much lower update frequency compared with conventional designs based on
instantaneous CSI in the sub-6 GHz bands. Therefore, the proposed optimization
criterion has a reduced computational complexity.

###  VI-A Single-user Case

Before tackling the general optimization problem, we would like to first
justify the setting of $\left|f_{k}(\bm{\Phi})\right|=N$ in Section IV-C. To
this end, this subsection aims to solve the phase shift optimization problem
in the special single-user case without spatial correlation and EMI.

In the single-user case, only the user $k$ is present. We need to find the
phase shifts matrix $\bf\Phi$ that maximizes the lower bound of the ergodic
rate $\underline{R}_{k}$, which is obtained by setting
$I_{ki}\left(\mathbf{\Phi}\right)=0$ in Theorem 2. Note that only the
scenarios with $N>1$, $\delta>0$ and $\varepsilon_{k}>0,\forall k$ are
considered (if not, $\bf\Phi$ can be set arbitrarily). We can then find that
the rate depends on the phase shifts matrix $\bf\Phi$ only through the term
$\left|f_{k}\left(\mathbf{\Phi}\right)\right|^{2}$. By defining
$x\triangleq\left|f_{k}\left(\mathbf{\Phi}\right)\right|^{2}$, we can rewrite
the rate $\underline{R}_{k}$ in Theorem 2 as follows

$\displaystyle\underline{R}_{k}$
$\displaystyle=\tau^{o}\log_{2}\left(1+\mathrm{SNR}_{k}\left(x\right)\right)$
$\displaystyle=\tau^{o}\log_{2}\left(1+\frac{E_{k}^{\rm
signal}\left(x\right)}{E_{k}^{\rm
leak}\left(x\right)+\frac{\sigma^{2}}{p}E_{k}^{\rm{noise}}\left(x\right)}\right)$
$\displaystyle\triangleq\tau^{o}\log_{2}\left(1+\frac{\left(s_{1}x+s_{2}\right)^{2}}{t_{1}x+t_{2}}\right),$
(104)

where
$\mathrm{SNR}_{k}\left(x\right)\triangleq\frac{\left(s_{1}x+s_{2}\right)^{2}}{t_{1}x+t_{2}}$,
$E_{k}^{\text{signal}}\left(x\right)\triangleq\left(s_{1}x+s_{2}\right)^{2}$
and
$E_{k}^{\text{leak}}\left(x\right)+\frac{\sigma^{2}}{p}E_{k}^{\text{noise}}\left(x\right)\triangleq
t_{1}x+t_{2}$. The exact expressions of $s_{1}$, $s_{2}$, $t_{1}$, and $t_{2}$
are omitted due to their cumbersome formulation and because they are
irrelevant for solving the optimization problem of interest. They can be
retrieved by direct inspection of Theorem 2. The only important property to
highlight is that $s_{1},s_{2},t_{1},t_{2}>0$.

Note that we have assumed that $N>1$ in this subsection. From Lemma 3, we know
that the domain of the variable $x$ is $0\leq x\leq N^{2}$. Based on (VI-A),
therefore, the optimization problem can be formulated as follows

$\displaystyle\max_{x}\;\;\mathrm{SNR}_{k}\left({x}\right)=\frac{\left(s_{1}x+s_{2}\right)^{2}}{t_{1}x+t_{2}},$
(105a) $\displaystyle\text{ s.t. }\quad 0\leq x\leq N^{2}.$ (105b)

To solve the problem in (105), we compute the first-order derivative of
$\mathrm{SNR}_{k}\left(x\right)$ with respect to $x$, as follows

$\displaystyle\frac{\partial\mathrm{SNR}_{k}\left(x\right)}{\partial
x}=\frac{\left(s_{1}x+s_{2}\right)\left(s_{1}t_{1}x+2s_{1}t_{2}-s_{2}t_{1}\right)}{\left(t_{1}x+t_{2}\right)^{2}}.$
(106)

The first-order derivative of $\mathrm{SNR}_{k}\left(x\right)$ is positive or
negative depending on the numerator in (106), which is a quadratic function of
$x$, i.e., a parabola opening upward, with two roots. The two roots can be
obtained by setting (106) equal to zero, which yields

$\displaystyle x_{0}^{L}=\frac{-s_{2}}{s_{1}},\quad
x_{0}^{R}=\frac{s_{2}t_{1}-2s_{1}t_{2}}{s_{1}t_{1}},$ (107)

where $x_{0}^{L}<0$ while $x_{0}^{R}$ can be positive.

We can design the optimal configuration of $\bf\Phi$ by analyzing the
derivative $\frac{\partial\mathrm{SNR}_{k}\left(x\right)}{\partial x}$ in the
domain of $x$, i.e., (105b), which can be identified by discussing the value
of $x_{0}^{R}$. For example, if $x_{0}^{R}\leq 0$, for a parabola opening
upward, we obtain $\frac{\partial\mathrm{SNR}_{k}\left(x\right)}{\partial
x}\geq 0$ in the domain $0\leq x\leq N^{2}$. The design framework is
summarized in the following theorem.

######  Theorem 5.

For RIS-aided single-user systems subject to imperfect CSI, the optimal phase
shift matrix $\bf\Phi$ in terms of the UatF bound of the achievable rate is as
follows

  * •

Case 1: If $x_{0}^{R}\leq 0$, the optimal design for $\bf\Phi$ corresponds to
setting $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$.

  * •

Case 2: If $x_{0}^{R}\geq N^{2}$, the optimal design for $\bf\Phi$ corresponds
to setting $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=0$.

  * •

Case 3: If $0<x_{0}^{R}<N^{2}$, the optimal design for $\bf\Phi$ corresponds
to setting $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=0$ when ${\rm
SNR}_{k}\left(0\right)>{\rm SNR}_{k}\left(N^{2}\right)$ and
$\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$ when ${\rm
SNR}_{k}\left(0\right)\leq{\rm SNR}_{k}\left(N^{2}\right)$.

  * •

Case 4: If $N\to\infty$, the optimal design for $\bf\Phi$ corresponds to
setting $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$.

Proof: In Case 1, we obtain
$\frac{\partial\mathrm{SNR}_{k}\left(x\right)}{\partial x}\geq 0$ in the
domain $0\leq x\leq N^{2}$. Thus, the SNR is an increasing function of $x$ in
its domain, which implies that the maximum SNR is reached at the endpoint
$x=N^{2}$. Therefore, it is optimal to set
$\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$. In Case 2, we obtain
$\frac{\partial\mathrm{SNR}_{k}\left(x\right)}{\partial x}\leq 0$ in the
domain of $x$. Thus, the SNR is a decreasing function of $x$, which implies
that the maximum SNR is reached at the endpoint $x=0$. Therefore, it is
optimal to set $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=0$. In Case 3,
the SNR first decreases for $x<x_{0}^{R}$, and then increases for
$x>x_{0}^{R}$. Therefore, the maximum SNR is obtained either at $x=0$ or at
$x=N^{2}$. By comparing ${\rm SNR}_{k}\left(0\right)$ with ${\rm
SNR}_{k}\left(N^{2}\right)$, we can identify the optimal design. In Case 4,
from Theorem 2, we have the following: (1) if
$\left|f_{k}\left(\mathbf{\Phi}\right)\right|=0$, the SNR is bounded if
$N\to\infty$; and (2) if $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$, the
SNR increases without bound if $N\to\infty$. Therefore, the optimum SNR is
attained for $\left|f_{k}\left(\mathbf{\Phi}\right)\right|=N$ if $N\to\infty$.
$\blacksquare$

Finally, we note that the optimal design obtained in Case 4 substantiates the
analysis reported in Section IV-C for large $N$.

###  VI-B Multi-user Case

In this subsection, we focus on the design of the RIS phase shifts in the
general multi-user scenario with $K>1$. In the multi-user case, as mentioned
in Remark 2, it is necessary to guarantee some fairness among different users.
To this end, we aim to maximize the minimum rate of the users. As a result,
the optimization problem can be formulated as follows

$\displaystyle\max_{\bf\Phi}\;\min_{k\in\mathcal{K}}\;\;\underline{R}_{k}\left({\bf\Phi}\right)\text{
or }\underline{R}_{c,k}\left({\bf\Phi}\right),$ (108a) $\displaystyle\text{
s.t. }\quad\left|\left[\mathbf{\Phi}\right]_{n,n}\right|=1,\forall n,$ (108b)

where $\underline{R}_{k}\left({\bf\Phi}\right)$ is given by (39) in Theorem 2
and $\underline{R}_{c,k}\left({\bf\Phi}\right)$ is given by (96) in Theorem 4.
Constraint (108b) is the unit modulus constraint for the RIS phase shift
matrix.

For tractability, define vector
$\bm{\theta}=[\theta_{1},\theta_{2},\ldots,\theta_{N}]^{T}$ and
$\bm{c}=[e^{j\theta_{1}},e^{j\theta_{2}},\ldots,e^{j\theta_{N}}]^{T}$ so that
$\bm{c}=e^{j\bm{\theta}}$ and
$\mathbf{\Phi}=\mathrm{diag}\left(\bm{c}\right)$. Then, the problem in (108)
can be solved effectively based on the gradient ascent method with respect to
real variable $\bm{\theta}$. It is worth noting that our method is different
from the existing work which adopted the projected gradient ascent method with
respect to complex variable $\bm{c}$[46]. To be specific, after updating
$\bm{c}$, a projection operation is needed to enforce the updated variable
$\bm{c}_{new}$ satisfying the unit modulus constraint
$\left|\bm{c}_{new}\right|=\mathbf{1}$. By contrast, our method can avoid the
rate loss in the projection operation since the sine and cosine functions are
periodic with $\bm{\theta}$ and the unit modulus constraint always holds for
all angles. Besides, the gradient ascent method is sensitive to the setting of
the step size, and the method of adjusting step sizes is more robust in the
real field.

Next, we need to calculate the gradient with respect to $\bm{\theta}$.
However, the objective function in (108) includes a minimal function, which is
non-differential. To tackle this challenge, we first approximate the objective
function in (108) as

$\displaystyle\min_{k}\underline{R}_{k}(\bm{\theta})\approx-\frac{1}{\mu}\ln\left\\{\sum_{k=1}^{K}\exp\left\\{-\mu\underline{R}_{k}(\mathbf{\bm{\theta}})\right\\}\right\\}\triangleq
f(\mathbf{\bm{\theta}}),$ (109)
$\displaystyle\min_{k}\underline{R}_{c,k}(\bm{\theta})\approx-\frac{1}{\mu}\ln\left\\{\sum_{k=1}^{K}\exp\left\\{-\mu\underline{R}_{c,k}(\mathbf{\bm{\theta}})\right\\}\right\\}\triangleq
f_{c}(\mathbf{\bm{\theta}}),$ (110)

where $\mu$ is a constant value for controlling the approximation accuracy. It
can be proved that the approximation error is smaller than $\frac{\ln K}{\mu}$
following a similar method in [72]. Thus, the problem in (108) can be recast
as

$\displaystyle\max_{\bm{\theta}}\;f(\mathbf{\bm{\theta}})\text{ or
}f_{c}(\mathbf{\bm{\theta}}),$ (111a) $\displaystyle\text{ s.t. }\quad
0\leq\theta_{n}<2\pi,\forall n.$ (111b)

Note that constraint (111b) can be neglected since the objective functions
$f(\mathbf{\bm{\theta}})$ and $f_{c}(\mathbf{\bm{\theta}})$ including sine and
cosine functions are periodic with respect to $\bm{\theta}$. Therefore, there
is no need to perform the projection operation after updating variable
$\bm{\theta}$. Then, we need to calculate the gradient of
$f(\mathbf{\bm{\theta}})$ and $f_{c}(\mathbf{\bm{\theta}})$. Since these two
gradients can be calculated in a similar way, we only provide the detailed
process of deriving $\frac{\partial f_{c}(\bm{\theta})}{\partial\bm{\theta}}$
for brevity. Based on the chain rule, we have

$\displaystyle\frac{\partial
f_{c}(\bm{\theta})}{\partial\bm{\theta}}=\frac{\tau^{o}\sum_{k=1}^{K}\left\\{\frac{\exp\left\\{-\mu\underline{R}_{c,k}(\bm{\theta})\right\\}}{1+\operatorname{SINR}_{c,k}(\bm{\theta})}\frac{\partial\operatorname{SINR}_{c,k}(\bm{\theta})}{\partial\bm{\theta}}\right\\}}{(\ln
2)\left(\sum_{k=1}^{K}\exp\left\\{-\mu\underline{R}_{c,k}(\bm{\theta})\right\\}\right)},$
(112)

and

$\displaystyle\begin{aligned}
&\frac{\partial\operatorname{SINR}_{c,k}(\bm{\theta})}{\partial\bm{\theta}}=\frac{p\frac{\partial
E_{c,k}^{\mathrm{signal}}}{\partial\bm{\theta}}}{pE_{c,k}^{\mathrm{leak}}+p\sum_{i=1,i\neq
k}^{K}I_{c,ki}+\sigma_{e}^{2}E_{c,k}^{\mathrm{emi}}+\sigma^{2}E_{c,k}^{\mathrm{noise}}}\\\
&-pE_{c,k}^{\mathrm{signal}}\frac{p\frac{\partial
E_{c,k}^{\mathrm{leak}}}{\partial\bm{\theta}}+p\sum_{i=1,i\neq
k}^{K}\frac{\partial
I_{c,ki}}{\partial\bm{\theta}}+\sigma_{e}^{2}\frac{\partial
E_{c,k}^{\mathrm{emi}}}{\partial\bm{\theta}}+\sigma^{2}\frac{\partial
E_{c,k}^{\mathrm{noise}}}{\partial\bm{\theta}}}{\left(pE_{c,k}^{\mathrm{leak}}+p\sum_{i=1,i\neq
k}^{K}I_{c,ki}+\sigma_{e}^{2}E_{c,k}^{\mathrm{emi}}+\sigma^{2}E_{c,k}^{\mathrm{noise}}\right)^{2}}.\end{aligned}$
(113)

Therefore, the gradient of $f_{c}(\bm{\theta})$ can be obtained after
calculating $\frac{\partial E_{c,k}^{\mathrm{signal}}}{\partial\bm{\theta}}$,
$\frac{\partial E_{c,k}^{\mathrm{leak}}}{\partial\bm{\theta}}$,
$\frac{\partial I_{c,ki}}{\partial\bm{\theta}}$, $\frac{\partial
E_{c,k}^{\mathrm{emi}}}{\partial\bm{\theta}}$, and $\frac{\partial
E_{c,k}^{\mathrm{noise}}}{\partial\bm{\theta}}$ in (113). Recall Theorem 4. It
can be seen that $E_{c,k}^{\mathrm{signal}}$, $I_{c,ki}$,
$E_{c,k}^{\mathrm{leak}}$, $E_{c,k}^{\mathrm{emi}}$ and
$E_{c,k}^{\mathrm{noise}}$ are nearly comprised of terms listed in (95).
Therefore, we first provide two useful lemmas and then use them to calculate
the gradient of terms in (95).

######  Lemma 4.

For deterministic matrices $\mathbf{A}$ and $\mathbf{B}$, the gradient of
$\operatorname{Tr}\left\\{\mathbf{A}\bm{\Phi}\mathbf{B}\bm{\Phi}^{H}\right\\}$
with respect to $\bm{\theta}$ is given by

$\displaystyle\begin{aligned}
\frac{\partial\operatorname{Tr}\left\\{\mathbf{A}\bm{\Phi}\mathbf{B}\bm{\Phi}^{H}\right\\}}{\partial\bm{\theta}}&=j\bm{\Phi}^{T}\left(\mathbf{A}^{T}\odot\mathbf{B}\right)\bm{c}^{*}-j\bm{\Phi}^{H}\left(\mathbf{A}\odot\mathbf{B}^{T}\right)\bm{c}\\\
&\triangleq\bm{f}_{d}(\mathbf{A},\mathbf{B}).\end{aligned}$ (114)

If $\mathbf{A}=\mathbf{A}^{H},\mathbf{B}=\mathbf{B}^{H}$, we further have

$\displaystyle\frac{\partial\operatorname{Tr}\left\\{\mathbf{A}\bm{\Phi}\mathbf{B}\bm{\Phi}^{H}\right\\}}{\partial\bm{\theta}}=2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{A}\odot\mathbf{B}^{T}\right)\bm{c}\right\\}.$
(115)

Proof: See Appendix I. $\blacksquare$

######  Lemma 5.

Define
$\psi_{k}^{1}=\widehat{c}_{k}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\mathbf{\Phi}\overline{\mathbf{h}}_{k}+\gamma_{k}$
and $\mathbf{\Upsilon}_{k}=\psi_{k}^{1}\mathbf{\Upsilon}_{k}^{1}$. Then, for
deterministic matrix $\bf T$, the gradient of
$\operatorname{Tr}\left\\{\mathbf{T}\mathbf{\Upsilon}_{k}\right\\}$ with
respect to $\bm{\theta}$ is given by

$\displaystyle\begin{aligned}
&\frac{\partial\operatorname{Tr}\left\\{\mathbf{T}\mathbf{\Upsilon}_{k}\right\\}}{\partial\bm{\theta}}\\\
&=2\widehat{c}_{k}\left\\{\operatorname{Tr}\left\\{\mathbf{T}\mathbf{\Upsilon}_{k}^{1}\right\\}-\psi_{k}^{1}\operatorname{Tr}\left\\{\mathbf{T}\left(\mathbf{\Upsilon}_{k}^{1}\right)^{2}\right\\}\right\\}\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{R}_{ris}\odot\left(\overline{\mathbf{h}}_{k}\overline{\mathbf{h}}_{k}^{H}\right)^{T}\right)\bm{c}\right\\}\\\
&-\frac{2\sigma_{e}^{2}\beta}{\tau
p(\delta+1)}\psi_{k}^{1}\operatorname{Tr}\left\\{\mathbf{T}\left(\mathbf{\Upsilon}_{k}^{1}\right)^{2}\right\\}\operatorname{Im}\left\\{\mathbf{\Phi}^{H}\left(\mathbf{R}_{ris}\odot\mathbf{R}_{emi}\right)\bm{c}\right\\}\\\
&-\frac{\sigma_{e}^{2}\beta\delta}{\tau
p(\delta+1)}\psi_{k}^{1}\bm{f}_{d}\left(\overline{\mathbf{H}}_{2}^{H}\mathbf{\Upsilon}_{k}^{1}\mathbf{T}\mathbf{\Upsilon}_{k}^{1}\overline{\mathbf{H}}_{2},\mathbf{R}_{emi}\right)\\\
&\triangleq\mathbf{z}_{k}(\mathbf{T})\end{aligned}$ (116)

Proof: It follows a similar procedure to the proof of Lemma 4 after using the
chain rule of inverse matrix
$\partial\left(\mathbf{X}^{-1}\right)=-\mathbf{X}^{-1}(\partial\mathbf{X})\mathbf{X}^{-1}$.
$\blacksquare$

Next, using the chain rule and applying Lemma 4 and 5, we are able to obtain
the following lemma.

######  Lemma 6.

Gradients of terms in (95) with respect to $\bm{\theta}$ are given by

$\displaystyle\begin{aligned}
\bm{f}_{c,1}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,1}(\bm{\Phi})}{\partial\bm{\theta}}=2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{R}_{ris}\odot\mathbf{R}_{emi}\right)\bm{c}\right\\},\\\
\bm{f}_{c,k,2}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,2}(\bm{\Phi})}{\partial\bm{\theta}}=2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{R}_{ris}\odot\left(\overline{\mathbf{h}}_{k}\overline{\mathbf{h}}_{k}^{H}\right)^{T}\right)\bm{c}\right\\},\\\
\bm{f}_{c,k,3}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,3}(\bm{\Phi})}{\partial\bm{\theta}}\\\
&=\mathbf{z}_{k}\left(\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\right)+2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{2}\overline{\mathbf{H}}_{2}\odot\mathbf{R}_{emi}\right)\bm{c}\right\\}\\\
&+\mathbf{z}_{k}\left(\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right),\end{aligned}$
(117) $\displaystyle\begin{aligned}
\bm{f}_{c,k,4}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,4}(\bm{\Phi})}{\partial\bm{\theta}}=2\mathbf{z}_{k}\left(\bm{\Upsilon}_{k}\right),\\\
\bm{f}_{c,k,5}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,5}(\bm{\Phi})}{\partial\bm{\theta}}=2\operatorname{Tr}\left\\{\mathbf{\Upsilon}_{k}\right\\}\mathbf{z}_{k}\left(\mathbf{I}_{M}\right),\\\
\bm{f}_{c,k,6}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,6}(\bm{\Phi})}{\partial\bm{\theta}}\\\
&=2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{R}_{ris}\bm{\Phi}\overline{\mathbf{h}}_{k}\overline{\mathbf{h}}_{k}^{H}\bm{\Phi}^{H}\mathbf{R}_{ris}\odot\mathbf{R}_{emi}\right)\bm{c}\right\\}\\\
&+2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{R}_{ris}\mathbf{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\mathbf{R}_{ris}\odot\left(\overline{\mathbf{h}}_{k}\overline{\mathbf{h}}_{k}^{H}\right)^{T}\right)\bm{c}\right\\},\quad\qquad\qquad\end{aligned}$
(118) $\displaystyle\begin{aligned}
\bm{f}_{c,k,7}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,k,7}(\bm{\Phi})}{\partial\bm{\theta}}=2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\mathbf{a}_{N}\mathbf{a}_{N}^{H}\odot\left(\overline{\mathbf{h}}_{k}\overline{\mathbf{h}}_{k}^{H}\right)^{T}\right)\bm{c}\right\\},\\\
\bm{f}_{c,ki,8}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,ki,8}(\bm{\Phi})}{\partial\bm{\theta}}\\\
&=\mathbf{z}_{k}\left(\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right)+\mathbf{z}_{k}\left(\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\right)\\\
&+2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{2}\overline{\mathbf{H}}_{2}\odot\left(\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\right)^{T}\right)\bm{c}\right\\},\quad\qquad\qquad\qquad\qquad\end{aligned}$
(119) $\displaystyle\begin{aligned}
\bm{f}_{c,ki,9}^{\prime}(\bm{\theta})&=\frac{\partial
f_{c,ki,9}(\bm{\Phi})}{\partial\bm{\theta}}\\\
&=\mathbf{z}_{k}\left(\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right)\\\
&+2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\odot\mathbf{R}_{emi}\right)\bm{c}\right\\}\\\
&+\mathbf{z}_{k}\left(\overline{\mathbf{H}}_{2}\bm{\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right)\\\
&+2\operatorname{Im}\left\\{\bm{\Phi}^{H}\left(\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\mathbf{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{H}\overline{\mathbf{H}}_{2}\odot\left(\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\right)^{T}\right)\bm{c}\right\\}.\end{aligned}$
(120)

Proof: It follows by applying the chain rules and using Lemma 4 and 5. Take
$\bm{f}_{c,ki,3}^{\prime}(\bm{\theta})$ for example. Apply the chain rule, we
have

$\displaystyle\begin{aligned}
&\bm{f}_{c,k,3}^{\prime}(\bm{\theta})=\frac{\partial\operatorname{Tr}\left\\{\bm{\Upsilon}_{k}\bm{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right\\}}{\partial\bm{\theta}}=\left.\frac{\partial\operatorname{Tr}\left\\{\mathbf{T}\bm{\Upsilon}_{k}\right\\}}{\partial\bm{\theta}}\right|_{\mathbf{T}=\mathbf{\Upsilon}_{k}\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}}\\\
&+\left.\frac{\partial\operatorname{Tr}\left\\{\mathbf{T}\bm{\Upsilon}_{k}\right\\}}{\partial\bm{\theta}}\right|_{\mathbf{T}=\overline{\mathbf{H}}_{2}\bm{\Phi}\mathbf{R}_{emi}\bm{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}}+\left.\frac{\partial\operatorname{Tr}\left\\{\mathbf{A}\bm{\Phi}\mathbf{B}\bm{\Phi}^{H}\right\\}}{\partial\bm{\theta}}\right|_{\mathbf{A}=\overline{\mathbf{H}}_{2}^{H}\bm{\Upsilon}_{k}^{2}\overline{\mathbf{H}}_{2},\mathbf{B}=\mathbf{R}_{emi}}.\end{aligned}$
(121)

Then, utilizing Lemma 4 and 5, the proof can be completed. The detailed proofs
for other terms are similar but cumbersome and therefore are omitted for
brevity. $\blacksquare$

Now, we can readily obtain the gradient vectors $\frac{\partial
f_{c}(\bm{\theta})}{\partial\bm{\theta}}$ using (112), (113), Lemmas 4, 5, 6
and the chain rule. For example, we have

$\displaystyle\frac{\partial
E_{c,k}^{\mathrm{signal}}}{\partial\bm{\theta}}=\frac{\partial\left\\{\left(E_{c,k}^{\mathrm{noise}}\right)^{2}\right\\}}{\partial\bm{\theta}}=2E_{c,k}^{\mathrm{noise}}\frac{\partial
E_{c,k}^{\mathrm{noise}}}{\partial\bm{\theta}},$ (122)

and

$\displaystyle\frac{\partial
E_{c,k}^{\mathrm{noise}}}{\partial\bm{\theta}}=M\widehat{c}_{k}\delta\bm{f}_{c,k,7}^{\prime}(\bm{\theta})+\left\\{\widehat{c}_{k}f_{c,k,2}(\bm{\Phi})+\gamma_{k}\right\\}\mathbf{z}_{k}\left(\mathbf{I}_{M}\right)+\widehat{c}_{k}\operatorname{Tr}\left\\{\bm{\Upsilon}_{k}\right\\}\bm{f}_{c,k,2}^{\prime}(\bm{\theta}).$
(123)

Algorithm 1 Gradient Ascent Algorithm

1: Initialize $\bm{\theta}_{0}$ randomly, $i=0$, $\kappa_{b}=0.5$;

2: while 1 do

3: Calculate gradient vector
$\bm{f}_{c}^{\prime}(\bm{\theta}_{i})=\left.\frac{\partial
f_{c}(\bm{\theta})}{\partial\bm{\theta}}\right|_{\bm{\theta}=\bm{\theta}_{i}}$,
$\kappa_{a}=1$;

4: while
$f_{c}(\bm{\theta}_{i}+\kappa_{a}\bm{f}_{c}^{\prime}(\bm{\theta}_{i}))<f_{c}(\bm{\theta}_{i})+\kappa_{b}\kappa_{a}\bm{f}_{c}^{\prime}(\bm{\theta}_{i})^{T}\bm{f}_{c}^{\prime}(\bm{\theta}_{i})$
do

5: $\kappa_{a}=0.8\times\kappa_{a}$;

6: end while

7:
$\bm{\theta}_{i+1}=\bm{\theta}_{i}+\kappa_{a}\bm{f}_{c}^{\prime}(\bm{\theta}_{i})$;

8: if $f_{c}(\bm{\theta}_{i+1})-f_{c}(\bm{\theta}_{i})<10^{-5}$ then

9: $\bm{\theta}^{*}=\bm{\theta}_{i+1}$, break;

10: end if

11: $i=i+1$;

12: end while

Following a similar process as (122) and after some algebraic manipulations,