# Power Scaling Law Analysis and Phase Shift Optimization of RIS-aided Massive
MIMO Systems with Statistical CSI

Kangda Zhi, Cunhua Pan, Hong Ren and Kezhi Wang (Corresponding author: Cunhua
Pan). K. Zhi, C. Pan 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@qmul.ac.uk). H. Ren is with the National Mobile
Communications Research Laboratory, Southeast University, Nanjing 210096,
China. (hren@seu.edu.cn). K. Wang is with Department of Computer and
Information Sciences, Northumbria University, UK. (e-mail:
kezhi.wang@northumbria.ac.uk).

###### Abstract

This paper considers an uplink reconfigurable intelligent surface (RIS)-aided
massive multiple-input multiple-output (MIMO) system with statistical channel
state information (CSI). The RIS is deployed to help conventional massive MIMO
networks serve the users in the dead zone. We consider the Rician channel
model and exploit the long-time statistical CSI to design the phase shifts of
the RIS, while the maximum ratio combination (MRC) technique is applied for
the active beamforming at the base station (BS) relying on the instantaneous
CSI. Firstly, we derive the closed-form expressions for the uplink achievable
rate which holds for arbitrary numbers of base station (BS) antennas. Based on
the theoretical expressions, we reveal the power scaling laws, provide the
average asymptotic rate when using random phase shifts and discuss the rate
performance under some special cases. Then, we consider the sum-rate
maximization and the minimum user rate maximization problems by optimizing the
phase shifts at the RIS. However, these two optimization problems are
challenging to solve due to the complicated data rate expression. To solve
these problems, we propose a novel genetic algorithm (GA) with low complexity
but can achieve considerable performance. Finally, extensive simulations are
provided to validate the benefits by integrating RIS into conventional massive
MIMO systems. Besides, our simulations demonstrate the feasibility of
deploying large-size but low-resolution RIS in massive MIMO systems.

###### Index Terms:

Intelligent reflecting surface (IRS), reconfigurable intelligent surface
(RIS), massive MIMO, Rician fading channels, uplink achievable rate,
statistical CSI.

##  I Introduction

The massive multiple-input multiple-output (MIMO) technology is an essential
technique to provide the extremely high network throughput in current and
future communication systems[1]. However, to achieve such high throughput,
hundreds of antennas should be equipped at the base station (BS), which raises
the issues of high cost and energy consumption. Besides, to provide seamless
coverage in the urban environment with dense obstacles, active relay or small
BSs should be densely deployed, which also sharply increases the cost. On the
other hand, reconfigurable intelligent surface (RIS), also known as
intelligent reflecting surface (IRS), has been proposed as a revolutionary
technology to support high data rate while maintaining at a low cost and
energy consumption[2, 3, 4]. Specifically, RIS can constructively reflect the
signal from the multi-antenna BS to multiple users which cannot directly
communicate with the BS due to blockages, and it can also operate in a full-
duplex (FD) mode without self-interference. Therefore, RIS is an efficient and
cost-effective solution for the blockage problem of conventional massive MIMO
systems.

Due to the above advantages, RIS-aided communication systems have been widely
investigated in various scenarios[5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24, 25]. Specifically, for single-cell multi-
antenna systems, the authors in [5] jointly considered the active and passive
beamforming optimizations to demonstrate the potential of RIS, while a deep
reinforcement learning-based method was proposed in [6]. Zhang _et al._ [7]
characterized the fundamental capacity limit of RIS-aided MIMO systems with
the narrowband and broadband transmission. Downlink multigroup multicast
communication systems were presented in [8] and the RIS-aided simultaneous
wireless information and power transfer (SWIPT) systems were studied in [9].
The benefits of using RIS in multi-user FD two-way communication systems were
demonstrated in [10]. Meanwhile, an energy efficiency maximization problem was
considered in [11]. To investigate the performance of RIS-aided multi-cell
MIMO networks, the authors in [12] proposed to deploy an RIS at the cell edge
and demonstrated the benefits of the RIS to mitigate the inter-cell
interference. Furthermore, RIS-aided mobile edge computing (MEC) systems were
studied in [13], which showed that significant latency can be reduced by
integrating RIS into conventional MEC systems. The authors in [14] further
investigated the wireless powered orthogonal-frequency-division-multiplexing
(OFDM) MEC systems under the assistance of an RIS. Meanwhile, RIS-aided
unmanned aerial vehicle (UAV) networks were studied in [15, 16, 17].
Specifically, the work in [15] considered the joint optimization of UAV’s
trajectory and RIS’s phase shifts in a single-user network, and a novel
symbiotic UAV-aided multiple RIS radio system was studied in [16]. Wang _et
al._ [17] further investigated the UAV-aided multi-RIS multi-user systems
using a deep reinforcement learning approach. Taking into consideration the
impact of imperfect cascaded channels, the authors in [18] firstly studied the
robust active and passive beamforming optimization problem to minimize the
total transmit power. Besides, RIS-aided space shift keying and RIS-aided
spatial modulation schemes were investigated in [19, 20]. Considering the
secure communication scenarios, the authors in[21] studied the performance of
artificial noise-aided MIMO systems with the aid of an RIS. RIS-aided secure
communications with imperfect RIS-eavesdropper channels were considered in
[22], while the authors in [23] further investigated the robust transmission
design in RIS-aided secure communications with cascaded channel error.
Furthermore, RIS-aided MIMO and FD cognitive radio systems were respectively
studied in [24] and [25].

However, all of the above contributions considered to design the phase shifts
of the RIS based on instantaneous channel state information (CSI). Those
schemes are suitable for the scenarios with a fixed location or low mobility,
which enable the BS to carry out the channel estimation, design the optimal
RIS phase shifts and adjust the phase shifts of the RIS in each channel
coherence time. However, for the scenarios with high mobility and short
channel coherence time, it is more practical to design and tune the phase
shifts of the RIS relying on statistical CSI. Furthermore, this statistical
CSI-based strategy can effectively reduce the feedback overhead required for
RIS[4], reduce the power consumed by RIS’s controller and release the capacity
requirement for the RIS’s control link. In addition, significant computational
complexity can be reduced at the BS since the phase shift matrix is only
needed to be updated when the statistical CSI varies, which occurs in a much
larger time scale than the instantaneous CSI.

Due to the above benefits, some researchers have exploited the statistical CSI
to design the RIS-aided communication systems [26, 27, 28, 29, 30, 31, 32, 33,
34]. For the single-user systems, Han _et al._ [26] first presented the
optimal RIS phase shift design based on the derived ergodic capacity
expression under the Rician channel model. The authors in [27] further
designed the RIS-aided systems with a serving BS and an interfered BS. For the
multi-user case, Peng _et al._ [28] investigated the performance of RIS-aided
multi-pair communication systems and verified the effectiveness of applying
genetic algorithm (GA) in the optimization of the phase shifts of the RIS. The
performance of RIS-aided Internet of Things under correlated Rayleigh channels
was evaluated in [29]. The authors in [30] proposed a location information-
aided multi-RIS system, where a low-complexity BS-RIS maximum-ratio
transmission beamforming scheme was proposed. By resorting to random matrix
theory, You _et al._[31] considered the energy efficiency maximization problem
in MIMO networks under the correlated Rayleigh channel model, and Nadeem _et
al._ [32] considered the minimum signal-to-interference-plus-noise ratio
(SINR) maximization problem with line-of-sight (LoS) BS-RIS channel matrix. A
novel two-timescale beamforming optimization scheme was proposed in [33],
where the passive beamforming was first optimized based on statistical CSI and
then the active beamforming was designed based on instantaneous CSI. Besides,
the authors in [34] studied the IRS-aided opportunistic beamforming scheme
with statistical CSI.

However, based on the statistical CSI, the RIS-aided massive MIMO systems
under the Rician channel model have not been investigated. On one hand, since
the RIS is often deployed on the facade of tall buildings, the RIS-related
channels may possess the LoS channel components. Therefore, the more general
Rician fading model should be adopted. On the other hand, it is crucial to
characterize the interplay between the promising RIS technology and the
existing massive MIMO technology, and evaluate the potential of RIS-aided
massive MIMO systems. To the best of our knowledge, only [35] studied the RIS-
aided massive MIMO networks. However, in [35], the correlated Rayleigh channel
was considered and the phase shifts of RIS are simply set as an identity
matrix.

Against the above background, in this paper, we theoretically analyze and
optimize the uplink RIS-aided massive MIMO systems with the Rician channel
model and statistical CSI. Specifically, the low-complexity maximum-ratio
combination (MRC) technique is employed for the active beamforming based on
the instantaneous CSI, while the phase shifts of the RIS are designed and
adjusted by exploiting the statistical CSI. The Rician channel model is
applied in this paper to capture the achievable spatial multiplexing gain of
RIS-aided massive MIMO systems. We present the closed-form analytical
expression for the uplink achievable rate which holds for arbitrary numbers of
antennas at the BS. Our main contributions are summarized as follows:

  * •

First, we derive the closed-form expression of the uplink achievable rate
using the Rician channel model that holds for any finite number of antennas at
the BS, and this analytical expression only depends on the locations and
angles information and Rician factors. Based on the derived expressions, we
reveal the scaling laws with respect to the number of RIS’s elements and the
number of BSs’ antennas. We also evaluate the average asymptotic rate achieved
by random phase shifts.

  * •

Then, by using the derived expression, we utilize the GA-based method to solve
the sum-rate maximization problem and the minimum user rate maximization
problem, by taking into consideration the impact of discrete phase shifts.

  * •

Finally, extensive simulations are carried out to characterize the gains by
employing RIS into massive MIMO networks. Our results reveal the trade-off
between the increase of spatial multiplexing gain and the decrease of path
loss in the RIS-aided massive MIMO systems. Meanwhile, we validate the
feasibility of deploying large-size RIS with low-resolution hardware into
existing massive MIMO systems.

The remainder of this paper is organized as follows. Section II describes the
model of uplink RIS-aided massive MIMO systems with Rician channel. Section
III derives the closed-form analytical expressions for the uplink achievable
rate with arbitrary numbers of BS antennas, and discusses the power scaling
laws and some special cases. Section IV presents the GA-based method to solve
the sum-rate maximization and the minimum user rate maximization problems.
Section V provides extensive simulation results to characterize the achievable
spatial multiplexing gain and other benefits brought by RIS. Finally, Section
VI concludes this paper.

_Notations_ : The vectors and the matrices are respectively expressed in
lowercase blodface and uppercase blodface letters. $\mathbf{A}^{H}$,
$\mathbf{A}^{T}$ and $\mathbf{A}^{*}$ represent the conjugate transpose,
transpose and conjugate operators, respectively. $\left|\mathbf{a}\right|$
denotes the modulus of the complex number and $\left\|\mathbf{a}\right\|$
denotes $l_{2}$-norm of the vector. $\left[\mathbf{A}\right]_{m,n}$ denotes
the $\left(m,n\right)$-th entry of the matrix. $\left[\mathbf{a}\right]_{m}$
denotes the $m$-th entry of the vector. Re$\left\\{\cdot\right\\}$ represents
the real part. $\mathbb{E}\left\\{\cdot\right\\}$ and
$\mathrm{Tr}\left\\{\cdot\right\\}$ denote the expectation and trace operator,
respectively. $\mathbf{I}_{N}$ is the identity matrix with $N$ dimension.
$\mathbb{C}^{M\times N}$ represents the $M\times N$ complex-valued matrix.
Besides, $x\sim{\cal C}{\cal N}\left({a,b}\right)$ denotes that random
variable $x$ follows the complex Gaussian distribution with mean $a$ and
variance $b$. Operation $\left\lfloor n\right\rfloor$ means rounding $n$
toward the negative infinity and operation $\bmod$ means taking the remainder
after division.

##  II System Model

Figure 1: An RIS-assisted uplink massive MIMO communication system.

Consider a typical uplink RIS-aided communication system with a base station
(BS), an RIS and $K$ single-antenna users. The BS and RIS are equipped with
$M$ antennas and $N$ reflecting elements, respectively. The RIS is connected
with the BS with a dedicated transmission link and its phase shifts can be
controlled by the BS.

Since the ground communication links can be obstructed by buildings, humans
and trees, we assume that the direct links between the BS and users are
blocked similar to [32, 30, 31]. To assist the communications for users, an
RIS is deployed on the building and helps users communicate with the BS, as
shown in Fig. 1. Considering the fact that the RIS can be deployed on the wall
of tall buildings, it can assist in creating channels dominated by LoS
propagation along with a few scatters. Therefore, we adopt the Rician fading
model, and the channels between users and the RIS and that between the RIS and
the BS can be respectively expressed as:

$\displaystyle{{\bf{H}}_{1}}=\left[{{{\bf h}_{1}},{{\bf h}_{2}},...,{{\bf
h}_{K}}}\right],{{\bf
h}_{k}}=\sqrt{{\alpha_{k}}}\left({\sqrt{\frac{{{\varepsilon_{k}}}}{{{\varepsilon_{k}}+1}}}{{\bf\bar{h}}_{k}}+\sqrt{\frac{1}{{{\varepsilon_{k}}+1}}}{{\bf\tilde{h}}_{k}}}\right),$
(1) $\displaystyle{{\bf
H}_{2}}=\sqrt{\beta}\left({\sqrt{\frac{{{\delta}}}{{{\delta}+1}}}{{\bf\bar{H}}_{2}}+\sqrt{\frac{1}{{\delta+1}}}{{\bf\tilde{H}}_{2}}}\right),$
(2)

where ${{\bf H}_{1}}\in{\mathbb{C}^{N\times K}}$, ${{\bf
H}_{2}}\in{\mathbb{C}^{M\times N}}$, ${\alpha_{k}}$ and $\beta$ are the large-
scale path loss of the user $k$-RIS link and RIS-BS link, respectively.
${\varepsilon_{k}}$ and ${\delta}$ are the Rician factors of the user $k$-RIS
link and RIS-BS link, respectively.
${{\bf\bar{h}}_{k}}\in{{\mathbb{C}}^{N\times 1}}$ and
${{\bf\bar{H}}_{2}}\in{{\mathbb{C}}^{M\times N}}$ denote the LoS components of
the user $k$-RIS link and RIS-BS link.
${{\bf\tilde{h}}_{k}}\in{{\mathbb{C}}^{N\times 1}}$ and
${{\bf\tilde{H}}_{2}}\in{{\mathbb{C}}^{M\times N}}$ represent the non-line-of-
sight (NLoS) components, whose elements are independently and identically
distributed (i.i.d) complex Gaussian random variables following the
distribution of ${\cal C}{\cal N}\left({0,1}\right)$.

Assume that the BS and RIS are equipped with uniform square planar array
(USPA) with size of $\sqrt{M}\times\sqrt{M}$ and $\sqrt{N}\times\sqrt{N}$,
respectively. Therefore, the LoS components
${{\bf\bar{h}}_{k}}\in{{\mathbb{C}}^{N\times 1}}$ and
${{\bf\bar{H}}_{2}}\in{{\mathbb{C}}^{M\times N}}$ can be respectively
expressed as:

