Impact of Channel Aging on Reconfigurable Intelligent
Surface Aided Massive MIMO Systems with Statistical CSI
Anastasios Papazafeiropoulos, Ioannis Krikidis, Pandelis Kourtessis A. Papazafeiropoulos is with the Communications and Intelligent Systems Research Group, University of Hertfordshire, Hatfield AL10 9AB, U. K., and with the SnT at the University of Luxembourg, Luxembourg. I. Krikidis is with the IRIDA Research Centre for Communication Technologies, Department of Electrical and Computer Engineering, University of Cyprus, Cyprus. P. Kourtessis is with the Communications and Intelligent Systems Research Group, University of Hertfordshire, Hatfield AL10 9AB, U. K. A. Papazafeiropoulos was supported  by the University of Hertfordshire's 5-year Vice Chancellor's Research Fellowship. Also, this work was co-funded by the European Regional Development Fund and the Republic of Cyprus through the Research and Innovation Foundation under the project INFRASTRUCTURES/1216/0017 (IRIDA). It has also received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 819819).
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
The incorporation of reconfigurable intelligent surface (RIS) into massive multiple-input-multiple-output (mMIMO) systems can unleash the potential of next-generation networks by improving the performance of  user equipments (UEs) in service dead zones. However, their requirement for accurate channel state information (CSI) is critical, and especially,  applications with UE mobility  that induce channel aging make challenging the achievement of adequate quality of service. Hence, in this work, we investigate the impact of channel aging on the performance of RIS-assisted mMIMO systems under both spatial correlation and imperfect CSI conditions. Specifically, by accounting for channel aging during both uplink training and downlink data transmission phases, we first perform  minimum mean square error (MMSE) channel estimation to obtain the UE effective channels with low overhead similar to conventional systems without RIS. Next, we derive the downlink achievable sum spectral efficiency (SE) with  regularized zero-forcing (RZF) precoding in closed-form  being dependent only on large-scale statistics by using the deterministic equivalent (DE) analysis. Subsequently, we present the attractive optimization of the achievable sum SE with respect to the phase shifts and the total transmit power that can be performed every several coherence intervals due to the slow variation of the large-scale statistics. Numerical results validate the analytical expressions and demonstrate the performance while allowing the extraction of insightful design conclusions for common scenarios including UE mobility. In particular,  channel aging degrades the performance but its impact can be controlled by choosing appropriately the frame duration or by increasing the number of RIS elements.

Reconfigurable intelligent surface (RIS), channel aging, channel estimation, achievable spectral efficiency, beyond 5G networks.

§ INTRODUCTION

Emerging applications bring challenges that demand ever-higher data rates and increased connectivity/coverage together with ultra-reliable and low-latency wireless communication (URLLC) requirements. Especially, these applications have led to the development of  disruptive technologies such as massive multiple-input multiple-output (mMIMO) systems and millimeter-wave (mmWave) communications [1].  Unfortunately, existing techniques incur additional power and hardware costs while they cannot guarantee an adequate quality of service (QoS) in dead zones due to obstacles. For example, mMIMO exhibits poor performance in low scattering conditions, and  the large number of active elements might result in prohibitive energy usage. In particular, they focus on improvements regarding the transmission and reception, while the wireless propagation environment is left uncontrollable.   Furthermore, the time-varying and random nature of the wireless channel constitute the ultimate impediment to achieving the URLLC and rate targets.

