# Multi-RIS-aided Wireless Systems: Statistical Characterization and
Performance Analysis

Tri Nhu Do, Georges Kaddoum, Thanh Luan Nguyen, Daniel Benevides da Costa, and
Zygmunt J. Haas T. N. Do and G. Kaddoum are with the Department of Electrical
Engineering, the École de Technologie Supérieure (ÉTS), Université du Québec,
Montréal, QC H3C 1K3, Canada (emails<EMAIL_ADDRESS>georges.kaddoum@etsmtl.ca).T. L. Nguyen is with the Faculty of Electronics
Technology, Industrial University of Ho Chi Minh City, Ho Chi Minh City
700000, Vietnam (e-mail: thanhluan.nguyen.iuh@outlook.com).D. B. da Costa is
with the Intelligent Wireless Communications (IWiCom) Research Group - Future
Technology Research Center, National Yunlin University of Science and
Technology, Douliou, Yunlin 64002, Taiwan, R.O.C, and with the Department of
Computer Engineering, Federal University of Ceará (UFC), Sobral 62010-560,
Brazil (e-mail: danielbcosta@ieee.org).Z. J. Haas is with the Department of
Computer Science, University of Texas at Dallas, Richardson, TX 75080, USA,
and also with the School of Electrical and Computer Engineering, Cornell
University, Ithaca, NY 14853, USA (e-mail: zhaas@cornell.edu).

###### Abstract

In this paper, we study the statistical characterization and modeling of
distributed multi-reconfigurable intelligent surface (RIS)-aided wireless
systems. Specifically, we consider a practical system model where the RISs
with different geometric sizes are distributively deployed, and wireless
channels associated to different RISs are assumed to be independent but not
identically distributed (i.n.i.d.). We propose two purpose-oriented multi-RIS-
aided schemes, namely, the exhaustive RIS-aided (ERA) and opportunistic RIS-
aided (ORA) schemes. A mathematical framework, which relies on the method of
moments, is proposed to statistically characterize the end-to-end (e2e)
channels of these schemes. It is shown that either a Gamma distribution or a
Log-Normal distribution can be used to approximate the distribution of the
magnitude of the e2e channel coefficients in both schemes. With these
findings, we evaluate the performance of the two schemes in terms of outage
probability (OP) and ergodic capacity (EC), where tight approximate closed-
form expressions for the OP and EC are derived. Representative results show
that the ERA scheme outperforms the ORA scheme in terms of OP and EC. In
addition, under i.n.i.d. fading channels, the reflecting element settings and
location settings of RISs have a significant impact on the system performance
of both the ERA or ORA schemes.

###### Index Terms:

Reconfigurable Intelligent Surfaces (RIS), Multi-RIS, Gamma distribution, Log-
Normal distribution, method of moments, ergodic capacity, outage probability.

##  I Introduction

Reconfigurable intelligent surfaces (RISs) have recently been prospected as
one of the key technologies to achieve smart radio environment (SRE) for the
sixth generation (6G) of wireless communication systems [1]. Specifically, in
order to create a truly SRE, one of the key considered approaches is to make
the radio medium controllable [1]. To this end, RISs, also known as
intelligent reflecting surfaces (IRSs) [2] or software-defined metasurfaces
(SDMs) [3], have been developed. In particular, a RIS consists of nearly-
passive reflecting elements, which are programmable and controllable via a RIS
controller, allowing it to reflect and steer impinging signals toward desired
directions. In order to realize this goal, the phase-shifts of the reflecting
elements are adjusted, such that the directionality of the beam of scattered
signals can be controlled. With this phase-tuning capability, coherent signal
combining of reflecting signals from different elements of different RISs can
be properly implemented such that the multi-path received signals
constructively add at the receiver side [4, 3, 5].

Single-RIS-aided systems have been extensively studied in the literature [6,
7, 8, 9, 10, 11, 12]. In [6], the authors presented an analytical
representation of the ideal phase-shift configuration that maximizes the
signal-to-noise ratio (SNR) of the synthesized received signal. In [7],
considering a large intelligent surface (LIS)-aided single-input single-output
(SISO) system without the direct link, the authors showed that the magnitude
of the end-to-end (e2e) channel coefficient under Rayleigh fading can
approximately follow a Gamma distribution. In [8], similar to [7], considering
a large reflecting surface (LRS)-aided SISO system under Rayleigh fading, the
authors showed that the magnitude of the e2e channel coefficient approximately
follows a Nakagami distribution. In [9], considering a similar system setting
as in [7] and [8], the authors derived approximate closed-form expressions for
the probability density function (PDF) and cumulative distribution function
(CDF) of the magnitude of the e2e channel coefficient. In [10], considering a
single-RIS-aided system under Rician fading, the authors relied on the central
limit theorem (CLT) to show that the magnitude of the e2e channel coefficient
(without the direct channel coefficient) approximately follows a Gaussian
distribution. The aforementioned works made a common assumption that channels
associated with different reflecting elements of the same RIS are independent
and identically distributed (i.i.d.). In [11], the authors showed that i.i.d.
Rayleigh fading does not physically occur when using a single RIS in an
isotropic scattering environment. In [12, 13, 14], relying on the method of
moments, the authors showed that the distribution of the e2e channel
coefficient of single-RIS-aided systems can be approximated by a Gamma
distribution.

In this work, we focus on a more general system setup and consider multi-RIS-
aided systems. Next, we discuss some recent papers which assume such a system
configuration. In [4], considering a multi-RIS-aided SISO system with direct
link under i.i.d. Nakagami-$m$ fading channels, the authors relied on the CLT
to show that the magnitude of the e2e channel coefficient can approximately
follow a Gaussian distribution. In [15], considering an LIS-aided system under
Rician fading, where the LIS consists of RIS units and each RIS unit has a
large number of reflecting elements, the authors relied on the law of large
numbers and the CLT to approximate the distributions of random variables
(RVs), which are functions of the squared magnitude of the channel
coefficient, by Gaussian distributions. In [16], the authors considered a
multi-RIS-aided point-to-point transmission assuming that the direct link does
not exist and only the best RIS is selected to participate in the
transmission. Assuming that all the Rayleigh fading channels between different
RISs are i.i.d., the authors showed that the e2e SNR approximately follows a
non-central chi-square distribution. In [17], considering a multi-RIS-aided
system, in which the number of elements of each RIS can be arbitrarily
adjusted, the authors proposed RIS selection strategies based on the RISs’
location information. The performance analysis was carried out based on the
assumption that the magnitudes of the channel coefficients associated with
different RISs are i.i.d. RVs. In [18], considering multi-hop multi-RIS-aided
systems, the authors proposed a RIS selection strategy that maximizes the e2e
SNR. However, the authors only took into account the impact of path-loss and
ignored the impact of fading. In [19], the authors discussed the impact of the
centralized and distributed RIS deployment strategies on the capacity region
of multi-RIS-aided systems. However, in the distributed RIS deployment, the
authors assumed that the channels associated with the different distributed
RISs underwent i.i.d. Rayleigh fading. In [20], the authors considered multi-
RIS-aided systems for both indoor and outdoor communications, where a direct
link between a source and a destination is unavailable; aiming for low-
complexity transmission, the authors proposed a RIS selection strategy that
selects the RIS with the highest SNR to assist the communication. However,
small-scale fading was ignored and the performance analysis of the RIS
selection strategy was not carried out. In [21] and [22], the authors
evaluated the spatial throughput of single-cell and multi-cell systems,
respectively, in which multiple RISs are used to assist multiple users under
intra-cell and/or inter-cell interference.

Although previous works have provided important contributions to multi-RIS-
aided systems, the accurate characterization of the fading model is still an
open problem. Specifically, to facilitate the performance analysis, some
existing works only considered path-loss effects and/or ignored small-scale
fading effects, as in [18, 19, 20]. On the other hand, when taking into
account small-scale fading, existing works relied on the i.i.d. fading channel
model [4, 16], and [19] or deterministic fading channel, as in [15]. Channels
associated with different elements of the same RIS can be reasonably assumed
to be i.i.d. as aforementioned since the elements are typically in sub-
wavelength size [1] and are installed closely to each other on the same panel.
However, channels between different RISs cannot be assumed to be i.i.d.,
because in distributed multi-RIS-aided systems, the RISs are installed
significantly far apart, e.g., tens of meters.

Motivated by the aforementioned observations, in this work, we propose two
distributed multi-RIS-aided wireless schemes, namely, the exhaustive RIS-aided
(ERA) and opportunistic RIS-aided (ORA) schemes. Specifically, in the former,
all the RISs participate in the transmission, whereas in the later, only one
scheduled RIS is employed. It is shown that each approach has its own
advantages and drawbacks. In particular, the ERA scheme may give better
performance at the expense of additional complexity and fronthaul/backhaul
load. On the other hand, by using only one RIS, as in the ORA scheme, the
other RISs can be employed for other purposes (assisting other pairs of users,
for instance), consequently providing a more efficient use of the available
resources. Thus, depending on the needs and goals of the application, one
approach can be more useful than the other. Considering such multi-RIS-aided
schemes, we raise a research question that still needs to be properly answered
in the literature: what distributions should be used to statistically model
the magnitude of the e2e channel coefficient of distributed multi-RIS-aided
systems, which is a combination of direct and all reflecting channel
coefficients? Thus, to answer the raised research question comprehensively, we
consider a practical multi-RIS-aided system setting, in which channels
associated with reflecting elements of the same RIS are assumed to be i.i.d.,
whereas channels associated with different RISs are assumed to be independent
but not identically distributed (i.n.i.d.), and the system undergoes
Nakagami-$m$ fading.

The key contributions of this work are summarized as follows:

  * •

The unique technical contribution of our paper lies in our comprehensive
analysis approach. Specifically, we introduce a framework for implementing the
method of moments for statistically characterizing the e2e channel of any
wireless system, including distributed multi-RIS-aided wireless systems. Under
the considered fading environment, we prove that the distribution of the
magnitude of the e2e channel coefficient of the ERA scheme can be approximated
by either a Gamma or Log-Normal distribution. On the other hand, the magnitude
of the e2e channel coefficient of the ORA scheme can be approximated by a Log-
Normal distribution. It is noted that these findings have not yet been
reported in the literature.

  * •

For the system performance analysis, invoking these findings, tight
approximate closed-form expressions for the outage probability (OP) and
ergodic capacity (EC) of the ERA and the ORA schemes are derived, based on
which insightful discussions are drawn. In addition, we provide an in-depth
analysis of the accuracy of the proposed distribution approximation, i.e., the
accuracy in using Gamma and Log-Normal distributions to approximate the true
distribution of the e2e channel coefficient, based on the lower tail
characteristics of the Gamma and Log-Normal distributions, the squared
Gaussian RVs-based physical origin, the Kullback-Leibler (KL) divergence, and
the Kolmogorov-Smirnov (KS) goodness-of-fit test.

  * •

In our simulation, we consider realistic simulation settings, as shown in
Table I. Notably, our simulation results are in agreement with the analytical
results in all simulation settings. From an engineering (and system design)
perspective, our results reveal that the ERA scheme outperforms the ORA scheme
in terms of OP and EC; nevertheless, the ORA scheme achieves a higher energy
efficiency (EE) in a specific range of the target spectral efficiency (SE)
under a certain network setting. From the system infrastructure viewpoint, the
number of elements installed at each RIS and the locations of the RISs have a
significant impact on the performance of the proposed schemes, in which, the
ERA scheme is more robust against changes in the system infrastructure than
the ORA scheme. We show that both centralized and distributed multi-RIS-aided
systems have their own advantages and drawbacks, which depends on the pre-
planed locations of the RISs.

Notations: $\Gamma(\cdot)$ denotes the Gamma function [23, Eq. (8.310.1)],
$\gamma(\cdot,\cdot)$ denotes the lower incomplete Gamma function [23, Eq.
(8.350.1)], $K_{\nu}(z)$ denotes the modified Bessel function of the second
kind [23, Eq. (8.407.1)], $\ln(x)=\log_{e}(x)$, $\mathrm{erf}(\cdot)$ denotes
the error function [24, Eq. (7.1.1)], $G^{m,n}_{p,q}[\cdot]$ denotes the
Meijer-G function [25, Eq. (8.2.1.1)], and $\mathrm{erfc}(\cdot)$ denotes the
complementary error function [23, Eq. (8.250.4)]. $X\overset{\rm
approx.}{\sim}\mathcal{D}(\cdot)$ means $X$ approximately follows the
distribution $\mathcal{D}(\cdot)$. For parameters of distributions, e.g.,
${\alpha_{X}}$ means the parameter $\alpha$ in the distribution of random
variable $X$, $\mu_{X}(k)$ denotes the $k$-th moment of $X$, and
$F_{A}(\cdot)$ denotes the Lauricell function Type-A [26, Eq.(1.4.1)].
$\mathbb{E}[X]$ and $\mu_{X}$ denote the expectation of $X$ and the
corresponding expected value.

##  II System Model

We consider a distributed multi-RIS-aided system, consisting of one source,
$\mathrm{S}$, that communicates with one destination, $\mathrm{D}$, through a
direct link with the assistance of $N$ distributed RISs,
$\mathrm{R}_{n},n=1,\ldots,N$, as depicted in Fig. 1. It is assumed that
$\mathrm{S}$ and $\mathrm{D}$ are single-antenna nodes, whereas the $n$-th RIS
is equipped with $L_{n}$ passive reflecting elements. It is noted that the
RISs may have different geometric sizes, i.e., $L_{i}\neq
L_{j},i,j\in\\{1,\ldots,N\\}$. Let
$\boldsymbol{\mathrm{\Theta}}_{n}=\mathrm{diag}([\kappa_{n1}e^{j\theta_{n1}},\ldots,\kappa_{nl}e^{j\theta_{nl}},\ldots,\kappa_{nL_{n}}e^{j\theta_{nL_{n}}}])$
denote the phase-shift matrix [27], where $\kappa_{nl}\in(0,1]$ and
$\theta_{nl}\in[0,2\pi)$ stand for, respectively, the amplitude reflection
coefficient and the phase-shift of the $l$-th reflecting element of the $n$-th
RIS. Each RIS is programmed by a controller, where the controllers can be
connected, for example, via a high-speed backhaul to synchronize the phase-
shift configuration over the entire system. To provide a more general system
model and comprehensive performance analysis, the direct link from
$\mathrm{S}$ to $\mathrm{D}$ is assumed to be available. Accordingly, the use
of RISs aims to increase the reliability of the direct transmission by
creating additional propagation paths. Furthermore, we assume that perfect
global channel state information (CSI) is available at the RIS controllers and
the source. In order to obtain the CSI, we resort to channel estimation
techniques proposed for (multi)-RIS-aided systems in the literature, such as
in [28, 29]. It is noted that a possible deployment of the distributed multi-
RIS-aided systems is over quasi-static wireless channels in low mobility
environments, such as walking pedestrian environments. In addition, we assume
interference-free signal reflection between RISs. For instance, one possible
scenario is that RISs are deployed as uniform rectangular arrays (URAs)
implanted on a wall, which leads to no interference between the RISs.

Figure 1: Illustrations of a distributed multi-RIS-aided wireless system,
where RISs can be installed on houses and/or buildings.