$\displaystyle{{\bf\bar{h}}_{k}}={{\bf
a}_{N}}\left({\varphi_{kr}^{a},\varphi_{kr}^{e}}\right),$ (3)
$\displaystyle{{\bf\bar{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),$ (4)

with array response vector as

$\displaystyle{{\bf
a}_{X}}\left({\vartheta^{a},\vartheta^{e}}\right)={\left[{1,...,{e^{j2\pi\frac{d}{\lambda}\left({x\sin\vartheta^{a}\sin\vartheta^{e}+y\cos\vartheta^{e}}\right)}},...,{e^{j2\pi\frac{d}{\lambda}\left({\left({\sqrt{X}-1}\right)\sin\vartheta^{a}\sin\vartheta^{e}+\left({\sqrt{X}-1}\right)\cos\vartheta^{e}}\right)}}}\right]^{T}},$
(5)

where ${\rm{0}}\leq x,y\leq\sqrt{X}-1$, $d$ and $\lambda$ are the element
spacing and carrier wavelength, $\varphi_{kr}^{a}$ and $\varphi_{kr}^{e}$ are
respectively the azimuth and elevation angles of arrival (AoA) at the RIS from
user $k$. $\varphi_{t}^{a}$ and $\varphi_{t}^{e}$ respectively denote the
azimuth and elevation angles of departure (AoD) from the RIS towards the BS.
$\phi_{r}^{a}$ and $\phi_{r}^{e}$ respectively represent the AoA at the BS
from the RIS. Note that ${{\bf\bar{h}}_{k}}$ and ${{\bf\bar{H}}_{2}}$ only
rely on the AoA and AoD, which could keep invariant within the considered time
period. Besides, we assume that these angles are known based on some
technologies. For example, it can be calculated by the locations obtained from
the global position system (GPS).

With the help of RIS, the received signal at the BS can be written as:

$\displaystyle{\bf y}={\bf
G\mathbf{P}x+n}={\bf{H_{2}}\Phi{H_{1}}\mathbf{P}x+n},$ (6)

where
$\mathbf{P}={\rm{diag}}\left\\{\sqrt{p_{1}},\sqrt{p_{2}},...,\sqrt{p_{K}}\right\\}$,
$p_{k}$ is the transmit power of user $k$. ${{{\bf
G}\buildrel\Delta\over{=}\bf{H_{2}}\Phi{H_{1}}}}\in{{\mathbb{C}}^{M\times K}}$
represents the cascaded user-RIS-BS channel.
${\bf\Phi}={\rm{diag}}\left\\{{{e^{j{\theta_{1}}}},{e^{j{\theta_{2}}}},...,{e^{j{\theta_{N}}}}}\right\\}$
is the reflection matrix of RIS and ${\theta_{n}}\in[0,2\pi)$ is the phase
shift introduced by the RIS reflector $n$. ${\bf
x}={\left[{{x_{1}},{x_{2}},...,{x_{K}}}\right]^{T}}\in{{\mathbb{C}}^{K\times
1}}$ denotes the information symbols from $K$ users, where
$\mathbb{E}\left\\{\left|x_{k}\right|^{2}\right\\}=1$. ${\bf n}\sim{\cal
C}{\cal N}\left({0,{\sigma^{2}}{\bf I}}\right)$ is the additional white
Gaussian noise (AWGN).

Adopting the maximal-ratio-combining (MRC) technique, the received signal at
the BS can be written as

$\displaystyle{{\bf r}={{\bf G}^{H}}{\bf y}}={{{\bf G}^{H}}{\bf
G\mathbf{P}x}+{{\bf G}^{H}}\bf n},$ (7)

and the signal of user $k$ can be expressed as

$\displaystyle{r_{k}}=\sqrt{p}_{k}{{\bf g}_{k}^{H}}{{\bf
g}_{k}}{x_{k}}+\sum\limits_{i=1,i\neq k}^{K}{\sqrt{p}_{i}{{{\bf
g}_{k}^{H}}}{{\bf g}_{i}}{x_{i}}}+{{{\bf g}_{k}^{H}}}{\bf n}.$ (8)

where ${{\bf g}_{k}}\buildrel\Delta\over{=}{{\bf H}_{2}}{\bf\Phi}{{\bf
h}_{k}}\in{{\mathbb{C}}^{M\times 1}}$ is the $k$-th column of matrix $\bf G$
representing the cascaded user $k$-RIS-BS channel.

Considering the ergodic channel, the uplink achievable rate of user $k$ can be
expressed as

$\displaystyle
R_{k}=\mathbb{E}\left\\{\log_{2}\left(1+\frac{{p_{k}}\left\|\mathbf{g}_{k}\right\|^{4}}{\sum\limits_{i=1,i\neq
k}^{K}{p_{i}}\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}+\sigma^{2}\left\|\mathbf{g}_{k}\right\|^{2}}\right)\right\\},$
(9)

and the sum rate is as

$\displaystyle R=\sum\limits_{i=1}^{K}{{R_{k}}}.$ (10)

##  III Uplink Achievable Rate Analysis

In this section, we derive the closed-form expression of the achievable rate
in the uplink RIS-aided multi-user system. The theoretical results can capture
the impacts of various variables, including the number of antennas at the BS,
the number of reflecting elements at the RIS, the transmit power and Rician
factors. We will also present asymptotic expressions in some special cases.

###  III-A Preliminary Results

We first give a key Lemma which will be used in further derivations.

######  Lemma 1.

The expectation of ${{{\left\|{{{\bf g}_{k}}}\right\|}^{2}}}$,
${{{\left\|{{{\bf g}_{k}}}\right\|}^{4}}}$ and ${{{\left|{{\bf
g}_{{}_{k}}^{H}{{\bf g}_{i}}}\right|}^{2}}}$ are respectively given by

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}=\frac{M\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\left(\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}+\left(\delta+\varepsilon_{k}+1\right)N\right),\qquad\qquad\qquad\qquad\quad$
(11)
$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}=M\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\times\\\
\left\\{M\left(\delta\varepsilon_{k}\right)^{2}\left|f_{k}({\bf\Phi})\right|^{4}+2\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}\left(2MN\delta+MN\varepsilon_{k}+MN+2M+N\varepsilon_{k}+N+2\right)\right.\\\
+MN^{2}\left(2\delta^{2}+\varepsilon_{k}^{2}+2\delta\varepsilon_{k}+2\delta+2\varepsilon_{k}+1\right)+N^{2}\left(\varepsilon_{k}^{2}+2\delta\varepsilon_{k}+2\delta+2\varepsilon_{k}+1\right)\\\
\left.+MN\left(2\delta+2\varepsilon_{k}+1\right)+N\left(2\delta+2\varepsilon_{k}+1\right)\right\\},\end{array}$
(16)

and

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}=M\frac{\beta^{2}\alpha_{i}\alpha_{k}}{(\delta+1)^{2}\left(\varepsilon_{i}+1\right)\left(\varepsilon_{k}+1\right)}\times\left\\{M\delta^{2}\varepsilon_{k}\varepsilon_{i}\left|f_{k}({\bf\Phi})\right|^{2}\left|f_{i}({\bf\Phi})\right|^{2}\right.\\\
+\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}\left(\delta
MN+N\varepsilon_{i}+N+2M\right)+\delta\varepsilon_{i}\left|f_{i}({\bf\Phi})\right|^{2}\left(\delta
MN+N\varepsilon_{k}+N+2M\right)\\\
+N^{2}\left(M\delta^{2}+\delta\left(\varepsilon_{i}+\varepsilon_{k}+2\right)+\left(\varepsilon_{k}+1\right)\left(\varepsilon_{i}+1\right)\right)+MN\left(2\delta+\varepsilon_{i}+\varepsilon_{k}+1\right)\\\
\left.+M\varepsilon_{k}\varepsilon_{i}\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}+2M\delta\varepsilon_{k}\varepsilon_{i}\operatorname{Re}\left\\{f_{k}^{H}({\bf\Phi})f_{i}({\bf\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\}\right\\},\end{array}$
(21)

where ${f_{c}}\left({\bf\Phi}\right)\in{{\mathbb{C}}^{1\times
1}},c\in\left\\{k,i\right\\}$ is defined as

$\displaystyle{f_{c}}\left({\bf\Phi}\right)\buildrel\Delta\over{=}{\bf
a}_{N}^{H}\left({\varphi_{t}^{a},\varphi_{t}^{e}}\right){\bf\Phi}\overline{\bf
h}_{c}=\sum\limits_{n=1}^{N}{{e^{j2\pi\frac{d}{\lambda}\left({x{p_{c}}+y{q_{c}}}\right)+j{\theta_{n}}}}},$
(22)

with $x=\left\lfloor{\left({n-1}\right)/\sqrt{N}}\right\rfloor$,
$y=\left({n-1}\right)\bmod\sqrt{N}$,
${p_{c}}=\sin\varphi_{cr}^{a}\sin\varphi_{cr}^{e}-\sin\varphi_{t}^{a}\sin\varphi_{t}^{e}$,
${q_{c}}=\cos\varphi_{cr}^{e}-\cos\varphi_{t}^{e}$.

Proof: Please refer to Appendix A. $\blacksquare$

Note that $\left|{f_{c}}\left({\bf\Phi}\right)\right|\leq N$ and the equality
holds when
$\theta_{n}=-2\pi\frac{d}{\lambda}\left({x{p_{c}}+y{q_{c}}}\right),\forall n$.
In this paper, we refer to the phase shift solution that maximizes
$\left|{f_{k}}\left({\bf\Phi}\right)\right|$ as “phase aligned to user $k$”.
In this setting, when $N\rightarrow\infty$,
$\left|{f_{k}}\left({\bf\Phi}\right)\right|$ can grow without bound. However,
$\left|{f_{i}}\left({\bf\Phi}\right)\right|,i\neq k$ will be bounded unless
user $i$ has nearly the same azimuth and elevation AoA with user $k$. Note
that we ignore this rare situation in this section.

Lemma 1 shows that both
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$ and
$\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$
are on the order of $\mathcal{O}\left(M^{2}\right)$. However, their scaling
laws with respect to $N$ depends on the value of $\bf\Phi$. For example, when
the phase shifts of RIS are aligned to user $k$,
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$ is on the order
of $\mathcal{O}\left(N^{4}\right)$ whereas
$\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$
is on the order of $\mathcal{O}\left(N^{3}\right)$.

###  III-B Main Results

Next, with the above results (11)$\sim$ (21), we provide the closed-form
expression of the uplink achievable rate under the general case with any
number of antennas.

######  Theorem 1.

In the RIS-aided massive MIMO systems, the uplink achievable rate of user $k$
can be approximated as

$\displaystyle
R_{k}\approx\log_{2}\left(1+\frac{p_{k}{E}_{k}^{(signal)}({\bf\Phi})}{\sum\limits_{i=1,i\neq
k}^{K}p_{i}I_{ki}({\bf\Phi})+\sigma^{2}E_{k}^{(noise)}({\bf\Phi})}\right),$
(23)

where
$E_{k}^{(signal)}({\bf\Phi})\triangleq\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$,
$I_{ki}({\bf\Phi})\triangleq\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$,
and
$E_{k}^{(noise)}({\bf\Phi})\triangleq\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}$.

Proof: It can be readily proved by using Jensen’s inequality as in [36, Lemma
1]. $\blacksquare$

Rate expression (23) characterizes the impacts of $\bf\Phi$, $M$, $N$, $P$,
different kinds of AoA and AoD, path-loss parameters and Rician factors on the
data rate performance. We can see that this theoretical expression is only
determined by locations, AoA and AoD of the BS, the RIS and users, which could
keep invariant for a long time. Therefore, designing the phase shifts of RIS
based on statistical CSI can significantly reduce the computational complexity
and channel estimation overhead in practical systems.

######  Corollary 1.

In the RIS-aided single user systems, i.e., without the multi-user
interference, the achievable rate of user $k$ is

$\displaystyle
R_{k}\approx\log_{2}\left(1+\frac{p_{k}E_{k}^{(signal)}({\bf\Phi})}{\sigma^{2}E_{k}^{(noise)}({\bf\Phi})}\right),$
(24)

which can achieve the gain of
$\mathcal{O}\left(\mathrm{log_{2}}(MN^{2})\right)$.

It is well known that this performance gain comes from the active beamforming
gain of multi-antenna, passive beamforming gain of RIS and the inherent
aperture gain of RIS[5]. However, when considering the multi-user
interference, this performance gain cannot be obtained. We can see that both
$E_{k}^{(signal)}({\bf\Phi})$ and $I_{ki}({\bf\Phi})$ in (23) increase on the
order of $\mathcal{O}\left(M^{2}\right)$, therefore it cannot achieve a gain
of $\mathcal{O}\left(\mathrm{log_{2}}(M)\right)$, i.e., the rate in (23)
cannot grow without bound when $M\rightarrow\infty$. This is because the
channels of different users share the same RIS-BS channel $\mathbf{H}_{2}$.
For example, recalling (165) and (182) in Appendix A, we can see that the
common term $\overline{\mathbf{H}}_{2}$ and $\tilde{\mathbf{H}}_{2}$ bring the
factor $M^{2}$ to the interference term $I_{ki}({\bf\Phi})$.

Meanwhile, rate expression (23) shows that the order of magnitude with respect
to $N$ is determined by the setting of $\mathbf{\Phi}$. If we adjust RIS’s
phase shifts to maximize the desired signal power of user $k$, i.e.,
$\left|{f_{k}}\left({\bf\Phi}\right)\right|=N$, we can find that the rate
$R_{k}$ is on the order of $\mathcal{O}\left(\mathrm{log_{2}}(N)\right)$. This
means that the sum achievable rate can grow without bound when
$N\rightarrow\infty$. Although the rate does not have the order of
$\mathcal{O}\left(\mathrm{log_{2}}(N^{2})\right)$ in this simple case, the
rate performance could be improved by properly designing the phase shifts of
RIS to increase the desired signal as well as mitigating the interference, and
RIS’s interference mitigation capability has been validated in the previous
contributions [27, 12]. This discussion emphasizes the importance of the phase
shift design in RIS-aided massive MIMO systems.

To better understand Theorem 1, we will present the asymptotic results under
some special cases. Firstly, we reveal the power-scaling law of the uplink
achievable rate in the RIS-aided massive MIMO systems in the following.

######  Corollary 2.

Assume that the transmit power of each user is scaled with the number of
antennas at the BS according to $p_{k}={E_{u}}/{M}$, $\forall k$, where
${E_{u}}$ is fixed. When $M\rightarrow\infty$, we have

$\displaystyle
R_{k}\to\log_{2}\left(1+\frac{E_{u}\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\mathcal{A}_{k}^{(1)}({\bf\Phi})}{\sum_{i=1,i\neq
k}^{K}E_{u}\frac{\beta\alpha_{i}}{(\delta+1)\left(\varepsilon_{i}+1\right)}{\mathcal{A}^{(2)}_{ki}}({\bf\Phi})+\sigma^{2}\mathcal{A}_{k}^{(3)}({\bf\Phi})}\right),$
(25)

where

$\displaystyle\mathcal{A}_{k}^{(1)}({\bf\Phi})$
$\displaystyle=\left(\mathcal{A}_{k}^{(3)}({\bf\Phi})\right)^{2}+2\delta\varepsilon_{k}\left|f_{k}(\mathbf{\Phi})\right|^{2}(N\delta+2)+N\left(N\delta^{2}+2\delta+2\varepsilon_{k}+1\right),$
(26) $\displaystyle\mathcal{A}_{ki}^{(2)}({\bf\Phi})$
$\displaystyle=\varepsilon_{k}\varepsilon_{i}\left|\delta
f_{k}^{H}(\mathbf{\Phi})f_{i}(\mathbf{\Phi})+\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}$
$\displaystyle\;\;\;+\left(\delta^{2}N+2\delta\right)\left(\varepsilon_{k}\left|f_{k}(\mathbf{\Phi})\right|^{2}+\varepsilon_{i}\left|f_{i}(\mathbf{\Phi})\right|^{2}\right)+N\left(N\delta^{2}+2\delta+\varepsilon_{i}+\varepsilon_{k}+1\right),$
(27) $\displaystyle\mathcal{A}^{(3)}_{k}({\bf\Phi})$
$\displaystyle=\delta\varepsilon_{k}\left|f_{k}(\mathbf{\Phi})\right|^{2}+\left(\delta+\varepsilon_{k}+1\right)N.$
(28)

Proof: By substituting $p_{k}=E_{u}/M,\forall k$ into rate expression (23),
when $M\to\infty$, we can ignore the insignificant terms which don’t scale
with $M$. Then, after some simplifications, we can complete the proof.
$\blacksquare$

From Corollary 2, we can see that similar to traditional massive MIMO systems,
users in RIS-aided systems can cut down their transmit power by a factor $1/M$
while the rate will converge to a non-zero value as $M\to\infty$. However,
different from the traditional systems, both the signal, interference and
noise terms in rate (25) are related with $\bf\Phi$. To clearly show the
difference, we consider a special case where the RIS is deployed in the
environment with pure NLoS channels, i.e., $\delta=\varepsilon_{k}=0,\forall
k$. Then, the power scaling law in Corollary 2 becomes

$\displaystyle\tilde{R}_{k}\rightarrow\log_{2}\left(1+\frac{E_{u}\beta\alpha_{k}(N+1)}{\sum_{i=1,i\neq
k}^{K}E_{u}\beta\alpha_{i}+\sigma^{2}}\right),\text{as }M\to\infty.$ (29)

By contrast, in traditional non-RIS massive MIMO systems with large-scale
path-loss $\gamma_{k}$, when scaling the power by $p_{k}=E_{u}/M$, the rate
can be written as [36, Theorem 1]:

$\displaystyle{R}^{(w/o)}_{k}\rightarrow\log_{2}\left(1+\frac{E_{u}\gamma_{k}}{\sigma^{2}}\right),\text{as
}M\to\infty.$ (30)

Comparing Eq. (29) with Eq. (30), we can see that the rate can reap
significant benefits by deploying RIS with large number of elements in massive
MIMO systems.

######  Corollary 3.

If the phase shifts of RIS are aligned to user $k$, the transmit power of user
$k$ is scaled down by $p_{k}=\frac{E_{u}}{MN^{2}}$, while the transmit power
of other users are scaled down by $p_{i}=\frac{E_{u}}{MN},\forall i\neq k$.
When both $M$ and $N$ are large, we have

$\displaystyle
R_{k}\rightarrow\log_{2}\left(1+\frac{E_{u}\frac{\varepsilon_{k}}{\left(\varepsilon_{k}+1\right)}}{E_{u}\sum_{i=1,i\neq
k}^{K}\frac{\alpha_{i}}{\left(\varepsilon_{i}+1\right)\alpha_{k}}+\left(1+\frac{1}{\delta}\right)\frac{\sigma^{2}}{\beta\alpha_{k}}}\right),$
(31) $\displaystyle R_{i}\rightarrow 0,\forall i\neq k,$ (32)

Proof: Please refer to Appendix B. $\blacksquare$

Corollary 3 means that with large $M$ and $N$, we can further cut down user
$k$’s transmission power to $E_{u}/(MN^{2})$ while keeping the data rate as a
non-zero value. Meanwhile, this rate will be improved if the environment has
few scatters, i.e., with a larger $\alpha_{k}$, $\beta$ and $\delta$.

######  Corollary 4.

For both the ideal RIS with continuous phase shifts and non-ideal RIS with
$b>1$ bits discrete phase shifts, if the phase shift matrix $\bf\Phi$ is
randomly adjusted in each time block, when $N\rightarrow\infty$ and
$M\rightarrow\infty$, we have

$\displaystyle
R_{k}\rightarrow\log_{2}\left(1+\frac{p_{k}\alpha_{k}\left(2\delta^{2}+2\delta+1\right)}{\sum_{i=1,i\neq
k}^{K}p_{i}\alpha_{i}\delta^{2}}\right).$ (33)

Proof: Please refer to Appendix C. $\blacksquare$

Corollary 4 shows that with a large number of antennas at the BS and a large
number of reflecting elements at the RIS, the sum achievable rate is still
bounded if the phase shifts are randomly adjusted. This conclusion shows the
necessity of optimizing the phase shifts of RIS in the RIS-aided massive MIMO
systems. Besides, we can see that the data rate in (33) decreases when
$\delta$ increases, which has a different trend from Corollary 3. The reason
lies in that when the phase shifts are adjusted randomly in each time block,
it tends to equally allocate the passive beamforming gain to all the users.
However, when $\delta\rightarrow\infty$, the channel with unit rank will be
unable to support the multi-user communications.

######  Corollary 5.

If $\delta=\varepsilon_{k}=0,\forall k$, i.e., only NLoS paths exist in the
environment, we have

$\displaystyle\tilde{R}_{k}\rightarrow\log_{2}\left(1+\frac{p_{k}\beta\alpha_{k}(MN+M+N+1)}{\sum_{i=1,i\neq
k}^{K}p_{i}\beta\alpha_{i}(M+N)+\sigma^{2}}\right).$ (34)

Proof: The proof can be completed by removing the terms with zero values when
setting $\delta=\varepsilon_{k}=0,\forall k$. $\blacksquare$

Corollary 5 represents the environment where rich scatters exist and the
Rician channel degrades to the Rayleigh channel. We can see that with
uncorrelated Rayleigh channel, there is no need to design the phase shifts of
RIS. Therefore, in the environment with rich scatters, the phase shifts of RIS
can be set arbitrarily. Besides, with a large number of antennas or a large
number of reflecting elements, the rate in (34) will converge to

$\displaystyle\tilde{R}_{k}\rightarrow\log_{2}\left(1+\frac{p_{k}\alpha_{k}(N+1)}{\sum_{i=1,i\neq
k}^{K}p_{i}\alpha_{i}}\right),\mathrm{as}\;M\rightarrow\infty,$ (35)
$\displaystyle\tilde{R}_{k}\rightarrow\log_{2}\left(1+\frac{p_{k}\alpha_{k}(M+1)}{\sum_{i=1,i\neq
k}^{K}p_{i}\alpha_{i}}\right),\mathrm{as}\;N\rightarrow\infty.$ (36)

Therefore, even the LoS link does not exist, significant performance gain can
be achieved by deploying RIS with large numbers of elements in the massive
MIMO systems.

######  Corollary 6.

When $\delta=\varepsilon_{k}\rightarrow\infty,\forall k$, i.e., only LoS paths
exist, we have

$\displaystyle\bar{R}_{k}\rightarrow\log_{2}\left(1+\frac{p_{k}\beta\alpha_{k}M\left|f_{k}({\bf\Phi})\right|^{2}}{\sum_{i=1,i\neq
k}^{K}p_{i}{\beta\alpha_{i}M\left|f_{i}({\bf\Phi})\right|^{2}}+\sigma^{2}}\right).$
(37)

By contrast, in the conventional massive MIMO systems without RIS, the rate
under LoS channel $\sqrt{\gamma_{k}}\mathbf{\bar{h}}_{k}^{w/o}$ is

$\displaystyle\bar{R}_{k}^{w/o}=\log_{2}\left(1+\frac{p_{k}\gamma_{k}M}{\sum_{i=1,i\neq
k}^{K}p_{i}\gamma_{i}\frac{\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}}{M}+\sigma^{2}}\right).$
(38)

Proof: Please refer to Appendix D. $\blacksquare$

Corollary 6 clearly presents the difference between RIS-aided massive MIMO
systems and conventional non-RIS massive MIMO systems. We can see that in
conventional uplink massive MIMO systems without fast fading, when the number
of antennas is large, the multi-user interference term will be zero compared
with the useful signal power. However, this property does not hold in the RIS-
aided massive MIMO systems with a low-complexity MRC scheme. However, this
rate degradation can be compensated by properly designing phase shifts
$\mathbf{\Phi}$. For example, when the phase shifts are aligned to user $k$,
the sum inter-user interference suffered by user $k$ will become negligible
compared with the desired signal received by user $k$. This observation
emphasizes the importance of the optimization of $\bf\Phi$.

######  Corollary 7.

In the RIS-aided massive MIMO systems, if the RIS has discrete phase shifts
with $b$ bits resolution, the sum achievable rate can still achieve a gain of
$\mathcal{O}\left(\mathrm{log_{2}}(N)\right)$.

Proof: Please refer to Appendix E. $\blacksquare$

Corollary 7 states that the phase noise of RIS does not impact the scaling
laws, and the rate can still grow without bound when $N\rightarrow\infty$.
Corollary 7 indicates that the negative effect brought by RIS’s low-resolution
elements can be easily compensated by increasing the size of RIS. Therefore,
this result demonstrates the feasibility of deploying low-resolution RIS with
a large number of reflecting elements in the massive MIMO systems.

##  IV Phase Shifts Design

In this section, we will design the phase shifts of RIS based on the long-term
statistical CSI, which could effectively reduce the training overhead and the
frequency of updating phase shifts at the RIS. For the RIS-aided massive MIMO
systems, we respectively formulate two optimization problems with different
objective functions, and both continuous and discrete phase shifts of RIS are
considered. The sum user rate-oriented optimization problem which can
characterize the system capacity limitation is formulated as

$\max\limits_{\mathbf{\Phi}}\;\;\sum\limits_{k=1}^{K}R_{k},\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
(39a) $\text{s.t. }\;\;\theta_{n}\in[0,2\pi),\forall
n,\text{or}\qquad\quad\qquad\qquad\qquad$ (39b)
$\qquad\qquad{\theta}_{n}\in\left\\{0,\frac{2\pi}{2^{b}},2\times\frac{2\pi}{2^{b}},\ldots,\left(2^{b}-1\right)\frac{2\pi}{2^{b}}\right\\},\forall
n,$ (39c)

where $R_{k}$ is given by (23). Note that constraint (39b) corresponds to the
continuous phase shift case, while constraint (39c) corresponds to the
discrete phase shift case with $b$ bits precision.

Next, the minimum user rate-oriented optimization problem which could
guarantee fairness and characterize networks spatial multiplexing is
formulated as follows

$\max\limits_{\mathbf{\Phi}}\;\;\min\limits_{k}\;\;R_{k},\qquad\;$ (40a)
$\text{s.t.
}\;\;(\text{\ref{constraint1}})\;\text{or}\;\text{(\ref{constraint2})}.$

###  IV-A Special Cases

To begin with, we will discuss phase shifts design in some special cases.

######  Proposition 1.

For problem (39) and (40), if $N=1$, any phase shift satisfying (39b) or (39c)
is optimal.

Proof: Recalling (22), if $N=1$ we have $x=y=0$ and
${f_{c}}\left({\bf\Phi}\right)=e^{j{\theta_{1}}}$. Therefore, any $\theta_{1}$
will have the same results of $\left|{f_{c}}\left({\bf\Phi}\right)\right|=1$
and $f_{k}^{H}\left({\bf\Phi}\right)f_{i}\left({\bf\Phi}\right)=1$.
$\blacksquare$

######  Proposition 2.

For problem (39) and (40), if $\delta=0$ or $\varepsilon_{k}=0,\forall k$, any
phase shift satisfying (39b) or (39c) is optimal.

Proof: Recalling (23), if $\delta=0$ or $\varepsilon_{k}=0,\forall k$, all
terms related to $\bf\Phi$ become zero. $\blacksquare$

This result indicates that if the environment between the BS and the RIS or
that between the RIS and all users has rich scatters, there is no need to
design the phase shifts of RIS based on the statistical CSI.

######  Proposition 3.

For problem (39), if $p_{i}=0,\forall i\neq k$, aligning $\bf\Phi$ to user $k$
is optimal.

Proof: If $p_{i}=0,\forall i\neq k$, it becomes RIS-aided single-user systems.
In this case, we can directly apply the results in [26]. $\blacksquare$

Proposition 3 also indicates that if a user $k$ locates very close to the RIS,
aligning the phase shifts of RIS to this user will almost yield an optimal sum
rate.

###  IV-B General Case

Next, we consider the optimization problems (39) and (40) in the general case.
Since the expression of the rate has a complicated form and the active and
passive beamforming are closely coupled, it is hard to obtain a globally
optimal solution in general. Therefore, we propose a GA-based method to solve
the two optimization problems.

GA simulates the evolution of a population in the nature[37, 38], and its main
steps are summarized in Fig. 2. Next, we describe the implementation details
of the GA-based optimization method for both the problem (39) and (40).

Figure 2: The outline of GA-based method.

####  IV-B1 Initial population

GA is initialized by generating a population with $N_{t}$ individuals. Each
individual contains $N$ chromosomes, and the $n$-th chromosome corresponds to
RIS’s phase shift $\theta_{n}$. For continuous phase shifts, we randomly
generate the initial chromosomes of individuals in the population in
$\left[0,2\pi\right)$. For discrete phase shifts, we randomly generate the
initial chromosomes of individuals from the set in (39c). Then, we tend to
evolve this initial population to the next generation following the steps
below.

####  IV-B2 Fitness evaluation and scaling

We first evaluate the fitness of each individual in the current population.
The fitness evaluation function is the objective function in (39a) or (40a),
respectively. This operation means that an individual with better fitness in
the population corresponds to a better solution for optimization problems (39)
or (40). Next, we need to scale the raw fitness value of individuals based on
their rank in the population. We sort the raw fitness of individuals and
compute their scaled fitness as follows

$\displaystyle\begin{aligned}
&f_{i}=\frac{1}{\sqrt{\operatorname{rank}_{i}}},\operatorname{rank}_{i}\in\left[1,\ldots,i,\ldots,N_{t}\right],\\\
&f_{i}^{\text{scaled}}=2N_{c}\frac{f_{i}}{\sum_{i=1}^{N_{\text{t
}}}f_{i}},\end{aligned}$ (41)

where $\operatorname{rank}_{i}$ is the index of raw fitness of individual $i$
after descending sort, $f_{i}^{\text{scaled}}$ is the scaled fitness of
individual $i$, $N_{c}$ is a parameter used in the selection operation. This
scaling operation can restrict the effect of individuals with large fitness
which may reproduce their chromosomes too frequently and cause prematurity.
After the adjustment of raw fitness values, raw fitness values are converted
to a more suitable range and we can carry out the selection operation better.

####  IV-B3 Selection

Here we will select some individuals from current population, and some of them
are selected as elites, some of them are chosen as parents which could
generate offspring. First, $N_{e}$ individuals with larger
$f_{i}^{\text{scaled}}$ are selected as elites, and they will be directly
passed to the next generation. Then we will select $2N_{c}$ parents based on
stochastic universal sampling, which has a fast and accurate selection
mechanism. To perform stochastic universal sampling, we first form a roulette
wheel which has $2N_{c}$ slots and the size of slot $i$ is proportional to
$f_{i}^{\text{scaled}}$ as follows

$\displaystyle\operatorname{slot}_{i}=\frac{f_{i}^{\text{scaled}}}{2N_{c}},$
(42)

where we have $\sum_{i=1}^{N_{t}}\operatorname{slot}_{i}=1$. Therefore, each
slot corresponds to an individual. Then we rotate the roulette wheel $2N_{c}$
times, each time forwarding with an equal step $\frac{1}{2N_{c}}$. After each
time rotation, we find where the wheel pointer falls and select the
corresponding individual as a parent. After $2N_{c}$ times rotation, we can
select $2N_{c}$ parents which will be used for crossover operation. Note that
one individual may appear multiple times in this $2N_{c}$ combination, and its
appearance probability is proportional to its scaled fitness value. Finally,
the remaining $N_{m}=N_{t}-N_{e}-N_{c}$ individuals will be used for mutation
operation.

####  IV-B4 Crossover

We will use previously selected $2N_{c}$ parents to perform crossover and
generate $N_{c}$ offspring. Crossover operation can extract the best
chromosome from different parents and recombine them into potentially superior
offspring. When $N\leq 2$, we will use the single point crossover method.
Otherwise, two points crossover method is adopted in this paper. The pseudo
codes of crossover operation are shown in Algorithm 1.

Algorithm 1 Crossover Algorithm

1: Set $c_{1}=1$, $c_{2}=2$;

2: if $N>2$ then

3: for $i=1:N_{c}$ do

4: Select the $c_{1}$-th and the $c_{2}$-th parents in the $2N_{c}$
combination;

5: Generate different integers $i_{1}$ and $i_{2}$ randomly from
$\left[1,N-1\right]$;

6: if $i_{1}>i_{2}$ then

7: Swap $i_{1}$ and $i_{2}$;

8: Swap parents $c_{1}$ and $c_{2}$;

9: end if

10: Generate the $i$-th offspring by $\left[\text{parent
}c_{1}(1:i_{1}),\text{parent }c_{2}(i_{1}+1,i_{2}),\text{parent
}c_{1}(i_{2}+1,N)\right]$;

11: $c_{1}=c_{1}+2$, $c_{2}=c_{2}+2$;

12: end for

13: else

14: Generate a random integer $i_{1}$ and perform single point crossover;

15: end if

####  IV-B5 Mutation

$N_{m}$ parents will experience mutation operation with probability $p_{m}$
and produce $N_{m}$ offspring. Mutation operation can increase the diversity
of the population and bring the likelihood that offspring with better fitness
is generated. We use the uniform mutation method and its pseudo codes are
shown in Algorithm 2.

Algorithm 2 Mutation Algorithm

1: for $i=1:N_{m}$ do

2: for $n=1:N$ do

3: if $\text{rand}\left(1\right)<p_{m}$ then

4: if RIS has continous phase shifts then

5: the $n$-th chromosome $\theta_{n}$ of parent $i$ mutates to
$2\pi\times\text{rand}(1)$;

6: else

7: the $n$-th chromosome ${\theta}_{n}$ of parent $i$ mutates to a value
randomly selected from the set in (39c);

8: end if

9: end if

10: end for

11: end for

After the above operations, we combine the $N_{e}$ elite children, $N_{c}$
children coming from crossover and $N_{m}$ children coming from mutation to
form the next generation population. The GA will stop if the number of
generations is larger than $N_{max}$ or the change of the average fitness
value is less than $\xi$.