In this direction, sixth generation (6G) networks, aiming at covering the higher rate demands and more stringent constraints, have appeared with the reconfigurable intelligent surface (RIS) being among its proposed promising technologies. Actually, RIS   has attracted significant attention since it overcomes the aforementioned issues [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. Specifically, a RIS is a software-defined surface that is usually attached  to existing infrastructure to alleviate blockage effects. It consists of a large number of  individually-controlled, low-cost, and nearly passive elements. A RIS achieves to adapt to changes in the propagation environment and modify the radio waves since each of its elements can induce an adjustable phase shift to each incident signal, which enables a dynamic control over the wireless propagation channel. For instance, in [3], a minimization of the transmit power at the base station (BS) with signal-to-interference-plus-noise ratio (SINR) constraints took place in a RIS-assisted multi-user (MU) multiple-input single-output (MISO) communication system by  jointly optimizing the precoding and reflecting beamforming matrices (RBMs). In [4], the sum rate was maximized subject to a transmit power constraint, while in [5], the sum rate was optimized by accounting also for correlated Rayleigh fading and inevitable hardware impairments at both the transceiver and the RIS. Similarly, in [6], the maximization of the minimum UE rate was studied in the case of a large number of antennas,  in [7], the impact of hardware impairments was evaluated, and in [11], the impact of imperfect CSI on the outage probability was investigated. Notably, given that RIS-assisted MIMO systems, having a reduced number of active  radio frequency (RF) chains,  can achieve similar performance to mMIMO without RIS [13], it is indicated that a more cost and energy-efficient implementation of mMIMO is possible.

To reap the benefits of RIS and mMIMO and arrive at realistic conclusions, the acquisition of accurate channel state information (CSI) is of paramount importance [14, 15, 8, 16, 17, 18]. [Note that many previous works assumed perfect CSI, which is a highly unrealistic assumption.] However, channel estimation (CE) in RIS-aided systems is quite challenging because of two main reasons. First, although its passive elements render RIS energy-efficient, they make infeasible conventional CE through transmitting and receiving pilots, which require active elements. Second, RIS generally includes a large number of elements, which require a prohibitively high training overhead that severely reduces  the achievable rate. For instance, in  [14], an ON/OFF CE scheme was proposed, where   the least-squares  estimates of all RIS-assisted MISO channels   with  a single user were calculated one by one. Moreover, other works  such as [15, 8] do not provide analytical expressions for the estimated channel that could be exploited for the derivation of the spectral efficiency (SE).  In [16], all RIS elements were assumed active during training but a number of sub-phases equal at least to the number of RIS elements are required, which results in a lower rate because the overhead on the coherence time for CE is larger. Moreover, this method does not provide the covariance of the estimated channel vector from all RIS elements to a specific UE but estimates of the individual channels while leaving the correlation among them unknown. An  effective method with low overhead compared to previous works is to estimate the cascaded BS-RIS-UE channel  in a single phase based on  minimum mean square error (MMSE) as in [5].  Notably, another method, which reduces the overhead has been presented by exploiting RIS partitioning into subgroups, e.g., see [18]. Also, therein, an insightful categorization of the various CE approaches concerning RIS-assisted systems has been provided.

In practice, CSI is not only imperfect but can also be outdated because of channel aging [19, 20, 21, 22, 23]. The cause of channel aging is the UE mobility, which renders the channel time-varying, i.e., contrary to the standard block fading model, the channel evolves with time and is different during each symbol. Thus, a mismatch appears between the current channel and the estimated channel used for detection or precoding. Interestingly, in  [22], channel aging was also considered during the training phase for a more realistic study. Several works have studied the impact of channel aging in mMIMO systems as mentioned but little attention has been given to its effect in  RIS-assisted systems despite its great significance [23].

In principle, the phase shifts  optimization lies on two methodologies with respect to CSI, namely, instantaneous CSI (I-CSI) [3, 4] and statistical CSI (S-CSI) [24, 6, 5, 25, 26, 23, 27, 28, 29]. The first approach suggests the optimization of the phases at every coherence interval because the related expressions depend on small-scale fading, while the second approach concerns expressions that depend on  large-scale statistics, which vary every several coherence intervals. Hence, the latter approach enables considerably the reduction of the signal overhead and the computational complexity, which can become excessively high in the case of a large number of RIS elements and BS antennas. Moreover, for the same reasons, the S-CSI approach is more energy-efficient. Notably, in high mobility scenarios, which are faster time-varying, the I-CSI method would be very challenging to be implemented since the tuning of the RIS parameters should be repeated very frequently. On the contrary, the application of the  S-CSI appears to be more practical.

§.§ Motivation/Contributions

Faced with these challenges, the motivation of this work is to conduct a realistic characterization of the downlink achievable sum SE of RIS-assisted mMIMO systems accounting for UE mobility  and imperfect CSI under correlated Rayleigh fading conditions, when  regularized zero-forcing (RZF) precoding is applied.

* Contrary to the majority of existing works on RIS-assisted systems such as [3, 4, 24, 6, 17, 18, 5, 25, 26, 27, 28], which  considered static UEs, we account for channel aging due to UE mobility. To the best of our knowledge, the only  previous works considering channel aging are [23] and [29]. The former focused on  mmWave communications with LoS links and with a finite number of BS antennas, while we consider mMIMO systems, correlated Rayleigh fading, and channel aging during the training phase too.  The latter did not account for correlated fading  mMIMO, and ZF precoding, while we have taken correlation into account and we have focused on mMIMO. Note that no optimization took place in [29]. In addition, we have resorted to the deterministic equivalent analysis to provide the rate for RZF. Also, despite that many previous works have relied on independent Rayleigh fading e.g., [3, 4], we consider correlated Rayleigh fading [30], which appears unavoidable in practice. Moreover, we account for S-CSI instead of I-CSI since the former is more suitable for studying time-varying channels. In particular,  compared to other works, which are based on statistical CSI such as [24, 6, 5, 25, 26, 23, 27, 28, 29], our work is the only one that has studied the impact of channel aging by taking into account correlated fading, imperfect CSI, and RZF being a more advanced precoder, which increases the difficulty for the derivation of closed-form expressions. For example, compared to [5] focusing on the uplink and maximal ratio combining (MRC), we have assumed a more suitable model for RIS correlation, the downlink, RZF, and we have focused on the impact of channel aging. Similarly, compared to [6], we have assumed imperfect CSI, have focused on the sum rate instead of the max-min rate, and have studied channel aging. Also, although channel aging has been studied in [29], no correlation has been considered, ZF instead of RZF has been applied, and no optimization has been performed.

* We introduce channel aging not only in the downlink data transmission phase but also during the uplink training phase as in [22]. In particular, based on [5], we perform MMSE  estimation and obtain the effective channel estimate that ages with time. The proposed approach provides the estimated channel with low overhead and in closed-form that enables further manipulations to derive the achievable SE. Previous works, e.g., [15, 8] do not provide analytical expressions or have other disadvantages such as high overhead and unknown inter-element correlation [16].

* Exploiting the deterministic equivalent (DE) analysis, we obtain the DE of the downlink sum SE of RIS-assisted mMIMO systems with RZF precoding under UE mobility and correlated Rayleigh fading conditions.  The DE results are of great importance because they provide closed-form expressions in terms of a convergent system of fixed-point 	equations that allow efficient optimization.

* We formulate the maximization problem regarding the sum SE with respect to RBM and total transmit power constraints. Notably, given that the sum SE depends only on large-scale statistics, the proposed optimization  can be performed every several coherence intervals, and thus, reduce significantly the signal overhead, which is large in time-varying channels.

* We verify the analytical results with Monte Carlo (MC) simulations, and we  shed light on the impact of channel aging on the downlink sum SE of a RIS-assisted mMIMO system  due to correlated fading and channel aging. For comparison, we depict results corresponding to no mobility to show the degradation due to channel aging and the inferior performance of maximum ratio transmission (MRT) precoding.

§.§ Paper Outline

The remainder of this paper is organized as follows. Section <ref> presents the system model of a RIS-assisted mMIMO system  with imperfect CSI under correlated Rayleigh fading and  channel aging conditions. Section <ref> describes the CE accounting for channel aging. Section <ref> presents the downlink sum SE, while Section <ref> provides the optimization regarding the  RBM and the transmit power.
The numerical results are discussed in Section <ref>, and Section <ref> concludes the paper.

A downlink RIS-assisted mMIMO communication system with $ M $ BS antennas, $ L $ RIS elements, and moving $ K $ UEs.

§.§ Notation
Vectors and matrices are denoted by boldface lower and upper case symbols, respectively. The notations $(\cdot)^\T$, $(\cdot)^\H$, and $( · )$ represent the transpose, Hermitian transpose, and trace operators, respectively. The expectation operator is denoted by $[·]$  while $ () $ represents an $ n×n $ diagonal matrix with diagonal elements being the elements of vector $ $.
%	In the case of a matrix $ \bA $, $\diag\left(\bA\right) $ denotes  a vector with elements the diagonal elements of $ \bA $.
Also, the notations $$ and $a_n≍b_n$ with $a_n$ and $b_n$ being two infinite sequences denote almost sure convergence as $ M →∞$. The notation $ lim_x→c  f(x) $ denotes the limit of $ f $ of $ x $ as $ x $ approaches $ c $, and the notation $  f(x)x   $ denotes the partial derivative of $ f $ with respect to  $ x $.  Finally, $∼(0̱,Σ)$ represents a circularly symmetric complex Gaussian vector with {zero mean} and covariance matrix $Σ$.
\section{System Model}\label{SystemModel}
We consider a RIS-assisted downlink mMIMO  system, where a BS, equipped with $ M $ antennas communicates with $ K $  single-antenna noncooperative UEs behind obstacles. A RIS, consisting of $ L $ passive reflecting elements is located in the LoS of the BS to assist the communication with the UEs, e.g., imagine the common scenario where   both the BS and RIS are deployed at high altitude and their locations are fixed, as shown in Fig. \ref{Fig0}. The RIS can  dynamically adjust the phase shift induced by each reflecting element on the impinging electromagnetic waves through a perfect  smart controller that is connected to  the BS in terms of a perfect backhaul link. The size of each RIS element is $ d_H×d_V $, where $d_V$ and $d_H$ express its vertical height and its horizontal width, respectively. The proposed model considers also the potential presence of direct links between the BS and the UEs. However, these could be also neglected in the cases of high penetration losses and/or big signal blockages.

\subsection{Channel Model}\label{ChannelModel}
We account for a quasi-static fading model with coherence bandwidth much larger than the channel bandwidth. We employ the standard block fading model with each coherence interval/block including $τ_c= B_cT_c$ channel uses, where $ B_c $ and $ T_c $ are the coherence bandwidth and the coherence time in $Hz$ and $s$, respectively.

Within the transmission in each coherence block and during the $ n $th time slot, let  $ _1=[_11…,_1L ] ∈ℂ^M ×L$ and $ _2k,n ∈ℂ^L ×1$ be  the LoS channel between the BS and the RIS and the channel between the RIS and UE $ k $ at the $ n $th time instant. Note that  $ _1i $ for $ i=1,…,L $ denotes the $ i $th column vector of $ _1 $.	Similarly, let $ _k,n ∈ℂ^M ×1 $ be the direct channel between the BS and UE $ k $ at the $ n $th time instant.  Despite that  the majority of existing works, e.g., [3, 4], assumed independent Rayleigh model, in practice, correlated fading appears, which affects the performance  [30].		Thus, $ _2k,n $ and  $ _k,n $ are described in terms of correlated Rayleigh fading distributions as
\begin{align}
\bh_{2k,n}&=\sqrt{\beta_{\mathrm{h}_{2},k}}\bR_{\mathrm{RIS},k}^{1/2}\bq_{2k,n},\\
\bg_{k,n}&=\sqrt{\beta_{\mathrm{g},k}}\bR_{\mathrm{BS},k}^{1/2}\bq_{k,n},
\end{align}
where $ _RIS,k ∈ℂ^L ×L $ and $ _BS,k ∈ℂ^M ×M $ express the deterministic Hermitian-symmetric positive semi-definite correlation
matrices at the RIS and the BS respectively with $ (_RIS,k )=L $ and $ (_BS,k )=M $.  The correlation matrices $ _RIS,k $ and $ _BS,k ∀k$ are assumed to be known by the network since they can be obtained by  existing estimation methods [31].\footnote{The correlation matrices and the path-losses are independent of  $n$ because these represent effects that vary with time in a much slower pace than the coherence time.} Another way of practical calculation of
the covariance matrices follows. Especially, as can be seen by the expression of the covariance matrices, they depend on
the distances $ d_BS $ and $ d_IRS $ and the angles. The distances are based on the construction of the BS and the IRS.
Moreover, the angles can be calculated when the locations are given. Moreover,  $ β_h_2,k $ and $ β_g,k $ describe the path-losses  of the RIS-UE $ k $ and BS-UE $ k $ links, respectively. Especially, $ β_g,k $ is expected to be small because of the blockages between the BS and the UEs. Also, $ _2k,n∼𝒞𝒩(0̱,_L) $ and $ _k,n ∼𝒞𝒩(0̱,_M) $ denote the corresponding fast-fading vectors at the $ n $th time instant. Note that fast fading vectors change within each coherence block, while the correlation matrices are assumed constant for a large number of coherence blocks.

%		Given that the RIS is designed to be installed in a location providing an LoS channel with the BS, the channel matrix $ \bH_{1} $ will likely have rank one, which results in performance gains only when $ K=1 $ [6]. However, an MU scenario requires $ \mathrm{rank}\left(\bH_{1}\right)\ge K $. 				The higher rank could be achieved by placing the RIS close to the BS or by 	deterministic scattering between the BS and the RIS.
The high rank LoS channel $ _1 $ is described as
\begin{align}
[\bH_{1}]_{m,l}=\sqrt{\beta_{1}} \exp\big(j \frac{2 \pi }{\lambda}\left(m-1\right)d_{\mathrm{BS}}\sin \theta_{1,l}\sin \phi_{1,l}\nn\\
+\left(l-1\right)d_{\mathrm{RIS}}\sin \theta_{2,m}\sin \phi_{2,m}\big)\!,
\end{align}
where  $ β_1 $ is the path-loss between the BS and RIS, $ λ$ is the carrier wavelength, while $ d_BS $ and $ d_RIS $ are the inter-antenna separation at the BS and inter-element separation at the RIS, respectively [16]. Also, $ θ_1,l $ and $ ϕ_1,l $ denote the elevation and azimuth LoS angles of departure (AoD) at the BS with respect to RIS element $ l $, and $ θ_2,l $ and $ ϕ_2,l $ denote the elevation and azimuth LoS angles of arrival (AoA) at the RIS. It is worthwhile to mention that  $ _1 $ can be obtained similarly to the covariance matrices
since the dependence of their expressions on the distances and the angles is similar.

The response of the $ L $ elements is described by the diagonal RBM    $=diag( μ_1e^j θ_1, …, μ_Le^j θ_L )∈ℂ^L×L$, where $ θ_l ∈[0,2π]$ and $ μ_l∈[0,1] $ are the phase  and amplitude coefficient for RIS element $ l $, respectively. Herein, we assume maximum signal reflection, i.e., $ μ_l=1  ∀l$ [3].\footnote{Recently, it was shown that the amplitude and phase responses are 	intertwined in practice, while the assumption of  independence between the amplitude and the phase shift or even a unity amplitude is unrealistic, e.g., see [32]. However, this 	assumption regarding independence still allows revealing fundamental properties of the channel aging of the proposed model, while the consideration of the phase shift model in [32] is an interesting idea for extension of the current work, i.e., to study the impact of  channel aging on RIS-assisted systems by accounting for this intertwinement.    } For the sake of exposition, the overall channel vector $ _k,n=_k,n+ _1_2k,n $, conditioned on $ $ is distributed as $ _k,n∼( 0, _k ) $, where $ _k= β_g,k_BS,k+ β_h_2,k_1 _RIS,k^_1^$. Given that $R_k$ depends on the path-losses, the correlation matrices, and $ _1 $, which are all assumed to known as explained previously, $R_k$ can also be assumed known by the network.

\begin{remark}\label{rem1}
Although, it is uncommon to meet  independent Rayleigh fading in practice [30], in such a case, we have $ {\bR}_{\mathrm{RIS},k}=\Id_{L}$. Then, the overall covariance becomes $ \bR_{k} =\beta_{\mathrm{d},k}\bR_{\mathrm{BS},k}+\beta_{2,k}\bH_{1}\bH_{1}^{\H}$. Obviously,  $ \bR_{k} $ does not depend on the RBM and cannot be optimized. Hence, the RIS cannot be exploited. 	Nevertheless, note that even under these conditions, the  RIS enhances  the communication with an additional signal to the receiver.
\end{remark}

\subsection{Channel Aging}
In practice, the relative movement between the UEs and the RIS, i.e., the UE mobility causes a phenomenon, known as channel aging~[19, 20, 21].\footnote{Normally, all RIS elements have the same relative movement comparing to a specific UE.} In particular, this movement results in a Doppler shift that makes the channel change with time. Hence, contrary to the conventional block fading channel model, the channel coefficients, exhibiting flat fading, vary from symbol to symbol. However, they are constant within one symbol.  The symbol duration is assumed smaller than or equal to the coherence time of all UEs. This assumption is common in works studying the impact of channel aging such as [19]-[22]. The channel use is denoted by $ n ∈{ 1, …, τ_c} $.

Mathematically,  the channel realization $ _k,n $ at the $ n $th time instant is modeled as a function of its
initial state $ _k,0 $ and an innovation component as [22]
\begin{align}
\bh_{k,n}=\al_{k,n}\bh_{k,0}+\bar{\al}_{k,n}\bee_{k,n},\label{GaussMarkoModel}
\end{align}
where $ _k,n∼𝒞𝒩(0̱,_k)$  denotes  the independent innovation component at the $ n $th time instant and $ _k,n= J_0(2 πf_DT_sn) $ is the temporal correlation coefficient of UE $ k $  between the channel realizations at time $ 0 $ and $ n $ with  $J_0(·)$ being the zeroth-order Bessel function of the first kind, $T_s$ being the channel sampling duration,  and $f_D=vf_c/c$ being the maximum Doppler shift.\footnote{The second-order statistics of the channel including the path-losses and $\al_{k,n}$ are estimated during the connection establishment when the BS estimates the location and the velocity of the UE. While the path-losses are estimated using the average pilot power, $ \al_{k,n} $ can be estimated using the temporal correlation among the same set of pilots.} Also, we denote $ _k,n=√(1-^2_k,n)$. Note that $v$ is the velocity of the UE, $c = 3 ×10^8$ $m/s$ is the speed of light, and $f_c$ is the carrier frequency. As can be seen, a higher UE velocity of the UE  or higher delay result in decrease of $ _k,n $ though not monotonically, since there are some ripples. It is worthwhile to mention that the model in \eqref{GaussMarkoModel}  is not  autoregressive of first-order as in previous works [19, 20] since the current channel is not determined in terms of its state   at the previous time instant,  but it depends on its state at an initial time $ _k,0 $. The advantage is that the statistics of the model exactly match with that of the Jakes’ model [22].
%	\begin{remark}
%		Generally, the Doppler shifts observed at the BS and the RIS are different due to their different locations. By writing the channel aging model for each individual channel, and adding the two signal models we result in an expression similar to \eqref{GaussMarkoModel},  In other words, \eqref{GaussMarkoModel}  describes the overall channel aging effect with $ \al_{k,n} $ expressing the sum of the temporal correlations of the individual channels.
%	\end{remark}

\section{Channel Estimation with Channel Aging}\label{ChannelEstimation}
Perfect CSI is not available in practice but the BS needs to estimate the channel. On this ground,  we consider the standard time-division-duplex (TDD) protocol, where each block consists of $τ$ channel uses for the uplink training phase and $τ_d=τ_c-τ$ channel uses for the downlink data transmission phase~[33]. The disadvantage of RIS-assisted systems is that the RIS, which consists of passive elements, cannot  process the received pilot symbols from the UEs to obtain the  estimated channels and cannot  send pilots to the BS for CE. Herein, contrary to ON/OFF channel estimation schemes such as  [14] and [16] that require $ L+1 $ phases,  we perform the CE in a single phase.  Notably, the consideration of  the cascaded channel CE instead of the individual channels has already been applied in several works such as [5, 6]. Actually, it is more beneficial to consider the overall channel because this allows computing the correlation among inter-element links, while in the case of individual channels, this correlation remains unknown. Also, our CE is accompanied by reduced feedback and allows higher achievable SE  due to the larger pre-log factor since the training overhead is much lower.

During the uplink training phase, 	each UE transmits a $τ$-length
mutually orthogonal training sequences, i.e., $ψ_k= [ψ_k,1,…,ψ_k,τ ]^∈^τ×1$. We assume that the pilot sequence consists of $ K $  pilot symbols, which is the  minimal number for channel estimation, i.e., $ τ≥K $ [34]. Also, $ $ is assumed fixed. Thus, the received signal by the BS at time $ n $ is  given by
\begin{align}
\bY_{n}=\sqrt{p_{\mathrm{p}}}\sum_{i=1}^{K}\bh_{i,n}\bpsi_{i}^{\H}+\bW_{n},\label{TrainingReceived}
\end{align}
where $ p_p≥0 $ is the common pilot transmit power for all UEs   and $ _n∈ℂ^M×τ $ is spatially white additive Gaussian noise matrix at the BS   during this phase. After correlation of the received signal with the training sequence of UE $ k $ $1/√(p_)_k$, we obtain
\begin{align}
\tilde{\by}_{k,n}=\bh_{k,n}+\frac{1}{\sqrt{p_{\rp}}}\tilde{\bw}_{k,n},\label{TrainingReceived1}
\end{align}
where $ _k,n=_n_k∼𝒞𝒩(0̱,σ̃^2_M) $.

Although this received signal can be used to estimate the channel at any time slot of the block, the estimated channel  will deteriorate as the time interval between training and transmission increases. On this ground, we consider the channel estimate at $ n=K+1 $ since the estimate will be worse at a later instant. Based on \eqref{GaussMarkoModel}, the channel at the $ n $th instant ($ n≤K $) can be described in terms of the channel at time $ K+1 $ as
\begin{align}
\bh_{k,n}=\al_{k,\zeta-n}\bh_{k,\zeta}+\bar{\al}_{k,\zeta-n}\tilde{\bee}_{k,n},\label{innov1}
\end{align}
where $  _k,n∼𝒞𝒩(0̱,_k) $ is the  independent innovation vector, which relates $ 	_k,n $ and $ _k,ζ $. Also, we have defined $ ζ=K+1 $ to simplify the notation. Inserting \eqref{innov1} into \eqref{TrainingReceived1}, we obtain
\begin{align}
\tilde{\by}_{k,n}&=\al_{k,\zeta-n}\bh_{k,\zeta}+\bar{\al}_{k,\zeta-n}\tilde{\bee}_{k,n}+\frac{1}{\sqrt{p_{\rp}}}\tilde{\bw}_{k,n}\label{TrainingReceived}.
\end{align}

By applying  the standard minimum mean square error (MMSE) estimation [33], the BS obtains the channel estimate of $ _k,ζ $ as
\begin{align}
\hat{\bh}_{k,\zeta}=\al_{k,\zeta-n}\bR_{k} \bQ			\tilde{\by}_{k,n},\label{estimatedChannel}
\end{align}
where $ =(_k+σ̃^2/p__M)^-1 $. The estimate  	$ _k,ζ $  is distributed as $ 𝒞𝒩(0̱,_k) $, where $ _k=^2_k,ζ-n_k _k$.  According to the orthogonality property of MMSE estimation, the independent channel estimation error vector is $	_k,ζ=	_k,ζ-	_k,ζ  $ and distributed as $ 𝒞𝒩(0̱,_k) $, where $ _k=_k-_k $. Notably, the channel estimate 	$ _k,ζ $  in \eqref{estimatedChannel} includes the  degradation due to the  channel aging.
\begin{remark}
In the case of no channel aging, i.e., when $ \al_{k,\zeta-n}=1 $, we reduce to the conventional block-fading model. Moreover, as can be seen, the estimation error takes values between $\bR_{k}(\Id_{M}- \al^{2}_{k,1} \bQ \bR_{k})  $ and  $\bR_{k}(\Id_{M}- \al^{2}_{k,K} \bQ \bR_{k})  $. Also, we observe that in the case of no mobility, the estimation error vanishes as the pilot signal-to-noise ratio (SNR) $ \gamma=\frac{p_{\rp}}{\tilde{\sigma}^{2}} $ increases, while it saturates as $\displaystyle\lim_{\gamma\to \infty}\bPsi_{k}= (1-\al^{2}_{k,\zeta-n})\bR_{k} $ in the case of channel aging. The latter shows that the estimation error increases as the UE moves with higher velocity and as the number of UEs increases since $ \zeta  $ increases.
\end{remark}

In Sec. \ref{Numerical}, we illustrate the normalized mean square error (NMSE) defined as
\begin{align}
\mathrm{NMSE}_{k}&=\frac{\tr\big(\EE[(\hat{\bh}_{k}-{\bh}_{k})(\hat{\bh}_{k}-{\bh}_{k})^{\H}]\big)}{\tr\left(\EE[{\bh}_{k}{\bh}_{k}^{\H}]\right)}\\
\end{align}

According to \eqref{nmse1},  an increase in channel aging results in the increase of the $ NMSE_k $. Overall, channel aging has a detrimental on channel estimation. Below, we elaborate on its impact during the downlink transmission.

\section{Downlink Transmission}\label{PerformanceAnalysis}
The downlink transmission of data from the BS to all UEs consists of a broadcast channel that has to make use of a certain precoding strategy in terms of a precoding vector $_k,n ∈^M ×1$. In parallel, taking advantage of TDD and its channel reciprocity, the downlink channel is the Hermitian transpose of the uplink channel. Thus, the received signal $r_k,n∈$ by UE $ k $ during the data transmission phase ($n=K+1,…, τ_c$)  can be written as
\begin{align}
\end{align}
where  $_n=∑_i=1^K√(p_i)_i,nx_i,n$ describes the transmit signal vector  by the BS,  $p_i≥0$ is the transmit power to UE $ i $, and $z_k,n ∼(0,σ^2)$ is complex Gaussian noise at UE $k$. Note that   $_i,n ∈^M ×K$ and
% $x_{i,n}\in \sim \cC\cN(0,1)$
$ x_i,n $ are  the linear precoding vector  and the data symbol with $ {|x_i,n|^2}=1 $, respectively.\footnote{If we assume that mmWwave communication takes place, hybrid beamforming can be introduced as the best solution that achieves a good trade-off between cost and complexity as usually adopted in the literature. However, the study of the impact of channel aging in the mmWwave  region is left for future research   due to limited space.} The precoding vector is normalized based on the average total power constraint
\begin{align}
\EE\{\|\bs_{n}\|^{2}\}=\tr(\bP\bF_{n}^{\H}\bF_{n})\le  P_{\mathrm{max}}\label{PowerConstraint},
\end{align}
where $ _n=[_1,n, …, _K,n] ∈ℂ^M ×K$, $ =(p_1, …, p_K)$, and $  P_max>0$ is the total transmit power. However, the channel $ _k,n $ in \eqref{DLreceivedSignal} can be expressed as
\begin{align}
\bh_{k,n}=&\al_{k,n-K}\hat{\bh}_{k}+\al_{k,n-K}\tilde{\bh}_{k}+\bar{\al}_{k,n-K}\tilde{\bee}_{k,n}\nn\\
\end{align}
where $ _k $ expresses the channel vector at the beginning of the data transmission phase, and $ _k $ is the corresponding channel estimation error.

Substitution of \eqref{innov2} into  \eqref{DLreceivedSignal} provides
\begin{align}
r_{k,n}&=\al_{k,n-K}\sqrt{ p_{k}}{\bh}_{k}^{\H}\bff_{k,n}x_{k,n}
+\bar{\al}_{k,n-K}\sqrt{ p_{k}}\tilde{\bee}_{k,n}^{\H}\bff_{k,n}x_{k,n}\nn\\
&		+\sum_{i\ne k}^{K}\sqrt{ p_{i}}\bh^\H_{k,n}\bff_{i,n}x_{i,n}+z_{k,n}.\label{DLreceivedSignal1}
\end{align}

Although  UEs do not have instantaneous CSI, we can assume that UE $ k $ has access to $ {_k^_k,nx_k,n} $. Then, by using the technique in~[35], where UE $ k $ is aware of only the statistical CSI, the received signal is written as
\begin{align}
&r_{k,n}\!=\!\al_{k,n-K}\sqrt{p_{k}}\EE\{{\bh}_{k}^{\H}\bff_{k,n}\}x_{k,n}	\!+\!\!\bar{\al}_{k,n-K}\sqrt{p_{k}}\tilde{\bee}_{k,n}^{\H}\bff_{k,n}x_{k,n}\nn\\
&+\sum_{i\ne k}^{K}\sqrt{p_{i}}\bh^\H_{k,n}\bff_{i,n}x_{i,n}+z_{k,n}\label{DLreceivedSignal1}.
\end{align}

\begin{proposition}\label{LowerBound1}
The downlink average SE for UE $ k$ of an RIS-assisted mMIMO system, accounting for imperfect CSI and channel aging,  is lower bounded by
\begin{align}
\mathrm{SE}_{k}	=\frac{1}{\tau_{\mathrm{c}}}\sum_{n=K+1}^{\tau_{\mathrm{c}}}\log_{2}\left ( 1+\gamma_{k,n}\right)\!,\label{LowerBound}
\end{align}
where $ \gamma_{k,n}$ is the achievable SINR at time $ n $  given  by \eqref{sinr}.
\begin{figure*}
\begin{align}
\gamma_{k,n}=\frac{\al_{k,n-K}^{2}p_{k}|\EE\{{\bh}_{k}^{\H}\bff_{k,n}\}|^{2}}{\al_{k,n-K}^{2}p_{k}\mathrm{Var}\left\{{\bh}_{k}^{\H}\bff_{k,n}\right\}		+
\sum_{i\ne k}^{K}\al_{i,n-K}p_{i}\EE\left\{|\bh^\H_{k,n}\bff_{i,n}|^{2}\right\} 		+\bar{\al}_{k,n-K}^{2}p_{k}\EE\left\{|	\tilde{\bee}_{k,n}^{\H}\bff_{k,n}|^{2}\right\}+
\sigma^{2}}.\label{sinr}
\end{align}
\line(1,0){490}
\end{figure*}
\end{proposition}
\begin{proof}
First, based on a similar approach to~[36, 37],     the average achievable SE including the achievable SINR $ \gamma_{k,n} $ in the transmission phase is computed  for each $ n $. Next,   the average over these SEs is obtained as in \eqref{LowerBound}.

Regarding the achievable SINR $ \gamma_{k,n} $,  given  the Gaussianity of the input symbols, it  is obtained by accounting for    a worst-case assumption for the computation of the mutual information~\cite[Theorem $1$]{Hassibi2003}. In particular, except for $ \EE\{\bh^\H_{k,n}\bff_{k,n}\} $ being the deterministic desired signal detected by UE $ k $, all others terms are treated as independent Gaussian noise with zero mean and  variance equal to the variance of interference plus noise.
\end{proof}

The following analysis requires $ M $, $ L $, and $ K $ increase  but with a given bounded ratio as $ 0 < liminfK
/M< ∞$  and $ 0 < liminfL
/M< ∞$. Henceforth, this notation is denoted as $ M →∞$. Taking into account that, according to \eqref{innov2}, the available CSI at time $ n $ is $ _k,n-K_k $, the BS designs its RZF precoder as\footnote{Despite that RZF precoding is generally suboptimal,  it has been applied selected in the case of mMIMO in many works due to reasons of complexity and to provide closed-form expressions. In other words, it is common in mMIMO to consider a linear precoder between the MRT and RZF precoders. Hence, for these reasons, in this work, we have selected the more optimal precoder, i.e., RZF despite its complexity. Note that  RZF precoding is a very good choice compared to MRT, and the corresponding derivation of the rate together with its optimization require delicate manipulations, which raise the novelty of this work.}
\begin{align}
\bff_{k,n}=
\al_{k,n-K}\sqrt{\lambda} \bSigma\hat{\bh}_{k}, \label{precoderRZF}
\end{align}
where $ λ$ is   a normalization parameter, which is obtained due to \eqref{PowerConstraint} as
\begin{align}
\lambda=\frac{P_{\mathrm{max}}}{	\al_{k,n-K}^{2}\tr(\bP\hat{\bH}^{\H}\bSigma^{2}\hat{\bH})}. \label{eq:lamda}
\end{align}
Note that $≜(_k,n-K^2^ + + M _M)^-1$, where $ ∈ℂ^M ×M $ is  an arbitrary Hermitian non negative definite matrix, and $$ is a regularization scaled by $M$ to make expressions converge to a constant as $M →∞$.  As a result, by denoting $ ρ=P_max/	σ^2 $ the SNR at the downlink transmission phase, the SINR of UE $ k $ under RZF precoding can be written as \eqref{SINRwithPrecoding} at the top of the next page. For the sake of convenience, we denote $ S_k,n $ and $ I_k,n $ the numerator and denominator of \eqref{SINRwithPrecoding}, respectively, i.e., we have  $ γ_k,n=S_k,n/I_k,n$.
\begin{figure*}
\begin{align}
\gamma_{k,n}=\frac{p_{k}|\EE\{{\bh}_{k}^{\H}\bSigma\hat{\bh}_{k}	 \}|^{2}}{p_{k}\mathrm{Var}\left\{{\bh}_{k}^{\H}\bSigma\hat{\bh}_{k}	 \right\}		+
\sum_{i\ne k}^{K}p_{i}\frac{\al_{i,n-K}^{4}}{\al_{k,n-K}^{4}}\EE\left\{|\bh^\H_{k,n}\bSigma\hat{\bh}_{i}	 |^{2}\right\} 		+p_{k}\frac{\bar{\al}_{k,n-K}^{2}}{\al_{k,n-K}^{2}}\EE\left\{|	\tilde{\bee}_{k,n}^{\H}\bSigma\hat{\bh}_{k}	 |^{2}\right\}+
\frac{\tr(\bP\hat{\bH}^{\H}\bSigma^{2}\hat{\bH})}{\al_{k,n-K}^{2}\rho}}\label{SINRwithPrecoding}.
\end{align}
\line(1,0){490}
\end{figure*}

Based on similar assumptions to \cite[Assump. A1-A3]{Hoydis2013} regarding the covariance matrices under study, we rely on the DE analysis to obtain the DE downlink SE of UE $ k $.\footnote{The DE analysis results in deterministic expressions, which make lengthy Monte-Carlo simulations unnecessary. Also, its results are tight approximations even for conventional systems with moderate  dimensions, e.g., an $ 8 \times 8 $ matrix~[39]. Hence, the following DE SINR in \eqref{DLdelayedCSIetaRZF} is of great practical importance.} 	The DE SINR  obeys to $γ_k,n-γ̅_k,n0$ while the deterministic SE of UE $k$ obeys to
\begin{align}
\mathrm{SE}_{k}-\mathrm{\overline{SE}}_{k}\xrightarrow[M \rightarrow \infty]{\mbox{a.s.}}0,\label{DeterministicSumrate}
\end{align}
where $SE_k= 1/T_c∑_n=K+1^T_clog_2(1 + γ̅_k,n) $ based on the dominated convergence and the continuous mapping theorem~[39].

\begin{Theorem}\label{theorem:DLagedCSIRZF}
The downlink DE of the SINR of UE $k$ with RZF precoding at time $n$, accounting for correlated Rayleigh fading, imperfect CSI, and channel aging due to  UE mobility  is given by \eqref{DLdelayedCSIetaRZF}, where
\begin{figure*}
\begin{align}
\!	\!\bar{\gamma}_{k,n} \!=\! \frac{ p_{k}\delta_{k}^{2}	 }{\displaystyle p_{k}\tilde{\delta}_{k}\!+\!	p_{k}\frac{\bar{\al}_{k,n-K}^{2}}{\al_{k,n-K}^{2}}\delta_{k}^{\mathrm{e}}\!+\!\sum_{i\ne k}\!p_{i}\frac{\al_{i,n-K}^{4}(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}}{\al_{k,n-K}^{4}(1\!+\!\al_{i,n-K}^{2}{\delta_{i}})^{2}}Q_{ik}\!+\!	\frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}	\delta_{i}^{\lambda}	 }{\rho(1+	\al_{i,n-K}^{2}\delta_{i})^{2}}}\label{DLdelayedCSIetaRZF}
\end{align}
\line(1,0){490}
\end{figure*}
\begin{align}
\bar{\lambda}&=\frac{P_{\mathrm{max}}}{  \frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{\al_{k,n-K}^{2}	\delta_{i}^{\lambda}	 }{(1+	\al_{i,n-K}^{2}\delta_{i})^{2}}},
\end{align}
$\delta_{k}= \frac{1}{M}\tr(\bPhi_{k} \bT ),
{\tilde{\delta}}_{k}=\frac{1}{M^{2}}\tr \big((\bR_{k}-\bPhi_{k})\tilde{\bT}(\bPhi_{k})\big),
\delta_{k}^{\mathrm{e}}=\frac{1}{M}\tr\big(\bR_{k}\tilde{\bT}(\bPhi_{k})\big)$, 	$\delta_{i}^{\lambda}=\frac{1}{M}\tr\big(\bPhi_{i}\tilde{\bT}(\Id_{M})\big)		$,
$ \zeta_{ki}=\frac{1}{M^{2}}\tr\big(\bR_{k}\tilde{\bT}(\bPhi_{i})\big) $, 	 $ 	\mu_{ki}=\frac{1}{M^{2}}\tr\big(\bPhi_{k}\tilde{\bT}(\bPhi_{i})\big) $,
$Q_{ik}= \zeta_{ki}\!+\!\frac{|\delta_{k}|^{2}\mu_{ki}}{\left(1\!+\!\al_{k,n-K}^{2}\delta_{k}\right)^{2}}\!-\!2\mathrm{Re}\left\{\!\!
\frac{ \delta_{k}^{*}\mu_{ki}}{1\!+\!\al_{k,n-K}^{2}\delta_{k}}\right\}$,
\begin{itemize}
\renewcommand{\labelitemi}{$\ast$}
\item $\bT\!=\!\Big(\!\displaystyle\sum_{i=1 }^{K} \frac{1}{\!M\left(1+\al_{i,n-K}^{2}\delta_{i}\right)}\bPhi_{i}\!+\!\al \Id_{M}\Big)^{-1}$,
\item $ \tilde{\bT}(\bL)=\bT\bL\bPhi_{i}\bT+\sum_{i=1}^{K}\frac{\tilde{\delta}_{i}\bT\bPhi_{i}\bT}{M\left(1+\al_{i,n-K}^{2}\delta_{i}\right)^{2}}$, $ \bL=\bPhi_{i}, \Id_{M} $
\item $ \tilde{\deltav}(\bL) \!\!=\left(\Id_{K}\!-\!\bF\right)^{\!-1}\bf$ with $ 	\left[\bF\right]_{k,i}\!=\!\frac{1}{M^{2}\left(1+ \al_{i,n-K}^{2}\delta_{i}\right)}\!\tr\!\left(\bPhi_{k}\bT\bL\bT\right) $, $ \bL=\bPhi_{i}, \Id_{M} $.
%			 and $ \left[\bff\right]_{k}=\frac{1}{W}\tr\left(\bD_{k}\bPhi_{k}\bD_{k}\bT\bD_{k}\bPhi_{k}\bD_{k}\bT\right) $.
\end{itemize}
\end{Theorem}
\proof: Please see Appendix~\ref{theorem3}.\endproof

