Secrecy Analysis and Error Probability of LIS-Aided Communication Systems under Nakagami-m Fading

Abstract

Large intelligent surfaces (LIS) are a new trend to achieve higher spectral efficiency and signal-to-noise ratio in mobile communications. For this reason, this paper proposes metrics to analyze the performance of systems with multiple antennas aided by LIS and derive the spectral efficiency, secrecy outage probability, and bit error probability in an environment with Nakagami-m distributed fading. In addition to an eavesdropper, there is a single-antenna user, an array of antennas at the transmitter side and the possibility of a direct link between transmitter and receiver. This study assumes that the LIS performs non-ideal phase cancellation leading to a residual phase error that follows a Von Mises distribution, and shows that the resulting channel can be accurately approximated by a Gamma distributed SNR whose parameters are analytically derived. From these formulas, it is possible to evaluate the effect of the strength of the line-of-sight link by varying the Nakagami parameter, m.

1. Introduction

Large Intelligent surfaces (LIS) are a promising technology for beyond fifth-generation (B5G) systems, given the number of papers emphasizing their advantages, whether compared to relays [1] or even when used to enhance the power of millimeter wave technologies [2]. Furthermore, reflecting signals with extreme precision and without power consumption can reduce the interference and improve the signal-to-noise ratio at the receiver, especially when the direct path between transmitter and destination is weak and needs to be strengthened.

In addition, known as large reflecting surfaces, they have recently been studied as a solution for different modulation schemes and communication channels. Their performance metrics show their significant potential for mobile communications. For example, Yang et al. [3] proposed a transmission protocol to reduce the channel estimation overhead when adjacent cells share the same reflection coefficients. In addition, optimization methods are used to allocate the transmit power and maximize the achievable rate in an orthogonal frequency division multiplexing (OFDM) scheme under frequency-selective channels.

In [4], Basar presented a mathematical framework to obtain the signal-to-noise ratio and derive the symbol error probability of an LIS-aided communication system, with or without knowledge of the channel phases. The author also proposed an access point sending signals directly to the users aided by a LIS system. Wymeersch et al. [2] emphasized that, although there are already other techniques for high frequencies (0.1 to 1 THz), these technologies are limited by multipath propagation and obstacles presented in the environment. In this case, LIS can control the physical propagation environment, decrease energy consumption, and simplify location and mapping systems, creating a line-of-sight (LoS) path between transmitter and receiver.

In [5], the authors presented solutions for the adjustment of the LIS elements’ phases, which optimizes the channel capacity and the precoder applied on the transmitter side. Elbir et al. [6] developed a deep learning framework to obtain the channel state information (CSI) in a massive multiuser MIMO system aided by a LIS. The authors estimated each user’s composite channel and the direct path through a convolutional neural network whose inputs are the received pilot signals. Lin et al. [7] performed channel estimation by applying Lagrange multipliers and a dual ascent-based scheme iteratively. They also found a closed-form solution for Cramer–Rao lower bounds and proposed a method that improves the accuracy of the classical least-square method. Taha et al. [8] presented an energy-efficient architecture where all the LIS’s elements are passive except for a few distributed active elements that are arranged in a non-uniform manner. The reflector array applies deep learning models to obtain the optimal matrices of phase shifts.

Although an LIS is usually a panel of reflectors physically organized in planar shapes, Hu et al. [9] proposed alternative structures with a three-dimensional spatial configuration with spherical surfaces. In addition to broader coverage, they have a more straightforward positioning system when compared to the conventional planar arrays.

LIS must be large in far-field communications to compete with classic massive MIMO systems and compensate for multipath propagation and electromagnetic interference. Besides that, optimizing the phase shifts associated with each element of the LIS is a great challenge. Therefore, in [10], Najafi et al. proposed an optimization method based on the physical modeling of the propagation and clusterization of a thousand reflectors into small subsets, also known as tiles. Based on concepts from radar communications, they modeled the impact of each tile on the overall channel, calculated the associated electric and magnetic fields, and showed that it is possible to optimize the operation of the LIS to maximize some quality of service (QoS) criteria.

On the other hand, Garcia et al. [11] focused on near-field environments and established a relation between the array size and the Fresnel zones. The punctual approximation of the scattering characterization presented dependence on the second and third-order moments of the distance. On the contrary, for far-field, the dependence is given for the fourth power. Kishk et al. [12] employed some stochastic geometry tools to analyze the effect of the large-scale deployment of LIS on the performance of cellular networks in the presence of blockages surfaces. They established a relation between the density of LIS panels and blockages.

Mukherjee [13] explores the idea of integrating LIS with mobile edge computing (MEC) technology that intends to leave computing involved in processing the received signal to a cloud server and describes how these technologies can mutually benefit and create a framework competitive for 6G. Finally, Malandrino et al. [14] analyze the possible benefits of using intelligent and reflective surfaces to increase the privacy and security of mobile communications through secrecy rate, considering that passive eavesdroppers are involved in the system, in addition to legitimate users.

In addition to the works related to optimal estimation and power control in transmission systems aided by LIS, it has become a trend to compute the channel’s capacity in the face of eavesdroppers. The question to answer is: “Does such a system offer the physical layer security that prevents an intruder from receiving a signal not intended for him?” The secrecy outage probability metric can answer this question since it means the probability that the instantaneous secrecy capacity is less than or equal to a given capacity threshold. Below are some works that, like ours, are concerned with information security in systems assisted by LIS.

Related Works

For the case of Gaussian distributed channels and considering parameters such as the distances between devices and the number of LIS elements, Yang et al. [15] derived closed-form expression for the secrecy outage probability (SOP) assuming that the LIS uses CSI to implement the phase shifting perfectly. In its turn, Trigui et al. [16] assumed a more realistic model in which there are errors caused by phase quantization. By leveraging Fox’s H transforms, they obtained exact SOP expression under the assumption that many reconfigurable elements of LIS and channels were distributed according to the Rayleigh distribution.

On the other hand, Ai et al. [17] demonstrated the potential of improving secrecy with LIS aid under different scenarios where a passive eavesdropper is attempting to retrieve the transmitted information: a vehicular-to-vehicular and a vehicular-to-infrastructure. Makarfi et al. [18] showed how the source power, eavesdropper distance, the number of LIS elements, the source-to-relay distance, and the secrecy threshold affect the secrecy capacity and SOP when the vehicular source uses an LIS as an access point.

Following the perspective of the physical layer security, this paper analyzes the secrecy outage probability of a LIS-assisted system in which K antennas at BS transmit simultaneous signals to only one user. As shown in [19], the overall fading coefficient is approximately gamma distributed, even for small values of N and K, but only when the Nakagami-m fading channels have $m=1$ (Rayleigh distribution). The reasoning is extended here to more general scenarios in which m assumes arbitrary values. To the best of the authors’ knowledge, this is the first analysis covering both channels with and without a line of sight. The derived closed-expressions for bit error probability (BER) and SOP allow us to conclude that it is possible to evaluate the system performance and design it without performing several Monte Carlo simulations that would be computationally costly in a scenario with multiple antennas and multiple reflectors. The use of the gamma approximation is investigated for a more general scenario, in which it is possible for a line of sight path to exist or not in each one of the intermediate channels (i.e., paths between the transmitter and the LIS, and between the LIS and the user).

In contrary to our previous work [19], this study focuses on secrecy analysis and extends the system model to near-field scenarios. The presence of an unwanted eavesdropper link is a realistic consideration since the information leakage becomes increasingly worrisome, especially for banking, corporate, and government communications in addition to demonstrating the validity of the proposed bit error probability approach when analyzing environments with Nakagami-m fading.

The paper is organized as follows: Section 2 presents the system model and the initial equations that based the formulation of the problem, while Section 3 presents the closed-form expressions for spectral efficiency, BER, upper bound, and SOP. Finally, Section 4 demonstrates the validity of the proposed analytical expressions through Monte Carlo simulations and Section 5 presents the final considerations. Demonstrations and mathematical deductions are presented in Appendix A.

2. System Model

As shown in Figure 1, this study considers a base station (BS) equipped with an antenna array of K antennas transmitting the same signal to a unique single antenna user, the destination. Additionally, a large intelligent surface system with N reflecting elements aids the system. Both channels BS to LIS and LIS to the user are modeled by the Nakagami-m distribution. There is a direct link between the user and the BS and between an eavesdropper and the BS whose channel is also Nakagami-m distributed. The signal that arrives at the destination antenna is given by

