A Simulation Framework for Cooperative Reconfigurable Intelligent Surface-Based Systems

We present a simulation framework for evaluating the performance of cooperative reconfigurable intelligent surface (RIS) based systems, which may ultimately deploy an arbitrary number of RISs to overcome adverse propagation-related effects, such as cascaded fading. The physical model underlying the proposed framework considers the (optional) presence of a dominant signal path between the source and RIS, and then between each subsequent stage of the communication link to the destination. Accompanying the dominant signal component is a non-isotropic scattered signal contribution, which accounts for angular selectivity within the cascaded RIS stages between the source and destination. The simulation of the time-correlated scattered signal, reflected by the illuminated reflective elements, is achieved using autoregressive modelling. As a by-product of our analysis, significant insights are drawn which enable us to characterize the amplitude and phase properties of the received signal, and the associated complex autocorrelation functions (ACFs) for the product of multiple Rician channels. For both single and cooperative RIS systems, the outage probability (OP), and important second-order statistics, such as the level crossing rate (LCR) and average outage duration (AOD), are analyzed for a variety of system configurations, accounting for practical limitations, such as phase errors. It is shown that by using multiple RISs cooperatively, the AOD is reduced at a lower signal-to-noise-ratio (SNR) compared to single RIS-assisted transmission under the same operating conditions. Lastly, increased channel variations (i.e., higher maximum Doppler frequencies) are shown to decrease the AOD in the case of absent phase errors; yet, this improvement is not observed when phase errors are present.


