Large Intelligent Surfaces Communicating Through Massive MIMO Rayleigh Fading Channels

Large intelligent surfaces (LIS) promises not only to improve the signal to noise ratio, and spectral efficiency but also to reduce the energy consumption during the transmission. We consider a base station equipped with an antenna array using the maximum ratio transmission (MRT), and a large reflector array sending signals to a single user. Each subchannel is affected by the Rayleigh flat fading, and the reflecting elements perform non-perfect phase correction which introduces a Von Mises distributed phase error. Based on the central limit theorem (CLT), we conclude that the overall channel has an equivalent Gamma fading whose parameters are derived from the moments of the channel fading between the antenna array and LIS, and also from the LIS to the single user. Assuming that the equivalent channel can be modeled as a Gamma distribution, we propose very accurate closed-form expressions for the bit error probability and a very tight upper bound. For the case where the LIS is not able to perform perfect phase cancellation, that is, under phase errors, it is possible to analyze the system performance considering the analytical approximations and the simulated results obtained using the well known Monte Carlo method. The analytical expressions for the parameters of the Gamma distribution are very difficult to be obtained due to the complexity of the nonlinear transformations of random variables with non-zero mean and correlated terms. Even with perfect phase cancellation, all the fading coefficients are complex due to the link between the user and the base station that is not neglected in this paper.