##  V Numerical Results

In this section, we validate our analysis and demonstrate the benefits brought
by deploying RIS into massive MIMO systems. Our simulation parameters are set
as in [5, 12]. We assume the locations of the BS and the RIS are $(0,0,25)$
and $(5,100,30)$, respectively. Similar to [5], we assume users are located on
a half circle centered at $(5,100)$ with radius of $5$m and height of $1.6$m.
The AoA and AoD of BS, RIS and users are generated randomly from $[0,2\pi]$[9,
12] and these angles will be fixed after initial generation. Unless otherwise
stated, our simulation parameters are set as follows: element spacing of
$d=\frac{\lambda}{2}$, number of users of $K=4$, number of reflecting elements
of $N=16$, number of antennas of $M=64$, transmit power of $p_{k}=30$ dBm,
$\forall k$, noise power of $\sigma^{2}=-104$ dBm and Rician factor of
$\delta=1$, $\varepsilon_{k}=10,\forall k$. Large-scale path-loss is
calculated as $\alpha_{k}=\frac{1}{1000{d_{k}}^{\alpha_{kUR}}},\forall k$ and
$\beta=\frac{1}{1000d_{0}^{\beta_{RB}}}$ [12] where $d_{k}$ and $d_{0}$ are
respectively the distances of user $k$-RIS and RIS-BS, and the path-loss
exponents are $\alpha_{kUR}=\beta_{RB}=2.8,\forall k$ [5]. GA population
$N_{t}$ is 200, elite number is $N_{e}=10$, number of crossover parents is
$N_{c}=152$, number of mutation parents is $N_{m}=38$. The following
simulation results are obtained by averaging over $10000$ random channel
generations.

Figure 3: Desired signal power and sum interference power of user $1$ under
random RIS phase shifts.

To begin with, we validate the correctness of our key derivation in Lemma 1.
In Fig. 3, we show the desired signal
$\mathbb{E}\left\\{\left\|\mathbf{g}_{1}\right\|^{4}\right\\}$ and sum inter-
user interference
$\sum_{i=2}^{4}\mathbb{E}\left\\{\left|\mathbf{g}_{1}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$
for user $1$ under two independent random realizations of $\mathbf{\Phi}$.
Fig. 3 shows that our derived expressions perfectly match the Monte Carlo
simulation, which verifies the accuracy of our results.

Next, we evaluate the impact of various system parameters on the data rate of
the RIS-aided massive MIMO system. To this end, two kinds of optimal phase
shifts $\mathbf{\Phi^{*}}$ are obtained by respectively solving optimization
problem (39) and problem (40), and the obtained $\mathbf{\Phi^{*}}$ will be
used to calculate two different performance metrics, i.e., the sum user rate
$\sum_{k=1}^{K}R_{k}\left(\mathbf{\Phi}^{*}\right)$ and the minimum user rate
$\min\limits_{k}R_{k}\left(\mathbf{\Phi^{*}}\right)$. We refer to the sum user
rate calculated by $\mathbf{\Phi}^{*}$ obtained from problem (39) as “sum rate
by max-sum”, refer to the minimum user rate calculated by $\mathbf{\Phi}^{*}$
obtained from problem (39) as “min rate by max-sum”, refer to the sum rate
calculated by $\mathbf{\Phi}^{*}$ obtained from problem (40) as “sum rate by
max-min” and refer to the minimum user rate calculated by $\mathbf{\Phi}^{*}$
obtained from problem (40) as “min rate by max-min”, respectively. Besides, we
will calculate the sum rate and minimum user rate under random RIS phase
shifts setting by averaging over 1000 random phase shifts generations.

###  V-A Trade-off between path-loss and spatial multiplexing

Figure 4: Sum rate and minimum user rate vesus the Rician factor of RIS-BS
channel.

Fig. 4 shows the four different kinds of rate versus the Rician factor
$\delta$ of RIS-BS channel $\mathbf{H}_{2}$. Results show that our
approximated analytical rate expression (23) matches well with the simulation
result, which verifies the correctness of the derived results. We can see that
when $\delta$ is small, both the sum rate maximization (39) and minimum rate
maximization (40) lead to a similarly good performance (both in terms of sum
rate and minimum rate). This means that in the rich scattering environment,
one can simultaneously achieve a large system capacity while guaranteeing user
fairness. However, when $\delta$ increases, it becomes impossible to balance
the system capacity and fairness. If we maximize the sum rate, the minimum
user rate will approach zero. Conversely, if we want to maintain the minimum
rate, the sum rate will be severely degraded, which nearly equals the rate
achieved by random phases. This result is totally different from the RIS-aided
single-user system with statistical CSI[26, 27], whose rate performance will
be improved by increasing the Rician factor. The reason lies in that when
$\delta$ increases, channel LoS components will become more dominant, which
increases the channel correlation between different users, as well as
increases the inter-user interference and reduces the spatial multiplexing
gain. Specifically, when $\delta\rightarrow\infty$, the rank of the cascaded
channel $\mathbf{G}$ will approach one, and the system will be incapable of
supporting the communication of multiple users.

According to the above discussion, we know it is better to deploy the RIS in
the environment with relatively rich scatters to support multi-user
communications. However, to ensure the rich scatters, the distance between the
BS and RIS should be increased, yielding an increased path loss and a larger
path-loss exponent. Therefore, we present Fig. 5 to show the impacts of RIS-BS
channel path-loss exponent $\beta_{RB}$.

Figure 5: Sum rate and minimum rate versus the path-loss exponent $\beta_{RB}$
of RIS-BS channel.

Firstly, we can see that when $\beta_{RB}$ is small, both the max-sum problem
(39) and max-min problem (40) can achieve similarly good performance.
Secondly, we can see that as $\beta_{RB}$ keeps increasing, if we want to
maintain fairness, the rate performance (sum-rate and min-rate) will decrease
and eventually approach the rate achieved by random phases. These observations
indicate that if we want to simultaneously achieve high system throughput and
guarantee fairness, the path-loss exponent should be as small as possible,
which corresponds to short distances and high value of Rician factors.
Therefore, there exists a trade-off between the achievable spatial
multiplexing gain and the unwanted channel path-loss.

###  V-B The interplay between RIS and massive MIMO

In this subsection, we aim to answer the question about what benefits are
brought by deploying RIS in massive MIMO systems. Note that to guarantee
fairness, only the minimum user rate maximization (40) is considered in the
following simulation.

Figure 6: Average condition numbers of the cascaded channel $\mathbf{G}$.

Fig. 6 shows the standard condition number (i.e., the ratio between the
largest to the smallest eigenvalue [39]) of the cascaded channel $\mathbf{G}$
versus the number of RIS elements $N$, and the result is obtained from Monte
Carlo simulation. It is well known that channel matrix with a lower condition
number can achieve better performance in the high signal-to-noise ratio (SNR)
regime[40], and the channel matrix with condition number $1$ is referred to as
“well-conditioned”. Fig. 6 shows that the condition number of the cascaded
channel decreases quickly as $N$ increases. Besides, after the optimization of
RIS’s phase shifts, we can see that the channel will become nearly well-
conditioned. This finding indicates that RIS can reshape the channel in
massive MIMO systems, reduce the disparity among the channel singular values
and achieve a higher capacity in the high SNR regime.

Figure 7: Sum rate and minimum user rate versus the number of BS antennas $M$.

Fig. 7 shows the data rate performance of RIS-aided massive MIMO systems with
the simple MRC technique. We can see that although the inter-user interference
makes the minimum rate and sum rate approach saturation when
$M\rightarrow\infty$, it still has some promising features. Firstly, by
increasing the number of RIS elements $N$, the data rate can be significantly
improved, which demonstrates the benefits of integrating RIS into massive MIMO
networks. By contrast, in the conventional massive MIMO networks without RIS,
the number of the antennas should be extremely large to serve excessive number
of users. However, the increase of the number of active antennas requires a
large-sized array, high power consumption and high hardware cost. By observing
Fig. 7, we can find that thanks to RIS’s passive beamforming gain, only a
moderate number of antennas are enough to bring promising throughput. For
example, $100$ antennas with $64$ RIS elements can outperform $400$ antennas
with $16$ RIS elements. Therefore, RIS-aided massive MIMO systems are
promising to be applied in future communication systems with much reduced
hardware cost and power consumption, while still maintaining the network
capacity requirement.

Figure 8: Achievable rate versus the number of served users $K$.

Fig. 8 further examines the capability of supporting multiple users in RIS-
aided massive MIMO systems. Here we increase the number of users which are
located on the same circle centered at the IRS with a radius of $5$m, and six
users’ angles are randomly generated. We can see that the minimum user rate
decreases with the increase of the number of users, but the sum rate increases
with $K$. This result is reasonable since we consider the minimum user rate
maximization problem. Although the minimum user rate decreases by increasing
the number of users served simultaneously, it can be significantly promoted by
increasing the number of RIS’s elements and carefully designing RIS’s phase
shifts.

In Fig. 9, we examine the power scaling laws in the RIS-aided massive MIMO
systems, where the transmit power of each user is scaled down as
$p_{k}=100/M,\forall k$.

Figure 9: Minimum user rate versus the number of BS antennas, with scaled
transmission power $p_{k}=100/M,\forall k$.

It has been proved that the massive MIMO technique can help users decrease
their uplink transmitting power while maintaining the data rate performance
[41]. Besides, in the RIS-aided massive MIMO systems, the transmit power of
each user can be further cut down by carefully designing the phase shifts of
RIS relying on statistical CSI. Meanwhile, the increase of RIS’s size also has
a positive impact on saving power consumption.

###  V-C The impacts brought by RIS limited precision

Finally, in Fig. 10, we assess the performance degradation brought by RIS’s
discrete phase shifts in massive MIMO systems.

Figure 10: Minimum user rate under continuous and discrete phase shifts.

Firstly, we can see that both random continuous phase shifts and random
discrete phase shifts lead to the same rate performance, which is consistent
with our derivation in Corollary 4. Secondly, we can see that in the RIS-aided
massive MIMO systems, the degradation due to low-resolution reflecting
elements is marginal which does not enlarge when $N$ increases. Hence, it will
not be an implementation bottleneck in practical systems. Meanwhile, the
degradation can be easily compensated by increasing $N$, and the degradation
does not enlarge when increasing the number of antennas at BS. We conjecture
that the reason for this phenomenon lies in that the robustness of data rate
against the low-resolution of individual reflecting elements is increased by
means of large $N$ and large $M$. Therefore, it is feasible to deploy RIS with
low hardware cost but large size in massive MIMO systems.

##  VI Conclusion

This paper has analyzed and optimized the uplink achievable rate performance
in the uplink RIS-aided massive MIMO systems with statistical CSI. We have
applied the RIS to provide an additional communication link to the user in the
dead zone of conventional massive MIMO systems. We have designed the phase
shifts of the RIS based on statistical CSI, which could reduce the
implementation complexity and the signaling overhead. To this end, first, we
have derived the closed-form expressions for the uplink achievable rate which
hold for any finite numbers of BS antennas. We have then investigated the
power scaling laws, analyzed the rate under some special cases and presented
the average asymptotic rate achieved by the random phase shift setting. Then,
we have studied the optimal phase shifts in some special cases and used the
GA-based method to solve the sum-rate maximization and the minimum user rate
maximization problems in the general case. Finally, we have provided the
numerical results to validate the potential of integrating RIS into existing
massive MIMO systems. Our results have revealed the trade-off between the
achievable spatial multiplexing gain and unwanted path-loss. Besides, we have
demonstrated that it is promising to use RIS with low-resolution hardware to
enhance the coverage in massive MIMO systems.

## Appendix A

To begin with, we present some definition and properties which will be
utilized in the following derivation.

According to the definition of Rician channels in (1) and (2), we can rewrite
the cascaded channels ${\bf g}_{k}$ for user $k$ and ${\bf g}_{i}$ for user
$i$ as follows

$\displaystyle\mathbf{g}_{k}=\\!\mathbf{H}_{2}{\bf\Phi}\mathbf{h}_{k}=\\!\sqrt{\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}}(\underbrace{\sqrt{\delta\varepsilon_{k}}\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{k}}_{\mathbf{g}_{k}^{1}}+\underbrace{\sqrt{\delta}\overline{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{k}}_{\mathbf{g}_{k}^{2}}+\underbrace{\sqrt{\varepsilon_{k}}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{k}}_{\mathbf{g}_{k}^{3}}+\underbrace{\tilde{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{k}}_{\mathbf{g}_{k}^{4}}),$
(43)
$\displaystyle\mathbf{g}_{i}=\mathbf{H}_{2}{\bf\Phi}\mathbf{h}_{i}=\sqrt{\frac{\beta\alpha_{i}}{(\delta+1)\left(\varepsilon_{i}+1\right)}}(\underbrace{\sqrt{\delta\varepsilon_{i}}\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}}_{\mathbf{g}_{i}^{1}}+\underbrace{\sqrt{\delta}\overline{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{i}}_{\mathbf{g}_{i}^{2}}+\underbrace{\sqrt{\varepsilon_{i}}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}}_{\mathbf{g}_{i}^{3}}+\underbrace{\tilde{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{i}}_{\mathbf{g}_{i}^{4}}).$
(44)

Note that $\tilde{\mathbf{H}}_{2}$, $\tilde{\mathbf{h}}_{k}$ and
$\tilde{\mathbf{h}}_{i}$ are independent with each other, and
$\tilde{\mathbf{H}}_{2}$, $\tilde{\mathbf{h}}_{k}$ and
$\tilde{\mathbf{h}}_{i}$ are composed of independent and identically
distributed random variables following $\mathcal{CN}\left(0,1\right)$.
Therefore, for arbitrary $m$ and $n$, we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right\\}=\mathbb{E}\left\\{\tilde{\mathrm{h}}_{km}\right\\}=\mathbb{E}\left\\{\tilde{\mathrm{h}}_{im}\right\\}=0,\\\
\mathbb{E}\left\\{\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\tilde{\mathrm{h}}_{km}\tilde{\mathrm{h}}_{im}\right\\}=\mathbb{E}\left\\{\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right\\}\mathbb{E}\left\\{\tilde{\mathrm{h}}_{km}\right\\}\mathbb{E}\left\\{\tilde{\mathrm{h}}_{im}\right\\}=0,\\\
\mathbb{E}\left\\{\tilde{\mathrm{h}}_{km}\tilde{\mathrm{h}}^{*}_{kn}\right\\}=\mathbb{E}\left\\{\tilde{\mathrm{h}}_{km}\right\\}\mathbb{E}\left\\{\tilde{\mathrm{h}}^{*}_{kn}\right\\}=0,\forall
m\neq n\\\
\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{km}\right|^{2}\right\\}=\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{im}\right|^{2}\right\\}=\mathbb{E}\left\\{\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\right\\}=1,\end{array}$
(49)

where $\left[{\mathbf{H}}\right]_{mn}$ denotes the $(m,n)$-th entry of matrix
${\mathbf{H}}$ and $\left[\mathbf{h}_{c}\right]_{m}\triangleq{\rm{h}}_{cm}$
represents the $m$-th element of column vector ${\bf h}_{c}$.

Next, we will derive
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}$,
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$ and
$\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$,
respectively.

###  A-A Derivation of
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}$

Using the definition in (43),
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}$ can be written
as

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}=\mathbb{E}\left\\{\mathbf{g}^{H}_{k}\mathbf{g}_{k}\right\\}=\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\mathbb{E}\left\\{\sum_{\omega=1}^{4}\left(\mathbf{g}_{k}^{\omega}\right)^{H}\sum_{\psi=1}^{4}\mathbf{g}_{k}^{\psi}\right\\}.$
(50)

Based on (49), we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left(\mathbf{g}_{k}^{\omega}\right)^{H}\mathbf{g}_{k}^{\psi}\right\\}=0,\forall\omega\neq\psi.\end{array}$
(52)