\section{Sum SE Maximization}\label{SumSEMaximizationDesign}
In this section, we focus on the  optimization of the sum SE of a RIS-assisted time-varying mMIMO  system with channel aging. Specifically, the sum rate maximization problem is described as
\begin{subequations}
\begin{align}
(\mathcal{P}1)~~&\max_{\bTheta,\bp\ge 0} 	\;	\mathrm{\overline{SE}}=\sum_{i=1}^{K}\mathrm{\overline{SE}}_{k}
%		~~~~\mathrm{for}~i=\mathrm{JGP}, \mathrm{PGP}
\label{Maximization1} \\
&~	\mathrm{s.t}~~~\;\!	|\tilde{\theta}_{l}|=1,~~ l=1,\dots,L,\label{Maximization2} \\
&\;\quad\;\;\;\;\;\!\!~\!\sum_{i=1}^{K}p_{i}\le P_{\mathrm{max}},
\label{Maximization3}
\end{align}
\end{subequations}
where  we have denoted the elements of $ $ as $θ̃_l= exp(j θ_l) $ for all $ l $ and the vector $ =[p_1, …,p_K]^ $. The first constraint in \eqref{Maximization2} means that each RIS element induces a phase shift without any change on the amplitude of the incoming signal, while  constraint \eqref{Maximization3}  ensures that the BS transmit power is kept below the maximum power $ P_max $.\footnote{The considered optimization problem allocates the available resources (power vector, RIS configuration) to maximize the total spectral efficiency (i.e, sum-rate). This is a well-known objective function when the total/aggregate rate of the network is the key design priority. It is worth also noting that the  sum capacity/rate is a fundamental performance metric (one  dimension) in multi-user networks from information theoretic perspective. To consider QoS per individual user, the consideration of other objective functions that take into account fairness (e.g. max-min rate) and/or minimum individual rates is required, which are beyond the scope of this paper and can be considered for future work.}
%, obtained by replacing the constraint \eqref{PowerConstraint}  with its asymptotically equivalent,  ensures that the BS transmit power is kept below the maximum power $ P_{\mathrm{max}} $.