Introduction
The future of mobile digital communications in the age of the internet of things (IoT) requires to optimize the energy consumption for transmission, improve the signal to noise ratio (SNR) at the receiver, increase the spectral efficiency, and propose communication protocols, channel estimation methods and beamforming strategies suitable for the adopted system model.
Many solutions have been proposed as alternatives to the sixth generation of mobile communications. Zhang et al. [1] make an excellent review of the literature on these emerging techniques, citing among them large intelligent surfaces (LIS), holographic beamforming (HBF), angular orbital momentum (OAM) multiplexing, laser and visible-light communications (VLC) [2] and the advent of quantum computing which is increasingly present in large technology companies like Google and allows unmatched performance and security for quantum communication systems. Nawaz et al. [3] investigate the use of quantum machine learning strategies to improve the performance of the processes involved in the network structure since we mess with many parallel operations involving large arrays and tensors with data loaded and that through quantum computing can be mapped into large tensors product spaces where operations are handled by quantum processors that take advantage of the phenomenon of quantum superposition to achieve large communication rate and encryption security.
The proposal to make use of large intelligent surfaces (LIS) to improve transmission quality in massive MIMO systems is recent and has been gaining visibility in the literature as a concrete solution for the sixth generation (6G) of mobile communications and has presented a competitive performance in comparison with classic methods like relaying switches.
One of the great challenges for the implementation of the LIS is to estimate the channel and obtain the distribution of the fading coefficient. Wang et al. [4] propose channel estimation methods for multiuser massive MIMO systems assisted by LIS and present alternatives to decrease the training time necessary to have complete knowledge of the channel coefficients. Tataria et al. [5] discuss practical aspects of real-time implementation of LIS, especially in terms of processing and applications in radio frequency (RF) communications. Elbir et al. [6] present a deep learning framework for channel estimation, considering the massive MIMO scenario using mm-Wave.
Yu et al. [7] propose the use of LIS to improve the coverage of a cellular IoT in the so-called beyond fifth-generation (B5G). The LIS project aims to minimize the energy consumption and study the impact of channel parameters on spectral efficiency.
Ye et al. [8] propose techniques to minimize the symbol error rate (SER) by optimizing the phase shifts and the precoder for a MIMO reconfigurable intelligent surface (RIS) considering a finite alphabet of symbols. Among the strategies proposed are to fix the phase shifts and obtain the optimal precoder or to fix the precoder and find the phase shifts, that solution is useful to reduce the dimensionality of the optimization task and also the performance of the proposed RIS strategy is compared with a relay system and we see the advantage of using the techniques proposed by the authors.
Wu et al. [9] develope a mmWave point-to-point communication system assisted by multiple intelligent subsurfaces with passive reflecting elements, and antenna arrays on the transmitter and receiver. The authors derived the system achievable rate and have found the optimal precoding and power allocation for the LIS phase shift design. He et al. [10] investigate the theoretical limits and Cramér-Rao bounds for the LIS performance in a MIMO 5G system dealing with mmWave and considering the existence of a direct path (NLoS and LoS).
Dardari [11] derives analytical expressions for the channel gain and the spatial degrees-of-freedom (DoF) for the optimal LIS design considering MIMO systems. The analysis is based on electromagnetic theory and employs only geometric arguments. Jung et al. [12] consider that a MIMO system assisted by LIS can be modeled as an LoS after phase cancellation. The authors also analyze the theoretical limitations of the practical system's performance considering spatially correlated Rician channels and demonstrate that the NLoS component can be neglected when the number of antennas increases.
Yan et al. [13] present a multiuser MIMO (Mu-MIMO) system in which intelligent electromagnetic reflectors perform passive beamforming. The authors also propose to design a receiver with two estimation modules. One for the signal transmitted by the base station and the other, to estimate the additional On/Off information associated with the reflectors that modulate the digital signal arriving at them.
Badiu et al. [14] shows that the perfect estimation of the reflection angles at the LIS array is unfeasible, so we have to model the phase errors due to the estimation and discretization errors. The authors claim that the overall channel, including the LIS, can be modeled as Nakagami-m distributed, for phase errors having a generic distribution.
Cavers [15] defines maximal ratio transmission (MRT), establishing that the base station applies a vector of complex weights to compensate the downlink channel by canceling the phase and perform a signal reinforcement. He also shows a generalization for the effects of fading when the system has multiple users, although there is no exact generic solution for the optimal precoder in this scenario.
Makarfi et al. [16] propose to apply reconfigurable intelligent surfaces to expand coverage and improve the signal-to-noise ratio of a vehicular network, which can be seen as a case study of the IoT area using this new massive MIMO solution. The authors explore the idea of using smart radio environments for IoT problems and discuss some relevant aspects beyond 5G to establish communication between vehicles.
Qian et al. [17] present a MIMO system that uses LIS and has an array of antennas on the transmitter and receiver. The signals suffer uncorrelated Rayleigh fading in each channel. The authors obtain good approximations and performance studies based on analytical derivations of the statistical moments associated with the largest eigenvalues of the Wishart matrices related to the LoS and NLoS component. Without losing generality, they assume that the largest eigenvalues have a Gamma distribution and their moments are a function of the number of LIS elements and the number of antennas in the array.
Björnson et al. [18] discuss how the correlation matrix of the LIS elements can be computed under certain conditions, considering Rayleigh fading channels with a direct path between the signal and the final user. In this study, the authors make a more geometric analysis of the problem considering a rectangular panel formed by several reflectors and their constructive parameters, thereby establishing a relationship between the degrees of freedom of the LIS and the rank of the autocorrelation matrix of the reflector panel. Asymptotic analyzes of the SNR variation and channel hardening are also performed when the number of antennas and reflectors increase.
In this work, we investigate the performance of system employing LIS, also known as large reflective surfaces (LRS), taking into account Rayleigh channels and phase errors due to imperfect channel phase cancellation. This work is very general since it considers a direct link between the base station with multiple antennas and the single user. We investigate the system performance and quality of the proposed approximations for channel distribution in terms of the Kullback-Leibler divergence metric. We also present analytical expressions for the bit error probability and a very tight upper bounds for different scenarios in terms of the Von Misses parameter.
Note that the phase estimation errors is modeled as zero mean Von Mises distribution [14], which has a concentration parameter, κ, that helps us model the accuracy of the estimation. Large values of the Von Mises κ implies small errors, when κ → ∞ the zero mean Von Mises probability density function is impulsive at zero, and for κ = 0 the probability distribution is the uniform distribution.
Analytically obtaining the equivalent fading distribution and the bit error probability is quite intricate due to the considerable sums and transformations of random variables with distinct distributions.
This paper tries to cover some gaps in the literature, with regard to obtaining analytical solutions in the form of simple algebraic expressions involving the channel parameters and the Von Mises parameter, in this way we were able to simplify the analysis of a very generalist system model that may include the direct link with the user, may have one or more antennas at the base station and we maintain the validity of the analysis even when the phase correction algorithm is not efficient. Badiu et al. [14] use the central limit theorem to solve a simpler scenario with one Rayleigh NLoS channel and a rician LoS channel, but considering a single antenna transmitter. Björnson et al. [18] propose asymptotic approaches in a scenario that takes into account the correlation of the LIS elements but is also restricted to the context with a single antenna transmitter which simplifies the analytical solutions.
For organizational purposes, we summarize what is covered in each section. In Section 2, we define the mathematical notation used in the next chapters, in Section 3, the problem is well described mathematically. Section 4 discusses about the effect of phase correction errors. In Section 5, we propose an approximation for the SNR distribution and the error probability. In Section 7, we present the Monte Carlo simulations and the analytical results. Finally, in Section 8, we make our final considerations about the system performance. For easier reading of this paper, the analytical calculations of the mean and variance of the fading coefficient and the Von Mises trigonometric moments are left for the Appendix C.