Therefore, we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\mathbf{g}_{k}^{H}\mathbf{g}_{k}\right\\}\\\
=\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\mathbb{E}\left\\{\sum_{\omega=1}^{4}\left(\mathbf{g}_{k}^{\omega}\right)^{H}\mathbf{g}_{k}^{\omega}\right\\}\\\
=\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\left(\delta\varepsilon_{k}\left\|\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{k}\right\|^{2}+\delta\mathbb{E}\left\\{\left\|\overline{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{k}\right\|^{2}\right\\}+\varepsilon_{k}\mathbb{E}\left\\{\left\|\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{k}\right\|^{2}\right\\}+\mathbb{E}\left\\{\left\|\tilde{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{k}\right\|^{2}\right\\}\right)\\\
{\mathop{=}\limits^{\left(a\right)}}\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\left(\delta\varepsilon_{k}M\left|f_{k}({\bf\Phi})\right|^{2}+\delta
MN+\varepsilon_{k}MN+MN\right)\\\
=M\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\left(\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}+\left(\delta+\varepsilon_{k}+1\right)N\right),\end{array}$
(58)

where $(a)$ utilizes the following results

$\displaystyle\begin{array}[]{l}\left\|\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{k}\right\|^{2}=\left\|\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\right\|^{2}\left\|\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\overline{\mathbf{h}}_{k}\right\|^{2}=M\left|f_{k}({\bf\Phi})\right|^{2},\\\
\mathbb{E}\left\\{\tilde{\mathbf{h}}_{k}\tilde{\mathbf{h}}_{k}^{H}\right\\}=\mathbf{I}_{N},{\bf\Phi}{\bf\Phi}^{H}={\bf
I}_{N},\\\
\mathbb{E}\left\\{\tilde{\mathbf{h}}_{k}^{H}\tilde{\mathbf{h}}_{k}\right\\}=\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{k}=N,\\\
\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}=M\mathbf{I}_{N},\\\
\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}\tilde{\mathbf{H}}_{2}^{H}\right\\}=N\mathbf{I}_{M}.\end{array}$
(64)

###  A-B Derivation of
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$

We can divide $\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$
into the following two parts

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}\\\
=\mathbb{E}\left\\{\left(\sum\limits_{m=1}^{M}\left|{\rm{g}}_{km}\right|^{2}\right)^{2}\right\\}\\\
=\sum\limits_{m=1}^{M}\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{4}\right\\}+2\sum\limits_{m=1}^{M-1}\sum\limits_{h=m+1}^{M}\mathbb{E}\left\\{\left|{{\rm{g}}}_{km}\right|^{2}\left|{\rm{g}}_{kh}\right|^{2}\right\\},\end{array}$
(68)

where ${\rm{g}}_{km}$ is the $m$-th entry of $\mathbf{g}_{k}$.

Next, we will calculate
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{4}\right\\}$ and
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{2}\left|{\rm{g}}_{kh}\right|^{2}\right\\}$,
respectively.

####  A-B1 Calculate
$\mathbb{E}\left\\{\left|{\mathrm{g}}_{km}\right|^{4}\right\\}$

Recalling (1) $\sim$ (4), we can rewrite $\mathrm{g}_{km}$ in the following
form

$\displaystyle\begin{aligned}
&\mathrm{g}_{km}=\sqrt{\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}}\times\left(\underbrace{\sqrt{\delta\varepsilon_{k}}\mathrm{a}_{Mm}\left(\phi_{r}^{a},\phi_{r}^{e}\right)f_{k}({\bf\Phi})}_{\mathrm{g}_{km}^{1}}+\underbrace{\sqrt{\delta}\mathrm{a}_{Mm}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}}_{\mathrm{g}_{km}^{2}}\right.\\\
&\qquad+\left.\underbrace{\sqrt{\varepsilon_{k}}\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{mn}e^{j\theta_{n}}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)}_{{\rm{g}}_{km}^{3}}+\underbrace{\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{mn}e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}}_{\mathrm{g}_{km}^{4}}\right),\end{aligned}$
(69)

where $\mathrm{a}_{Xi}\left({\vartheta^{a},\vartheta^{e}}\right)$ is the
$i$-th element of ${\bf a}_{X}\left({\vartheta^{a},\vartheta^{e}}\right)$.

Therefore, $\mathbb{E}\left\\{\left|{\mathrm{g}}_{km}\right|^{4}\right\\}$ can
be calculated as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathrm{g}_{km}\right|^{4}\right\\}\\\
=\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}+\mathrm{g}_{km}^{2}+\mathrm{g}_{km}^{3}+\mathrm{g}_{km}^{4}\right|^{4}\right\\}\\\
{\mathop{=}\limits^{\left(b\right)}}\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\left(\mathbb{E}\left\\{\sum\limits_{\omega=1}^{4}\left|\mathrm{g}_{km}^{\omega}\right|^{4}\right\\}+2\mathbb{E}\left\\{\sum\limits_{\omega=1}^{3}\sum\limits_{\psi=\omega+1}^{4}\left|\mathrm{g}_{km}^{\omega}\right|^{2}\left|\mathrm{g}_{km}^{\psi}\right|^{2}\right\\}\right.\\\
\qquad\qquad\qquad\qquad+\left.4\mathbb{E}\left\\{\sum\limits_{\omega=1}^{3}\sum\limits_{\psi=\omega+1}^{4}\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{\omega}\right)^{*}\mathrm{g}_{km}^{\psi}\right\\}\right)^{2}\right\\}\right.\\\
\qquad\qquad\qquad\qquad+8\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{2}\right\\}\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{3}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right\\}\\\
\qquad\qquad\qquad\qquad+8\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{3}\right\\}\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{2}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right\\}\Bigg{)}\end{array}$
(76)

where $(b)$ is obtained by removing the zero terms. Since each element in
$\tilde{\mathbf{H}}_{2}$ and $\tilde{\mathbf{h}}_{k}$ is composed of
independent real and imaginary parts following
$\mathcal{N}\left(0,\frac{1}{2}\right)$, we can filter the zero items based on
the property that the $k$-order raw moment
$\mathbb{E}\left\\{s^{k}\right\\}=0$, when $k$ is odd and $s$ is a normal
distribution variable with zero mean[42].

Next, we will calculate the above terms in (76) one by one.

Firstly, we calculate
$\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{\omega}\right|^{4}\right\\},1\leq\omega\leq
4$. When $\omega=1$ we have

$\displaystyle\mathbb{E}\left\\{\left|{\rm{g}}_{km}^{1}\right|^{4}\right\\}=\left|{\rm{g}}_{km}^{1}\right|^{4}=\left(\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}\right)^{2}.$
(77)

When $\omega=2$, we have

$\displaystyle\begin{aligned}
\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{4}\right\\}&=\mathbb{E}\left\\{\left|\sqrt{\delta}\mathrm{a}_{Mm}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}\\\
&=\delta^{2}\mathbb{E}\left\\{\left|\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}=\delta^{2}\mathbb{E}\left\\{\left(\left|\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\right)^{2}\right\\}\\\
&=\delta^{2}\mathbb{E}\left\\{\left(\sum_{n=1}^{N}\left|\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\right.\right.\\\
&\qquad\qquad\left.\left.+2\sum_{n_{1}=1}^{N-1}\sum_{n_{2}=n_{1}+1}^{N}\operatorname{Re}\left\\{\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\tilde{\mathrm{h}}_{kn_{2}}^{*}e^{-j\theta_{n_{2}}}\mathrm{a}_{Nn_{2}}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\right\\}\right)^{2}\right\\}\\\
&{\mathop{=}\limits^{\left(c\right)}}\,\delta^{2}\mathbb{E}\left\\{\left(\sum_{n=1}^{N}\left|\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\right)^{2}\right\\}\\\
&\quad+4\delta^{2}\mathbb{E}\left\\{\left(\sum_{n_{1}=1}^{N-1}\sum_{n_{2}=n_{1}+1}^{N}\operatorname{Re}\left\\{\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\tilde{\mathrm{h}}_{kn_{2}}^{*}e^{-j\theta_{n_{2}}}\mathrm{a}_{Nn_{2}}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\right\\}\right)^{2}\right\\}\\\
&{\mathop{=}\limits^{\left(d\right)}}\delta^{2}\sum_{n=1}^{N}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}\\\
&\quad+2\delta^{2}\sum_{n_{1}=1}^{N-1}\sum_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\right|^{2}\right\\}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn_{2}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{2}}}\tilde{\mathrm{h}}_{kn_{2}}\right|^{2}\right\\}\\\
&\quad+4\delta^{2}\sum_{n_{1}=1}^{N-1}\sum_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\tilde{\mathrm{h}}_{kn_{2}}^{*}e^{-j\theta_{n_{2}}}\mathrm{a}_{Nn_{2}}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\right\\}\right)^{2}\right\\},\end{aligned}$
(78)

where $(c)$ and $(d)$ are obtained by removing the zero expectation terms in
binomial expansion.

Assume $\tilde{\mathrm{h}}_{kn}=s+jt$, where
$s\sim\mathcal{N}\left(0,1/2\right)$ and $t\sim\mathcal{N}\left(0,1/2\right)$.
Then, we have

$\displaystyle\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}=\mathbb{E}\left\\{|s+jt|^{4}\right\\}=\mathbb{E}\left\\{s^{4}+t^{4}+2s^{2}t^{2}\right\\}=2,$
(79)

according to the fact that
$\mathbb{E}\left\\{s^{4}\right\\}=\mathbb{E}\left\\{t^{4}\right\\}=3/4$ and
$\mathbb{E}\left\\{s^{2}\right\\}=\mathbb{E}\left\\{t^{2}\right\\}=1/2$ .

Thus, we can derive

$\displaystyle\begin{array}[]{l}\sum\limits_{n=1}^{N}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}=\sum\limits_{n=1}^{N}\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{kn}\right|^{4}\right\\}=2N.\end{array}$
(81)

Likewise, we have

$\displaystyle\begin{array}[]{l}2\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\right|^{2}\right\\}\mathbb{E}\left\\{\left|\mathrm{a}_{Nn_{2}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{2}}}\tilde{\mathrm{h}}_{kn_{2}}\right|^{2}\right\\}\\\
=2\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{kn_{1}}\right|^{2}\right\\}\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{kn_{2}}\right|^{2}\right\\}\\\
=N(N-1).\end{array}$ (85)

Assume that
$\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}e^{-j\theta_{n_{2}}}\mathrm{a}_{Nn_{2}}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)=\sigma_{n}^{c}+j\sigma_{n}^{s}$,
where $\left(\sigma_{n}^{c}\right)^{2}+\left(\sigma_{n}^{s}\right)^{2}=1$.
Besides, assume that $\tilde{\mathrm{h}}_{kn_{1}}=s_{kn_{1}}+jt_{kn_{1}}$ and
$\tilde{\mathrm{h}}_{kn_{2}}=s_{kn_{2}}+jt_{kn_{2}}$, then we have

$\displaystyle\begin{array}[]{l}4\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\mathrm{a}_{Nn_{1}}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n_{1}}}\tilde{\mathrm{h}}_{kn_{1}}\tilde{\mathrm{h}}_{kn_{2}}^{*}e^{-j\theta_{n_{2}}}\mathrm{a}_{Nn_{2}}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\right\\}\right)^{2}\right\\}\\\
=4\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\mathbb{E}\left\\{\left(\sigma_{n}^{c}s_{kn_{1}}s_{kn_{2}}-\sigma_{n}^{s}t_{kn_{1}}s_{kn_{2}}+\sigma_{n}^{c}t_{kn_{1}}t_{kn_{2}}+\sigma_{n}^{s}s_{kn_{1}}t_{kn_{2}}\right)^{2}\right\\}\\\
=4\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\left(\mathbb{E}\left\\{\left(\sigma_{n}^{c}s_{kn_{1}}s_{kn_{2}}\right)^{2}\right\\}+\mathbb{E}\left\\{\left(\sigma_{n}^{s}t_{kn_{1}}s_{kn_{2}}\right)^{2}\right\\}\right.\\\
\left.\qquad\quad\qquad\qquad+\mathbb{E}\left\\{\left(\sigma_{n}^{c}t_{kn_{1}}t_{kn_{2}}\right)^{2}\right\\}+\mathbb{E}\left\\{\left(\sigma_{n}^{s}s_{kn_{1}}t_{kn_{2}}\right)^{2}\right\\}\right)\\\
=4\sum\limits_{n_{1}=1}^{N-1}\sum\limits_{n_{2}=n_{1}+1}^{N}\left(\left(\sigma_{n}^{c}\right)^{2}+\left(\sigma_{n}^{s}\right)^{2}\right)\frac{1}{4}\times
2\\\ =N(N-1).\end{array}$ (92)

Substituting (81), (85) and (92) into (78), we complete the calculation of
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}^{2}\right|^{4}\right\\}$ as follows

$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{4}\right\\}=\delta^{2}(2N+2N(N-1))=2\delta^{2}N^{2}.$
(93)

When $\omega=3,4$, similarly, we have

$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{4}\right\\}={\varepsilon^{2}_{k}}(2N+2N(N-1))=2{\varepsilon^{2}_{k}}N^{2},$
(94)

and

$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{4}\right\\}=4N+2N(N-1)=2N(N+1).$
(95)

Secondly, we focus on
$\sum\limits_{\omega=1}^{3}\sum\limits_{\psi=\omega+1}^{4}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{\omega}\right|^{2}\left|\mathrm{g}_{km}^{\psi}\right|^{2}\right\\}$.

When $\omega=1$, we can derive

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{km}^{2}\right|^{2}\right\\}=\left|\mathrm{g}_{km}^{1}\right|^{2}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\right\\}=\delta^{2}\varepsilon_{k}\left|f_{k}(\mathbf{\Phi})\right|^{2}N,\\\
\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{km}^{3}\right|^{2}\right\\}=\left|\mathrm{g}_{km}^{1}\right|^{2}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\right\\}=\delta\varepsilon_{k}^{2}\left|f_{k}(\mathbf{\Phi})\right|^{2}N,\\\
\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{km}^{4}\right|^{2}\right\\}=\left|\mathrm{g}_{km}^{1}\right|^{2}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{2}\right\\}=\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}N.\end{array}$
(99)

When $\omega=2$, by utilizing the property of independence and removing the
terms with zero expectation, we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{km}^{3}\right|^{2}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\left|\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\left|\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{mn}e^{j\theta_{n}}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\right|^{2}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\sum_{n=1}^{N}\left|\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}e^{j\theta_{n}}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\right|^{2}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\sum_{n=1}^{N}\left|\tilde{\mathrm{h}}_{kn}\right|^{2}\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\right\\}\\\
=\delta\varepsilon_{k}\sum_{n=1}^{N}\mathbb{E}\left\\{\left|\tilde{\mathrm{h}}_{kn}\right|^{2}\right\\}\sum_{n=1}^{N}\mathbb{E}\left\\{\left.\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\right\\}\\\
=\delta\varepsilon_{k}N^{2},\end{array}$ (106)

and

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{km}^{4}\right|^{2}\right\\}\\\
=\delta\mathbb{E}\left\\{\left|\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\left|\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{mn}e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}\right|^{2}\right\\}\\\
=\delta\mathbb{E}\left\\{\sum_{n=1}^{N}\left|\tilde{\mathrm{h}}_{kn}\right|^{2}\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\left|\tilde{\mathrm{h}}_{kn}\right|^{2}\right\\}\\\
=\delta\mathbb{E}\left\\{\sum_{n_{1}=1}^{N}\sum_{n_{2}=1,n_{2}\neq
n_{1}}^{N}\left|\tilde{\mathrm{h}}_{kn_{1}}\right|^{2}\left|\tilde{\mathrm{h}}_{kn_{2}}\right|^{2}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn_{2}}\right|^{2}+\sum_{n_{1}=1}^{N}\left|\tilde{\mathrm{h}}_{kn_{1}}\right|^{4}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn_{1}}\right|^{2}\right\\}\\\
=\delta N(N+1).\end{array}$ (112)

When $\omega=3$, similarly, we have

$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\left|\mathrm{g}_{km}^{4}\right|^{2}\right\\}=\varepsilon_{k}N(N+1).$
(113)

Thirdly, we calculate
$\sum\limits_{\omega=1}^{3}\sum\limits_{\psi=\omega+1}^{4}\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{\omega}\right)^{*}\mathrm{g}_{km}^{\psi}\right\\}\right)^{2}\right\\}$.
Using the similar methods in (92), we can extract the real parts and then
calculate the expectation of their square. Then we can obtain the following
results after some straightforward simplifications:

$\displaystyle\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{2}\right\\}\right)^{2}\right\\}=\frac{\delta^{2}\varepsilon_{k}}{2}N\left|f_{k}({\bf\Phi})\right|^{2},\quad\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{3}\right\\}\right)^{2}\right\\}=\frac{\delta{\varepsilon_{k}^{2}}}{2}N\left|f_{k}({\bf\Phi})\right|^{2},$
$\displaystyle\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right)^{2}\right\\}=\frac{\delta{\varepsilon_{k}}}{2}N\left|f_{k}({\bf\Phi})\right|^{2},\;\,\quad\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{2}\right)^{*}\mathrm{g}_{km}^{3}\right\\}\right)^{2}\right\\}=\frac{\delta{\varepsilon_{k}}}{2}N^{2},$
$\displaystyle\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{2}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right)^{2}\right\\}={\frac{\delta}{2}N\left(N+1\right)},\quad\quad\;\;\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{3}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right)^{2}\right\\}={\frac{{\varepsilon_{k}}}{2}N\left(N+1\right)}.$
(114)

The remaining two terms in (76) can also be derived similarly, which is
presented as follows

$\displaystyle\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{2}\right\\}\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{3}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right\\}=\frac{1}{2}\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2},$
(115)
$\displaystyle\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{1}\right)^{*}\mathrm{g}_{km}^{3}\right\\}\operatorname{Re}\left\\{\left(\mathrm{g}_{km}^{2}\right)^{*}\mathrm{g}_{km}^{4}\right\\}\right\\}=\frac{1}{2}\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}.$
(116)

Substituting the above intermediate results (77) and (93) $\sim$ (116) into
(76), we complete the calculation of
$\mathbb{E}\left\\{\left|{\mathrm{g}}_{km}\right|^{4}\right\\}$, which is not
related with its subscript $m$. Here we omit its detailed expression since it
is straightforward.

####  A-B2 Calculate
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{2}\left|{\rm{g}}_{kh}\right|^{2}\right\\}$

Similar to (69), we can express $\mathrm{g}_{kh}$ as follows

$\displaystyle\begin{aligned}
&\mathrm{g}_{kh}=\sqrt{\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}}\times\left(\underbrace{\sqrt{\delta\varepsilon_{k}}\mathrm{a}_{Mh}\left(\phi_{r}^{a},\phi_{r}^{e}\right)f_{k}({\bf\Phi})}_{\mathrm{g}_{kh}^{1}}+\underbrace{\sqrt{\delta}\mathrm{a}_{Mh}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\sum_{n=1}^{N}\mathrm{a}_{Nn}^{*}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}}_{\mathrm{g}_{kh}^{2}}\right.\\\
&\qquad+\left.\underbrace{\sqrt{\varepsilon_{k}}\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{hn}e^{j\theta_{n}}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)}_{{\rm{g}}_{kh}^{3}}+\underbrace{\sum_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{hn}e^{j\theta_{n}}\tilde{\mathrm{h}}_{kn}}_{\mathrm{g}_{kh}^{4}}\right),\end{aligned}$
(117)

Note that $\left[\tilde{\mathbf{H}}_{2}\right]_{mn}$ is independent to
$\left[\tilde{\mathbf{H}}_{2}\right]_{hn}$ and both of them have zero mean. We
can extract the terms with non-zero expectation after the binomial expansion
as follows

$\displaystyle\begin{aligned}
\mathbb{E}\left\\{\left|\mathrm{g}_{km}\right|^{2}\left|\mathrm{g}_{kh}\right|^{2}\right\\}&=\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\mathbb{E}\left\\{\left|\sum\limits_{\omega=1}^{4}\mathrm{g}_{km}^{\omega}\right|^{2}\left|\sum\limits_{\psi=1}^{4}\mathrm{g}_{kh}^{\psi}\right|^{2}\right\\}\\\
&=\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\times\left(\sum\limits_{\omega=1}^{4}\sum\limits_{\psi=1}^{4}\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{\omega}\right|^{2}\left|\mathrm{g}_{kh}^{\psi}\right|^{2}\right\\}\right.\\\
&\quad+4\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{1}\left(\mathrm{g}_{kh}^{2}\right)^{*}\right\\}\right\\}\\\
&\quad+4\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\\\
&\quad+4\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{3}\left(\mathrm{g}_{km}^{4}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{1}\left(\mathrm{g}_{kh}^{2}\right)^{*}\right\\}\right\\}\\\
&\quad+\left.4\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{3}\left(\mathrm{g}_{km}^{4}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\right).\end{aligned}$
(118)

Next, we will calculate the above terms in (118) one by one.

First, we focus on
$\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{\omega}\right|^{2}\left|\mathrm{g}_{kh}^{\psi}\right|^{2}\right\\},1\leq\omega,\psi\leq
4$. These terms can be derived following the similar process in the
calculation of
$\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{\omega}\right|^{2}\left|\mathrm{g}_{km}^{\psi}\right|^{2}\right\\}$.
Therefore, we can directly obtain the following results

$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{kh}^{1}\right|^{2}\right\\}=\left(\delta\varepsilon_{k}\right)^{2}\left|f_{k}(\mathbf{\Phi})\right|^{4},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{kh}^{2}\right|^{2}\right\\}=\delta^{2}\varepsilon_{k}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
(119)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{kh}^{3}\right|^{2}\right\\}=\delta\varepsilon_{k}^{2}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{1}\right|^{2}\left|\mathrm{g}_{kh}^{4}\right|^{2}\right\\}=\delta\varepsilon_{k}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
(120)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{kh}^{1}\right|^{2}\right\\}=\delta^{2}\varepsilon_{k}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{kh}^{2}\right|^{2}\right\\}={2\delta^{2}N^{2}},$
(121)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{kh}^{3}\right|^{2}\right\\}=\delta\varepsilon_{k}N^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{2}\right|^{2}\left|\mathrm{g}_{kh}^{4}\right|^{2}\right\\}={\delta(N^{2}+N)},$
(122)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\left|\mathrm{g}_{kh}^{1}\right|^{2}\right\\}=\delta\varepsilon_{k}^{2}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\left|\mathrm{g}_{kh}^{2}\right|^{2}\right\\}=\delta\varepsilon_{k}N^{2},$
(123)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\left|\mathrm{g}_{kh}^{3}\right|^{2}\right\\}=\varepsilon_{k}^{2}N^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{3}\right|^{2}\left|\mathrm{g}_{kh}^{4}\right|^{2}\right\\}=\varepsilon_{k}N^{2},$
(124)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{2}\left|\mathrm{g}_{kh}^{1}\right|^{2}\right\\}=\delta\varepsilon_{k}N\left|f_{k}(\mathbf{\Phi})\right|^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{2}\left|\mathrm{g}_{kh}^{2}\right|^{2}\right\\}={\delta(N^{2}+N)},$
(125)
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{2}\left|\mathrm{g}_{kh}^{3}\right|^{2}\right\\}=\varepsilon_{k}N^{2},$
$\displaystyle\mathbb{E}\left\\{\left|\mathrm{g}_{km}^{4}\right|^{2}\left|\mathrm{g}_{kh}^{4}\right|^{2}\right\\}{=N^{2}+N}.$
(126)

Next, we will derive the remaining four parts in (118). To begin with, the
first one is

$\displaystyle\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{1}\left(\mathrm{g}_{kh}^{2}\right)^{*}\right\\}\right\\}=\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\right)^{2}\right\\}=\frac{\delta^{2}\varepsilon_{k}}{2}N\left|f_{k}({\bf\Phi})\right|^{2}.$
(127)

The second one is

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\operatorname{Re}\left\\{f_{k}({\bf\Phi})\sum\limits_{n=1}^{N}\mathrm{a}_{Nn}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{-j\theta_{n}}\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right.\\\
\left.\qquad\qquad\times\operatorname{Re}\left\\{\sum\limits_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{hn}e^{j\theta_{n}}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\sum\limits_{n=1}^{N}\left[\tilde{\mathbf{H}}_{2}\right]_{hn}^{*}e^{-j\theta_{n}}\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\operatorname{Re}\left\\{\sum\limits_{n=1}^{N}\left(f_{k}({\bf\Phi})\mathrm{a}_{Nn}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{-j\theta_{n}}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right.\\\
\left.\qquad\qquad\times\operatorname{Re}\left\\{\sum\limits_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{hn}\right|^{2}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right\\}.\end{array}$
(133)

Assume that

$\displaystyle\begin{array}[]{l}f_{k}({\bf\Phi})\mathrm{a}_{Nn}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{-j\theta_{n}}=\sigma_{c}^{tn}+j\sigma_{s}^{tn},\\\
\tilde{\mathrm{h}}_{kn}=s^{n}+jt^{n},\\\
\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)=\sigma_{c}^{krn}+j\sigma_{s}^{krn},\end{array}$
(137)

and after some algebraic simplifications, we can obtain

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\\\
=\delta\varepsilon_{k}\mathbb{E}\left\\{\sum\limits_{n=1}^{N}\sigma_{c}^{tn}\sigma_{c}^{krn}\left(s^{n}\right)^{2}+\sigma_{s}^{tn}\sigma_{s}^{krn}\left(t^{n}\right)^{2}\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\sum\limits_{n=1}^{N}\left(\sigma_{c}^{tn}\sigma_{c}^{krn}+\sigma_{s}^{tn}\sigma_{s}^{krn}\right)=\frac{\delta\varepsilon_{k}}{2}\sum\limits_{n=1}^{N}\operatorname{Re}\left\\{\left(\sigma_{c}^{tn}+j\sigma_{s}^{tn}\right)\left(\sigma_{c}^{krn}-j\sigma_{s}^{krn}\right)\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\sum\limits_{n=1}^{N}\operatorname{Re}\left\\{f_{k}({\bf\Phi})\mathrm{a}_{Nn}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)e^{-j\theta_{n}}\mathrm{a}_{Nn}^{*}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\operatorname{Re}\left\\{f_{k}(\mathbf{\Phi})\mathrm{a}_{N}^{H}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\mathbf{\Phi}^{H}\mathrm{a}_{N}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\operatorname{Re}\left\\{f_{k}(\mathbf{\Phi})f_{k}^{H}(\mathbf{\Phi})\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\left|f_{k}({\bf\Phi})\right|^{2}.\end{array}$
(145)

Then, we can easily find that

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{3}\left(\mathrm{g}_{km}^{4}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{1}\left(\mathrm{g}_{kh}^{2}\right)^{*}\right\\}\right\\}\\\
=\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{1}\left(\mathrm{g}_{km}^{2}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\\\
=\frac{\delta\varepsilon_{k}}{2}\left|f_{k}({\bf\Phi})\right|^{2}.\end{array}$
(149)

The last one can be derived as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\mathrm{g}_{km}^{3}\left(\mathrm{g}_{km}^{4}\right)^{*}\right\\}\operatorname{Re}\left\\{\mathrm{g}_{kh}^{3}\left(\mathrm{g}_{kh}^{4}\right)^{*}\right\\}\right\\}\\\
=\varepsilon_{k}\mathbb{E}\left\\{\operatorname{Re}\left\\{\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\operatorname{Re}\left\\{\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{hn}\right|^{2}\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right\\}\\\
=\varepsilon_{k}\mathbb{E}\left\\{\left(\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\operatorname{Re}\left\\{\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right)\left(\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{hn}\right|^{2}\operatorname{Re}\left\\{\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right)\right\\}\\\
=\varepsilon_{k}\mathbb{E}\left\\{\sum_{n=1}^{N}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{mn}\right|^{2}\left|\left[\tilde{\mathbf{H}}_{2}\right]_{hn}\right|^{2}\left(\operatorname{Re}\left\\{\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right)^{2}\right\\}\\\
=\varepsilon_{k}\sum_{n=1}^{N}\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{\mathrm{a}_{Nn}\left(\varphi_{kr}^{a},\varphi_{kr}^{e}\right)\tilde{\mathrm{h}}_{kn}^{*}\right\\}\right)^{2}\right\\}\\\
=\frac{1}{2}\varepsilon_{k}N.\end{array}$ (156)

Substituting (119) $\sim$ (127) and (145) $\sim$ (156) into (118), we can
obtain the expression of
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{2}\left|{\rm{g}}_{kh}\right|^{2}\right\\}$,
which is not related with its subscript $m$ and $h$. Since we have obtained
the expressions of
$\mathbb{E}\left\\{\left|{\mathrm{g}}_{km}\right|^{4}\right\\}$ and
$\mathbb{E}\left\\{\left|{\rm{g}}_{km}\right|^{2}\left|{\rm{g}}_{kh}\right|^{2}\right\\}$
, we can directly obtain
$\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}$ by using

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}=M\mathbb{E}\left\\{\left|\mathrm{g}_{km}\right|^{4}\right\\}+M(M-1)\mathbb{E}\left\\{\left|\mathrm{g}_{km}\right|^{2}\left|\mathrm{g}_{kh}\right|^{2}\right\\}.$
(157)

###  A-C Derivation of
$\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$

Before the proof, we first provide an important property as follows

$\displaystyle\mathbb{E}\left\\{\operatorname{Re}\left\\{\tilde{\mathbf{H}}_{2}\mathbf{A}\tilde{\mathbf{H}}_{2}\right\\}\right\\}=\mathbf{0},$
(158)

where $\mathbf{A}\in\mathbb{C}^{N\times M}$ is an arbitrary deterministic
matrix. This conclusion can be readily proved by firstly considering the case
of one dimension and then generalizing it to high dimensions by mathematical
induction.

Note that since the communication of different users goes through the same
RIS-BS channel $\mathbf{H}_{2}$, $\mathbf{g}_{k}$ is no longer independent to
$\mathbf{g}_{i}$, which is different from the scenario without RIS. Recalling
(43) and (44), when calculating
$\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}$,
we can ignore the terms with zero expectation based on (49) and (158), and
then we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}\\\
=\frac{\beta^{2}\alpha_{k}\alpha_{i}}{(\delta+1)^{2}\left(\varepsilon_{k}+1\right)\left(\varepsilon_{i}+1\right)}\mathbb{E}\left\\{\left|\sum\limits_{\omega=1}^{4}\sum\limits_{\psi=1}^{4}\left(\mathbf{g}_{k}^{\omega}\right)^{H}\mathbf{g}_{i}^{\psi}\right|^{2}\right\\}\\\
=\frac{\beta^{2}\alpha_{k}\alpha_{i}}{(\delta+1)^{2}\left(\varepsilon_{k}+1\right)\left(\varepsilon_{i}+1\right)}\times\left(\mathbb{E}\left\\{\sum\limits_{\omega=1}^{4}\sum\limits_{\psi=1}^{4}\left|\left(\mathbf{g}_{k}^{\omega}\right)^{H}\mathbf{g}_{i}^{\psi}\right|^{2}\right\\}\right.\\\
+2\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{1}\left(\mathbf{g}_{i}^{3}\right)^{H}\mathbf{g}_{k}^{3}\right\\}\right\\}+2\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{2}\left(\mathbf{g}_{i}^{4}\right)^{H}\mathbf{g}_{k}^{3}\right\\}\right\\}\\\
\left.+2\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{1}\left(\mathbf{g}_{i}^{3}\right)^{H}\mathbf{g}_{k}^{4}\right\\}\right\\}+2\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{2}\left(\mathbf{g}_{i}^{4}\right)^{H}\mathbf{g}_{k}^{4}\right\\}\right\\}\right).\end{array}$
(164)

Then we will calculate the above terms in (164) one by one.

First, we focus on
$\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{\omega}\right)^{H}\mathbf{g}_{i}^{\psi}\right|^{2}\right\\},1\leq\omega,\psi\leq
4$.

When $\omega=1$, we have

$\displaystyle\begin{aligned}
&\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{1}\right|^{2}\right\\}=\left|\sqrt{\delta\varepsilon_{k}}\sqrt{\delta\varepsilon_{i}}\overline{\mathbf{h}}_{k}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\right|^{2}=\delta^{2}\varepsilon_{k}\varepsilon_{i}M^{2}\left|f_{k}(\mathbf{\Phi})\right|^{2}\left|f_{i}(\mathbf{\Phi})\right|^{2},\\\
&\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{2}\right|^{2}\right\\}\\\
&=\delta^{2}\varepsilon_{k}\mathbb{E}\left\\{\left|\overline{\mathbf{h}}_{k}^{H}\mathbf{\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{i}\right|^{2}\right\\}=\delta^{2}\varepsilon_{k}M^{2}\left|f_{k}({\bf\Phi})\right|^{2}\mathbb{E}\left\\{\left|\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\tilde{\mathbf{h}}_{i}\right|^{2}\right\\}\\\
&=\delta^{2}\varepsilon_{k}M^{2}N\left|f_{k}({\bf\Phi})\right|^{2},\\\
&\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{3}\right|^{2}\right\\}\\\
&=\delta\varepsilon_{k}\varepsilon_{i}\left|f_{k}({\bf\Phi})\right|^{2}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\mathbf{a}_{M}^{H}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\overline{\mathbf{h}}_{i}\\\
&=\delta\varepsilon_{k}\varepsilon_{i}\left|f_{k}({\bf\Phi})\right|^{2}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}M\mathbf{I}_{N}{\bf\Phi}\overline{\mathbf{h}}_{i}\\\
&=\delta\varepsilon_{k}\varepsilon_{i}\left|f_{k}({\bf\Phi})\right|^{2}MN,\end{aligned}$
(165)