The solution of the optimization problem $ 	(𝒫1) $ is challenging  due to its non-convexity and the  unit-modulus constraint concerning $ θ̃_l $. To tackle this difficulty, we consider the alternating optimization (AO) technique, where $ $ 	and  $ $ are going to be solved separately and iteratively. In particular,  first,  we focus on finding  the optimum  $ $ given  fixed $ $. Next, we solve for $ $ with  fixed $ $. The iteration of this procedure, where the sum rate increases  at  each iteration step, continues until convergence  to  the optimum value  since the sum rate is upper-bounded  subject to the power constraint \eqref{Maximization3}.  Note that all computations take place at the BS.

\subsection{RIS Configuration}
The exploitation of the RIS potentials implies the optimization of the RBM towards maximum sum SE. The presence of the RBM appears inside the covariance matrices  in the DE achievable SINR in \eqref{DLdelayedCSIetaRZF}. Given that the logarithm function is monotonic, it is sufficient to maximize $ γ̅_k,n $ instead of $SE_k $. Hence, by assuming infinite resolution phase shifters, the RBM optimization problem   is formulated as
\begin{align}\begin{split}
(\mathcal{P}2)~~~~~~~\max_{\bTheta} ~~~	&	\;	\mathrm{\overline{SE}}\\
%		~~~~\mathrm{for}~i=\mathrm{JGP}, \mathrm{PGP}\\
\mathrm{s.t}~~~&|\tilde{\theta}_{l}|=1,~~ l=1,\dots,L.
\end{split}\label{Maximization}
\end{align}

Although the  problem  $ 	(𝒫2) $ is  non-convex in terms of $ $ and it is subject to  a unit-modulus constraint regarding $ θ̃_l $, application of  the projected gradient 	ascent algorithm  can achieve a local optimal solution  by projecting the solution onto the closest feasible point at every step 	  until converging to a 	stationary point. Specifically, let the vector $ _t =[θ̃_t,1, …, θ̃_t,L]^$ include the phases at step  $ t $. The next step of the algorithm is described by
\begin{align}
\tilde{\bc}_{t+1}&=\bc_{t}+\tilde{\mu} \bq_{t},~
\bc_{t+1}=\exp\left(j \arg \left(\tilde{\bc}_{t+1}\right)\right),\label{sol2}
\end{align}
where the parameter $ μ̃$ denotes the step size and  $ _t=  γ̅_k,n_t^* $  describes the  ascent direction at step $ t $, given below by Proposition \ref{Prop:optimPhase}.\footnote{Given that the feasible set is the unit circle, then any point  $ x $ should be   projected on this circle, i.e., it should be $ x/|x| $, which is equal to $ \exp(j \angle(x)) $.} Note that  the suitable step size  at each iteration is selected based on  the backtracking line search [40]. The solution is obtained by the  projection problem $ min_|θ_l |=1, l=1,…,L-^2 $ under the unit-modulus constraint. 	Algorithm \ref{Algoa1} presents the overview of this procedure.  For the sake of convenience, we denote  the partial derivative with respect to $ _t^* $  by $ (·)' $.

\begin{algorithm}
\caption{Projected Gradient Ascent Algorithm for the RIS Design}
1.				 \textbf{Initialisation}: $ \bc_{0} =\exp\left(j\pi/2\right)\one_{L}$, $ \bTheta_{0}=\diag\left(\bc_{0}\right) $, $\bar{\gamma}_{k,n}=f\left(\bTheta_{0}\right) $ given by \eqref{DLdelayedCSIetaRZF}; $ \epsilon>0 $\\
2. \textbf{Iteration} $ t $: \textbf{for} $ t=0,1,\dots, $ do\\
3. $\bq_{t}=   \pdv{\bar{\gamma}_{k,n}}{\bc_{t}^{*}}  $, where $\pdv{\bar{\gamma}_{k,n}}{\bc_{t}^{*}} $ is given by Proposition \ref{Prop:optimPhase};\\
4. \textbf{Find} $ \tilde{\mu} $ by backtrack line search $( f\left(\bTheta_{0}\right),\bq_{t},\bc_{t})$ [40];\\
5. $ \tilde{\bc}_{t+1}=\bc_{t}+\tilde{\mu} \bq_{t} $;\\
6. 	$ \bc_{t+1}=\exp\left(j \arg \left(\tilde{\bc}_{l+1}\right)\right) $; $ \bTheta_{l+1}= \diag\left(\bc_{l+1}\right) $;\\
7. $ \bar{\gamma}_{k,n}^{t+1}=f\left(\bTheta_{t+1}\right) $;\\
8. \textbf{Until} $ \| \bar{\gamma}_{k,n}^{t+1}- \bar{\gamma}_{k,n}^{t}\|^{2} <\epsilon$; \textbf{Obtain} $ \bTheta^{\star}=\bTheta_{t+1}$;\\
9. \textbf{end for}\label{Algoa1}
\end{algorithm}