I. INTRODUCTION
T O TRULY enable the next generation of wireless tech- nologies, such as enhanced mobile broadband (eMBB), ultra-reliable low-latency communication (uRLLC), and massive machine-type communication (mMTC), calls for a major shift in the network infrastructure paradigm [2], [3].One potential approach is the utilization of intelligent meta-surfaces [4], which are two-dimensional planar arrays composed of a high number of passive reflective metaatoms/antennae (hereby referred to as elements) that are on the order of sub-wavelengths of the operating frequency [5], [6].The reflective elements within a meta-surface can be controlled to alter their reflective properties in order to manipulate the propagation of any signals which are impinged upon them [7].As a result, this property allows the control of phase shifts applied by each of the reflective elements [8].By using these smart surfaces distributed throughout a network, operators will be able to overcome unfavorable propagation phenomena in cellular environments [9], [10], in a cost effective manner [11].In the literature, the same technology has been referred to as software-controlled metasurfaces [9], intelligent reflective surfaces [12], [13], [14] or reconfigurable intelligent surfaces (RISs) [15], [16], [17].In the present contribution, we adopt the latter term, RISs, which are assumed to be composed of multiple reflective elements, that receive and re-radiate signals [4], [10], [12], [13], [14], [18], [19], [20].
A number of studies have investigated the theoretical performance of RISs in recent years [14], [19], [21], [22], [23], [24], [25], [26].To this end, the authors in [19] used the phase-shifts controlled by the illuminated reflective elements to fully optimize the received signal-to-noise-ratio (SNR).This work also showed that the outage probability (OP), the average symbol error rate, and the achievable rate, all are directly affected by the number of illuminated reflective elements, with greater numbers improving the end-to-end communications performance.It was shown in [21] that if two RISs work cooperatively to overcome a severe blockage, they can provide enhanced SNR compared to a single RIS enabled system.The authors in [27] showed that the OP and ergodic capacity improve with successive RISs working cooperatively to assist transmission, within Nakagami-m environments.While beyond the scope of this paper, we note that obtaining channel state information (CSI) is important for the successful operation of RISs and will be made more difficult for transmission via cascaded RISs.One possible solution is to incorporate passive beam steering between RISs working cooperatively, though it would increase the training overhead as more reflection coefficients need to be estimated [28].It was shown in [29] that the cascaded CSI is sufficient for designing the cooperative beam steering at two RISs to maximize the data transmission rate.Furthermore, in [28], the CSI estimation for two RISs working cooperatively is significantly improved (in terms of overhead), while having high accuracy.This is essential to maximize the time resources available for data transmission given a limited channel coherence interval.The authors in [24] and [25] showed that practical implementations of RISs will surely be subject to errors in the desired phase shifts due to hardware imperfections in the reflective elements, or imperfect estimation.Hence, as the desired phase shifts from the illuminated reflective elements are often used to optimize RIS operation, accounting for their errors is essential for realistic assessment of the achievable performance [13], [14], [15], [19], [26].While providing an opportunity to improve the endto-end performance of wireless links, increasing the number of RIS stages will also be subject to practical considerations (such as, e.g., reliable CSI acquisition).Not only will it present difficulties with the estimation of CSI as discussed above, it will also entail significant challenges, for example when localizing users, which has already been shown to be a challenging issue in single RIS systems [30].
It is recalled that the theoretical performance analysis of wireless communication systems often makes use of first-order statistics, such as the probability density functions (PDFs) of the received signal envelope and received signal phase.Although these are useful for interpreting the overall distribution of metrics, such as the SNR, OP, etc., crucially, they provide no insight on how these metrics are distributed with respect to time.On the contrary, the second-order statistics, such as the level crossing rate (LCR), do convey this information and can be used in conjunction with the OP to determine the corresponding average outage duration (AOD) [31], [32].For example, a primary use case of these concepts for RISs is in delay sensitive applications because of their potential to increase reliability and improve latency [22], [23].Thus, acknowledging the investigation of the AOD in RIS systems is critical and is a key motivator for this work.
Unlike the theoretical treatment of performance analysis [14], [19], [21], [22], [23], [24], [25], [26], the simulation of RIS systems has received less attention despite its importance for providing insights into the system performance in the presence of spatio-temporal variations.A simulator was implemented in [33] to model the RIS operation based on the 3GPP 5G channel study, discussed in [34].This simulator provided useful insights into how RIS placement is crucial for improving the achievable rates; however, it neither considered spatial correlation effects arising from the close proximity of the reflective elements in the RIS, nor did it incorporate the impact of Doppler related effects.In this contribution, we develop a simulation framework, which is built upon underlying Gaussian random variables (RVs), for wireless transmission from a source to destination through multiple cooperative RISs.Through our analysis, we have been able to obtain expressions for a number of the first-and second-order statistics associated with product (or equivalently cascaded) complex Rician channels, including the complex autocorrelation function (ACF).Importantly, through the complex ACF, we are able to account for important channel properties through each of the cascaded links, such as potential Doppler effects, a dominant signal component, and a non-isotropic scattered signal contribution which caters for angular selectivity.By further exploiting the complex ACF and the properties of the fading model, we harness the power of autoregressive (AR) modeling [35] to implement our simulator.We recall that AR modeling has previously been successfully used for channel prediction [36], and the simulation of correlated Rayleigh and Nakagami-m fading channels [35], [37].
Following from the above, the main contributions of this paper are summarized as follows: 1) We present a new comprehensive simulation framework that enables the generation of realizations of the complex received signal, encapsulating the autocorrelation properties of the signal reflected by each illuminated reflective element, and also the spatial correlations between them.2) Novel expressions for the distribution of the received signal phase for both single and cooperative RIS systems are derived, based upon the physical models representing the cascaded fading in each stage of the end-to-end communication link.3) We obtain a closed-form expression for the complex ACF of the product of two Rician channels, which encapsulates the temporal variations of a dominant signal component, a non-isotropic scattered signal component, and Doppler effects.The non-isotropic scattered signal component caters for angular selectivity of the direction of departure (DoD) and direction of arrival (DoA) of the signals at the transmitting source, a RIS, and the receiving destination, respectively.This is subsequently extended to the product of multiple Rician channels.4) We extend the system model of single RIS-assisted transmission including phase errors, to multiple cooperative RISs assisting transmission, which renders the proposed simulator particularly versatile.5) Capitalizing on the above contributions, we indicatively simulate RIS-assisted device-to-device (D2D) communications (making use of either a single or multiple cooperative RISs) and evaluate the corresponding OP, LCR, and AOD while considering various realistic RIS configurations, phase errors, and critically, different practical signal propagation conditions.The remainder of the paper is organized as follows: In Section II, the physical models are presented for both the single and cooperative RIS systems, along with some of the associated first-and second-order statistics for the product of multiple complex Rician channels, which arise as a result of the analysis performed here.Next, the simulation technique for generating the signal at each stage of an RIS-assisted transmission, is detailed in Section III.Following this, Section IV presents the system models for both single and cooperative RIS systems with phase errors.In Section V, single and cooperative RIS systems are simulated using our proposed framework and insights into the performance of the considered systems under the prescribed conditions are provided.Lastly, Section VI outlines our concluding remarks.
Notation: The complex conjugate of x is given by x * , x T denotes the transpose of vector x, arg(•) returns the argument, while ⊗ represents the Kronecker product.The {•} ∈ Q notation is used to indicate that the elements in {•} are members of set Q. Also, the real numbers greater than zero are denoted as R >0 , non-negative real numbers including zero are represented by R ≥0 , and the set of natural numbers greater than zero is denoted as N 1 .The statistical expectation operator is represented as E [•] and the variance operator is denoted by V [•], while P[•] denotes the probability operator.Circularly-symmetric complex Gaussian distributed RVs are denoted as X ∼ CN (a, b) where a is the mean value and b > 0 is the variance.The gamma function is denoted as Γ(•) [38,Eq. (8.310

II. PHYSICAL MODEL
For the single RIS system illustrated in Fig. 1(a), we assume that there is no direct link between the source and destination due to severe blockage effects.We also consider that there are L illuminated reflective elements available, and that the received signal at the destination is the superposition of signals coming from the L illuminated reflective elements; each of these signals is subject to a cascaded channel, corresponding to the paths between the source and the l th illuminated reflective element (i.e., p l ), and between the l th element and the destination (i.e., g l ).The simulation of this system is further detailed in Fig. 2 (with Λ being the number of RISs assisting transmission, i.e., depicted in Fig. 1(b), Λ = 1), where the physical model introduced in the sequel is used to simulate the signal propagated through each stage of the transmission.We then apply knowledge of the RIS functionality (i.e., reflection coefficients and phase errors that will be introduced in Section IV) in the system model to obtain the received SNR at the destination.
In our physical model, for a single RIS case, we consider the optional presence of dominant signal paths between the source and the RIS and then between the RIS and the destination.Accompanying this is a scattered signal contribution which can be modeled as a zero-mean complex Gaussian random process [12].As we will show later, for RISs, the scattered signal may exhibit non-isotropic transmission and reception behavior.It is also noted that each stage of the cascaded channel can be denoted in vector form as p = [p 1 , . . ., p l , . . ., p L ] T and g = [g 1 , . . ., g l , . . ., g L ] T .Following from this, the l th illuminated reflective elements of both p and g are characterized as complex non-zero mean Gaussian RVs due to the superposition of the dominant and scattered signal contributions within each cascaded stage of the source-to-destination link.Therefore, the magnitudes of the l th illuminated reflective elements of both p and g follow the Rician distribution [40].
The complex signal envelope at the destination is denoted as S = R exp(jΘ), where R is the received signal envelope and Θ is the received signal phase.Also, the in-phase and quadrature components are represented by X and Y , respectively.Based on this, it follows that S = X + jY , R 2 = X 2 + Y 2 , Θ = arg(X + jY ), X = R cos(Θ) and Y = R sin(Θ).From Fig. 1(a), sub index λ = 1 refers to the path between the source and the RIS for each illuminated reflective element (i.e., p l ).Correspondingly, sub index λ = 2 refers to the path between the RIS and the destination for each illuminated reflective element (i.e., g l ).Letting the quadrature components, which make up the signal contribution from each illuminated reflected element, be represented by complex non-zero mean Gaussian RVs, the received signal power at the destination can now be expressed as where L is the number of illuminated reflective elements, and B λ,l and C λ,l are mutually independent Gaussian RVs with The time-varying amplitudes of the in-phase and quadrature components of the dominant signal component are represented by u λ,l and v λ,l , respectively, with the variation related to the embedded Doppler effect [41].The Rician k λ,l factor represents the ratio between the total power of the dominant signal component, δ 2 λ,l = u 2 λ,l + v 2 λ,l , and the total power of the scattered signal component 2σ 2 λ,l , with where rλ,l = E[R 2 λ,l ] denotes the root mean square (rms) of the signal envelope.Based on this, it follows that k λ,l = , while by defining ϖ λ,l = arg(u λ,l + jv λ,l ) as a phase parameter we have, Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.If only one reflective element is illuminated (i.e., L = 1 and the l subindex is omitted), the model in (1) simplifies to Let us now consider the case where there are Λ RISs working cooperatively, in order to trasnmit the signal from the source to destination.To this effect, the received signal power at the destination can now be expressed as ( The number of illuminated reflective elements at the first RIS is L and the channel vector from the last RIS to the destination is denoted by g T , with M being the number of illuminated reflective elements at the last RIS. 1 The channel matrices between the Λ cooperating RISs are represented by H λ whilst the matrix dimensions in each H λ correspond to the number of illuminated reflective elements at each RIS.Hereafter, when considering two RISs working cooperatively (i.e., Λ = 2) the λ subindex within H λ is omitted to become H.For instance, in Figs.1(b) and 2 (using Λ = 2), with p and g as previously defined, therefore , represents the channel matrix between the two cooperating RISs assisting transmission [21].Each element in all of the vectors and matrices are complex non-zero mean Gaussian RVs, corresponding to Rician fading.Also, when only one reflective element is illuminated at each RIS (i.e., L = M = 1 and therefore g, H λ and p each reduce to a single complex Rician RV, respectively), the model in ( 5) simplifies to

A. Products of Rician Channels
As a consequence of the physical models presented for RIS-assisted transmission, we now present new results for the product of Rician channels.This is predicated on the received signal models for the product of two and then Λ + 1 complex Rician channels, as shown in ( 4) and ( 6), respectively.The product of multiple RVs, and their corresponding received signal distributions, can be used to model cascaded fading channels, such as those encountered in multihop cooperative communications and multiple-input-multiple-output (MIMO) keyhole communication systems [42].
1) First-Order Statistics: Considering the mathematical model in (4), the PDF of the received signal envelope, R, for the product of two complex Rician channels is given by [43, Eq. ( 19)], which is fully convergent.Using this and following from the signal model given in (6), the PDF of the received signal envelope, R, for the product of Λ + 1 complex Rician channels can be obtained by Theorem 1 below.
, the PDF of the received signal envelope of the model in (6), for the product of Λ + 1 complex Rician channels is given by, where Proof: See Appendix A. It is worth highlighting that the PDF of the received signal phase for the product of two complex Rician channels has not been determined in the literature.Yet, knowledge of the received signal phase is essential for the accurate characterization of complex channels, given in [44] and [45] and the design of phase-based modulation schemes, as articulated in [46].Based on this, the corresponding phase PDF is given by Theorem 2, as shown below.
Theorem 2: For the PDF of the received signal phase of the model in ( 4), which considers the product of two complex Rician channels, is given by Proof: See Appendix B. The PDF of the received signal phase of the product of Λ + 1 complex Rician channels is now obtained by Theorem 3 as shown below.
Theorem 3: For {k 1 , . . ., k Λ+1 } ∈ R ≥0 , {r 1 , . . ., rΛ+1 } ∈ R >0 , {ϖ 1 , . . ., ϖ Λ+1 , θ} ∈ [−π, π], the PDF of the received signal phase of the model in (6), which considers the product of Λ + 1 complex Rician channels is given by, where f R1,...,RΛ+1 (r 1 , . . ., r Λ+1 ) is the product of the PDF of the individual signal envelopes of the Λ + 1 complex non-zero mean Gaussian RVs.Notably, the convolution of the marginal distribution of the individual signal phases over the interval [−π, π], given the signal envelope, f Θ|R1,...,RΛ+1 (θ|r 1 , . . ., r Λ+1 ), follows the Von Mises distribution [43, Eq. ( 4)] such that Proof: See Appendix C. Several propagation scenarios arise from the model in (4), which have a significant impact on the properties of the received signal phase.For illustrative purposes, and without loss of generality, ϖ λ is considered to be time-invariant for all examples presented in this section.Figure 3 illustrates the form of the PDF of the received signal phase for a number of propagation scenarios of interest.Firstly, as expected, we observe that when k 1 = k 2 = 0, the received signal phase follows the circular uniform distribution.Most strikingly, when only one of the Gaussian RVs has a zero mean, it is found that the resultant received signal phase is still uniformly distributed, which is consistent with the result presented in [47].This observation suggests that the clustering of the signal around a particular phase, which is introduced by the dominant signal component in one stage of the cascaded transmission, is diffused in the other stage of transmission.When both Gaussian RVs have non-zero means, both Rician k factors (i.e., k 1 and k 2 ) control the concentration of the distribution around the maximum, which occurs at θ = ϖ 1 + ϖ 2 , while the minima is located at θ = ϖ 1 + ϖ 2 + π.
2) Second-Order Statistics: The backbone of our proposed simulation framework is the use of AR modeling to induce the time variant properties of the fading into the complex Gaussian RVs underlying each of the cascaded stages of transmission.Central to the construction of the appropriate AR model for channel simulation in our proposed framework Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. is knowledge of the associated complex ACF.This requires expansion of the physical model to include other propagation related effects not apparent within the first-order statistics.Initially, let us consider the product of two complex Rician channels given in (4).This is representative of the case when only one reflective element is illuminated for a single RIS-assisted transmission.For ease of understanding, and to remove any potential ambiguity, we let δ λ (t) represent the complex amplitude of the dominant signal component, and ϱ λ (t) account for the scattered signal component in each of the complex Gaussian RVs in p and g. 2 The time dependent complex received signal for a single Rician channel is already know as in [48].We hereafter extend this result to the product of two Rician channels (in this analysis when a single reflective element is illuminated for the single RIS-assisted transmission case), which can be expressed as where rλ and k λ have been defined previously.When λ = 1, (13) represents the first stage of transmission which occurs between the source and the RIS (i.e., p), and when λ = 2, this denotes the second stage of transmission which occurs between the RIS and the destination (i.e., g).In addition, we let f δ λ and α δ λ represent the relative Doppler frequency and the DoA of the dominant signal component, respectively.
Using [49], the dominant signal component, δ λ (t), is expressed as exp (j (2πf δ λ t cos (α δ λ ) + φ 0λ )) where, φ 0λ , is the initial phase, uniformly distributed on the interval [−π, π].Now by adopting the notation provided in [50, p. 82] a mathematical reference model for the scattered signal contribution, ϱ λ (t), can be written as where Q is the number of propagation paths, f D λ and f Aλ are the maximum Doppler frequencies of the signal departing (denoted with subscript D) from the source and arriving (denoted with subscript A) at the RIS (i.e., f D 1 and f A1 ), and departing from the RIS and arriving at the destination (i.e., f D 2 and f A2 ), respectively.The Doppler shifts observed in an RIS system can occur for multiple reasons, such as the mobility of the source and destination, mobile scatterers in the local environment, or some combination thereof.The random DoD of the q th propagation path with reference to the respective velocity vector is α q,D λ , while α q,A λ is the random DoA.
The variable φ q λ is a random phase uniformly distributed on [−π, π] and is independent of α q,D λ and α q,A λ for all q.This reference model allows greater flexibility in being able to consider either, or both ends of the transmission, to be in motion [50, p. 82].
It is noted here that since the RISs can be deployed in urban environments to assist transmission to users in blind zones [10], the scattered signal component may be subject to spatial filtering and as a consequence be non-isotropic in nature [51].In this case, the distribution of the DoD or DoA of the scattered signal component can be modeled using the versatile Von Mises distribution [52].This distribution can be used to model propagation scenarios in which the DoD or DoA of multipath waves are either isotropic or non-isotropic.The Von Mises distribution for the DoD (or DoA) of the scattered signal contributions, f Aλ (α λ ), is given by [52] where α λ ∈ [−π, π], κ λ ≥ 0 controls the spread of the DoD (κ D λ ) and the DoA (κ Aλ ).The mean DoD or DoA are accounted for by ᾱDλ and ᾱAλ , respectively, and each can take values in the range [−π, π].When κ λ = 0, (15) reduces to the circular uniform distribution and the signal received is the result of isotropic scattering.As κ λ increases, the DoD or DoA becomes increasingly unidirectional with κ λ → ∞ describing extremely non-isotropic scattering.As discussed in [52], a useful estimator of the spread of the DoD or DoA in polar coordinates for increasing values of κ λ (i.e., κ λ > 1) can be obtained by considering the inflection points of f Aλ (α λ ) which are approximately equal to ±1/ √ κ λ .For highly nonisotropic scattering, e.g., for κ λ = 3 and 5, the spread of the DoD or DoA of these contributions can be estimated as 2/ √ κ λ , giving 66 • and 51 • , respectively.Using the previous definition of the dominant signal component, δ λ (t), and the scattered signal component, ϱ λ (t), the complex ACF with respect to the time lag, τ , of the model given in ( 13) is now presented by Theorem 4 below.
Theorem 4: For the complex ACF of the model for the product of two Rician channels in ( 13) Fig. 4. The normalized complex ACF of a single RIS system with one illuminated reflective element (lines) in ( 16) alongside simulation results using an AR model order of 200 (shapes) for varying k 1 and k 2 , with can be expressed as, Proof: The proof is provided in Appendix D. The derived complex ACF in ( 16) is a general expression for both isotropic (i.e., κ D λ = 0 or κ Aλ = 0) and non-isotropic (i.e., κ D λ > 0 or κ Aλ > 0) scattering conditions that can occur between the source and the RIS, and then between the RIS and the destination.It is also recalled that is often more convenient to use the complex ACF in its normalized form, given by φSS (τ ).This can be straightforwardly obtained by setting rλ = 1 in (16).
As an example, now let us consider the real and imaginary parts of the complex ACF for the case where a source transmits to an RIS, illuminating a single reflective element, and the RIS reflects the signal to the destination.The real part of the normalized ACF for this type of scenario is given in Fig. 4(a).It can be seen that when either k 1 or k 2 are decreased, there is weaker correlation, as compared to the case where there is at least one strong dominant signal component present (e.g., k 1 = 5 and k 2 = 0.8).Next, we compare the real part in Fig. 4(a) to the imaginary part in Fig. 4(b), which is a measure of the cross-correlation of the real and imaginary parts of the complex received signal [50].We observe that in the presence of stronger dominant signal components, there is greater correlation between the two quadrature components in this particular example, where the scattering conditions have remained identical.
Up until now we have only considered the product of two Rician channels.Nonetheless, this can readily be extended to the more general case of the product of Λ + 1 Rician channels by changing the upper limit of the product in ( 16).
The complex ACF in this case is given by Corollary 1 as shown below.
Corollary 1: For the complex ACF of the model for the product of Λ + 1 Rician channels in ( 6) can be expressed as, which reduces to (16) for the special case that Λ = 1.
A summary of the parameters used in ( 17) is given in Table I.

III. SIMULATION OF THE COMPLEX RECEIVED SIGNAL
Underlying our simulation framework is the well-known approach of AR modeling, which we use to introduce the necessary time-related variation in the complex received signal.AR modeling can be used to exploit the autocorrelation properties of a signal and has been successfully used to simulate wireless channels based on the theoretical model from [35], and the empirical model from [53] for the channel's complex ACF.In what follows, we provide an overview on how to generate p, H 1 , . . ., H Λ−1 , and g over time.To assist with the understanding of the procedure, a block diagram is provided in Fig. 5.The simulator framework developed can be found here [54].

A. AR-Based Simulation
An AR model is an infinite impulse response filter in which the current output is dependent on the previous inputs.In order to filter the complex white noise, the AR parameters must be estimated, whereas the number of AR parameters is equal to the AR model order.It can be seen in Fig. 5 that the complex ACF of the scattered signal contribution is used as the input Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.(17) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
for this procedure.To estimate the AR parameters, the complex ACF of the scattered signal contributions in ( 35) is used to establish the Yule-Walker equations, which are then solved (see [48] for further mathematical details of this process).For each of the scattered signal contributions reflected from each of the illuminated reflective elements in p, H 1 , . . ., H Λ−1 , or g, complex white noise must be filtered using the AR parameters estimated for the different complex ACFs of p, H 1 , . . ., H Λ−1 , and g.This technique ensures that the scattered signal reflected from each of the illuminated reflective elements at the RISs in p, H 1 , . . ., H Λ−1 , and g are all autocorrelated.

B. Spatial Correlation
Using a single RIS system as an example, we will now detail the approach used to generate spatial correlations between the scattered signals reflected by the illuminated reflective elements of the RIS(s).To achieve this, we create the departing and arriving spatial correlation matrices, denoted by Φ D λ and Φ Aλ , respectively.As shown in Fig. 5, these correlation matrices are used as inputs.In a single RIS system, when λ = 1, Φ D 1 dictates the spatial correlations introduced by the placement of antennas (or single antenna) at the transmitting source.Reciprocally, Φ A1 , represents the spatial correlations between the illuminated reflective elements at the RIS.The spatial correlation matrix, Φ DA1 can be expressed as the Kronecker product of both sets of spatial correlations such that Φ DA1 = Φ D 1 ⊗ Φ A1 .To generate the spatial correlations between the complex Gaussian RVs, a Cholesky lower decomposition is performed on Φ DA1 to generate w 1 .This is then multiplied on the vectorized Gaussian RVs produced by the AR-based simulation, as described in Section III-A.This procedure must be repeated to introduce the desired spatial correlation between the signal reflected from each of the illuminated reflective elements at the RIS and propagated towards the destination (i.e., λ = 2).The spatial correlation matrix of the arriving signals at an RIS will be the same as the spatial correlation matrix of the signals departing the same RIS (i.e., Φ A1 = Φ D 2 ), even though they correspond to separate stages of the transmission process.This is because the illuminated reflective elements used to redirect the signal are in the same position and, thereby, subject to the same spatial correlation properties.This assumption is based on the physical arrangement of the illuminated reflective elements in an RIS, although the simulation is not limited to this.In other words, the spatial correlation matrices of the arriving and departing signal can be set such that Φ A1 ̸ = Φ D 2 .As a result, the simulation framework proposed here can easily be adjusted for other applications, such as relays using adaptive antenna array structures, where different elements may be used for transmission and reception.

C. Complex Signal Formation
As shown in Fig. 5, the next stage of the simulation is the complex signal formation.The RVs representing the time-correlated scattered signal reflected from each of the illuminated elements now need to be altered to have the correct power and include a dominant component (if present).This is straightforwardly achieved by adding u λ and v λ from (3) to the zero mean in-phase and quadrature components respectively, which have been multiplied by σ λ in (2).This stage, and the two previous stages of the simulation framework are repeated until p, H 1 , . . ., H Λ−1 , and g are all generated, representing each stage of the communication link in RIS-assisted transmission.

D. Output
Lastly, depending on the RIS system (single or cooperative) p, H 1 , . . ., H Λ−1 , and g representing each stage of the RIS-assisted transmission are multiplied together to form the final complex received signal.

IV. AN RIS SYSTEM EXAMPLE
As an illustrative example of the use of our simulation framework within an RIS system, we now consider the cases of both single and cooperative RIS-assisted transmission.We first consider communication from a source to a destination through a single RIS (i.e., Λ = 1) having L illuminated reflective elements, as illustrated in Fig. 1(a).The signals before (i.e., p l ) and after (i.e., g l ) the RIS are assumed to be Rician distributed, as previously defined in Section II, and hence p l and g l are each modeled as complex non-zero mean Gaussian RVs.Based on this, the received signal at the destination using a single RIS with L illuminated reflective elements can be written as where x is the signal transmitted by the source satisfying E |x| 2 = P , with P being the transmit power and n ∼ CN 0, σ 2 N denoting the additive white Gaussian noise (AWGN), with zero mean and variance σ 2 N .Furthermore, η l represents the reflection coefficient introduced by the l th illuminated reflective element of the RIS, whilst ϑ l represents the phase shift introduced by the l th illuminated reflective element of the RIS.Note that the received signal in (18) can also be expressed as [12] where p and g are defined as previously, and is a diagonal matrix which captures the reflection functionality of the RIS elements.
Following from this, we now consider communication from a source to a destination through two RISs working cooperatively (i.e., Λ = 2), as illustrated in Figs.1(b) and 2. We assume that the first RIS has L illuminated reflective elements, while the second RIS has M illuminated reflective elements.As before, we assume that there is no direct link between the source and destination due to severe blockage effects.The signal propagating between individual reflective elements is assumed to be Rician distributed as detailed in Section II, and hence p l , h ml and g m are each modeled as complex non-zero mean Gaussian RVs.Based on this, the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.received signal at the destination through two cooperative RISs having L and M illuminated reflective elements, respectively, can be expressed as where g, H, Θ, p, x, and n are as defined previously, whilst is a diagonal matrix which captures the reflection properties of the M illuminated reflective elements of the second RIS. 3he reflection coefficients introduced by the l th and m th illuminated reflective element of the first and second RIS, are represented by η l and η ′ m , respectively.Likewise, ϑ l and ϑ ′ m represent the phase shift introduced by the l th and m th illuminated reflective element of the first and second RIS, respectively.

A. Phase Errors
Given the phases arg(p l ), arg(g l ), arg(g m ) and arg(h ml ), the received SNR can be maximized by intelligently controlling the reflective properties of each illuminated reflective element within the RIS.In other words, the optimal choice of ϑ l for the single RIS system is obtained when they are adjusted so as to cancel the overall phase shift arg(p l ) + arg(g l ), whilst the optimal choice of ϑ l and ϑ ′ m , for the cooperative RIS system, is obtained when they are adjusted so as to cancel the overall phase shift arg(p l ) + arg(h ml ) + arg(g m ).Now, assuming the phase shifts induced by the channels are not estimated perfectly and/or the desired phases cannot be set precisely, we consider the deviation of ϑ ′ m and ϑ l from the ideal setting through the phase noise θ′ m and θl terms, respectively, which are uniformly distributed between [−π, π].The received signal for the single and cooperative RIS systems can now be expressed as and respectively, where Following from ( 21) and ( 22), the received SNR for the single and cooperative RIS systems can now be written as and respectively, where γ = P/σ 2 N is the average transmit SNR.Alternatively, (23) can be expressed as whereas (24) can be expressed as Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where the complex non-zero mean Gaussian RV magnitudes, V. PERFORMANCE ANALYSIS AND RESULTS In this section, we use the proposed simulation framework with the system model described above to study the performance of D2D communications when assisted by single (i.e., Λ = 1) and then cooperative (i.e., Λ = 2) RIS systems (see Fig. 1).As well as presenting the OP for both systems, we also demonstrate the utility of the simulation framework by reporting important second-order statistics, such as the AOD and LCR.To the best of the authors' knowledge, the offered results have not been previously reported in the open literature.We now provide details of the simulations for both communication systems.
Without loss of generality, it is assumed that the transmitting and receiving devices are equipped with a single directional antenna.Also, both systems operated in an urban environment, resulting in non-isotropic behavior of the scattered signals.
At each stage of the communication link for both systems, we considered a dominant signal component with f δ λ = 0 Hz and α δ λ = 0 rad.The rest of the parameters are summarized in Table II.It should be noted that f D 1 , f A1 and f A2 have values which correspond to both the users in the D2D communication link being in motion.All other f D λ , f Aλ Doppler frequencies, which are below 1 correspond to the stationary RISs.The proximity of the illuminated reflective elements in the single and cooperative RIS systems results in spatial correlations, such that each entry in the associated systems spatial correlation matrices is equal to 0.9 (i.e., Φ A1 and Φ D 2 for the single RIS system, and Φ A1 , Φ D 2 , Φ A2 and Φ D 3 for the cooperative RIS system).In all simulations, an AR model order of 200 was used, with a sampling frequency of 1 kHz to generate 2 × 10 6 samples using a bias of 1 × 10 −3 whilst the noise power of the AWGN was set to σ 2 N = 1.

A. Outage Probability
The OP, P OP (γ th ), is defined as the probability that the instantaneous SNR, γ, falls below a certain threshold, γ th , namely For the figures presented in this subsection, we assume γ th = 5 dB. Figure 6 depicts the OP versus γ for both scenarios with and without phase errors, and for different values of η l and η m .As expected, the OP is lower for higher Here, ϑe = θl for the single RIS system and ϑe = θl + θ′ m for the cooperative RIS system.Fig. 7. OP versus γ for the single RIS system, for an increasing number of illuminated reflective elements, with and without phase errors for η l = 0.4.
values of the reflection coefficients and the OP is greater when there were phase errors present.Figure 7 shows the variation in the OP versus γ for the single RIS system, with an increasing numbers of illuminated reflective elements, L, with and without phase errors, for η l = 0.4.We observe that the OP decreases as the number of illuminated reflective elements increases, since there is a greater signal contribution re-radiated towards the receiving device.As expected, the OP is lowest in the absence of phase errors. 4

B. Average Outage Duration
The AOD, T(γ th ), measured in seconds, is used to describe how long on average a wireless system remains in outage [31], 4 While beyond the scope of this paper, we make the point that the proposed simulation framework can be used to test the performance of various signal processing techniques designed to be used with RIS systems.For example, approaches to phase optimization such as that proposed in [55].
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Hz for the cooperative RIS system.Here, ϑe = θl for the single RIS system and ϑe = θl + θ′ m for the cooperative RIS system.[32].Mathematically, the AOD is given by [56, Eq. (1.29)] where P OP (γ th ) is the OP and N (γ th ) is the LCR, defined as the rate at which the SNR crosses a given threshold level, γ th , in a positive (or negative) direction per unit time.
For both the single and cooperative RIS systems, when the average transmitted SNR is at −7.5 dB and below, and in the presence of phase errors with a reflection coefficient of 0.5, the AOD is higher compared to when there are no phase errors present.The effect of varying the maximum Doppler frequency of the scattered signal contributions is depicted in Fig. 9.It is evident that when no phase errors are present in either system, the increase in the Doppler frequencies causes a decrease in the AOD, as expected, due to the increased temporal variation.

VI. CONCLUSION
For the first time, this paper has provided a method for simulating a cooperative RIS system where the received signal is autocorrelated, spatial correlations can exist between the scattered signal reflected from each of the illuminated reflective elements at an RIS, and phase errors can be present.To this end, physical models for both single and cooperative RIS-assisted transmission have been presented, with the fluctuation of the signal transmitted between the reflective elements assumed to follow that described by the Rician fading model.Based on the physical models, novel formulations for the PDF of the received signal phase were provided for the product of two and multiple complex Rician channels.A number of examples have been included for the PDF of the received signal phase alongside respective simulated results to demonstrate its validity.To facilitate the simulation of cooperative RISassisted transmission, the complex ACF of the signal reflected by each illuminated reflective element was obtained, taking into account non-isotropic transmission and Doppler effects caused by mobility and/or the environment.Thanks to the physical model being fully defined in terms of underlying Gaussian RVs, introduction of temporal correlation within our framework was achieved using AR-based simulation, utilizing the derived complex ACF.As a byproduct of the physical placement of the illuminated reflective elements in an RIS, spatial correlations are likely to exist between the scattered signal at each stage of the cascaded link and these were implemented in our framework using the Kronecker product.
Moreover, example simulations using our new framework were conducted to gain insights into the OP, LCR, and AOD for RIS-assisted D2D communications.It was shown that the OP decreases for higher numbers of illuminated reflective elements in a single RIS system, and the OP further decreases when two RISs cooperatively assist transmission.The importance of including Doppler effects into the simulations of RIS systems was examined, where an increase in the maximum Doppler frequency of the scattered signal contribution was shown to noticeably decrease the AOD when no phase errors were present.This enhancement was further increased in cooperative RIS-assisted transmission.However, quite surprisingly, such effect was not observed when phase errors were present.Finally, it is worth remarking that the proposed model and simulation framework is not just limited to RIS systems.It can be used to simulate a wide variety of wireless systems such as relay or metasurface assisted communications, D2D and pointto-point communications which make use of antenna arrays and cascaded links.

APPENDIX A PROOF OF THEOREM 1
We first consider the mathematical model given in (6), which is the product of Λ + 1 complex non-zero mean Gaussian RVs.The PDF of the received signal envelope can be determined by extending the principle used in [43,Theorem IV.3] to Λ+1 complex non-zero mean Gaussian RVs.The resultant expression of the PDF of the received signal envelope for the product of Λ + 1 complex Rician channels is shown in (8).

APPENDIX B PROOF OF THEOREM 2
The PDF of the received signal phase for the product of two complex Rician channels is found by integrating the joint envelope-phase distribution with respect to R, the received signal envelope.The joint envelope-phase distribution of the product of two independent non-zero mean complex Gaussians is already known as [43, Eq. ( 6 (from (2) with the l subscript now removed) into (31) and then integrating the result with respect to the received signal envelope by sequentially using the identities [39, Eq. (03.02.06.0037.01)]and [38,Eq. (6.561.16)],results in (10).

APPENDIX C PROOF OF THEOREM 3
The PDF of the received signal phase for the product of Λ + 1 complex Rician channels is obtained by following the same procedure as presented in [43,Appendix A].Now it is extended to Λ + 1 complex non-zero mean Gaussian RVs and integrated with respect to the received signal envelope to obtain (11).This completes the proof.

APPENDIX D PROOF OF THEOREM 4
The ACF of a wide sense stationary process may be written as ϕ SS (τ ) = E [S(t)S * (t + τ )], where t is the current time and τ is the time lag [50].Using (13) as a model for the complex received signal when one reflective element is illuminated and after some mathematical manipulations the complex ACF may be expressed as Evaluating the autocorrelation of the scattered signal component separately for each complex Gaussian RV gives (33), as shown at the bottom of the page.This conveniently reduces to (34), as shown at the top of the next page.Knowing that α D λ and α Aλ are independent RVs (E[α D λ α Aλ ] = E[α D λ ]E[α Aλ ]), followed by substituting the PDF given in (15) into (34) along with [38,Eq. (3.338.4)],yields (35), as shown at the top of the next page.Now substituting (35) into (32), we obtain the complex ACF given in (16), completing the proof.

Fig. 3 .
Fig.3.The PDF of the received signal phase of the product of two complex Rician channels in(10) for varying cases.The lines represent analytical results, with a maximum number of terms of b = c = d = 15, while the shapes represent the results of the Monte-Carlo simulations using 1 × 10 6 samples for each stage of the transmission.

Fig. 5 .
Fig. 5. Block diagram of the framework used to simulate signal reception in each stage of the communication link in RIS-assisted transmission.

Fig. 6 .
Fig. 6.OP versus γ for single and cooperative RIS systems in RIS-assisted D2D communication, with and without phase errors for different values of η l and ηm when L = M = 4. Grey solid and dashed lines indicate η l = ηm = 0.8, while black solid and dashed lines indicate η l = ηm = 0.5.Here, ϑe = θl for the single RIS system and ϑe = θl + θ′ m for the cooperative RIS system.

Fig. 8 .
Fig. 8. LCR (a) and AOD (b) versus γ for single and cooperative RIS systems in RIS-assisted D2D communication, with and without phase errors for different values of η l and ηm when L = M = 4. Grey solid and dashed lines indicate η l = ηm = 0.8, while black solid and dashed lines indicate η l = ηm = 0.5.Here, ϑe = θl for the single RIS system and ϑe = θl + θ′ m for the cooperative RIS system.

Fig. 9 .
Fig. 9. AOD versus γ for single and cooperative RIS systems in RIS-assisted D2D communication, with and without phase errors and η l = ηm = 0.6 when L = M = 4. Black lines represent the single RIS system and grey lines indicate the cooperative RIS system.The solid lines used f D λ and f Aλ as provided in Table II and the dashed lines used f D 1 = f A2 = 20 Hz, f A1 = f D 2 = 15 Hz for the single RIS system and f D 1= f A3 = 20 Hz, f A1 = f D 2 = f A2 = f D 3 = 15Hz for the cooperative RIS system.Here, ϑe = θl for the single RIS system and ϑe = θl + θ′ m for the cooperative RIS system.

TABLE II SUMMARY
OF PARAMETERS USED FOR SIMULATIONS