and

$\displaystyle\begin{aligned}
&\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{4}\right|^{2}\right\\}\\\
&=\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}\mathbf{a}_{M}^{H}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}{\bf\Phi}\tilde{\mathbf{h}}_{i}\tilde{\mathbf{h}}_{i}^{H}{\bf{\bf\Phi}}^{H}\tilde{\mathbf{H}}_{2}^{H}\right\\}\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\\\
&=\delta\varepsilon_{k}\left|f_{k}({\bf\Phi})\right|^{2}MN.\end{aligned}$
(166)

Similarly, when $\omega=2$, we have

$\displaystyle\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{1}\right|^{2}\right\\}=\delta^{2}\varepsilon_{i}\left|f_{i}({\bf\Phi})\right|^{2}M^{2}N.$
(167)

Next we have

$\displaystyle\begin{aligned}
&\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{2}\right|^{2}\right\\}\\\
&=M^{2}\delta^{2}\mathbb{E}\left\\{\tilde{\mathbf{h}}_{k}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\tilde{\mathbf{h}}_{i}\tilde{\mathbf{h}}_{i}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\tilde{\mathbf{h}}_{k}\right\\}\\\
&{\mathop{=}\limits^{\left(e\right)}}M^{2}\delta^{2}\mathbb{E}\left\\{\tilde{\mathbf{h}}_{k}^{H}\mathbf{\Phi}^{H}\mathbf{a}_{N}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\mathbb{E}\left\\{\tilde{\mathbf{h}}_{i}\tilde{\mathbf{h}}_{i}^{H}\right\\}{\bf\Phi}^{H}\mathbf{a}_{N}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right)\mathbf{a}_{N}^{H}\left(\varphi_{t}^{a},\varphi_{t}^{e}\right){\bf\Phi}\tilde{\mathbf{h}}_{k}\right\\}\\\
&=\delta^{2}M^{2}N^{2},\end{aligned}$ (168)

where $(e)$ is due to the independence between $\tilde{\mathbf{h}}_{k}$ and
$\tilde{\mathbf{h}}_{i}$.

Similarly, we have

$\displaystyle\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{3}\right|^{2}\right\\}=\delta\varepsilon_{i}MN^{2},$
(169)
$\displaystyle\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{4}\right|^{2}\right\\}=\delta
MN^{2}.$ (170)

When $\omega=3$, we can readily obtain the first two terms as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{3}\right)^{H}\mathbf{g}_{i}^{1}\right|^{2}\right\\}=\delta\varepsilon_{i}\varepsilon_{k}\left|f_{i}({\bf\Phi})\right|^{2}MN,\\\
\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{3}\right)^{H}\mathbf{g}_{i}^{2}\right|^{2}\right\\}=\delta\varepsilon_{k}MN^{2}.\end{array}$
(173)

The third term can be derived as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{3}\right)^{H}\mathbf{g}_{i}^{3}\right|^{2}\right\\}=\varepsilon_{k}\varepsilon_{i}\mathbb{E}\left\\{\left|\overline{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\right|^{2}\right\\}\\\
=\varepsilon_{k}\varepsilon_{i}\overline{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\overline{\mathbf{h}}_{k}.\end{array}$
(176)

Assume that
$\tilde{\mathbf{H}}_{2}=\left[\mathbf{J}_{1},\ldots,\mathbf{J}_{i},\ldots,\mathbf{J}_{N}\right]$
and
$\left[{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\right]_{mn}=\alpha_{mn}$,
we can rewrite the $\left(n1,n2\right)$-th entry of
$\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}$
as follows

$\displaystyle\begin{array}[]{l}{\left[\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right]_{n1,n2}}=\sum\limits_{h=1}^{N}\sum\limits_{m=1}^{N}\mathbf{J}_{n1}^{H}\mathbf{J}_{m}\alpha_{mh}\mathbf{J}_{h}^{H}\mathbf{J}_{n2},\end{array}$
(178)

which can be calculated by discussing the values of $h$ and $m$ under
different situations. After some algebraic simplifications, we can obtain the
following results

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left[\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right]_{n1,n2}\right\\}=\alpha_{n1n2}M^{2},\\\
\mathbb{E}\left\\{\left[\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right]_{n1,n1}\right\\}=M(M+N),\\\
\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}=M^{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}+MN\mathbf{I}_{N}.\end{array}$
(182)

Substituting (182) into (176), we have

$\displaystyle\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{3}\right)^{H}\mathbf{g}_{i}^{3}\right|^{2}\right\\}=\varepsilon_{k}\varepsilon_{i}M\left(N^{2}+M\left|\overline{\mathbf{h}}_{k}^{H}\overline{\mathbf{h}}_{i}\right|^{2}\right).$
(183)

Using (182), we know that
$\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\mathbf{I}_{N}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}=M\left(M+N\right)\mathbf{I}_{N}$.
Therefore, we can obtain the fourth term as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{3}\right)^{H}\mathbf{g}_{i}^{4}\right|^{2}\right\\}\\\
=\varepsilon_{k}\overline{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\overline{\mathbf{h}}_{k}\\\
=\varepsilon_{k}MN(M+N).\end{array}$ (187)

When $\omega=4$, similarly, we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{4}\right)^{H}\mathbf{g}_{i}^{1}\right|^{2}\right\\}=\delta\varepsilon_{i}\left|f_{i}(\mathbf{\Phi})\right|^{2}MN,\\\
\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{4}\right)^{H}\mathbf{g}_{i}^{2}\right|^{2}\right\\}=\delta
MN^{2},\\\
\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{4}\right)^{H}\mathbf{g}_{i}^{3}\right|^{2}\right\\}=\varepsilon_{i}MN(M+N),\\\
\mathbb{E}\left\\{\left|\left(\mathbf{g}_{k}^{4}\right)^{H}\mathbf{g}_{i}^{4}\right|^{2}\right\\}=MN(M+N).\end{array}$
(192)

Similar to the above derivation, the remaining four parts in (164) can be
readily derived as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{1}\left(\mathbf{g}_{i}^{3}\right)^{H}\mathbf{g}_{k}^{3}\right\\}\right\\}\\\
=\delta\varepsilon_{k}\varepsilon_{i}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\mathbf{\Phi}^{H}\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\overline{\mathbf{h}}_{k}\right\\}\\\
=\delta\varepsilon_{k}\varepsilon_{i}M^{2}\operatorname{Re}\left\\{f_{k}^{H}({\bf\Phi})f_{i}({\bf\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\},\end{array}$
(196)
$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{1}\right)^{H}\mathbf{g}_{i}^{2}\left(\mathbf{g}_{i}^{4}\right)^{H}\mathbf{g}_{k}^{3}\right\\}\right\\}\\\
=\delta\varepsilon_{k}\operatorname{Re}\left\\{\overline{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\mathbb{E}\left\\{\tilde{\mathbf{h}}_{i}\tilde{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\overline{\mathbf{h}}_{k}\right\\}\\\
=\delta\varepsilon_{k}M^{2}\left|f_{k}({\bf\Phi})\right|^{2},\end{array}$
(200)
$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{1}\left(\mathbf{g}_{i}^{3}\right)^{H}\mathbf{g}_{k}^{4}\right\\}\right\\}\\\
=\delta\varepsilon_{i}\mathbb{E}\left\\{\operatorname{Re}\left\\{\tilde{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\mathbb{E}\left\\{\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\tilde{\mathbf{h}}_{k}\right\\}\right\\}\\\
=\delta\varepsilon_{i}M\operatorname{Re}\left\\{\mathrm{Tr}\left\\{\overline{\mathbf{H}}_{2}{\bf\Phi}\overline{\mathbf{h}}_{i}\overline{\mathbf{h}}_{i}^{H}\mathbb{E}\left\\{{\tilde{\mathbf{h}}_{k}\tilde{\mathbf{h}}_{k}^{H}}\right\\}{\bf\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\right\\}\right\\}\\\
=\delta\varepsilon_{i}M\operatorname{Re}\left\\{\mathrm{Tr}\left\\{\left|f_{i}({\bf\Phi})\right|^{2}\mathbf{a}_{M}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\mathbf{a}_{M}^{H}\left(\phi_{r}^{a},\phi_{r}^{e}\right)\right\\}\right\\}\\\
=\delta\varepsilon_{i}M^{2}\left|f_{i}({\bf\Phi})\right|^{2},\end{array}$
(206)

and

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\mathbf{g}_{k}^{2}\right)^{H}\mathbf{g}_{i}^{2}\left(\mathbf{g}_{i}^{4}\right)^{H}\mathbf{g}_{k}^{4}\right\\}\right\\}\\\
=\delta\mathbb{E}\left\\{\operatorname{Re}\left\\{\tilde{\mathbf{h}}_{k}^{H}{\bf\Phi}^{H}\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\mathbb{E}\left\\{\tilde{\mathbf{h}}_{i}\tilde{\mathbf{h}}_{i}^{H}{\bf\Phi}^{H}\tilde{\mathbf{H}}_{2}^{H}\tilde{\mathbf{H}}_{2}\right\\}{\bf\Phi}\tilde{\mathbf{h}}_{k}\right\\}\right\\}\\\
=\delta
M\operatorname{Re}\left\\{\operatorname{Tr}\left\\{\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}{\bf\Phi}\mathbb{E}\left\\{\tilde{\mathbf{h}_{k}}\tilde{\mathbf{h}}_{k}^{H}\right\\}{\bf\Phi}^{H}\right\\}\right\\}\\\
=\delta
M\operatorname{Re}\left\\{\operatorname{Tr}\left\\{\overline{\mathbf{H}}_{2}^{H}\overline{\mathbf{H}}_{2}\right\\}\right\\}\\\
=\delta M^{2}N.\end{array}$ (212)

Substituting (165) $\sim$ (173) and (183) $\sim$ (212) into (164), we can
complete the proof of Lemma 1 after some trivial simplifications.

## Appendix B

When the phase shifts of RIS are aligned to user $k$, we have
$f_{k}\left(\mathbf{\Phi}\right)=N$ but
$\left|f_{i}\left(\mathbf{\Phi}\right)\right|$ is bounded when
$N\rightarrow\infty$. Therefore, when $N\rightarrow\infty$ and
$M\rightarrow\infty$, we can obtain the order of magnitude as follows

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}=\mathcal{O}\left(M^{2}N^{4}\right),$
(213)
$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{i}\right\|^{4}\right\\}=\mathcal{O}\left(M^{2}N^{2}\right),$
(214)
$\displaystyle\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}=\mathcal{O}\left(M^{2}N^{3}\right),$
(215)
$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}=\mathcal{O}\left(MN^{2}\right).$
(216)

Therefore, when $M\rightarrow\infty$ and $N\rightarrow\infty$, user $k$’s rate
$R_{k}$ can maintain a non-zero value when we cut the transmission power of
user $k$ as $p_{k}=E_{u}/(MN^{2})$ and cut the transmission power of other
users as $p_{i}=E_{u}/(MN),\forall i\neq k$. However, at the same time, since
$\left|f_{i}\left(\mathbf{\Phi}\right)\right|$ is bounded, the rate of user
$i$ will be zero. At this time, when $N\rightarrow\infty$ and
$M\rightarrow\infty$, the dominant terms in rate expression (23) are those
terms which are on the order of $MN^{2}$:

$\displaystyle\begin{array}[]{l}\frac{E_{u}}{MN^{2}}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}\rightarrow
E_{u}\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}\left(\delta\varepsilon_{k}\right)^{2}MN^{2},\\\
\frac{E_{u}}{MN}\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}\rightarrow
E_{u}\frac{\beta^{2}\alpha_{i}\alpha_{k}}{(\delta+1)^{2}\left(\varepsilon_{i}+1\right)\left(\varepsilon_{k}+1\right)}\delta^{2}\varepsilon_{k}MN^{2},\\\
\sigma^{2}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}\rightarrow\sigma^{2}\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\delta\varepsilon_{k}MN^{2}.\end{array}$
(220)

Thus, after some simplification, the rate can be calculated as follows

$\displaystyle R_{k}$
$\displaystyle=\log_{2}\left(1+\frac{\frac{E_{u}}{MN^{2}}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}}{\sum_{i=1,i\neq
k}^{K}\frac{E_{u}}{MN}\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}+\sigma^{2}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{2}\right\\}}\right)$
$\displaystyle\to\log_{2}\left(1+\frac{E_{u}\frac{\varepsilon_{k}}{\left(\varepsilon_{k}+1\right)}}{E_{u}\sum_{i=1,i\neq
k}^{K}\frac{\alpha_{i}}{\left(\varepsilon_{i}+1\right)\alpha_{k}}+\left(1+\frac{1}{\delta}\right)\frac{\sigma^{2}}{\beta\alpha_{k}}}\right),\text{as
}M,N\to\infty.$ (221)

Besides, we can see that (B) is an increasing function with respect to
$\alpha_{k}$, $\beta$ and $\delta$.

## Appendix C

To begin with, we need to provide some necessary preliminary results. Firstly,
for the ideal RIS with continuous phase shifts, we assume that the phase shift
of each reflecting element $\theta_{n}$ is randomly and independently adjusted
at each fading block following the uniform distribution of
$\mathcal{U}\left[0,2\pi\right]$. Then, for $k_{1}$ with arbitrary values, we
have

$\displaystyle\mathbb{E}\left\\{\cos\left(k_{1}+\theta_{n}\right)\right\\}=\frac{1}{2\pi}\int_{0}^{2\pi}\cos\left(k_{1}+\theta_{n}\right)d\theta_{n}=0,$
(222)
$\displaystyle\mathbb{E}\left\\{\cos^{2}\left(k_{1}+\theta_{n}\right)\right\\}=\frac{1}{2}\left(1+\mathbb{E}\left\\{\cos\left(2k_{1}+2\theta_{n}\right)\right\\}\right)=\frac{1}{2}\left(1+\frac{1}{2\pi}\int_{0}^{2\pi}\cos\left(2k_{1}+2\theta_{n}\right)d\theta_{n}\right)=\frac{1}{2}.$
(223)

Next, for the non-ideal RIS with finite $b$ bits discrete phase shifts, we
assume that each $\theta_{n}$ is randomly and independently adjusted from
$\left\\{0,\frac{2\pi}{2^{b}},2\times\frac{2\pi}{2^{b}},\ldots,\left(2^{b}-1\right)\frac{2\pi}{2^{b}}\right\\}$.
When $b>1$, for $k_{1}$ with arbitrary values, we have

$\displaystyle\mathbb{E}\left\\{\cos\left(k_{1}+\theta_{n}\right)\right\\}$
$\displaystyle=\frac{1}{2^{b}}\sum\nolimits_{t=0}^{2^{b}-1}\cos\left(k_{1}+t\frac{2\pi}{2^{b}}\right)$
$\displaystyle=\frac{1}{2^{b}}\sum\nolimits_{t=0}^{2^{(b-1)}-1}\left(\cos\left(k_{1}+t\frac{2\pi}{2^{b}}\right)+\cos\left(k_{1}+\left(t+2^{(b-1)}\right)\frac{2\pi}{2^{b}}\right)\right)$
$\displaystyle\;{\mathop{=}\limits^{\left(f\right)}}\;0,$ (224)

and