(1)$r=({\mathbf{h}}_{SL}^H{\mathsf{\Phi }}^H{\mathbf{h}}_{LD}^H+{\mathbf{h}}_{SD}^H)\mathsf{\Psi }+\eta ,$

where ${\mathbf{h}}_{SL}\in {\mathbb{C}}^{N\times 1}$ is the link between the source and the LIS, ${\mathbf{h}}_{LD}\in {\mathbb{C}}^{K\times N}$ is the link between the LIS and the destination and ${\mathbf{h}}_{SD}\in {\mathbb{C}}^{K\times 1}$ is the direct link between the source and the destination. The term $\mathsf{\Phi }\in {\mathbb{C}}^{N\times N}$ is a diagonal matrix, whose elements are the phase shifts $e^{-j{\varphi }_1}\dots e^{-j{\varphi }_N}$ applied by the LIS to the incident electromagnetic waves. The LIS’s phases, ${\varphi }_n\forall n$, are assumed continuous in the interval of 0 to $2\pi$ radians. The term $\mathsf{\Psi }=\mathbf{v}s$ represents the precoded signal, where the data symbol is $s\sim \mathcal{CN}(0,1)$ and the optimal precoding vector is applied by BS, according to the MRT (maximum ratio transmission) criterion, i.e.,

(2)$\mathbf{v}=\frac{{\mathbf{h}}^H}{∥\mathbf{h}∥}.$

Finally, the term $\eta \sim \mathcal{CN}(0,1)$ is additive white Gaussian noise (AWGN) with zero mean and unit variance. Suppose that there is no LoS in the direct link and that it is modeled as a complex normal random variable, with zero mean and variance ${\sigma }_{SD}^2$. Additionally, the magnitude of the channels ${\mathbf{h}}_i=\left|{\mathbf{h}}_i\right|e^{j{\varphi }_i}$ with $i\in \{SL,LD\}$ are Nakagami-m distributed with probability density function (PDF) given by

(3)$f_X(x)=\frac{2m_i^{m_i}}{\Gamma (m_i){\Omega }_i^{m_i}}{x}^{2m_i-1}{e}^{-\frac{m_i}{{\Omega }_i}{x}^2}.$

In this work, the parameters $m_i$ and ${\Omega }_i$ refer to the shape and spread of the Nakagami-m PDF, respectively. The distribution of the phases is not specified since, for this model, these phases are not relevant. Then, the overall channel, including the LIS and the antenna array, can be defined as

(4)$\mathbf{h}={\mathbf{h}}_{SL}^H{\mathsf{\Phi }}^H{\mathbf{h}}_{LD}^H+{\mathbf{h}}_{SD}^H,$

whose representation in scalar form is

(5)$h_k=\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|e^{j({\varphi }_i-{\varphi }_i^{SL}-{\varphi }_{ki}^{LD})}+{h_k}^{SD}.k\in \mathbb{N}$

Perfect phase cancelling occurs when ${\varphi }_i-{\varphi }_i^{SL}-{\varphi }_{ki}^{LD}=0$. Therefore, in this scenario, it follows that ${\varphi }_i={\varphi }_i^{SL}+{\varphi }_{ki}^{LD}$. However, the task of removing the overall channel phase is unfeasible. Some residual phase noise is left behind, in this case, ${\varphi }_i-{\varphi }_i^{SL}-{\varphi }_{ki}^{LD}={\theta }_{ki}$, where ${\theta }_{ki}$ is the phase noise, which, in this work, is modeled as a Von Mises random variable with concentration parameter $\kappa$. Therefore, the overall channel can be written as

(6)$h_k=\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|e^{j{\theta }_{ki}}+h_k^{SD}.$

It is expected that there is no phase error in the best case analysis, but this situation is entirely unfeasible. However, it is possible to estimate an optimal phase adjustment matrix that provides a performance as good as possible, so it is expected that, on average, the phase errors are zero. The zero mean Von Mises circular distribution can be proper to model the phases of each antenna’s fading coefficients [20]. It has nonzero support in the interval −$\pi$ and $\pi$ and a concentration parameter $\kappa$ associated with the quality of the phase adjustment promoted by the LIS and the efficiency of the channel estimation method.

The moment-generating function (MGF) of the Von Mises distribution is useful since a complex exponential represents the phase adjustments. With the MGF, it is possible to calculate the statistical moments associated with the channel coefficients.

Let X be a Von Mises random variable; therefore, its MGF is given by ${\phi }_p=E[e^{-jpX}]={\alpha }_p+j{\beta }_p$. Since the zero mean Von Mises distribution is symmetric about zero, then the imaginary part of the MGF ${\beta }_p=E[sinpX]=0$, and the real part is ${\alpha }_p=\frac{I_p(\kappa )}{I_0(\kappa )}$, where $I_p(\kappa )$ is the modified Bessel function of first kind and order p.

Considering that the precoder is the normalized hermitian of the overall channel, the SNR of the desired link is

(7)${\gamma }_D=({\mathbf{h}}_{SL}^H{\mathsf{\Phi }}^H{\mathbf{h}}_{LD}^H+{\mathbf{h}}_{SD}^H)\mathbf{v}={∥\mathbf{h}∥}^2.$

Assuming, as an approximation that ${\gamma }_D$ is Gamma distributed, then its statistical moments, $\alpha$ and $\beta$ can be estimated as

(8)$\alpha =\frac{\mathbb{E}{[{\gamma }_D]}^2}{\mathrm{var}({\gamma }_D)},\beta =\frac{\mathbb{E}[{\gamma }_D]}{\mathrm{var}({\gamma }_D)},$

where $\alpha$ and $\beta$ are the shape and rate parameters, while $\mathbb{E}[{\gamma }_D]$ and $\mathrm{var}({\gamma }_D)$ are the expected value and variance of ${\gamma }_D$, respectively, as shown throughout Appendix A.

The assumption that the distribution of ${\gamma }_D$ is Gamma distributed can be assessed using the Hellinger distance. According to Beran [21], the Hellinger distance between two arbitrary discrete probability distributions $p_k$ and $q_k$ can be obtained as

(9)$D_{HL}=\frac{1}{\sqrt{2}}\sqrt{\sum _{k=0}^{N_p-1}{(\sqrt{p_k}-\sqrt{q_k})}^2},$

where $N_p$ is the number of samples available to calculate the distance. The Hellinger distance is limited in the interval $0\le D_{HL}\le 1$ and can be considered as an absolute metric.

In Figure 2, the realizations of a Monte Carlo simulation of the channels involved in the system are used to compose a histogram that approximates the PDF of the overall channel that is compared to the Gamma distribution predicted by the approximation proposed in this study. To perform the analysis, ${10}^6$ iterations were performed with unit variance for all channels, Von Mises concentration parameter $\kappa =2$ and the Nakagami parameter $m=2$. The results show that the Hellinger distance decreases when N and K increase. In the last case, for $K=16$, the decrease is even more pronounced. Therefore, this accurate approximation motivates us to formulate the problem further.

3. Problem Formulation

Knowing that the SNR can be approximated by a Gamma random variable, closed-form expressions for spectral efficiency, BER and SOP are derived in the following subsections.

3.1. Spectral Efficiency

The average spectral efficiency of the system can be defined as

(10)$\begin{array}{c}C=\mathbb{E}[{log}_2(1+\gamma )]=\underset{0}{\overset{\infty }{\int }}{log}_2(1+\gamma )\frac{{\beta }^{\alpha }}{\Gamma (\alpha )}{\gamma }^{\alpha -1}{e}^{-\beta \gamma }d\gamma ,\end{array}$

whose approximated solution is given by

(11)$\begin{array}{c}C=\frac{{(-1)}^{-\alpha }\beta }{log(2)\Gamma (\alpha )}[{(-1)}^{\alpha }(\Gamma (\alpha -1){}_2{F}_2(1,1;2,2-\alpha ;\beta )+\frac{\Gamma (\alpha )({\psi }^{(0)}(\alpha )-log(\beta ))}{\beta })+\pi \beta csc(\pi \alpha )(\Gamma (\alpha )-\Gamma (\alpha ,-\beta ))],\end{array}$

