# Reconfigurable Intelligent Surface-Assisted Cell-Free Massive MIMO Systems Over Spatially-Correlated Channels

Trinh van Chien, Hien Quoc Ngo, Symeon Chatzinotas, Marco Di Renzo, Bjorn Ottersten

► To cite this version:

Trinh van Chien, Hien Quoc Ngo, Symeon Chatzinotas, Marco Di Renzo, Bjorn Ottersten. Reconfigurable Intelligent Surface-Assisted Cell-Free Massive MIMO Systems Over Spatially-Correlated Channels. IEEE Transactions on Wireless Communications, 2021, pp.1-1. 10.1109/TWC.2021.3136925. hal-03506040v2

HAL Id: hal-03506040

https://centralesupelec.hal.science/hal-03506040v2

Submitted on 1 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Reconfigurable Intelligent Surface-Assisted Cell-Free Massive
MIMO Systems Over Spatially-Correlated Channels

Trinh Van Chien, Member, IEEE, Hien Quoc Ngo, Senior Member, IEEE, Symeon Chatzinotas, Senior Member, IEEE, Marco Di Renzo, Fellow, IEEE, and Björn Ottersten, Fellow, IEEE

**Abstract**—Cell-Free Massive multiple-input multiple-output (MIMO) and reconfigurable intelligent surface (RIS) are two promising technologies for application to beyond-5G networks. This paper considers Cell-Free Massive MIMO systems with the assistance of an RIS for enhancing the system performance under the presence of spatial correlation among the engineered scattering elements of the RIS. Distributed maximum-ratio processing is considered at the access points (APs). We introduce an aggregated channel estimation approach that provides sufficient information for data processing with the main benefit of reducing the overhead required for channel estimation. The considered system is studied by using asymptotic analysis which lets the number of APs and/or the number of RIS elements grow large. A lower bound for the channel capacity is obtained for a finite number of APs and engineered scattering elements of the RIS, and closed-form expressions for the uplink and downlink ergodic net throughput are formulated in terms of only the channel statistics. Based on the obtained analytical frameworks, we unveil the impact of channel correlation, the number of RIS elements, and the pilot contamination on the net throughput of each user. In addition, a simple control scheme for optimizing the configuration of the engineered scattering elements of the RIS is proposed, which is shown to increase the channel estimation quality, and, hence, the system performance. Numerical results demonstrate the effectiveness of the proposed system design and performance analysis. In particular, the performance benefits of using RISs in Cell-Free Massive MIMO systems are confirmed, especially if the direct links between the APs and the users are of insufficient quality with high probability.

***Index Terms***—Cell-free Massive MIMO, reconfigurable intelligent surface, maximum ratio processing, ergodic net throughput.

I. INTRODUCTION

In the last few decades, we have witnessed an exponential
growth of the demand for wireless communication systems
that provide reliable communications and ensure ubiquitous
coverage, high spectral efficiency and low latency [2]. To meet

Manuscript received April 12, 2021; revised July 30, 2021 and October 28, 2021; accepted December 07, 2021. Date of publication XXX; date of current version XXX. The work of T. V. Chien, S. Chatzinotas, and B. Ottersten was supported by RISOTTI - Reconfigurable Intelligent Surfaces for Smart Cities under project FNR/C20/IS/14773976/RISOTTI. The work of H. Q. Ngo was supported by the UK Research and Innovation Future Leaders Fellowships under Grant MR/S017666/1. The work of M. Di Renzo was supported in part by the European Commission through the H2020 ARIADNE project under grant agreement number 871464 and through the H2020 RISE-6G project under grant agreement number 101017011. Parts of this paper were presented at IEEE GLOBECOM 2021 [1]. The associate editor coordinating the review of this paper and approving it for publication was Z. Zhang. (Corresponding author: Trinh Van Chien.)

T. V. Chien, S. Chatzinotas, B. Ottersten are with the Interdisciplinary Centre for Security, Reliability and Trust (SnT), University of Luxembourg, L-1855 Luxembourg, Luxembourg (email: vanchien.trinh@uni.lu, symeon.chatzinotas@uni.lu, and bjorn.ottersten@uni.lu).

H. Q. Ngo is with the School of Electronics, Electrical Engineering and Computer Science, Queen's University Belfast, Belfast BT7 1NN, United Kingdom (email: hien.ngo@qub.ac.uk).

M. Di Renzo is with Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire des Signaux et Systèmes, 3 Rue Joliot-Curie, 91192 Gif-sur-Yvette, France (email: marco.di-renzo@universite-paris-saclay.fr).

these requirements, several new technologies have been incor-
porated in 5G communication standards, which include Mas-
sive multiple-input multiple-output (MIMO) [3], millimeter-
wave communications [4], and network densification [5].
Among them, Massive MIMO has gained significant attention
since it can offer a good service to many users in the network.
Moreover, the net throughput offered by a Massive MIMO
system is close to the Shannon capacity, in many scenarios,
by only employing simple linear processing techniques, such
as maximum ratio (MR) or zero forcing (ZF) processing.
Since the net throughput can be computed in a closed-
form expression that only depends on the channel statistics,
the optimized designs are applicable for a long period of
time [6]. The colocated Massive MIMO architecture has the
advantage of low backhaul requirements since the base station
antennas are installed in a compact array. Conventional cellular
networks, however, are impaired by intercell interference. In
particular, the users at the cell boundaries are impaired by
high intercell interference and path loss, and hence, they may
experience insufficient performance. More advanced signal
processing methods are necessary to overcome the inherent
intercell interference that characterizes conventional cellular
network deployments.

Cell-Free Massive MIMO has recently been introduced to
reduce the intercell interference that characterizes colocated
Massive MIMO architectures. Cell-Free Massive MIMO is a
network deployment where a large number of access points
(APs) are located in a given coverage area to serve a small
number of users [7]–[10]. All APs collaborate with each other
via a backhaul network and serve all the users in the absence of
cell boundaries. The system performance is enhanced in Cell-
Free Massive MIMO systems because they inherit the benefits
of the distributed MIMO and network MIMO architectures, but
the users are also close to the APs. When each AP is equipped
with a single antenna, MR processing results in a good net
throughput for every user, while ensuring a low computational
complexity and offering a distributed implementation that is
convenient for scalability purposes [11]. However, Cell-Free
Massive MIMO cannot guarantee a good quality of service
under harsh propagation conditions, such as in the presence
of poor scattering environments or high attenuation due to the
presence of large obstacles.

Reconfigurable intelligent surface (RIS) is an emerging
technology that is capable of shaping the radio waves at the
electromagnetic level without applying digital signal process-
ing methods and without requiring power amplifiers [12]–
[14]. Each element of the RIS scatters (e.g., reflects) the
incident signal without using radio frequency chains and power
amplification [15]. Integrating an RIS into wireless networks
introduces digitally controllable links that scale up with the
number of engineered scattering elements of the RIS, whose
estimation is, however, challenged by the lack of digital signal
processing units at the RIS [16]–[20]. For simplicity, the
main attention has so far been concentrated on designing the
phase shifts under the assumption of perfect channel state
information (CSI) [16], [21] and the references therein. In [20],
the authors have recently discussed the fundamental issues
of performing channel estimation in RIS-assisted wireless
systems. The impact of the channel estimation overhead and
reporting on the spectral efficiency, energy efficiency, and
their tradeoff has recently been investigated in [17]. In [16]
and [18], to reduce the impact of the channel estimation
overhead, the authors have investigated the design of RIS-
assisted communications in the presence of statistical CSI. As
far as the integration of Cell-Free Massive MIMO and RIS
is concerned, recent works have formulated and solved op-
timization problems with different communication objectives
under the assumption of perfect (and instantaneous) CSI [22]–
[26]. Recent results in the context of single-input single-output
(SISO) and multi-user MIMO systems have, however, shown
that designs for the engineered scattering elements of the RIS
that are based on statistical CSI may be of practical interest
and provide good performance [18], [27]–[29].

In the depicted context, no prior work has analyzed the
performance of RIS-assisted Cell-Free Massive MIMO sys-
tems in the presence of spatially-correlated channels. In this
work, motivated by these considerations, we introduce an
analytical framework for analyzing and optimizing the uplink
and downlink transmissions of RIS-assisted Cell-Free Massive
MIMO systems under spatially correlated channels and in the
presence of direct links subject to the presence of blockages.
In particular, the main contributions made by this paper can
be summarized as follows:

* We consider an RIS-assisted Cell-Free Massive MIMO under spatially correlated channels. All APs estimate the instantaneous channels in the uplink pilot training phase. We exploit a channel estimation scheme that estimates the aggregated channels including both the direct and indirect links, instead of every individual channel coefficient as in previous works [20], [24]. For generality, the pilot contamination is assumed to originate from an arbitrary pilot reuse pattern.

* We analytically show that, even by using a low complex-
ity MR technique, the non-coherent interference, small-
scale fading effects, and additive noise are averaged out
when the number of APs and RIS elements increases.
The received signal includes, hence, only the desired
signal and the coherent interference. In addition, we show
that the indirect links become dominant if the number of
engineered scattering elements of the RIS increases.

* We derive a closed-form expression of the net throughput for both the uplink and downlink data transmissions. The impact of the array gain, coherent joint transmission, channel estimation errors, pilot contamination, spatial correlation, and phase shifts of the RIS, which determine the system performance, are explicitly observable in the obtained analytical expressions.

* With the aid of numerical simulations, we verify the

<figure><img src="image_2.png" /><figcaption>Fig. 1. An RIS-assisted Cell-Free Massive MIMO system where *M* APs collaborate with each other to serve *K* distant users.</figcaption></figure>

effectiveness of the proposed channel estimation scheme
and the accuracy of the closed-form expressions of the
net throughput. The obtained numerical results show that
the use of RISs significantly enhances the net throughput
per user, especially when the direct links are blocked with
high probability.

The rest of this paper is organized as follows: Section II presents the system model, the channel model, and the channel estimation protocol. The uplink data transmission protocol and the asymptotic analysis by assuming a very large number of APs and engineered scattering elements at the RIS are discussed in Section III. A similar analysis for the downlink data transmission is reported in Section IV. Finally, Section V illustrates several numerical results, while the main conclusions are drawn in Section VI.

*Notation:* Upper and lower bold letters are used to denote matrices and vectors, respectively. The identity matrix of size *N* × *N* is denoted by *I*<sub>*N*</sub>. The imaginary unit of a complex number is denoted by *j* with $\sqrt{j} = -1$. The superscripts (·)*, (·)<sup>T</sup>, and (·)<sup>H</sup> denote the complex conjugate, transpose, and Hermitian transpose, respectively. $\mathbb{E}\{\cdot\}$ and $\mathrm{Var}\{\cdot\}$ denote the expectation and variance of a random variable. The circularly symmetric Gaussian distribution is denoted by $CN(\cdot, \cdot)$ and $diag(\mathbf{x})$ is the diagonal matrix whose main diagonal is given by $\mathbf{x}$. tr(·) is the trace operator. The Euclidean norm of vector $\mathbf{x}$ is $||\mathbf{x}||$, and $||\mathbf{X}||_2$ is the spectral norm of matrix $\mathbf{X}$. Finally, mod(·, ·) is the modulus operation and $[\cdot]_j$ denotes the truncated argument.

## II. SYSTEM MODEL, CHANNEL ESTIMATION, AND RIS
PHASE SHIFT CONTROL

We consider an RIS-assisted Cell-Free Massive MIMO
system, where *M* APs connected to a central processing unit
(CPU) serve *K* users on the same time and frequency resource,
as schematically illustrated in Fig. 1. All APs and users are
equipped with a single antenna and they are randomly located
in the coverage area. Since the considered users are far away
from the APs, the communication is assisted by an RIS that
comprises *N* engineered scattering elements that can modify
the phases of the incident signals.<fnref n="1" /> The matrix of phase shifts

<footnote n="1">¹In general, a completed system should include many users randomly located in the coverage area. There are some favorable users where their links to some APs are strong. But there are also some unfavorable users where the their links to the APs are weak. This may come from the large path loss (long distances) or heavy shadowing. In our work, we consider the cases where RISs are deployed to improve the performance of these unfavorable users, and hence, the coverage of the whole system can be increased.</footnote>
of the RIS is denoted by $\mathbf{\Phi} = \text{diag}([e^{j\theta_1}, \dots, e^{j\theta_N}]^T)$, where $\theta_n \in [-\pi, \pi]$ is the phase shift applied by the n-th element of the RIS. The phase shifts are adjusted by a controller which exchanges information with the APs via a backhaul link (see Fig. 1). As a canonical form of Cell-Free Massive MIMO systems, we assume that the system operates in time-division duplexing (TDD) mode. Thus, we assume that channel reciprocity holds in the consisted system model.

## A. Channel Model

We assume a quasi-static block fading model where the channels are static and frequency flat in each coherence interval comprising $\tau_c$ symbols. We assume that the APs estimate the channel during the uplink pilot training phase. Thus, $\tau_p$ symbols ($\tau_p < \tau_c$) in each coherence interval are dedicated to the channel estimation phase and the remaining ($\tau_c - \tau_p$) symbols are utilized for the uplink and downlink data transmission phases.

The following notation is used: $g_{mk}$ is the channel between the user $k$ and the AP $m$, which is the direct link [12]; $\mathbf{h}_m \in \mathbb{C}^N$ is the channel between the AP $m$ and the RIS; and $\mathbf{z}_k \in \mathbb{C}^N$ is the channel between the RIS and the user $k$. The pair $\mathbf{h}_m$ and $\mathbf{z}_k$, which constitutes the cascaded channel, results in an indirect link (virtual line-of-sight link), which enhances the communication reliability between the AP $m$ and the user $k$ [30]. The majority of existing works assume that the wireless channels undergo uncorrelated Rayleigh fading. In this paper, we consider a more realistic channel model by taking into account the spatial correlation among the engineered scattering elements of the RIS, which is due to their sub-wavelength size, sub-wavelength inter-distance, and geometric layout. In an isotropic propagation environment, in particular, $g_{mk}$, $\mathbf{h}_m$, and $\mathbf{z}_k$ can be modeled as follows

$$g_{mk} \sim \mathcal{CN}(0, \beta_{mk}), \mathbf{h}_m \sim \mathcal{CN}(0, \mathbf{R}_m), \mathbf{z}_k \sim \mathcal{CN}(0, \tilde{\mathbf{R}}_k), \quad (1)$$

where $\beta_{mk}$ is the large-scale fading coefficient; $\mathbf{R}_m \in \mathbb{C}^{N \times N}$ and $\tilde{\mathbf{R}}_k \in \mathbb{C}^{N \times N}$ are the covariance matrices that characterize the spatial correlation among the channels of the RIS elements. The covariance matrices in (1) correspond to a general model, which can be further particularized for application to typical RIS designs and propagation environments. For example, a simple exponential model was used to describe the spatial correlation among the engineered scattering elements of the RIS in [31]. Another recent model that is applicable to isotropic scattering with uniformly distributed multipath components in the half-space in front of the RIS was recently reported in [32], whose covariance matrices are<fnref n="2" />

$$\mathbf{R}_m = \alpha_m d_H d_V \mathbf{R} \text{ and } \tilde{\mathbf{R}}_k = \tilde{\alpha}_k d_H d_V \mathbf{R}, \quad (2)$$

where $\alpha_m, \tilde{\alpha}_k \in \mathbb{C}$ are the large-scale channel coefficients, which, for example, model the signal attenuation due to large objects and due to the transmission distance. The matrices in (2) assume that the size of each element of the RIS is

$d_H \times d_V$, with $d_H$ being the horizontal width and $d_V$ being the vertical height of each RIS element. In particular, the $(m', n')$-th element of the spatial correlation matrix $\mathbf{R} \in \mathbb{C}^{N \times N}$ in (2) is $[\mathbf{R}]_{m'n'} = \text{sinc}(2\|\mathbf{u}_{m'} - \mathbf{u}_{n'}\|/\lambda)$, where $\lambda$ is the wavelength and $\text{sinc}(x) = \sin(\pi x)/(\pi x)$ is the sinc function. The vector $\mathbf{u}_x, x \in \{m', n'\}$ is given by $\mathbf{u}_x = [0, \text{mod}(x-1, N_H)d_H, \lfloor(x-1)/N_H\rfloor d_V]^T$, where $N_H$ and $N_V$ denote the total number of RIS elements in each row and column, respectively. The channel model in (1) is significantly distinct from related works since the small-scale fading and the spatial correlation matrices are included in both links of the virtual line-of-sight link that comprises the RIS. In [31], by contrast, the channels between the transmitters and the RIS are assumed to be deterministic, for analytical tractability.

## B. Uplink Pilot Training Phase

The channels are independently estimated from the $\tau_p$ pilot sequences transmitted by the $K$ users. All the users share the same $\tau_p$ pilot sequences. In particular, $\boldsymbol{\phi}_k \in \mathbb{C}^{\tau_p}$ with $\|\boldsymbol{\phi}_k\|^2 = 1$ is defined as the pilot sequence allocated to the user $k$. We denote by $\mathcal{P}_k$ the set of indices of the users (including the user $k$) that share the same pilot sequence as the user $k$. The pilot sequences are assumed to be mutually orthogonal such that the pilot reuse pattern is