Let $\tilde{h}_{nl}$ and $\tilde{g}_{nl}$ denote the complex channel
coefficients from $\mathrm{S}$ to the $l$-th reflecting element of the $n$-th
RIS, and from that reflecting element to $\mathrm{D}$, respectively; and let
$\tilde{h}_{0}$ denote channel coefficient of the direct
$\mathrm{S}\to\mathrm{D}$ link. Knowing that, for a complex number $z$, one
can represent $z=r(\cos\theta+j\sin\theta)=re^{j\theta}$, where $\theta=\angle
z$ and $r=|z|$, and $|re^{j\theta}|=|re^{-j\theta}|,\forall r\in\mathbb{R}$,
the channels can be expressed in their polar representations as
$\tilde{h}_{0}=h_{0}e^{j\phi_{0}}$, $\tilde{h}_{nl}=h_{nl}e^{j\phi_{nl}}$,
$\tilde{g}_{nl}=g_{nl}e^{j\psi_{nl}}$, where $h_{0}$, $h_{nl}$, and $g_{nl}$
are the magnitudes of the channel coefficients, i.e., $h_{0}=|\tilde{h}_{0}|$,
$h_{nl}=|\tilde{h}_{nl}|$, and $g_{nl}=|\tilde{g}_{nl}|$, and $\phi_{0}$,
$\phi_{nl}$, and $\psi_{nl}$ are the phases of $\tilde{h}_{0}$,
$\tilde{h}_{nl}$, and $\tilde{g}_{nl}$, respectively, with
$\\{\phi_{0},\phi_{nl},\psi_{nl}\\}\in[0,2\pi]$. Let $x_{\mathrm{S}}$, with
$\mathbb{E}[|x_{\mathrm{S}}|^{2}]=1$, be the transmit symbol at $\mathrm{S}$,
$P_{\mathrm{S}}$ be the transmit power (in [dBm]) of $\mathrm{S}$, and
$w_{\mathrm{D}}$ be the additive white Gaussian noise (AWGN) at $\mathrm{D}$
with zero mean and variance $\sigma_{\mathrm{D}}^{2}$, i.e.,
$w_{\mathrm{D}}\sim\mathcal{CN}(0,\sigma_{\mathrm{D}}^{2})$.

###  II-A The Exhaustive RIS-aided Scheme

In the ERA scheme, all the RISs assist the transmission between $\mathrm{S}$
and $\mathrm{D}$, i.e., when $\mathrm{S}$ transmits its signal
$x_{\mathrm{S}}$ to $\mathrm{D}$, all $N$ RISs are controlled to reflect
replicas of $x_{\mathrm{S}}$ to $\mathrm{D}$ over the same time-frequency
channel. Herein, we assume that all RISs are controlled to steer their beams
to the destination, and therefore they cannot interfere with each other, as in
[4, 16, 17], and [20]. As a result, $D$ receives a superposition (combination)
of the all received multi-path signals. Thus, using coherent combination [4]
and [12], the received signal at $\mathrm{D}$, which is synthesized from the
direct signal and reflected signals from all $N$ RISs, can be written as

$\displaystyle
y^{\mathrm{ERA}}=\sqrt{P_{\mathrm{S}}}\bigg{(}\tilde{h}_{0}+\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}\tilde{g}_{nl}\kappa_{nl}e^{j\theta_{nl}}\tilde{h}_{nl}\bigg{)}x_{\mathrm{S}}+w_{\mathrm{D}}.$
(1)

Using polar representation of the complex channel coefficients, the received
SNR at $\mathrm{D}$ can be expressed as

$\displaystyle\mathrm{SNR}^{\mathrm{ERA}}$
$\displaystyle=\bar{\rho}\bigg{|}h_{0}e^{j\phi_{0}}+\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}g_{nl}\kappa_{nl}h_{nl}e^{j(\theta_{nl}+\phi_{nl}+\psi_{nl})}\bigg{|}^{2}$
$\displaystyle=\bar{\rho}\underbrace{|e^{j\phi_{0}}|^{2}}_{=1}\bigg{|}h_{0}+\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}g_{nl}\kappa_{nl}h_{nl}e^{j\delta_{nl}}\bigg{|}^{2},$
(2)

where $\bar{\rho}=P_{\mathrm{S}}/\sigma_{\mathrm{D}}^{2}$ denotes the average
transmit SNR [dB] and
$\delta_{nl}\triangleq\theta_{nl}+\phi_{nl}+\psi_{nl}-\phi_{0}$ is the phase
error of the $l$-th reflecting element of the $n$-th RIS.

In order to overcome the destructive effect of multi-path fading, the phase-
shifts of the RISs are reconfigured such that all the received signals are
constructively added to achieve the largest SNR. Mathematically, the ideal
phase-shift configuration of the $l$-th reflecting element of the $n$-th RIS
can be expressed as

$\displaystyle\theta_{nl}^{*}=\arg\max_{\theta_{nl}\in\mathcal{Q}_{n}}\mathrm{SNR}^{\mathrm{ERA}}(\theta_{nl}),\forall
l,\forall n.$ (3)

In practice, considering discrete phase-shifts, the number of phase-shifts is
limited and constrained by a phase-shift resolution, which is the number of
quantization bits, $b_{n}$, for the $n$-th RIS. Let $Q_{n}\triangleq
2^{b_{n}}$ denote the phase-shift resolution of the $n$-th RIS. Thus, the
value of a phase-shift can be picked from the following set
$\mathcal{Q}_{n}=\left\\{0,\frac{2\pi}{Q_{n}},\frac{4\pi}{Q_{n}},\ldots,\frac{2\pi(Q_{n}-1)}{Q_{n}}\right\\}$
[30]. Assuming a high phase-shift resolution, i.e., $2^{b_{n}}\gg 1$, and
perfect channel state information (CSI), the $n$-th RIS is able to generate
ideal phase-shifts such that the phase errors can be zero, i.e.,
$\delta_{nl}=0,\forall l,\forall n$. It is noted that the ideal phase-shift
assumption has been extensively used for single-RIS-aided systems [6, 9, 31,
10], and [11], multi-RIS-aided systems [4] and [16], and multi-hop multi-RIS-
aided networks [18]. As a result,
$\theta_{nl}^{*}=\phi_{0}-(\phi_{nl}+\psi_{nl})$, which yields a synthesized
received signal at $\mathrm{D}$ having the largest amplitude. Note that since
$e^{j(\theta+2k\pi)}=e^{j\theta}$, $\theta_{nl}^{*}\in[0,2\pi)$. Thus, the
received SNR at $\mathrm{D}$ can be re-expressed as

$\displaystyle\mathrm{SNR}^{\mathrm{ERA}}=\bar{\rho}\bigg{|}h_{0}+\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}\kappa_{nl}h_{nl}g_{nl}\bigg{|}^{2}.$
(4)

###  II-B The Opportunistic RIS-aided Scheme

Instead of using all the RISs, we propose the ORA scheme with the aim to
reduce the resource usage. In this case, only the most appropriate RIS
participates in assisting the direct transmission. The benefit of the ORA
scheme is that it provides a low-complexity and energy-efficient transmission
protocol at the receiver side. Consequently, the receiver in the ORA scheme
handles a significantly lower number of reflecting signals compared to the ERA
scheme. Specifically, assuming that the $n$-th RIS is scheduled to assist the
direct transmission, the received signal at $\mathrm{D}$ can be written as

$\displaystyle
y^{\mathrm{ORA}}_{n}=\sqrt{P_{\mathrm{S}}}\bigg{(}\tilde{h}_{0}+\sum_{l=1}^{L_{n}}\tilde{g}_{nl}\kappa_{nl}e^{j\theta_{nl}}\tilde{h}_{nl}\bigg{)}x_{\mathrm{S}}+w_{\mathrm{D}}.$
(5)

The received SNR associated with the $n$-th RIS can be expressed as

$\displaystyle\mathrm{SNR}^{\mathrm{ORA}}_{n}=\bar{\rho}\underbrace{|e^{j\phi_{0}}|^{2}}_{=1}\bigg{|}h_{0}+\sum_{l=1}^{L_{n}}g_{nl}\kappa_{nl}h_{nl}e^{j\delta_{nl}}\bigg{|}^{2},$
(6)

where $\delta_{nl}=\theta_{nl}+\phi_{nl}+\psi_{nl}-\phi_{0}$. In the ORA
scheme, the selected RIS is the one that potentially provides the highest e2e
received SNR at $\mathrm{D}$. Specifically, the opportunistic scheduling
criterion to choose the most appropriate RIS can be mathematically expressed
as

$\displaystyle n^{*}=\arg\max_{1\leq n\leq
N}\max_{\theta_{nl}\in\mathcal{Q}_{n}}\mathrm{SNR}^{\mathrm{ORA}}_{n},\forall
l,\forall n.$ (7)

In order to meet (7), the ideal phase-shift configuration is first estimated
at each reflecting element of each RIS, i.e.,
$\theta_{nl}^{*}=\phi_{0}-(\phi_{nl}+\psi_{nl}),\forall l,\forall n$. Thus,
the criterion (7) can be rewritten as $n^{*}=\arg\max_{1\leq n\leq
N}\mathrm{SNR}^{\mathrm{ORA}}_{n}$, i.e.,

$\displaystyle\mathrm{SNR}^{\mathrm{ORA}}_{n^{*}}$ $\displaystyle=\max_{1\leq
n\leq N}\mathrm{SNR}^{\mathrm{ORA}}_{n}$ $\displaystyle=\max_{1\leq n\leq
N}\bigg{\\{}\bar{\rho}\bigg{|}h_{0}+\sum_{l=1}^{L_{n}}g_{nl}\kappa_{nl}h_{nl}\bigg{|}^{2}\bigg{\\}}.$
(8)

As can be observed, the operation of the ORA scheme requires the same number
of CSI estimations as required by the ERA scheme. Nevertheless, the
destination node in the ORA scheme just needs to process $(L_{n}+1)$ received
signals compared to $(N\times L_{n}+1)$ signals in the ERA scheme. This
difference in the participating passive elements results in a significant gap
between the EE of the two schemes.

##  III Performance Analysis of the ERA Scheme

For notational simplicity, along the analysis, let
$U_{nl}\triangleq\kappa_{nl}h_{nl}g_{nl}$,
$V_{n}\triangleq\sum_{l=1}^{L_{n}}U_{nl}$, $M_{V}\triangleq\max_{1\leq n\leq
N}V_{n}$, $T\triangleq\sum_{n=1}^{N}V_{n}$, $Z\triangleq h_{0}+T$, and
$R\triangleq h_{0}+M_{V}$. The magnitudes of the e2e channel coefficients in
the ERA and the ORA schemes can be represented by the RVs $Z$ and $R$,
respectively. Directly deriving true distributions of these RVs is infeasible
because of their complicated structures. To deal with this technical problem,
we propose a methodology to determine an approximate version of the true
distribution.

###  III-A The Proposed Framework for Statistical Characterization of the
Magnitude of the e2e Channel

The 3-step framework is proposed as follows:

  * •

Step 1: We first determine which parametric distribution is the best candidate
to approximate the true distribution. Specifically, we heuristically and
numerically match the simulated true distribution to some known candidate
distributions, e.g., Gamma distribution, Log-Normal distribution, Gaussian
distribution, Burr distribution, Weibull distribution, just to name of few, by
using the $\mathrm{fitdist}$ function of Matlab [32]. We then select the best
candidate distributions, which give a high goodness of fit and accurate
simulation results of performance metrics of interests, e.g., outage
probability (OP) and/or ergodic capacity (EC), as will be shown in Theorems 1,
2, and 4. If no appropriate candidate distribution can be found, we focus on
matching the simulated true distribution of a key component of the considered
RV to the known candidate distribution, as will be shown in Theorem 3.

  * •

Step 2: We next determine the statistical characteristics, i.e., parameters,
of the candidate parametric distribution using the method of moments [33]. In
particular, based on the knowledge of the statistical characteristics of the
true distribution of individual components, e.g., the $k$-th moment of the
Nakagami-$m$ distribution of each individual channel, we determine the
probability distribution and the associated parameters of the candidate
distribution that best fits the statistical characteristics of the magnitude
of the e2e channel coefficient. The key derivation is to determine the method
of moments estimators of the parameters of the candidate distribution. More
specifically, to derive the estimators, we match the first $k$ moments of the
candidate distribution, which are unknown, with that of the simulated true
distribution, which are known as population moments and can be derived, to
form a system of equations that can be solved to find the representation of
the estimators. For instance, if the candidate distribution of the magnitude
of the e2e channel coefficient is postulated to be a parametric Gamma
distribution with the PDF given in (11), the tricky part is how to determine
the estimators for the parameters $\alpha$ and $\beta$ of this Gamma
distribution, e.g., solve the system of equations and find the $k$-th moment
of the simulated true distribution, such that this distribution can accurately
characterize the true distribution of the magnitude of the e2e channel
coefficient.

  * •

Step 3: We verify the accuracy of the obtained approximate distribution.
Specifically, in order to evaluate the accuracy of the obtained statistical
model, we rely on the corroboration between the simulated true distribution
and the matched distribution.

We elaborate more on the advantages of the proposed framework as follows.
According to the definition of i.i.d. RVs, as in [34, pp. 239], the i.i.d. RVs
have the same distribution with the same statistical characteristics. For
instance, let $X_{1}$ and $X_{2}$ are Nakagami-$m$ independent RVs, i.e.,
$X_{1}\sim\mathrm{Nakagami}(m_{1},\Omega_{1})$ and
$X_{2}\sim\mathrm{Nakagami}(m_{2},\Omega_{2})$. If $X_{1}$ and $X_{2}$ are
i.i.d. RVs, then $m_{1}=m_{2}=m$ and $\Omega_{1}=\Omega_{2}=\Omega$. If
$X_{1}$ and $X_{2}$ are i.n.i.d. RVs, then there is no restriction between
$m_{1}$ and $m_{2}$, meaning that $m_{1}$ and $m_{2}$ can be identical
($m_{1}=m_{2}$) or non-identical ($m_{1}\neq m_{2}$). This relation also
applies to $\Omega_{1}$ and $\Omega_{2}$. Directly applying the method of
moments to either $Z$ or $R$ is infeasible because of their complicated
structures. For instance, in the ERA scheme, $Z$ is the sum of $N$ i.n.i.d.
RVs, where each of these RVs consists of $2L_{n}$ i.i.d. RVs. Our proposal
framework is more flexible in the sense that it can apply the method of
moments to any small parts of $Z$ or $R$, and gradually leads to an ultimate
approximate distribution of $Z$ or $R$.

It is noteworthy that, for large-scale systems with co-channel inter-cell
interference, our proposed framework can be used as follows. Instead of
dealing directly with the e2e signal-to-interference-plus-noise ratio (SINR),
we deal with individual RVs in the numerator and denominator of the SINR
expression. Specifically, we use the framework to find the approximate
distributions of the true distributions of these individual RVs. The method of
moments can be iteratively applied to the approximate distributions, leading
to the ultimate approximate distribution, which facilitates the performance
analysis of interest. It is noted that the more distribution fittings are
performed, the more inaccurate the approximate distribution is. Nevertheless,
this inaccuracy can be compensated by increasing the degree of freedom of the
e2e channel, e.g., by installing more reflecting elements.

Next, we present some distributions, which will be used along the performance
analysis. Let $X$ be a RV following a Nakagami-$m$ distribution with its PDF
and CDF, parameterized by $m$ and $\Omega$, given by [34]

$\displaystyle f_{X}(x;m,\Omega)$
$\displaystyle=\frac{2{m}^{m}}{\Gamma(m){\Omega}^{m}}x^{2m-1}e^{-\frac{m}{\Omega}x^{2}},$
(9) $\displaystyle F_{X}(x;m,\Omega)$
$\displaystyle=\frac{\gamma\left(m,\frac{m}{\Omega}x^{2}\right)}{\Gamma(m)}.$
(10)