\begin{proposition}\label{Prop:optimPhase}
The derivative of $ \bar{\gamma}_{k,n} $ with respect to $ \bc_{t}^{*} $ is given  by
\begin{align}
\!\bar{\gamma}_{k,n}' =\frac{S_{k}'I_{k}-S_{k}I_{k}'}{I_{k}^{2}},\label{gam1}
\end{align}
\begin{align}
S_{k}'&= 2p_{k}\delta_{k}\delta_{k}', \label{sg11}\\
I_{k}'&=\displaystyle\tilde{\delta}_{k}'\!+\!(\delta_{k}^{\mathrm{e}})'\!+\!\sum_{i\ne k}\!\frac{\al_{i,n-K}^{2}(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}}{\al_{k,n-K}^{2}(1\!+\!\al_{i,n-K}^{2}{\delta_{i}})^{2}}Q_{ik}'	\nn\\
\!&+\frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{2(1\!+\!\al_{k,n-K}^{2}\delta_{k})\delta_{k'}	\delta_{i}^{\lambda}+(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}	(\delta_{i}^{\lambda})'	 }{(1+	\al_{i,n-K}^{2}\delta_{i})^{2}}\nn\\
\!&+\frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{\al_{i,n-K}^{2}(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}\rho\delta_{i}'	\delta_{i}^{\lambda}	 }{(1+	\al_{i,n-K}^{2}\delta_{i})^{3}}\nn\\
\!&+\sum_{i\ne k}\!\frac{2 \al_{i,n-K}^{2}\al_{k,n-K}^{2}\delta_{k}'(1\!+\!\al_{k,n-K}^{2}\delta_{k})}{\al_{k,n-K}^{2}\rho(1\!+\!\al_{i,n-K}^{2}{\delta_{i}})^{3}}Q_{ik}\!\nn\\
&-\!\sum_{i\ne k}\!\frac{3\al_{i,n-K}^{4}\delta_{i}'(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}}{\al_{k,n-K}^{2}(1\!+\!\al_{i,n-K}^{2}{\delta_{i}})^{5}}Q_{ik}\label{intDeriv}
\end{align}
with the auxiliary variables  given by
\begin{align}
\!\!	\!\!\!\! &\tilde{\delta}_{k}'\!=\!\frac{1}{M}\tr\!\big(\!(\bR_{k}'-\bPhi_{k}')\tilde{\bT}(\bPhi_{k})\!+\!(\bR_{k}-\bPhi_{k})\tilde{\bT}'(\bPhi_{k})\!\big),\label{sg41}\\
( &\delta_{k}^{\mathrm{e}})'\!=\!\frac{1}{M}\tr\big(\bR_{k}'\tilde{\bT}(\bPhi_{k})+\bR_{k}\tilde{\bT}'(\bPhi_{k})\big),\label{sg51}\\
& \bT'=-\bT(\bT^{-1})'\bT,	\\
&	(\bT^{-1})'=\frac{1}{M}\displaystyle \sum_{i=1 }^{K}\frac{\bPhi_{i}'\left(1+\al_{i,n-K}^{2}\delta_{i}\right)-\al_{i,n-K}^{2}\bPhi_{i}\delta_{i}'}{(1+\al_{i,n-K}^{2}\delta_{i})^{2}},\label{sg31}\\
&-\frac{1}{M}\sum_{i=1}^{K}\frac{2\al_{i,n-K}^{2}\tilde{\delta}_{i}\bT\bPhi_{i}\bT \delta_{i}'}{(1+\al_{i,n-K}^{2}\delta_{i})^{}},\\
& 	Q_{ik}'= \zeta_{ki}'\!+\!\frac{\delta_{k}\delta_{k}'\mu_{ki}+|\delta_{k}|^{2}\mu_{ki}}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!2}}\!+\!\frac{2\al_{k,n-K}^{2}|\delta_{k}|^{2}\mu_{ki}\delta_{k}'}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!3}}\nn\\
\frac{ (\delta_{k}^{*})'\mu_{ki}+\delta_{k}^{*}\mu_{ki}'}{1\!+\!\al_{k,n-K}^{2}\delta_{k}}+\!\!\!
\frac{ \al_{k,n-K}^{2}\delta_{k}^{*}\mu_{ki}\delta_{k}'}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!2}}\right\},\\
& 	\zeta_{ki}'=\frac{1}{M^{2}}\tr\big(\bR_{k}'\tilde{\bT}(\bPhi_{i})+\bR_{k}\tilde{\bT}'(\bPhi_{i})\big), \\
& 	\bar{\lambda}'=\frac{K \big( \frac{1}{M} \tr \left(\frac{\Zm}{M} +\al \Id_M\right) \tilde{\bT}'(\Id_{M})-\frac{1}{M}\tr\Tm'\big)}{ \big(\frac{1}{M}\tr\Tm - \frac{1}{M} \tr \left(\frac{\Zm}{M} +\al \Id_M\right) \tilde{\bT}(\Id_{M})\big)^{2}}.
\end{align}

\end{proposition}
\proof The proof  is given in Appendix~\ref{optimPhase}.\endproof

The RBM beamforming design is based on the gradient ascent , which results in a significant advantage because the gradient ascent is derived in a closed-form. This method comes with low computational complexity because it consists of simple matrix operations. Specifically, the complexity of Algorithm \ref{Algoa1} is $ 𝒪(MN^2+N+M) $, which consists of the fundamental system parameters  $ M $ and $N $ with the number of RIS elements having the higher (square) impact.

\subsection{Power Allocation}
Given  a fixed RBM $ $, the objective is the maximization of the sum SE with respect to $ $. Specifically, we have
\begin{align}\begin{split}
(\mathcal{P}3)~~~~~~~\max_{\bp\ge 0} ~~~	&			\mathrm{\overline{SE}}\\
\mathrm{s.t}~~\;\!&\sum_{i=1}^{K}p_{i}\le P_{\mathrm{max}},
\end{split}\label{Maximization}
\end{align}
where  	$SE$ is given by \eqref{Maximization1}. 	This problem is not convex but a local optimal solution can be obtained by using a weighted minimum mean square error (MMSE) reformulation of the sum SE maximization. The SINR $γ̅_k,n $ in \eqref{DLdelayedCSIetaRZF} can be described as a function of the downlink power coefficients given by the vector $ $ as
\begin{align}
\bar{\gamma}_{k,n} =\frac{p_{k} q_{k}}{\bc^{\T}\bp},\label{SINR1}
\end{align}
where  $=[c_1, …,c_K ]^  $ with
\begin{align}
q_{k}&\!=\!\delta_{k}^{2},\! \\
c_{k}&\!=\!\tilde{\delta}_{k}\!+\!	\frac{\bar{\al}_{k,n-K}^{2}}{\al_{k,n-K}^{2}}\delta_{k}^{\mathrm{e}}+\!\frac{1}{M}	\delta_{i}^{\lambda}\!,~\forall g \\
&+\frac{1}{M}\frac{(1\!+\!\al_{k,n-K}^{2}\delta_{k})^{2}	\delta_{i}^{\lambda}	 }{(1+	\al_{i,n-K}^{2}\delta_{i})^{2}},~\forall g, ~\forall i\ne g.
\end{align}

The MMSE reformulation includes writing  the SINR in \eqref{SINR1} in terms of a single‐input and single‐output (SISO) channel that is described as
\begin{align}
\tilde{y}_{k}=\sqrt{p_{k} q_{k}}s_{k}+\sum_{i=1}^{K}\sqrt{p_{i}c_{i}}s_{i},
\end{align}
where $ ỹ_k $ is the received signal, and $ s_i ∈ℂ $ expresses the normalized and independent random data signal with $ {|s_i|^2=1} $. The receiver can compute an estimate $ŝ_k=v_k^* 	ỹ_k $ of the desired signal $ s_k $ by minimizing the  MSE 	$ e_k(,v_k) =[|ŝ_k-s_k|^2]$, where  $ v_k $ as a scalar combining coefficient. Specifically, the MSE is written as
\begin{align}
e_{k}(\bp,v_{k})=v_{k}^{2}\left(p_{k}q_{k}+\bc^{\T}\bp\right)-2 v_{k}\sqrt{q_{k}P_k}+1.\label{mse1}
\end{align}

The coefficient $ 	v_k $, which minimizes the MSE 	$ e_k(,v_k)$ for a given $ $, is given by
\begin{align}
v_{k}=\frac{\sqrt{p_{k} q_{k}}}{p_{k}q_{k}+\sum_{i=1}^{K}p_{i}c_{i}}.
\end{align}

Substituting $ v_k $ into \eqref{mse1}, the MSE becomes $ 1/(1+	γ̅_k,n ) $. According to the  weighted MMSE method, we introduce the auxiliary weight parameter $ d_k≥0 $ for the MSE $  e_k$ and focus on the solution of the following  optimization problem
\begin{align}\begin{split}
(\mathcal{P}4)~\min_{\substack{\bp\ge 0,\\
\{v_{k}, d_{k}\ge 0: k=1,\ldots,K\}}} 	&		\!\sum_{i=1}^{K}\!d_{i} e_{i}(\bp,\bv_{i})-\ln(d_{1})\\
\mathrm{s.t}~~\;\!&\sum_{i=1}^{K}p_{i}\le P_{\mathrm{max}}.
\end{split}\label{Maximization10}
\end{align}

It is worthwhile to mention that  problems  $ 	(𝒫3) $ and $ 	(𝒫4) $ are equivalent because they have the same global optimal solution. The equivalence  relies on the fact that the optimal $ d_k $ in \eqref{Maximization10} is $ 1/ e_k=(1+		γ̅_k,n)$. The benefit of  the reformulation results in the following lemma, which is  adapted based on \cite[Th. 3]{Shi2011}.

\begin{lemma}
The block descent coordinate algorithm, described as Algorithm \ref{Algoa2} below, converges to a local optimum of $ 	(\mathcal{P}4) $ by means of AO among three blocks of variables being $ \{v_{k}:  k=1,\ldots,K\} $, $ \{d_{k}:  k=1,\ldots,K\} $, and $ \bp $.
\end{lemma}

\begin{algorithm}
\caption{Block coordinate descent algorithm for solving $ 	(\mathcal{P}4) $}
%		the sum		SE maximization problem in \eqref{Maximization10}}
1.				 \textbf{Initialisation}: Set $ \bp=\frac{P_{\mathrm{max}}}{K}\one_{K} $ (arbitrary value) and the solution accuracy $ \epsilon >0 $, \\
2. \textbf{while} the objective function in \eqref{Maximization10} is not improved more than $ \epsilon $ \textbf{do}\\
3. $ 	v_{k}=\frac{\sqrt{p_{k} q_{k}}}{p_{k}q_{k}+\sum_{i=1}^{K}p_{i}c_{i}}, $ $ i=1,\ldots, K $\\
4. $ d_{k}=1/e_{k}(\bp,v_{k}), $ $ k=1,\ldots, K $\\
5. Solve the following problem for the current values of $ v_{k} $ and $ d_{k} $:\\\vskip-12mm
\begin{align}\begin{split}
(\mathcal{P}5)~~~~~~~\min_{\bp \ge 0} ~~~	&		\sum_{i=1}^{K}d_{i} e_{i}(\bp,\bv_{i})\\
\mathrm{s.t}~~\;\!&\sum_{i=1}^{K}p_{i}\le P_{\mathrm{max}},\\
\end{split}\label{Maximization11}
\end{align}	\vskip-2mm
6. Update $ \bp $ by the obtained solution to \eqref{Maximization11}\\
7. \textbf{end while}\\
8. \textbf{Output:} $ \bp^{\star} $\label{Algoa2}
\end{algorithm}

Algorithm \ref{Algoa2}  delineates the whole operation for the power allocation.  Step $ 5 $  includes a subproblem that has to be solved in every iteration. Fortunately, its solution can be obtained in closed-form with    $ √(p_k), k=1…, K $  treated as optimization variables by   decomposing  the problem into $ K$ independent subproblems. In particular, the solution is given by
\begin{align}
\end{align}

%	Actually, this subproblem is a convex quadratically-constrained quadratic programming problem with respect to $ \sqrt{P_{g}}, g=1\ldots, G $, which does not admit a closed-form solution but can be solved by using the solver SDPT3 with CVX.

\begin{remark}
The proposed algorithms, i.e., Algorithms \ref{Algoa1} and \ref{Algoa2} converge quickly and have low computation complexity. Moreover, given that both algorithms achieve  a local optimum,  and that the overall algorithm is based on AO,  the final solution corresponds to a local optimum, which means that different initializations  will result  in different solutions,  as will be shown below in Sec. \ref{Numerical}.
\end{remark}

The power allocation presents a similar complexity to the RBM design since similar matrix operations take place in Algorithm 2, i.e., its complexity is $ 𝒪(MN^2+N+M) $.

\section{Numerical Results}\label{Numerical}
In this section, we elaborate on the numerical results concerning the downlink sum SE of RIS-assisted mMIMO systems with  channel aging, imperfect CSI, and correlated Rayleigh fading. MC simulations, portrayed by "\ding{53}" marks, verify our analysis for $n →∞$ even for finite (conventional) system dimensions.  This property has been already observed in previous works relying on DE analysis~[39, 38, 20]. For the sake of comparison, in certain cases, we also depict the scenarios with MRT precoding, static UEs, and the absence of RIS.

The simulation setup consists of a a single cell with radius  $R = 1000$ meters. In the center of the cell, there is a BS, which consists of a  uniform linear array (ULA) with $ M =100$ antennas that serve   $ K = 20 $ moving UEs while the communication is assisted  by an RIS assembled by a uniform planar array (UPA) of $ L=100 $ elements.  The length of the uplink training duration is $τ=K$ symbols and the pilot transmit power is $p_p=6 dBm$. Assuming that the coherence time and bandwidth are $T_c=2 ms$ and $B_c=100 kHz$, respectively, i.e., the coherence block consists of $ 200 $ channel uses. Also,   $ σ^2=-174+10log_10B_c $. For the sake of exposition, we assume that the temporal correlation coefficient $ _k,n $ is the same across all UEs.

The correlation matrices $ _BS,k $ and $ _RIS,k $ are  computed according to [38] and [30], respectively. The size of each RIS element, required by the latter, is given by $ d_H=d_V=λ/4 $. In addition, the  path-losses corresponding to  the BS-to-RIS and RIS-to-UE $ k $ links are given by [3, 6]
\begin{align}
\beta_{1}=\frac{C_{1}}{r_{1}^{\al_{1}}},~~~\beta_{2,k}=\frac{C_{2}}{r_{2,k}^{\al_{2}}},
\end{align}
where $ _i $ and $ r_i $ for $ i=1,2 $ are the path-loss exponent and distance for link $ i $. Note that  $ C_1=26 $ dB and $ C_2=28 $ dB, which describe  the path-losses at a reference distance of $ 1$ m  [12]. In the case of the LoS link, we consider the same value for $ β_d,k $ as for $ β_2,k $ but with an additional penetration loss of $ 15 dB $. The choice of these values  has been made based on the 3GPP Urban Micro (UMi) scenario from TR36.814 for a carrier frequency of $ 2 $ GHz and a noise level of $ -80 $ dBm. In particular,  the path-losses for links $ 1 $ and $ 2 $ are generated based on the LOS and NLOS versions [42].  Specifically, we have $ _1 =2.2$ and $ _2 =3.67$ while $ r_1 =8 m$ and $ r_2 =60 m$. The Doppler spread, corresponding to a relative velocity of  $135$ km/h between the BS and the UEs, is $f_D=250 Hz$.    Moreover, if we assume that the   bandwidth is $W=20MHz$, the symbol time is $T_s=1/(2W)=0.025 μs$.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph1.eps}
\caption{\footnotesize{NMSE of UE $ k $ versus the SNR of an RIS-assisted MIMO system with imperfect CSI ($ M=100 $, $ L=100 $, $ K=20 $) for different UE velocities in the cases of uncorrelated fading at the RIS (Analytical results and MC simulations). }}
\label{Fig1}
\end{center}
\end{figure}

Fig. \ref{Fig1} depicts the relative estimation error per channel element, i.e., the NMSE versus the downlink SNR for different velocities. Certain noise floors appear  as $ p →∞$.
Moreover, we illustrate the result corresponding to no mobility, which decreases without bound. In addition, we  observe that the error floors take larger values with increasing  mobility (channel aging). It is shown that the NMSE saturates after $ 25 dB$ for high velocity while it saturates much later for slow velocity. Note that these results correspond to uncorrelated fading at the RIS  to avoid any RBM optimization because, in this case, the NMSE is independent of $ $
(Remark~\ref{rem1}). However, below, in the case of the sum SE, the optimization  with respect to $ $ is applied.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph2.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus the normalized Doppler shift $ f_{\mathrm{D}}T_{\mathrm{s}}  $ for different
number of RIS elements $ L $ ($ M=100 $, $ K=20 $)  (Analytical results and MC simulations). }}
\label{Fig2}
\end{center}
\end{figure}