Notation
The mathematical notation adopted in this paper considers E|X|, var(X), and cov(X) as the expected value, variance, and covariance of the random variable X, respectively. The function I p (κ) is the modified Bessel function of first kind and order p, and Q(.) is the Gaussian error function. The operation A • B is the Haddamard tensor product between the tensors A and B. The term X H represents the Hermitian of the complex matrix X (transposed conjugate) and z ∈ C is an element of the set of complex numbers.

System Model
In this section, we describe the mathematical model adopted in this paper and present the rationale to justify the models used for the channel, the distribution of the fading coefficients in the direct (LoS) and indirect (NLoS) links, in addition to the probability distribution associated with the error of phase accomplished by the LIS when performing the beamforming.
This paper considers a multiple-input single-output (MISO) system between a base station (BS) equipped with an antenna array composed of M antennas and a single-antenna user as shown in Figure 1. The signal path passes through the LIS environment dividing the system fading in an LoS component between the BS and the user. There are two indirect paths, between each antenna and the LIS reflector, and between each reflector and the user, these indirect links form a composite channel between the base station and the user.We suppose that the LIS is far from the BS, and the user is also far from the LIS. So, the fading coefficients are modeled as uncorrelated Rayleigh. In this formulation the signal received by the user can be written as where h LIS ∈ C N×1 is the Rayleigh channel between the LIS and the user, G = [g 1 . . . g N ] ∈ C M×N is the Rayleigh channel between the BS and the LIS, h BS ∈ C M×1 is the complex normal fading of the direct path between the antenna array and the user (LoS component), x ∈ C M×1 is the transmitted symbol after precoding and Φ = diag( e −jφ 1 . . . e −jφ N ) ∈ C N×N is a diagonal matrix representing the response of the LIS where φ n ∈ [0, 2π], ∀n is the adjustable phase-shift produced by the nth LIS's element. The variable ζ ∼ CN (0, 1) is the additive white Gaussian noise (AWGN) term. The Tx signal x is defined as x = us where u ∈ C M×1 is the precoding vector and s ∼ CN (0, 1) is the data symbol. The precoding vector u is applied by the antenna array at the BS before the transmission. Considering the MRT criterion, the optimal precoder is given as [19] where √ p is the precoder gain and w is the overall channel defined as Let η ki and θ i be the phases of g ki and h LIS i , respectively. Therefore, we can rewrite each channel fading coefficient as From (4), we see that the best situation occurs when the composite fading coefficient is perfectly corrected by the LIS and we can state that φ i = θ i + η ki . But this scenario is unfeasible because perfect channel state information is not a very realistic assumption. Therefore, both cases are approached: (i) the case where the LIS is able to perform perfect phase cancellation, and; (ii) the case where imperfect cancellation is assumed.
Considering the first case, the composite channel can be written as Let C BS,k = Re{h BS k } and = S BS,k = Im{h BS k }. Then, the square of the fading vector norm can be written as w 2 = w H w and To evaluate the system performance and understand the relationship between the bit error rate and the energy per bit applied by the transmitter, we need to know how the fading coefficients of the overall channel are distributed. Therefore, we need to obtain the statistical moments and the distribution of w 2 .

