Spatially Correlated Dual Hop RIS Aided Next Generation Wireless Systems: An Outage Perspective

Design and analysis of Reconfigurable Intelligent Surfaces (RIS) assisted new wireless communication systems have attracted much attention among the 6G technology researchers. Recently, studies were made to compare the use of RIS and Decode and Forward (DF) relaying when a weak direct path exists between the transmitter and receiver. Amalgamating the merits of both RIS and DF relay technologies, an RIS assisted DF relay-based cooperative wireless system with energy harvesting is proposed in this paper. A Lambertian model of reflectance is employed in the analysis of the reflecting nature of RIS. Further, the impact of spatial correlation which arises due to the tightly-packed planar RIS array geometry is also considered for performance analysis. By varying inter-element distances and transmit power, the performance of the proposed system is analyzed in terms of outage probability and energy harvesting capability. Numerical expressions are derived for the outage probability of the proposed system and verified using Monte-Carlo simulations. Based on the high performance, the concept of joint use of RIS and DF Relay is extended to the Index Modulation (IM) based wireless system.


I. INTRODUCTION
Futuristic 6G wireless technologies will supposedly accomplish the expectations not met with 5G. In infrastructure aspect, the key enablers of 6G are envisioned to be Reconfigurable Intelligent Surfaces (RIS), ultra-massive Multiple Input Multiple Output (MIMO) communications, and user-centric networking [1]- [3]. Currently, the research focus is on investigating the customization aspect of the wireless propagation environment [4], [5]. Controllable intelligent channels, enabled by RIS, are utilized to adjust the phase of incident signals such that the signals are added constructively at the end-user. In RIS, arbitrary incident waves are phase tuned by the discrete reflecting elements, and the reflected waves out of RIS, reach the end-user as a coherent wave [6]. In this way, operators can have certain control over the random channel behavior. In [7], it is proved that the RIS aided The associate editor coordinating the review of this manuscript and approving it for publication was Barbara Masini .
communications would provide a more favorable path than the free space when the direct link is obstructed. The utilization of RIS for assisting wireless communication has numerous benefits such as facilitation of ultra-reliable link even at very low Signal to Noise Ratio (SNR), low-cost implementation, and energy-efficient hardware. Energy Harvesting (EH) is the enabling factor to improve the energy efficiency in 5G as well as in 6G wireless networks [8]. A realistic power consumption model for RIS based communication is presented in [9]. [10] investigates Simultaneous Wireless Information and Power Transfer (SWIPT) of RIS aided system. It is shown through simulations that RIS aided SWIPT improves the energy efficiency. Wireless energy transfer for RIS aided OFDM system is investigated in [11]. It is shown through theoretical results that it demonstrates ultra energy efficient operation due to RIS assist.
The performance of RIS based communication system is analyzed in detail in [12]. It is also extended to RIS assisted MIMO Index Modulation in [13]. Low power and highly VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ secure back-scatter RIS based spatial modulated system is explored in [14]. In [15], hybrid transmitting schemes are proposed for the joint operation of both RIS and Decode and Forward (DF) relaying. Achievable rates of the proposed systems are extensively analyzed and significant performance improvement is observed. In [16], a system model is proposed with two RIS placed side by side to the relay to improve full-duplex performance. A generalized hardware impairment model of RIS aided system is investigated in [17]. When comparing RIS assisted relaying system with conventional MIMO DF relaying system, RIS assisted system has the merit of using only passive elements to aid communication while the latter uses multiple energy-hungry RF chains. Performance between RIS aided communication and relaying are examined in [18]- [21]. All these analysis are based on the assumption of independent and identically distributed (i.i.d) Rayleigh fading channel. However, as the spatial samples at RIS are correlated, the number of independent paths and signal diversity is reduced [22], [23]. Spatial correlation in RIS-aided systems is a very critical factor in analyzing the performance of the systems. An accurate spatial correlation model for RIS systems is presented in [24] as a function of the array geometry of reflecting units. The proposed model in [24] is more accurate than the conventional Kronecker model and in accordance with the 3D Clarke's model [25].
In this paper, the major contributions are summarized as follows: • A novel co-operative communication system is proposed by combining the merits of RIS and EH in DF relaying system in a Non-Line of Sight (NLOS) wireless environment • The analytical expressions are derived for the outage performance of the proposed system at both hops, in the presence of spatial correlation at RIS, using the spatial correlation model in [24].
• Alternative approaches are presented to generate correlated Rayleigh envelopes to carry out the outage analysis using simulation.
• Further, the proposed system model is extended with IM at the receiver side. The outage performance of the extended system is also analyzed in the presence of spatial correlation.
A. NOTATIONS z denotes the closest integer smaller or equal to z (floor operation). V{z} is the variance of a random variable z. E{z} is the mean of a random variable z. COV{z 1 , z 2 } is the covariance between random variables z 1 and z 2 . |z| denotes the magnitude of a complex random variable z. N (µ, σ 2 ) denotes a real Gaussian random variable with mean µ and variance σ 2 .