Here, $m>0$ is the shape parameter, indicating the severity of fading and
$\Omega>0$ is the spread parameter of the distribution. Next, we use the
alternative representation for denoting a Nakagami-$m$ RV:
$X\sim\mathrm{Nakagami}(m,\Omega)$. It is noted that $\Omega$ denotes the mean
square value of $X$, i.e., $\Omega=\mathbb{E}[X^{2}]$ [35], which is
equivalent to the average channel (power) gain. The distribution of the
magnitude of each individual channel can be expressed as
$h_{0}\sim\mathrm{Nakagami}(m_{0},\Omega_{0})$,
$h_{nl}\sim\mathrm{Nakagami}(m_{\mathrm{h}_{n}},\Omega_{\mathrm{h}_{n}})$, and
$g_{nl}\sim\mathrm{Nakagami}(m_{\mathrm{g}_{n}},\Omega_{\mathrm{g}_{n}})$,
where $l=1,\ldots,L_{n}$ and $n=1,\ldots,N$.

Let $Y$ be a RV following a Gamma distribution, whose PDF and CDF,
parameterized by $\alpha$ and $\beta$, respectively given by [34]

$\displaystyle f_{Y}(y;\alpha,\beta)$
$\displaystyle=\frac{{\beta}^{\alpha}}{\Gamma(\alpha)}y^{\alpha-1}e^{-\beta
y},y\geq 0,$ (11) $\displaystyle F_{Y}(y;\alpha,\beta)$
$\displaystyle=\frac{\gamma(\alpha,\beta y)}{\Gamma(\alpha)},y\geq 0.$ (12)

Here, $\alpha>0$ is the shape parameter and $\beta>0$ is the rate parameter of
the distribution. Hereafter, we use the following representation to denote a
Gamma RV: $Y\sim\mathrm{Gamma}(\alpha,\beta)$.

Let $W$ be a RV following a Log-Normal distribution, whose PDF and CDF are
given by [34]

$\displaystyle f_{W}\left(w;\nu,\zeta\right)$
$\displaystyle=\frac{1}{w\sqrt{2\pi\zeta^{2}}}e^{-\frac{\left(\ln
w-\nu\right)^{2}}{2\zeta^{2}}},$ (13) $\displaystyle
F_{W}\left(w;\nu,\zeta\right)$
$\displaystyle=\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{\ln
w-\nu}{\sqrt{2\zeta^{2}}}\right),$ (14)

respectively, where $\nu$ and $\zeta^{2}$, with $\zeta>0$, are the mean and
the variance of the distribution of $W$. Hereafter, we use the following
representation to denote a Log-Normal RV:
$W\sim\mathrm{LogNormal}(\nu,\zeta)$.

###  III-B Statistical Channel Characterization of the ERA Scheme Based on
Gamma Distribution

Relying on the proposed distribution estimation framework, we show that the
true distribution of $Z$ can be accurately approximated by the Gamma
distribution, as shown in Theorem 1.

######  Theorem 1.

The true distribution of $Z$ can be approximated by the Gamma distribution,
which is characterized by two parameters $\alpha_{Z}$ and $\beta_{Z}$, i.e.,
$Z\overset{\rm approx.}{\sim}\mathrm{Gamma}({\alpha_{Z}},{\beta_{Z}})$, where
the estimators of $\alpha_{Z}$ and $\beta_{Z}$ can be expressed as

$\displaystyle{\alpha_{Z}}$
$\displaystyle=\frac{(\mathbb{E}[Z])^{2}}{\mathrm{Var}[Z]}=\frac{[\mu_{Z}(1)]^{2}}{\mu_{Z}(2)-[\mu_{Z}(1)]^{2}},$
(15) $\displaystyle{\beta_{Z}}$
$\displaystyle=\frac{\mathbb{E}[Z]}{\mathrm{Var}[Z]}=\frac{\mu_{Z}(1)}{\mu_{Z}(2)-[\mu_{Z}(1)]^{2}},$
(16)

respectively, where $\mu_{Z}(1)$ and $\mu_{Z}(2)$ are presented in (85) and
(86), respectively. Thus, the approximate PDF and CDF of $Z$, i.e.,
$f_{Z}(z;{\alpha_{Z}},{\beta_{Z}})$ and $F_{Z}(z;{\alpha_{Z}},{\beta_{Z}})$,
can be expressed using (11) and (12), respectively.

###### Proof:

The proof is provided in Appendix A. ∎

It is worth noting that the tricky part is to determine the statistical
characteristics of the Gamma distribution, i.e., $\mu_{Z}(1)$ and
$\mu_{Z}(2)$. In addition, knowing that for arbitrary $X$ and $Y$, where
$Y=cX^{2}$, we have $F_{Y}(y)=F_{X}(\sqrt{y/c})$ and
$f_{Y}(y)=\frac{1}{2\sqrt{cy}}f_{X}\left(\sqrt{\frac{y}{c}}\right)$, the PDF
and CDF of $Z^{2}$ can be obtained as

$\displaystyle f_{Z^{2}}(x)$
$\displaystyle\approx\frac{1}{2\sqrt{x}}\frac{{\beta_{Z}}^{\alpha_{Z}}}{\Gamma(\alpha_{Z})}\\!\left(\sqrt{x}\right)^{\alpha_{Z}-1}\\!\\!e^{-\beta_{Z}\sqrt{x}},$
(17) $\displaystyle F_{Z^{2}}(x)$
$\displaystyle\approx\frac{\gamma\left(\alpha_{Z},\beta_{Z}\sqrt{x}\right)}{\Gamma(\alpha_{Z})},$
(18)

respectively. On the other hand, the PDF and CDF of the Generalized Gamma (GG)
distribution are, respectively, given by [36]

$\displaystyle f_{\mathrm{GG}}(x;a,d,p)$
$\displaystyle=\frac{\left(p/a^{d}\right)x^{d-1}e^{-\left(\frac{x}{a}\right)^{p}}}{\Gamma\left(d/p\right)},$
(19) $\displaystyle F_{\mathrm{GG}}(x;a,d,p)$
$\displaystyle=\frac{\gamma\left(d/p,\left(x/a\right)^{p}\right)}{\Gamma(d/p)},$
(20)

where $a,d,p$ are parameters of the distribution. Thus, (17) can be rewritten
using (19) as

$\displaystyle f_{Z^{2}}(x)$
$\displaystyle=\frac{[1/2]/\left[\left(\beta_{Z}^{-2}\right)^{\frac{\alpha_{Z}}{2}}\right]}{\Gamma([\alpha_{Z}/2]/[1/2])}x^{\frac{\alpha_{Z}}{2}-1}e^{-\left(\frac{x}{\beta_{Z}^{-2}}\right)^{\frac{1}{2}}}.$
(21)

Next, we present the statistical characterization of the e2e channel power
gain of the ERA scheme. By comparing (19) and (21), the distribution functions
of $Z^{2}$ can be represented in the form of the Generalized Gamma
distribution as in the following Remark 1.

######  Remark 1.

From Theorem 1, the true distribution of $Z^{2}$ can be approximated by using
the Generalized Gamma distribution, where the PDF and CDF of $Z^{2}$ can be
expressed as

$\displaystyle f_{Z^{2}}(x)$
$\displaystyle=f_{\mathrm{GG}}\left(x;a_{Z^{2}},d_{Z^{2}},p_{Z^{2}}\right),$
(22) $\displaystyle F_{Z^{2}}(x)$
$\displaystyle=F_{\mathrm{GG}}\left(x;a_{Z^{2}},d_{Z^{2}},p_{Z^{2}}\right),$
(23)

respectively, where $f_{\mathrm{GG}}(\cdot)$ and $F_{\mathrm{GG}}(\cdot)$ are
presented in (19) and (20), respectively, $a_{Z^{2}}=\beta_{Z}^{-2}$,
$d_{Z^{2}}={\alpha_{Z}}/2$, and $p_{Z^{2}}=1/2$, where ${\alpha_{Z}}$ and
${\beta_{Z}}$ are expressed in (15) and (16), respectively.

####  III-B1 Outage Probability

The OP can be defined as the probability that the instantaneous mutual
information of the ERA scheme drops below a target spectral efficiency (SE)
threshold $R_{\rm th}$ [b/s/Hz], i.e., $P_{\rm
out}^{\mathrm{ERA}}=\Pr(\log_{2}(\mathrm{SNR}^{\mathrm{ERA}}+1)<R_{\rm th})$.
From (4), and invoking Theorem 1, an approximate closed-form expression for
the OP of the ERA scheme can be obtained as

$\displaystyle P_{\rm
out}^{\mathrm{ERA},\mathrm{Gam}}=\Pr\left(\bar{\rho}Z^{2}\leq\rho_{\rm
th}\right)\overset{\rm(a)}{\approx}F_{Z}\left(\sqrt{\rho_{\rm
th}/\bar{\rho}}\right),$ (24)

where $\rho_{\rm th}=2^{R_{\rm th}}-1$ and step (a) in (24) is due to the fact
that $F_{X^{2}}(x)=F_{X}(\sqrt{x}),x>0$.

####  III-B2 Ergodic Capacity

The system ergodic capacity can be determined by averaging the system
instantaneous capacity over a large number of channel realizations. Invoking
Theorem 1, an approximate closed-form expression for the EC of the ERA scheme
can be expressed as

$\displaystyle\bar{C}^{\mathrm{ERA},\mathrm{Gam}}=\mathbb{E}\left[\log_{2}\left(1+\mathrm{SNR}^{\mathrm{ERA}}\right)\right]$
$\displaystyle\quad=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}Z^{2}\right)\right]$
$\displaystyle\quad\approx\frac{1}{\Gamma({\alpha_{Z}})\ln
2}\frac{2^{{\alpha_{Z}}-1}}{\sqrt{\pi}}G^{5,1}_{3,5}\left(\frac{({\beta_{Z}})^{2}}{4\bar{\rho}}\middle|\begin{matrix}0,\frac{1}{2},1\\\
\frac{{\alpha_{Z}}}{2},\frac{{\alpha_{Z}}+1}{2},0,\frac{1}{2},0\end{matrix}\right).$
(25)

The detailed derivation of (25) is provided in Appendix B.

###  III-C Statistical Channel Characterization of the ERA Scheme Based on
Log-Normal Distribution

To give more insights into the ERA scheme’s fading model, we show that the
Log-Normal distribution can alternatively be used to approximate the true
distribution of the magnitude of the e2e channel coefficient of the ERA
scheme.

######  Theorem 2.

The true distribution of $Z$ can be approximated by the Log-Normal
distribution, which is characterized by two parameters $\nu_{Z}$ and
$\zeta_{Z}$, i.e., $Z\overset{\rm
approx.}{\sim}\mathrm{LogNormal}(\nu_{Z},\zeta_{Z})$, where the estimators of
$\nu_{Z}$ and $\zeta_{Z}$ are given in (26) and (27), respectively.

###### Proof:

Since the Log-Normal distribution is characterized by two parameters $\nu$ and
$\zeta$, following Step 2 of our framework, we need to match the first two
moments of $Z$, i.e., solving the following system of equations:
$\mathbb{E}[Z]=\mathbb{E}[X]$ and $\mathbb{E}[Z^{2}]=\mathbb{E}[X^{2}]$, where
$X$ is a Log-Normal RV. As a result, the estimators of $\nu_{Z}$ and
$\zeta_{Z}$ can be written as [37]

$\displaystyle\nu_{Z}$
$\displaystyle=\ln\left(\frac{\left(\mathbb{E}[Z]\right)^{2}}{\sqrt{\mathbb{E}[(Z)^{2}]}}\right)=\ln\left(\frac{\left(\mu_{Z}(1)\right)^{2}}{\sqrt{\mu_{Z}(2)}}\right),$
(26) $\displaystyle\zeta_{Z}$
$\displaystyle=\sqrt{\ln\left(\frac{\mathbb{E}[Z^{2}]}{(\mathbb{E}[Z])^{2}}\right)}=\sqrt{\ln\left(\frac{\mu_{Z}(2)}{(\mu_{Z}(1))^{2}}\right)}.$
(27)

Next, in order to determine the exact values of $\nu_{Z}$ and $\zeta_{Z}$, we
derive a general expression for the $k$-th moment of $Z$, i.e.,
$\mu_{Z}(k)=\mathbb{E}[Z^{k}]$. After some mathematical manipulations,
$\mu_{Z}(k)$ can be obtained as

$\displaystyle\mu_{Z}(k)$
$\displaystyle=\sum_{l=0}^{k}\binom{k}{l}\mathbb{E}[(h_{0})^{l}]\mathbb{E}[T^{k-l}]$
$\displaystyle=\sum_{l=0}^{k}\binom{k}{l}\mu_{h_{0}}(l)\mu_{T}(k-l),$ (28)

where $\mu_{h_{0}}(k)$ and $\mu_{T}(k)$ are derived in (69) and (A),
respectively. Consequently, the PDF, $f_{Z}(x;\nu_{Z},\zeta_{Z})$, and the
CDF, $F_{Z}(x;,\nu_{Z},\zeta_{Z})$, can be expressed using (13) and (14),
respectively. This completes the proof of Theorem 2. ∎

####  III-C1 Outage Probability

An approximate closed-form expression for the OP of the ERA scheme can be
obtained as

$\displaystyle P_{\rm out}^{\mathrm{ERA},\rm LN}$
$\displaystyle=\Pr\left(\bar{\rho}Z^{2}\leq\rho_{\rm th}\right)$
$\displaystyle=F_{Z}\left(\sqrt{\frac{\rho_{\rm th}}{\bar{\rho}}}\right)$
$\displaystyle\approx\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{\ln\left(\sqrt{\frac{\rho_{\rm
th}}{\bar{\rho}}}\right)-\nu_{Z}}{\sqrt{2\zeta_{Z}^{2}}}\right).$ (29)

####  III-C2 Ergodic Capacity

The EC of the ERA scheme can be determined as
$\bar{C}^{\mathrm{ERA},\mathrm{LN}}=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}Z^{2}\right)\right]$.

######  Lemma 1.

Given $X\sim\mathrm{LogNormal}(\nu_{X},\zeta_{X})$, and let
$\bar{C}_{X}=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}X\right)\right]$, it
can be shown that $\bar{C}_{X}=\Upsilon(\nu_{X},\zeta_{X})$, where

$\displaystyle\Upsilon(\nu_{X},\zeta_{X})$
$\displaystyle\triangleq\frac{1}{\ln(2)}\bigg{[}\Xi\left(\frac{1}{\zeta\sqrt{2}},\frac{\ln(\bar{\rho})+\nu}{\zeta\sqrt{2}}\right)$
$\displaystyle\quad+\Xi\left(\frac{1}{\zeta\sqrt{2}},-\frac{\ln(\bar{\rho})+\nu}{\zeta\sqrt{2}}\right)+\frac{\zeta}{\sqrt{2\pi}}e^{-\frac{(\ln(\bar{\rho})+\nu)^{2}}{2\zeta^{2}}}$
$\displaystyle\quad+\frac{\ln(\bar{\rho})+\nu}{2}\mathrm{erfc}\left(-\frac{\ln(\bar{\rho})+\nu}{\zeta\sqrt{2}}\right)\bigg{]}.$
(30)

Here, $\Xi(\cdot,\cdot)$ is presented in (34) and derived in (35).

###### Proof:

We have

$\displaystyle\bar{C}_{X}$
$\displaystyle=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}X\right)\right]=\int_{0}^{\infty}\log_{2}(1+x)\frac{1}{\bar{\rho}}f_{X}\left(\frac{x}{\bar{\rho}}\right)dx.$
(31)

Substituting (13) into (31) yields