Von Mises Distributed Continuous Phase Estimation Errors
Since the phase adjustments performed by the intelligent reflectors are imperfect and cannot completely cancel the channel phase, a term associated with the phase error appears in the equation of the composite channel phase.
Consider that φ i = θ k + η ki + δ ki is the phase correction performed by the LIS, so the fading coefficients for each antenna is where the term δ ki is the phase error, here supposed as Von Mises distributed with probability density function Therefore, we have that This error model considers a matrix ∆ ∈ C M×N in which we have the Von Mises phase errors and the Haddamard product is an elementwise product between the phase errors and each channel fading magnitude. In this case, the moment generating function (MGF) of the Von Mises distribution is useful to obtain the trigonometric moments that are needed to obtain the mean and variance of the fading coefficients. For a random variable δ Von Mises distributed, the MGF can be calculated by where α p = I p (κ) I 0 (κ) and β p = 0 are defined in terms of the modified Bessel function of first kind.

Approximated Gamma Fading Distribuition
Since each fading coefficient w k is the summation of independent and equally probable random variables, we can apply the central limit theorem (CLT). So, for large values of N each w k is approximately complex Gaussian. The term w 2 is the sum of squared Gaussian random variables whose generalized distribution is the Gamma distribution. Let V ∼ Γ(α, β) be a Gamma-distributed variable with shape parameter α and rate parameter β, therefore its probability density function is given as Note that the mean of the Gamma random variable V is given by E[V] = α β , and the variance as var[V] = α β 2 [20]. We can compute the mean and variance of w 2 , denoted here as µ w 2 and σ 2 w 2 , respectively, and match with E[V] and var [V]. Using this rationale, the following can be written: Solving this linear equation system, we get that therefore, with the mean and variance of w 2 , we can generate its Gamma approximated probability density function. Although the idea might seem very simple, the mean and variance of the channel norm are very difficult to be obtained. For the sake of clarity, we detail these calculations in Appendix A, for the case where there are no phase errors. For this case, the mean µ w 2 and variance σ 2 w 2 , are given in (A17) and (A45), respectively. In the same way, for the case where phase error occurs, Appendix B presents the mean µ w 2 and variance σ 2 w 2 as in (A76) and (A91), respectively. The trigonometric moments needed to perform the calculations are in the Appendix C.

Kullback-Leibler Divergence
To evaluate the accuracy of approximating w 2 as a Gamma random variable, we can use Kullback-Leibler divergence [21]. The Gamma distribution will be compared with the simulation obtained by the Monte Carlo method.
Although the distribution of the fading coefficient is continuous, for purposes of numerical calculation, we estimate the PDF with a finite number of points and thus we also sample the Gamma distribution and calculate the Kullback-Leibler divergence in its discrete form [21] where D 1 and D 2 are the simulated and the theoretical distributions, respectively, whose probability distribution functions are d 1 (x) and d 2 (x) respectively and χ is the set of points available to represent the distributions.

Error Probability Calculations
The error probability for the M-QAM modulation can be approximately obtained by [22] P QAM where M is the size of the M-QAM constellation. Under Gamma fading, the mean error probability can be calculated asP where γ = pγ 0 and γ 0 is the SNR at the receiver whileP QAM e is the mean error probability considering the fading coefficient v and the Gamma pdf f w 2 (v).
Therefore the error probability can be expressed as where α = α w 2 and β = β w 2 are calculated by (12). In [23], we have a useful approximation for the bit error probability on an M-QAM schema. Considering coherent detection, we can state that By using the approximation in (17), we propose an upper bound for the bit error probability for the transmission of M-QAM symbols under Gamma fading by applying the Chernoff bound Q(x) < − 1 2 e − 1 2 x 2 and solving the integral formula in (16). The proposed bound for error probability can be calculated asP The gamma approximation for the resulting fading coefficient is adequate and works even for small values of N and M when we consider the scenario without phase errors, as in Figure 2, and in the case where we have the Von Mises distributed phase errors as shown in Figure 3.   As we can see in Figure 4 the Kullback-Leibler divergence decays with the number of reflectors at the LIS and the distribution is well represented by the proposed gamma approximation.