where ${\psi }^{(0)}(.)$ is the digamma function, $\Gamma (.)$ is the gamma function, $\Gamma (.,.)$ is the incomplete gamma function, and ${}_2{F}_2(a,b;c,d;e)$ is the generalized hypergeometric function. It is noteworthy that, although this study did not find an explicit solution for the spectral efficiency, it does present a more generic solution for the integral in the reference [22].

3.2. Bit Error Probability

The error probability for the M-QAM modulation can be approximately obtained by [23]

(12)$P_e^{QAM}(\gamma )=1-{(1-2(1-\frac{1}{\sqrt{\mathcal{M}}})Q[\sqrt{\frac{3\gamma {log}_2\mathcal{M}}{(\mathcal{M}-1)}}])}^2.$

Assuming, as an approximation that $\gamma$ is Gamma distributed, the mean bit error probability ${\overline{P}}_e^{QAM}$ can be calculated by

(13)${\overline{P}}_e^{QAM}(\gamma )=\underset{0}{\overset{\infty }{\int }}{P}_e^{QAM}(\gamma v)f_{{∥\mathbf{w}∥}^2}(v)dv,$

where $f_{{∥\mathbf{w}∥}^2}(v)$ is the pdf of ${∥\mathbf{w}∥}^2$ as a function of an independent variable v.

Ferreira et al. [19] derived a close upper bound for the mean error probability of an M-QAM schema under Gamma fading by using the approximation

(14)${\overline{P}}_e^{QAM}(\gamma )\approx \frac{4}{{log}_{2}M}Q(\sqrt{\frac{3\gamma {log}_{2}M}{M-1}}).$

From the Chernoff bound $Q(x)\le \frac{1}{2}{e}^{-\frac{1}{2}{x}^2}$, they obtained the following upper bound for BER

(15)${\overline{P}}_e^{QAM}(\gamma )<\frac{1.38629{(\frac{2.16404\gamma log(M)}{(M-1)\beta }+1)}^{-\alpha }}{log(M)},$

which is close to the exact solution.

3.3. Secrecy Outage Probability

Considering that an eavesdropper has access to the signal provided by the source and according to [24], the secrecy capacity associated with the two fading channels can be obtained as