$\displaystyle\bar{C}_{X}$
$\displaystyle=\frac{1}{\bar{\rho}\zeta\sqrt{2\pi}}\int_{0}^{\infty}\frac{\log_{2}(1+x)}{x/\bar{\rho}}e^{-\frac{(\ln(x/\bar{\rho})-\nu)^{2}}{2\zeta^{2}}}dx.$
(32)

Since $\log_{2}(1+x)=\ln(1+x)/\ln(2)$, we have

$\displaystyle\bar{C}_{X}$
$\displaystyle=\frac{1}{\zeta\sqrt{2\pi}\ln(2)}\int_{0}^{\infty}\frac{\ln(1+x)}{x}e^{-\frac{[\ln(x)-(\ln(\bar{\rho})+\nu)]^{2}}{2\zeta^{2}}}dx.$
(33)

Let $\Xi(\cdot,\cdot)$ be defined as

$\displaystyle\Xi(a,b)\triangleq\frac{a}{\sqrt{\pi}}\int_{0}^{1}\frac{\ln(1+x)}{x}e^{-[a\ln(x)-b]^{2}}dx.$
(34)

From the fact that $\ln(1+x)=\sum_{k=1}^{8}{a_{k}x^{k}+\epsilon(x)},~{}0<x<1$
[38], next, using the identity [23, Eq. (8.252.6)], and then using the fact
that $1-{\rm erf}(x)={\rm erfc}(x)$ [23, Eq. (8.250.4)], after some
mathematical manipulations, $\Xi(a,b)$ in (34) can be approximated by

$\displaystyle\Xi(a,b)\approx\frac{e^{b^{2}}}{2}\sum_{k=1}^{8}a_{k}\exp\left(\left(b+\frac{k}{2a}\right)^{2}\right){\rm
erfc}\left(b+\frac{k}{2a}\right).$ (35)

Using (35), and following [38, Subsection II-A], one can arrive at the
approximate closed-form expression for $\Upsilon(\nu_{X},\zeta_{X})$ in (1).
This completes the proof of Lemma 1. ∎

Based on $Z\overset{\rm approx.}{\sim}\mathrm{LogNormal}(\nu_{Z},\zeta_{Z})$
as stated in Theorem 2, we can derive the PDF of $Z^{2}$, however, the
resultant PDF of $Z^{2}$ using this method leads to a non-closed-form
expression of the integral in (31). To circumvent this problem, we perform
Step $2$ of the proposed statistic characterization framework again, and after
some mathematical manipulation steps, we obtain $Z^{2}\overset{\rm
approx.}{\sim}\mathrm{LogNormal}(\nu_{Z^{2}},\zeta_{Z^{2}})$, where
$\nu_{Z^{2}}=\ln\left(\frac{\left(\mu_{Z}(2)\right)^{2}}{\sqrt{\mu_{Z}(4)}}\right)$,
$\zeta_{Z^{2}}=\sqrt{\ln\left(\frac{\mu_{Z}(4)}{(\mu_{Z}(2))^{2}}\right)}$.
From (31), an approximate closed-form expression for the EC of the ERA scheme
can be obtained as

$\displaystyle\bar{C}^{\mathrm{ERA},\rm
LN}=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}Z^{2}\right)\right]\approx\Upsilon(\nu_{Z^{2}},\zeta_{Z^{2}}).$
(36)

##  IV Performance Analysis of the ORA Scheme

###  IV-A Statistical Channel Characterization of the ORA Scheme

We first rely on the proposed framework to match the true distribution of the
key component of $R$, i.e., $V_{n}$, to a Gamma distribution, as stated in
(77), and then we derive approximate closed-form expressions for the PDF and
CDF of $R$. Specifically, since $h_{0}$ and $V_{n}$ are independent RVs and
$h_{0},V_{n}\geq 0$, the opportunistic scheduling criterion in (7) can be re-
expressed as

$\displaystyle n^{*}=\arg\max_{1\leq n\leq
N}\left\\{\bar{\rho}\left[h_{0}+V_{n}\right]^{2}\right\\}=\arg\max_{1\leq
n\leq N}V_{n}.$ (37)

Thus, the e2e SNR of the ORA scheme in (II-B) can be rewritten as

$\displaystyle\mathrm{SNR}^{\mathrm{ORA}}_{n^{*}}=\bar{\rho}\left[h_{0}+M_{V}\right]^{2}.$
(38)

The statistical characterization of the magnitude of the e2e channel
coefficient of the ORA scheme is presented in the following theorem.

######  Theorem 3.

Approximate closed-form expressions for the CDF and PDF of $R$ can be,
respectively, obtained as

$\displaystyle F_{R}(x)\approx\Phi(x)$
$\displaystyle\triangleq\sum_{m=1}^{M}\left[F_{h_{0}}\left(\frac{m}{M}x\right)-F_{h_{0}}\left(\frac{m-1}{M}x\right)\right]$
$\displaystyle\quad\times F_{M_{V}}\left(\frac{M-m+1}{M}x\right).$ (39)
$\displaystyle
f_{R}(x)\\!\approx\\!\Delta(x)\\!\triangleq\\!\sum_{m=1}^{M}\bigg{[}\bigg{[}\frac{m}{M}f_{h_{0}}\left(\frac{m}{M}x\right)\\!-\\!\frac{m-1}{M}f_{h_{0}}\left(\frac{m-1}{M}x\right)\\!\bigg{]}$
$\displaystyle\times F_{M_{V}}\left(\frac{M-m+1}{M}x\right)\frac{M-m+1}{M}$
$\displaystyle\times\left(F_{h_{0}}\left(\frac{m}{M}x\right)-F_{h_{0}}\left(\frac{m-1}{M}x\right)\right)f_{M_{V}}\left(\frac{M-m+1}{M}x\right)\bigg{]}.$
(40)

where $f_{h_{0}}(\cdot)$, $F_{h_{0}}(\cdot)$, and $F_{M_{V}}(\cdot)$ are
presented in (9), (10), and (49), respectively.

###### Proof:

Since the $\max\\{\cdot\\}$ operation on $V_{n}$ changes the statistical
characteristics of its argument, $V_{n}$, we first rely on the Total
Probability Theorem [34] and the definition of conditional probability [34] to
determine the probability distribution of $R$. The CDF of $R$, $F_{R}(x)$, can
be expressed as

$\displaystyle F_{R}(x)=\Pr(R\leq x)=\Pr((h_{0}+M_{V})\leq x).$ (41)

From (37) and (38), (41) can be further expressed as

$\displaystyle F_{R}(x)$ $\displaystyle=\Pr((h_{0}+V_{n^{*}})\leq x)$
$\displaystyle=\sum_{n=1}^{N}\Pr(n^{*}=n)\Pr\left((h_{0}+V_{n^{*}})\leq x\mid
n^{*}=n\right).$ (42)

Based on the conditional probability definition [34], we have

$\displaystyle\Pr\left((h_{0}+V_{n^{*}})\leq x\mid n^{*}=n\right)$
$\displaystyle\quad=\frac{\Pr(\\{(h_{0}+V_{n^{*}})\leq
x\\}\cap\\{n^{*}=n\\})}{\Pr(n^{*}=n)}.$ (43)

Thus, (42) can be rewritten as

$\displaystyle F_{R}(x)=\sum_{n=1}^{N}\Pr(\\{(h_{0}+V_{n^{*}})\leq
x\\}\cap\\{n^{*}=n\\}).$ (44)

From (II-B) and (37), (44) can be further expressed as

$\displaystyle F_{R}(x)=\sum_{n=1}^{N}\Pr(\\{(h_{0}+V_{n})\leq
x\\}\cap\\{V_{n}=\max_{1\leq k\leq N}V_{k}\\}).$ (45)

It is noted that $V_{n}=\max_{1\leq k\leq N}V_{k},n=1,\ldots,N,$ are pairwise
disjoint events, and

$\displaystyle\sum_{n=1}^{N}\Pr\left(V_{n}=\max_{1\leq k\leq
N}V_{k}\right)=1.$ (46)

Thus, for an arbitrary event $A$, relying on the Total Probability Theorem
[34], we have

$\displaystyle\sum_{n=1}^{N}\Pr\left(\left\\{V_{n}=\max_{1\leq k\leq
N}V_{k}\right\\}\cap A\right)=\Pr(A).$ (47)

Thus, (45) can be further expressed as

$\displaystyle F_{R}(x)$ $\displaystyle=\Pr\left(\big{\\{}\max_{1\leq k\leq
N}V_{k}\leq(x-h_{0})\big{\\}}\cap\\{h_{0}<x\\}\right)$
$\displaystyle=\int_{0}^{x}F_{M_{V}}(x-y)f_{h_{0}}(y)dy.$ (48)

Since $V_{k}$, $k=1,\ldots,N$, are i.n.i.d. RVs, from (78), the CDF of $M_{V}$
can be derived as

$\displaystyle F_{M_{V}}(x)$ $\displaystyle=\Pr\left(\max_{1\leq k\leq
N}V_{k}\leq x\right)$
$\displaystyle=\prod_{k=1}^{N}F_{V_{k}}(x)\approx\prod_{k=1}^{N}\frac{\gamma(L_{k}{\alpha_{U}}_{k},{\beta_{U}}_{k}z)}{\Gamma\left(L_{k}{\alpha_{U}}_{k}\right)}.$
(49)

Inserting (49) into (48) results in an integral that cannot be expressed in a
closed-form. To address this problem, we further express (48) as

$\displaystyle
F_{R}(x)=\int_{0}^{x}\int_{0}^{x-y}f_{M_{V}}(z)f_{h_{0}}(y)dzdy.$ (50)

Using the $M$-staircase approximation, (50) can be derived as

$\displaystyle F_{R}(x)$
$\displaystyle\overset{M\to\infty}{\approx}\sum_{m=1}^{M}\int_{\frac{m-1}{M}x}^{\frac{m}{M}x}f_{h_{0}}(y)\int_{0}^{\frac{M-m+1}{M}x}f_{M_{V}}(z)dzdy.$
(51)

After some mathematical manipulations, (51) can be derived as (39). Next, by
taking the derivative of (39), i.e., $f_{h_{0}}(x)=\frac{dF_{h_{0}}(x)}{dx}$,
the PDF of $R$ is obtained as in (40), while the PDF of $M_{V}$,
$f_{M_{V}}(x)$, is derived by taking the derivative of (49). This completes
the proof of Theorem 3. ∎

####  IV-A1 Outage Probability

Invoking Theorem 3, an approximate closed-form expression for the OP of the
ORA scheme can be obtained as

$\displaystyle P_{\rm out}^{\mathrm{ORA},\rm Gam}$
$\displaystyle=\Pr(\bar{\rho}R^{2}\leq\rho_{\rm th})=F_{R^{2}}\left(\rho_{\rm
th}/\bar{\rho}\right)$ $\displaystyle=F_{R}\left(\sqrt{\rho_{\rm
th}/\bar{\rho}}\right)\approx\Phi\left(\sqrt{\rho_{\rm
th}/\bar{\rho}}\right),$ (52)

where $\Phi(\cdot)$ is presented in (39).

####  IV-A2 Ergodic Capacity

An analytical closed-form expression for the EC of the ORA scheme can be
expressed as

$\displaystyle\bar{C}^{\mathrm{ORA},\mathrm{Gam}}$
$\displaystyle=\mathbb{E}[\log_{2}(1+\bar{\rho}R^{2})]$
$\displaystyle=\frac{1}{\ln
2}\int_{0}^{\infty}\frac{1}{z+1}\left[1-F_{R}\left(\sqrt{\frac{z}{\bar{\rho}}}\right)\right]dz.$
(53)

With the PDF of $R$ in (39), the integral in (53) does not admit a closed-form
expression. Nevertheless, an approximate closed-form expression for the EC of
the ORA scheme is obtained in subsection IV-B2.

###  IV-B Statistical Channel Characterization of the ORA Scheme Based on Log-
Normal Distribution

Using the Gamma distribution approximation presented in subsection IV-A, the
closed-form expression for the EC of the ORA scheme cannot be derived. It
motivates us to alternatively use a Log-Normal distribution to statistically
characterize the e2e channel coefficient.

######  Theorem 4.

The true distribution of $R$ can be approximated by the Log-Normal
distribution, characterized by two parameters $\nu_{R}$ and $\zeta_{R}$, i.e.,
$R\overset{\rm approx.}{\sim}\mathrm{LogNormal}(\nu_{R},\zeta_{R})$, where the
estimators of $\nu_{R}$ and $\zeta_{R}$ are presented in (54) and (55),
respectively.

###### Proof:

By using the same mathematical reasoning as in Subsection III-C, the
parameters characterizing the Log-Normal distribution of $R$ can be obtained
as

$\displaystyle\nu_{R}$
$\displaystyle=\ln\left(\frac{\left(\mathbb{E}[R]\right)^{2}}{\sqrt{\mathbb{E}[R^{2}]}}\right)=\ln\left(\frac{\left(\mu_{R}(1)\right)^{2}}{\sqrt{\mu_{R}(2)}}\right),$
(54) $\displaystyle\zeta_{R}$
$\displaystyle=\sqrt{\ln\left(\frac{\mathbb{E}[R^{2}]}{(\mathbb{E}[R])^{2}}\right)}\\!=\\!\sqrt{\ln\left(\frac{\mu_{R}(2)}{(\mu_{R}(1))^{2}}\right)}.$
(55)

To obtain the value of $\nu_{R}$ and $\zeta_{R}$, we derive the $k$-th moment
of $R$, which can be obtained as

$\displaystyle\mu_{R}(k)=\sum_{v=0}^{k}\binom{k}{v}\mu_{h_{0}}(v)\mu_{M_{V}}(k-v),$
(56)

From the CDF of $M_{V}$ in (49), an expression for $\mu_{M_{V}}(k)$ is derived
in the following Lemma.

######  Lemma 2.

A tight approximate closed-form expression for the $k$-th moment of $M_{V}$
can be obtained as in (57) at the top of the next page.

$\displaystyle\mu_{M_{V}}(k)$
$\displaystyle=\left[\sum_{t=1}^{N}{\beta_{U}}_{t}\right]^{-k}\left[\prod_{t=1}^{N}\chi_{t}^{L_{t}\alpha_{U_{t}}}\right]\sum_{n=1}^{N}\frac{\Gamma(\Lambda+k)}{\Gamma(L_{n}\alpha_{U_{n}})}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}\frac{1}{\Gamma(L_{t}{\alpha_{U}}_{t}+1)}$
$\displaystyle\quad\times
F_{A}^{(N-1)}\bigg{[}\Lambda+k;\underbrace{1,\ldots,1}_{(N-1)\text{
terms}};\underbrace{L_{1}{\alpha_{U}}_{1}+1,\ldots,L_{N}{\alpha_{U}}_{N}+1}_{\text{exclude
}(L_{n}\alpha_{U_{n}}+1)};\underbrace{\chi_{1},\ldots,\chi_{N}}_{\text{exclude
}\chi_{n}}\bigg{]}.$ (57)