II. SYSTEM MODEL
A three node co-operative wireless communication system is shown in Fig. 1. The proposed RIS-DF-EH system has RIS at both source and destination blocks. RIS at source block (RIS-S) is located in close proximity to the source node S and is employed as an access point as described in [12, section III]. The advantage of placing RIS in close proximity to the RF source is that the phase compensation is to be performed for the channel between RIS and desired destination only. It is also proved that RIS position must be at either closer to the transmitter or receiver for efficient operation [26], [27]. Similar technique is also used in commercial holographic beamformer [28]. Similarly, RIS at destination block (RIS-D) is placed near the destination node D. The basic motivation of the proposed system is to improve energy harvesting at the relay node for increased operational time and data rate. One of the physical scenarios in which the proposed system can be applied shown in Fig. 2, with Coordinated Multi-Point (CoMP) transmission and reception [29].
Let N A and N B be the number of elements of the planar RIS-S and RIS-D respectively. Let where N H A is number of horizontal elements and N V A is the number of vertical elements in RIS-S. Similarly N B = N H B × N V B . Further, the area of a single RIS reflecting unit is given by d H d V , where d H and d V are horizontal and vertical inter-element distances respectively. (On the account of RIS reflecting elements being deployed in an edge to edge manner, it is apparent that inter-element distances can be swapped for sides of the reflecting unit).
A plane wave with wavelength λ c corresponding to the carrier frequency f c is reflected from RIS-S with the azimuth angle of ϕ and elevation angle of θ. As the RIS element's surface reflectance behaves as isotropic and diffusely reflecting, the azimuth electromagnetic intensity obeys Lambert's cosine law. The joint Probability Density Function (PDF) of ϕ and θ is given as Angular distribution of planar RIS geometry is shown in Fig. 3. It is observed that the resultant scatter plot is  highly directional compared to i.i.d. spread. With the angular constraints from (1), the array response vector is given by [24] a(ϕ, θ) = e jk(ϕ,θ) T p 1 , . . . , e jk(ϕ,θ) where k(ϕ, θ) is the wavenumber vector. It is defined as The position vector of the elements in RIS-S is given by With isotropic scattering in the half-space in front of the RIS, (n, m) th element of correlation matrix where sinc(x) = sin(πx)/(πx) is the sinc function. The covariance matrix (R G ) is defined as where ρ gi,j are the pairwise covariance values between each complex Gaussian random variables. The real and imaginary part of individual complex Gaussian samples are i.i.d, but the correlation occurs only between the successive pairs of Gaussian random variables. There are infinitely many multipath components in an isotropic scattering environment between RIS-S and relay node R. By considering L plane waves, the channel between RIS-S and relay node R is modeled as where c l ∈ C is the complex signal attenuation of the lth component, ϕ l and θ l are the azimuth and elevation angles of l th plane wave. The coefficients c 1 , . . . , c L are i.i.d. with zero mean and variance d H d V . Using the covariance matrix R G from (5), h ∈ N C (0, d H d V R G ). It is noted that channel vector h is the independently distributed but it is spatially correlated due to RIS geometry. The complex channel coefficient h i in polar form is represented as h i = α i e −jθ i , α i is the magnitude of the channel coefficient h i and modeled as Rayleigh The relay node R is assumed to be operating in DF mode with EH capability. The M-ary digital modulated RF signal from the source node S is transmitted to the relay node R via RIS-S. Then, the relay node R decodes and re-transmits the signal to the destination node D through RIS-D. Let φ i be the phase shift introduced by the i th reflecting element of the RIS-S. The active phase shifting vector of the RIS-S is defined as The received signal at the relay node R is given by P T is the transmitted symbol power and large scale channel gain β PL1 is a function of propagation distance between RIS-S to R. For the NLOS wireless propagation scenario, β PL1 defined as [30, at the relay node R. G t and G r in dB are the antenna gains of source node S and relay node R.
The relay node R is equipped with a power splitter to co-ordinate EH and information processing. More specifically, the power splitting factor λ PS is utilized for EH and (1−λ PS ) for information processing. After energy harvesting, the relay node forwards the decoded symbol to destination node. Hence, from (8) the instantaneous SNR at relay node R during the information processing time is given by, Considering ideal RIS operation i.e total elimination of the phase component of the channel coefficients at the relay. RIS-S is tuned in a way such that for every Then (10) can be written as, α is the sum of correlated Rayleigh variables with covariance matrix R R . Following Central Limit Theorem (CLT), the mean ofα is expressed as . The convergence is tested and elaborated subsequently. According to CLT,α approximates to N (µα, σ 2 α ) for inter-element distances greater than or equal to 0.5λ c and for higher values of N A . This is confirmed by depicting Quantile-Quantile(QQ) plot in Fig. 4 for the inter-element spacing Fig. 5. It is observed that the left tail doesn't follow the standard normal distribution quantile line and it is confirmed that it does not converge to normal distribution.  As the channel coefficients are spatially correlated, the variance ofα is determined from R R . Using [31, (9)], the covariance matrix R R is given by Substituting (5) in (12), the elements of covariance matrix (R R ) is defined as n, m = 1, . . . , N A . In a matrix form R R is written as Proposition 1: Variance of sum of correlated random vari- Since the covariance matrix R R is square semi-definite and symmetric over its principal diagonal, the trace of the covariance matrix R R is equal to V (X i ) and sum of the upper triangular elements is equal to COV X i , X j . Similarly, the sum of the lower triangular matrix also equals to COV X i , X j . The overall sum of the elements of the covariance matrix is equal to Thus, the variance of sum of correlated random variables In vector form, it is simply written as Corollary 1: In an ideal isotropic scattering environment, all elements are uncorrelated and R G becomes an identity matrix I. Applying CLT, N (µα, This case holds only when RIS elements are deployed in a Uniform Linear Array (ULA) geometry with spacing in multiples of half wavelength. In the proposed system, both the RISs are deployed in rectangular array geometry. As the spacings between the diagonal elements are not exact multiples of half wavelength, spatial correlation always exists. The resultant variance is typically lesser than tr( Corollary 2: When the elements are fully correlated, R G is equal to ones square matrix. This is practically impossible due to the fact that (5) For notational convenience, letα 2 = A χ 2 . The PDF of A χ 2 is a non-central chi-squared distribution with one degree of freedom and non-centrality param- . Using R R non-centrality parameter 1 is written as [32], In order to represent spatial correlation, it is modified as ρ ri,j = η sc . It represents spatial correlation term. For completely uncorrelated i.i.d, η sc = 0. Substituting in [33, (1.6)], the PDF of A χ 2 is written as The outage occurs due to uncertainty of the channel conditions and lower received SNR than the predefined threshold (γ th ). Using (11), the outage probability at relay node R is given by, (20) Using [34], the outage probability is expressed as where Q M (., .) is the Marcum-Q-function. Using [35], it is expressed as Time switching policy of energy harvesting and information processing is shown in Fig. 6. ρ is the time proportion for VOLUME 9, 2021 Time Switching Relay(TSR) protocol and T is the block duration for both EH and information processing. During the time slot ρ 2 T , the harvested energy at relay node R is given by, where η is the energy efficiency factor. The power available for information transmission at relay node R is given by The received signal at the destination node D is given by, x is the decoded and re-transmitted symbol from the relay node R. Large scale channel gain β PL2 is a function of propagation distance between R to RIS-D. g = [g 1 , . . . , g k , . . . , g N B ], where g k is the channel coefficient from relay node to the k th reflecting element of RIS-D. ψ k = [ψ 1 , . . . , ψ k , . . . , ψ N B ], ψ k is the phase shift introduced by the RIS-D. The magnitude component is given by |g k | = β k , also for ideal RIS operation the phase component is nullified i.e { g k } = 0, ∀ k and n D is the receiver AWGN at destination node. The instantaneous SNR at destination node D is given by . B χ 2 follows a non-central chi-squared distribution with one degree of freedom and non centrality parameter 2 = η sc +0.25d H d V N B (π (N B − 1) + 4). Substituting (24) in (26), the instantaneous SNR at destination node is expressed as The outage probability at the destination node D is given by, Substituting (27) in (28), the outage probability at the destination node P D out is expressed as, Grouping the product of RVs A χ 2 , B χ 2 , (29) is rewritten as, The PDF of the product of two non central chi-square distributions A χ 2 and B χ 2 , is given by [36], [37] κ 1 and κ 2 are the degrees of freedom of the A χ 2 and B χ 2 respectively. K g (z) is the modified Bessel function of the second kind. The exponentially decaying function is denoted by If equal number of reflectors are allocated for source and destination RISs, then N A = N B = N , 1 = 2 = and κ 1 = κ 2 = 1. The expression for the PDF is simplified as The outage probability at destination node D is given by where β PLi ] A practical approach to solve (49) numerically using simulation software is given in Appendix A.

C. END TO END P S→D out
An outage occurs if the end-to-end SNR of any link i.e.,γ S→R or γ R→D falls below an outage threshold γ th . The overall outage probability at the destination node P S→D out (γ th ) is written as,

IV. RIS AIDED DF RELAYING INDEX MODULATED SYSTEM WITH EH (RIS-DF-IM-EH)
In this section, RIS Aided DF Relaying Index Modulated System with EH (RIS-DF-IM-EH) is proposed as an extension to improve the spectral efficiency of the proposed system in section III. In conventional Index Modulation (IM) system, the indices of transmit antenna are used to provide additional information along with an M-ary modulated carrier to increase spectral efficiency [38], [39]. The proposed IM system is inspired from RIS-IM scheme involving receive antenna index to transmit additional spatial bits as proposed in [13, Section II-B]. The system model is shown in Fig. 7. One of the physical scenarios in which the proposed IM-DF-RIS relay system can be applied is shown in Fig. 8 with CoMP transmission and reception. In the first hop of the proposed RIS-DF-IM-EH, the information bits from the source node S are clustered into two groups log 2 N R (Spatial bits) and log 2 M (Data bits). M denotes to the number of digital modulation symbols per transmission. The first group of log 2 N R adjusts the RIS-S phases according to the selected receive antenna at relay node R with index p.
The spatial bits are sent to RIS-S controller system via the bit-splitter and the reflection parameters are appropriately changed such that the signal reaches the desired antenna at relay node R. The second group of log 2 M is transmitted using single RF chain from source node S.
The received signal at the p th receive antenna of the relay node is written as h p,i is the channel coefficient from i th reflecting element of RIS-S to the p th receive antenna at the relay node. The complex channel coefficient in polar form is, . φ i is the phase shift introduced by the i th reflecting element of the RIS-S, x is the transmitted symbol and n R is the received AWGN with N (0, N 0 ) at the relay node. Similarly to the previous case, the relay is equipped with a power splitter to co-ordinate EH and information processing. More specifically, the power splitting factor λ PS is utilized for EH and (1 − λ PS ) for information processing. For further clarity, Maximum Likelihood (ML) detection for the proposed RIS-DF-IM-EH operation involves joint decoding of both antenna index p and the symbol x at relay node R is given as From (36), instantaneous SNR at the RIS-DF-IM-EH relay node R is given by p denotes the selected antenna index.

A. OUTAGE PROBABILITY AT RELAY (P R out (γ th ))
The index modulated channel between the source nodes and relay node R consists of N R independent paths. Ideal operation of RIS-S is considered to eliminate all the phase components of the channel. Hence, the PDF of the IM channel is modeled by substituting κ = N R degrees of freedom in the PDF from (19) of non-central chi-squared distribution for IM case is given as, where The probability of outage at the RIS-DF-IM-EH relay node is given by the P R out (γ th ) = Pr γ R inst ≤ γ th . It is determined as

B. OUTAGE PROBABILITY AT DESTINATION NODE D
In the second hop of the proposed RIS-DF-EH-IM, the relay node R forwards the decoded bits both spatial and data bits to the destination node D using the harvested energy. Similar to the first hop, the first group of log 2 N D adjusts the RIS-R phases according to the selected receive antenna at destination node D with index q. The spatial bits are sent to RIS-R controller system via the bit-splitter and the reflection parameters are appropriately changed such that the signal reaches the desired antenna at destination node D. The second group of log 2 M is transmitted using single RF chain from relay node R. RIS-DF-IM-EH relay node R harvests energy and then re-transmits the signal to destination with the aid of decoded bit splitter to tune the RIS-R. Single transmit antenna at relay node R is used in transmission of the modulated signal to the RIS-R elements placed in the near proximity of the relay node R (similar to Wireless Access Point (WAP) operation) [12], [13]. RIS-R reflects the beam to the desired destination antenna index q. The received signal at destination node D is given by The instantaneous SNR at the destination node D is expressed as As described in previous section, substituting (24) in (43) the instantaneous SNR is written as and B q χ 2 are non-central chi-squared distribution with degrees of freedom N R and N D and non-centrality parameters 1 and 2 respectively. The PDF of the product of A p From (45), the probability of outage P D out at destination node D is written as

V. RESULTS AND DISCUSSIONS
In this section, the performance of the proposed systems is analyzed using derived analytical expressions and simulations parameters listed in Table 1. It is assumed that d H = d V and same at both RIS-S and RIS-D. Using (5), the variations in spatial correlation between the first and second elements are shown in Fig. 9 with respect to the change in inter-element spacing. For simulation purpose, two algorithms are given  in Appendix B to generate correlated Rayleigh envelopes. It is observed that both the algorithms exhibit perfect match between the theoretical and simulation values of spatial correlation. It is also verified using numerical integration of E a(ϕ, θ)a(ϕ, θ) H [24, (9)]. The noise variance N 0 in dBm is modeled through [40] N 0 = −174 + 10 log 10 ( and P T is assumed to be unity. The outage performance of the proposed system at the relay node R is shown in Fig. 10 with respect to SNR (in dB) and spatial correlation coefficients. It is inferred that outage probability increases as spatial correlation increases since it introduces severe degradation in signal diversity. The outage performance at the relay node R of the proposed RIS-DF-EH system is shown in Fig. 11 with respect to SNR (in dB) and inter-element distances. From the color map plot, it is inferred that simultaneous increase in SNR and the inter-element   spacing substantially improves the outage performance of the proposed system. The dark blue quadrant refers to no outage zone, i.e., the capacity exceeds the threshold of 10 bits/s/Hz and the yellow area in the bottom left depicts a full outage region. As expected, tight packing (high spatial correlation) and low SNR level cause severe degradation in outage performance. The ergodic capacity at the relay node R is shown in Fig. 12. A maximum of 19.6 bits/s/Hz is achieved at 25dB with d H = d V = 0.7λ c .  The outage performance and ergodic capacity at the destination node of the proposed system are shown in Fig. 13 and Fig. 14 with respect to SNR (in dB) and inter-element distances. The maximum capacity is around 30 bits/s/Hz at 25dB SNR and 0.7 inter-element distance. This due to fact that the signal experiences two independent RIS phasing operations. The energy harvesting capability of the proposed systems is shown in Fig. 16. When the power splitting factor  and d H = d V are increased, the average harvested energy is also increased (23).
Using Fig. 15, the relationship between the correlation coefficient and outage probability is analyzed for three different cases: i) without IM, ii) IM with N R = 2 and iii) IM with N R = 4. It is inferred that in 'without IM' the outage probability of 10 −1 is at correlation coefficient 0.92, whereas in 'IM with N R = 2' and IM with N R = 4, outage probability of 10 −1 is at correlation coefficient 0.95 and 0.97 respectively. As expected, IM improves the degrees of freedom and provides enhancement in maximum capacity of the RIS-DF-EH systems by the excess sum of log 2 (N R ).
The harvested energy at relay node R, as a function of SNR is shown in Fig 17,    2) With RIS and no compensation for phase distortion.
Only aperture gain (Nd H d V ) is obtained in this case. 3) Without the use of RIS in the system. It is observed that the harvested energy in case 1 is approximately 20dB more than the system without RIS.