(16)$C=\left\{\begin{array}{cc}ln(1+{\gamma }_D)-ln(1+{\gamma }_E)\hfill & {\gamma }_D>{\gamma }_E\hfill \\ 0\hfill & {\gamma }_D\le {\gamma }_E\hfill \end{array}\right.,$

where ${\gamma }_E$ is the SNR of the link between the source and the eavesdropper. Therefore, the SOP is defined as the probability that the instantaneous secrecy capacity, C, be less than or equal to a given capacity threshold, $ln(1+{\gamma }_{th})$, which is expressed as

(17)$\mathrm{SOP}=Pr[ln\frac{(1+{\gamma }_D)}{(1+{\gamma }_E)}\le ln(1+{\gamma }_{th})]={\int }_0^{\infty }{\int }_0^{(1+{\gamma }_E)(1+{\gamma }_D)-1}{f}_{{\gamma }_E}(w)f_{{\gamma }_D}(u)dudw,$

where $Pr[.]$ denotes the probability of a random event.

Considering a Nakagami-m distributed eavesdropper channel, the SOP can be obtained as follows:

(18)$\mathrm{SOP}={\int }_0^{\infty }{\int }_0^{({\gamma }_{th}+1)(x+1)-1}\frac{{\beta }^{\alpha }(2m^m{x}^{2m-1})y^{\alpha -1}exp(-\beta y)exp(-\frac{m{x}^2}{\Omega })}{\Gamma (\alpha )({\Omega }^m\Gamma (m))}dydx$

Solving the first integral, the remaining expression becomes

(19)$\mathrm{SOP}=\underset{0}{\overset{\infty }{\int }}\frac{2e^{-\frac{m{x}^2}{\Omega }}{m}^m{x}^{2m-1}{\Omega }^{-m}[\Gamma (\alpha )-\Gamma (\alpha ,\beta (x+{\gamma }_{th}(1+x)))]}{\Gamma (m)\Gamma (\alpha )}dx,$

where the term $\Gamma (\alpha )-\Gamma (\alpha ,\beta (x+{\gamma }_{th}(1+x)))$ can be rewritten as function of the lower incomplete gamma function

(20)$\gamma (s,x)\triangleq \underset{0}{\overset{x}{\int }}{t}^{s-1}{e}^{-t}dt,$

considering that $\Gamma (s)=\gamma (s,x)+\Gamma (s,x)$.

Representing the exponential in terms of power series, it follows that the incomplete Gamma function can be written as

(21)$\underset{0}{\overset{x}{\int }}{t}^{s-1}{e}^{-t}dt=\underset{0}{\overset{x}{\int }}\sum _{k=0}^{\infty }{(-1)}^k\frac{t^{s+k-1}}{k!}dt=x^s\sum _{k=0}^{\infty }\frac{{(-x)}^k}{k!(s+k)}.$

Applying the expansion (21) in (19), and using the result obtained by the reference [22] in its table of integrals for integrands of type $x^a{e}^{-p{x}^2}\gamma (\nu ,cx)$, the SOP can be rewritten as (22), where ${}_p{\tilde{F}}_q$ is the regularized ${}_p{F}_q$ hypergeometric function and $v=\frac{1+{\gamma }_{th}}{{\gamma }_{th}}$. Although this expression is an infinite sum, it is possible to verify, in Section 4, that the error is small if only the first term is considered to compute the SOP:

(22)$\begin{array}{c}\mathrm{SOP}={\displaystyle \sum _{k=0}^{\infty }}\frac{{(-1)}^k{\beta }^{\alpha +k}{\gamma }_{th}^{\alpha +k}}{\Gamma (\alpha )\Gamma (k+1)}\\ \times (\frac{\pi m^m{2}^{-\alpha -k}{v}^{-2m}{\Omega }^{-m}\Gamma (m+\frac{1}{2})csc(\pi (\alpha +k+2m)){}_2{\tilde{F}}_2(m,m+\frac{1}{2};\frac{1}{2}(k+2m+\alpha +1),\frac{1}{2}(k+2m+\alpha +2);-\frac{m}{v^2\Omega })}{\Gamma (-k-\alpha +1)}\\ +\frac{{\pi }^{3/2}{\Omega }^{\frac{\alpha +k}{2}}{m}^{\frac{1}{2}(-\alpha -k)}{v}^{\alpha +k}}{2\Gamma (m)}(\frac{2csc(\frac{1}{2}\pi (\alpha +k+2m)){}_2{\tilde{F}}_2(\frac{1}{2}(-k-\alpha ),\frac{1}{2}(-k-\alpha +1);\frac{1}{2},\frac{1}{2}(-k-2m-\alpha +2);-\frac{m}{v^2\Omega })}{\alpha +k}\\ -\frac{\sqrt{m}sec(\frac{1}{2}\pi (\alpha +k+2m)){}_2{\tilde{F}}_2(\frac{1}{2}(-k-\alpha +1),\frac{1}{2}(-k-\alpha +2);\frac{3}{2},\frac{1}{2}(-k-2m-\alpha +3);-\frac{m}{v^2\Omega })}{v\sqrt{\Omega }})),\end{array}$

4. Numerical Results

This section analyzes the accuracy of the proposed approximations and discusses the improvements in capacity provided by LIS. In unspecified cases, this study adopts, by default, the Nakagami-m shape parameters $m_{SL}=m_{LD}=m=2$. The spread parameters ${\Omega }_{SL}={\Omega }_{LD}=\Omega$ were chosen to make the variances ${\sigma }_{SL}^2={\sigma }_{LD}^2=1$, the Von Mises concentration parameter $\kappa =2$, $K=16$ antennas at the source, the size of the M-QAM constellation is $M=16$, and the number of iterations is ${10}^6$ for each Monte Carlo simulation.

For each iteration of the Monte Carlo method, we generate the coefficients $h_k$ of (6), the magnitudes are generated using the Nakagami-m distribution and the phase errors with the Von Mises distribution. Given the coefficients, it is easy to estimate the bit error rate, spectral efficiency, and the SOP in each realization of the random variables and approximating the simulated results by the mean value of these quantities. We compare each of the simulated results in several iterations with the theoretical formulas described in terms of channel parameters.

Figure 3 shows the simulated and theoretical BERs considering the Von Mises and uniformly distributed phase errors. The theoretical BER is obtained assuming that the overall fading channel has a Gamma distribution. Note that the larger the number of reflectors, N, the smaller the error probability for any SNR value. When the phase errors are uniformly distributed ($\kappa =0$), the error probability is higher than in the Von Mises scenario. This result shows the importance of accurately estimating the phases and channel gains and choosing the optimization method to find the best LIS phase shifts. Uniformly distributed phase noise indicates that the algorithm has equal chances to present significant phase errors (close to $\pm \pi$) or small phase errors (close to zero). That implies greater bit error probabilities, which can be compensated only with a large number of antennas at the transmitter or with a large number of reflectors at the LIS.

In its turn, Figure 4 confirms that large reflecting surfaces can produce an LoS link between the transmitter and the user even in a far-field Rayleigh fading channel. However, in a near-field scenario, a stronger LoS link (higher Nakagami-$m$ parameter) implies a lower probability of error.

Even in weak LoS scenarios, LIS can decrease the probability of bit error by creating an LoS that is the result of beamforming toward the target user. Figure 3 shows that, for an environment with a fixed value of m, the increase in the number of reflectors (N), or the improvement of the phase adjustment performed by the LIS (related to the concentration parameter $\kappa$) can reduce the bit error rate in an aided LIS system.

The upper bound (15) for the error probability proposed by Ferreira et al. [19] is very close to the bit error rate as shown in Figure 5, even when the fading coefficients are Nakagami-m distributed.

Regarding spectral efficiency, this study considers two scenarios. The first one has uniformly distributed phase noise, and the second one has a Von Mises distributed phase, as shown in Figure 6. Notably, the spectral efficiency increases when the LIS has a more significant number of reflectors, thus indicating a better sharing of the spectrum for the transmission of signals for a multiuser scenario. Moreover, the efficiency is higher for the case in which the phase errors have a Von Mises distribution and lower when the phase errors are uniformly distributed, which means that the phase adjustment of the LIS is a highly relevant factor in improving the spectral efficiency, reinforcing the importance of channel estimation and choice of the phase correction applied to reflectors. When the phase error distribution is more concentrated around zero (higher $\kappa$ values), then the spectral efficiency is higher. Using an array of antennas on the base station can be a good choice to achieve better spectrum sharing in diverse scenarios. It is also remarkable that the result predicted by the formula proposed for the spectral efficiency is very close to the results obtained by the Monte Carlo simulation.

In Figure 7, the secrecy outage probability for a Nakagami-$m$ eavesdropper link with $\Omega =1$ and $m=1.4$ is shown. The sum was truncated up to the index 1000, and the number of iterations used was ${10}^6$ to generate the Gamma distributed random SNR with parameters $\alpha$ and $\beta$, the Von Mises concentration parameter $\kappa =2$, $K=2$ antennas, unity variance, and Nakagami-m fading distribution for all channels between the antennas, the LIS, and the user. The larger the number of reflectors or the SNR, then the greater is the SOP.

The first-order approximation of the SOP, considering that the Nakagami-m parameters are $m=2.5$ and $\Omega =0.1$ for all the channels in the system model, is also close to the simulated result as shown in Figure 8.

5. Final Considerations

This work has presented an in-depth analysis of the performance of systems aided by large intelligent surfaces considering the existence of an eavesdropper link in generic scenarios that contemplate channels with and without LoS links, employing the Nakagami-m distribution, and channels with or without a direct link to the transmitter and the user. This study derives very accurate analytical expressions from computing the secrecy outage probability, bit error probability, and secrecy capacity, in addition to reasonable approximations for estimating the equivalent channel parameters based on the central limit theorem.

Author Contributions

All authors contributed equally to this study working in the bibliographic review, calculating the mathematical expressions, executing the simulations, and analysis. All authors have read and agreed to the published version of the manuscript.

Funding

This work was partially supported by RNP, with resources from MCTIC, Grant No. 01245.010604/2020-14, under the 6G Mobile Communications Systems project of the Radiocommunication Reference Center (Centro de Referência em Radiocomunicações—CRR) of the National Institute of Telecommunications (Instituto Nacional de Telecomunicações—Inatel), Brazil, R&D ANEEL—PROJECT COPEL 2866-0366/2013 and EMBRAPII—Empresa Brasileira de Pesquisa e Inovação Industrial, CNPq under Grant 301777/2019-5 (Gustavo Fraidenraich) and 313239/2017-7 (Eduardo R. de Lima).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Acknowledgments

We are thankful for the financial support of the Eldorado Research Institute that paid the charges involved in the publication of this study.

Conflicts of Interest

The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. System model with eavesdropper link.

Figure 2. Hellinger distance.

Figure 3. Bit error probability varying the SNR.

Figure 4. Bit error probability varying the LoS strength.

Figure 5. Bit error probability upper bound.

Figure 6. Spectral Efficiency varying M and N.

Figure 7. SOP for [Forumla omitted. See PDF.] antennas and different number of reflectors (N).

Figure 8. First order approximation of SOP for [Forumla omitted. See PDF.] antennas and different number of reflectors (N).

Appendix A. Parameters of γ_D

It is possible to obtain the parameters of the $${\gamma }_D$$ distribution, following the steps below.

Appendix A.1. Expected Value of Each Fading Coefficient

Since the expected value is a linear operator, then (A1)$$\mathbb{E}[h_k]=\mathbb{E}[\sum _{i=1}^N|h_{ki}^{LD}|\left|h_i^{SL}\right|e^{j{\theta }_{ki}}+h_k^{SD}]=N\mathbb{E}[\left|h_{ki}^{LD}{h}_i^{SL}\right|e^{j{\theta }_{ki}}]+\mathbb{E}[h_k^{SD}].$$

Since the channels are independent and identically distributed, $$\mathbb{E}[h_k^{SD}]=0$$, and $$\mathbb{E}[e^{j{\theta }_{ki}}]={\alpha }_1$$, thus (A2)$$\mathbb{E}[h_k]=N{\mu }_{LD}\times {\mu }_{SL}\times {\alpha }_1,$$ where (A3)$${\mu }_{LD}=\mathbb{E}[\left|h_{ki}^{LD}\right|]=\frac{\Gamma (m_{LD}+0.5)}{\Gamma (m_{LD})}{(\frac{{\Omega }_{LD}}{m_{LD}})}^{\frac{1}{2}},$$ and (A4)$${\mu }_{SL}=\mathbb{E}[\left|h_i^{SL}\right|]=\frac{\Gamma (m_{SL}+0.5)}{\Gamma (m_{SL})}{(\frac{{\Omega }_{SL}}{m_{SL}})}^{\frac{1}{2}},$$ are the expected values of each Nakagami-m channel. To obtain the variance of the overall channel fading coefficient, the mean of $$c_k=Re\left\{h_k\right\}$$ and $$s_k=Im\left\{h_k\right\}$$ needs to be calculated, i.e., the in-phase and quadrature components of the fading coefficient, respectively. The in-phase component can be written as (A5)$$c_k=\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki}+Re\left\{h_k^{SD}\right\},$$ while the quadrature component is (A6)$$s_k=\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki}+Im\left\{h_k^{SD}\right\}.$$

Then, the expected value of $$c_k$$ is (A7)$$\mathbb{E}[c_k]=\sum _{i=1}^N\mathbb{E}[\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki}]+\mathbb{E}[Re\left\{h_k^{SD}\right\}].$$

Since $$\mathbb{E}[Re\left\{h_k^{SD}\right\}]=0$$ and all the summation terms are independent, (A8)$$\mathbb{E}[c_k]=N\times \mathbb{E}[\left|h_{ki}^{LD}\right|]\mathbb{E}[\left|h_i^{SL}\right|]\mathbb{E}[cos{\theta }_{ki}]N\times {\mu }_{LD}{\mu }_{SL}{\alpha }_1,$$ as well as $$\mathbb{E}[h_k]$$.

In its turn, the expected value of $$s_k$$ is given by (A9)$$\mathbb{E}[s_k]=\sum _{i=1}^N\mathbb{E}[\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki}]+\mathbb{E}[Im\left\{h_k^{SD}\right\}].$$