TABLE I: Simulation Parameters Parameters |  Values | Parameters | Values
---|---|---|---
$\mathrm{S}\to\mathrm{D}$ distance, $d_{\mathrm{S}\mathrm{D}}$ [m] |  100 | Avg. height of RISs’ location, $H$ [m] | 10
Number of RISs, $N$ |  $5$ | Total reflecting elements, $\sum_{n}L_{n}$ | $[125,200]$
Amplitude reflection coef., $\kappa_{nl}$ |  $1$ [27] | Bandwidth, $\rm BW$ [MHz] |  $10$ [27]
Carrier frequency, $f_{c}$ [GHz] |  $3$ [27] | Nakagami shape parameter, $m$ | $\sim\mathcal{U}[2,3]$
CDP at RIS, $\tilde{P}_{nl}$ [mW] |  $7.8$ [30] | Target SE, $R_{\rm th}$ [b/s/Hz] | $1$
CDP at $\mathrm{S}$, $\tilde{P}_{\mathrm{S}}$, and $\mathrm{D}$, $\tilde{P}_{\mathrm{D}}$ [mW] |  10 [30] | Thermal noise power density [dBm/Hz] |  $-174$ [27]
Transmit power, $P_{\mathrm{S}}$ [dBm] |  $[0,30]$ | RIS’s coordinate [m] |  $x_{\mathrm{R}_{n}}\sim\mathcal{U}[5,95]$, $y_{\mathrm{R}_{n}}\sim\mathcal{U}[1,9]$
Antenna gain $G_{\mathrm{S}}$, $G_{\mathrm{D}}$, $G_{\mathrm{R}_{n}}$ [dB] |  5 [27] | Noise figure, $\rm NF$ [dBm] |  $10$ [27]
|  |  |

###### Proof:

The proof is provided in Appendix C. ∎

Invoking Lemma 2, and then substituting (69) and (57) into (56), we obtain the
closed-form expression for $\mu_{R}(k)$. This completes the proof of Theorem
4. ∎

####  IV-B1 Outage Probability

An approximate closed-form expression for the OP of the ORA scheme can be
obtained as

$\displaystyle P_{\rm out}^{\mathrm{ORA},\mathrm{LN}}$
$\displaystyle=\Pr(\bar{\rho}R^{2}\leq\rho_{\rm
th})=F_{R}\left(\sqrt{\frac{\rho_{\rm th}}{\bar{\rho}}}\right)$
$\displaystyle\approx\frac{1}{2}+\frac{1}{2}\mathrm{erf}\left(\frac{\ln\left(\sqrt{\frac{\rho_{\rm
th}}{\bar{\rho}}}\right)-\nu_{R}}{\sqrt{2\zeta_{R}^{2}}}\right).$ (58)

####  IV-B2 Ergodic Capacity

Using the same argument made in subsection III-C2, and after some mathematical
manipulations, we obtain $R^{2}\overset{\rm
approx.}{\sim}\mathrm{LogNormal}(\nu_{R^{2}},\zeta_{R^{2}})$,
$\nu_{R^{2}}=\ln\left(\frac{\left(\mu_{R}(2)\right)^{2}}{\sqrt{\mu_{R}(4)}}\right)$,
$\zeta_{R^{2}}=\sqrt{\ln\left(\frac{\mu_{R}(4)}{(\mu_{R}(2))^{2}}\right)}$.
Invoking Lemma 1, an approximate closed-form expression for the EC of the ORA
scheme can be expressed as

$\displaystyle\bar{C}^{\mathrm{ORA},\mathrm{LN}}=\mathbb{E}\left[\log_{2}\left(1+\bar{\rho}R^{2}\right)\right]\approx\Upsilon(\nu_{R^{2}},\zeta_{R^{2}}).$
(59)

##  V Results and Discussions

In this section, we provide representative results on the PDFs and CDFs of $Z$
and $R$, as well as the OP, EC, and energy efficiency (EE) of the ERA and the
ORA schemes, not only to verify the correctness of the developed analysis, but
also to give insights into the performance of the two schemes. It is noted
that the results in this section are reproducible using our Matlab code111The
source code is available at https://github.com/trinhudo/Multi-RIS, in which we
detail the coding of the analytical results, the obtained CDFs and PDFs, and
especially the Lauricell function Type-A, $F_{A}(\cdot)$, used in (57). One of
the technical contributions of our code is that we use symbolic computations
in Matlab [32] to obtain highly accurate results, e.g., using variable-
precision arithmetic ($\mathrm{vpa}$ functions) to obtain analytical results
with $32$ significant decimal digits of precision. The reason of using
symbolic computations is that the analysis of multi-RIS-aided systems under
i.n.i.d. fading channels involves a high number of RVs, especially, when the
number of elements of each RIS is large, some non-symbolic operations of
floating-point arithmetic in Matlab are invalid and return not-a-number (NaN)
results.

Unless otherwise specified, the simulation parameters are set to the values
provided in Table I, where the equivalent noise power at $\mathrm{D}$ is
$\sigma_{\mathrm{D}}^{2}=N_{0}+10\log(\mathrm{BW})+\mathrm{NF}$ [dBm].
Considering non-line-of-sight (NLoS) condition in the 3GPP Urban Micro (UMi)
path-loss model [39, 20], and [27], the path-loss is integrated in the spread
parameter of the Nakagami-$m$ distribution as $\Omega_{\rm
XY}=G_{\mathrm{X}}+G_{\mathrm{Y}}-22.7-26\log(f_{c})-36.7\log(d_{\rm
XY}/d_{0})$ [dB], where $d_{\rm XY}$ [m] is the distance in two-dimensional
Cartesian coordinate system, where $\rm
X\in\\{\mathrm{S},\mathrm{R}_{1},\ldots,\mathrm{R}_{N}\\}$, $\rm
Y\in\\{\mathrm{R}_{1},\ldots,\mathrm{R}_{N},\mathrm{D}\\}$, and $d_{0}$ [m] is
the reference distance, herein, $d_{0}=1$ m. Let
$\boldsymbol{\mathrm{L}}\in\mathbb{R}^{1\times N}$ denote a vector that
indicates the number of reflecting elements of each one of the $N$ RISs, i.e.,
$\boldsymbol{\mathrm{L}}=[L_{n}:n=1,\ldots,N]$. Specifically, we consider four
reflecting element setting as follows:
$\boldsymbol{\mathrm{L}}_{1}=[25,25,25,25,25]$, is used in Figs. 2-4,
$\boldsymbol{\mathrm{L}}_{2}=[40,40,40,40,40]$,
$\boldsymbol{\mathrm{L}}_{3}=[20,30,40,50,60]$, and
$\boldsymbol{\mathrm{L}}_{4}=[60,50,40,30,20]$, which are used in Figs. 5-7.
Let $\boldsymbol{\mathrm{D}}\in\mathbb{R}^{1\times
2N}=[(x_{\mathrm{R}_{n}},y_{\mathrm{R}_{n}}):n=1,\ldots,N]$ represent the
coordinates of the RISs. The locations of the source and the destination are
set as $(0,0)$ and $(100,0)$, respectively, while we consider two location
settings for RISs, namely
$\boldsymbol{\mathrm{D}}_{1}=[(7,2),(13,6),(41,8),(75,4),(93,3)]$, used in
Figs. 2-7, and
$\boldsymbol{\mathrm{D}}_{2}=[(5,2),(13,7),(37,6),(69,1),(91,3)]$, used in
Fig. 7. The number of steps, $M$, used in the $M$-staircase approximation, is
set as $M=100$.

(a)

(b)

Figure 2: Graphical demonstration of the estimated true distributions and the
obtained approximate distributions in (a) Theorems 1 and 3 and (b) Theorems 2
and 4.

In Figs. 2-4, we verify the correctness of our developed analysis.
Specifically, in Fig. 2, we plot the approximate Gamma and Log-Normal
distributions in the ERA and ORA schemes, which are presented in Theorems 1-4.
As can be seen in Fig. 2, the approximate CDFs and PDFs are in agreement with
the true CDF and PDF, which are estimated numerically based on simulation
data. Next, in Fig. 3, we plot the OP and EC results of the ERA and the ORA
schemes, i.e., $P_{\rm out}^{\mathrm{ERA},\mathrm{Gam}}$,
$\bar{C}^{\mathrm{ERA},\mathrm{Gam}}$, $P_{\rm out}^{\mathrm{ERA},\rm LN}$,
$\bar{C}^{\mathrm{ERA},\rm LN}$, $P_{\rm out}^{\mathrm{ORA},\rm Gam}$,
$\bar{C}^{\mathrm{ORA},\mathrm{Gam}}$, $P_{\rm
out}^{\mathrm{ORA},\mathrm{LN}}$, and $\bar{C}^{\mathrm{ORA},\mathrm{LN}}$
presented in (24), (25), (29), (36), (52), (53), (58), and (59), respectively.
As can be observed in Fig. 3, the theoretical and simulation results are well
corroborated, which validates our developed analysis.

(a)

(b)

(c)

(d)

Figure 3: The OP and EC of the ERA and the ORA schemes as a function of
transmit power of $\mathrm{S}$, $P_{\mathrm{S}}$ [dBm], in particular, (a)
$P_{\rm out}^{\mathrm{ERA},\mathrm{Gam}}$ and $P_{\rm out}^{\mathrm{ORA},\rm
Gam}$, (b) $\bar{C}^{\mathrm{ERA},\mathrm{Gam}}$ and
$\bar{C}^{\mathrm{ORA},\mathrm{Gam}}$, (c) $P_{\rm out}^{\mathrm{ERA},\rm LN}$
and $P_{\rm out}^{\mathrm{ORA},\mathrm{LN}}$, and (d)
$\bar{C}^{\mathrm{ERA},\rm LN}$ and $\bar{C}^{\mathrm{ORA},\mathrm{LN}}$.

(a)

(b)

(c)

Figure 4: Comparison between the Gamma, Log-Normal, and Gaussian distribution
approximations in terms of a) PDF in the ERA scheme, b) the discrepancy in OP,
and c) the discrepancy in EC.

Furthermore, in Fig. 4, we discuss the performance behavior, and compare the
impact of the Gamma, Log-Normal, and Gaussian distributions on the OP and EC.
First, as can be seen in Fig. (4a), the Gamma and Log-Normal distributions
provide better approximations compared to the Gaussian distribution. Moreover,
from Fig. (4b), one can see that, compared to the Gamma and Log-Normal
distributions, the OP associated with the Gaussian distribution shows a
tangible gap between the approximate and true values at the low OP range,
i.e., when the OP is lower than $10^{-2}$. The reason behind this OP
discrepancy is the inaccurate approximation occurring at the left-hand tail of
the distributions. Fig. (4c) shows that the analytical and simulation results
of the EC obtained by the Gamma, Log-Normal, or Gaussian distributions are
well corroborated even in the high transmission power regime. This EC behavior
is in line with the results reported in [40]. In particular, the EC is not as
sensitive to the numerical inaccuracy as the OP. In short, the Gamma and Log-
Normal distribution approximations cause a small OP discrepancy in the large
SNR region; nevertheless, this discrepancy is negligible compared to the OP
discrepancy caused by the Gaussian distribution approximation. Moreover, the
discrepancy in the EC caused by the Gaussian, Gamma, or Log-Normal
distribution approximations is very small in the large SNR region and can be
negligible.

###  V-A The Difference in Lower Tail Between Gamma and Log-Normal
Distributions

In this subsection, we elaborate more on difference between the Gamma and Log-
Normal distributions. As shown in Fig. (4b), when we compare the OP computed
based on the Gamma and Log-Normal distributions, there is a very small
discrepancy, which can be negligible. Rigorously, this slight difference in OP
is due to the difference at the left tail (the lower tail) between the two
distributions. Specifically, let $\mathcal{N}(0,1)$ denote a standard Gaussian
(Normal) RV. Considering a Gamma RV $Y$ with the PDF in (11), mean
$\mathbb{E}[Y]=\alpha/\beta$ and variance $\mathrm{VAR}[Y]=\alpha/\beta^{2}$.
Note that $\beta$ is the rate parameter. Following the central-limit theorem
(CLT), when $\alpha$ is sufficiently large, $Y$ can be remodeled as

$\displaystyle
Y\sim\frac{\alpha}{\beta}\left[1+\frac{\mathcal{N}(0,1)}{\sqrt{\alpha}}\right],~{}\alpha\to\infty.$
(60)

Considering a Log-Normal RV $W$ with the PDF in (13), mean
$\mathbb{E}[W]=e^{\nu+\zeta^{2}/2}$ and variance
$\mathrm{VAR}[W]=(e^{\zeta^{2}}-1)e^{2\nu+\zeta^{2}}$, $W$ can be re-expressed
as $W=\exp(\nu+\zeta\mathcal{N}(0,1))$. Following the method in [41], to show
the relationship between the Gamma and Log-Normal distributions, we match
their means and variances as follows:

$\displaystyle\left\\{\begin{array}[]{l}\exp\left(\nu+\frac{\zeta^{2}}{2}\right)=\frac{\alpha}{\beta},\\\
\left(\exp\left(\zeta^{2}\right)-1\right)\exp\left(2\nu+\zeta^{2}\right)=\frac{\alpha}{\beta^{2}}.\end{array}\right.$
(63)

Solving (63) with respect to $\alpha$ and $\beta$ yields

$\displaystyle\zeta^{2}=\log\left(\frac{\alpha+1}{\alpha}\right),~{}\nu=\log\left(\frac{\alpha}{\beta}\frac{\sqrt{\alpha}}{\sqrt{\alpha+1}}\right).$
(64)

As $\alpha$ is sufficiently large, we have $\zeta\to\frac{1}{\sqrt{\alpha}}$
and $\nu\to\log\big{(}\frac{\alpha}{\beta}\big{)}$. Hence, $W$ can be
remodeled as
$W\sim\frac{\alpha}{\beta}\left[1+\frac{\mathcal{N}(0,1)}{\sqrt{\alpha}}\right],~{}\alpha\to\infty$.

This result is consistent with (60), which shows the similarity between the
Gamma and Log-Normal distributions. In Fig. 6, it can be observed that when
$\alpha$ is sufficiently large, e.g., $\alpha\geq 30$, the lower tails of the
Gamma and Log-Normal distributions are similar. Note that in our distributed
multi-RIS-aided system, the Nakagami shape parameter $m$ of each individual
channel is uniformly distributed in the range of $[2,3]$ as in Table I,
consequently, the corresponding shape parameter $\alpha_{Z}$ in (15) in the
ERA scheme is in the range of $[77,120]$.

###  V-B The Decomposition to The Sum of Squared Normal Random Variables

To give in-depth insight into the physical origin of the e2e channel
coefficient of the considered system, we can express the e2e channel
coefficients of either the ERA and ORA schemes as the sum of squared Gaussian
RVs as follows.

Let $U$, $V$, and $W$ be a Nakagami-$m$ RV, Gamma RV, and Gaussian RV,
respectively, specifically, $U\sim\mathrm{Nakagami}(m,\Omega)$,
$V\sim\mathrm{Gamma}(\alpha,\zeta)$, and $W\sim\mathcal{N}(0,\sigma^{2})$. On
one hand, knowing that if $U\sim\mathrm{Nakagami}(m,\Omega)$ and $U=\sqrt{V}$,
then $V\sim\mathrm{Gamma}(m,\Omega/m)$. On the other hand, the sum of $M$
squared Gaussian RVs, $\sum_{m=1}^{M}W_{m}^{2}$, follows a Gamma distribution,
i.e., $\sum_{m=1}^{M}W_{m}^{2}\sim\mathrm{Gamma}(M/2,2\sigma^{2})$. Based on
the above relationship between Nakagami-$m$, Gamma, and Gaussian RVs, assuming
that $m_{0}$, $m_{\mathrm{h}_{n}}$ and $m_{\mathrm{g}_{n}}$ are integers, we
can express the e2e channel coefficient of the ERA and ORA schemes as

$\displaystyle
Z=\sqrt{\sum_{i=1}^{m_{0}}\Big{(}X_{0,i}^{2}+Y_{0,i}^{2}\Big{)}}$
$\displaystyle\\!\\!+\\!\\!\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}\\!\kappa_{nl}\\!\sqrt{\sum_{j=1}^{m_{{\rm
h}_{n}}}\Big{(}X_{{\rm h}_{n},j}^{2}+Y_{{\rm
h}_{n},j}^{2}\Big{)}}\\!\sqrt{\sum_{k=1}^{m_{{\rm g}_{n}}}\Big{(}X_{{\rm
g}_{n},k}^{2}+Y_{{\rm g}_{n},k}^{2}\Big{)}},$ (65) $\displaystyle
R=\sqrt{\sum_{i=1}^{m_{0}}\\!\Big{(}X_{0,i}^{2}+Y_{0,i}^{2}\Big{)}}$
$\displaystyle\\!\\!+\\!\\!\max_{1\leq n\leq
N}\sum_{l=1}^{L_{n}}\kappa_{nl}\sqrt{\sum_{j=1}^{m_{{\rm
h}_{n}}}\Big{(}X_{{\rm h}_{n},j}^{2}\\!+\\!Y_{{\rm
h}_{n},j}^{2}\Big{)}}\sqrt{\sum_{k=1}^{m_{{\rm g}_{n}}}\\!\\!\Big{(}X_{{\rm
g}_{n},k}^{2}\\!+\\!Y_{{\rm g}_{n},k}^{2}\Big{)}},$ (66)