$$\boldsymbol{\phi}_{k'}^{H} \boldsymbol{\phi}_{k} = \begin{cases} 1, & \text{if } k' \in \mathcal{P}_{k}, \\ 0, & \text{if } k' \notin \mathcal{P}_{k}. \end{cases} \quad (3)$$

During the pilot training phase, all the $K$ users transmit the pilot sequences to the $M$ APs simultaneously. In particular, the user $k$ transmits the pilot sequence $\sqrt{\tau_p} \boldsymbol{\phi}_k$. The received training signal at the AP $m$, $\mathbf{y}_{pm} \in \mathbb{C}^{\tau_p}$, can be written as

$$\mathbf{y}_{pm} = \sum_{k=1}^{K} \sqrt{p \tau_p} g_{mk} \boldsymbol{\phi}_k + \sum_{k=1}^{K} \sqrt{p \tau_p} \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k \boldsymbol{\phi}_k + \mathbf{w}_{pm}, \quad (4)$$

where $p$ is the normalized signal-to-noise ratio (SNR) of each pilot symbol, and $\mathbf{w}_{pm} \in \mathbb{C}^{\tau_p}$ is the additive noise at the AP $m$, which is distributed as $\mathbf{w}_{pm} \sim \mathcal{CN}(0, \mathbf{I}_{\tau_p})$. In order for the AP $m$ to estimate the desired channels from the user $k$, the received training signal in (4) is projected on $\boldsymbol{\phi}_k^H$ as

$$y_{pmk} = \boldsymbol{\phi}_k^H \mathbf{y}_{pm} = \sqrt{p \tau_p} (g_{mk} + \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k) + \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sqrt{p \tau_p} (g_{mk'} + \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_{k'}) + \mathbf{w}_{pmk}, \quad (5)$$

where $w_{pmk} = \boldsymbol{\phi}_k^H \mathbf{w}_{pm} \sim \mathcal{CN}(0, 1)$. We emphasize that the co-existence of the direct and indirect channels due to the presence of the RIS results in a complicated channel estimation process. In particular, the cascaded channel in (5) results in a nontrivial procedure for applying the minimum mean-square error (MMSE) estimation method, as reported in previous works, for processing the projected signals [8], [9]. Based on the specific signal structure in (5), we denote the channel between the AP $m$ and the user $k$ through the RIS as

$$u_{mk} = g_{mk} + \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k, \quad (6)$$

which is referred to as the *aggregated channel* that comprises the direct and indirect link between the user $k$ and the AP $m$.

<footnote n="2">²This paper considers a spatial correlation model between engineered scattering elements in the far-field scenarios that is applicable and when the covariance matrices among the users differ only in terms of large-scale channel coefficients. Other scenarios where the users have covariance matrices with different phases [33] are of interest and are left for future work.</footnote>
In contrast to previous works, where the matrix **Φ** of the RIS phase shifts is optimized based on instantaneous CSI, in this paper, **Φ** is optimized based on statistical CSI. This is detailed in Section II-C. By capitalizing on the definition of the aggregated channel in (6), the required channels can be estimated in an effective manner even in the presence of the RIS. In particular, the aggregated channel in (6) is given by the product of weighted complex Gaussian and spatially correlated random variables, as given in (1). Despite the complex analytical form, the following lemma gives information on the statistics of the aggregated channels.

**Lemma 1.** The second and fourth moments of the aggregated channel $u_{mk}$ can be formulated as follows

$$ \mathbb{E}\{|u_{mk}|^2\} = \delta_{mk}, \quad (7) $$

$$ \mathbb{E}\{|u_{mk}|^4\} = 2\delta_{mk}^2 + 2\text{tr}(\mathbf{\Theta}_{mk}^2), \quad (8) $$

where $\mathbf{\Theta}_{mk} = \mathbf{\Phi}^H \mathbf{R}_m \mathbf{\Phi} \tilde{\mathbf{R}}_k$ and $\delta_{mk} = \beta_{mk} + \text{tr}(\mathbf{\Theta}_{mk})$. Moreover, the aggregated channels are mutually independent for $m \neq m'$ and $k \neq k'$, i.e.,

$$ \mathbb{E}\{u_{mk}u_{m'k'}^*\} = \mathbb{E}\{u_{mk}\}\mathbb{E}\{u_{m'k'}^*\} = 0. \quad (9) $$

In addition, the aggregated channels $u_{mk}$ and $u_{m'k}$, $\forall m \neq m'$, and the aggregated channels $u_{mk}$ and $u_{m'k}$, $\forall k \neq k'$, are mutually uncorrelated, i.e.,

$$ \mathbb{E}\{u_{mk}u_{m'k}^*\} = 0 \text{ and } \mathbb{E}\{u_{mk'}u_{mk}^*\} = 0. \quad (10) $$

Besides, the aggregated channels $u_{mk}$, $u_{mk'}$, $u_{m'k}$, and $u_{m'k'}$ fulfill the following conditions

$$ \mathbb{E}\{|u_{mk}u_{m'k'}^*|^2\} = \delta_{mk}\delta_{m'k'}, \quad m \neq m', k \neq k', \quad (11) $$

$$ \mathbb{E}\{u_{mk}^*u_{mk'}u_{m'k'}u_{m'k}\} = \text{tr}(\mathbf{\Theta}_{mk'}\mathbf{\Theta}_{m'k}), \quad m \neq m', k \neq k', \quad (12) $$

$$ \mathbb{E}\{|u_{mk}u_{m'k}^*|^2\} = \delta_{mk}\delta_{m'k} + \text{tr}(\mathbf{\Theta}_{mk}\mathbf{\Theta}_{m'k}), \quad m \neq m', \quad (13) $$

$$ \mathbb{E}\{|u_{mk}u_{mk'}^*|^2\} = \delta_{mk}\delta_{mk'} + \text{tr}(\mathbf{\Theta}_{mk}\mathbf{\Theta}_{mk'}), \quad k \neq k'. \quad (14) $$

*Proof.* See Appendix B. □

The moments in Lemma 1 are employed next for analyzing the channel estimation and the net throughput performance. We note, in addition, that the odd moments of $u_{mk}$, e.g., the first and third moments, are equal to zero. Conditioned on the phase shifts, we employ the linear MMSE method for estimating $u_{mk}$ at the AP. In spite of the complex structure of the RIS-assisted cascaded channel, Lemma 2 provides analytical expressions of the estimated channels.

**Lemma 2.** By assuming that the AP $m$ employs the linear MMSE estimation method based on the signal observation in (5), the estimate of the aggregated channel $u_{mk}$ can be formulated as

$$ \hat{u}_{mk} = \frac{\mathbb{E}\{y_{pmk}^* u_{mk}\} y_{pmk}}{\mathbb{E}\{|y_{pmk}|^2\}} = c_{mk} y_{pmk}, \quad (15) $$

where $c_{mk} = \mathbb{E}\{y_{pmk}^* u_{mk}\}/\mathbb{E}\{|y_{pmk}|^2\}$ has the following closed-form expression

$$ c_{mk} = \frac{\sqrt{p \tau_p} \delta_{mk}}{p \tau_p \sum_{k' \in \mathcal{P}_k} \delta_{mk'} + 1}. \quad (16) $$

The estimated channel in (15) has zero mean and variance $\gamma_{mk}$ equal to

$$ \gamma_{mk} = \mathbb{E}\{|\hat{u}_{mk}|^2\} = \sqrt{p \tau_p} \delta_{mk} c_{mk}. \quad (17) $$

Also, the channel estimation error $e_{mk} = u_{mk} - \hat{u}_{mk}$ and the channel estimate $\hat{u}_{mk}$ are uncorrelated. The channel estimation error has zero mean and variance equal to

$$ \mathbb{E}\{|e_{mk}|^2\} = \delta_{mk} - \gamma_{mk}. \quad (18) $$

*Proof.* It is similar to the proof in [34], and is obtained by applying similar analytical steps to the received signal in (5) and by taking into account the structure of the RIS-assisted cascaded channel and the spatial correlation matrices in (1). □

Lemma 2 shows that, by assuming **Φ** fixed, the aggregated channel in (6) can be estimated without increasing the pilot training overhead, i.e., only $\tau_p$ symbols in each coherence interval are needed for channel estimation, which is the same as for conventional Cell-Free Massive MIMO systems. If the user $k'$ uses the same pilot sequence as the user $k$ does, then $\hat{u}_{mk'} = c_{mk'} y_{pmk}$ from (15). Consequently, we obtain the relation $\hat{u}_{mk'} = \frac{c_{mk'} y_{pmk}}{c_{mk}} \hat{u}_{mk}$, which implies that, because of pilot contamination, it may be difficult to distinguish the signals of two generic users. In that regard it is worth noting that, to get rid of pilot contamination, one can assign mutually orthogonal pilot signals to all the users in the network (if the coherence time is long enough so that $\tau_p \ge K$). Under mutually orthogonal pilot sequences, $c_{mk}$ and $\gamma_{mk}$ simplify to $c_{mk}^o$ and $\gamma_{mk}^o$, respectively, as follows

$$ c_{mk}^o = \frac{\sqrt{p \tau_p} \delta_{mk}}{p \tau_p \delta_{mk} + 1}, \quad \gamma_{mk}^o = \sqrt{p \tau_p} \delta_{mk} c_{mk}^o. \quad (19) $$

This implies that, in the absence of pilot contamination, we have $\gamma_{mk}^o \to \delta_{mk}$ as $\tau_p \to \infty$, i.e., the variance of the channel estimation error in (17) is equal to zero. The channel estimates given in Lemma 2 can be applied to an arbitrary set of phase shifts and covariance matrices. To facilitate the performance analysis presented next, we introduce the following corollary that characterizes the correlation between the aggregated channels and their estimates.

**Corollary 1.** Let us consider the two aggregated channels $u_{mk}$ and $u_{m'k}$ with $m \neq m'$, and let us denote

$$ o_{mk} = \sqrt{\omega_{mk} \hat{u}_{mk}^* u_{mk}} - \sqrt{\omega_{mk}} \mathbb{E}\{\hat{u}_{mk}^* u_{mk}\}, \quad (20) $$

$$ o_{m'k} = \sqrt{\omega_{m'k} \hat{u}_{m'k}^* u_{m'k}} - \sqrt{\omega_{mk}} \mathbb{E}\{\hat{u}_{mk}^* u_{mk}\}, \quad (21) $$

where $\omega_{mk}$ and $\omega_{m'k}$ are two non-negative deterministic scalars. Then, the following holds

$$ \mathbb{E}\{o_{mk}o_{m'k}^*\} = p \tau_p \sqrt{\omega_{mk}\omega_{m'k}} c_{mk}c_{m'k} \sum_{k' \in \mathcal{P}_k} \text{tr}(\mathbf{\Theta}_{mk}\mathbf{\Theta}_{m'k}). \quad (22) $$

*Proof.* See Appendix C. □

Both Lemma 1 and Corollary 1 indicate that the presence of an RIS makes the channel statistics more complex, compared to a conventional Cell-Free Massive MIMO system, due to the correlation among the propagation channels. In the next
sections, the analytical expression of the channel estimates in Lemma 2 and the properties in Corollary 1 and Lemma 1 are employed for signal detection in the uplink and for beamforming in the downlink. Also, the same results are used to optimize the phase shifts of the RIS in order to minimize the channel estimation error and to evaluate the corresponding ergodic net throughput.

### C. RIS Phase Shift Control and Optimization

Channel estimation is a critical aspect in Cell-Free Massive MIMO. As discussed in previous text, in many scenarios, non-orthogonal pilots have to be used. This causes pilot contamination, which may significantly reduce the system performance. In this section, we design an RIS-assisted phase shift control scheme that is aimed to improve the quality of channel estimation. To this end, we introduce the normalized mean square error (NMSE) of the channel estimate of the user $k$ at the AP $m$, as follows

$$ \mathrm{NMSE}_{mk} = \frac{\mathbb{E}\{|e_{mk}|^2\}}{\mathbb{E}\{|u_{mk}|^2\}} = 1 - \frac{p \tau_p \delta_{mk}}{p \tau_p \sum_{k' \in \mathcal{P}_k} \delta_{mk'} + 1}, \quad (23) $$

where the last equality is obtained from (7) and (18). The NMSE is a suitable metric to evaluate the channel estimation quality and to measure the relative channel estimation error per AP. By definition, the NMSE lies in the range [0, 1]. In particular, the NMSE tends to zero if orthogonal pilot signals are used for every user and the pilot power is sufficiently large. In general, however, the NMSE is greater than zero if the $K$ users reuse the pilot signals, i.e., $\tau_p < K$, since $\mathrm{NMSE}_{mk} \to 1 - \delta_{mk}/(\sum_{k' \in \mathcal{P}_k} \delta_{mk'})$ as $p \to \infty$. We propose to optimize the phase shift matrix $\Phi$ of the RIS so as to minimize the total NMSE obtained from all the users and all the APs, as follows

$$ \begin{aligned} \underset{\{\theta_n\}}{\text{minimize}} & \quad \sum_{m=1}^{M} \sum_{k=1}^{K} \mathrm{NMSE}_{mk} \quad (24) \\ \text{subject to} & \quad -\pi \le \theta_n \le \pi, \forall n. \end{aligned} $$

We emphasize that the optimal phase shifts obtained by solving the problem in (24) are independent of the instantaneous CSI and depend only on the statistical CSI, i.e., the large-scale fading coefficients and the channel covariance matrices. Problem (24) is a fractional program, whose globally-optimal solution is not simple to be obtained for an RIS with a large number of independently tunable elements. Nonetheless, in the special network setup where the direct links from the APs to the users are weak enough to be negligible with respect to the RIS-assisted links, the optimal solution to problem (24) is available in a closed-form expression as in Corollary 2.

**Corollary 2.** If the direct links are weak enough to be negligible and the RIS-assisted channels are spatially correlated as formulated in (2), the optimal minimizer of the optimization problem in (24) is $\theta_1 = \dots = \theta_N$, i.e., the equal phase shift design is optimal.

*Proof.* See Appendix D. □

If the direct links are completely blocked and the spatial correlation model in (2) holds, Corollary 2 provides a simple but effective option to design the phase shifts of the

RIS while ensuring the optimal estimation of the aggregated channels according to the sum-NMSE minimization criterion. Therefore, an efficient channel estimation protocol can be designed even in the presence of an RIS equipped with a large number of engineered scattering elements. The numerical results in Section V show that the phase shifts design obtained in Corollary 2 offers good gains in terms of net throughput even if the direct links cannot be completely ignored.

**Remark 1.** The proposed optimization method of the phase shifts of the RIS is based on the minimization of the sum-NMSE, and it is, therefore, based on improving the channel estimation quality. This is a critical objective in Massive MIMO systems, since improving the accuracy of channel estimation results in a noticeable enhancement of the uplink and downlink net throughput [35], [36]. If the direct links are not weak enough, the equal phase shift design is not optimal anymore, and the optimal solution to problem (24) may be obtained numerically. For example, one can compute the first-order derivative of the sum-NMSE in (24) with respect to each reflecting element, and the gradient descent algorithm may be utilized to obtain a locally-optimal solution of (24) in an iterative manner. Another option would be to optimize the phase shifts of the RIS based on the maximization of the uplink or downlink ergodic net throughput. The solution of the corresponding optimization problem is, however, challenging and depends on whether the uplink or the downlink transmission phases are considered. Due to space limitations, therefore, we postpone this latter criterion for optimizing the phase shifts of the RIS to a future research work.

### III. UPLINK DATA TRANSMISSION AND PERFORMANCE ANALYSIS WITH MR COMBINING

In this section, we first introduce a procedure to detect the uplink transmitted signals by capitalizing on the channel estimation method introduced in the previous section. Then, we derive an asymptotic closed-form expression of the ergodic net throughput.

#### A. Uplink Data Transmission Phase

In the uplink, all the $K$ users transmit their data to the $M$ APs simultaneously. Specifically, the user $k$ transmits a modulated symbol $s_k$ with $\mathbb{E}\{|s_k|^2\} = 1$. This symbol is weighted by a power control factor $\sqrt{\eta_k}$, $0 \le \eta_k \le 1$, which enhances the spectral efficiency by, for example, compensating the near-far effects and mitigating the mutual interference among the users. Then, the received baseband signal, $y_{um} \in \mathbb{C}$, at the AP $m$ is

$$ \begin{aligned} y_{um} &= \sqrt{\rho_u} \sum_{k=1}^{K} \sqrt{\eta_k} (g_{mk} + \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k) s_k + w_{um} \quad (25) \\ &= \sqrt{\rho_u} \sum_{k=1}^{K} \sqrt{\eta_k} u_{mk} s_k + w_{um}, \end{aligned} $$

where $\rho_u$ is the normalized uplink SNR of each data symbol, which is defined as the maximum transmit power divided by the noise variance, and $\rho_u \eta_k$ is the corresponding SNR of the user $k$; and $w_{um} \sim CN(0, 1)$ is the normalized additive
noise. For data detection, the MR combining method is used at the CPU, i.e., $\hat{u}_{mk}, \forall m, k$, in (15) is employed to detect the data transmitted by the user $k$. In mathematical terms, the corresponding decision statistic is

$$ 
\begin{aligned}
r_{uk} &= \sum_{m=1}^{M} \hat{u}_{mk}^{*} y_{um} \\
&= \sqrt{\rho_{u}} \sum_{m=1}^{M} \sum_{k'=1}^{K} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'sk'} + \sum_{m=1}^{M} \hat{u}_{mk}^{*} w_{um}.
\end{aligned} 
\quad (26) $$

Based on the observation $r_{uk}$, the uplink ergodic net throughput of the user $k$ is analyzed in the next subsection.

## B. Asymptotic Analysis of the Uplink Received Signal

In the considered system model, the number of APs, $M$, and the number of tunable elements of the RIS, $N$, can be large. Therefore, we analyze the performance of two case studies: (i) $N$ is fixed and $M$ is large; and (ii) both $N$ and $M$ are large. The asymptotic analysis is conditioned upon a given setup of the CSI, which includes the large-scale fading coefficients, the covariance matrices, and the power utilized for the pilot and data transmission phases. To this end, the uplink weighted signal in (26) is split into three terms based on the pilot reuse set $\mathcal{P}_k$, as follows

$$ 
\begin{aligned}
r_{uk} = & \sqrt{\rho_{u}} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'sk'} + \\ & \underbrace{\sqrt{\rho_{u}} \sum_{k' \notin \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'sk'} + \sum_{m=1}^{M} \hat{u}_{mk}^{*} w_{um}}_{\mathcal{T}_{k2}} \\ & \underbrace{\sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'sk'} + \sum_{m=1}^{M} \hat{u}_{mk}^{*} w_{um}}_{\mathcal{T}_{k3}}
\end{aligned} 
\quad (27) $$

where $\mathcal{T}_{k1}$ accounts for the signals received from all the users in $\mathcal{P}_k$, and $\mathcal{T}_{k2}$ accounts for the mutual interference from the users that are assigned orthogonal pilot sequences. The impact of the additive noise obtained after applying MR combining is given by $\mathcal{T}_{k3}$. From (5), (6), and (15), we obtain the following identity

$$ 
\begin{aligned}
& \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'} \\
&= \sum_{m=1}^{M} \sqrt{\eta_{k'}} c_{mk} u_{mk'} \left( \sum_{k'' \in \mathcal{P}_k} \sqrt{p_{\tau p}} u_{mk''}^{*} + w_{pmk}^{*} \right) \\
&= \sum_{m=1}^{M} \sqrt{p_{\tau p} \eta_{k'}} c_{mk} |u_{mk'}|^2 + \sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} \sum_{m=1}^{M} \sqrt{p_{\tau p} \eta_{k''}} \\
& \quad \times c_{mk} u_{mk'} u_{mk''}^{*} + \sum_{m=1}^{M} \sqrt{\eta_{k'}} c_{mk} u_{mk'} w_{pmk}^{*}.
\end{aligned} 
\quad (28) $$

1) Case I: $N$ is fixed and $M$ is large, i.e., $M \to \infty$. In this case, we divide both sides of (28) by $M$ and exploits

Tchebyshev's theorem [37]$^3$ and (7) to obtain

$$ \frac{1}{M} \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'} \xrightarrow[M \to \infty]{P} \frac{1}{M} \sum_{m=1}^{M} \sqrt{p_{\tau p} \eta_{k'}} c_{mk} \delta_{mk'}, \quad (29) $$

where $P \xrightarrow{}$ denotes the convergence in probability.$^4$ Specifically, the second and third terms in (28) converge to zero due to the so-called favorable propagation conditions and since the aggregated channel and the noise are mutually independent [38]. By inserting (29) into the decision variable in (27), we obtain the following deterministic value

$$ \frac{1}{M} r_{uk} \xrightarrow[M \to \infty]{P} \frac{1}{M} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{p_{\tau p} \eta_{k'}} c_{mk} \delta_{mk'sk'}, \quad (30) $$

because $\mathcal{T}_{k2}/M \to 0$ and $\mathcal{T}_{k3}/M \to 0$ as $M \to \infty$. The result in (30) unveils that, for a fixed $N$, the channels become asymptotically orthogonal. In particular, the small-scale fading, the non-coherent interference, and the additive noise vanish. The only residual impairment is the pilot contamination caused by the users that employ the same pilot sequence. This result is the evidence that, due to pilot contamination, the system performance cannot be improved by adding more APs if MR combining is used. The contributions of both the direct and RIS-assisted indirect channels explicitly appear in (30) through the terms $\beta_{mk'}$ and $\text{tr}(\Theta_{mk'})$, respectively.

2) Case II: Both $N$ and $M$ are large, i.e., $N \to \infty$ and $M \to \infty$. To analyze this case study, we need some assumptions on the covariance matrices $\mathbf{R}_m$ and $\tilde{\mathbf{R}}_k$, as summarized as follows.

**Assumption 1.** For $m = 1, \dots, M$ and $k = 1, \dots, K$, the covariance matrices $\mathbf{R}_m$ and $\tilde{\mathbf{R}}_k$ are assumed to fulfill the following properties

$$ \limsup_N \|\mathbf{R}_m\|_2 < \infty, \liminf_N \frac{1}{N} \text{tr}(\mathbf{R}_m) > 0, \quad (31) $$

$$ \limsup_N \|\tilde{\mathbf{R}}_k\|_2 < \infty, \liminf_N \frac{1}{N} \text{tr}(\tilde{\mathbf{R}}_k) > 0. \quad (32) $$

The assumptions in (31) and (32) imply that the largest singular value and the sum of the eigenvalues (counted with their multiplicity) of the $N \times N$ covariance matrices that characterize the spatial correlation among the channels of the RIS elements are finite and positive. Dividing both sides of (28) by $MN$ and then applying Tchebyshev's theorem and (7), we obtain

$$ \frac{1}{MN} \sum_{m=1}^{M} \sqrt{\eta_{k'}} \hat{u}_{mk}^{*} u_{mk'} \xrightarrow[M \to \infty]{P} \frac{1}{MN} \sum_{m=1}^{M} \sqrt{p_{\tau p} \eta_{k'}} c_{mk} \text{tr}(\Theta_{mk'}). \quad (33) $$

We first observe that $\Theta_{mk'}$ is similar to $\tilde{\mathbf{R}}_{k'}^{1/2} \Phi \mathbf{R}_m \Phi^H \tilde{\mathbf{R}}_{k'}^{1/2}$, which is a positive semi-definite matrix.$^5$ Since similar matrix-

<footnote>$^3$Let $X_1, \dots, X_n$ be independent random variables such that $\mathbb{E}\{X_i\} = \bar{x}_i$ and $\text{Var}\{X_i\} \le c < \infty$. Then, Tchebyshev's theorem states $\frac{1}{n} \sum_{n'=1}^{n} \frac{P}{n} \xrightarrow[n \to \infty]{P} \frac{1}{n} \sum_{n'=1}^{n} \bar{x}_{n'}$.</footnote>

<footnote>$^4$A sequence $\{X_n\}$ of random variables converges in probability to the random variable $X$ if, for all $\epsilon > 0$, it holds that $\lim_{n \to \infty} \text{Pr}(|X_n - X| > \epsilon) = 0$, where $\text{Pr}(\cdot)$ denotes the probability of an event.</footnote>

<footnote>$^5$Two matrices $\mathbf{A}$ and $\mathbf{B}$ of size $N \times N$ are similar if there exists an invertible $N \times N$ matrix $\mathbf{U}$ such that $\mathbf{B} = \mathbf{U}^{-1}\mathbf{A}\mathbf{U}$.</footnote>
ces have the same eigenvalues, it follows that $\text{tr}(\Theta_{mk'}) > 0$. Based on Assumption 1, we obtain the following inequalities

$$ 
\begin{align} 
\frac{1}{N} \text{tr}(\Theta_{mk'}) &\stackrel{(a)}{\le} \frac{1}{N} \|\Phi\|_2 \text{tr}(\mathbf{R}_m \mathbf{\Phi} \tilde{\mathbf{R}}_{k'}) \tag{34} \\ 
&\stackrel{(b)}{=} \frac{1}{N} \text{tr}(\mathbf{\Phi} \tilde{\mathbf{R}}_{k'} \mathbf{R}_m) \stackrel{(c)}{\le} \frac{1}{N} \|\tilde{\mathbf{R}}_{k'}\|_2 \text{tr}(\mathbf{R}_m), 
\end{align} 
$$

where (a) is obtained from Lemma 3 in Appendix A; (b) follows because $\|\Phi\|_2 = 1$; and (c) is obtained by applying again Lemma 3. Based on Assumption 1, the last inequality in (34) is bounded by a positive constant. From (33), therefore, the decision variable in (27) can be formulated as

$$ \frac{1}{MN} r_{uk} \xrightarrow[M \to \infty]{P} \frac{1}{MN} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{p \tau_p \rho_u \eta_{k'} c_{mk} \text{tr}(\Theta_{mk'} s_{k'})} \quad (35) $$

which is bounded from above thanks to (34). The expression obtained in (35) reveals that, as $M, N \to \infty$, the post-processed signal at the CPU consists of the desired signal of the intended user $k$ and the interference from the other users in $\mathcal{P}_k$. Compared with (30), we observe that (35) is independent of the direct links and depends only on the RIS-assisted indirect links. This highlights the potentially promising contribution of the RIS, in the limiting regime $M, N \to \infty$, for enhancing the system performance.

### C. Uplink Ergodic Net Throughput with a Finite Number of APs and RIS Elements

In this section, we focus our attention on the practical setup in which $M$ and $N$ are both finite. By utilizing the user-and-then forget channel capacity bounding method [39], the uplink ergodic net throughput of the user $k$ can be written as follows

$$ R_{uk} = B\nu_u (1 - \tau_p / \tau_c) \log_2 (1 + \text{SINR}_{uk}), [\text{Mbps}], \quad (36) $$

where $B$ is the system bandwidth measured in MHz and $0 \le \nu_u \le 1$ is the portion of each coherence time interval that is dedicated to the uplink data transmission phase. The effective uplink signal-to-noise-plus-interference ratio (SINR), which is denoted by $\text{SINR}_{uk}$, is defined as follows

$$ \text{SINR}_{uk} = \frac{|\text{DS}_{uk}|^2}{\mathbb{E}\{|\text{BU}_{uk}|^2\} + \sum_{k'=1, k' \neq k}^{K} \mathbb{E}\{|\text{UI}_{uk'k}|^2\} + \mathbb{E}\{|\text{NO}_{uk}|^2\}}, \quad (37) $$

where the following definitions hold

$$ \text{DS}_{uk} = \sqrt{\rho_u \eta_k} \mathbb{E} \left\{ \sum_{m=1}^{M} \hat{u}_{mk}^* u_{mk} \right\}, \quad (38) $$

$$ \text{BU}_{uk} = \sqrt{\rho_u \eta_k} \left( \sum_{m=1}^{M} \hat{u}_{mk}^* u_{mk} - \mathbb{E} \left\{ \sum_{m=1}^{M} \hat{u}_{mk}^* u_{mk} \right\} \right), \quad (39) $$

$$ \text{UI}_{uk'k} = \sqrt{\rho_u \eta_{k'}} \sum_{m=1}^{M} \hat{u}_{mk}^* u_{mk'}, \quad (40) $$

$$ \text{NO}_{uk} = \sum_{m=1}^{M} \hat{u}_{mk}^* w_{um}. \quad (41) $$

In particular, $\text{DS}_{uk}$ denotes the (average) strength of the desired signal, $\text{BU}_{uk}$ denotes the beamforming uncertainty,

which represents the randomness of the effective channel gain for a given beamforming method, $\text{UI}_{uk'k}$ denotes the interference caused by the user $k'$ to the user $k$, and $\text{NO}_{uk}$ denotes the additive noise. We emphasize that the net throughput in (36) is achievable since it is a lower bound of the channel capacity. A closed-form expression for (36) is given in Theorem 1.

**Theorem 1.** If the CPU utilizes the MR combining method, a lower-bound closed-form expression for the uplink net throughput of the user $k$ is given by (36), where the SINR is

$$ \text{SINR}_{uk} = \frac{\rho_u \eta_k \left( \sum_{m=1}^{M} \gamma_{mk} \right)^2}{\text{MI}_{uk} + \text{NO}_{uk}}, \quad (42) $$

where $\text{MI}_{uk}$ is the mutual interference and $\text{NO}_{uk}$ is the noise, which are formulated as follows

$$ 
\begin{align} 
\text{MI}_{uk} &= \rho_u \sum_{k'=1}^{K} \sum_{m=1}^{M} \eta_{k'} \gamma_{mk} \delta_{mk'} \nonumber \\ 
&\quad + p \tau_p \rho_u \sum_{k'=1}^{K} \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} c_{mk} c_{m'k} \text{tr}(\Theta_{mk'} \Theta_{m'k''}) \tag{43} \\ 
&\quad + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \eta_{k'} c_{mk}^2 \text{tr}(\Theta_{mk'}^2) \nonumber \\ 
&\quad + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \eta_{k'} \left( \sum_{m=1}^{M} c_{mk} \delta_{mk'} \right)^2 \nonumber \\ 
\text{NO}_{uk} &= \sum_{m=1}^{M} \gamma_{mk}, \tag{44} 
\end{align} 
$$

with $\delta_{mk'} = \beta_{mk'} + \text{tr}(\Theta_{mk'}), c_{mk}$ given in (16), and $\gamma_{mk}$ given in (17).

*Proof.* See Appendix E. □

By direct inspection of the SINR in (42), we observe that the numerator increases with the square of the sum of the variances of the channel estimates, $\gamma_{mk}, \forall m$ thanks to the joint signal processing, i.e., the received signals form the $M$ APs are sent to the CPU for centralized data detection. On the other hand, the first term in (43) represents the power of the mutual interference. The use of an RIS to support multiple users introduce the extra interference shown in the second and third terms in (43). Due to the limited and finite number of orthogonal pilot sequences being used, the fourth term in (43) dominates the impact of pilot contamination. The second term in the denominator in (42) is the additive noise. If the coherence time is sufficiently large that every user can utilize its own orthogonal pilot sequence, the uplink net throughput of the user $k$ can still be obtained from (36), but the effective SINR simplifies to (45). The SINR in (42) is a multivariate function of the matrix of phase shifts of the RIS and of the channel statistics, i.e., the channel covariance matrices. Table I gives a comparison of the obtained uplink SINR of the user $k$ with and without the presence of the RIS. By direct inspection of $|\text{DS}_{uk}|^2$, we evince that the strength of the desired signal increases thanks to the assistance of the RIS. However, the mutual interference becomes more severe as
$$ \text{SINR}_{uk} = \frac{\rho_u \eta_k \left( \sum_{m=1}^{M} \gamma_{mk}^o \right)^2}{\sum_{m=1}^{M} \gamma_{mk}^o + \rho_u \sum_{k'=1}^{K} \sum_{m=1}^{M} \eta_{k'} \gamma_{mk}^o \delta_{mk'} + p \tau_p \rho_u \sum_{k'=1}^{K} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} c_{mk}^o c_{m'k}^o \text{tr}(\Theta_{mk'} \Theta_{m'k}) + p \tau_p \rho_u \eta_k \sum_{m=1}^{M} c_{mk}^2 \text{tr}(\Theta_{mk}^2)} \quad (45) $$

well, due to the need of estimating both the direct and indirect links in the presence of the RIS. By assigning orthogonal pilot signals to all the K users, the coherent interference can be completely suppressed. In Section V, the performance of Cell-Free Massive MIMO and RIS-assisted Cell-Free Massive MIMO is compared with the aid of numerical simulations.

## IV. DOWNLINK DATA TRANSMISSION AND PERFORMANCE ANALYSIS WITH MR PRECODING

In this section, we consider the downlink data transmission phase and analyze the received signal at the users when the number of APs is large or when the numbers of RIS elements and APs are both large. A closed-form expression of the downlink ergodic net throughput that is attainable with MR precoding and for an arbitrary phase shift matrix of the RIS elements is provided.

### A. Downlink Data Transmission Phase

By exploiting channel reciprocity, the AP m treats the channel estimates obtained in the uplink as the true channels in order to construct the beamforming coefficients. Accordingly, the downlink signal transmitted from the AP m is⁶

$$ x_m = \sqrt{\rho_d} \sum_{k=1}^{K} \sqrt{\eta_{mk}} \hat{u}_{mk}^* q_k, \quad (46) $$

where $\rho_d$ is the normalized SNR in the downlink; $q_k$ is the complex data symbol that is to be sent (cooperatively by virtue of the considered coherent joint transmission scheme) by all the M APs to the user k, with $\mathbb{E}\{|q_k|^2\} = 1$; and $\eta_{mk}$ is the power control coefficient of the AP m, which satisfies the power budget constraint as follows

$$ \mathbb{E}\{|x_m|^2\} \le \rho_d \Rightarrow \sum_{k=1}^{K} \eta_{mk} \gamma_{mk} \le 1. \quad (47) $$

The cooperation among the M APs for jointly transmitting the same data symbol to a particular user creates the major distinction between the downlink and uplink data transmission phases. Based on (46), the received signal at the user k is the superposition of the signals transmitted by the M APs as

$$ \begin{aligned} r_{dk} &= \sum_{m=1}^{M} u_{mk} x_m + w_{dk} \\ &= \sqrt{\rho_d} \sum_{m=1}^{M} \sum_{k'=1}^{K} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* q_{k'} + w_{dk}. \end{aligned} \quad (48) $$

where $w_{dk} \sim CN(0, 1)$ is the additive noise at the user k. The user k decodes the desired data symbol based on the observation in (48).

⁶In this paper, the downlink data transmission is conducted based on the uplink channel estimates, which depend on the RIS phase shift matrix. Since the MR processing is used for downlink data transmission based on the uplink channel estimates, closed-form analytical expressions for the downlink can be obtained if the same phase shift matrix is utilized in the uplink and in the downlink.

### B. Asymptotic Analysis of the Downlink Received Signal

In contrast with the uplink data processing where the CPU needs only the channel estimate $\hat{u}_{mk}$ for detecting the data of the user k, as displayed in (27), the received signal in (48) depends on the channel estimates of the K users in the network, since the channel estimates from the K users are used for MR precoding. Therefore, the analysis of the uplink and downlink data transmission phases are different. First, we split (48) into three terms, by virtue of the pilot reuse pattern $\mathcal{P}_k$, as follows

$$ r_{dk} = \sqrt{\rho_d} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* q_{k'} + \sqrt{\rho_d} \sum_{k' \notin \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* q_{k'} + w_{dk}. \quad (49) $$

Then, we investigate the two asymptotic regimes for $M \to \infty$ and $M, N \to \infty$. In particular, the first term in (49) can be rewritten as

$$ \begin{aligned} & \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* \\ & \stackrel{(a)}{=} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} c_{mk'} u_{mk} \left( \sum_{k'' \in \mathcal{P}_k} \sqrt{p \tau_p} u_{mk''}^* + w_{pmk'}^* \right) \\ & \stackrel{(b)}{=} \sum_{m=1}^{M} \sqrt{\eta_{mk'} p \tau_p} c_{mk'} |u_{mk}|^2 + \sum_{k'' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \sqrt{\eta_{mk'} p \tau_p} \\ & \quad \times c_{mk'} u_{mk} u_{mk''}^* + \sum_{m=1}^{M} \sqrt{\eta_{mk'}} c_{mk'} u_{mk} w_{pmk'}^* \end{aligned} \quad (50) $$

where (a) is obtained by utilizing the channel estimates in (15) and (b) is obtained by extracting the aggregated channel of the user k from the summation. By letting M and/or N be large, similar to the uplink data transmission phase, we obtain the following asymptotic results

$$ \frac{1}{M} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* \xrightarrow{M \to \infty} \frac{P}{M} \sum_{m=1}^{M} \sqrt{\eta_{mk'} p \tau_p} c_{mk'} \delta_{mk}, \quad (51) $$

$$ \frac{1}{MN} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* \xrightarrow{M \to \infty} \frac{P}{N} \sum_{m=1}^{M} \sqrt{\eta_{mk'} p \tau_p} c_{mk'} \text{tr}(\Theta_{mk}), \quad (52) $$

which are bounded from above based on Assumption 1. Consequently, the received signal at the user k converges to a deterministic equivalent as the number of APs is large, i.e., $M \to \infty$, and as the number of APs and RIS elements are
<table><caption>TABLE I
COMPARISON OF THE UPLINK SINR BETWEEN CELL-FREE MASSIVE MIMO AND RIS-ASSISTED CELL-FREE MASSIVE MIMO</caption><thead><tr><th></th><th>Uplink SINR</th><th>Cell-Free Massive MIMO</th><th>RIS-Assisted Cell-Free Massive MIMO</th></tr></thead><tbody><tr><td rowspan="6">(42)</td><td>δ<sub>mk</sub></td><td>β<sub>mk</sub></td><td>β<sub>mk</sub> + tr(Θ<sup>H</sup>R<sub>m</sub>Φ<sub>R</sub>k)</td></tr><tr><td>c<sub>mk</sub></td><td><sup>√pτ<sub>p</sub>β<sub>mk</sub></sup>/<sub>pτ<sub>p</sub>Σ<sub>k'∈P<sub>k</sub></sub>β<sub>mk'</sub>+1</sub></td><td><sup>√pτ<sub>p</sub>δ<sub>mk</sub></sup>/<sub>pτ<sub>p</sub>Σ<sub>k'∈P<sub>k</sub></sub>δ<sub>mk'</sub>+1</sub></td></tr><tr><td>γ<sub>mk</sub></td><td>√pτ<sub>p</sub>β<sub>mk</sub>c<sub>mk</sub></td><td>√pτ<sub>p</sub>δ<sub>mk</sub>c<sub>mk</sub></td></tr><tr><td>|DS<sub>uk</sub>|<sup>2</sup></td><td>ρ<sub>u</sub>η<sub>k</sub>(Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub>)<sup>2</sup></td><td>ρ<sub>u</sub>η<sub>k</sub>(Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub>)<sup>2</sup></td></tr><tr><td>MI<sub>uk</sub></td><td>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>γ<sub>mk</sub>δ<sub>mk'</sub> +<br>pτ<sub>p</sub>ρ<sub>u</sub>Σ<sub>k'∈P<sub>k</sub>\{k</sub></sub>η<sub>k'</sub>(Σ<sub>m=1</sub><sup>M</sup>c<sub>mk</sub>δ<sub>mk'</sub>)<sup>2</sup></td><td>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>γ<sub>mk</sub>δ<sub>mk'</sub> +<br>pτ<sub>p</sub>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>c<sub>mk</sub>c<sub>m'</sub>tr(Θ<sub>mk'</sub>Θ<sub>m'</sub>k'') +<br>pτ<sub>p</sub>ρ<sub>u</sub>Σ<sub>k'∈P<sub>k</sub>\{k</sub></sub>η<sub>k'</sub>(Σ<sub>m=1</sub><sup>M</sup>c<sub>mk</sub>δ<sub>mk'</sub>)<sup>2</sup></td></tr><tr><td>NO<sub>uk</sub></td><td>Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub></td><td>Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub></td></tr><tr><td rowspan="6">(45)</td><td>c<sub>mk</sub><sup>o</sup></td><td>√pτ<sub>p</sub>β<sub>mk</sub></td><td>√pτ<sub>p</sub>δ<sub>mk</sub></td></tr><tr><td>γ<sub>mk</sub><sup>o</sup></td><td>√pτ<sub>p</sub>β<sub>mk</sub>c<sub>mk</sub><sup>o</sup></td><td>√pτ<sub>p</sub>δ<sub>mk</sub>c<sub>mk</sub><sup>o</sup></td></tr><tr><td>|DS<sub>uk</sub>|<sup>2</sup></td><td>ρ<sub>u</sub>η<sub>k</sub>(Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub><sup>o</sup>)<sup>2</sup></td><td>ρ<sub>u</sub>η<sub>k</sub>(Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub><sup>o</sup>)<sup>2</sup></td></tr><tr><td>MI<sub>uk</sub></td><td>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>γ<sub>mk</sub><sup>o</sup>δ<sub>mk'</sub> + Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub><sup>o</sup></td><td>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>γ<sub>mk</sub><sup>o</sup>δ<sub>mk'</sub> + pτ<sub>p</sub>ρ<sub>u</sub>Σ<sub>k'=1</sub><sup>K</sup>Σ<sub>m=1</sub><sup>M</sup>η<sub>k'</sub>c<sub>mk</sub><sup>o</sup>c<sub>m'</sub>tr(Θ<sub>mk'</sub>Θ<sub>m'</sub>k') +<br>pτ<sub>p</sub>ρ<sub>u</sub>η<sub>k'</sub>Σ<sub>m=1</sub><sup>M</sup>c<sub>mk</sub><sup>o</sup>tr(Θ<sub>mk</sub><sup>2</sup>)</td></tr><tr><td>NO<sub>uk</sub></td><td>Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub><sup>o</sup></td><td>Σ<sub>m=1</sub><sup>M</sup>γ<sub>mk</sub><sup>o</sup></td></tr></tbody></table>

large, i.e., $M, N \to \infty$. More precisely, the received signal converges (asymptotically) to

$$ \frac{1}{M} r_{dk} \xrightarrow{P} \frac{1}{M} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{p \tau_p \rho_d \eta_{mk'} c_{mk'} \delta_{mk} s_{k'}} \quad (53) $$

$$ \frac{1}{MN} r_{dk} \xrightarrow{P} \frac{1}{MN} \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sqrt{p \tau_p \rho_d \eta_{mk'} c_{mk'} \text{tr}(\Theta_{mk}) s_{k'}} \quad (54) $$

which indicates the inherent coexistence of the users in $\mathcal{P}_k$. The deterministic equivalents in (53) and (54) unveil that the impact of the channel estimation accuracy and the channel statistics is different between the uplink and the downlink. In particular, the asymptotic received signal in the uplink only depends on the channel estimation quality of each individual user, which is manifested by the coefficient $c_{mk}$. The asymptotic received signal in the downlink depends, on the other hand, on the channel estimation quality of all the users that share the same orthogonal pilot sequences, i.e., $c_{m,k'}, \forall k' \in \mathcal{P}_k$.

### C. Downlink Ergodic Net Throughput with a Finite Number of APs and RIS Elements

By utilizing the channel capacity bounding technique [39], similar to the analysis of the uplink data transmission phase, the downlink ergodic net throughput of the user $k$ can be written as follows

$$ R_{dk} = Bv_d (1 - \tau_p / \tau_c) \log_2 (1 + \text{SINR}_{dk}) \quad [\text{Mbps}] \quad (55) $$

where $0 \le v_d \le 1$ is the portion of each coherence time interval dedicated to the downlink data transmission phase,

with $v_u + v_d = 1$, and the effective downlink SINR is defined as

$$ \text{SINR}_{dk} = \frac{|\text{DS}_{dk}|^2}{\mathbb{E}\{|\text{BU}_{dk}|^2\} + \sum_{k'=1, k' \neq k}^{K} \mathbb{E}\{|\text{UI}_{dk'k}|^2\} + 1} \quad (56) $$

where the following definitions hold

$$ \text{DS}_{dk} = \sqrt{\rho_d} \mathbb{E} \left\{ \sum_{m=1}^{M} \sqrt{\eta_{mk}} u_{mk} \hat{u}_{mk}^* \right\} \quad (57) $$

$$ \text{UI}_{dk'k} = \sqrt{\rho_d} \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk}^* \quad (58) $$

$$ \text{BU}_{dk} = \sqrt{\rho_d} \left( \sum_{m=1}^{M} \sqrt{\eta_{mk}} u_{mk} \hat{u}_{mk}^* - \mathbb{E} \left\{ \sum_{m=1}^{M} \sqrt{\eta_{mk}} u_{mk} \hat{u}_{mk}^* \right\} \right) \quad (59) $$

In particular, $\text{DS}_{dk}$ denotes the (average) strength of the desired signal received by the user $k$, $\text{BU}_{dk}$ denotes the beam-forming uncertainty, $\text{UI}_{dk'k}$ denotes the interference caused to the user $k$ by the signal intended to the user $k'$. The downlink ergodic net throughput in (55) is achievable since it is a lower bound of the channel capacity, similar to the uplink data transmission phase. In contrast to the uplink ergodic net throughput, which only depends on the combining coefficients of each individual user, the downlink net throughput of the user $k$ depends on the precoding coefficients of all the $K$ users. A closed-form expression for (55) is given in Theorem 2.

**Theorem 2.** If the CPU utilizes the MR precoding method, a lower-bound closed-form expression for the downlink net throughput of the user $k$ is given by (55), where the SINR is

$$ \text{SINR}_{dk} = \frac{\rho_d (\sum_{m=1}^{M} \sqrt{\eta_{mk}} \gamma_{mk})^2}{\text{MI}_{dk} + 1} \quad (60) $$
$$ \text{SINR}_{dk} = \frac{\rho_d \left( \sum_{m=1}^{M} \sqrt{\eta_{mk} \gamma_{mk}} \right)^2}{1 + \rho_d \sum_{k'=1}^{K} \sum_{m'=1}^{M} \eta_{mk'} \gamma_{mk'}^o \delta_{mk} + p \tau_p \rho_d \sum_{k'=1}^{K} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \text{tr}(\Theta_{mk} \Theta_{m'k'}) + p \tau_p \rho_d \sum_{m=1}^{M} \eta_{mk} c_{mk}^2 \text{tr}(\Theta_{mk}^2)}} \quad (62) $$

where $Ml_{dk}$ is the mutual interference, which is defined as follows

$$ \begin{aligned} Ml_{dk} ={}& \rho_d \sum_{k'=1}^{K} \sum_{m=1}^{M} \eta_{mk'} \gamma_{mk'} \delta_{mk} + p \tau_p \rho_d \times \\ & \sum_{k'=1}^{K} \sum_{k'' \in \mathcal{P}_{k'}} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \text{tr}(\Theta_{mk} \Theta_{m'k''})} \\ & + p \tau_p \rho_d \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \text{tr}(\Theta_{mk}^2) + \\ & p \tau_p \rho_d \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \left( \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} \delta_{mk}} \right)^2. \end{aligned} \quad (61) $$

*Proof.* The main steps of the proof are similar to those of the proof of Theorem 1. However, there are also major differences that are due to the coherent joint data transmission scheme among the APs. The details of the proof are available in Appendix F. □

From Theorem 2, we observe that the effective downlink SINR has some similarities and differences as compared with its uplink counterpart in Theorem 1. Similar to the uplink, the numerator of (60) is a quadratic function that depends on the channel estimation quality and the coherent joint transmission processing. Differently from the uplink, the transmit power coefficients appear explicitly in the numerator of (60) as a result of the cooperation among the APs. The impact of pilot contamination in (61), in contrast to the uplink, depends on all the transmit power coefficients. In addition, we observe that the impact of pilot contamination scales up with the number of APs and with the number of elements of the RIS. When the K users employ orthogonal pilot sequences, the downlink SINR expression of the user k simplifies to (62).

A comparison of the obtained analytical expressions of the downlink SINR for Cell-Free Massive MIMO and RIS-assisted Cell-Free Massive MIMO systems is given in Table II. By comparing Table I and Table II, the difference and similarities between the uplink and downlink transmission phases can be identified as well. With the aid of numerical results, in Section V, we will illustrate the advantages of RIS-assisted Cell-Free Massive MIMO especially if the direct links are not sufficiently reliable (e.g., they are blocked) with high probability.

**Remark 2.** We observe that the ergodic net throughput in the uplink (Theorem 1) and downlink (Theorem 2) data transmission phases depend only on the large-scale fading statistics and on the channel covariance matrices, while they are independent of the instantaneous CSI. This simplifies the deployment and optimization of RIS-assisted Cell-Free Massive MIMO systems. As anticipated in Remark 1, in fact,

the phase shifts of the RIS can be optimized based on the (long-term) analytical expressions of the ergodic net throughputs in Theorem 1 and Theorem 2, which are independent of the instantaneous CSI. In this paper, we have opted for optimizing the phase shifts of the RIS in order to minimize the channel estimation error, which determines the performance of both the uplink and downlink transmission phases. The optimization of the phase shifts of the RISs based on the closed-form expressions in Theorem 1 and Theorem 2 is postponed to a future research work.

## V. NUMERICAL RESULTS

In this section, we report some numerical results in order to illustrate the performance of the RIS-assisted Cell-Free Massive MIMO system introduced in the previous sections. We consider a geographic area of size $1.5 \times 1.5$ km$^2$, where the locations of the APs and users are given in terms of $(x, y)$ coordinates. The four vertices of the considered region are $[-0.75, -0.75]$ km, $[-0.75, 0.75]$ km, $[0.75, 0.75]$ km, $[0.75, -0.75]$ km. To simulate a harsh communication environment, the M APs are uniformly distributed in the sub-region $x, y \in [-0.75, -0.5]$ km, while the K users are uniformly distributed in the sub-region $x, y \in [0.375, 0.75]$ km. The RIS is located at the origin, i.e., $(x, y) = (0, 0)$. The carrier frequency is 1.9 GHz and the system bandwidth is 20 MHz. Each coherence interval comprises $\tau_c = 200$ symbols, which may correspond to a coherence bandwidth equal to $B_c = 200$ KHz and a coherence time equal to $T_c = 1$ ms, except in Fig. 12(b) where $\tau_c = 5000$ symbols. We assume $\tau_p = 5$ orthonormal pilot sequences that are shared by all the users. The large-scale fading coefficients in dB are generated according to the three-slope propagation model in [8, (51)–(53)], where the path loss exponent depends on the distance between the transmitter and the receiver. The shadow fading has a log-normal distribution with standard deviation equal to 8 dB. The distance thresholds for the three slopes are 10 m and 50 m. The height of the APs, RIS, and users is 15 m, 30 m, and 1.65 m, respectively. The direct links, $g_{mk}, \forall m, k$, are assumed to be unblocked with a given probability. More specifically, the large-scale fading coefficient $\beta_{mk}$ is formulated as follows

$$ \beta_{mk} = \bar{\beta}_{mk} a_{mk}, \quad (63) $$

where $\bar{\beta}_{mk}$ accounts for the path loss due to the transmission distance and the shadow fading according to the three-slope propagation model in [8]. The binary variables $a_{mk}$ accounts for the probability that the direct links are unblocked, and it is defined as

$$ a_{mk} = \begin{cases} 1, & \text{with a probability } \tilde{p}, \\ 0, & \text{with a probability } 1 - \tilde{p}, \end{cases} \quad (64) $$
<table><caption>TABLE II
COMPARISON OF THE DOWNLINK SINR BETWEEN CELL-FREE MASSIVE MIMO AND RIS-ASSISTED CELL-FREE MASSIVE MIMO (SOME PARAMETERS ARE DEFINED IN TABLE I)</caption><thead><tr><th>Downlink SINR</th><th>Cell-Free Massive MIMO</th><th>RIS-Assisted Cell-Free Massive MIMO</th></tr></thead><tbody><tr><td rowspan="3">(60)</td><td>|DS<sub>dk</sub>|<sup>2</sup></td><td>ρ<sub>d</sub> (M/Σ<sub>m=1</sub> √η<sub>mk</sub>γ<sub>mk</sub>)<sup>2</sup></td><td>ρ<sub>d</sub> (M/Σ<sub>m=1</sub> √η<sub>mk</sub>γ<sub>mk</sub>)<sup>2</sup></td></tr><tr><td>M<sub>ldk</sub></td><td>ρ<sub>d</sub> (Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>γ<sub>mk'</sub> δ<sub>mk</sub> + pτ<sub>p</sub>ρ<sub>d</sub> Σ<sub>k'=∈P<sub>k</sub>\{k}</sub> (M/Σ<sub>m=1</sub> √η<sub>mk'</sub>c<sub>mk'</sub> δ<sub>mk</sub>)<sup>2</sup></td><td>ρ<sub>d</sub> (Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>γ<sub>mk'</sub> δ<sub>mk</sub> + pτ<sub>p</sub>ρ<sub>d</sub> Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>η<sub>m'k'</sub>c<sub>mk'</sub>c<sub>m'k'</sub>tr(Θ<sub>mk</sub>Θ<sub>m'k'</sub>) + pτ<sub>p</sub>ρ<sub>d</sub> Σ<sub>k'=∈P<sub>k</sub>\{k}</sub> (M/Σ<sub>m=1</sub> √η<sub>mk'</sub>c<sub>mk'</sub> δ<sub>mk</sub>)<sup>2</sup></td></tr><tr><td>NO<sub>dk</sub></td><td>1</td><td>1</td></tr><tr><td rowspan="3">(62)</td><td>|DS<sub>dk</sub>|<sup>2</sup></td><td>ρ<sub>d</sub> (M/Σ<sub>m=1</sub> √η<sub>mk</sub>γ<sub>mk</sub><sup>o</sup>)<sup>2</sup></td><td>ρ<sub>d</sub> (M/Σ<sub>m=1</sub> √η<sub>mk</sub>γ<sub>mk</sub><sup>o</sup>)<sup>2</sup></td></tr><tr><td>M<sub>ldk</sub></td><td>ρ<sub>d</sub> (Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>γ<sub>mk'</sub><sup>o</sup> δ<sub>mk</sub> + 1</td><td>ρ<sub>d</sub> (Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>γ<sub>mk'</sub><sup>o</sup> δ<sub>mk</sub> + pτ<sub>p</sub>ρ<sub>d</sub> Σ<sub>k'=1</sub><sup>K</sup> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk'</sub>η<sub>m'k'</sub>c<sub>mk'</sub>c<sub>m'k'</sub>tr(Θ<sub>mk</sub>Θ<sub>m'k'</sub>) + pτ<sub>p</sub>ρ<sub>d</sub> Σ<sub>m=1</sub><sup>M</sup> η<sub>mk</sub>c<sub>mk</sub><sup>2</sup>tr(Θ<sub>mk</sub><sup>2</sup>)</td></tr><tr><td>NO<sub>dk</sub></td><td>1</td><td>1</td></tr></tbody></table>

where $\tilde{p} \in [0, 1]$ is the probability that the direct link is not blocked. The noise variance is $-92$ dBm, which corresponds to a noise figure of 9 dB. The covariance matrices are generated according to the spatial correlation model in (2). The power of the pilot sequences is 100 mW and the power budget of each AP is 200 mW. The time intervals of the data transmission phase, in each coherence time, that are allocated to the uplink and downlink transmissions are, $v_u = v_d = 0.5$. The uplink and downlink power control coefficients are $\eta_k = 1, \forall k$, and $\eta_{mk} = (\sum_{k'=1}^K \gamma_{mk'})^{-1}, \forall m, k$, which is directly obtained from (47) to satisfy the limited power budget per AP.<fnref n="7" /> As far as the optimization of the phase shifts of the RIS elements are concerned, we assume that they are optimized according to the sum-NMSE minimization criterion in the absence of direct links, according to Corollary 2 and Remark 1. Without loss of generality, in particular, the $N$ phase shifts in $\Phi$ are all set equal to $\pi/4$, except in Figs. 8 and 9 where different phase shifts are considered for comparison. In order to evaluate the advantages and limitations of RIS-assisted Cell-Free Massive MIMO systems, three system configurations are considered for comparison:

i) **RIS-Assisted Cell-Free Massive MIMO:** This is the proposed system model in which the direct links are blocked according to (64). This setup is denoted by "RIS-CellFree".

ii) **(Conventional) Cell-Free Massive MIMO:** This setup is the same as the previous one with the only exception that the RIS is not deployed in the network. This setup is denoted by "CellFree".

iii) **RIS-Assisted Cell-Free Massive MIMO with blocked direct links:** This is the worst case study in which the direct links are blocked with unit probability and the uplink and downlink transmission phases are ensured only

through the RIS. This setup is denoted by "RIS-CellFree-NoLOS".

In Fig. 2, we illustrate the cumulative distribution function (CDF) of the net throughput by using Monte Carlo simulations and the proposed analytical framework. The CDF is computed with respect to the locations of the APs and users in the considered area. The Monte Carlo simulation results are obtained by using (36) and (55) with the SINRs given in (37) and (56), while the analytical results are obtained by using Theorem 1 and Theorem 2. We observe a very good overlap between the numerical simulations and the obtained analytical expressions. From Fig. 2, we evince that the downlink net throughput per user is about $2.6 \times$ better than the uplink net throughput. This is due to the higher transmission power of the APs and the gain of the joint processing of the APs. Since the Monte Carlo simulations are not simple to obtain for larger values of the simulation parameters, the rest of the figures are obtained by using the closed-form expressions of the net throughput derived in Theorem 1 and Theorem 2.

In Fig. 3, we illustrate the average sum net throughput as a function of the probability $\tilde{p}$ in (64). In particular, the average uplink sum net throughput is defined as $\sum_{k=1}^K \mathbb{E}\{R_{uk}\}$ and the downlink sum net throughput is defined as $\sum_{k=1}^K \mathbb{E}\{R_{dk}\}$, where the uplink and downlink SINRs are obtained by using Theorem 1 and Theorem 2. In particular, the expectation is computed with respect to the locations of the APs and users in the considered area. From the obtained results, we evince that Cell-Free Massive MIMO provides the worst performance if the blocking probability is large ($\tilde{p}$ is small). As expected, the average net throughput offered by Cell-Free Massive MIMO tends to zero if $\tilde{p} \to 0$ (the direct links are unreliable). For example, at $\tilde{p} = 0.1$, the average sum net throughput of Cell-Free Massive MIMO is approximately $2.0 \times$ and $1.4 \times$ smaller, in the uplink and downlink, respectively, than the average sum net throughput of the worst-case RIS-assisted Cell-Free Massive MIMO setup (i.e., RIS-CellFree-NoLOS). In the considered case study, in addition, we note that the

<footnote n="7">⁷In the worst case, if all the direct links are blocked, we introduce a damping constant when Cell-Free Massive MIMO systems in the absence of the RIS are considered, since in those cases we have $\sum_{m=1}^M \gamma_{mk'} = 0$.</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 2. Monte Carlo simulations versus the analytical frameworks with $M = 20$, $K = 5$, $N = 64$, $\tau_p = 2$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 1.0$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 4. CDF of the sum net throughput [Mbps] with $M = 100$, $K = 10$, $N = 900$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 6. Average sum net throughput [Mbps] versus the number of users with $M = 100$, $N = 900$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 3. Average sum net throughput [Mbps] versus the unblocked probability of the direct links $\bar{p}$ with $M = 100$, $K = 10$, $N = 900$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Average sum net throughput [Mbps] versus the number of APs with $K = 10$, $N = 900$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Average sum net throughput [Mbps] versus the number of pilot sequences with $M = 100$, $K = 10$, $N = 900$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

proposed RIS-assisted Cell-Free Massive MIMO setup offers the best average net throughput, since it can overcome the unreliability of the direct links thanks to the presence of the RIS. The presence of the RIS is particularly useful if $\bar{p}$ is small, i.e., $\bar{p} < 0.2$ in Fig. 3, since the direct links are not able to support a high throughput. In this case, the combination of Cell-Free Massive MIMO and RIS is capable of providing a high throughput and signal reliability.

In Fig. 4, we compare the three considered systems in terms of average sum net throughput when $\bar{p} = 0.2$. We observe the net advantage of the proposed RIS-assisted Cell-Free Massive

MIMO system, especially in the downlink. In the uplink, in addition, even the worst-case RIS-assisted Cell-Free Massive MIMO system setup (i.e., $\bar{p} = 0$) outperforms the Cell-Free Massive MIMO setup in the absence of an RIS. In Figs. 5–7, we show the average sum net throughput as a function of the number of APs, the number of users, and the number of orthonormal pilot signals, respectively. We observe that the average sum net throughput increases with the number of APs, with the number of users, and with the number of pilot sequences, and that the RIS-assisted Cell-Free Massive MIMO setup outperforms, especially in the downlink, the
<figure><img src="image_9.png" /><figcaption>Fig. 8. Average uplink sum net throughput [Mbps] versus the number of engineered scattering elements with $M = 100$, $K = 10$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Average downlink sum net throughput [Mbps] versus the number of engineered scattering elements with $M = 100$, $K = 10$, $\tau_p = 5$, and $d_H = d_V = \lambda/4$. The unblocked probability of the direct links is $\bar{p} = 0.2$.</figcaption></figure>

other benchmark schemes. Gains of the order of 1.7× and 2.6× are obtained in the considered setups.

In the following figures, we focus our attention only on the RIS-assisted Cell-Free Massive MIMO setup, since it provides the best performance in the analyzed setups. In Figs. 8 and 9, we report the uplink and downlink average sum net throughput as a function of number of engineered scattering elements of the RIS. In particular, we compare the average sum net throughput when the phase shifts of the RIS are randomly chosen and are optimized according to Corollary 2 over spatially-independent and spatially-correlated fading channels according to (2). In the presence of spatial correlation, the channel correlation matrices are $\mathbf{R}_m = \alpha_m d_H d_V \mathbf{I}_N$ and $\tilde{\mathbf{R}}_{mk} = \tilde{\alpha}_{mk} d_H d_V \mathbf{I}_N, \forall m, k$. We see different performance trends over spatially-independent and spatially-correlated fading channels. If the spatial correlation is not considered, we observe that there is no significant difference between the random and uniform phase shifts setup. In the presence of spatial uncor- relation on the other hand, the uniform phase shift design obtained from Corollary 2 provides a much higher average throughput. This result highlights the relevance of using even simple optimization designs for RIS-assisted communications over spatially-correlated fading channels.

In Fig. 10, we analyze the impact of the size of the engineered scattering elements of the RISs on the uplink and downlink average net throughput, while keeping the total number of RIS elements $N$ fixed. The size of the considered RIS, which is a compact surface, is $Nd_H d_V$, which implies that it increases as the size $d_H d_V$ of each element of the RIS increases. In this setup, we observe that the average net throughput increases as the physical size of each element of the RIS increases. In Fig. 11, on the other hand, we analyze a setup in which the total size of the RIS is kept constant and equal to $Nd_H d_V = 10\lambda \times 10\lambda$ while the triplet $(N, d = d_H = d_V)$ is changed accordingly. With the considered fading spatial correlation model and for a size of the RIS elements no smaller than $\lambda/3$, we do not observe a significant difference on the average net throughput. Further studies are, however, necessary for deep sub-wavelength RIS structures, for different optimization criteria of the phase shifts of the RIS, and in the presence of mutual coupling in addition to the fading spatial correlation [40], [41].

In Fig. 12, we plot the ergodic net throughput per user by utilizing different channel capacity bounding techniques, by assuming $\tau_c = 200$ symbols and $\tau_c = 5000$ symbols. The main benefit of the use-and-then-forget bound in (36) and the hardening bound in (55) is the possibility of obtaining a closed-form expression for the net throughput by capitalizing on the fundamentals properties of Massive MIMO communi- cations, as demonstrated in Theorems 1 and 2. However, the channel hardening capability reduces in the presence of an RIS [38]. Thus, other channel capacity bounding techniques may result in a better estimate of the net throughput as compared to the actual channel capacity. The statistical bound in [42], [43] was originally derived for application to the downlink data transmission when no instantaneous CSI is available at the users and the channel vectors may be less hardened [10]. In order to apply the same bound to the uplink data transmission, we assume that only the channel statistics are available at the CPU. Even though the statistical bound results in a better ergodic net throughput per user, some realizations of user locations and shadow fading may lead to negative values of the ergodic net throughput. This may occur in high mobility scenarios, which correspond to small values of $\tau_c$, as illustrated in Fig. 12(a). The statistical bound provides consistent values of the ergodic net throughput in low mobility scenarios when $\tau_c$ is large, as displayed in Fig. 12(b). The numerical results unveils the need of developing different and more accurate channel bounding methods for evaluating the net throughput of RIS-assisted Cell-Free Massive MIMO systems.

## VI. CONCLUSION

Cell-Free Massive MIMO and RIS are two disruptive technologies for boosting the system performance of future wireless networks. These two technologies are not competing with each other, but have complementary features that can be integrated and leveraged for enhancing the system per- formance in harsh communication environments. Therefore, we have considered an RIS-assisted Cell-Free Massive MIMO system that operates according to the TDD mode. An efficient channel estimation scheme has been introduced to overcome the high overhead that may be associated with the estimation of the individual channels of the RIS elements. Based on the
<figure><img src="image_11.png" /><figcaption>Fig. 10. Average sum net throughput [Mbps] versus the size of engineered scattering elements {$d_H, d_V$}, but for a different size of the RIS with $M = 100, K = 10, N = 900$, and $\tau_p = 5$. The unblocked probability of the direct links is $\tilde{p} = 0.2$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 12. A comparison of different channel capacity bounding techniques with $M = 20, K = 10, N = 64, \tau_p = 10, d_H = d_V = \lambda/4, \tilde{p} = 1.0$, and two coherence intervals: (a) $\tau_c = 200$ symbols; (b) $\tau_c = 5000$ symbols.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. Average sum net throughput [Mbps] for a fixed total size of the RIS but for a different number of RIS elements with $M = 100, K = 10, \tau_p = 5$, and $d = d_H = d_V$. The unblocked probability of the direct links is $\tilde{p} = 0.2$.</figcaption></figure>

<figure><img src="image_14.png" /></figure>

proposed channel estimation scheme, an optimal design for the phase shifts of the RIS that minimizes the channel estimation error has been devised and has been used for system analysis. Based on the proposed channel estimation method, closed-form expressions of the ergodic net throughput for the uplink and downlink data transmission phases have been proposed. Based on them, the performance of RIS-assisted Cell-Free Massive MIMO has been analyzed as a function of the fading spatial correlation and the blocking probability of the direct AP-user links. The numerical results have shown that the presence of an RIS is particularly useful if the AP-user links are unreliable with high probability.

Possible generalizations of the results illustrated in this paper include the optimization of the phase shifts of the RIS that maximize the uplink or downlink throughput, the analysis of the impact of the fading spatial correlation for non-compact and deep sub-wavelength RIS structures, and the analysis and optimization of RIS-assisted systems in the presence of mutual coupling.

# APPENDIX

## A. Useful Lemmas

This section reports three useful lemmas that are utilized for asymptotic analysis.

**Lemma 3.** [44, Lemma B.7] For an arbitrary matrix $\mathbf{X} \in \mathbb{C}^{N \times N}$ and a positive semi-definite matrix $\mathbf{Y} \in \mathbb{C}^{N \times N}$, it holds

that $|\text{tr}(\mathbf{X}\mathbf{Y})| \le \|\mathbf{X}\|_2\text{tr}(\mathbf{Y})$. If $\mathbf{X}$ is also a positive semi-definite matrix, then $\text{tr}(\mathbf{X}\mathbf{Y}) \le \|\mathbf{X}\|_2\text{tr}(\mathbf{Y})$.

**Lemma 4.** [45, Lemma 9] For a random variable $\mathbf{x} \in \mathbb{C}^N$ distributed as $CN(0, \bar{\mathbf{R}})$ with $\bar{\mathbf{R}} \in \mathbb{C}^{N \times N}$ and a given deterministic matrix $\mathbf{M} \in \mathbb{C}^{N \times N}$, it holds that $\mathbb{E}\{|\mathbf{x}^H\mathbf{M}\mathbf{x}|^2\} = |\text{tr}(\bar{\mathbf{R}}\mathbf{M})|^2 + \text{tr}(\bar{\mathbf{R}}\mathbf{M}\bar{\mathbf{R}}\mathbf{M}^H)$.

**Lemma 5.** For a random vector $\mathbf{x} \in \mathbb{C}^N$ distributed as $\mathbf{x} \sim CN(0, \bar{\mathbf{R}})$ with $\bar{\mathbf{R}} \in \mathbb{C}^{N \times N}$ and two deterministic matrices $\mathbf{M}, \mathbf{N} \in \mathbb{C}^{N \times N}$, it holds that

$$ \mathbb{E}\{\mathbf{x}^H\mathbf{M}\mathbf{x}^H\mathbf{N}\mathbf{x}\} = \text{tr}(\bar{\mathbf{R}}\mathbf{M}\bar{\mathbf{R}}\mathbf{N}) + \text{tr}(\bar{\mathbf{R}}\mathbf{M})\text{tr}(\bar{\mathbf{R}}\mathbf{N}). \quad (65) $$

*Proof.* Consider $\mathbf{x} = \bar{\mathbf{R}}^{1/2}\tilde{\mathbf{x}}$ with $\tilde{\mathbf{x}} \sim CN(0, \mathbf{I}_N)$. Let us further denote $\tilde{\mathbf{M}} = \bar{\mathbf{R}}^{1/2}\mathbf{M}\bar{\mathbf{R}}^{1/2}$ and $\tilde{\mathbf{N}} = \bar{\mathbf{R}}^{1/2}\mathbf{N}\bar{\mathbf{R}}^{1/2}$, where $[\bar{\mathbf{M}}]_{mn}$ and $[\tilde{\mathbf{N}}]_{mn}$ are the $(m, n)$-th elements of matrix $\bar{\mathbf{M}}$ and $\tilde{\mathbf{N}}$, respectively. Then, the expectation on the left-hand side of (65) is

$$ \begin{aligned} & \mathbb{E}\{\mathbf{x}^H\mathbf{M}\mathbf{x}^H\mathbf{N}\mathbf{x}\} \\ &= \mathbb{E}\left\{\left(\sum_{m=1}^N \sum_{n=1}^N \tilde{x}_m^*[\bar{\mathbf{M}}]_{mn}\tilde{x}_n\right)\left(\sum_{m=1}^N \sum_{n=1}^N \tilde{x}_m^*[\tilde{\mathbf{N}}]_{mn}\tilde{x}_n\right)\right\} \quad (66) \\ &= \sum_{m=1}^N \sum_{n=1}^N \sum_{m'=1}^N \sum_{n'=1}^N [\bar{\mathbf{M}}]_{mn} [\tilde{\mathbf{N}}]_{m'n'} \mathbb{E}\{\tilde{x}_m^*[\tilde{x}_n]\tilde{x}_{m'}^*[\tilde{x}_{n'}]\}. \end{aligned} $$
where $\tilde{x}_n$ is the n-th element of vector $\tilde{x}$. By noting that $\mathbb{E}\{|\tilde{x}_m|^4\} = 2$ from Lemma 4 and $\mathbb{E}\{|\tilde{x}_m|^2|\tilde{x}_n|^2\} = 1$ if $m \neq n$, we obtain the following

$$ \mathbb{E}\{\tilde{x}_m^* \tilde{x}_n \tilde{x}_{m'}^* \tilde{x}_{n'}\} = \begin{cases} 2, & \text{if } m = n = m' = n', \\ 1, & \text{if } (m = n) \neq (m' = n'), \\ 1, & \text{if } (m = n') \neq (m' = n), \\ 0, & \text{otherwise.} \end{cases} \quad (67) $$

Consequently, (66) can be further simplified as

$$ \mathbb{E}\{\mathbf{x}^H \mathbf{M}\mathbf{x}^H \mathbf{N}\mathbf{x}\} = \sum_{m=1}^{N} \sum_{n=1}^{N} [\tilde{\mathbf{M}}]_{mm} [\tilde{\mathbf{N}}]_{nn} + \sum_{m=1}^{N} \sum_{n=1}^{N} [\tilde{\mathbf{M}}]_{mn} [\tilde{\mathbf{N}}]_{nm}, \quad (68) $$

which, with the aid of some algebraic manipulations, coincides with (65). □

## B. Proof of Lemma 1

We first compute the second moment of the aggregated channel $u_{mk}$ by capitalizing on the statistical independence of the direct and indirect channels, as follows

$$ \begin{aligned} \mathbb{E}\{|u_{mk}|^2\} &= \mathbb{E}\{|g_{mk}|^2\} + \mathbb{E}\{|\mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k|^2\} \\ &\stackrel{(a)}{=} \beta_{mk} + \mathbb{E}\{\text{tr}(\mathbf{\Phi}^H \mathbf{h}_m \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k \mathbf{z}_k^H)\} \\ &\stackrel{(b)}{=} \beta_{mk} + \text{tr}(\mathbf{\Phi}^H \mathbb{E}\{\mathbf{h}_m \mathbf{h}_m^H\} \mathbf{\Phi} \mathbb{E}\{\mathbf{z}_k \mathbf{z}_k^H\}) \\ &= \beta_{mk} + \text{tr}(\mathbf{\Theta}_{mk}) = \delta_{mk}, \end{aligned} \quad (69) $$

where (a) follows by applying the trace of product property $\text{tr}(\mathbf{X}\mathbf{Y}) = \text{tr}(\mathbf{Y}\mathbf{X})$ for some given size-matched matrices $\mathbf{X}$ and $\mathbf{Y}$; and (b) is obtained thanks to the independence of the cascaded channels $\mathbf{h}_m$ and $\mathbf{z}_k$. The fourth moment of the aggregated channel can be written, from (6), as follows

$$ \mathbb{E}\{|u_{mk}|^4\} = \mathbb{E}\left\{\left|g_{mk}\right|^2 + g_{mk}^* \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k + g_{mk} \mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_m + |\mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k|^2\right|^2\} \quad (70) $$

By setting $a = |g_{mk}|^2$, $b = g_{mk}^* \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k$, $c = g_{mk} \mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_m$, and $d = |\mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k|^2$, (70) can be equivalently written as follows

$$ \mathbb{E}\{|u_{mk}|^4\} = \mathbb{E}\{|a|^2\} + \mathbb{E}\{|b|^2\} + \mathbb{E}\{|c|^2\} + 2\mathbb{E}\{ad\} + \mathbb{E}\{|d|^2\}. \quad (71) $$

By applying Lemma 4 with $g_{mk} \sim \mathbf{CN}(0, \beta_{mk})$, the first expectation on the right-hand side of (71) is equal to

$$ \mathbb{E}\{|a|^2\} = 2\beta_{mk}^2. \quad (72) $$

By exploiting the independence between the direct and RIS-assisted links, the next three expectations on the right-hand side of (71) are equal to

$$ \mathbb{E}\{|b|^2\} = \mathbb{E}\{|c|^2\} = \mathbb{E}\{ad\} = \beta_{mk} \text{tr}(\mathbf{\Theta}_{mk}). \quad (73) $$

By introducing the normalized variable $\tilde{z} = \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k / \|\mathbf{R}_m^{1/2} \mathbf{\Phi} \mathbf{z}_k\|$ with $\tilde{z} \sim \mathbf{CN}(0, 1)$, the last expectation on the right-hand side of (71) is equal to

$$ \begin{aligned} \mathbb{E}\{|d|^2\} &= \mathbb{E}\left\{\left\|\mathbf{R}_m^{1/2} \mathbf{\Phi} \mathbf{z}_k\right\|^4 |\tilde{z}|^4\right\} \\ &\stackrel{(a)}{=} \mathbb{E}\left\{\left\|\mathbf{R}_m^{1/2} \mathbf{\Phi} \mathbf{z}_k\right\|^4\right\} \mathbb{E}\{|\tilde{z}|^4\} = 2(\text{tr}(\mathbf{\Theta}_{mk}))^2 + 2\text{tr}(\mathbf{\Theta}_{mk}^2), \end{aligned} \quad (74) $$

where (a) follows because $\tilde{z}$ is independent of the remaining random variables; and (b) is obtained by virtue of Lemma 4. Inserting (72)–(74) into (71), the proof follows with the aid of some algebraic manipulations.

Also, by exploiting the second moment in (7), the expectation of the two independent aggregated channels can be formulated as shown in (11). The correlation between the two aggregated channels $u_{mk}$ and $u_{m'k}$, $m \neq m'$, is, by definition, as follows

$$ \begin{aligned} \mathbb{E}\{u_{mk} u_{m'k}^*\} &= \mathbb{E}\{g_{mk} g_{m'k}^*\} + \mathbb{E}\{g_{mk}(\mathbf{h}_{m'}^H \mathbf{\Phi} \mathbf{z}_k)^*\} \\ &\quad + \mathbb{E}\{\mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k g_{m'k}^*\} + \mathbb{E}\{\mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k (\mathbf{h}_{m'}^H \mathbf{\Phi} \mathbf{z}_k)^*\} \stackrel{(a)}{=} 0, \end{aligned} \quad (75) $$

where (a) follows because the propagation channels are independent.

The expectation in (13) can be written as follows

$$ \begin{aligned} \mathbb{E}\{|u_{mk} u_{m'k}^*|^2\} &= \mathbb{E}\{|u_{mk}|^2 |u_{m'k}|^2\} \\ &= \mathbb{E}\{|g_{mk}|^2 |g_{m'k}|^2\} + \mathbb{E}\{|g_{mk}|^2 |\mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_{m'} \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k|^2\} \\ &\quad + \mathbb{E}\{|g_{m'k}|^2 |\mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_m \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k|^2\} \\ &= \mathbb{E}\{\mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_m \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k \mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_{m'} \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k\} \quad (76) \\ &= \beta_{mk} \beta_{m'k} + \beta_{mk} \text{tr}(\mathbf{\Theta}_{m'k}) + \beta_{m'k} \text{tr}(\mathbf{\Theta}_{mk}) \\ &\quad + \mathbb{E}\{\mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{R}_m \mathbf{\Phi} \mathbf{z}_k \mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{R}_{m'} \mathbf{\Phi} \mathbf{z}_k\} \\ &\stackrel{(a)}{=} (\beta_{mk} + \text{tr}(\mathbf{\Theta}_{mk})) (\beta_{m'k} + \text{tr}(\mathbf{\Theta}_{m'k})) + \text{tr}(\mathbf{\Theta}_{mk} \mathbf{\Theta}_{m'k}), \end{aligned} $$

where (a) is obtained by utilizing Lemma 5. Similar steps can be applied to the two aggregated channels $u_{mk}$ and $u_{m'k}$ with $k \neq k'$.

The expectation in (12) can be written as follows

$$ \begin{aligned} \mathbb{E}\{u_{mk}^* u_{m'k}^* u_{m'k'} u_{m'k}\} \\ &\stackrel{(a)}{=} \mathbb{E}\{\mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_m \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k' \mathbf{z}_k^H \mathbf{\Phi}^H \mathbf{h}_{m'} \mathbf{h}_m^H \mathbf{\Phi} \mathbf{z}_k\} \\ &= \text{tr}(\mathbf{\Phi}^H \mathbf{R}_m \mathbf{\Phi} \bar{\mathbf{R}}_{k'} \mathbf{\Phi}^H \mathbf{R}_{m'} \mathbf{\Phi} \bar{\mathbf{R}}_k), \end{aligned} \quad (77) $$

where (a) follows from the independence of the direct links.

## C. Proof of Corollary 1

By utilizing the identities $\mathbb{E}\{(X - \mathbb{E}\{X\})(Y - \mathbb{E}\{Y\})\} = \mathbb{E}\{XY\} - \mathbb{E}\{X\}\mathbb{E}\{Y\}$ and $\mathbb{E}\{|X - \mathbb{E}\{X\}|^2\} = \mathbb{E}\{|X|^2\} - |\mathbb{E}\{X\}|^2$, the expectation in (22) can be formulated as follows

$$ \begin{aligned} \mathbb{E}\{o_{mk} o_{m'k}^*\} &= \sqrt{\alpha_{mk} \alpha_{m'k}} \underbrace{\mathbb{E}\{\hat{u}_{mk}^* u_{mk} \hat{u}_{m'k} u_{m'k}^*\}}_{=Q_{mm'}} \\ &= \sqrt{\alpha_{mk} \alpha_{m'k} \gamma_{mk} \gamma_{m'k}}. \end{aligned} \quad (78) $$

By using the analytical expressions of the projected training signal in (5) and the channel estimate in (15), $Q_{mm'k}$ defined in
(78) can be formulated in a closed-form expression, as follows

$$
\begin{align}
Q_{mm'k} &\stackrel{(a)}{=} c_{mk} c_{m'k} \mathbb{E} \left\{ \left( \sqrt{p \tau_p} \sum_{k' \in \mathcal{P}_k} u_{mk'}^* + w_{pmk} \right) u_{mk} \times \right. \nonumber \\
&\qquad \left. \left( \sqrt{p \tau_p} \sum_{k' \in \mathcal{P}_k} u_{m'k'} + w_{pm'k} \right) u_{m'k'}^* \right\} \nonumber \\
&\stackrel{(b)}{=} c_{mk} c_{m'k} p \tau_p \mathbb{E} \{ |u_{mk}|^2 |u_{m'k'}|^2 \} + c_{mk} c_{m'k} p \tau_p \times \nonumber \\
&\qquad \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \mathbb{E} \{ u_{mk'}^* u_{mk} u_{m'k'}^* u_{m'k'} \} \nonumber \\
&= c_{mk} c_{m'k} p \tau_p \delta_{mk} \delta_{m'k} + c_{mk} c_{m'k} p \tau_p \sum_{k' \in \mathcal{P}_k} \operatorname{tr}(\Theta_{mk} \Theta_{m'k'}), \tag{79}
\end{align}
$$

where (a) is obtained by retaining only the terms whose expectation is not zero based on (10) and (b) follows by utilizing (13). From (17), finally, we obtain

$$
Q_{mm'k} = \gamma_{mk}\gamma_{m'k} + c_{mk}c_{m'k}p\tau_p \sum_{k' \in \mathcal{P}_k} \operatorname{tr}(\Theta_{mk}\Theta_{m'k'}). \quad (80)
$$

The proof follows by inserting (78) in (80).

D. *Proof of Corollary 2*

We introduce the shorthand notation $b_{mk} = p\tau_p\alpha_m\tilde{\alpha}_k d_H^2 d_V^2$ and $d_{mk} = p\tau_p d_H^2 d_V^2 \alpha_m \sum_{k' \in \mathcal{P}_k} \tilde{\alpha}_{k'}$. When the direct links are weak enough to be negligible, the NMSE of the channel estimate of the user $k$ at the AP $m$ can be reformulated as follows

$$
\mathrm{NMSE}_{mk} = 1 - \frac{b_{mk} \mathrm{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})}{1 + d_{mk} \mathrm{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})}. \quad (81)
$$

Let us denote by $f(\operatorname{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})) = \sum_{m=1}^{M} \sum_{k=1}^{K} \mathrm{NMSE}_{mk}$ the objective function of the problem in (24), the first-order derivative of $\mathrm{NMSE}_{mk}$ with respect to $\operatorname{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})$ is

$$
\frac{df(\operatorname{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R}))}{d\operatorname{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})} = - \sum_{m=1}^{M} \sum_{k=1}^{K} \frac{b_{mk}}{(1 + d_{mk} \operatorname{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R}))^2} < 0, \quad (82)
$$

which implies that the objective function is a monotonically
decreasing function of tr(Φ^H RΦR) since b_{mk} ≥ 0, ∀m, k.
Moreover, Φ^H RΦR is similar to R^{1/2}Φ^H RΦR^{1/2}, which is
a positive semidefinite matrix. Thus, we obtain

$$
\mathrm{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R}) = |\mathrm{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R})| = \left| \sum_{n=1}^{N} \sum_{n'=1}^{N} (r^{nn'})^2 e^{j(\theta_n - \theta_{n'})} \right|. \quad (83)
$$

Let us introduce the two vectors **a**, **b** ∈ ℂ<sup>N<sup>2</sup></sup> defined as follows

$$
\mathbf{a} = [r^{11}, \dots, r^{nn'}, \dots, r^{NN}]^T, \quad (84)
$$

$$
\mathbf{b} = [r^{11} e^{j(\theta_1 - \theta_1)}, \dots, r^{nn'} e^{j(\theta_n - \theta_{n'})}, \dots, r^{NN} e^{j(\theta_N - \theta_N)}]^T. \quad (85)
$$

With the aid of Cauchy-Schwarz’s inequality, we obtain the
following upper bound for (83)

$$
\mathrm{tr}(\mathbf{\Phi}^H \mathbf{R} \mathbf{\Phi} \mathbf{R}) = |\mathbf{a}^H \mathbf{b}| \le \|\mathbf{a}\| \|\mathbf{b}\| = \sum_{n=1}^{N} \sum_{n'=1}^{N} (r^{nn'})^2, \quad (86)
$$

which holds with equality if and only if the two vectors **a** and
**b** in (84) and (85) are parallel. This implies θn = θn', ∀n, n'.
By combining (82) and (86), the proof is concluded.

E. *Proof of Theorem 1*

To obtain the closed-form expression of the uplink SINR in (37), we first compute $|DS_k|^2$ by using the definition of $u_{mk}$ in (6) as

$$
\begin{align}
|DS_{uk}|^2 &\stackrel{(a)}{=} \rho_u \eta_k \left| \mathbb{E} \left\{ \sum_{m=1}^{M} \hat{u}_{mk}^* (\hat{u}_{mk} + e_{mk}) \right\} \right|^2 \tag{87} \\
&\stackrel{(b)}{=} \rho_u \eta_k \left| \sum_{m=1}^{M} \mathbb{E} \{|\hat{u}_{mk}|^2\} \right|^2 = \rho_u \eta_k \left( \sum_{m=1}^{M} \gamma_{mk} \right)^2, \nonumber
\end{align}
$$

where (a) is obtained by expressing the original channel $u_{mk}$ into the summation of its channel estimate and its estimation error as stated in Lemma 2; and (b) follows because the channel estimate and the channel estimation error are uncorrelated. Since the aggregated channels sharing the same AP index are correlated, the first expectation in the denominator of (37) is equal to

$$
\begin{align}
\mathbb{E}\{|\mathbf{B}u_{uk}|^2\} &= \rho_u \eta_k \mathbb{E}\left\{\left|\sum_{m=1}^{M} o_{umk}\right|^2\right\} \nonumber \\
&= \rho_u \eta_k \sum_{m=1}^{M} \sum_{m'=1, m' \neq m}^{M} \mathbb{E}\{o_{umk} o_{um'k}^*\} + \rho_u \eta_k \sum_{m=1}^{M} \mathbb{E}\{|o_{umk}|^2\}, \nonumber \\
&\qquad \underbrace{=}_{= \bar{T}_{u0}} \nonumber \\
&= \rho_u \eta_k \sum_{m=1}^{M} \sum_{m'=1, m' \neq m}^{M} \mathbb{E}\{|o_{umk}|^2\}, \nonumber \\
&\qquad \underbrace{=}_{= \bar{T}_{u1}} \tag{88}
\end{align}
$$

where $o_{umk} = \hat{u}_{mk}^* u_{mk} - \mathbb{E}\{\hat{u}_{mk}^* u_{mk}\}$ with $\omega_{mk} = 1$. The closed-form expression of $\bar{T}_{u0}$ defined in (88) is as follows

$$
\bar{T}_{u0} = p \tau_p \rho_u \eta_k \sum_{k' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1, m' \neq m}^{M} c_{mk} c_{m'k} \operatorname{tr}(\Theta_{mk} \Theta_{m'k'}), \quad (89)
$$

thanks to (22) in Corollary 1. The expectation of $\bar{T}_{u1}$ defined
in (88) can be rewritten as follows

$$
\begin{align}
\bar{T}_{u1} &= \rho_u \eta_k \sum_{m=1}^{M} \mathbb{E} \left\{ |\hat{u}_{mk}^* u_{mk} - \mathbb{E}\{\hat{u}_{mk}^* u_{mk}\}|^2 \right\} \nonumber \\
&\stackrel{(a)}{=} \rho_u \eta_k \sum_{m=1}^{M} \mathbb{E} \left\{ |\hat{u}_{mk}^* u_{mk}|^2 \right\} - \rho_u \eta_k \sum_{m=1}^{M} \mathbb{E} \left\{ |\hat{u}_{mk}^* u_{mk}|^2 \right\} \tag{90} \\
&\stackrel{(b)}{=} \rho_u \eta_k \sum_{m=1}^{M} \mathbb{E} \left\{ |\hat{u}_{mk}^* u_{mk}|^2 \right\} - \rho_u \eta_k \sum_{m=1}^{M} \gamma_{mk}^2, \nonumber
\end{align}
$$

where (a) is obtained by using the identity $\mathbb{E}\{|X - \mathbb{E}\{X\}|^2\} = \mathbb{E}\{|X|^2\} - |\mathbb{E}\{X\}|^2$; and (b) is obtained from $u_{mk}^* = \hat{u}_{mk}^* + e_{mk}^*$, by taking into account that the channel estimate and the channel estimation error are uncorrelated random variables as stated in Lemma 2. By replacing $\hat{u}_{mk}^* = c_{mk} \gamma_{mk}^*$ in the first
expectation of (88), we obtain

$$
\begin{align}
\mathbb{E}\{|\hat{u}_{mk}^* u_{mk}|^2\} &= c_{mk}^2 \mathbb{E}\{|y_{pmk}^* u_{mk}|^2\} \nonumber \\
&= c_{mk}^2 \mathbb{E}\left\{\left|\left(\sum_{k' \in \mathcal{P}_k} \sqrt{p\tau_p} u_{mk'}^* + w_{pmk}^*\right) u_{mk}\right|^2\right\} \nonumber \\
&= \underbrace{c_{mk}^2 p \tau_p \mathbb{E}\{|u_{mk}|^4\} + c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \mathbb{E}\{|u_{mk'}^* u_{mk}|^2\}}_{=T_{11}} \nonumber \\
&\quad + \underbrace{c_{mk}^2 \mathbb{E}\{|w_{pmk}^* u_{mk}|^2\}}_{=T_{12}}. \tag{91}
\end{align}
$$

Let us analyze the three terms, $T_{11}$, $T_{12}$, and $T_{13}$, in the last equality of (91). The first term can be computed by exploiting the fourth moment given in (8), as follows

$$
\begin{align}
T_{11} &= 2c_{mk}^2 p \tau_p \delta_{mk}^2 + 2c_{mk}^2 p \tau_p \text{tr}(\mathbf{\Theta}_{mk}^2) \tag{92} \\
&\stackrel{(a)}{=} \gamma_{mk}^2 + c_{mk}^2 p \tau_p \delta_{mk}^2 + 2c_{mk}^2 p \tau_p \text{tr}(\mathbf{\Theta}_{mk}^2), \nonumber
\end{align}
$$

where (a) is obtained by using the variance of the channel estimate in (17), with $\delta_{mk}$ and $\xi_{mk}$ that are defined in the statement of the theorem. The second term can be computed by exploiting the uncorrelation of the two cascaded channels, as follows

$$
\begin{align}
T_{12} &= c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \mathbb{E}\{|u_{mk'}|^2|u_{mk}|^2\} \nonumber \\
&= c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \delta_{mk'} \delta_{mk} + c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{mk}). \tag{93}
\end{align}
$$

The last term can be computed by exploiting the independence of the channel and noise, as follows

$$
\begin{align}
T_{13} &= c_{mk}^2 \mathbb{E}\{|w_{pmk}|^2\} \mathbb{E}\{|u_{mk}|^2\} \tag{94} \\
&= c_{mk}^2 \delta_{mk}. \nonumber
\end{align}
$$

By inserting (92)–(94) into (91) and with the aid of some algebraic steps, we obtain

$$
\begin{align}
\mathbb{E}\{|\hat{u}_{mk}^* u_{mk}|^2\} &= c_{mk}^2 p \tau_p \text{tr}(\mathbf{\Theta}_{mk}^2) + \gamma_{mk}^2 \nonumber \\
&\quad + c_{mk} \sqrt{p \tau_p} \delta_{mk}^2 + c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{mk}) \tag{95} \\
&= c_{mk}^2 p \tau_p \text{tr}(\mathbf{\Theta}_{mk}^2) + \gamma_{mk}^2 + \gamma_{mk} \delta_{mk} \nonumber \\
&\quad + c_{mk}^2 p \tau_p \sum_{k' \in \mathcal{P}_k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{mk}), \nonumber
\end{align}
$$

where the final identity is obtained by using the relationship between $\gamma_{mk}$ and $c_{mk}$ in (17). Combining (88) and (95), the first term in the denominator of (37) simplifies to

$$
\begin{align}
\mathbb{E}\{|\mathbf{B}U_{uk}|^2\} ={}& \rho_u \eta_k \sum_{m=1}^{M} \gamma_{mk} \delta_{mk} + p \tau_p \rho_u \eta_k \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} c_{mk} \nonumber \\
& \times c_{m'k} \text{tr}(\mathbf{\Theta}_{mk} \mathbf{\Theta}_{m'k''}) + p \tau_p \rho_u \eta_k \sum_{m=1}^{M} c_{mk}^2 \text{tr}(\mathbf{\Theta}_{mk}^2). \tag{96}
\end{align}
$$

The second term in the denominator of (37) can be split into the two terms based on the pilot reuse pattern defined in $\mathcal{P}_k$, as follows

$$
\sum_{k'=1, k' \neq k}^{K} \mathbb{E}\{|\mathbf{U}_{uk'k}|^2\} = \sum_{k' \notin \mathcal{P}_k} \mathbb{E}\{|\mathbf{U}_{uk'k}|^2\} \\
\qquad + \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \mathbb{E}\{|\mathbf{U}_{uk'k}|^2\}. \quad (97)
$$

The first term on the right-hand side of (97) is the non-coherent interference, which is equal to

$$
\begin{align}
\sum_{k' \notin \mathcal{P}_k} \mathbb{E}\{|\mathbf{U}_{uk'k}|^2\} &= \rho_u \sum_{k' \notin \mathcal{P}_k} \eta_{k'} \mathbb{E}\left\{\left|\sum_{m=1}^{M} \hat{u}_{mk}^* u_{mk'}\right|^2\right\} \tag{98} \\
&= \rho_u \sum_{k' \notin \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} \underbrace{\mathbb{E}\{\hat{u}_{mk}^* u_{mk'} \hat{u}_{m'k} u_{m'k'}^*\}}_{=T_{mk'm'k}}. \nonumber
\end{align}
$$

We compute each expectation $T_{mk'm'k}$ by utilizing the channel estimate in (15) with the aid of some algebraic manipulations, as follows

$$
\begin{align}
T_{mk'm'k} ={}& p \tau_p c_{mk} c_{m'k} \sum_{k'' \in \mathcal{P}_k} \mathbb{E}\{u_{mk''}^* u_{mk'} u_{m'k''} u_{m'k'}^*\} \nonumber \\
& + c_{mk} c_{m'k} \mathbb{E}\{w_{pmk}^* u_{mk'} w_{pm'k} u_{m'k'}^*\} \tag{99} \\
={}& \begin{cases}
p \tau_p c_{mk} c_{m'k} \sum_{k'' \in \mathcal{P}_k} \bar{T}_{mk'm'k''}, & \text{if } m \neq m', \\
\bar{T}_{mk'mk}, & \text{if } m = m'.
\end{cases} \nonumber
\end{align}
$$

where the following definitions hold

$$
\begin{align}
\bar{T}_{mk'm'k''} &= \mathbb{E}\{u_{mk''}^* u_{mk'} u_{m'k''} u_{m'k'}^*\}, \tag{100} \\
\bar{T}_{mk'mk} &= p \tau_p c_{mk}^2 \sum_{k'' \in \mathcal{P}_k} \mathbb{E}\{|u_{mk''}^* u_{mk'}|^2\} + c_{mk}^2 \mathbb{E}\{|u_{mk'}|^2\}. \tag{101}
\end{align}
$$

The closed-form expression of $\bar{T}_{mk'm'k''}$ is obtained by utilizing the uncorrelated property of the quadruple aggregated channels in (12), as follows

$$
\bar{T}_{mk'm'k''} = \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k''}). \quad (102)
$$

In addition, the closed-form expression of $\bar{T}_{mk'mk}$ can be computed by utilizing the results in Lemma 1, as follows

$$
\begin{align}
\bar{T}_{mk'mk} &= c_{mk}^2 p \tau_p \sum_{k'' \in \mathcal{P}_k} \delta_{mk''} \delta_{mk'} + c_{mk}^2 p \tau_p \times \nonumber \\
&\quad \sum_{k'' \in \mathcal{P}_k} \text{tr}(\mathbf{\Theta}_{mk''} \mathbf{\Theta}_{mk'}) + c_{mk}^2 \delta_{mk'} \tag{103} \\
&= \gamma_{mk} \delta_{mk'} + c_{mk}^2 p \tau_p \sum_{k'' \in \mathcal{P}_k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{mk''}). \nonumber
\end{align}
$$

Inserting (102) and (103) into (99), and with the aid of some algebraic manipulations, (98) can be equivalently formulated as follows

$$
\sum_{k' \notin \mathcal{P}_k} \mathbb{E}\{|\mathbf{U}_{uk'k}|^2\} = \rho_u \sum_{k' \notin \mathcal{P}_k} \sum_{m=1}^{M} \eta_{k'} \gamma_{mk} \delta_{mk'} + p \tau_p \rho_u \times \\
\sum_{k' \notin \mathcal{P}_k} \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} c_{mk} c_{m'k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k''}). \quad (104)
$$
The second term on the right-hand side of (97) is the coherent
interference. By using (5) and (15), it simplifies as follows

$$
\begin{align}
\mathbb{E}\{|\mathrm{UI}_{uk'k}|^2\} &= \mathbb{E}\left\{\left|\sum_{m=1}^{M} c_{mk} y_{pmk}^* u_{mk'}\right|^2\right\} \nonumber \\
&= \mathbb{E}\left\{\left|\sum_{m=1}^{M} c_{mk} \left(\sum_{k'' \in \mathcal{P}_k} \sqrt{p_{\tau p}} u_{mk''}^* + w_{pmk}^*\right) u_{mk'}\right|^2\right\} \nonumber \\
&= \rho_u \eta_{k'} \mathbb{E}\left\{\left|\sum_{m=1}^{M} c_{mk} w_{pmk}^* u_{mk'}\right|^2\right\} \nonumber \\
&\quad + \rho_u \eta_{k'} p_{\tau p} \mathbb{E}\left\{\underbrace{\left|\sum_{m=1}^{M} c_{mk} \left(\sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} u_{mk''}^* u_{mk'}\right)\right|^2}_{=T_{u21}}\right\} \nonumber \\
&\quad + \rho_u \eta_{k'} p_{\tau p} \mathbb{E}\left\{\underbrace{\left|\sum_{m=1}^{M} c_{mk} |u_{mk'}|^2\right|^2}_{=T_{u22}}\right\} \nonumber \\
&= \rho_u \eta_{k'} \sum_{m=1}^{M} c_{mk}^2 \delta_{mk'} + \rho_u \eta_{k'} p_{\tau p} (T_{u21} + T_{u22}), \tag{105}
\end{align}
$$

which is obtained by using the identity $\mathbb{E}\{|X+Y|^2\} = \mathbb{E}\{|X|^2\} + \mathbb{E}\{|Y|^2\}$ that holds true for zero-mean and uncorrelated random variables. The expectation $T_{u21}$ can be simplified as follows

$$
\begin{align*}
T_{u21} &= \sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} \sum_{m=1}^{M} c_{mk}^2 \mathbb{E}\{|u_{mk''}|^2|u_{mk'}|^2\} + \\
&\quad \sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} \sum_{m'=1, m' \neq m}^{M} \sum_{m=1}^{M} c_{mk} c_{m'k} \mathbb{E}\{u_{mk''}^* u_{mk'} u_{m'k''} u_{m'k'}^*\} \\
&= \sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} \sum_{m=1}^{M} c_{mk}^2 \delta_{mk''} \delta_{mk'} + \\
&\quad \sum_{k'' \in \mathcal{P}_k \setminus \{k'\}} \sum_{m=1}^{M} \sum_{m'=1}^{M} c_{mk} c_{m'k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k''}).
\end{align*}
$$

(106)

Similarly, the expectation $T_{u22}$ in (105) can be simplified as
follows

$$
\begin{align}
T_{u22} &= \sum_{m=1}^{M} \sum_{m'=1}^{M} c_{mk} c_{m'k} \mathbb{E}\{|u_{mk'}|^2|u_{m'k'}|^2\} \nonumber \\
&= \sum_{m=1}^{M} c_{mk}^2 \mathbb{E}\{|u_{mk'}|^4\} + \sum_{m=1}^{M} \sum_{m'=1, m' \neq m} \sum_{m=1}^{M} c_{mk} c_{m'k} \mathbb{E}\{|u_{mk'}|^2|u_{m'k'}|^2\} \nonumber \\
&= 2 \sum_{m=1}^{M} c_{mk}^2 \delta_{mk'}^2 + 2 \sum_{m=1}^{M} c_{mk}^2 \text{tr}(\mathbf{\Theta}_{mk'}^2) + \sum_{m=1}^{M} \sum_{m'=1, m' \neq m} \sum_{m=1}^{M} c_{mk} \times \nonumber \\
&\quad c_{m'k} \delta_{mk'} \delta_{m'k'} + \sum_{m=1}^{M} \sum_{m'=1, m' \neq m} \sum_{m=1}^{M} c_{mk} c_{m'k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k'}) \nonumber \\
&= \sum_{m=1}^{M} c_{mk}^2 \delta_{mk'}^2 + \sum_{m=1}^{M} c_{mk}^2 \text{tr}(\mathbf{\Theta}_{mk'}^2) + \left(\sum_{m=1}^{M} c_{mk} \delta_{mk'}\right)^2 \nonumber \\
&\quad + \sum_{m=1}^{M} \sum_{m'=1}^{M} c_{mk} c_{m'k} \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k'}). \tag{107}
\end{align}
$$

By inserting (106) and (107) into (105), and with the aid of
some algebraic manipulations, we obtain

$$
\begin{align}
& \rho_u \eta_{k'} \sum_{m=1}^{M} c_{mk}^2 \delta_{mk'} + \rho_u \eta_{k'} p_{\tau p} \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} c_{mk}^2 \delta_{mk''} \delta_{mk'} \nonumber \\
&= \rho_u \eta_{k'} \sum_{m=1}^{M} c_{mk}^2 \delta_{mk'} \left(1 + p_{\tau p} \sum_{k'' \in \mathcal{P}_k} \delta_{m''k}\right) \nonumber \\
&\stackrel{(a)}{=} \rho_u \eta_{k'} \sum_{m=1}^{M} c_{mk} \sqrt{p_{\tau p}} \delta_{mk} \delta_{mk'} \stackrel{(b)}{=} \rho_u \eta_{k'} \sum_{m=1}^{M} \gamma_{mk} \delta_{mk'}, \tag{108}
\end{align}
$$

where (a) is obtained by using (16) and (b) by using (17).
Therefore, the total mutual interference between the user $k'$
and the user $k$ who share the same pilot sequence can be
written as follows

$$
\begin{align}
\sum_{k' \in \mathcal{P}_k \setminus \{k\}} \mathbb{E}\{|\mathrm{UI}_{uk'k'}|^2\} &= \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \eta_{k'} \gamma_{mk} \delta_{mk'} \nonumber \\
&\quad + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \eta_{k'} c_{mk}^2 \text{tr}(\mathbf{\Theta}_{mk'}^2) \tag{109} \\
&\quad + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} c_{mk} c_{m'k} \times \nonumber \\
&\quad \text{tr}(\mathbf{\Theta}_{mk'} \mathbf{\Theta}_{m'k''}) + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \eta_{k'} \left(\sum_{m=1}^{M} c_{mk} \delta_{mk'}\right)^2. \nonumber
\end{align}
$$

Let us denote $l_{uk} = \sum_{k'=1, k' \neq k}^{K} \mathbb{E}\{|\mathrm{UI}_{uk'k'}|^2\}$. By combining (104) and (109), the mutual interference at the user $k$ can
be formulated in a closed-form expression as follows

$$
\begin{align}
I_{uk} = \rho_u & \sum_{k'=1, k' \neq k}^{K} \sum_{m=1}^{M} \eta_{k'} \gamma_{mk} \delta_{mk'} + p \tau_p \rho_u \times \nonumber \\
& \sum_{k'=1, k' \neq k}^{K} \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} \eta_{k'} c_{mk} c_{m'k} \text{tr}(\Theta_{mk'} \Theta_{m'k''}) + \nonumber \\
& + p \tau_p \rho_u \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \eta_{k'} c_{mk}^2 \text{tr}(\Theta_{mk'}^2) + p \tau_p \rho_u \times \nonumber \\
& \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \eta_{k'} \left( \sum_{m=1}^{M} c_{mk} \delta_{mk'} \right)^2 . \tag{110}
\end{align}
$$

Finally, the expectation of the additive noise after MR pro-
cessing can be written as follows

$$
\begin{align}
\mathbb{E}\{|\text{NO}_{uk}|^2\} &= \sum_{m=1}^{M} \mathbb{E}\{|\hat{u}_{mk}^* w_{um}|^2\} \tag{111} \\
&= \sum_{m=1}^{M} \mathbb{E}\{|\hat{u}_{mk}|^2\} \mathbb{E}\{|w_{um}|^2\} = \sum_{m=1}^{M} \gamma_{mk}, \nonumber
\end{align}
$$

thanks to the independence between the channel estimate and
the noise. The proof follows by inserting (87), (91), (109), and
(111) into the SINR in (37) with the aid of some algebraic
manipulations.

F. *Proof of Theorem 2*

Consider the downlink SINR in (56). Thanks to the uncorre-
lation between the channel estimate and the channel estimation
error, the numerator simplifies to

$$
|DS_{dk}|^2 = \rho_d \left| \sum_{m=1}^{M} \sqrt{\eta_{mk}} \mathbb{E}\{|\hat{u}_{mk}|^2\} \right|^2 = \rho_d \left( \sum_{m=1}^{M} \sqrt{\eta_{mk}} \gamma_{mk} \right)^2 . \quad (112)
$$

The beamforming uncertainty term in the denominator of (56)
can be simplified by using (15), as follows

$$
\begin{align}
\mathbb{E}\{|\mathbf{BU}_{dk}|^2\} &= \rho_d \mathbb{E}\left\{\left|\sum_{m=1}^{M} o_{dkm}\right|^2\right\} \nonumber \\
&= \rho_d \sum_{m=1}^{M} \underbrace{\sum_{m'=1, m' \neq m}^{M} \mathbb{E}\{o_{dkm} o_{dm'k}^*\}}_{=\hat{T}_{d0}} + \rho_d \sum_{m=1}^{M} \underbrace{\sum_{m=1}^{M} \mathbb{E}\{|o_{dkm}|^2\}}, \nonumber \\
&= \hat{T}_{d1} \tag{113}
\end{align}
$$

where $o_{dkm} = \sqrt{\eta_{mk}}u_{mk}u_{mk}^* - \sqrt{\eta_{mk}}\mathbb{E}\{u_{mk}u_{mk}^*\}$ with $\omega_{mk} = \eta_{mk}$ (see Corollary 1 for $\omega_{mk}$). The closed-form expression of $\bar{T}_{d0}$ in (113) can be formulated as follows

$$
\begin{align}
\bar{T}_{d0} ={}& \nonumber \\
& p \tau_p \rho_d \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1, m' \neq m}^{M} \sqrt{\eta_{mk} \eta_{m'k} c_{mk} c_{m'k} \operatorname{tr}(\Theta_{mk} \Theta_{m'k''})}, \tag{114}
\end{align}
$$

by utilizing (22) in Corollary 1. The expectation $\bar{T}_{d1}$ in (114)
can be rewritten as follows

$$
\begin{align}
\bar{T}_{d1} ={}& \rho_d \sum_{m=1}^{M} \eta_{mk} \mathbb{E}\{|u_{mk} \hat{u}_{mk}^*|^2\} - \rho_d \sum_{m=1}^{M} \eta_{mk} |\mathbb{E}\{|\hat{u}_{mk}|^2\}|^2 \nonumber \\
& \underbrace{=}_{T_{d1}} T_{d1} \nonumber \\
& = T_{d1} - \rho_d \sum_{m=1}^{M} \eta_{mk} \gamma_{mk}^2, \tag{115}
\end{align}
$$

where we have used the identities $\mathbb{E}\{|X-\mathbb{E}\{X\}|^2\} = \mathbb{E}\{|X|^2\}-$
$|\mathbb{E}\{X\}|^2$ and $\mathbb{E}\{|X+Y|^2\} = \mathbb{E}\{|X|^2\} + \mathbb{E}\{|Y|^2\}$ for zero-mean
uncorrelated random variables. By using the identities in (91)
and (95), $T_{d1}$ in (115) can be formulated as follows

$$
\begin{align}
T_{d1} ={}& p \tau_p \rho_d \sum_{m=1}^{M} \eta_{mk} c_{mk}^2 \operatorname{tr}(\Theta_{mk}^2) + \rho_d \sum_{m=1}^{M} \eta_{mk} \gamma_{mk}^2 + \rho_d \times \nonumber \\
& \sum_{m=1}^{M} \eta_{mk} \gamma_{mk} \delta_{mk} + p \tau_p \rho_d \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \eta_{mk} c_{mk}^2 \operatorname{tr}(\Theta_{mk} \Theta_{mk''}), \tag{116}
\end{align}
$$

by using (5). Inserting (114) and (116) into (113), we obtain

$$
\begin{align}
\mathbb{E}\{|\mathbf{BU}_{dk}|^2\} ={}& \rho_d \sum_{m=1}^{M} \eta_{mk} \gamma_{mk} \delta_{mk} + p \tau_p \rho_d \times \nonumber \\
& \sum_{k'' \in \mathcal{P}_k} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk} \eta_{m'k} c_{mk} c_{m'k} \operatorname{tr}(\Theta_{mk} \Theta_{m'k''})} \tag{117} \\
& + p \tau_p \rho_d \sum_{m=1}^{M} \eta_{mk} c_{mk}^2 \operatorname{tr}(\Theta_{mk}^2). \nonumber
\end{align}
$$

The mutual interference term in the denominator of (56) can
be rewritten as follows

$$
\begin{align}
\sum_{k'=1, k' \neq k}^{K} \mathbb{E}\{|\mathrm{UI}_{dk'k}|^2\} ={}& \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \underbrace{\rho_d \mathbb{E}\{|\mathrm{UI}_{dk'k}|^2\} +}_{=T_d} \nonumber \\
& \rho_d \sum_{k' \notin \mathcal{P}_k} \sum_{m=1}^{M} \eta_{mk'} \gamma_{mk'} \delta_{mk} + p \tau_p \rho_d \times \nonumber \\
& \sum_{k' \notin \mathcal{P}_k} \sum_{k'' \in \mathcal{P}_{k'}} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \operatorname{tr}(\Theta_{mk} \Theta_{m'k''})}, \tag{118}
\end{align}
$$

where the first term on the right-hand side of (118) is obtained
by exploiting the orthogonality of the pilot sequences. In
the second summation of (118), $T_d = \mathbb{E}\{|\mathrm{UI}_{dkk'}|^2\}$ can be
rewritten as follows

$$
\begin{align}
T_d &= \rho_d \mathbb{E} \left\{ \left| \sum_{m=1}^{M} \sqrt{\eta_{mk'}} u_{mk} \hat{u}_{mk'}^* \right|^2 \right\} \nonumber \\
&= \rho_d \mathbb{E} \left\{ \left| \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} u_{mk} y_{pmk'}^*} \right|^2 \right\} \nonumber \\
&= \rho_d \mathbb{E} \left\{ \left| \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} u_{mk} \left( \sum_{k'' \in \mathcal{P}_{k'}} \sqrt{p \tau_p} u_{mk''}^* + w_{pmk'}^* \right)} \right|^2 \right\} \nonumber \\
&= \rho_d \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \mathbb{E} \{|u_{mk} w_{pmk'}^*|^2\} \nonumber \\
&\quad + \rho_d p \tau_p \mathbb{E} \left\{ \left| \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} u_{mk} \left( \sum_{k'' \in \mathcal{P}_{k'}} u_{mk''}^* \right)} \right|^2 \right\} + \nonumber \\
&\quad \underbrace{\rho_d p \tau_p \mathbb{E} \left\{ \left| \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} |u_{mk}|^2} \right|^2 \right\}}_{=T_{d21}} \nonumber \\
&= \rho_d \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \delta_{mk} + p \tau_p \rho_d (T_{d21} + T_{d22}), \tag{119}
\end{align}
$$

which is obtained by taking into account that the noise is
circularly symmetric. The term $T_{d21}$ depends on the non-
coherent interference and can be simplified as follows

$$
\begin{align}
T_{d21} &= \sum_{k'' \in \mathcal{P}_{k'}} \sum_{\{k\} m=1}^{M} \eta_{mk'} c_{mk'}^2 \mathbb{E}\{|u_{mk} u_{mk''}^*|^2\} \nonumber \\
&\quad + \sum_{k'' \in \mathcal{P}_{k'}} \sum_{m=1}^{M} \sum_{m'=1, m' \neq m}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'}} \nonumber \\
&\quad \times \mathbb{E}\{u_{mk''}^* u_{mk} u_{m'k}^* u_{m'k''}\} \nonumber \\
&= \sum_{k'' \in \mathcal{P}_{k'}} \sum_{\{k\} m=1}^{M} \eta_{mk'} c_{mk'}^2 \delta_{mk} \delta_{mk''} + \nonumber \\
&\quad \sum_{k'' \in \mathcal{P}_{k'}} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk} \eta_{m'k'} c_{mk'} c_{m'k'} \text{tr}(\Theta_{mk} \Theta_{m'k'})}, \tag{120}
\end{align}
$$

by using the second moment in (7) and the uncorrelation
among the aggregated channels. The last term $T_{d22}$ in (119)

can be simplified as follows

$$
\begin{align}
T_{d22} &= \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \mathbb{E}\{|u_{mk}|^2 |u_{m'k}|^2\}} \nonumber \\
&= \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \mathbb{E}\{|u_{mk}|^4\} + \nonumber \\
&\qquad \sum_{m=1}^{M} \sum_{m'=1, m' \neq m} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \mathbb{E}\{|u_{mk}|^2 |u_{m'k}|^2\}} \tag{121} \\
&= \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \delta_{mk}^2 + \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \text{tr}(\Theta_{mk}^2) \nonumber \\
&\quad + \left(\sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} \delta_{mk}}\right)^2 \nonumber \\
&\quad + \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \text{tr}(\Theta_{mk} \Theta_{m'k})} \nonumber
\end{align}
$$

where the last equality follows from Lemma 1. By inserting (120) and (121) into (119), and by denoting $l_{dk} = \sum_{k'=1,k' \neq k}^{K} \mathbb{E}\{|U_{dk'k}|^2\}$, (118) can be rewritten as follows

$$
\begin{align}
l_{dk} ={}& \rho_d \sum_{k'=1,k' \neq k}^{K} \sum_{m=1}^{M} \eta_{mk'} \gamma_{mk'} \delta_{mk} + p \tau_p \rho_d \times \nonumber \\
& \sum_{k'=1,k' \neq k}^{K} \sum_{k'' \in \mathcal{P}_{k'}} \sum_{m=1}^{M} \sum_{m'=1}^{M} \sqrt{\eta_{mk'} \eta_{m'k'} c_{mk'} c_{m'k'} \text{tr}(\Theta_{mk}^2 \Theta_{m'k''})} \nonumber \\
& + p \tau_p \rho_d \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \sum_{m=1}^{M} \eta_{mk'} c_{mk'}^2 \text{tr}(\Theta_{mk}^2) + \nonumber \\
& p \tau_p \rho_d \sum_{k' \in \mathcal{P}_k \setminus \{k\}} \left( \sum_{m=1}^{M} \sqrt{\eta_{mk'} c_{mk'} \delta_{mk}} \right)^2, \tag{122}
\end{align}
$$

The proof follows, by inserting (112), (113), and (122) into
(56) and by using some algebraic manipulations.

REFERENCES

[1] T. V. Chien, H. Q. Ngo, S. Chatzinotas, M. D. Renzo, and B. Ottersten, "RIS and Cell-Free Massive MIMO: A marriage for harsh propagation environments," in *Proc. IEEE GLOBECOM*, 2021.

[2] M. Giordani, M. Polese, M. Mezzavilla, S. Rangan, and M. Zorzi, "Toward 6G networks: Use cases and technologies," *IEEE Commun. Mag.*, vol. 58, no. 3, pp. 55–61, 2020.

[3] T. V. Chien and E. Björnson, *Massive MIMO Communications*. Springer International Publishing, 2017, pp. 77–116.

[4] T. S. Rappaport, Y. Xing, G. R. MacCartney, A. F. Molisch, E. Melli os, and J. Zhang, “Overview of millimeter wave communications for fifth-generation (5G) wireless networks—With a focus on propagation models,” *IEEE Trans. Antennas Propag.*, vol. 65, no. 12, pp. 6213–6230, 2017.

[5] E. Björnson, L. Sanguinetti, and M. Kountouris, “Deploying dense networks for maximal energy efficiency: Small cells meet massive MIMO,” *IEEE J. Sel. Areas Commun.*, vol. 34, no. 4, pp. 832–847, 2016.

[6] T. Van Chien, T. N. Canh, E. Björnson, and E. G. Larsson, “Power control in cellular massive MIMO with varying user activity: A deep learning solution,” *IEEE Trans. Wireless Commun.*, vol. 19, no. 9, pp. 5732–5748, 2020.

[7] G. Interdonato, E. Björnson, H. Q. Ngo, P. Frenger, and E. G. Larson, “Ubiquitous cell-free massive MIMO communications,” *EURASIP Journal on Wireless Communications and Networking*, vol. 2019, no. 1, pp. 1–13, 2019.
[8] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, "Cell-free massive MIMO versus small cells," *IEEE Trans. Wireless Commun.*, vol. 16, no. 3, pp. 1834-1850, 2017.

[9] T. Van Chien, E. Björnson, and E. G. Larsson, "Joint power allocation and load balancing optimization for energy-efficient cell-free Massive MIMO networks," *IEEE Trans. Wireless Commun.*, vol. 19, no. 10, pp. 6798-6812, 2020.

[10] G. Interdonato, H. Q. Ngo, and E. G. Larsson, "Enhanced normalized conjugate beamforming for Cell-Free Massive MIMO," *IEEE Trans. Commun.*, vol. 69, no. 5, pp. 2863-2877, 2021.

[11] E. Björnson and L. Sanguinetti, "Scalable cell-free massive MIMO systems," *IEEE Trans. Commun.*, vol. 68, no. 7, pp. 4247-4261, 2020.

[12] Q. Wu and R. Zhang, "Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming," *IEEE Trans. Wireless Commun.*, vol. 18, no. 11, pp. 5394-5409, 2019.

[13] T. A. Le, T. Van Chien, and M. Di Renzo, "Robust probabilistic-constrained optimization for IRS-aided MISO communication systems," *IEEE Wireless Commun. Lett.*, vol. 10, no. 1, pp. 1-5, 2021.

[14] M. Di Renzo, A. Zappone, M. Debbah, M. S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, "Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead," *IEEE J. Sel. Areas Commun.*, vol. 38, no. 11, pp. 2450-2525, 2020.

[15] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, "Intelligent reflecting surface aided wireless communications: A tutorial," *IEEE Trans. Commun.*, vol. 69, no. 5, pp. 3313-3351, 2021.

[16] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, "Intelligent reflecting surface enhanced wireless networks: Two-timescale beamforming optimization," *IEEE Trans. Wireless Commun.*, vol. 20, no. 1, pp. 2-17, 2021.

[17] A. Zappone, M. Di Renzo, F. Shams, X. Qian, and M. Debbah, "Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments," *IEEE Trans. Wireless Commun.*, vol. 20, no. 1, pp. 126-141, 2021.

[18] A. Abrardo, D. Dardari, and M. Di Renzo, "Intelligent reflecting surfaces: Sum-rate optimization based on statistical position information," *IEEE Trans. Commun.*, vol. 69, no. 10, pp. 7121-7136, 2021.

[19] B. Zheng and R. Zhang, "Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization," *IEEE Wireless Commun. Lett.*, vol. 9, no. 4, pp. 518-522, 2020.

[20] X. Wei, D. Shen, and L. Dai, "Channel estimation for RIS assisted wireless communications: Part I-Fundamentals, solutions, and future opportunities," *IEEE Commun. Lett.*, vol. 25, no. 5, pp. 1398-1402, 2021.

[21] N. S. Perović, L.-N. Tran, M. Di Renzo, and M. F. Flanagan, "Achievable rate optimization for MIMO systems with reconfigurable intelligent surfaces," *IEEE Trans. Wireless Commun.*, vol. 20, no. 6, pp. 3865-3882, 2021.

[22] T. Zhou, K. Xu, X. Xia, W. Xie, and J. Xu, "Achievable rate optimization for aerial intelligent reflecting surface-aided Cell-Free Massive MIMO system," *IEEE Access*, vol. 9, pp. 3828-3837, 2020.

[23] Z. Zhang and L. Dai, "Capacity improvement in wideband reconfigurable intelligent surface-aided cell-free network," in *Proc. IEEE SPAWC*, 2020, pp. 1-5.

[24] M. Bashar, K. Cumanan, A. G. Burr, P. Xiao, and M. Di Renzo, "On the performance of reconfigurable intelligent surface-aided cell-free massive MIMO uplink," in *Proc. IEEE GLOBECOM*, 2020, pp. 1-6.

[25] Y. Zhang, B. Di, H. Zhang, J. Lin, Y. Li, and L. Song, "Reconfigurable intelligent surface aided cell-free MIMO communications," *IEEE Wireless Commun. Lett.*, vol. 10, no. 4, pp. 775-779, 2021.

[26] Y. Zhang, B. Di, H. Zhang, J. Lin, C. Xu, D. Zhang, Y. Li, and L. Song, "Beyond cell-free MIMO: Energy efficient reconfigurable intelligent surface aided cell-free MIMO communications," *IEEE Trans. on Cog. Comm. and Net.*, vol. 7, no. 2, pp. 412-426, 2021.

[27] T. Van Chien, L. T. Tu, S. Chatzinotas, and B. Ottersten, "Coverage probability and ergodic capacity of intelligent reflecting surface-enhanced communication systems," *IEEE Commun. Lett.*, vol. 25, no. 1, pp. 69-73, 2021.

[28] T. Hou, Y. Liu, Z. Song, X. Sun, Y. Chen, and L. Hanzo, "Reconfigurable intelligent surface aided NOMA networks," *IEEE J. Sel. Areas Commun.*, vol. 38, no. 11, pp. 2575-2588, 2020.

[29] T. Van Chien, A. K. Papazafeiropoulos, L. T. Tu, R. Chopra, S. Chatzinotas, and B. Ottersten, "Outage probability analysis of IRS-assisted systems under spatially correlated channels," *IEEE Wireless Commun. Lett.*, vol. 10, no. 8, pp. 1815-1819, 2021.

[30] M. Di Renzo, M. Debbah, D.-T. Phan-Huy, A. Zappone, M.-S. Alouini, C. Yuen, V. Sciancalepore, G. C. Alexandropoulos, J. Hoydis, H. Gacanin et al., "Smart radio environments empowered by reconfigurable AI meta-surfaces: An idea whose time has come," *EURASIP Journal on Wireless Communications and Networking*, vol. 2019, no. 1, pp. 1-20, 2019.

[31] H. Alwazani, A. Kammoun, A. Chaaban, M. Debbah, and M.-S. Alouini, "Intelligent reflecting surface-assisted multi-user MISO communication: Channel estimation and beamforming design," *IEEE Open Jour. of the Commun. Soc.*, vol. 1, pp. 661-680, 2020.

[32] E. Björnson and L. Sanguinetti, "Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces," *IEEE Wireless Commun. Lett.*, vol. 10, no. 4, pp. 830-834, 2021.

[33] B. Clerckx, G. Kim, and S. Kim, "Correlated fading in broadcast MIMO channels: Curse or blessing?" in *Proc. IEEE GLOBECOM*, 2008.

[34] S. Kay, *Fundamentals of Statistical Signal Processing: Estimation Theory*. Prentice Hall, 1993.

[35] T. C. Mai, H. Q. Ngo, M. Egan, and T. Q. Duong, "Pilot power control for Cell-Free Massive MIMO," *IEEE Trans. Veh. Technol.*, vol. 67, no. 11, pp. 11264-11268, 2018.

[36] T. Van Chien, E. Björnson, and E. G. Larsson, "Joint pilot design and uplink power allocation in multi-cell Massive MIMO systems," *IEEE Trans. Wireless Commun.*, vol. 17, no. 3, pp. 2000-2015, 2018.

[37] H. Cramér, *Random variables and probability distributions*. Cambridge University Press, 2004, vol. 36.

[38] T. Van Chien, H. Q. Ngo, S. Chatzinotas, and B. Ottersten, "Reconfigurable intelligent surface-assisted Massive MIMO: Favorable propagation, channel hardening, and rank deficiency," *IEEE Signal Process. Mag.*, 2021, accepted for publication. [Online]. Available: arXiv preprint arXiv:2107.03434

[39] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, *Fundamentals of Massive MIMO*. Cambridge University Press, 2016.

[40] G. Gradoni and M. Di Renzo, "End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances," *IEEE Wireless Commun. Lett.*, vol. 10, no. 5, pp. 938-942, 2021.

[41] X. Qian and M. Di Renzo, "Mutual coupling and unit cell aware optimization for reconfigurable intelligent surfaces," *IEEE Wireless Commun. Lett.*, vol. 10, no. 6, pp. 1183-1187, 2021.

[42] G. Caire, "On the ergodic rate lower bounds with applications to Massive MIMO," *IEEE Trans. Wireless Commun.*, vol. 17, no. 5, pp. 3258-3268, 2018.

[43] G. Interdonato, H. Q. Ngo, P. Frenger, and E. G. Larsson, "Downlink training in cell-free Massive MIMO: A blessing in disguise," *IEEE Trans. Wireless Commun.*, vol. 18, no. 11, pp. 5153-5169, 2019.

[44] E. Björnson, J. Hoydis, and L. Sanguinetti, "Massive MIMO networks: Spectral, energy, and hardware efficiency," *Foundations and Trends in Signal Processing*, vol. 11, no. 3-4, pp. 154-655, 2017.

[45] T. V. Chien and H. Q. Ngo, *Massive MIMO Channels*. IET Publishers, 2020, ch. 11.

<figure><img src="image_15.png" /></figure>

Trinh Van Chien (S'16-M'20) received the B.S. degree in Electronics and Telecommunications from Hanoi University of Science and Technology (HUST), Vietnam, in 2012. He then received the M.S. degree in Electrical and Computer Engineering from Sungkyunkwan University (SKKU), Korea, in 2014 and the Ph.D. degree in Communication Systems from Linköping University (LiU), Sweden, in 2020. He is now a research associate at University of Luxembourg. His interest lies in convex optimization problems and machine learning applications for

wireless communications and image & video processing. He was an IEEE wireless communications letters exemplary reviewer for 2016 and 2017. He also received the award of scientific excellence in the first year of the 5G wireless project funded by European Union Horizon's 2020.
<figure><img src="image_16.png" /></figure>

**Hien Quoc Ngo** received the B.S. degree in electrical engineering from the Ho Chi Minh City University of Technology, Vietnam, in 2007, the M.S. degree in electronics and radio engineering from Kyung Hee University, South Korea, in 2010, and the Ph.D. degree in communication systems from Linköping University (LiU), Sweden, in 2015. In 2014, he visited the Nokia Bell Labs, Murray Hill, New Jersey, USA. From January 2016 to April 2017, Hien Quoc Ngo was a VR researcher at the Department of Electrical Engineering (ISY), LiU.

He was also a Visiting Research Fellow at the School of Electronics, Electrical Engineering and Computer Science, Queen's University Belfast, UK, funded by the Swedish Research Council.

Hien Quoc Ngo is currently a Reader (Associate Professor) at Queen's University Belfast, UK. His main research interests include massive (large-scale) MIMO systems, cell-free massive MIMO, physical layer security, and cooperative communications. He has co-authored many research papers in wireless communications and co-authored the Cambridge University Press textbook *Fundamentals of Massive MIMO* (2016).

Dr. Hien Quoc Ngo received the IEEE Com Soc Stephen O. Rice Prize in Communications Theory in 2015, the IEEE Com Soc Leonard G. Abraham Prize in 2017, and the Best PhD Award from EURASIP in 2018. He also received the IEEE Sweden VT-COM-IT Joint Chapter Best Student Journal Paper Award in 2015. He was an *IEEE Communications Letters* exemplary reviewer for 2014, an *IEEE Transactions on Communications* exemplary reviewer for 2015, and an *IEEE Wireless Communications Letters* exemplary reviewer for 2016. He was awarded the UKRI Future Leaders Fellowship in 2019. Dr. Hien Quoc Ngo currently serves as an Editor for the IEEE Transactions on Wireless Communications, IEEE Wireless Communications Letters, Digital Signal Processing, Elsevier Physical Communication (PHYCOM), and IEICE Transactions on Fundamentals of Electronics, Communications and Computer Sciences. He was a Guest Editor of IET Communications, special issue on “Recent Advances on 5G Communications” and a Guest Editor of IEEE Access, special issue on “Modelling, Analysis, and Design of 5G Ultra-Dense Networks”, in 2017. He has been a member of Technical Program Committees for several IEEE conferences such as ICC, GLOBECOM, WCNC, and VTC.

<figure><img src="image_17.png" /></figure>

**Marco Di Renzo** (Fellow, IEEE) received the Laurea (cum laude) and Ph.D. degrees in electrical engineering from the University of L'Aquila, Italy, in 2003 and 2007, respectively, and the Habilitation à Diriger des Recherches (Doctor of Science) degree from University Paris-Sud (now Paris-Saclay University), France, in 2013. Since 2010, he has been with the French National Center for Scientific Research (CNRS), where he is a CNRS Research Director (Professor) with the Laboratory of Signals and Systems (L2S) of Paris-Saclay University –

CNRS and CentraleSupelec, Paris, France. In Paris-Saclay University, he serves as the Coordinator of the Communications and Networks Research Area of the Laboratory of Excellence DigiCosme, and as a Member of the Admission and Evaluation Committee of the Ph.D. School on Information and Communication Technologies. He is the Editor-in-Chief of IEEE Communications Letters and a Distinguished Speaker of the IEEE Vehicular Technology Society. In 2017-2020, he was a Distinguished Lecturer of the IEEE Vehicular Technology Society and IEEE Communications Society. He has received several research distinctions, which include the SEE-IEEE Alain Glavieux Award, the IEEE Jack Neubauer Memorial Best Systems Paper Award, the Royal Academy of Engineering Distinguished Visiting Fellowship, the Nokia Foundation Visiting Professorship, the Fulbright Fellowship, and the 2021 EURASIP Journal on Wireless Communications and Networking Best Paper Award. He is a Fellow of the UK Institution of Engineering and Technology (IET), a Fellow of the Asia-Pacific Artificial Intelligence Association (AAIA), an Ordinary Member of the European Academy of Sciences and Arts (EASA), and an Ordinary Member of the Academia Europaea (AE). Also, he is a Highly Cited Researcher.

<figure><img src="image_18.png" /></figure>

**Symeon Chatzinotas** is Full Professor and Head of the SIGCOM Research Group at SnT, University of Luxembourg. He is coordinating the research activities on communications and networking, acting as a PI for more than 20 projects and main representative for 3GPP, ETSI, DVB. He is currently serving in the editorial board of the IEEE Transactions on Communications, IEEE Open Journal of Vehicular Technology and the International Journal of Satellite Communications and Networking.

In the past, he has been a Visiting Professor at the University of Parma, Italy and was involved in numerous R&D projects for NCSR Demokritos, CERTH Hellas and CCSR, University of Surrey.

He was the co-recipient of the 2014 IEEE Distinguished Contributions to Satellite Communications Award and Best Paper Awards at EURASIP WJCN, CROWNCOM, ICSSC. He has (co-)authored more than 500 technical papers in refereed international journals, conferences and scientific books.

<figure><img src="image_19.png" /></figure>

**Björn Ottersten** (S'87–M'89–SM'99–F'04) received the M.S. degree in electrical engineering and applied physics from Linköping University, Linköping, Sweden, in 1986, and the Ph.D. degree in electrical engineering from Stanford University, Stanford, CA, USA, in 1990. He has held research positions with the Department of Electrical Engineering, Linköping University, the Information Systems Laboratory, Stanford University, the Katholieke Universiteit Leuven, Leuven, Belgium, and the University of Luxembourg, Luxembourg. From 1996 to

1997, he was the Director of Research with ArrayComm, Inc., a start-up in San Jose, CA, USA, based on his patented technology. In 1991, he was appointed Professor of signal processing with the Royal Institute of Technology (KTH), Stockholm, Sweden. Dr. Ottersten has been Head of the Department for Signals, Sensors, and Systems, KTH, and Dean of the School of Electrical Engineering, KTH. He is currently the Director for the Interdisciplinary Centre for Security, Reliability and Trust, University of Luxembourg. He is a recipient of the IEEE Signal Processing Society Technical Achievement Award, the EURASIP Group Technical Achievement Award, and the European Research Council advanced research grant twice. He has co-authored journal papers that received the IEEE Signal Processing Society Best Paper Award in 1993, 2001, 2006, 2013, and 2019, and 8 IEEE conference papers best paper awards. He has been a board member of IEEE Signal Processing Society, the Swedish Research Council and currently serves of the boards of EURASIP and the Swedish Foundation for Strategic Research. Dr. Ottersten has served as Editor in Chief of EURASIP Signal Processing, and acted on the editorial boards of IEEE Transactions on Signal Processing, IEEE Signal Processing Magazine, IEEE Open Journal for Signal Processing, EURASIP Journal of Advances in Signal Processing and Foundations and Trends in Signal Processing. He is a fellow of EURASIP.