In Fig \ref{Fig2}, we depict the achievable average  $ SE_n $ versus the normalized Doppler shift $ f_DT_s $ for $ n=160 $ while varying the number of RIS elements, i.e., $ L=20, 60, 120,180 $. The ripples  are caused because of the dependence from the correlation coefficient,  which takes the form of the Bessel function of zeroth order, i.e., $_k,n-K^2=J_0^2(2 πf_DT_sn)  $. It is shown that by increasing $f_DT_s$,  the average sum SE decreases. An increase in $ f_DT_s $ is equivalent to an increase in velocity or duration of the symbol time.  Especially, when $ f_DT_s≈0.39 $, the sum rate becomes almost zero for the first time, and then the magnitude fluctuates while tending to zero. Moreover, a larger RIS in terms of number of elements results in larger average sum SE but the curves keep the same shape. Even the zeros of the sum SE appear at the same normalized Doppler shifts. Hence, the presence of a RIS improves the performance even in channel aging conditions.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph3.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus the time instant index $ n $ for different
velocities $ v $ ($ M=100 $, $ L=100 $, $ K=20 $)  (Analytical results and MC simulations).  }}
\label{Fig3}
\end{center}
\end{figure}

Fig. \ref{Fig3} shows the achievable average  $ SE_n $ versus the time instant index $ n $ for different UEs velocities. The time index starts at $ n=K+1 $, where the data transmission begins. As can be seen, at a given time instant, $ SE_n $ decreases with UE mobility.  Moreover, as the velocity increases, the first zero position moves to the left, which suggests that the length of  coherence time should be designed to account for the impact of channel aging. Also, we depict the scenario with no RIS in the case $ v=40 km/h $, which appears lower SE.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph4.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus  the number of RIS elements $ L $ for different
velocities $ v $ ($ M=100 $,  $ K=20 $)  (Analytical results and MC simulations).  }}
\label{Fig4}
\end{center}
\end{figure}

Fig. \ref{Fig4} illustrates the achievable downlink sum  $ SE $ versus the number of RIS elements $ L $ for varying UE mobility in terms of the relative velocity $ v $. Obviously, $ SE $ increases with $ L $ as expected, but an increase in velocity degrades the performance. The larger the velocity becomes, the larger the loss in the sum SE is observed. However, an increase with respect to the number of elements still allows for improving the SE despite the degradation due to mobility.  Also, we show the results corresponding to static UEs and no RIS that correspond to the horizontal lines, which are parallel to the $ x $-axis. For the sake of further comparison, we illustrate the performance of the simpler MRT precoding under channel aging conditions (dotted lines), which increases with $ L $ but is inferior compared to RZF.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph5.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus  the number BS antennas $ M $ for different
velocities $ v $ ($ L=100 $, $ K=20 $)  (Analytical results and MC simulations). }}
\label{Fig5}
\end{center}
\end{figure}

Fig. \ref{Fig5} shows also the impact of channel aging in terms of the achievable sum SE with respect to the number of BS antennas $ M $ for different velocities. Generally, the downlink sum rate
presents  a rise with $ M $. Notably, we observe that when the velocity increases, the sum rate decreases while an increase in the number of BS antennas increases the performance. For the sake of comparison, we have included the solid lines in terms of simulation corresponding to the SE with CE performed according to [16]   when $  v=10 km/h$ and $  v=100 km/h $. Although the CE in [16] does not account for channel aging, the achievable sum SE is much lower than the proposed method due to its high training overhead.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph6.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus  the  number of iterations for different
velocities $ v $ and number of RIS elements $ L $ ($ M=100 $, $ K=20 $)  (Analytical results and MC simulations).}}
\label{Fig6}
\end{center}
\end{figure}

In Fig. \ref{Fig6}, we depict the convergence of the proposed entire algorithm, which is founded on AO and includes the subproblems of RBM and transmit power optimizations. In particular, the horizontal axis corresponds to the number of iterations. As can be observed, the convergence is achieved quite fast for the cases under study, i.e., about $ 6 $ iterations are required. In the same figure, we have shown how the increase in the number of RIS elements and the channel aging affect the convergence. As $ L $ increases, more iterations are required to reach convergence because the number of optimization variables is larger. However, channel aging does not affect the rate of convergence but only the rate since the number of optimization variables  does not change.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=0.9\linewidth]{Graph7.eps}
\caption{\footnotesize{Downlink achievable sum SE with RZF precoding of an RIS-assisted MIMO system with channel aging versus $ 30 $ channel realizations ($ M=100 $, $ L=100 $, $ K=20 $)  (Analytical results and MC simulations).}}
\label{Fig7}
\end{center}
\end{figure}

In addition, Fig. \ref{Fig7} assesses the dependence on the choice of the  initial points. Specifically,  given that the optimization problem \ref{Maximization3} is non-convex, its solution is contingent on the initialization. Hence, we consider  $30 $ channel realizations and $ v=50 km/h $. The overall algorithm, consisted of Algorithms $ 1 $ and $ 2 $ is initialized with $ _0 =exp(jπ/2)_L$ and $ =P_max/K_K  $. The line corresponding to "Alg. 1-Test"  considers the best initial point out of $ 100 $ random initial points for each channel instance. As observed, different initial points lead to different solutions. However,  the sum SE in both cases is almost identical, which indicates that the selected initialization regarding the phase shifts and the transmit power  is a good choice.
\section{Conclusion} \label{Conclusion}
In this paper, we studied the impact of channel aging on RIS-assisted mMIMO systems by taking additionally into account the effects of spatial correlation and imperfect CSI. Having introduced channel aging and correlated Rayleigh fading not only during the data transmission phase but also in the uplink training phase, we obtained the effective channel estimates and derived the DE achievable sum SE with RZF in closed form. Next, we provided its maximization with respect to the RIS phase shifts and power budget constraints by applying an efficient AO algorithm based on statistical CSI that reduces  both the computational complexity and the feedback overhead. Hence, we illustrated the impact of channel aging and how its interconnection with other fundamental parameters affects performance. For instance, a suitable selection of the number of RIS elements and the length of the frame duration can mitigate channel aging. Notably, this work suggests interesting directions for future research with the study of wideband systems being of prominent significance.
\begin{appendices}
\section{Proof of Theorem~\ref{theorem:DLagedCSIRZF}}\label{theorem3}
Regarding the  term in the desired signal power given by $ S_{k,n} $, we have
\begin{align}
\EE\{{\bh}_{k}^{\H}\bff_{k,n}\}&=	\al_{k,n-K} \EE\{(\hat{\bh}_{k}^{\H}+\tilde{\bh}_{k}^{\H}) \bSigma_{k}\hat{\bh}_{k}\}\label{ds0}\\
&\asymp	\al_{k,n-K}\EE\left\{\frac{\hat{\bh}_{k}^{\H}	 \bSigma_{k}\hat{\bh}_{k}}{1+	\al_{k,n-K}^{2}\hat{\bh}_{k}^{\H} \bSigma_{k}\hat{\bh}_{k}}\right\}\label{d1}\\
&\asymp	\al_{k,n-K}\frac{\frac{1}{M}\tr(\bPhi_{k} \bT )	 }{1+	\al_{k,n-K}^{2}\frac{1}{M}\tr(\bPhi_{k} \bT )}\label{d3}\\
&=\frac{	\al_{k,n-K}\delta_{k}	 }{1+	\al_{k,n-K}^{2}\delta_{k}},\label{d2}
\end{align}
where $ {\bSigma}_k $ is defined as
\begin{align}
{\bSigma}_k &=\big(\al_{k,n-K}^{2}\hat{\bH}_{k}\hat{\bH}_{k}^{\H} -	\al_{k,n-K}^{2}\hat{\bh}_{k}\hat{\bh}_{k}^{\H}+ \bZ + \al M\Id_M\big)^{-1}.\nn
\end{align}
In \eqref{d1}, we have applied the matrix inversion lemma. In \eqref{d3}, we have used  \cite[Lem. 14.3]{Bai2010} known as rank-1 perturbation lemma, \cite[Lem. B.26]{Bai2010}, and \cite[Theorem 1]{Wagner2012}. Also, we have set $\delta_{k}= \frac{1}{M}\tr(\bPhi_{k} \bT )	 $.

For the derivation of the DE of the normalization parameter $\lambda$, we focus on the denominator, and we result in
\begin{align}
&	\tr(\bP\hat{\bH}^{\H}\bSigma^{2}\hat{\bH})=\sum_{i=1}^{K}p_{i}\hat{\bh}_{i}^{\H}\bSigma^{2}\hat{\bh}_{i}  \\
&= \frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{\hat{\bh}_{i}^{\H}	 \bSigma_{i}^{2}\hat{\bh}_{i}}{(1+	\al_{i,n-K}^{2}\hat{\bh}_{i}^{\H} \bSigma_{i}\hat{\bh}_{i})^{2}}\label{lamda1}\\
&\asymp \frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{\frac{1}{M}\tr\big(\bPhi_{i} \tilde{\bT}(\Id_M)\big )	 }{(1+	\al_{i,n-K}^{2}\frac{1}{M}\tr(\bPhi_{i} \bT ))^{2}}\label{lamda2}\\
&=\frac{1}{M}\sum_{i=1}^{K}p_{i}\frac{	\delta_{i}^{\lambda}	 }{(1+	\al_{i,n-K}^{2}\delta_{i})^{2}},\label{desired1}
\end{align}
where, in \eqref{lamda1}, we have applied the matrix inversion lemma twice. In \eqref{lamda2}, we have applied the rank-1 perturbation lemma, \cite[Lem. B.26]{Bai2010}, and \cite[Theorem 1]{Wagner2012} together with  \cite[Theorem 2]{Hoydis2013} for $\bK=\Id_M$. In \eqref{desired1}, we have denoted $ \delta_{i}^{\lambda}=\frac{1}{M}\tr\big(\bPhi_{i} \tilde{\bT}(\Id_M) \big)$.

%	 The second and third terms in $ I_{k,n} $ vanish as $ M\to \infty $ as can be seen after applying the matrix inversion lemma, and then \cite[Lem. B.26]{Bai2010} since the corresponding vectors are independent.

The term in $ I_{k,n} $, which includes the deviation from the average effective
channel gain  is written as
\begin{align}
\mathrm{Var}\left\{{\bh}_{k}^{\H}\bff_{k,n}\right\}	&\asymp\frac{\al_{k,n-K}^{2}\EE\{|\tilde{\bh}_{k}^{\H}\bSigma_{k}\hat{\bh}_{k}|^{2}\}}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{var1}\\
&	\asymp \frac{\al_{k,n-K}^{2}\tr(\bSigma_{k}\bPhi_{k}\bSigma_{k}(\bR_{k}-\bPhi_{k}))}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{var2}\\
&	\asymp \frac{\frac{1}{M^{2}}\al_{k,n-K}^{2}\tr\big((\bR_{k}-\bPhi_{k})\tilde{\bT}(\bPhi_{k})\big)}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{var3},
\end{align}
where in \eqref{var1}, we have applied the matrix inversion lemma, \cite[Lem. B.26]{Bai2010}, and \cite[Theorem 1]{Wagner2012}. Next step includes application of the rank-1 perturbation lemma, \cite[Lem. B.26]{Bai2010}, and \cite[Lem. 10]{Krishnan2015}. and \cite[Lem. B.26]{Bai2010} again. In   \eqref{var3}, we have used \cite[Theorem 2]{Hoydis2013}.

Similar steps are followed to derive the DE of the term including the innovation error. Specifically, we have
\begin{align}
\EE\left\{|	\tilde{\bee}_{k,n}^{\H}\bff_{k,n}|^{2}\right\}&\asymp\frac{\al_{k,n-K}^{2}\EE\{|\tilde{\bee}_{k,n}^{\H}\bSigma_{k}\hat{\bh}_{k}|^{2}\}}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{error1}\\
&	\asymp \frac{\al_{k,n-K}^{2}\tr(\bSigma_{k}\bPhi_{k}\bSigma_{k}\bR_{k})}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{error2}\\
&	\asymp \frac{\frac{1}{M^{2}}\al_{k,n-K}^{2}\tr\big(\bR_{k}\tilde{\bT}(\bPhi_{k})\big)}{(1+\al_{k,n-K}^{2}\delta_{k})^{2}}\label{error3}.
\end{align}