respectively, where $X_{{\rm c},i}$, $Y_{{\rm c},i}$, where ${\rm
c}\in\\{0,{\rm h}_{n},{\rm g}_{n}\\}$, are Gaussian RVs with zero means and
variances of $\frac{{\Omega_{\rm c}}}{2m_{\rm c}}$, i.e., $X_{{\rm
c},i},Y_{{\rm c},i}\sim\mathcal{N}\left(0,\frac{{\Omega_{\rm c}}}{2m_{\rm
c}}\right)$. In Figs. (6a) and (6b), we demonstrate the expressions in (65)
and (66), respectively. It is noted that the expressions in (65) and (66) are
valid when the shape parameters of the Nakagami-$m$ distributions are integer.

Figure 5: Illustrations of the Gamma PDF (dashed curves) and the Log-Normal
PDF (solid curves) with different values of $\alpha$.

Figure 6: Comparison between Gamma distribution and the distribution based on
the sum of squared Gaussian RVs in (a) the ERA scheme and (b) the ORA scheme.

In Fig. 7, to demonstrate the impact of i.n.i.d. fading channels on the system
performance, we consider different element settings, i.e.,
$\boldsymbol{\mathrm{L}}_{2}$, $\boldsymbol{\mathrm{L}}_{3}$, and
$\boldsymbol{\mathrm{L}}_{4}$. As can be observed in Fig. 7, the number of
passive elements installed at the RISs has a significant impact on both the OP
and EC. It is because under i.n.i.d. fading channels, each RIS is subject to
different fading severity, i.e., the Nakagami-$m$ distributions between
different RISs have different shape parameters and spread parameters. When
comparing the performance gaps between different element settings, the ORA
scheme is more sensitive than the ERA scheme to the changes of element
settings. For instance, as shown in Fig. 7a, when changing from element
setting $\boldsymbol{\mathrm{L}}_{4}$ to $\boldsymbol{\mathrm{L}}_{2}$, the
transmit power of the ERA scheme needs to compensate $1.8$ dBm, whereas the
ORA scheme needs to compensate for $6.3$ dBm to achieve the same OP of
$10^{-4}$. In addition, it is evidenced from Fig. 7 that RIS-aided systems are
superior to non-RIS-aided systems, and that the ERA scheme outperforms the ORA
scheme. The reason behind this is that the more of the reflecting elements
that are installed, the better the performance that can be achieved. For the
case of large-scale networks, as in [22], we consider that there are two
adjacent sources located at $(-99,0)$ and $(255,0)$, respectively, causing co-
channel inter-cell interference to the system. As can be observed in Fig. 7,
the interference decreases the performance of the considered multi-RIS-aided
system.

(a)

(b)

Figure 7: Numerical results demonstrating the impact of the number of
reflecting elements on (a) the OP and (b) the EC of the ERA and the ORA
schemes, respectively, with co-channel inter-cell interference and $R_{\rm
th}=3$ b/s/Hz.

Figure 8: Energy efficiency [Mbit/Joule] of the ERA and the ORA schemes as a
function of target SE, $R_{\rm th}$ [b/s/Hz].

Figure 9: Performance comparison between centralized versus distributed RIS-
aided systems, with $P_{\mathrm{S}}=23$ dBm.

In Fig. 9, taking into account the circuit dissipation power (CDP) in the
transceivers and the RIS hardware, we evaluate the EE of the ERA and the ORA
schemes. The total power consumed, $P_{\rm tol}$, by a RIS-aided system can be
expressed as $P_{\rm
tol}=P_{\mathrm{S}}+\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}\tilde{P}_{nl}+\tilde{P}_{\mathrm{S}}+\tilde{P}_{\mathrm{D}}$
[30], where $\tilde{P}_{nl}$, $\tilde{P}_{\mathrm{S}}$,
$\tilde{P}_{\mathrm{D}}$ are the CDP at the $l$-th element of the $n$-th RIS,
$\mathrm{S}$, and $\mathrm{D}$, respectively. The EE of a RIS-aided system can
be defined as $\mathrm{EE}=\mathrm{BW}\times R_{\rm th}/P_{\rm tol}$
[Mbit/Joule] [11], where $\mathrm{BW}$ and $R_{\rm th}$ denote the bandwidth
and the target SE. As shown in Fig. 9, for element setting
$\boldsymbol{\mathrm{L}}_{3}$, when $R_{\rm th}\leq\check{R}=12$ b/s/Hz, where
$\check{R}$ denotes the crossing point as illustrated in Fig. 9, the ORA
scheme achieves higher EE than the ERA scheme, and vice-versa when $R_{\rm
th}\geq\check{R}=12$ b/s/Hz. In addition, the element setting strongly impacts
the EE, i.e., when the element setting is $\boldsymbol{\mathrm{L}}_{2}$ and
$\boldsymbol{\mathrm{L}}_{4}$, the crossing points are $\check{R}=12.5$ and
$15.5$ b/s/Hz, respectively.

In Fig. 9, we compare the performance of the distributed multi-RIS-aided
scheme, i.e., the ERA scheme, versus a centralized RIS (C-RIS)-aided system.
Specifically, in the centralized-ERA scheme, we assume that all elements are
installed at one large C-RIS, and we let the C-RIS moves along a horizontal
line from $\mathrm{S}\to\mathrm{D}$ as depicted in Fig. 9, where the unfilled
red squares are the possible locations of the C-RIS. From various location and
reflecting element settings, we can observe a common behavior, i.e., as the
C-RIS is located near either $\mathrm{S}$ or $\mathrm{D}$, the centralized-ERA
scheme yields a larger EC than the distributed scheme. Meanwhile, when the
C-RIS moves further away from either $\mathrm{S}$ or $\mathrm{D}$, the ERA
scheme achieves a better EC.

###  V-C The Kullback-Leibler (KL) Divergence and The Kolmogorov-Smirnov (KS)
Test

To gain more insights into the accuracy of the distribution approximation, we
carry out the numerical analysis using the Kullback-Leibler (KL) divergence
and the Kolmogorov–Smirnov (KS) goodness-of-fit test. Let $f_{X}(x)$ and
$F_{X}(x)$ be the true PDF and CDF of an arbitrary random variable $X$. Let
$\tilde{f}_{X}(x)$ and $\tilde{F}_{X}(x)$ be the approximate PDF and CDF of
$X$, respectively.

The Kullback-Leibler (KL) Divergence: The KL divergence is a metric that
measures the difference between two probability distributions. In our case,
the KL divergence of $\tilde{f}_{X}(x)$ from $f_{X}(x)$, denoted
$\mathrm{D_{KL}}(f_{X}(x)\|\tilde{f}_{X}(x))$, can be expressed as [42]

$\displaystyle\mathrm{D_{KL}}(f_{X}(x)\|\tilde{f}_{X}(x))$
$\displaystyle=\int_{-\infty}^{\infty}f_{X}(x)\ln\frac{f_{X}(x)}{\tilde{f}_{X}(x)}dx.$
(67)

It is noted that the lower the KL divergence value, the better the accuracy of
the approximate distribution one can achieve.

The Kolmogorov–Smirnov (KS) Goodness-of-Fit Test: We consider a null
hypothesis $H_{0}$, which states that the true distribution of the e2e channel
coefficient can be approximated by either the Gamma or Log-Normal
distributions; otherwise, it is an alternative hypothesis $H_{1}$. The
Kolmogorov-Smirnov (KS) goodness-of-fit test relies on the KS distance,
$\mathrm{D_{KS}}$, which is the maximum value of the absolute difference
between the true and approximate distributions, and the critical value,
$\mathrm{D}_{s,p}$, to decide whether to accept or reject $H_{0}$.
Mathematically, the KS test can be described as [42, 43]

$\displaystyle\max_{s}\left|F_{X}(\boldsymbol{\mathrm{x}})-\tilde{F}_{X}(\boldsymbol{\mathrm{x}})\right|\triangleq\mathrm{D_{KS}}\underset{H_{0}}{\overset{H_{1}}{\gtrless}}\mathrm{D}_{s,p},$
(68)

where the critical value $\mathrm{D}_{s,p}$ is given by
$D_{s,p}\approx\sqrt{-(1/2s)\ln(p/2)}$ [43], where $p$ is the significance
level, $s$ is the number of random samples of $X$, and
$\boldsymbol{\mathrm{x}}$ is the $s\times 1$ random sample vector.

(a)

(b)

Figure 10: Demonstration of (a) the Kullback-Leibler divergence,
$\mathrm{D_{KL}}(f_{X}(x)\|\tilde{f}_{X}(x))$, and (b) the Kolmogorov-Smirnov
distance, $\mathrm{D_{KS}}$, as a function of the total number of passive
reflecting elements, $L=\sum_{n=1}^{N}L_{n}$.

In Fig. 10, we demonstrate the KL divergence and the KS test, in which, with
the help of the Statistics and Machine Learning Toolbox provided by MATLAB
[32], we estimate $f_{X}(x)$ and $F_{X}(x)$ using the ${\tt ksdensity}$ and
${\tt ecdf}$ functions, respectively. On one hand, in our KS test, we set
$p=0.05$, $s=10^{3}$ samples, and as a result, $D_{s,p}=0.043$. From Fig. 10,
when the total number of the reflecting elements is greater than $100$, the KS
test accepts $H_{0}$ at the significant level of $p=0.05$. On the other hand,
as shown in Fig. 10, when the total number of the reflecting elements is
small, the Gamma distribution approximation results in a lower
$\mathrm{D_{KL}}$ and $\mathrm{D_{KS}}$ than the Log-Normal distribution
approximation, which implies that the Gamma distribution provides a higher
accuracy than the Log-Normal distribution. However, when the total number of
the reflecting elements is sufficiently large, the Gamma and Log-Normal
distributions provide comparable approximation accuracies.

##  VI Conclusions

In this paper, we proposed two multi-RIS-aided schemes, namely the ERA and the
ORA schemes. We focused on the statistical characterization and modeling of
the two schemes. Toward this end, we proposed a mathematical framework to
determine the distribution of the e2e fading channel in both schemes.
Specifically, for the ERA scheme, the framework determined that the true
distribution of the magnitude of the e2e channel coefficient can be
approximated by either the Gamma or Log-Normal distributions. For the ORA
scheme, we obtained approximate closed-form expressions for the CDF and PDF of
the magnitude of the e2e channel. Moreover, the framework determined that the
ORA scheme’s fading channel can also be modeled by a Log-Normal distribution.
Based on the resultant fading models, we evaluated the performance of the two
schemes in terms of OP and EC. Numerical results showed that the ERA scheme
outperforms the ORA scheme in terms of OP and EC. Nevertheless, the ORA scheme
achieves better EE than the ERA scheme for some specific values of the target
SE. We showed that for a given total number of passive elements, the element
allocation and RISs’ position have a significant impact on the system
performance. On the other hand, the centralized large-RIS-aided system can
yield a better EC than the ERA scheme when the centralized RIS is located near
either the source or the destination; otherwise, the ERA scheme outperforms
the centralized RIS-aided system.

##  Appendix A Proof of Theorem 1

Since $h_{0}\sim\mathrm{Nakagami}(m_{0},\Omega_{0})$, the $k$-th moment of
$h_{0}$, $\mu_{h_{0}}(k)\triangleq\mathbb{E}[h_{0}^{k}]$, can be obtained as

$\displaystyle\mu_{h_{0}}(k)=\frac{\Gamma(m_{0}+k/2)}{\Gamma(m_{0})}\bigg{(}\frac{m_{0}}{\Omega_{0}}\bigg{)}^{-k/2}.$
(69)

We now turn our focus to the $k$-th moment of $U_{nl}$. The exact PDF of
$U_{nl}$ is obtained as

$\displaystyle f_{U_{nl}}(z)$
$\displaystyle=\frac{4\lambda_{nl}^{m_{\mathrm{h}_{n}}+m_{\mathrm{g}_{n}}}}{\Gamma(m_{\mathrm{h}_{n}})\Gamma(m_{\mathrm{g}_{n}})}x^{m_{\mathrm{h}_{n}}+m_{\mathrm{g}_{n}}-1}K_{m_{\mathrm{g}_{n}}-m_{\mathrm{h}_{n}}}(2\lambda_{nl}z),$
(70)

where
$\lambda_{nl}=\sqrt{\frac{1}{\kappa_{nl}^{2}}\frac{m_{\mathrm{h}_{n}}}{\Omega_{\mathrm{h}_{n}}}\frac{m_{\mathrm{g}_{n}}}{\Omega_{\mathrm{g}_{n}}}}$.
The detailed derivation of (70) can be briefly presented as follows.

The aim of this derivation is not only to present the exact PDF of $U_{nl}$
but also to point out that the magnitude of the synthesized channel of the
individual dual-hop channel with respect to a single reflecting element of a
RIS in our considered system follows a Generalized-$K$ ($K_{G}$) distribution.
Recall that, for a given $n$, $U_{nl}=\kappa_{nl}h_{nl}g_{nl}$, where
$U_{n1},\ldots,U_{nL_{n}}$ are i.i.d. RVs, but $h_{nl}$ and $g_{nl}$ are
i.n.i.d. RVs, $\forall n,\forall l$. Specifically, knowing that
$f_{XY}(z)=\int_{0}^{\infty}\frac{1}{x}f_{Y}\left(\frac{z}{x}\right)f_{X}(x)dx$,
the PDF of $U_{nl}$ can be written as

$f_{U_{nl}}(z)=\frac{1}{\kappa_{nl}}\int_{0}^{\infty}\frac{1}{x}f_{h_{nl}}\left(\frac{z}{\kappa_{nl}x}\right)f_{g_{nl}}(x)dx.$
(71)

Recall that
$h_{nl}\sim\mathrm{Nakagami}(m_{\mathrm{h}_{n}},\Omega_{\mathrm{h}_{n}})$ and
$g_{nl}\sim\mathrm{Nakagami}(m_{\mathrm{g}_{n}},\Omega_{\mathrm{g}_{n}})$,
from (9), we have

$\displaystyle
f_{U_{nl}}(z)=\frac{4}{\Gamma(m_{\mathrm{h}_{n}})\Gamma(m_{\mathrm{g}_{n}})}\left(\frac{1}{\kappa_{nl}^{2}}\\!\\!\\!\\!\\!\frac{m_{\mathrm{h}_{n}}}{\Omega_{\mathrm{h}_{n}}}\right)^{m_{\mathrm{h}_{n}}}\\!\\!\\!\left(\frac{m_{\mathrm{g}_{n}}}{\Omega_{\mathrm{g}_{n}}}\right)^{m_{\mathrm{g}_{n}}}\\!\\!\\!\\!z^{2m_{\mathrm{h}_{n}}-1}$
$\displaystyle\quad\times\int_{0}^{\infty}x^{2m_{\mathrm{g}_{n}}-2m_{\mathrm{h}_{n}}-1}e^{\left(-\frac{z^{2}}{\kappa_{nl}^{2}}\frac{m_{\mathrm{h}_{n}}}{\Omega_{\mathrm{h}_{n}}}\frac{1}{x^{2}}-\frac{m_{\mathrm{g}_{n}}}{\Omega_{\mathrm{g}_{n}}}x^{2}\right)}dx.$
(72)