Since $$\mathbb{E}[Im\left\{h_k^{SD}\right\}]=0$$ and all the summation terms are independent, (A10)$$\mathbb{E}[s_k]=N\times \mathbb{E}[\left|h_{ki}^{LD}\right|]\mathbb{E}[\left|h_i^{SL}\right|]\mathbb{E}[sin{\theta }_{ki}]=0,$$ since $$\mathbb{E}[sin{\theta }_{ki}]=0$$.

Appendix A.2. Variance of the In-Phase and Quadrature Components of Each Fading Coefficient

The variance of the in-phase component is written as (A11)$$\mathrm{var}(c_k)=\mathrm{var}(\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki}+Re\left\{h_k^{SD}\right\}).$$

Since $$\mathrm{var}(h_k^{SD})={\sigma }_{SD}^2$$, the summation terms and $$Re\left\{h_k^{SD}\right\}$$ are independent and the $$h_k^{SD}$$ coefficient is zero mean, then (A12)$$\mathrm{var}(c_k)=N\times \mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki})+\frac{{\sigma }_{SD}^2}{2}.$$

Next, the variance of the term $$\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki}$$ is needed, considering that the variance of the product of two random variables X and Y is $$\mathrm{var}(XY)=\mathrm{var}(X)\mathrm{var}(Y)+\mathrm{var}(X)E{[Y]}^2+\mathrm{var}(Y)E{[X]}^2$$ and the phase noise is independent of the fading magnitudes. Thus, (A13)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki})=\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|)\mathrm{var}(cos{\theta }_{ki})+\dots \mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|){(\mathbb{E}[cos{\theta }_{ki}])}^2+\mathrm{var}(cos{\theta }_{ki}){(\mathbb{E}[\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|])}^2.$$

Since $$\left|h_{ki}^{LD}\right|$$ and $$\left|h_i^{SL}\right|$$ are independents, (A14)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|)=\mathrm{var}(\left|h_{ki}^{LD}\right|)\mathrm{var}(\left|h_i^{SL}\right|)+\mathrm{var}(\left|h_{ki}^{LD}\right|){(\mathbb{E}[\left|h_i^{SL}\right|])}^2+\dots \mathrm{var}(\left|h_i^{SL}\right|){(\mathbb{E}[\left|h_{ki}^{LD}\right|])}^2$$ and (A15)$$\mathbb{E}[\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|]={\mu }_{LD}\times {\mu }_{SL}.$$

Next, considering that $$\mathrm{var}(\left|h_i^{SL}\right|)={\sigma }_{SL}^2$$ and $$\mathrm{var}(\left|h_{ki}^{LD}\right|)={\sigma }_{LD}^2$$, then (A16)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|)={\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2.$$

By using (A16), (A13) can thus be evaluated as (A17)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki})=({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)\times \frac{1}{2}(1+{\alpha }_2)+({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2){\alpha }_1^2+\frac{1}{2}(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2,$$ which can be rewritten as (A18)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|cos{\theta }_{ki})=\frac{1}{2}(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+\dots \frac{1}{2}({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2).$$

Therefore, the variance of the in-phase component is (A19)$$\mathrm{var}(c_k)=\frac{{\sigma }_{SD}^2}{2}+\frac{N}{2}\times [(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+\dots ({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)].$$

On other hand, the variance of the quadrature component of each fading coefficient is given by (A20)$$\mathrm{var}(s_k)=\mathrm{var}(\sum _{i=1}^N\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki}+Im\left\{h_k^{SD}\right\}).$$

Considering that $$\mathrm{var}(h_k^{SD})={\sigma }_{SD}^2$$, the summation terms and $$Im\left\{h_k^{SD}\right\}$$ are independent and the term $$h_k^{SD}$$ is zero mean, then (A21)$$\mathrm{var}(s_k)=N\times \mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki})+\frac{{\sigma }_{SD}^2}{2}.$$

Next, for that reason, (A22)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki})=\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|)\mathrm{var}(sin{\theta }_{ki})+\dots \mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|){(\mathbb{E}[sin{\theta }_{ki}])}^2+\mathrm{var}(sin{\theta }_{ki}){(\mathbb{E}[\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|])}^2.$$

Since $$\mathbb{E}[sin{\theta }_{ki}]=0$$, applying (A16) therefore (A23)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki})=\frac{1}{2}(1-{\alpha }_2){\mu }_{SL}^2{\mu }_{LD}^2+\dots \frac{1}{2}(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2),$$ which can be rewritten as (A24)$$\mathrm{var}(\left|h_{ki}^{LD}\right|\left|h_i^{SL}\right|sin{\theta }_{ki})=\frac{1}{2}(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2).$$

Therefore, the variance of the quadrature component is expressed by (A25)$$\mathrm{var}(s_k)=\frac{{\sigma }_{SD}^2}{2}+\frac{N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2)}{2}.$$

Appendix A.3. Expected Value of γ_D

All the fading coefficients are independent and identically distributed and $$\mathbb{E}[{\gamma }_D]=\mathbb{E}[{∥\mathbf{h}∥}^2]$$. Therefore, (A26)$$\mathbb{E}[{\gamma }_D]=\mathbb{E}[\sum _{k=1}^K{\left|h_k\right|}^2]=K\times \mathbb{E}[{\left|h_k\right|}^2].$$