The second term, concerning the multi-user interference, is obtained as
\begin{align}
\EE\{|\bh^\H_{k,n}\bff_{i,n}|^{2}\}&=\al_{i,n-K}^{2}\EE\bigg\{\!\left	|\frac{\bh^\H_{k,n}\bSigma_{i}\hat{\bh}_{i}}{1+\al_{i,n-K}^{2}\hat{\bh}_{i}^{\H}\bSigma_{i}\hat{\bh}_{i}}\right|^{2}\!\bigg\}\label{int1}\\
&=\al_{i,n-K}^{2}	\EE\bigg\{\!\frac{\bh^\H_{k,n}\bSigma_{i}\hat{\bh}_{i}\hat{\bh}_{i}^{\H}\bSigma_{i}\bh_{k,n}}{(1+\al_{i,n-K}^{2}\hat{\bh}_{i}^{\H}\bSigma_{i}\hat{\bh}_{i})^{2}}\!\bigg\}\label{int2}\\
&\asymp\al_{i,n-K}^{2}	\EE\bigg\{\!\frac{\bh^\H_{k,n}\bSigma_{i}\bPhi_{i}\bSigma_{i}\bh_{k,n}}{(1+\al_{i,n-K}^{2}\delta_{i})^{2}}\!\bigg\},\label{int3}
\end{align}
where the matrix inversion lemma has been applied in \eqref{int1}. In \eqref{int3}, the mutual independence between $ \hat{\bh}_{i} $ and $\bh_{k,n}  $, the rank-1 perturbation lemma, \cite[Lem. B.26]{Bai2010}, and \cite[Theorem 1]{Wagner2012} are taken into account. However, $ \bSigma_{i} $ is not independent of $ \bh_{k,n} $. For this reason, application of \cite[Lemma~2]{Hoydis2013} gives
\begin{align}
\bSigma_i={\bSigma}_{ik}-\al_{k,n-K}^{2}\frac{{\bSigma}_{ik}\hat{\bh}_{k}\hat{\bh}_{k}^{\H}{\bSigma}_{ik}}{1+\al_{k,n-K}^{2}\hat{\bh}_{k}^{\H} {\bSigma}_{ik}\hat{\bh}_{k} },\label{eq:theorem2.I.51}
\end{align}
where the new matrix ${\bSigma}_{ik}$ is  defined as
\begin{align}
{\bSigma}_{ik}=&\big(\al_{k,n-K}^{2}\hat{\bH}_{k}\hat{\bH}_{k}^{\H} \!-\al_{k,n-K}^{2}\hat{\bh}_{k}\hat{\bh}_{k}^{\H}-\al_{i,n-K}^{2}\hat{\bh}_{i}\hat{\bh}_{i}^{\H}	\! \nn\\
&+\bZ + M \al \Id_M \big)^{-1}.
\end{align}