Making use of [23, Eq. (3.478.4)], the exact PDF of $U_{nl}$ is obtained as in
(70). Thus, we can conclude that $U_{nl}$ follows a $K_{G}$ distribution with
the shaping parameters $m_{\mathrm{g}_{n}}$ and $m_{\mathrm{h}_{n}}$. Based on
the exact PDF of $U_{nl}$, we are going to derive the $k$-th moment of
$U_{nl}$, which is defined as
$\mu_{U_{nl}}(k)\triangleq\mathbb{E}[U_{nl}^{k}]=\int_{0}^{\infty}z^{k}f_{U_{nl}}(z)dz$.
From (70) and after some mathematical manipulations, $\mu_{U_{nl}}(k)$ can be
obtained as

$\displaystyle\mu_{U_{nl}}(k)=\lambda_{nl}^{-k}\frac{\Gamma(m_{\mathrm{h}_{n}}+k/2)\Gamma(m_{\mathrm{g}_{n}}+k/2)}{\Gamma(m_{\mathrm{h}_{n}})\Gamma(m_{\mathrm{g}_{n}})}.$
(73)

The exact PDF of $U_{nl}$ in (70) helps derive the $k$-th moment of $U_{nl}$,
however, it makes the derivation of closed-form expressions for the CDF and
PDF of $Z$ intractable. To circumvent this problem, based on the obtained
$k$-th moment of $U_{nl}$ in (73), we fit $U_{nl}$ to a $\rm Gamma$
distribution. Note that for a given $n$, the $U_{nl}$ RVs are i.i.d. with
respect to $l=1,\ldots,L_{n}$. Relying on Step 2 in our framework, we have

$\displaystyle U_{nl}\overset{\rm
approx.}{\sim}\mathrm{Gamma}(\alpha_{U_{n}},\beta_{U_{n}}),$ (74)

where

$\displaystyle\alpha_{U_{n}}=\frac{\left(\mathbb{E}[U_{nl}]\right)^{2}}{\mathrm{Var}[U_{nl}]}=\frac{[\mu_{U_{nl}}(1)]^{2}}{\mu_{U_{nl}}(2)-[\mu_{U_{nl}}(1)]^{2}},$
(75)
$\displaystyle\beta_{U_{n}}=\frac{\mathbb{E}[U_{nl}]}{\mathrm{Var}[U_{nl}]}=\frac{\mu_{U_{nl}}(1)}{\mu_{U_{nl}}(2)-[\mu_{U_{nl}}(1)]^{2}},$
(76)

and $f_{U_{nl}}(z;\alpha_{U_{n}},\beta_{U_{n}})$ is expressed as in (11). The
analysis will rely on the approximate PDF of $U_{nl}$. In addition, for a
given $n$-th RIS, we consider fixed amplitude reflection coefficients [27],
i.e., $\kappa_{nl}=\kappa_{n},\forall l$.

The approximate distribution of $V_{n}$ can be obtained as

$\displaystyle V_{n}\overset{\rm
approx.}{\sim}\mathrm{Gamma}(L_{n}\alpha_{U_{n}},\beta_{U_{n}}).$ (77)

Consequently, the approximate CDF and PDF of $V_{n}$ are obtained as

$\displaystyle F_{V_{n}}(z)$
$\displaystyle\approx\frac{\gamma(L_{n}\alpha_{U_{n}},\beta_{U_{n}}z)}{\Gamma\left(L_{n}\alpha_{U_{n}}\right)},$
(78) $\displaystyle f_{V_{n}}(z)$ $\displaystyle\approx
f_{V_{n}}(z;L_{n}\alpha_{U_{n}},\beta_{U_{n}}),$ (79)

where $f_{V_{n}}(z;L_{n}\alpha_{U_{n}},\beta_{U_{n}})$ is expressed as in
(11).

Making use of a multinomial expansion [44], the $k$-th moment of $V_{n}$,
i.e., $\mu_{V_{n}}(k)\triangleq\mathbb{E}[V_{n}^{k}]$, can be obtained as

$\displaystyle\mu_{V_{n}}(k)$
$\displaystyle=\sum_{k_{1}=0}^{k}\sum_{k_{2}=0}^{k_{1}}\cdots\sum_{k_{L_{n}-1}=0}^{k_{L_{n}-2}}\binom{k}{k_{1}}\binom{k_{1}}{k_{2}}\cdots\binom{k_{L_{n}-2}}{k_{L_{n}-1}}$
$\displaystyle\quad\times\mu_{U_{n1}}(k-k_{1})\mu_{U_{n2}}(k_{1}-k_{2})\ldots\mu_{U_{nL_{n}}}(k_{L_{n}-1}).$
(80)

We turn our focus to the $k$-th moment of $T$, i.e.,
$\mu_{T}(k)\triangleq\mathbb{E}[T^{k}]$. Proceeding in a similar manner, the
$k$-th moment of $T$, i.e., $\mu_{T}(k)\triangleq\mathbb{E}[T^{k}]$, can be
obtained as

$\displaystyle\mu_{T}(k)$
$\displaystyle=\sum_{k_{1}=0}^{k}\sum_{k_{2}=0}^{k_{1}}\cdots\sum_{k_{N-1}=0}^{k_{N-2}}\binom{k}{k_{1}}\binom{k_{1}}{k_{2}}\cdots\binom{k_{N-2}}{k_{N-1}}$
$\displaystyle\quad\times\mu_{V_{1}}(k-k_{1})\mu_{V_{2}}(k_{1}-k_{2})\ldots\mu_{V_{n}}(k_{N-1}).$
(81)

From (73), (A), and (A), the first and second moments of $T$ can be expressed
as

$\displaystyle\mu_{T}(1)$
$\displaystyle=\sum_{n=1}^{N}\sum_{l=1}^{L_{n}}\mu_{U_{nl}}(1),$ (82)
$\displaystyle\mu_{T}(2)$
$\displaystyle=\sum_{n=1}^{N}\bigg{[}\sum_{l=1}^{L_{n}}\mu_{U_{nl}}(2)+2\sum_{l=1}^{L_{n}}\mu_{U_{nl}}(1)\sum_{l^{\prime}=l+1}^{L_{n}}\mu_{U_{nl^{\prime}}}(1)\bigg{]}$
$\displaystyle\quad+2\sum_{n=1}^{N}\bigg{[}\sum_{l=1}^{L_{n}}\mu_{U_{nl}}(1)\bigg{]}\sum_{n^{\prime}=n+1}^{N}\bigg{[}\sum_{l=1}^{L_{n^{\prime}}}\mu_{U_{n^{\prime}l}}(1)\bigg{]}.$
(83)

Since $h_{0}$ and $T$ are independent, the $k$-th moment of $Z$,
$\mu_{Z}(k)\triangleq\mathbb{E}[Z^{k}]$, can be obtained via the moments of
its summands, i.e., $h_{0}$ and $T$, by applying the binomial theorem. Thus,
$\mu_{Z}(k)$ can be obtained as

$\displaystyle\mu_{Z}(k)$
$\displaystyle=\mathbb{E}[(h_{0}+T)^{k}]=\mathbb{E}\bigg{[}\sum_{l=0}^{k}\binom{k}{l}h_{0}^{l}T^{k-l}\bigg{]}$
$\displaystyle=\sum_{l=0}^{k}\binom{k}{l}\mu_{h_{0}}(l)\mu_{T}(k-l).$ (84)

Thus, from (A), we have

$\displaystyle\mu_{Z}(1)$ $\displaystyle=\mu_{h_{0}}(1)+\mu_{T}(1),$ (85)
$\displaystyle\mu_{Z}(2)$
$\displaystyle=\mu_{h_{0}}(2)+\mu_{T}(2)+2\mu_{h_{0}}(1)\mu_{T}(1).$ (86)

From (69), we have

$\displaystyle\mu_{h_{0}}(1)$
$\displaystyle=\frac{\Gamma(m_{0}+1/2)}{\Gamma(m_{0})}\bigg{(}\frac{m_{0}}{\Omega_{0}}\bigg{)}^{-1/2},$
(87) $\displaystyle\mu_{h_{0}}(2)$
$\displaystyle=\frac{\Gamma(m_{0}+1)}{\Gamma(m_{0})}\frac{\Omega_{0}}{m_{0}}=\Omega_{0}.$
(88)

Plugging (87) and (82) into (85), and plugging (88) and (83) into (86), we
complete the proof of Theorem 1.

##  Appendix B Derivation of (25)

From (25), the EC of the ERA scheme can be rewritten as

$\displaystyle\bar{C}^{\mathrm{ERA},\rm Gam}$
$\displaystyle=\int_{0}^{\infty}\log_{2}\left(1+z\right)f_{\bar{\rho}Z^{2}}(z)dz$
$\displaystyle=\frac{1}{\ln
2}\int_{0}^{\infty}\frac{1}{z+1}\left[1-F_{Z}\left(\sqrt{\frac{z}{\bar{\rho}}}\right)\right]dz.$
(89)

Making use of the identity [25, Eq. (8.4.2.5)], i.e.,
$(1+x)^{-\xi}=\frac{1}{\Gamma(\xi)}G_{1,1}^{1,1}\left(x\middle|\begin{matrix}1-\xi\\\
0\end{matrix}\right)$, (89) can be further expressed as

$\displaystyle\bar{C}^{\mathrm{ERA},\rm Gam}\approx\frac{1}{\ln
2}\int_{0}^{\infty}G_{1,1}^{1,1}\left(z\middle|\begin{matrix}0\\\
0\end{matrix}\right)\frac{\Gamma\left({\alpha_{Z}},{\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\right)}{\Gamma({\alpha_{Z}})}dz.$
(90)

Based on [23, Eq. (8.356.3)], i.e.,
$\Gamma(\alpha,x)+\gamma(\alpha,x)=\Gamma(\alpha)$, we have

$\displaystyle
F_{Z}(z)=\frac{\gamma\left({\alpha_{Z}},{\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\right)}{\Gamma({\alpha_{Z}})}=1-\frac{\Gamma\left({\alpha_{Z}},{\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\right)}{\Gamma({\alpha_{Z}})}.$
(91)

Using the equivalence between the incomplete Gamma function and the Meijer-G
function [25, Eq. (8.4.16.2)], i.e.,
$\Gamma\left({\alpha_{Z}},{\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\right)=G^{2,0}_{1,2}\left({\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\middle|\begin{matrix}1\\\
{\alpha_{Z}},0\end{matrix}\right)$, (90) can be further expressed as

$\displaystyle\bar{C}^{\mathrm{ERA},\rm Gam}$
$\displaystyle\approx\frac{1}{\Gamma({\alpha_{Z}})\ln
2}\int_{0}^{\infty}G_{1,1}^{1,1}\left(z\middle|\begin{matrix}0\\\
0\end{matrix}\right)G^{2,0}_{1,2}\left({\beta_{Z}}\sqrt{\frac{z}{\bar{\rho}}}\middle|\begin{matrix}1\\\
{\alpha_{Z}},0\end{matrix}\right)dz.$ (92)

Next, we rely on the identity [25, Eq. (2.24.1.1)] for the general integral of

$\displaystyle\int_{0}^{\infty}x^{\alpha-1}G^{s,t}_{u,v}\left(\xi
x\middle|\begin{matrix}(c_{u})\\\
(d_{v})\end{matrix}\right)G^{m,n}_{p,q}\left(\omega
x^{l/k}\middle|\begin{matrix}(a_{p})\\\ (b_{q})\end{matrix}\right)dx.$

Thus, the EC of the ERA scheme can be obtained as

$\displaystyle\bar{C}^{\mathrm{ERA},\rm Gam}$
$\displaystyle\approx\frac{1}{\Gamma({\alpha_{Z}})\ln
2}\frac{2^{{\alpha_{Z}}-1}}{\sqrt{\pi}}G^{5,1}_{3,5}\left(\frac{({\beta_{Z}})^{2}}{4\bar{\rho}}\middle|\begin{matrix}0,\frac{1}{2},1\\\
\frac{{\alpha_{Z}}}{2},\frac{{\alpha_{Z}}+1}{2},0,\frac{1}{2},0\end{matrix}\right).$
(93)

##  Appendix C Proof of Lemma 2

The $k$-moment of $M_{V}$, $\mu_{M_{V}}(k)\triangleq\mathbb{E}[M_{V}^{k}]$,
can be expressed as

$\displaystyle\mu_{M_{V}}(k)=\sum_{n=1}^{N}\int_{0}^{\infty}x^{k}f_{V_{n}}(x)\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}F_{V_{t}}(x)dx.$ (94)

Thus, we have

$\displaystyle\mu_{M_{V}}(k)\approx\sum_{n=1}^{N}\int_{0}^{\infty}x^{k}f_{V_{n}}(x)\prod_{\begin{subarray}{c}t=1\\\
t\neq
n\end{subarray}}^{N}\frac{\gamma(L_{t}{\alpha_{U}}_{t},{\beta_{U}}_{t})}{\Gamma(L_{t}{\alpha_{U}}_{t})}dx.$
(95)

The $\gamma(\cdot,\cdot)$ in (78) can be re-expressed as [45, Eq. (6.5.29)]

$\displaystyle\frac{\gamma(L_{t}\alpha_{U_{n}},{\beta_{U}}_{t}x)}{\Gamma(L_{t}{\alpha_{U}}_{t})}=e^{-{\beta_{U}}_{t}x}\sum_{w_{t}=0}^{\infty}\frac{({\beta_{U}}_{t}x)^{L_{t}{\alpha_{U}}_{t}+w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}.$
(96)

Plugging (96) into $\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}F_{{V}_{t}}(x)$, after some mathematical manipulations, we
have

$\displaystyle\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{\gamma(L_{t}{\alpha_{U}}_{t},{\beta_{U}}_{t})}{\Gamma(L_{t}{\alpha_{U}}_{t})}$
$\displaystyle\quad=e^{-x\sum_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}{\beta_{U}}_{t}}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t}x)^{L_{t}{\alpha_{U}}_{t}+w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]},$
(97)

where

$\displaystyle\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\ t\neq
n\end{subarray}}}\triangleq\sum_{w_{1}=0}^{\infty}\ldots\sum_{w_{n-1}=0}^{\infty}\sum_{w_{n+1}=0}^{\infty}\ldots\sum_{w_{N}=0}^{\infty}.$

Thus, we have

$\displaystyle\mu_{M_{V}}(k)\approx\sum_{n=1}^{N}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t}x)^{L_{t}{\alpha_{U}}_{t}+w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}$
$\displaystyle\times\int_{0}^{\infty}\frac{(\beta_{U_{n}})^{L_{n}\alpha_{U_{n}}}}{\Gamma(L_{n}\alpha_{U_{n}})}x^{\sum_{j=1}^{N}L_{j}{\alpha_{U}}_{j}+\sum_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}w_{t}+k-1}e^{-x\sum_{i=1}^{N}{\beta_{U}}_{i}}dx.$ (98)

Let us denote $\Psi_{\check{n}}\triangleq\sum_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}w_{t}$ and
$\Lambda\triangleq\sum_{t=1}^{N}L_{t}{\alpha_{U}}_{t}$. After some
mathematical manipulations, $\mu_{M_{V}}(k)$ can be expressed as