VI. CONCLUSION
In this paper, the impact of spatial correlation among RIS elements in a rectangular array geometry is analyzed. A numerical approach for evaluating the product of two non-central chi-squared distribution is established as well. A novel co-operative RIS aided DF relaying with EH is proposed. The proposed system can enhance received signal power using constructive channel path utilization. Also, a substantial improvement is observed in harvested energy at the relay node. The proposed system model is extended to IM schemes for further improvement in outage performance. Both the theoretical and simulation results have perfect matching even at low SNRs and high spatial correlation. With the aforementioned characteristics, it is noted that the proposed systems provide a realistic analysis of the impact of spatial correlation on the tightly packed RIS array geometry.

APPENDIX-A
The joint PDF of The outage probability at destination node D is given by β PLi ] It describes a finite integral numerical solution for the probability of outage at destination P D out . Closed-form solution using Maclaurin approximations is not possible as it has a modified Bessel function of second kind. When γ → 0 then K g ( √ γ ) → ∞ and the exponentially decaying term inside the equation. In numerical computing environment such as MATLAB R , overloading it with twin infinite series is not recommended. Instead, an asymptotically equivalent finite sum solution for the nested infinite series in (49) is essential. The higher bound for i 1max increases the upper bound of the solution to the integral higher than unity and conversely higher bound for i 2max decreases the solution to very low values. In a numerical computation, the minimum value of upper limits of i 1 and i 2 are found by solving the following The outage probability at the destination node D is shown in Fig. 18 with d H = d V = 0.5λ c . It is observed that there is a tight match between the theoretical and simulation results.  Monte Carlo 3: for mc = 1 to 1000 do 4: for Npath = 1 to 100 do 5: θ = sin −1 (RAND) 6: φ = Phiall(RANDI(y)) 7: Generate k(θ, φ) vector 8: Compute array response a(θ, φ) vector through k(θ, φ) and pos n 9: h = multiply a(θ, φ) by i.i.d CN 10: end for 11: end for 12: return h Note: RAND is uniformly generated random variable (0 ≤ x ≤ 1). RANDI(u) is uniformly distributed random variable (0 ≤ u ≤ N c ). Properly simulating cosine PDF in (1), requires mapping an uniform random distribution (0 ≤ x ≤ 1) to cosine angular PDF, involves a simple change of variables, Simulating cos(x) with x ∈ [−π/2, π/2] would provide erroneous results.