By inserting \eqref{eq:theorem2.I.51} into \eqref{int3}, we result in
\begin{align}
\EE\{|\bh^\H_{k,n}\bff_{i,n}|^{2}\}=\frac{Q_{ik}}{\big(1+\al_{i,n-K}^{2}{\delta_{i}}\big)^{2}},\label{eq:theorem2.I.6}
\end{align}
where $Q_{ik}$ is given by \eqref{eq:theorem2.I.mu1} at the top of the next page.
\begin{figure*}
\begin{align}
\!\!\!\!\!\!\!\!\!\!\!\!Q_{ik}&\!= \!\bh^\H_{k,n}\bSigma_{ik} \bPhi_{i} \bSigma_{ik}  \bh_{k,n}\!+\!\al_{k,n-K}^{6}\frac{\big|  \bh^\H_{k,n}{\bSigma}_{ik} \hat{\bh}_{k}\big|^{2}\!\hat{\bh}^\H_{k}  {\bSigma}_{ik} \bPhi_{i}{\bSigma}_{ik}\hat{\bh}_{k}}{\big( 1+\al_{k,n-K}^{2}\hatvh^\H_{k} {\bSigma}_{ik} \hatvh_{k} \big)^{2}}\!-\!2\al_{k,n-K}^{4}\mathrm{Re}\!\left\{ \! \frac{\hatvh^\H_{k}{\bSigma}_{ik}\bh_{k,n}\bh_{k,n}^{\H}{\bSigma}_{ik}\bPhi_{i}{\bSigma}_{ik}\hatvh_{k}}{1+\al_{k,n-K}^{2}\hatvh^\H_{k} {\bSigma}_{ik} \hatvh_{k}}\!\right\}\!.
\label{eq:theorem2.I.mu1}
\end{align}
\line(1,0){470}
\end{figure*}
In \eqref{eq:theorem2.I.6}, we have applied again the matrix inversion lemma, and in the last step, we have used the rank-1 perturbation lemma, \cite[Lem. B.26]{Bai2010}, \cite[Theorem 1]{Wagner2012}, and \cite[Theorem 2]{Hoydis2013}.  The DE of each term in~\eqref{eq:theorem2.I.mu1} is obtained as
\begin{align}
\bh^\H_{k,n}\bSigma_{ik} \bPhi_{i} \bSigma_{ik}  \bh_{k,n}&\asymp  \frac{1}{M^{2}}\tr \bR_{k}\tilde{\bT}(\bPhi_{i})=\zeta_{ki}\label{q1},\\
% \begin{align}
% \left|  \bg^\H_{k,n} {\bSigma}_{ik}\bA_{n} \bg_{k,0}\right|^{2}\asymp \frac{1}{M^{2}}\tr \bA_{n}^{2}\bD_{k}\bT''= \frac{ \delta_{k}^{'}}{M}
% \end{align}
\hat{\bh}^\H_{k}  {\bSigma}_{ik} \bPhi_{i}{\bSigma}_{ik}\hat{\bh}_{k}&\asymp  \frac{1}{M^{2}}\tr \bPhi_{k}\tilde{\bT}(\bPhi_{i})=\mu_{ki},\\
\bh_{k,n}^{\H}{\bSigma}_{ik}\bPhi_{i}{\bSigma}_{ik}\hatvh_{k}&\asymp \al_{k,n-K} \frac{1}{M^{2}}\tr\bPhi_{k}\tilde{\bT}(\bPhi_{i}).\label{q3}
\end{align}
Hence, substitution of \eqref{q1}-\eqref{q3} into \eqref{eq:theorem2.I.mu1}  provides $ Q_{ik} $.

The definitions of the various parameters are given in the presentation of the theorem.  The DE SINR $\bar{\gamma}_{k,n} $ is obtained by substituting \eqref{d2}, \eqref{desired1}, \eqref{var3}, \eqref{error3}, and \eqref{eq:theorem2.I.6} into $ S_{k,n} $ and $ I_{k,n} $.

\section{Proof of Proposition~\ref{Prop:optimPhase}}\label{optimPhase}
For this proof, we are going to use the following lemma.
\begin{lemma}\label{traceProd}
Let 	 $ \bA \in\mathbb{C}^{M\times M} $ be independent of $ \bTheta$ and $\bR_{k}= \bH_{1} \bTheta$ $\bR_{\mathrm{RIS},k}\bTheta^{\H}\bH_{1}^{\H} $, then
\begin{align}
\tr\left( \!\!\bA\pdv{\bR_{k}}{\bc_{t}^{*}}\!\right) =\mu\diag\left(\bH_{1}^{\H}\bA\bH_{1} \bTheta\bR_{\mathrm{IRS},k}\right).
\end{align}\end{lemma}
\proof We have
\begin{align}
\!\tr\!\left(\! \!\bA\pdv{\bR_{k}}{\bc_{t}^{*}}\!\right) &\!=\!\pdv{\left(\diag\left(\bH_{1}^{\H}\bA\bH_{1} \bTheta\bR_{\mathrm{IRS},k}\right)\right)^{\T}\bc_{t}^{*}}{\bc_{t}^{*}}\\
&= \mu\diag\left(\bH_{1}^{\H}\bA\bH_{1} \bTheta\bR_{\mathrm{IRS},k}\right),\label{deriv4}
\end{align}
where we have used the property $ \tr\left(\bA \diag(\bc_{t})\right)=\left(\diag(A)\right)^{\T}\bc_{t} $.
\endproof
The  gradient of $ \!\bar{\gamma}_{k,n}  $  with respect to $ \bc_{t}^{*} $ is written as
\begin{align}
\!\bar{\gamma}_{k,n}' =\frac{S_{k}'I_{k}-S_{k}I_{k}'}{I_{k}^{2}},\label{gam1}
\end{align}
where the quotient rule derivative was simply applied. We continue with the computation of the partial derivatives $ S_{k}' $ and $ I_{k}' $. Specifically, $ S_{g}' $ is obtained as
\begin{align}
S_{k}'= 2p_{k}\delta_{k}\delta_{k}', \label{sg1}
\end{align}
\begin{align}
\delta_{k}'=\frac{1}{M}\tr\left(\bPhi_{k}'\bT+\bPhi_{k}\bT'\right).\label{sg2}
\end{align}
Note that $ \bPhi_{k}' $ is the derivative of $ \bPhi_{k} $  with respect to $\bc_{t}^{*}  $ whose expression is given by
\begin{align}
\bPhi_{k}'=\al^{2}_{k,\zeta-n}\left(\bR_{k}' \bQ \bR_{k}+\bR_{k} \bQ' \bR_{k}+\bR_{k} \bQ \bR_{k}'\right).
\end{align}
The computation of the trace of this expression requires Lemma \ref{traceProd} and that the derivative of the  inverse matrix $ \bQ $ is obtained according to \cite[Eq. 40]{Petersen2012} as  $ \bQ'=-\bQ(\bQ^{-1})'\bQ$, where the derivative of $ \bQ^{-1} $ is written as
\begin{align}
\end{align}
Also, $ \bT' $ concerns the derivative of an inverse matrix written as $ \bT'=-\bT(\bT^{-1})'\bT$, where the derivative of $ \bT^{-1} $ is written as
\begin{align}
(\bT^{-1})'=\frac{1}{M}\displaystyle \sum_{i=1 }^{K}\frac{\bPhi_{i}'\left(1+\al_{i,n-K}^{2}\delta_{i}\right)-\al_{i,n-K}^{2}\bPhi_{i}\delta_{i}'}{(1+\al_{i,n-K}^{2}\delta_{i})^{2}},\label{sg3}
\end{align}
After inserting \eqref{sg2}-\eqref{sg3} into \eqref{sg1}, we obtain $S_{k}'  $.

The derivative of $ I_{k} $ is based on simple derivative rules and it requires the computation of the derivatives of $ \delta_{k} $, $ \tilde{\delta}_{k} $,	$ \delta_{k}^{\mathrm{e}} $, $ Q_{ik} $, and $ 	\bar{\lambda} $ as shown in \eqref{intDeriv}. Regarding $ \delta_{k}' $, it is given by \eqref{sg2}, while $ \tilde{\delta}_{k} $ and	$ \delta_{k}^{\mathrm{e}}$ are obtained similarly as
\begin{align}
\!\!\!\! \tilde{\delta}_{k}'&\!=\!\frac{1}{M}\tr\!\big(\!(\bR_{k}'-\bPhi_{k}')\tilde{\bT}(\bPhi_{k})\!+\!(\bR_{k}-\bPhi_{k})\tilde{\bT}'(\bPhi_{k})\!\big),\label{sg4}\\
( \delta_{k}^{\mathrm{e}})'&\!=\!\frac{1}{M}\tr\big(\bR_{k}'\tilde{\bT}(\bPhi_{k})+\bR_{k}\tilde{\bT}'(\bPhi_{k})\big),\label{sg5}
\end{align}
\begin{align}
\tilde{\bT}'(\bPhi_{k})&=\bT'\bPhi_{k}\bT+\bT\bPhi_{k}'\bT+\bT\bPhi_{k}\bT'\nn\\
&-\frac{1}{M}\sum_{i=1}^{K}\frac{2\al_{i,n-K}^{2}\tilde{\delta}_{i}\bT\bPhi_{i}\bT \delta_{i}'}{(1+\al_{i,n-K}^{2}\delta_{i})^{}}.
\end{align}

The derivative of the normalization parameter is given by
\begin{align}
\bar{\lambda}'&=\frac{K \big( \frac{1}{M} \tr \left(\frac{\Zm}{M} +\al \Id_M\right) \tilde{\bT}'(\Id_{M})-\frac{1}{M}\tr\Tm'\big)}{ \big(\frac{1}{M}\tr\Tm - \frac{1}{M} \tr \left(\frac{\Zm}{M} +\al \Id_M\right) \tilde{\bT}(\Id_{M})\big)^{2}}.\label{lambdaderiv}
\end{align}

Moreover, the last required derivative, concerning $ Q_{ik} $, is obtained as
\begin{align}
Q_{ik}'&= \zeta_{ki}'\!+\!\frac{\delta_{k}\delta_{k}'\mu_{ki}+|\delta_{k}|^{2}\mu_{ki}}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!2}}\!+\!\frac{2\al_{k,n-K}^{2}|\delta_{k}|^{2}\mu_{ki}\delta_{k}'}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!3}}\nn\\
\frac{ (\delta_{k}^{*})'\mu_{ki}+\delta_{k}^{*}\mu_{ki}'}{1\!+\!\al_{k,n-K}^{2}\delta_{k}}+\!\!\!
\frac{ \al_{k,n-K}^{2}\delta_{k}^{*}\mu_{ki}\delta_{k}'}{\big(1\!+\!\al_{k,n-K}^{2}\delta_{k}\big)^{\!2}}\right\},\label{qik}
\end{align}
\begin{align}
\zeta_{ki}'&=\frac{1}{M^{2}}\tr\big(\bR_{k}'\tilde{\bT}(\bPhi_{i})+\bR_{k}\tilde{\bT}'(\bPhi_{i})\big), \\
\mu_{ki}'&=\frac{1}{M^{2}}\tr\big(\bPhi_{k}'\tilde{\bT}(\bPhi_{i})+\bPhi_{k}\tilde{\bT}'(\bPhi_{i})\big).
\end{align}
Having obtained \eqref{sg4}, \eqref{sg5}, \eqref{lambdaderiv}, and \eqref{qik}, the derivative of $ I_{k} $ is derived and the proof is concluded.
\end{appendices}
\bibliographystyle{IEEEtran}

% Generated by IEEEtran.bst, version: 1.12 (2007/01/11)
\begin{thebibliography}{10}
\providecommand{\url}[1]{#1}
\csname url@samestyle\endcsname
\providecommand{\newblock}{\relax}
\providecommand{\bibinfo}[2]{#2}
\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
\providecommand{\BIBentryALTinterwordstretchfactor}{4}
\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
\fontdimen4\font\relax}
\providecommand{\BIBforeignlanguage}[2]{{%
\expandafter\ifx\csname l@#1\endcsname\relax
\typeout{** WARNING: IEEEtran.bst: No hyphenation pattern has been}%
\typeout{** loaded for the language `#1'. Using the pattern for}%
\typeout{** the default language instead.}%
\else
\language=\csname l@#1\endcsname
\fi
\providecommand{\BIBdecl}{\relax}
\BIBdecl

[1]
F.~Boccardi \emph{et~al.}, ``Five disruptive technology directions for {5G},''
\emph{IEEE Commun. Mag.}, vol.~52, no.~2, pp. 74--80, 2014.

[2]
E.~Basar \emph{et~al.}, ``Wireless communications through reconfigurable
intelligent surfaces,'' \emph{IEEE Access}, vol.~7, pp. 116\,753--116\,773,

[3]
Q.~Wu and R.~Zhang, ``Intelligent reflecting surface enhanced wireless network
via joint active and passive beamforming,'' \emph{IEEE Trans. Wireless
Commun.}, vol.~18, no.~11, pp. 5394--5409, 2019.

[4]
C.~Pan \emph{et~al.}, ``Multicell {MIMO} communications relying on intelligent
reflecting surfaces,'' \emph{IEEE Trans. Wireless Commun.}, vol.~19, no.~8,
pp. 5218--5233, 2020.

[5]
A.~Papazafeiropoulos \emph{et~al.}, ``Intelligent reflecting surface-assisted
{MU-MISO} systems with imperfect hardware: {Channel} estimation, beamforming
design,'' \emph{IEEE Trans. Wireless Commun.}, 2021.

[6]
A.~Kammoun \emph{et~al.}, ``Asymptotic max-min {SINR} analysis of
reconfigurable intelligent surface assisted {MISO} systems,'' \emph{IEEE
Trans. Wireless Commun.}, vol.~19, no.~12, pp. 7748--7764, 2020.

[7]
A.~Papazafeiropoulos \emph{et~al.}, ``Asymptotic analysis of {Max-Min} weighted
{SINR} for {IRS}-assisted {MISO} systems with hardware impairments,''
\emph{IEEE Wireless Commun. Lett.}, pp. 1--1, 2021.

[8]
A.~M. {Elbir} \emph{et~al.}, ``Deep channel learning for large intelligent
surfaces aided {mm-Wave} massive {MIMO} systems,'' \emph{IEEE Wireless
Commun. Lett.}, vol.~9, no.~9, pp. 1447--1451, 2020.

[9]
H.~Guo \emph{et~al.}, ``Weighted sum-rate maximization for reconfigurable
intelligent surface aided wireless networks,'' \emph{IEEE Trans. Wireless
Commun.}, vol.~19, no.~5, pp. 3064--3076, 2020.

[10]
J.~Chen \emph{et~al.}, ``Intelligent reflecting surface: {A} programmable
wireless environment for physical layer security,'' \emph{IEEE Access},
vol.~7, pp. 82\,599--82\,612, 2019.

[11]
P.~Yang, L.~Yang, and S.~Wang, ``Performance analysis for {RIS-Aided Wireless
Systems With Imperfect CSI},'' \emph{IEEE Wireless Commun. Lett.}, vol.~11,
no.~3, pp. 588--592, 2021.

[12]
E.~Bj{\"o}rnson, {\"O}.~{\"O}zdogan, and E.~G. Larsson, ``Intelligent
reflecting surface versus decode-and-forward: {How} large surfaces are needed
to beat relaying?'' \emph{IEEE Wireless Commun. Lett.}, vol.~9, no.~2, pp.
244--248, 2019.

[13]
M.~Di~Renzo \emph{et~al.}, ``Smart radio environments empowered by
reconfigurable intelligent surfaces: {How} it works, state of research, and
the road ahead,'' \emph{IEEE J. Sel. Areas Commun.}, vol.~38, no.~11, pp.
2450--2525, 2020.

[14]
D.~Mishra and H.~Johansson, ``Channel estimation and low-complexity beamforming
design for passive intelligent surface assisted {MISO} wireless energy
transfer,'' in \emph{International Conference on Acoustics, Speech and Signal
Processing (ICASSP)}.\hskip 1em plus 0.5em minus 0.4em\relax IEEE, 2019, pp.

[15]
Z.-Q. He and X.~Yuan, ``Cascaded channel estimation for large intelligent
metasurface assisted massive {MIMO},'' \emph{IEEE Wireless Commun. Lett.},
vol.~9, no.~2, pp. 210--214, 2019.

[16]
Q.~{Nadeem} \emph{et~al.}, ``Intelligent reflecting surface-assisted multi-user
{MISO Communication: Channel} estimation and beamforming design,'' \emph{IEEE
Open J. Commun. Soc.}, vol.~1, pp. 661--680, 2020.

[17]
B.~Zheng and R.~Zhang, ``Intelligent reflecting surface-enhanced {OFDM:
C}hannel estimation and reflection optimization,'' \emph{IEEE Wireless
Commun. Lett.}, vol.~9, no.~4, pp. 518--522, 2019.

[18]
E.~Shtaiwi \emph{et~al.}, ``Channel estimation approach for {RIS assisted MIMO}
systems,'' \emph{IEEE Trans. Cogn. Commun. Net.}, vol.~7, no.~2, pp.
452--465, 2021.

[19]
K.~Truong and R.~Heath, ``Effects of channel aging in massive {MIMO} systems,''
\emph{IEEE/KICS J. Commun. Net., Special Issue on massive {MIMO}}, vol.~15,
no.~4, pp. 338--351, Aug 2013.

[20]
A.~K. Papazafeiropoulos and T.~Ratnarajah, ``Deterministic equivalent
performance analysis of time-varying massive {MIMO} systems,'' \emph{IEEE
Trans. Wireless Commun.}, vol.~14, no.~10, pp. 5795--5809, 2015.

[21]
A.~K. Papazafeiropoulos, ``Impact of general channel aging conditions on the
downlink performance of massive {MIMO},'' \emph{IEEE Trans. Veh. Tech.},
vol.~66, no.~2, pp. 1428--1442, Feb 2017.

[22]
R.~Chopra \emph{et~al.}, ``Performance analysis of {FDD} massive {MIMO} systems
under channel aging,'' \emph{IEEE Trans. Wireless Commun.}, vol.~17, no.~2,
pp. 1094--1108, 2017.

[23]
Y.~Chen, Y.~Wang, and L.~Jiao, ``Robust transmission for reconfigurable
intelligent surface aided millimeter wave vehicular communications with
statistical {CSI},'' \emph{IEEE Trans. Wireless Commun.}, pp. 1--1, 2021.

[24]
M.-M. Zhao \emph{et~al.}, ``Intelligent reflecting surface enhanced wireless
network: {Two}-timescale beamforming optimization,'' \emph{IEEE Trans.
Wireless Commun.}, vol.~19, no.~8, pp. 5218--5233, 2020.

[25]
T.~Van~Chien \emph{et~al.}, ``Outage probability analysis of {IRS}-assisted
systems under spatially correlated channels,'' \emph{IEEE Wireless Commun.
Lett.}, pp. 1--1, 2021.

[26]
A.~Papazafeiropoulos \emph{et~al.}, ``Coverage probability of distributed {IRS}
systems under spatially correlated channels,'' \emph{IEEE Wireless Commun
Lett.}, pp. 1--1.

[27]
L.~You \emph{et~al.}, ``Reconfigurable intelligent surfaces-assisted multiuser
{MIMO Uplink Transmission With Partial CSI},'' \emph{IEEE Trans. Wireless
Commun.}, vol.~20, no.~9, pp. 5613--5627, 2021.

[28]
K.~Zhi \emph{et~al.}, ``Power scaling law analysis and phase shift optimization
of {RIS-Aided Massive MIMO Systems With Statistical CSI},'' \emph{IEEE Trans.
Commun.}, vol.~70, no.~5, pp. 3558--3574, 2022.

[29]
Y.~Zhang \emph{et~al.}, ``Performance analysis of reconfigurable intelligent
surface assisted systems under channel aging,'' \emph{Intelligent and
Converged Networks}, vol.~3, no.~1, pp. 74--85, 2022.

[30]
E.~{Bj{\"o}rnson} and L.~{Sanguinetti}, ``Rayleigh fading modeling and channel
hardening for reconfigurable intelligent surfaces,'' \emph{IEEE Wireless
Commun. Lett.}, vol.~10, no.~4, pp. 830--834, 2021.

[31]
D.~Neumann, M.~Joham, and W.~Utschick, ``Covariance matrix estimation in
massive {MIMO},'' \emph{IEEE Signal Process. Lett.}, vol.~25, no.~6, pp.
863--867, 2018.

[32]
S.~{Abeywickrama} \emph{et~al.}, ``Intelligent reflecting surface: {Practical}
phase shift model and beamforming optimization,'' \emph{IEEE Trans. Commun.},
vol.~68, no.~9, pp. 5849--5863.

[33]
E.~Bj\"{o}rnson, J.~Hoydis, and L.~Sanguinetti, ``Massive {MIMO} networks:
{Spectral}, energy, and hardware efficiency,'' \emph{Foundations and
Trends{\textregistered} in Signal Processing}, vol.~11, no. 3-4, pp.
154--655, 2017.

[34]
B.~Hassibi and B.~Hochwald, ``How much training is needed in multiple-antenna
wireless links?'' \emph{IEEE Trans. Inform. Theory}, vol.~49, no.~4, pp.
951--963, April 2003.

[35]
M.~Medard, ``The effect upon channel capacity in wireless communications of
perfect and imperfect knowledge of the channel,'' \emph{IEEE Trans. Inf.
Theory}, vol.~46, no.~3, pp. 933--946, May 2000.

[36]
E.~Bj{\"o}rnson, M.~Matthaiou, and M.~Debbah, ``Massive {MIMO} with non-ideal
arbitrary arrays: {H}ardware scaling laws and circuit-aware design,''
\emph{IEEE Trans. Wireless Commun.}, vol.~14, no.~8, pp. 4353--4368, Aug

[37]
A.~Pitarokoilis, S.~Mohammed, and E.~Larsson, ``Uplink performance of
time-reversal {MRC} in massive {MIMO} systems subject to phase noise,''
\emph{IEEE Trans. Wireless Commun.}, vol.~14, no.~2, pp. 711--723, Feb 2015.

[38]
J.~Hoydis, S.~ten Brink, and M.~Debbah, ``Massive {MIMO} in the {UL/DL} of
cellular networks: How many antennas do we need?'' \emph{IEEE J. Select.
Areas Commun.}, vol.~31, no.~2, pp. 160--171, February 2013.

[39]
R.~Couillet and M.~Debbah, \emph{Random matrix methods for wireless
communications}.\hskip 1em plus 0.5em minus 0.4em\relax Cambridge University
Press, 2011.

[40]
S.~Boyd, S.~P. Boyd, and L.~Vandenberghe, \emph{Convex optimization}.\hskip 1em
plus 0.5em minus 0.4em\relax Cambridge university press, 2004.

[41]
Q.~Shi \emph{et~al.}, ``An iteratively weighted {MMSE} approach to distributed
sum-utility maximization for a {MIMO} interfering broadcast channel,''
vol.~59, no.~9, pp. 4331--4340, 2011.

[42]
E.~Access, ``Further advancements for {E-UTRA} physical layer aspects,''
\emph{3GPP Technical Specification TR}, vol.~36, p.~V2, 2010.

[43]
Z.~Bai and J.~W. Silverstein, \emph{Spectral analysis of large dimensional
random matrices}.\hskip 1em plus 0.5em minus 0.4em\relax Springer, 2010,

[44]
S.~Wagner \emph{et~al.}, ``Large system analysis of linear precoding in
correlated {MISO} broadcast channels under limited feedback,'' \emph{IEEE
Trans. Inf. Theory}, vol.~58, no.~7, pp. 4509--4537, July 2012.

[45]
R.~Krishnan \emph{et~al.}, ``Linear massive {MIMO} precoders in the presence of
phase noise-{A} large-scale analysis,'' \emph{IEEE Trans. Veh. Tech.},
vol.~65, no.~5, pp. 3057--3071, 2016.

[46]
K.~B. Petersen and M.~S. Pedersen, ``The matrix cookbook, nov 2012,'' \emph{URL
http://www2. imm. dtu. dk/pubdb/p. php}, vol. 3274, p.~14, 2012.

\end{thebibliography}

\end{document}