$\displaystyle\mu_{M_{V}}(k)$
$\displaystyle\approx\sum_{n=1}^{N}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t}x)^{L_{t}{\alpha_{U}}_{t}+w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}\frac{(\beta_{U_{n}})^{L_{n}\alpha_{U_{n}}}}{\Gamma(L_{n}\alpha_{U_{n}})}$
$\displaystyle\quad\times\int_{0}^{\infty}x^{\Lambda+\Psi_{\check{n}}+k-1}e^{-x\sum_{i=1}^{N}{\beta_{U}}_{i}}dx.$
(99)

Making use of the identity $\int_{0}^{\infty}x^{v-1}e^{-px}dx=\Gamma(v)/p^{v}$
[23, Eq. (3.381.4)], the integral in (99) can be derived as

$\displaystyle\mu_{M_{V}}(k)$
$\displaystyle\approx\sum_{n=1}^{N}\bigg{\\{}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t}x)^{L_{t}{\alpha_{U}}_{t}+w_{t}}}{}$
$\displaystyle\times{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}\frac{(\beta_{U_{n}})^{L_{n}\alpha_{U_{n}}}}{\Gamma(L_{n}\alpha_{U_{n}})}\Gamma(\Lambda+\Psi_{\check{n}}+k)$
$\displaystyle\times\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-(\Lambda+\Psi_{\check{n}}+k)}\bigg{\\}}.$
(100)

After some rearrangements, we have

$\displaystyle\mu_{M_{V}}(k)$
$\displaystyle\approx\sum_{n=1}^{N}\bigg{\\{}\frac{1}{\Gamma(L_{n}{\alpha_{U}}_{t})}\bigg{[}(\beta_{U_{n}})^{L_{n}\alpha_{U_{n}}}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}({\beta_{U}}_{t})^{L_{t}{\alpha_{U}}_{t}}\bigg{]}$
$\displaystyle\times\bigg{[}\sum_{t=i}^{N}{\beta_{U}}_{i}\bigg{]}^{-(\Lambda+k)}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\Gamma(\Lambda+\Psi_{\check{n}}+k)$
$\displaystyle\times\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t})^{w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-\Psi_{\check{n}}}\bigg{\\}}.$
(101)

Let $\chi_{t}={\beta_{U}}_{t}/\sum_{v=1}^{N}{\beta_{U}}_{v}$. It can be
observed that $\chi_{t}\in(0,1)$ and $\sum_{t=1}^{N}\chi_{t}=1$, which yields

$\displaystyle\bigg{[}(\beta_{U_{n}})^{L_{n}\alpha_{U_{n}}}\prod_{\begin{subarray}{c}t=1\\\
t\neq
n\end{subarray}}^{N}{\beta_{U}}_{t}^{L_{t}{\alpha_{U}}_{t}}\bigg{]}\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{(-\Lambda)}$
$\displaystyle=\prod_{j=1}^{N}(\chi_{j})^{L_{j}{\alpha_{U}}_{j}},$ (102)

By plugging (102) into (101), we have

$\displaystyle\mu_{M_{V}}(k)$
$\displaystyle\approx\sum_{n=1}^{N}\bigg{\\{}\frac{1}{\Gamma(L_{n}{\alpha_{U}}_{t})}\bigg{[}\prod_{j=1}^{N}(\chi_{j})^{L_{j}{\alpha_{U}}_{j}}\bigg{]}$
$\displaystyle\times\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-k}\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\
t\neq n\end{subarray}}}\Gamma(\Lambda+\Psi_{\check{n}}+k)$
$\displaystyle\times\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t})^{w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-\Psi_{\check{n}}}\bigg{\\}}.$
(103)

After some mathematical manipulations, it follows that

$\displaystyle\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{({\beta_{U}}_{t})^{w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\bigg{]}\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-\Psi_{\check{n}}}$
$\displaystyle\quad=\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}\frac{1}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}(\chi_{t})^{w_{t}}.$ (104)

Note that, in (104), $i$ is independent from $t$. Next, let $\langle
x\rangle_{n}\triangleq\Gamma(x+n)/\Gamma(x)$ denote a Pochhammer symbol, which
implies

$\displaystyle\Gamma(\Lambda+k+\Psi_{\check{n}})=\langle\Lambda+k\rangle_{\Psi_{\check{n}}}\Gamma(\Lambda+k).$
(105)

Thus, the right-hand-side (RHS) of (104) can be expressed as

$\displaystyle\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}\frac{(\chi_{t})^{w_{t}}}{\Gamma(L_{t}{\alpha_{U}}_{t}+w_{t}+1)}$
$\displaystyle=\bigg{[}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}\frac{1}{\Gamma(L_{t}{\alpha_{U}}_{t}+1)}\bigg{]}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}\frac{\langle 1\rangle_{w_{t}}}{\langle
L_{t}{\alpha_{U}}_{t}+1\rangle_{w_{t}}}\bigg{]}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}\frac{(\chi_{t})^{w_{t}}}{w_{t}!}\bigg{]}.$ (106)

After some mathematical manipulations, we have

$\displaystyle\mu_{M_{V}}(k)\approx\bigg{[}\sum_{i=1}^{N}{\beta_{U}}_{i}\bigg{]}^{-k}\bigg{[}\prod_{j=1}^{N}(\chi_{j})^{L_{j}{\alpha_{U}}_{j}}\bigg{]}$
$\displaystyle\times\sum_{n=1}^{N}\bigg{\\{}\frac{\Gamma(\Lambda+k)}{\Gamma(L_{n}\alpha_{U_{n}})}\bigg{[}\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}\frac{1}{\Gamma(L_{n}\alpha_{U_{n}}+1)}\bigg{]}$
$\displaystyle\times\widetilde{\sum_{\begin{subarray}{c}w_{t}\\\ t\neq
n\end{subarray}}}\bigg{[}\frac{\langle\Lambda+k\rangle_{\sum_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}w_{t}}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\langle 1\rangle_{w_{t}}}{\prod_{\begin{subarray}{c}t=1\\\
t\neq n\end{subarray}}^{N}\langle
L_{t}{\alpha_{U}}_{t}+1\rangle_{w_{t}}}\prod_{\begin{subarray}{c}t=1\\\ t\neq
n\end{subarray}}^{N}\frac{(\chi_{t})^{w_{t}}}{w_{t}!}\bigg{]}\bigg{\\}}.$
(107)

Using the series representation of the Lauricell function Type-A [26,
Eq.(1.4.1)], and relying on its integral representation [26, Eq.(9.4.35)], we
obtain the approximate closed-form expression for $\mu_{M_{V}}(k)$ in (57).
This completes the proof of Lemma 2.

## References

  * [1] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, Nov. 2020\.
  * [2] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [3] C. Liaskos, S. Nie, A. Tsioliaridou, A. Pitsillides, S. Ioannidis, and I. Akyildiz, “End-to-end wireless path deployment with intelligent surfaces using interpretable neural networks,” _IEEE Trans. Commun._ , vol. 68, no. 11, pp. 6792–6806, Nov. 2020.
  * [4] D. L. Galappaththige, D. Kudathanthirige, and G. A. A. Baduge, “Performance analysis of distributed intelligent reflective surfaces for wireless communications,” Oct. 2020. [Online]. Available: http://arxiv.org/abs/2010.12543
  * [5] B. Tahir, S. Schwarz, and M. Rupp, “Analysis of uplink IRS-assisted NOMA under Nakagami-m fading via moments matching,” _IEEE Wireless Commun. Lett._ , 2020, DOI: 10.1109/LWC.2020.3043810.
  * [6] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, Jun. 2019.
  * [7] F. A. P. de Figueiredo, M. S. P. Facina, R. C. Ferreira, Y. Ai, R. Ruby, Q.-V. Pham, and G. Fraidenraich, “Large intelligent surfaces with discrete set of phase-shifts communicating through double-rayleigh fading channels,” _IEEE Access_ , vol. 9, pp. 20 768–20 787, Feb. 2021.
  * [8] M.-A. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 184–188, Feb. 2020.
  * [9] A.-A. A. Boulogeorgos and A. Alexiou, “Performance analysis of reconfigurable intelligent surface-assisted wireless systems and comparison with relaying,” _IEEE Access_ , vol. 8, pp. 94 463–94 483, Jun. 2020.
  * [10] X. Gan, C. Zhong, Y. Zhu, and Z. Zhong, “User selection in reconfigurable intelligent surface assisted communication systems,” _IEEE Commun. Lett._ , 2021, DOI: 10.1109/LCOMM.2020.3048782.
  * [11] E. Bjornson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , 2020, DOI: 10.1109/LWC.2020.3046107.
  * [12] T. Van Chien, A. K. Papazafeiropoulos, L. T. Tu, R. Chopra, S. Chatzinotas, and B. Ottersten, “Outage Probability Analysis of IRS-Assisted Systems Under Spatially Correlated Channels,” _IEEE Wireless Communications Letters_ , Feb. 2021.
  * [13] H. Ibrahim, H. Tabassum, and U. T. Nguyen, “Exact coverage analysis of intelligent reflecting surfaces with nakagami-m channels,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 1, pp. 1072–1076, Jan. 2021.
  * [14] Z. Cui, K. Guan, J. Zhang, and Z. Zhong, “SNR coverage probability analysis of RIS-aided communication systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 4, pp. 3914–3919, Apr. 2021.
  * [15] M. Jung, W. Saad, Y. Jang, G. Kong, and S. Choi, “Performance analysis of large intelligent surfaces (LISs): Asymptotic data rate and channel hardening effects,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 3, pp. 2052–2065, Mar. 2020.
  * [16] L. Yang, Y. Yang, D. B. da Costa, and I. Trigui, “Outage probability and capacity scaling law of multiple RIS-aided networks,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 2, pp. 256–260, Feb. 2021.
  * [17] Y. Fang, S. Atapattu, H. Inaltekin, and J. Evans, “Optimum reconfigurable intelligent surface selection for indoor and outdoor communications,” Dec. 2020\. [Online]. Available: http://arxiv.org/abs/2012.11793
  * [18] W. Mei and R. Zhang, “Cooperative beam routing for multi-IRS aided communication,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 2, pp. 426–430, Feb. 2021.
  * [19] S. Zhang and R. Zhang, “Intelligent reflecting surface aided multiple access: Capacity region and deployment strategy,” Feb. 2020. [Online]. Available: http://arxiv.org/abs/2002.07091
  * [20] I. Yildirim, A. Uyrus, and E. Basar, “Modeling and analysis of reconfigurable intelligent surfaces for indoor and outdoor applications in future wireless networks,” _IEEE Trans. Commun._ , 2020, DOI: 10.1109/TCOMM.2020.3035391.
  * [21] J. Lyu and R. Zhang, “Spatial throughput characterization for intelligent reflecting surface aided multiuser system,” _IEEE Wireless Communications Letters_ , vol. 9, no. 6, pp. 834–838, Jun. 2020.
  * [22] ——, “Hybrid active/passive wireless network aided by intelligent reflecting surface: System modeling and performance analysis,” _IEEE Trans. Wirel. Commun._ , 2021, DOI:10.1109/TWC.2021.3081447.
  * [23] I. S. Gradshteyn and I. M. Ryzhik, _Tables of Integrals, Series, and Products_ , 7th ed. New York, NY, USA: Academic Press, 2007.
  * [24] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_. New York, NY, USA: Dover, 1965.
  * [25] A. Prudnikov, Y. Brychkov, and O. Marichev, _Integrals and Series, Volume 3: More Special Functions_. Boca Raton, FL, USA: CRC Press, 1999.
  * [26] H. M. Srivastava and P. W. Karlsson, _Multiple Gaussian Hypergeometric Series_. Hoboken, NJ, USA: Wiley, 1985\.
  * [27] E. Bjornson, O. Ozdogan, and E. G. Larsson, “Intelligent reflecting surface versus decode-and-forward: How large surfaces are needed to beat relaying?” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 244–248, Feb. 2020.
  * [28] B. Zheng and R. Zhang, “Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization,” _IEEE Wireless Communications Letters_ , vol. 9, no. 4, pp. 518–522, 2020.
  * [29] H. Alwazani, Q.-U.-A. Nadeem, and A. Chaaban, “Channel estimation for distributed intelligent reflecting surfaces assisted multi-user MISO systems,” in _2020 IEEE Globecom Workshops_. IEEE, Dec. 2020, pp. 1–6.
  * [30] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wirel. Commun._ , vol. 18, no. 8, pp. 4157–4170, Aug. 2019.
  * [31] L. Yang, F. Meng, J. Zhang, M. O. Hasna, and M. D. Renzo, “On the performance of RIS-assisted dual-hop UAV communication systems,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 9, pp. 10 385–10 390, Sep. 2020.
  * [32] _MATLAB version R2021a_ , The Mathworks, Inc., Natick, MA, USA, 2021.
  * [33] K. O. Bowman and L. R. Shenton, _Estimator: Method of Moments_. Hoboken, NJ, USA: Wiley, 1998.
  * [34] P. Z. Peebles, _Probability, Random Variables and Random Signal Principles_ , 4th ed. New York, NY, USA: McGraw-Hill Science, 2000.
  * [35] P. S. Bithas, V. Nikolaidis, A. G. Kanatas, and G. K. Karagiannidis, “UAV-to-ground communications: Channel modeling and UAV selection,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 5135–5144, Aug. 2020.
  * [36] E. W. Stacy, “A generalization of the gamma distribution,” _The Annals of Mathematical Statistics_ , vol. 33, no. 3, pp. 1187–1192, Sep. 1962.
  * [37] L. Fenton, “The sum of log-normal probability distributions in scatter transmission systems,” _IEEE Trans. Commun._ , vol. 8, no. 1, pp. 57–67, 1960.
  * [38] A. Laourine, A. Stéphenne, and S. Affes, “Estimating the ergodic capacity of Log-Normal channels,” _IEEE Commun. Lett._ , vol. 11, no. 7, pp. 568–570, Jul. 2007.
  * [39] “Evolved universal terrestrial radio access (E-UTRA); further advancements for e-utra physical layer aspects (release 9),” _Document 3GPP-TR-36.814 V9.0.0_ , Mar. 2010.
  * [40] S. Al-Ahmadi and H. Yanikomeroglu, “On the approximation of the generalized-K distribution by a Gamma distribution for modeling composite fading channels,” _IEEE Trans. Wirel. Commun._ , vol. 9, no. 2, pp. 706–713, Feb. 2010.
  * [41] I. Kosti, “Analytical approach to performance analysis for channel subject to shadowing and fading,” _IEE Proceedings - Communications_ , vol. 152, no. 6, p. 821, 2005.
  * [42] K. Peppas, “Accurate closed-form approximations to generalised-K sum distributions and applications in the performance analysis of equal-gain combining receivers,” _IET Commun._ , vol. 5, no. 7, pp. 982–989, May 2011\.
  * [43] L. Kong, Y. Ai, S. Chatzinotas, and B. Ottersten, “Effective Rate Evaluation of RIS-Assisted Communications Using the Sums of Cascaded $\alpha$-$\mu$ Random Variates,” _IEEE Access_ , vol. 9, pp. 5832–5844, 2021.
  * [44] D. B. da Costa, H. Ding, and J. Ge, “Interference-limited relaying transmissions in dual-pop cooperative networks over Nakagami-m fading,” _IEEE Commun. Lett._ , vol. 15, no. 5, pp. 503–505, May 2011.
  * [45] M. Abramowitz and I. A. Stegun, _Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables_. New York, NY, USA: Dover Books, 1972.