$\displaystyle\mathbb{E}\left\\{\cos^{2}\left(k_{1}+\theta_{n}\right)\right\\}=\frac{1}{2}\left(1+\mathbb{E}\left\\{\cos\left(2k_{1}+2\theta_{n}\right)\right\\}\right)=\frac{1}{2}\left(1+\frac{1}{2^{b}}\sum\nolimits_{t=0}^{2^{b}-1}\cos\left(2k_{1}+2t\frac{2\pi}{2^{b}}\right)\right)$
$\displaystyle=\frac{1}{2}\left(1+\frac{1}{2^{b}}\left(\sum\nolimits_{t=0}^{2^{(b-1)}-1}\cos\left(2k_{1}+t\frac{2\pi}{2^{b-1}}\right)+\sum\nolimits_{t=2^{(b-1)}}^{2^{b}-1}\cos\left(2k_{1}+t\frac{2\pi}{2^{b-1}}\right)\right)\right)$
$\displaystyle\;{\mathop{=}\limits^{\left(g\right)}}\frac{1}{2}\left(1+\frac{2}{2^{b}}\sum\nolimits_{t=0}^{2^{(b-1)}-1}\cos\left(2k_{1}+t\frac{2\pi}{2^{b-1}}\right)\right)$
$\displaystyle=\frac{1}{2}\left(1+\left.\mathbb{E}\left\\{\cos\left(2k_{1}+\theta_{n}\right)\right\\}\right|_{b=b-1}\right)$
$\displaystyle=\frac{1}{2},$ (225)

where $(f)$ and $(g)$ come from $\cos(\vartheta+\pi)=-\cos(\vartheta)$ and
$\cos(\vartheta+2\pi)=\cos(\vartheta)$, respectively.

(222) $\sim$ (C) prove that
$\mathbb{E}\left\\{\cos\left(k_{1}+\theta_{n}\right)\right\\}$ and
$\mathbb{E}\left\\{\cos^{2}\left(k_{1}+\theta_{n}\right)\right\\}$ have the
same values for both continuous and discrete phase shifts when $b>1$.
Therefore, we will not distinguish these two cases in the following
derivation. Besides, since the above equations hold for arbitrary $k_{1}$, we
can obtain the following results from trigonometric identities:

$\displaystyle\mathbb{E}\left\\{\sin\left(k_{1}+\theta_{n}\right)\right\\}=\mathbb{E}\left\\{\cos\left(\left(k_{1}-\frac{\pi}{2}\right)+\theta_{n}\right)\right\\}=0,$
(226)
$\displaystyle\mathbb{E}\left\\{\sin^{2}\left(k_{1}+\theta_{n}\right)\right\\}=\mathbb{E}\left\\{1-\cos^{2}\left(k_{1}+\theta_{n}\right)\right\\}=\frac{1}{2}.$
(227)

Then, using the above preliminary results can help us derive the asymptotic
average rate with random phase shifts. Since $\bf\Phi$ is independent with the
channel $\mathbf{H}_{2}$ and $\mathbf{h}_{k}$, rate $R_{k}$ can be calculated
by substituting the terms involving $\bf\Phi$ in (23) with their expectation.
Thus, we need to further calculate
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\right\\}$,
$\mathbb{E}\left\\{\left|f_{i}({\bf\Phi})\right|^{2}\right\\}$,
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{4}\right\\}$,
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\left|f_{i}({\bf\Phi})\right|^{2}\right\\}$
and
$\mathbb{E}\left\\{\operatorname{Re}\left\\{f_{k}^{H}({\bf\Phi})f_{i}({\bf\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\}\right\\}$,
respectively.

Recalling (22), we rewrite $f_{k}({\bf\Phi})$ and $f_{i}({\bf\Phi})$ as
follows

$\displaystyle\begin{array}[]{l}f_{k}({\bf\Phi})=\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{k}+\theta_{n}\right)},\\\
f_{i}({\bf\Phi})=\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{i}+\theta_{n}\right)}.\end{array}$
(230)

Then, using (222) $\sim$ (227) and the independence between
$\theta_{n1},\theta_{n2},\forall n1\neq n2$, we have

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{e^{j\left(k_{1}+\theta_{n1}\right)}\right\\}=\mathbb{E}\left\\{\cos\left(k_{1}+\theta_{n1}\right)\right\\}+j\mathbb{E}\left\\{\sin\left(k_{1}+\theta_{n1}\right)\right\\}=0,\\\
\mathbb{E}\left\\{e^{j\left(k_{1}+\theta_{n1}\right)}e^{j\left(k_{2}+\theta_{n2}\right)}\right\\}=\mathbb{E}\left\\{e^{j\left(k_{1}+\theta_{n1}\right)}\right\\}\mathbb{E}\left\\{e^{j\left(k_{2}+\theta_{n2}\right)}\right\\}=0,\\\
\mathbb{E}\left\\{\operatorname{Re}\left\\{e^{j\left(k_{1}+\theta_{n1}\right)}e^{-j\left(k_{2}+\theta_{n2}\right)}\right\\}\right\\}=\mathbb{E}\left\\{\cos\left(\left(k_{1}+\theta_{n1}\right)-\left(k_{2}+\theta_{n2}\right)\right)\right\\}=0,\\\
\mathbb{E}\left\\{\left(\operatorname{Re}\left\\{e^{j\left(k_{1}+\theta_{n1}\right)}e^{-j\left(k_{2}+\theta_{n2}\right)}\right\\}\right)^{2}\right\\}=\frac{1}{2}\left(1+\mathbb{E}\left\\{\cos\left(2\left(k_{1}+\theta_{n1}\right)-2\left(k_{2}+\theta_{n2}\right)\right)\right\\}\right)=\frac{1}{2}.\end{array}$
(235)

Utilizing (235), we will calculate the expectation of terms involving
$\bf\Phi$ one by one. Firstly, the term
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\right\\}$ can be
calculated as follows

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\right\\}\\\
=\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{k}+\theta_{n}\right)}\sum_{n=1}^{N}e^{-j\left(\zeta_{n}^{k}+\theta_{n}\right)}\\\
=\sum_{n=1}^{N}1+\mathbb{E}\left\\{\sum_{n1=1}^{N}\sum_{n_{2}=1,n_{2}\neq
n1}^{N}e^{-j\left(\zeta_{n1}^{k}+\theta_{n1}\right)}e^{j\left(\zeta_{n2}^{k}+\theta_{n2}\right)}\right\\}\\\
=N,\end{array}$ (240)

Similarly, for the term
$\mathbb{E}\left\\{\left|f_{i}({\bf\Phi})\right|^{2}\right\\}$ , we also have

$\displaystyle\mathbb{E}\left\\{\left|f_{i}({\bf\Phi})\right|^{2}\right\\}=N.$
(241)

Secondly, we focus on the term
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{4}\right\\}$ which can be
derived as:

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{4}\right\\}\\\
=\left|\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{k}+\theta_{n}\right)}\sum_{n=1}^{N}e^{-j\left(\zeta_{n}^{k}+\theta_{n}\right)}\right|^{2}\\\
=\left|N+2\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)}\right\\}\right|^{2}\\\
=N^{2}+4\mathbb{E}\left\\{\left|\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)}\right\\}\right|^{2}\right\\}\\\
=N^{2}+4\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\mathbb{E}\left\\{\left({\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)}\right\\}}\right)^{2}\right\\}\\\
=N^{2}+4\frac{N(N-1)}{2}\frac{1}{2}\\\ =2N^{2}-N.\end{array}$ (249)

Thirdly, the term
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\left|f_{i}({\bf\Phi})\right|^{2}\right\\}$
can be calculated as:

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\left|f_{i}({\bf\Phi})\right|^{2}\right\\}\\\
=\mathbb{E}\left\\{\left|\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{k}+\theta_{n}\right)}\right|^{2}\left|\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{i}+\theta_{n}\right)}\right|^{2}\right\\}\\\
=\mathbb{E}\left\\{\left(N+2\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)}\right\\}\right)\times\right.\\\
\qquad\left.\left(N+2\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{i}+\theta_{n1}-\zeta_{n2}^{i}-\theta_{n2}\right)}\right\\}\right)\right\\}\\\
=N^{2}+4\mathbb{E}\left\\{\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)}\right\\}\operatorname{Re}\left\\{e^{j\left(\zeta_{n1}^{i}+\theta_{n1}-\zeta_{n2}^{i}-\theta_{n2}\right)}\right\\}\right\\}\\\
=N^{2}+4\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\mathbb{E}\left\\{\cos\left(\zeta_{n1}^{k}+\theta_{n1}-\zeta_{n2}^{k}-\theta_{n2}\right)\cos\left(\zeta_{n1}^{i}+\theta_{n1}-\zeta_{n2}^{i}-\theta_{n2}\right)\right\\}\\\
{\mathop{=}\limits^{\left(h\right)}}N^{2}+2\sum_{n1=1}^{N-1}\sum_{n2=n1+1}^{N}\cos\left(\zeta_{n1}^{k}-\zeta_{n2}^{k}-\zeta_{n1}^{i}+\zeta_{n2}^{i}\right),\end{array}$
(257)

where $\left(h\right)$ is obtained by using prosthaphaeresis. Since the second
term in (257) is bounded, we have
$\mathbb{E}\left\\{\left|f_{k}({\bf\Phi})\right|^{2}\left|f_{i}({\bf\Phi})\right|^{2}\right\\}\rightarrow
N^{2}$ when $N\rightarrow\infty$.

The final term
$\mathbb{E}\left\\{\operatorname{Re}\left\\{f_{k}^{H}({\bf\Phi})f_{i}({\bf\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\}\right\\}$
is derived as:

$\displaystyle\begin{array}[]{l}\mathbb{E}\left\\{\operatorname{Re}\left\\{f_{k}^{H}({\bf\Phi})f_{i}({\bf\Phi})\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right\\}\right\\}\\\
=\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right)\sum_{n1=1}^{N}\sum_{n2=1}^{N}e^{-j\left(\zeta_{n1}^{k}+\theta_{n1}\right)}e^{j\left(\zeta_{n2}^{i}+\theta_{n2}\right)}\right\\}\right\\}\\\
=\mathbb{E}\left\\{\operatorname{Re}\left\\{\left(\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right)\sum_{n=1}^{N}e^{-j\left(\zeta_{n}^{k}+\theta_{n}\right)}e^{j\left(\zeta_{n}^{i}+\theta_{n}\right)}\right\\}\right\\}\\\
=\operatorname{Re}\left\\{\left(\overline{\mathbf{h}}_{i}^{H}\overline{\mathbf{h}}_{k}\right)\sum_{n=1}^{N}e^{j\left(\zeta_{n}^{i}-\zeta_{n}^{k}\right)}\right\\},\end{array}$
(262)

which is bounded when $N\rightarrow\infty$.

By substituting (240) $\sim$ (262) into the corresponding terms in rate
expression (23), we can see that when $M\rightarrow\infty$ and
$N\rightarrow\infty$, the dominant terms are those which have the order of
$\mathcal{O}\left(M^{2}N^{2}\right)$. Thus, when $M\rightarrow\infty$ and
$N\rightarrow\infty$, we have

$\displaystyle\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}\rightarrow$
$\displaystyle
M^{2}N^{2}\left(\frac{\beta\alpha_{k}}{(\delta+1)\left(\varepsilon_{k}+1\right)}\right)^{2}$
$\displaystyle\times\left(2\delta^{2}\varepsilon_{k}^{2}+2\delta\varepsilon_{k}\left(2\delta+\varepsilon_{k}+1\right)+\left(2\delta^{2}+\varepsilon_{k}^{2}+2\delta\varepsilon_{k}+2\delta+2\varepsilon_{k}+1\right)\right),$
(263)

and

$\displaystyle\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}\to
M^{2}N^{2}\frac{\beta^{2}\alpha_{i}\alpha_{k}}{(\delta+1)^{2}\left(\varepsilon_{i}+1\right)\left(\varepsilon_{k}+1\right)}\delta^{2}\left(\varepsilon_{k}\varepsilon_{i}+\varepsilon_{k}+\varepsilon_{i}+1\right).$
(264)

Then, we can complete the proof after some simple algebraic simplifications:

$\displaystyle R_{k}$
$\displaystyle\rightarrow\log_{2}\left(1+\frac{p_{k}\mathbb{E}\left\\{\left\|\mathbf{g}_{k}\right\|^{4}\right\\}}{\sum_{i=1,i\neq
k}^{K}p_{i}\mathbb{E}\left\\{\left|\mathbf{g}_{k}^{H}\mathbf{g}_{i}\right|^{2}\right\\}}\right)$
$\displaystyle\rightarrow\log_{2}\left(1+\frac{p_{k}\alpha_{k}\left(2\delta^{2}+2\delta+1\right)}{\sum_{i=1,i\neq
k}^{K}p_{i}\alpha_{i}\delta^{2}}\right),\text{as }M,N\to\infty.$ (265)

Besides, we can find that (C) is a decrease function with respect to $\delta$.

## Appendix D

Firstly, by selecting the non-zero terms when all the Rician factors grow to
infinity, we can complete the derivation of (37).

Secondly, we consider a conventional uplink non-RIS massive MIMO system with
one BS and $K$ users. Assume that the deterministic LoS channel between the BS
and user $k$ is $\sqrt{\gamma_{k}}\mathbf{\bar{h}}_{k}^{w/o}$, where
$\mathbf{\bar{h}}_{k}^{w/o}\in\mathbb{C}^{M\times 1}$. To facilitate the
analysis, we consider the uniform linear array (ULA), and the LoS channel
$\mathbf{\bar{h}}_{k}^{w/o}$ can be expressed as

$\displaystyle\overline{\mathbf{h}}_{k}^{w/o}=\left[1,e^{j2\pi\frac{d}{\lambda}\sin\vartheta_{k}},\ldots,e^{j2\pi\frac{d}{\lambda}(M-1)\sin\vartheta_{k}}\right]^{T},$
(266)

where $\vartheta_{k}$ is the AoA at the BS from user $k$. Besides, we have
$\left\|\mathbf{\bar{h}}_{k}^{w/o}\right\|^{2}=M$.

Thus, with the MRC technique, the rate of user $k$ is given by

$\displaystyle\bar{R}_{k}^{w/o}$
$\displaystyle=\log_{2}\left(1+\frac{p_{k}\left(\gamma_{k}\right)^{2}\left\|\overline{\mathbf{h}}_{k}^{w/o}\right\|^{4}}{\sum_{i=1,i\neq
k}^{K}p_{i}\gamma_{i}\gamma_{k}\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}+\sigma^{2}\gamma_{k}\left\|\overline{\mathbf{h}}_{k}^{w/o}\right\|^{2}}\right),$
$\displaystyle=\log_{2}\left(1+\frac{p_{k}\gamma_{k}M}{\sum_{i=1,i\neq
k}^{K}p_{i}\gamma_{i}\frac{\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}}{M}+\sigma^{2}}\right).$
(267)

Besides, with the ULA structure, we have

$\displaystyle\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}$
$\displaystyle=\left|\sum_{i=1}^{M}e^{-j2\pi\frac{d}{\lambda}(i-1)\sin\vartheta_{k}}e^{j2\pi\frac{d}{\lambda}(i-1)\sin\vartheta_{i}}\right|^{2}$
$\displaystyle=\left|\sum_{i=1}^{M}\left(e^{j2\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}\right)^{(i-1)}\right|^{2}=\left|\frac{1-e^{j2\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}}{1-e^{j2\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}}\right|^{2}$
$\displaystyle=\left|\frac{\left(e^{-j\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}-e^{j\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}\right)e^{j\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}}{\left(e^{-j\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}-e^{j\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}\right)e^{j\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}}\right|^{2}$
$\displaystyle=\left|\frac{\sin\left(\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)\right)}{\sin\left(\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)\right)}e^{j\pi\frac{d}{\lambda}(M-1)\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)}\right|^{2}$
$\displaystyle=\frac{\sin^{2}\left(\pi\frac{d}{\lambda}M\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)\right)}{\sin^{2}\left(\pi\frac{d}{\lambda}\left(\sin\vartheta_{i}-\sin\vartheta_{k}\right)\right)}.$
(268)

Therefore, when $\vartheta_{i}\neq\vartheta_{k}$, we know
$\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}$
is bounded. Thus, when $M\rightarrow\infty$, the inter-user interference terms
in (D) becomes

$\displaystyle\frac{\left|\left(\overline{\mathbf{h}}_{k}^{w/o}\right)^{H}\overline{\mathbf{h}}_{i}^{w/o}\right|^{2}}{M}\rightarrow
0,\text{as }M\to\infty.$ (269)

## Appendix E

If the RIS has discrete phase shifts with $b$ bits precision, the adjustable
phase shifts $\hat{\theta}_{n}$ can only be selected from
$\left\\{0,\frac{2\pi}{2^{b}},2\times\frac{2\pi}{2^{b}},\ldots,\left(2^{b}-1\right)\frac{2\pi}{2^{b}}\right\\}$.