Simulated Results
We have simulated the bit error probability and generated the fading coefficients of all LoS and NLoS channels using the Monte Carlo method with 10 6 iterations. We have assumed that the phase errors follow the uniform or the Von Mises distribution. To calculate the error probability, we have solved numerically the integral in (16) and compared it with the Monte Carlo simulation. Figure 5 shows the simulated bit error rate when there are no phase errors. As it can be observed, the simulation is very close to the analytical bit error probability. Moreover, as the number of LIS reflectors increases, the bit error probability decreases faster concerning the SNR. This result is valid even for small values of N, for example, for N = 8.
On the other hand, when the phase error is uniformly distributed, as shown in Figure 6, the bit error probability increases significantly compared to the scenario where there are no phase errors. Again, we can observe that the analytical curve matches perfectly the simulated results, even for a small number of antennas or a small number of elements at the LIS.
It is worth mentioning that a uniform phase error, as pointed out by [24], may mean that the LIS's channel estimation or phase correction was not so effective since large and small phase errors are equiprobable.   Comparing Figures 6 and 7, we can see that the error probability is smaller when κ > 0. The occurrence of small errors is more probable than large errors (±π).    Figure 8 shows how the bit error probability behaves as the concentration parameter varies for SNR of −25 dB. It is clear that the rate decreases as κ increases. Therefore, the κ parameter of the LIS can be considered a qualitative parameter of the phase correction performed by the reflectors for a specific channel estimation method. In Figure 9, we vary the size of the antenna array at BS and note that the bit error rate decreases significantly when M increases. Furthermore, our approximation is valid for both large and small values of M. Although the numerical calculation of the analytical expression of the bit error probability is computationally fast, it may still be interesting to use a direct expression that does not involve solving numerical integrals. The proposed upper bound is very close to the simulated results, as we can see in Figure 10.  Figure 11 shows the influence of the direct link in the bit error probability. In this figure, we have assumed that the direct link is 10 dB and 30 dB larger than the two indirect links. As it can be seen, when the direct link is strong, the error probability decays quickly with the increase of the SNR, on the other hand, when the direct link is weak, the bit error probability requires a larger SNR to decrease.  Figure 11. Effect of the direct link in the bit error probability.

Final Considerations
This paper has presented how some parameters like number of antennas at BS and electromagnetic reflectors at LIS, channel, and phase error distribution can influence on the performance of a massive MIMO system assisted by LIS.
Assuming that there is the direct link between the user and the base station and phase errors performed by the LIS, we have derived analytical expressions for the channel probability distribution function and bit error probability. As conclusion, all the sum of the fading coefficients and phase noise involved in the MIMO communication system assisted by LIS can be modeled as a Gamma random variable.
We have proved the accuracy of our approximation through the Kullback-Leibler divergence even when the phase error follows either the uniform or the Von Mises distribution with arbitrary concentration parameter. In the absence of phase error, the divergence between the simulated distribution and the proposed analytical approach decreases even faster with the increase of the number of reflectors at the LIS.
Further studies may include expanding this analysis to more general fading distributions such as Nakagami-m, which include the Rayleigh distribution as special case and can reasonably approximate Rician fading channels for large values of m.   (5) Departing from (5), the mean of the each total fading coefficient can be calculated as The mean of the product of two random variables X and Y is given by var [20], since |g ki | and |h LIS k | are independent and equally probable on the summation variable. So, The terms |g ki | and h LIS i are Rayleigh distributed, since the variables g ki and h LIS i are zero mean with variances σ 2 1 and σ 2 2 respectively. Therefore Let w k = c k + js k , where c k and s k are the in-phase and quadrature components of the fading coefficient with respect to the antenna k. Also, let µ c k = E[c k ] and µ s k = E[s k ] = 0. Then, Appendix A.2. Variance of w k Given in (5) The variance of w k is given by and can be expanded as The first term of (A8) can be written as Note in (A9) that it is necessary to compute the variance of the product of two random variables. Let X and Y, two independent random variables, then var(XY) = var(X)var(Y) + var(X)E and so, after some simplifications, we have that Since E [C BS,k ] = 0 and the terms C BS,k and therefore the variance of w k is given as The expected value of w 2 can be calculated by With E[w k ] and var(w k ), we can calculate E |w k | 2 as E |w k | 2 = var(w k ) + (E [w k ]) 2 . Therefore, we have so the mean of the overall channel fading coefficient is