Considering that $$\mathbb{E}[{\left|h_k\right|}^2]=\mathbb{E}[c_k^2+s_k^2]=\mathbb{E}[c_k^2]+\mathbb{E}[s_k^2]$$, then $$\mathbb{E}[c_k^2]=\mathrm{var}(c_k)+{(\mathbb{E}[c_k])}^2,$$ and, by using (A19) and (A25), (A27)$$\mathbb{E}[c_k^2]=\frac{{\sigma }_{SD}^2}{2}+N^2\times {\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+\frac{N}{2}\times [(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+\dots ({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)].$$

For the quadrature component, consider that $$\mathbb{E}[s_k]=0$$, and, consequently, $$\mathbb{E}[s_k^2]=\mathrm{var}(s_k)$$. Then, (A28)$$\mathbb{E}[s_k^2]=\frac{{\sigma }_{SD}^2}{2}+\frac{N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2)}{2}.$$ and the mean value of the overall fading coefficient magnitude is given by (A29)$$\begin{array}{c}\mathbb{E}[{\gamma }_D]=\frac{K}{2}\times (2{\sigma }_{SD}^2+2N^2\times {\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^{2}N[(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+\dots \\ ({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)]+\dots \\ N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2)).\end{array}$$

Appendix A.4. Variance of γ_D

Let $$Z_i={\left|h_i\right|}^2$$; then, (A30)$$\mathrm{var}({\gamma }_D)=\mathrm{var}({∥\mathbf{h}∥}^2)=\mathrm{var}(\sum _{i=1}^M{Z}_i).$$

Therefore, the variance of the sum of the terms $$Z_i$$ is given by (A31)$$\mathrm{var}(\sum _{i=1}^M{Z}_i)=\sum _{i=1}^M\mathrm{var}(Z_i)+2\sum _{1\le i<k\le M}\mathrm{cov}(Z_i,Z_k),$$ whose magnitudes are equally distributed. Therefore, (A32)$$\mathrm{var}(\sum _{i=1}^M{Z}_i)=M\mathrm{var}(Z_i)+M(M-1)\mathrm{cov}(Z_i,Z_k).$$

The covariance can be obtained by (A33)$$\mathrm{cov}(Z_i,Z_k)=\mathbb{E}[Z_i{Z}_k]-\mathbb{E}[Z_i]\mathbb{E}[Z_k],$$ where (A34)$$\mathbb{E}[Z_i{Z}_k]=\mathbb{E}[c_i^2{c}_k^2]+2\mathbb{E}[c_i^2{s}_k^2]+\mathbb{E}[s_i^2{s}_k^2].$$

The expected value of the product of two different in-phase coefficients can be written as (A35)$$\mathbb{E}[c_i{c}_k]=\mathbb{E}[(\sum _{l=1}^N\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il}+Re\left\{h_i^{SD}\right\})\times (\sum _{m=1}^N\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|cos{\theta }_{km}+Re\left\{h_k^{SD}\right\})].$$

By expanding the product, (A36)$$\mathbb{E}[c_i{c}_k]=\mathbb{E}[Re\left\{h_i^{SD}\right\}Re\left\{h_k^{SD}\right\}]+\sum _{l=1}^N\sum _{m=1}^N\mathbb{E}[\left|h_{il}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{km}],$$ where the independent terms can be separated as (A37)$$\mathbb{E}[\left|h_{il}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{km}]=\mathbb{E}{[\left|h_{km}^{LD}\right|]}^2\mathbb{E}{[\left|h_m^{SL}\right|]}^2\mathbb{E}{[cos{\theta }_{km}]}^2={\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2\forall l\ne m,$$ and (A38)$$\mathbb{E}[\left|h_{il}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{km}]=\mathbb{E}{[\left|h_{km}^{LD}\right|]}^2\mathbb{E}[{\left|h_m^{SL}\right|}^2]\mathbb{E}[cos{\theta }_{im}]\mathbb{E}[cos{\theta }_{km}]={\mu }_{LD}^2({\sigma }_{SL}^2+{\mu }_{SL}^2){\alpha }_1^2\forall i\ne k,l=m,$$ where the term $$\mathbb{E}[{\left|h_m^{SL}\right|}^2]={\sigma }_{SL}^2+{\mu }_{SL}^2$$ and the variance of the Nakagami-m distributed term is (A39)$${\sigma }_{SL}^2={\Omega }_{SL}(1-\frac{1}{m_{SL}}{(\frac{\Gamma (m_{SL}+0.5)}{\Gamma (m_{SL})})}^2).$$

For $$i=k$$, (A40)$$\mathbb{E}[\left|h_{il}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{km}]=\mathbb{E}[{\left|h_{km}^{LD}\right|}^2]\mathbb{E}[{\left|h_m^{SL}\right|}^2]\mathbb{E}[{cos}^2{\theta }_{km}]=\frac{1}{2}({\sigma }_{LD}^2+{\mu }_{LD}^2)({\sigma }_{SL}^2+{\mu }_{SL}^2)(1+{\alpha }_2)\forall l=m,i=k,$$ where (A41)$${\sigma }_{LD}^2={\Omega }_{LD}(1-\frac{1}{m_{LD}}{(\frac{\Gamma (m_{LD}+0.5)}{\Gamma (m_{LD})})}^2).$$

To obtain the variance according to (A34), the term $$\mathbb{E}[c_i^2{c}_k^2]$$ needs to be computed. The terms are approximately correlated Gaussian random variables by the central limit theorem (CLT), for large values of N, therefore (A42)$$\mathbb{E}[c_i^2{c}_k^2]=\underset{-\infty }{\overset{\infty }{\int }}\underset{-\infty }{\overset{\infty }{\int }}{x}^2{y}^2{f}_{c_i,c_k}(x,y)dxdy,$$ where $$f_{c_i,c_k}(x,y)$$ is the joint distribution of the two correlated Gaussian variables $$c_i$$ and $$c_k$$. The result of the integral is (A43)$$\mathbb{E}[c_i^2{c}_k^2]={\mu }_{c_k}^4+2{\mu }_{c_k}^2(1+2{\rho }_{c_i,c_k}){\sigma }_{c_k}^2+(1+2{\rho }_{c_i,c_k}^2){\sigma }_{c_k}^4,$$ where $${\mu }_{c_k}=\mathbb{E}[c_k]$$, $${\sigma }_{c_k}^2=\mathrm{var}(c_k)$$, and since $$\mathrm{var}(c_i)=\mathrm{var}(c_k)$$ and $$\mathbb{E}[c_i]=\mathbb{E}[c_k]$$, the correlation coefficient $${\rho }_{c_i,c_k}$$ can be obtained as (A44)$${\rho }_{c_i,c_k}=\frac{\mathbb{E}[c_i{c}_k]-{\mu }_{c_k}^2}{\mathrm{var}(c_k)},$$ where $$\mathbb{E}[c_i{c}_k]$$ can be calculated by (A45) in the next page. (A45)$$\mathbb{E}[c_i{c}_k]=\left\{\begin{array}{c}N(N-1)({\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2)+N({\mu }_{LD}^2({\sigma }_{SL}^2+{\mu }_{SL}^2){\alpha }_1^2)i\ne k\hfill \\ N(N-1)({\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2)+\frac{N}{2}(({\sigma }_{LD}^2+{\mu }_{LD}^2)({\sigma }_{SL}^2+{\mu }_{SL}^2)(1+{\alpha }_2))+\frac{1}{2}{\sigma }_{SD}^{2}i=k\hfill \end{array}\right.$$

Since the moments of $$c_k$$ were previously calculated, $${\rho }_{c_i,c_k}$$ can be obtained by (A46)$${\rho }_{c_i,c_k}=\frac{N(N-1)({\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2)+N({\mu }_{LD}^2({\sigma }_{SL}^2+{\mu }_{SL}^2){\alpha }_1^2)-{(N{\mu }_{LD}{\mu }_{SL}{\alpha }_1)}^2}{\frac{{\sigma }_{SD}^2}{2}+\frac{N}{2}[(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)]}i\ne k.$$

In its turn, applying (A46) for the correlation coefficient $${\rho }_{c_i,c_k}$$ in the definition of $$\mathbb{E}[c_i^2{c}_k^2]$$ in (A43), it is possible to find that (A47). (A47)$$\mathbb{E}[c_i^2{c}_k^2]=\left\{\begin{array}{ll}{a}_1^4(1-\frac{4}{N})+a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3-\frac{2}{N}+4a_5)2a_5+\frac{1}{4}{({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)}^2& i\ne k\\ a_1^4+3a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)+\frac{3}{4}({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)& i=k\end{array}\right.$$ where $$a_1=N{\mu }_{SL}{\mu }_{LD}{\alpha }_1,a_2={\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2,a_3=1+{\alpha }_2+2{\alpha }_1^2,a_4=N{\mu }_{SL}^2{\mu }_{LD}^2(1+{\alpha }_2),$$$$a_5=N{\mu }_{LD}^2{\alpha }_1^2({\sigma }_{SL}^2+{\mu }_{SL}^2)$$. The term $$\mathbb{E}[c_i^2{s}_k^2]$$ depends on two uncorrelated but not independent random variables $$c_i$$ and $$s_k$$, the correlation is zero and, to calculate the correlation of the squared product, it is necessary to expand the very definition of the two terms as shown in (A48)$$\mathbb{E}[c_i^2{s}_k^2]=\mathbb{E}[{(\sum _{l=1}^N\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il}+Re\left\{h_i^{SD}\right\})}^2\times {(\sum _{m=1}^N\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|sin{\theta }_{km}+Im\left\{h_k^{SD}\right\})}^2]$$

Expanding the product in (A48), it follows that (A49)$$\begin{array}{c}\mathbb{E}[{({\displaystyle \sum _{l=1}^N}\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il})}^2{({\displaystyle \sum _{m=1}^N}\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|sin{\theta }_{km})}^2]=\\ \mathbb{E}[{\displaystyle \sum _{l=1}^N}{\displaystyle \sum _{t=1}^N}{\displaystyle \sum _{d=1}^N}{\displaystyle \sum _{m=1}^N}\left|h_{il}^{LD}\right|\left|h_{it}^{LD}\right|\left|h_{kd}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|\left|h_d^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{it}sin{\theta }_{kd}sin{\theta }_{km}]\end{array}$$ where the general and simplified result given in Ferreira et al. [19] as (A50)$$\begin{array}{c}\mathbb{E}[c_i^2{s}_k^2]=\mathbb{E}[{({\displaystyle \sum _{l=1}^N}\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il})}^2{({\displaystyle \sum _{m=1}^N}\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|sin{\theta }_{km})}^2]\dots \\ +\mathbb{E}[{({\displaystyle \sum _{l=1}^N}\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il})}^2]\frac{{\sigma }_{SD}^2}{2}+\mathbb{E}[{({\displaystyle \sum _{m=1}^N}\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|sin{\theta }_{km})}^2]\frac{{\sigma }_{SD}^2}{2}+\frac{{\sigma }_{SD}^4}{4}\end{array}$$

The in-phase term is (A51)$$\mathbb{E}[{(\sum _{l=1}^N\left|h_{il}^{LD}\right|\left|h_l^{SL}\right|cos{\theta }_{il})}^2]=\mathbb{E}[\sum _{l=1}^N\sum _{t=1}^N\left|h_{il}^{LD}\right|\left|h_{it}^{LD}\right|\left|h_l^{SL}\right|\left|h_m^{SL}\right|cos{\theta }_{il}cos{\theta }_{it}]=\frac{N}{2}{\chi }_{LD}{\chi }_{SL}(1+{\alpha }_2)+N(N-1){\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2$$ and the quadrature term is (A52)$$\mathbb{E}[{(\sum _{m=1}^N\left|h_{km}^{LD}\right|\left|h_m^{SL}\right|sin{\theta }_{km})}^2]=\mathbb{E}[\sum _{d=1}^N\sum _{m=1}^N\left|h_{kd}^{LD}\right|\left|h_{km}^{LD}\right|\left|h_d^{SL}\right|\left|h_m^{SL}\right|sin{\theta }_{kd}sin{\theta }_{km}]=N({\chi }_{LD}{\chi }_{SL}\frac{1}{2}(1-{\alpha }_2))$$ where (A53)$${\chi }_{LD}=\mathbb{E}[{\left|h_{il}^{LD}\right|}^2]={\Omega }_{LD}$$ (A54)$${\chi }_{SL}=\mathbb{E}[{\left|h_d^{SL}\right|}^2]={\Omega }_{SL}$$ (A55)$${\xi }_{LD}=\mathbb{E}[{\left|h_{il}^{LD}\right|}^3]=\frac{\Gamma (1.5+m_{LD})}{(\frac{m_{LD}}{{\Omega }_{LD}})\Gamma (m_{LD})}$$ (A56)$${\xi }_{SL}=\mathbb{E}[{\left|h_d^{SL}\right|}^3]=\frac{\Gamma (1.5+m_{SL})}{(\frac{m_{SL}}{{\Omega }_{SL}})\Gamma (m_{SL})}$$ (A57)$${\tau }_{LD}=\mathbb{E}[{\left|h_{il}^{LD}\right|}^4]=\frac{(1+m_{LD}){\Omega }_{LD}^2}{m_{LD}}$$ (A58)$${\tau }_{SL}=\mathbb{E}[{\left|h_d^{SL}\right|}^4]=\frac{(1+m_{SL}){\Omega }_{SL}^2}{m_{SL}}$$

The term $$\mathbb{E}[s_i^2{s}_k^2]$$ can be easily obtained because $$s_i$$ and $$s_k$$ are uncorrelated and zero mean, then (A59)$$\mathbb{E}[s_i^2{s}_k^2]=\left\{\begin{array}{ll}{(\mathbb{E}[s_i^2])}^2& i\ne k\\ \mathbb{E}[s_k^4]=3{(\mathrm{var}(s_k))}^2& i=k\end{array}\right.$$ and, by substituting the variances; consequently, (A60). (A60)$$\mathbb{E}[s_i^2{s}_k^2]=\left\{\begin{array}{cc}\frac{1}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2\hfill & i\ne k\hfill \\ \\ \frac{3}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2\hfill & i=k\hfill \end{array}\right.$$

Finally, the term $$\mathbb{E}[Z_i{Z}_k]$$ can be calculated by using the formulas of $$\mathbb{E}[c_i^2{c}_k^2]$$, $$\mathbb{E}[c_i^2{s}_k^2]$$ and $$\mathbb{E}[s_i^2{s}_k^2]$$, resulting in (A61). (A61)$$\begin{array}{c}\mathbb{E}[Z_i{Z}_k]=\left\{\begin{array}{l}{a}_1^4(1-\frac{4}{N})+a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3-\frac{2}{N}+4a_5)2a_5+\frac{1}{4}{({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)}^2+\\ \frac{1}{2}(N^2-N)({\chi }_{LD}^2{\chi }_{SL}^2(1-{\alpha }_2^2)+4{\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}{\xi }_{SL}{\alpha }_1^2(1-{\alpha }_2))+\frac{N}{2}{\chi }_{LD}^2{\tau }_{SL}(1-{\alpha }_2^2)\\ +\frac{1}{2}{\sigma }_{SD}^2(N{\chi }_{LD}{\chi }_{SL}(1+{\alpha }_2)+2(N-1){\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+{\chi }_{LD}{\chi }_{SL}(1-{\alpha }_2))+\\ (N^3-3N^2+2N){\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}^2{\chi }_{SL}{\alpha }_1^2(1-{\alpha }_2)+\frac{1}{2}{\sigma }_{SD}^4+\\ \frac{1}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2& i\ne k\\ a_1^4+3a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)+\frac{3}{4}({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)+\\ \frac{1}{4}(N^2-N)({\chi }_{LD}^2{\chi }_{SL}^2(1-{\alpha }_4)+4{\mu }_{LD}{\xi }_{LD}{\mu }_{SL}{\xi }_{SL}({\alpha }_1-{\alpha }_3){\alpha }_1)+\frac{N}{2}{\tau }_{LD}{\tau }_{SL}(1-{\alpha }_2^2)+\\ (N^3-3N^2+2N){\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}^2{\chi }_{SL}{\alpha }_1^2(1-{\alpha }_2)+\frac{N}{2}{\sigma }_{SD}^2[{\chi }_{LD}{\chi }_{SL}(1+{\alpha }_2)+\\ 2(N-1){\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+{\chi }_{LD}{\chi }_{SL}(1-{\alpha }_2)]+\\ \frac{1}{2}{\sigma }_{SD}^4+\frac{3}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2& i=k\end{array}\right.\end{array}$$

The covariance $$\mathrm{cov}(Z_i{Z}_k)$$ can be calculated by (A62) using the result of (A61). Since $$\mathrm{var}(Z_i)=\mathrm{cov}(Z_i,Z_i)$$, the variance of the overall fading coefficient can be easily obtained by (A32) as shown in (A63): (A62)$$\begin{array}{c}\mathrm{cov}(Z_i{Z}_k)=\left\{\begin{array}{c}{a}_1^4(1-\frac{4}{N})+a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3-\frac{2}{N}+4a_5)2a_5+\frac{1}{4}{({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)}^2+\hfill \\ \frac{1}{2}(N^2-N)({\chi }_{LD}^2{\chi }_{SL}^2(1-{\alpha }_2^2)+4{\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}{\xi }_{SL}{\alpha }_1^2(1-{\alpha }_2))+\frac{N}{2}{\chi }_{LD}^2{\tau }_{SL}(1-{\alpha }_2^2)\hfill \\ +\frac{1}{2}{\sigma }_{SD}^2(N{\chi }_{LD}{\chi }_{SL}(1+{\alpha }_2)+2(N-1){\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+{\chi }_{LD}{\chi }_{SL}(1-{\alpha }_2))+\hfill \\ (N^3-3N^2+2N){\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}^2{\chi }_{SL}{\alpha }_1^2(1-{\alpha }_2)+\frac{1}{2}{\sigma }_{SD}^4+\hfill \\ \frac{1}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2\hfill \\ -\frac{1}{4}\times (2{\sigma }_{SD}^2+2N^2\times {\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+\hfill \\ N[(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)]+\hfill \\ N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2){)}^2\hfill & i\ne k\hfill \\ a_1^4+3a_1^2({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)+\frac{3}{4}({\sigma }_{SD}^2+a_4+N{a}_2{a}_3)+\hfill \\ \frac{1}{4}(N^2-N)({\chi }_{LD}^2{\chi }_{SL}^2(1-{\alpha }_4)+4{\mu }_{LD}{\xi }_{LD}{\mu }_{SL}{\xi }_{SL}({\alpha }_1-{\alpha }_3){\alpha }_1)+\frac{N}{2}{\tau }_{LD}{\tau }_{SL}(1-{\alpha }_2^2)+\hfill \\ (N^3-3N^2+2N){\mu }_{LD}^2{\chi }_{LD}{\mu }_{SL}^2{\chi }_{SL}{\alpha }_1^2(1-{\alpha }_2)+\frac{N}{2}{\sigma }_{SD}^2[{\chi }_{LD}{\chi }_{SL}(1+{\alpha }_2)+\hfill \\ 2(N-1){\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+{\chi }_{LD}{\chi }_{SL}(1-{\alpha }_2)]+\frac{1}{2}{\sigma }_{SD}^4+\hfill \\ \frac{3}{4}{({\sigma }_{SD}^2+N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2))}^2\hfill \\ -\frac{1}{4}\times (2{\sigma }_{SD}^2+2N^2\times {\mu }_{LD}^2{\mu }_{SL}^2{\alpha }_1^2+\dots \hfill \\ N[(1+{\alpha }_2)\times {\mu }_{SL}^2{\mu }_{LD}^2+\dots ({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2)(1+{\alpha }_2+2{\alpha }_1^2)]+\hfill \\ N(1-{\alpha }_2)({\sigma }_{LD}^2{\sigma }_{SL}^2+{\sigma }_{LD}^2{\mu }_{SL}^2+{\sigma }_{SL}^2{\mu }_{LD}^2+{\mu }_{SL}^2{\mu }_{LD}^2){)}^2\hfill & i=k\hfill \end{array}\right.\end{array}$$ (A63)$$\mathrm{var}({\gamma }_D)=M\mathrm{var}(Z_i)+M(M-1)\mathrm{cov}(Z_i,Z_k)$$

References

1. Boulogeorgos, A.A.A.; Alexiou, A. Performance Analysis of Reconfigurable Intelligent Surface-Assisted Wireless Systems and Comparison with Relaying. IEEE Access; 2020; 8, pp. 94463-94483. [DOI: https://dx.doi.org/10.1109/ACCESS.2020.2995435]

2. Wymeersch, H.; He, J.; Denis, B.; Clemente, A.; Juntti, M. Radio Localization and Mapping with Reconfigurable Intelligent Surfaces: Challenges, Opportunities, and Research Directions. IEEE Veh. Technol. Mag.; 2020; 15, pp. 52-61. [DOI: https://dx.doi.org/10.1109/MVT.2020.3023682]

3. Yang, Y.; Zheng, B.; Zhang, S.; Zhang, R. Intelligent Reflecting Surface Meets OFDM: Protocol Design and Rate Maximization. IEEE Trans. Commun.; 2020; 68, pp. 4522-4535. [DOI: https://dx.doi.org/10.1109/TCOMM.2020.2981458]

4. Basar, E. Transmission through large intelligent surfaces: A new frontier in wireless communications. Proceedings of the 2019 European Conference on Networks and Communications (EuCNC); Valencia, Spain, 18–21 June 2019; pp. 112-117.

5. Perović, N.S.; Renzo, M.D.; Flanagan, M.F. Channel Capacity Optimization Using Reconfigurable Intelligent Surfaces in Indoor mmWave Environments. Proceedings of the ICC 2020—2020 IEEE International Conference on Communications (ICC); Dublin, Ireland, 7–11 June 2020; pp. 1-7. [DOI: https://dx.doi.org/10.1109/ICC40277.2020.9148781]

6. Elbir, A.M.; Papazafeiropoulos, A.; Kourtessis, P.; Chatzinotas, S. Deep Channel Learning for Large Intelligent Surfaces Aided mm-Wave Massive MIMO Systems. IEEE Wirel. Commun. Lett.; 2020; 9, pp. 1447-1451. [DOI: https://dx.doi.org/10.1109/LWC.2020.2993699]

7. Lin, J.; Wang, G.; Fan, R.; Tsiftsis, T.A.; Tellambura, C. Channel Estimation for Wireless Communication Systems Assisted by Large Intelligent Surfaces. arXiv; 2019; arXiv: cs.IT/1911.02158

8. Taha, A.; Alrabeiah, M.; Alkhateeb, A. Deep learning for large intelligent surfaces in millimeter wave and massive MIMO systems. Proceedings of the 2019 IEEE Global Communications Conference (GLOBECOM); Island, HI, USA, 9–13 December 2019; pp. 1-6.

9. Hu, S.; Rusek, F. Spherical Large Intelligent Surfaces. Proceedings of the ICASSP 2020—2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP); Barcelona, Spain, 4–8 May 2020; pp. 8673-8677. [DOI: https://dx.doi.org/10.1109/ICASSP40776.2020.9054011]

10. Najafi, M.; Jamali, V.; Schober, R.; Poor, V.H. Physics-based Modeling and Scalable Optimization of Large Intelligent Reflecting Surfaces. arXiv; 2020; arXiv: cs.IT/2004.12957

11. Garcia, J.C.B.; Sibille, A.; Kamoun, M. Reconfigurable Intelligent Surfaces: Bridging the Gap Between Scattering and Reflection. IEEE J. Sel. Areas Commun.; 2020; 38, pp. 2538-2547. [DOI: https://dx.doi.org/10.1109/JSAC.2020.3007037]

12. Kishk, M.A.; Alouini, M.S. Exploiting Randomly-located Blockages for Large-Scale Deployment of Intelligent Surfaces. arXiv; 2020; arXiv: cs.IT/2001.10766[DOI: https://dx.doi.org/10.1109/JSAC.2020.3018808]

13. Mukherjee, M.; Kumar, V.; Guo, M.; da Costa, D.B.; Basar, E.; Ding, Z. The Interplay of Reconfigurable Intelligent Surfaces and Mobile Edge Computing in Future Wireless Networks: A Win-Win Strategy to 6G. arXiv; 2021; arXiv: 2106.11784

14. Malandrino, F.; Nordio, A.; Chiasserini, C.F. Eavesdropping with Intelligent Reflective Surfaces: Threats and Defense Strategies. arXiv; 2021; arXiv: 2108.00149

15. Yang, L.; Yang, J.; Xie, W.; Hasna, M.O.; Tsiftsis, T.; Renzo, M.D. Secrecy Performance Analysis of RIS-Aided Wireless Communication Systems. IEEE Trans. Veh. Technol.; 2020; 69, pp. 12296-12300. [DOI: https://dx.doi.org/10.1109/TVT.2020.3007521]

16. Trigui, I.; Agbogla, E.K.; Benjillali, M.; Ajib, W.; Zhu, W.P. Bit Error Rate Analysis for Reconfigurable Intelligent Surfaces with Phase Errors. arXiv; 2021; arXiv: 2101.03162

17. Ai, Y.; de Figueiredo, F.A.P.; Kong, L.; Cheffena, M.; Chatzinotas, S.; Ottersten, B. Secure Vehicular Communications Through Reconfigurable Intelligent Surfaces. IEEE Trans. Veh. Technol.; 2021; 70, pp. 7272-7276. [DOI: https://dx.doi.org/10.1109/TVT.2021.3088441]

18. Makarfi, A.U.; Rabie, K.M.; Kaiwartya, O.; Adhikari, K.; Li, X.; Quiroz-Castellanos, M.; Kharel, R. Reconfigurable Intelligent Surfaces-Enabled Vehicular Networks: A Physical Layer Security Perspective. arXiv; 2020; arXiv: eess.SP/2004.11288

19. Coelho Ferreira, R.; Facina, M.S.; de Figueiredo, F.A.; Fraidenraich, G.; de Lima, E.R. Large Intelligent Surfaces Communicating through Massive MIMO Rayleigh Fading Channels. Sensors; 2020; 20, 6679. [DOI: https://dx.doi.org/10.3390/s20226679] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33266451]

20. Best, D.; Fisher, N.I. Efficient simulation of the von Mises distribution. J. R. Stat. Soc. Ser. (Appl. Stat.); 1979; 28, pp. 152-157. [DOI: https://dx.doi.org/10.2307/2346732]

21. Beran, R. Minimum Hellinger distance estimates for parametric models. Ann. Stat.; 1977; 5, pp. 445-463. [DOI: https://dx.doi.org/10.1214/aos/1176343842]

22. Brychkov, Y.A.; Marichev, O.; Prudnikov, A. Integrals and Series: Special Functions; American Association of Physics Teachers Maryland: Botsford, MA, USA, 1988; Volume 2.

23. Proakis, J.G. Digital Communications; 5th ed. McGraw Hill: New York, NY, USA, 2007.

24. Lin, S.H.; Lu, R.R.; Fu, X.T.; Tong, A.L.; Wang, J.Y. Physical-Layer Security Analysis over M-Distributed Fading Channels. Entropy; 2019; 21, 998. [DOI: https://dx.doi.org/10.3390/e21100998]