. Correlation between the Fading Coefficients
Since the real and imaginary parts of the fading coefficients, without phase errors, are uncorrelated, we analyze only the correlation between the real parts as follows.
and simplifies to Since |g il | and |g km | are independent, therefore we only need to consider the case i = k, because for i = k we have that ρ c i ,c k = ρ c k ,c k = 1 therefore by performing algebraic simplifications we have that Appendix A.5. Variance of w 2 Let Z k = |w k | 2 , so the variance of the sum of the correlated random variables Z k will be given by Since the variables have the same distribution and parameters therefore so the pairwise covariance can be obtained by and simplifies to so, Considering that c k and c i are correlated Gaussian random variables with the same mean and variance, therefore we have that where f c i ,c k (x, y) is the joint probability density function of the correlated Gaussian random variables c i and c k .
Since c i and s k are independent and E[|s 2 k |] = var(s k ) so we have that where σ 2 123 = σ 2 1 σ 2 2 σ 2 3 . The fourth order moment E[s 4 k ] of a Gaussian random variable is well known in the literature and can be calculated by With the expected value E[Z i Z k ], the expected value of Z k = |w k | 2 and the consideration that E[Z i ] = E[Z k ] ∀i∀k, the covariance between Z i and Z k can be obtained by Note that var(Z k ) = cov(Z k , Z k ), so the variance of w 2 can be obtained by By substituting the terms E[Z i Z k ] and E[Z k ], given in (A43) and (A5), in the Equation (A45), we have finally have the analytical expression of the variance in (A46).

Appendix B. Mean and Variance of the Overall Channel Fading Coefficient with Von Mises Distributed Phase Errors
Appendix B.1. Mean ofw k Letw k =c k + js k be the fading coefficient with respect to the antenna k when phase errors occurs at the LIS.
To obtain the mean ofw k we need to calculate the mean ofc k ands k . The mean value of the in-phase fading componentc k is using the linearity of the expected value we can rewrite the equation as Therefore we can calculate the mean of the overall fading coefficient, for the antenna k by using

. Variance of the In-Phase and Quadrature Components
To obtain the variance ofw k we need to calculate the variance ofc k ands k . The variance of the quadrature component is The expected value of the cosine of a Von Mises variable can be calculated by or, using a trigonometric substitution, In terms of the the characteristic function, we have var(cos δ ki ) and simplifies to The variance of the in-phase fading coefficient will be and simplifies to var(c k ) = σ 2 3 + N 2σ 2 12 (1 + α 2 ) − To compute the variance of the total fading coefficient we need the mean value of the squared fading coefficients.
The mean of the squared in-phase component is given as The expected value of the squared quadrature component can be written as Therefore, the mean squared magnitude of the overall fading with respect to the antenna k is given as and simplifies to E[|w k | 2 ] = 4N(1 + α 2 )σ 2 12 + 2σ 2 3 (A67) Since the fading coefficients are identically distributed, then and Appendix B.6. Variance of w 2 The variance of the total fading coefficient is given as the expected value of the product of the squared magnitude of two diffent fading coefficients is The term E[c 2 where k 5 = N π 2 σ 12 α 1 , k 6 = 2Nσ 2 12 (1 + α 2 ) and a 3 = π 2 σ 12 α 1 . Since E[s k ] = 0, thus the forth order moment of the quadrature component is where k 7 = 2Nσ 2 12 (1 − α 2 ). We need to obtain the term E c 2 is 2 k without considering the CLT, because the correlation coefficient betweenc ands is zero. But the termsc 2 i ands 2 k are not independent and the approach that use the integral of the product to obtain the expected value is hard to solve analytically.
Considering the definition, we have that Expanding the power and the product of the two terms we can find the expression (A84).