# Performance Analysis of Distributed Intelligent Reflective Surfaces for
Wireless Communications

Diluka Loku Galappaththige, , Dhanushka Kudathanthirige, , and Gayan
Amarasuriya Aruma Baduge The authors are with the School of Electrical,
Computer, and Biomedical Engineering, Southern Illinois University,
Carbondale, IL, USA, Email:
<EMAIL_ADDRESS>This work in part
has been submitted to IEEE Global Communication Conference (Globecom), Dec.,
2020 [1].

###### Abstract

In this paper, a comprehensive performance analysis of a distributed
intelligent reflective surfaces (IRSs)-aided communication system is
presented. First, the optimal signal-to-noise ratio (SNR), which is attainable
through the direct and reflected channels, is quantified by controlling the
phase-shifts of the distributed IRSs. Next, this optimal SNR is statistically
characterized by deriving tight approximations to the exact probability
density function (PDF) and cumulative distribution function (CDF) for
Nakagami-$m$ fading. The accuracy/tightness of this statistical
characterization is investigated by deriving the Kullback-Leibler divergence.
Our PDF/CDF analysis is used to derive tight approximations/bounds for the
outage probability, achievable rate, and average symbol error rate (SER) in
closed-form. To obtain useful insights, the asymptotic outage probability and
average SER are derived for the high SNR regime. Thereby, the achievable
diversity order and array gains are quantified. Our asymptotic performance
analysis reveals that the diversity order can be boosted by using distributed
passive IRSs without generating additional electromagnetic (EM) waves via
active radio frequency chains. Our asymptotic rate analysis shows that the
lower and upper rate bounds converge to an asymptotic limit in large
reflective element regime. Our analysis is validated via Monte-Carlo
simulations. We present a rigorous set of numerical results to investigate the
performance gains of the proposed system model. Our analytical and numerical
results reveal that the performance of single-input single-output wireless
systems can be boosted by recycling the EM waves generated by a transmitter
through distributed passive IRS reflections to enable constructive signal
combining at a receiver.

##  I Introduction

Recently, a novel concept of coating physical objects such as building wall
and windows within a wireless propagation medium with intelligent reflective
surfaces (IRSs), which can interact with electromagnetic (EM) waves, has
emerged [1, 2, 3]. By invoking the recent research developments of meta-
surfaces [4], the feasibility of synthesizing the IRS has been explored [2].
An IRS consists of a large number of passive reflective elements, which can
introduce controllable phase-shifts to incident EM waves. These phase-shifts
can be intelligently controlled to enable a constructive addition of EM waves
at a desired destination and thereby boosting the signal-to-noise ratio (SNR)
of the end-to-end communication. Moreover, the IRS enables recycling of EM
waves without generating additional signals via radio-frequency (RF)
chains/amplifiers and thus enhancing the energy efficiency. The current state-
of-the-art wireless systems, including conventional communication systems such
as massive multiple-input multiple-output (MIMO), millimeter wave (mmWave)
communication, and ultra-dense deployments of small cells consume a high
amount of energy and a considerable hardware implementation cost [5]. The
emerging mmWave communications may suffer from high path-losses, penetration
losses, and signal blockages. To this end, IRS-aided communication systems can
potentially address the need of greener, more cost effective, and more
sustainable wireless connectivity by enabling control over the propagation
environment by recycling EM signals. The recent developments in physics of
meta-surfaces [4] may contribute to cost effective implementation of IRS-aided
wireless systems. Thus, by intelligently controlling the phase-shifts
introduced by the reflective elements, the IRS-assisted system can potentially
overcome unfavorable propagation conditions, increase the coverage area while
retaining a higher energy efficiency due to the passive reflective elements.

###  I-A Related Prior Research on Intelligent Reflective Surfaces

In [2, 4, 6, 7, 3], the initial designs and theoretical modeling aspects of
IRS-aided wireless communication are investigated. To this end, reference [2]
explores the core architectural designs and features of IRS-based wireless
communications. In [3], a new communication-theoretic model for the analysis
and optimization of IRS-assisted smart radio environments is proposed. In [4],
prototypes of meta-surfaces and meta-tiles that can be used to coat objects
embedded within a smart wireless propagation environment are already
developed. In [6], the ray-tracing techniques are used to design novel
propagation/path-loss models for IRS-assisted wireless communications. In [7],
the path-loss models for IRS channels are designed based on the EM properties,
and, these models are categorized into far-field, near-field beamforming, and
near-field broadcasting.

References [8, 9, 10, 11, 12, 13, 14, 15] provide performance analyses for
IRS-enable networks. Specifically, in [9], the fundamental performance metrics
of IRS-assisted communication systems operating over Rayleigh fading are
derived. In [8], novel IRS-aided techniques are explored to boost the physical
layer security aspects of wireless communications. In [10], the potential of
the large intelligent surface (LIS)-assisted communication in boosting the
performance is investigated. Reference [11] explores the feasibility of
exploiting statistical channel state information (CSI) to optimize phase-shift
designs at the IRS to maximize the achievable spectral efficiency. In [12],
the asymptotically achievable uplink sum rate of an IRS-aided system is
analyzed by deriving the rate distributions for Rician fading channels. In
[13], the impact of practically-viable/limited IRS control channels is
investigated, and a passive beamformer design, which can achieve the
asymptotic optimal received SNR under discrete reflection phase-shifts, is
presented. In [13], a resource allocation algorithm to maximize the asymptotic
sum rate under the passive beamforming is proposed. In [14], a performance
comparison between the spherical and two-dimensional (2D) LIS deployments is
presented, and it is shown that former can provide wider coverage while having
simpler positioning and flexible implementation over the 2D counterpart. In
[15], an analytical framework for random rotation based IRS-assisted systems
is presented, and four low-complexity schemes via selection-based and coding-
based approaches are proposed.

In [16, 17, 18, 19, 20, 21, 22, 23], the optimization-based algorithm designs
are developed to intelligently control the phases of passive reflective
elements of an IRS to enable a myriad of wireless communication applications.
In [16, 17], the joint optimization of precoders at the base-station (BS) and
phase-shifts at the IRS is investigated, and the proposed solutions are based
on the semidefinite relaxation and alternating optimization techniques.
Reference [18] presents two computationally efficient algorithms to maximize
the energy efficiency for the transmit power allocation based on alternating
maximization adopting gradient descent method and sequential fractional
programming approach. In [19], the beamforming designs are investigated with a
practical phase-shift at the IRS. In [20], the capacity regions of the IRS-
aided two user multiple-access channel are characterized. Thereby, the
algorithms to find the inner and outer bounds of the capacity region are
proposed for a centralized IRS deployment. Reference [21] investigates an IRS-
assisted wireless system operating in Rician fading channels by adopting
maximal ratio transmission (MRT) at the BS. In [22], the joint optimization of
active transmit and passive reflect beamforming is investigated to maximize
the secrecy rate for an IRS-assisted secure communication system. In [23], the
impact of artificial noise on IRS-assisted secure communication is
investigated.

This paper goes well beyond our related conference paper [1] by presenting the
average symbol error rate (SER), asymptotic outage/SER analysis in high SNR
regime, Kullback-Leibler (KL) divergence, diversity order, array/coding gain,
achievable rate limits in large reflective element regime, and transmit power
scaling laws for the distributed IRSs-aided communication systems.

###  I-B Motivation and Our Contribution

The main contribution of this work is to derive performance metrics for the
distributed IRSs deployment over Nakagami-$m$ fading in closed-form. A common
attribute of the above related prior researches [6, 16, 17, 9, 11, 8, 7, 10,
20, 13, 12, 21, 14, 15, 18, 19] is that their system models consist of a
single IRS. However, the concept of IRS was originally envisioned to coat
physical objects distributed within a wireless propagation channel. The
fundamental performance metrics of a distributed IRSs-aided set-up have not
yet been explored. To fill this important gap in IRS literature, we present a
statistical characterization of the end-to-end optimal SNR and the
corresponding performance metrics of a distributed IRSs-aided system over
Nakagami-$m$ fading. The optimal SNR is tightly approximated by using a
mathematically tractable counterpart via the central limit theorem (CLT) [24]
because the optimal counterpart does not seem to be amenable to closed-form
characterization of the performance metrics. Then, the exact probability
density function (PDF) and cumulative distribution function (CDF) of this
optimal SNR approximation are derived. We analytically validate the tightness
of our PDF/CDF approximations via the KL divergence. By using our PDF/CDF
analysis, the outage probability, the average SER, and achievable rate
bounds/approximations are derived in closed-form. In order to obtain further
insights, we report to an asymptotic analysis of the outage probability and
average SER in the high SNR regime. Thereby, the achievable diversity order
and array/coding gain are quantified. Thus, diversity gains can also be
achieved via passive reflections at the distributed IRSs deployment by
recycling existing EM waves without generating additional EM waves by active
RF chains. Our rate analysis in the large reflective element regime reveals
that our lower and upper bounds for the achievable rate are asymptotically
exact. The accuracy of our analysis is verified via Monte-Carlo simulations,
and a set of rigorous numerical results is presented to investigate the
performance of the distributed IRSs. Our results reveal that a distributed
deployment of IRSs can be exploited to boost the outage, average SER
performances, and achievable rate compared to the direct transmission and
single IRS-aided set-ups.

Notation: $\mathbf{x}^{T}$ denotes the transpose of $\mathbf{x}$.
$\mathbb{E}\\!\left[{X}\right]$ and $\mathbb{V}\mathrm{ar}\\!\left[{X}\right]$
represent the expectation and variance of a random variable (RV) $X$,
respectively. $X\sim\mathcal{N}\left(\mu_{X},\sigma_{X}^{2}\right)$ denotes
that $X$ is Gaussian distributed with $\mu_{X}$ mean and $\sigma_{X}^{2}$
variance. $\mathcal{O}(\cdot)$ denotes the higher order terms of the MacLaurin
series expansion [25, Eq. (0.318.2)]. $\mathcal{Q}(\cdot)$ is the
Gaussian-$\mathcal{Q}$ function [24], $\Gamma(t)$ is the Gamma function [25,
Eq. (8.310.1)], $\Gamma(\alpha,x)$ is the upper incomplete Gamma function [25,
Eq. (8.350.2)], and $\gamma(\alpha,x)$ denotes the lower incomplete Gamma
function [25, Eq. (8.350.1)]. Moreover, $\mathcal{K}_{v}(\cdot)$ represents
the modified Bessel function of the second kind [25, Eq. (8.407)] and
$F(\cdot,\cdot;\cdot;\cdot)$ is the Gauss Hypergeometric function [25, Eq.
(9.100)].

Destination $(D)$Source
$(S)$$\mathrm{IRS_{1}}$$\mathrm{IRS_{N}}$$\mathrm{IRS_{n}}$${\mathbf{h}_{1}}$${\mathbf{h}_{n}}$${\mathbf{h}_{N}}$${{u}}$${\mathbf{g}_{1}}$${\mathbf{g}_{n}}$${\mathbf{g}_{N}}$
Figure 1: System model - a distributed IRSs-aided communication set-up.

##  II System, Channel and Signal Models

###  II-A System and Channel Model

We consider a distributed IRSs-assisted wireless system having a single-
antenna source ($S$), a single-antenna destination ($D$), and $N$ distributed
IRSs with $L$ passive reflective elements at each IRS (see Fig. 1). The $S$
serves $D$ through $N$ distributed IRSs and via direct channel. It is assumed
that the phase-shifts at the IRS reflective elements can be intelligently
controlled such that the received signal at $D$ can be constructively added by
using an IRS controller, which connects the IRSs and the BS via a backhaul
connection [16, 17]. The channels between $S$ and the $l$th reflective element
in the $n$th IRS and the $l$th reflective element in the $n$th IRS and $D$ are
denoted by $h_{nl}$ and $g_{nl}$, respectively. Moreover, $u$ denotes the
direct channel between $S$ and $D$. All channel envelopes are assumed to be
independent Nakagami-$m$ distributed, where $m$ is the shape parameter. The
direct channel ($u$) can be written in its polar form as

$\displaystyle u=\alpha_{u}\mathrm{e}^{j\theta_{u}},$ (1)

where $\alpha_{u}$ is the envelope of $u$ and $\theta_{u}$ is the phase of
$u$. The PDF of $\alpha_{u}$ is given by [24]

$\displaystyle
f_{\alpha_{u}}(x)=\frac{2m_{u}^{m_{u}}x^{2m_{u}-1}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}\mathrm{exp}\left(\frac{-m_{u}x^{2}}{\xi_{u}}\right),$
(2)

where $m_{u}$ is the shape parameter and $\xi_{u}\\!=\\!m_{u}\zeta_{u}$ is the
scaling parameter. Here, $\zeta_{u}$ accounts for the large-scale fading/path-
loss. The $S$-to-IRS and IRS-to-$D$ channels can be defined as

$\displaystyle v_{nl}=\alpha_{v_{nl}}\mathrm{e}^{j\theta_{v_{nl}}},$ (3)

where $v\\!\in\\!\\{h,g\\}$, $\alpha_{v_{nl}}$ is the envelope of $v_{nl}$,
and $\theta_{v_{nl}}$ is the phase. The PDF of $\alpha_{v_{nl}}$ is given as

$\displaystyle
f_{\alpha_{v_{nl}}}(x)=\frac{2m_{v}^{m_{v}}x^{2m_{v}-1}}{\Gamma(m_{v})\xi_{v_{n}}^{m_{v}}}\mathrm{exp}\left(\frac{-m_{v}x^{2}}{\xi_{v_{n}}}\right),$
(4)

where $m_{v_{nl}}=m_{v}$ and
$\xi_{v_{nl}}=m_{v_{nl}}\zeta_{v_{nl}}=\xi_{v_{n}}$ are the shape and scaling
parameters, respectively. It is assumed that large-scale fading parameters are
the same for a given IRS because its reflective elements are co-located;
$\xi_{v_{nl}}=\xi_{v_{n}},\;\forall l$. However, $\xi_{v_{n}}$ depends on the
IRS index ($n$) because geographically distributed IRSs deployment.

###  II-B Signal Model

The signal transmitted by $S$ reaches at $D$ through the reflected channels of
$N$ distributed IRSs and the direct channel. The signal received at $D$ can be
written as

$\displaystyle
r=\sqrt{P}\left(u+\sum_{n=1}^{N}\mathbf{g}^{T}_{n}\mathbf{\Theta}_{n}\mathbf{h}_{n}\right)x+w,$
(5)

where $x$ is the transmitted signal at $S$ satisfying
$\mathbb{E}\\!\left[{|x|^{2}}\right]=1$, $P$ is the transmit power, and $w$ is
an additive white Gaussian noise (AWGN) at $D$ having zero mean and variance
of $\sigma_{w}^{2}$ such that $w\sim\mathcal{CN}(0,\sigma_{w}^{2})$. In (5),
$\mathbf{h}_{n}=[h_{n1},\cdots,h_{nl},\cdots,h_{nL}]^{T}\in\mathbb{C}^{L\times
1}$ is the channel vector between $S$ and the $n$th IRS and
$\mathbf{g}^{T}_{n}=[g_{n1},\cdots,g_{nl},\cdots,g_{nL}]\in\mathbb{C}^{1\times
L}$ accounts for the channel vector between the $n$th IRS and $D$. Moreover,
$\mathbf{\Theta}_{n}=\mathrm{diag}\left(\eta_{n1}\mathrm{e}^{j\theta_{n1}},\cdots,\eta_{nl}\mathrm{e}^{j\theta_{nl}},\cdots,\eta_{nL}\mathrm{e}^{j\theta_{nL}}\right)\in\mathbb{C}^{L\times
L}$ is a diagonal matrix, which captures the reflection properties of the
$n$th IRS. Here, $\eta_{nl}\mathrm{e}^{j\theta_{nl}}$ is a complex-valued
reflection coefficient at the $l$th reflective element of the $n$th IRS, where
$\eta_{nl}$ and $\theta_{nl}$ are the magnitude of attenuation and phase-
shift, respectively. Thus, the signal received at $D$ in (5) can be rewritten
as

$\displaystyle
r=\sqrt{P}\left(u+\sum_{n=1}^{N}\sum_{l=1}^{L}g_{nl}\eta_{nl}\mathrm{e}^{j\theta_{nl}}h_{nl}\right)x+w.$
(6)

The SNR at $D$ can be derived via (6) as

$\displaystyle\gamma=\bar{\gamma}\left|u+\sum_{n=1}^{N}\sum_{l=1}^{L}g_{nl}\eta_{nl}\mathrm{e}^{j\theta_{nl}}h_{nl}\right|^{2},$
(7)

where $\bar{\gamma}=P/\sigma^{2}_{w}$ is the transmit SNR. Then, by
substituting (1) and (3) into (7), this SNR can be written in terms of the
channel phases as

$\displaystyle\gamma=\bar{\gamma}\left|\alpha_{u}\mathrm{e}^{j\theta_{u}}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\alpha_{g_{nl}}\alpha_{h_{nl}}\mathrm{e}^{j\left(\theta_{nl}+\theta_{g_{nl}}+\theta_{h_{nl}}\right)}\right|^{2}.$
(8)

This analysis (8) reveals that the received SNR at $D$ can be maximized when
$NL$ signal terms inside the summation of (8) are constructively added to the
signal component received via the direct channel. By controlling the phase-
shifts at each IRS reflective element ($\theta_{nl}$), the phases inside the
double summation in (8) can be adjusted to enable a constructive addition of
the received signal components via the direct and reflected channels. Thus,
the optimal choice of $\theta_{nl}$ to maximize the received SNR at $D$ can be
given by [16]

$\displaystyle\theta_{nl}^{*}=\underset{-\pi\leq\theta_{nl}\leq\pi}{\mathrm{argmax}}\;{\gamma}=\theta_{u}-\left(\theta_{g_{nl}}+\theta_{h_{nl}}\right),$
(9)

for $n\in\\{1,\cdots,N\\}$ and $l\in\\{1,\cdots,L\\}$. By using (9), the
optimal SNR at $D$ can be derived as

$\displaystyle\gamma^{*}=\bar{\gamma}\left[\alpha_{u}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\alpha_{g_{nl}}\alpha_{h_{nl}}\right]^{2}.$
(10)

##  III Preliminary Analysis

###  III-A Statistical Characterization of the Optimal Received SNR

The $\alpha_{g_{nl}}$ and $\alpha_{h_{nl}}$ in (10) for $n\in\\{1,\cdots,N\\}$
and $l\in\\{1,\cdots,L\\}$ are independently distributed Nakagami RVs, and
hence, the exact derivations of the PDFs of

$\displaystyle
Y=\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\alpha_{g_{nl}}\alpha_{h_{nl}},\quad\text{and}\quad{\gamma}^{*}=\bar{\gamma}{[\underbrace{\alpha_{u}+Y}_{R}]}^{2}=\bar{\gamma}R^{2},$
(11)

seem mathematically involved and may not provide useful design insights.
Nevertheless, even for moderately large values of the product $NL$, $Y$ can be
tightly approximated by an one-sided Gaussian distributed RV ($\tilde{Y}$) by
invoking the CLT [24]. Then, an approximated PDF for $Y$ can be written as
(see Appendix A)

$\displaystyle f_{Y}(y)\approx
f_{\tilde{Y}}(y)=\frac{\psi}{\sqrt{2\pi\sigma_{Y}^{2}}}\mathrm{exp}\left(-\frac{(y-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}\right),\quad
y\geq 0,$ (12)

where $f_{\tilde{Y}}(y)=0$ for $y<0$, $\psi\triangleq
1/\mathcal{Q}\left(-\mu_{Y}/\sigma_{Y}\right)$ is a normalization factor to
satisfy $\int_{-\infty}^{\infty}f_{\tilde{Y}}(x)dx=1$, and
$\mathcal{Q}(\cdot)$ is the Gaussian-$\mathcal{Q}$ function [24]. In (12),
$\mu_{Y}$ and $\sigma_{Y}^{2}$ can be derived via the moment matching
technique as (see Appendix A)

$\displaystyle\mu_{Y}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\sqrt{\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\frac{\Gamma\left(m_{h}+1/2\right)\Gamma\left(m_{g}+1/2\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)},$
(13a) $\displaystyle\sigma_{Y}^{2}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}^{2}\left({\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\right)\frac{\Gamma\left(m_{h}+1\right)\Gamma\left(m_{g}+1\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)}-\mu_{Y}^{2},$
(13b)

where $\Gamma(t)=\int_{0}^{\infty}x^{t}\mathrm{e}^{-x}dx$ is the Gamma
function [25, Eq. (8.310.1)]. Then, a tight approximation of the PDF of $R$ in
(11) or the exact PDF of its approximation $\tilde{R}=\alpha_{u}+\tilde{Y}$
can be derived as (see Appendix B)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!f_{R}(x)\approx
f_{\tilde{R}}(x)=\lambda\mathrm{e}^{-\Delta\left(\frac{x-\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)^{\\!2}}\sum_{k=0}^{2m_{u}\\!-\\!1}\\!\\!\binom{2m_{u}\\!-\\!1}{k}\\!\\!\left(\frac{x\\!-\\!\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)^{\\!\\!2m_{u}\\!-\\!1\\!-\\!k}\\!\\!\left[\Gamma\\!\left(\\!\frac{k\\!+\\!1}{2},\\!\left(\\!\frac{x\\!-\\!\mu_{Y}}{2\sigma_{Y}^{2}\sqrt{a}}\right)^{\\!\\!2}\right)\right],$
(14)

where $\Gamma(\alpha,x)=\int_{x}^{\infty}\mathrm{e}^{-t}t^{\alpha-1}dt$ is the
upper incomplete Gamma function [25, Eq. (8.350.2)]. Here, $a$, $\lambda$, and
$\Delta$ are given by

$\displaystyle\\!\\!\\!\\!\\!\\!a$ $\displaystyle=$
$\displaystyle\frac{m_{u}}{\xi_{u}}+\frac{1}{2\sigma_{Y}^{2}},\quad\text{and}\quad\lambda=\frac{m_{u}^{m_{u}}\psi}{\Gamma(m_{u})\xi_{u}^{m_{u}}a^{m_{u}}\sqrt{2\pi\sigma_{Y}^{2}}},\quad\text{and}\quad\Delta=\left(\frac{1}{2\sigma_{Y}^{2}}-\frac{1}{4a\sigma_{Y}^{4}}\right)4\sigma_{Y}^{4}a.$
(15)

Via $\tilde{R}=\alpha_{u}+\tilde{Y}$ in (14), a tight approximation to the
optimal SNR ($\gamma^{*}$) can be written as

$\displaystyle\gamma^{*}\approx\tilde{\gamma}^{*}=\bar{\gamma}\tilde{R}^{2}.$
(16)

Then, an approximation for the PDF of $\gamma^{*}=\bar{\gamma}R^{2}$ or the
exact PDF of $\tilde{\gamma}^{*}=\bar{\gamma}\tilde{R}^{2}$ in (16) can be
derived by using (14) as [24]

$\displaystyle f_{\gamma^{*}}(y)$ $\displaystyle\approx$ $\displaystyle
f_{\tilde{R}}\left(\sqrt{{y}/{\bar{\gamma}}}\right)\Big{/}{2\sqrt{\bar{\gamma}y}}.$
(17)

Next, the CDF of $\tilde{R}$ can be derived as (see Appendix C)

$\displaystyle F_{\tilde{R}}(x)$ $\displaystyle=$ $\displaystyle
1-\int_{x}^{\infty}f_{\tilde{R}}(u)du=1-\frac{\lambda}{2a^{m_{u}}}\sum_{k=0}^{2m_{u}-1}\binom{2m_{u}-1}{k}I_{k},$
(18)

where $I_{k}$ is given as

$\displaystyle I_{k}=\begin{cases}I_{k}^{o},&\quad\text{odd $k$ },\\\
I_{k}^{e},&\quad\text{even $k$}.\end{cases}$ (19)

In (19), $I_{k}^{o}$ is given by

$\displaystyle\\!\\!\\!\\!\\!I_{k}^{o}\\!=\\!\\!\begin{cases}\\!\frac{q\left(\gamma_{o}\\!-\\!1\right)!}{2}\\!\\!\sum\limits_{i=0}^{\gamma_{o}\\!-\\!1}\\!\\!\frac{\left(\Delta+1\right)^{\frac{k}{2}-m_{u}-i}}{i!}\\!\left(2\Gamma\\!\left(m_{u}\\!+\\!i+\\!\frac{k}{2}\right)-\Gamma\\!\left(m_{u}\\!+\\!i\\!-\\!\frac{k}{2},\left(\Delta\\!+\\!1\right)l_{min}^{2}\right)\right),&\text{for}\quad
x<\mu_{Y},\\\
\\!\frac{q\left(\gamma_{o}-1\right)!}{2}\sum\limits_{i=0}^{\gamma_{o}-1}\frac{\left(\Delta+1\right)^{\frac{k}{2}-m_{u}-i}}{i!}\Gamma\left(m_{u}+i-\frac{k}{2},\left(\Delta+1\right)l_{min}^{2}\right),&\text{for}\quad
x>\mu_{Y},\end{cases}$ (20)

where $q=2\sigma_{Y}^{2}\sqrt{a}$, $\gamma_{o}=(k+1)/2$, and
$l_{min}=(x-\mu_{Y})/q$. Here, $I_{k}^{e}$ in (19) can be given as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!I_{k}^{e}$ $\displaystyle=$
$\displaystyle\frac{q(\gamma_{e}-1)!}{2}\sum\limits_{j=0}^{\gamma_{e}-1}\frac{\Delta^{j-\gamma_{e}}}{j!}\left(\Gamma\left(\frac{k+1}{2},l_{min}^{2}\right)l_{min}^{2j}\mathrm{e}^{-\Delta
l_{min}^{2}}-\frac{\Gamma\left(j+\frac{k}{2}+\frac{1}{2},(\Delta+1)l_{min}^{2}\right)}{(\Delta+1)^{j+\frac{k}{2}+\frac{1}{2}}}\right),$
(21)

where $\gamma_{e}=m-k/2$. Then, the CDF of $\gamma^{*}=\bar{\gamma}R^{2}$ can
be approximated as [24]

$\displaystyle F_{\gamma^{*}}(y)$ $\displaystyle=$
$\displaystyle\mathrm{Pr}\left(\gamma^{*}\leq y\right)\approx
F_{\tilde{R}}\left(\sqrt{y/\bar{\gamma}}\right).$ (22)

Figure 2: The PDF and CDF of the received SNR ($\gamma^{*}$) for
$\bar{\gamma}=-10$dB. The combinations of $N$ and $L$ for Case 1 to Case 5 are
set to $\\{N=2,L=32\\}$, $\\{N=6,L=32\\}$, $\\{N=2,L=64\\}$, $\\{N=8,L=32\\}$,
and $\\{N=2,L=128\\}$, respectively.

Remark 1: The accuracy of the approximated PDF and CDF of $\gamma^{*}$ is
verified by plotting (17), (22) and the Monte-Carlo simulations of the exact
counterparts in Fig. 2. It reveals that our analytical approximations for the
PDF (17) and CDF (22) of $\gamma^{*}$ are accurate even for moderately large
number of reflective elements $(L)$ at the IRSs. Since IRSs are typically made
of cost-effective reflective elements [26], a moderately large $N$ is
practically feasible.

###  III-B Kullback-Leibler (KL) Divergence

In this section, the accuracy of our statistical characterization of the
received SNR in Section III-A is investigated by using the KL divergence [27,
28]. Specifically, we analytically quantify the KL divergence between the
exact PDF of $z_{nl}=\alpha_{h_{nl}}\alpha_{g_{nl}}$ and its approximated PDF
as one-sided Gaussian distribution (12). Thus, the KL divergence can be
defined as [27, 28]

$\displaystyle\mathcal{D}_{KL}=\mathbb{E}\\!\left[{\mathrm{log}\left({f_{z_{nl}}^{exct}(x)}\big{/}{f_{z_{nl}}^{app}(x)}\right)}\right],$
(23)

where the expectation is taken with respect to the exact PDF of $z_{nl}$. The
exact PDF of $z_{nl}$ is given as (see Appendix A)

$\displaystyle f_{z_{nl}}^{exct}(x)$ $\displaystyle=$
$\displaystyle\alpha^{\prime}x^{m_{h}+m_{g}-1}\mathcal{K}_{m_{h}-m_{g}}\left(2x\sqrt{\frac{m_{h}m_{g}}{\xi_{h_{n}}\xi_{g_{n}}}}\right),\quad
x\geq 0,$ (24)

where $\alpha^{\prime}$ is given in (55b). The approximated PDF of $z_{nl}$
can be derived via (12) as

$\displaystyle
f_{z_{nl}}^{app}(x)=\frac{\psi_{z_{nl}}}{\sqrt{2\pi\sigma_{z_{nl}}^{2}}}\mathrm{exp}\left(-\frac{(x-\mu_{z_{nl}})^{2}}{2\sigma_{z_{nl}}^{2}}\right),\quad
x\geq 0,$ (25)

where $\psi_{z_{nl}}\triangleq
1/\mathcal{Q}\left(-\mu_{z_{nl}}/\sigma_{z_{nl}}\right)$ is a normalization
factor. Moreover, $\mu_{z_{nl}}$ and $\sigma_{z_{nl}}^{2}$ are given as (see
Appendix A)

$\displaystyle\mu_{z_{nl}}$ $\displaystyle=$
$\displaystyle\sqrt{\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\frac{\Gamma\left(m_{h}+1/2\right)\Gamma\left(m_{g}+1/2\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)},$
(26a) $\displaystyle\sigma_{z_{nl}}^{2}$ $\displaystyle=$
$\displaystyle\left({\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\right)\frac{\Gamma\left(m_{h}+1\right)\Gamma\left(m_{g}+1\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)}-\mu_{z_{nl}}^{2}.$
(26b)

Then, the KL divergence can be given as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mathcal{D}_{KL}=\digamma+\mathbb{E}\\!\left[{\mathrm{log}\left(x^{m_{h}+m_{g}-1}\right)}\right]+\mathbb{E}\\!\left[{\mathrm{log}\left(\mathcal{K}_{m_{h}-m_{g}}\left(2x\sqrt{\frac{m_{h}m_{g}}{\xi_{h_{n}}\xi_{g_{n}}}}\right)\right)}\right]+\mathbb{E}\\!\left[{\frac{(x-\mu_{z_{nl}})^{2}}{2\sigma_{z_{nl}}^{2}}}\right],$
(27)

where
$\digamma=\mathrm{log}\left(\alpha^{\prime}\right)-\mathrm{log}\left({\psi_{z_{nl}}}/{\sqrt{2\pi\sigma_{z_{nl}}^{2}}}\right)$.

Figure 3: The KL divergence of the approximated PDF of $z_{nl}$ for
$m_{h}=m_{g}=m\in\\{2,3,5\\}$.

Remark 2: In Fig. 3, the KL divergence of the approximated PDF of $z_{nl}$ is
plotted as a function of $\zeta_{c}=\xi_{h_{nl}}/m_{h}=\xi_{g_{nl}}/m_{g}$.
Since $\zeta_{c}$ accounts for the large-scale fading, the range of
$\zeta_{c}$ is set to $(0,1)$. Fig. 3 clearly shows that the KL divergence of
the approximated PDF for $m=3$ is as low as $\mathcal{D}_{KL}\approx 5.2\times
10^{-2}$ for the entire range of $\zeta_{c}$. This observation thus validates
the tightness of our approximated PDF of the received optimal SNR to the exact
counterpart.

##  IV Performance Analysis

###  IV-A Outage Probability

The probability that the instantaneous SNR falls below a threshold SNR
($\gamma_{th}$) is defined as the SNR outage probability, and a tight
approximation to it can be obtained from (22) as

$\displaystyle P_{out}=P_{r}\left(\gamma\leq\gamma_{th}\right)\approx
F_{\gamma^{*}}(\gamma_{th}).$ (28)

This definition for $P_{out}$ is used to quantify the asymptotic outage
probability and the achievable diversity order in Section IV-D.

###  IV-B Average Achievable Rate

The average achievable rate can be defined as follows:

$\displaystyle\mathcal{R}=\mathbb{E}\\!\left[{\mathrm{log}_{2}\left(1+\gamma^{*}\right)}\right].$
(29)

Since the exact derivation of $\mathcal{R}$ in (29) seems mathematically
intractable, we resort to tight upper and lower bounds by invoking the
Jensen’s inequality as [29]

$\displaystyle\mathcal{R}_{lb}\leq\mathcal{R}\leq\mathcal{R}_{ub}.$ (30)

In (30), $\mathcal{R}_{lb}$ and $\mathcal{R}_{ub}$ are defined as

$\displaystyle\mathcal{R}_{lb}$ $\displaystyle=$
$\displaystyle\mathrm{log}_{2}\left(1+\left(\mathbb{E}\\!\left[{1/\gamma^{*}}\right]\right)^{-1}\right)\approx\mathrm{log}_{2}\left(1+\left(\mathbb{E}\\!\left[{1/\tilde{\gamma}^{*}}\right]\right)^{-1}\right),$
(31a) $\displaystyle\mathcal{R}_{ub}$ $\displaystyle=$
$\displaystyle\mathrm{log}_{2}\left(1+\mathbb{E}\\!\left[{\gamma^{*}}\right]\right)\approx\mathrm{log}_{2}\left(1+\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right).$
(31b)

The expectation term in (31b) can be approximately derived as (see Appendix
D-A)

$\displaystyle\mathbb{E}\\!\left[{\gamma^{*}}\right]\approx\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]=\bar{\gamma}\left(\sigma_{u}^{2}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2}\right),$
(32)

where $\mu_{Y}$ and $\sigma_{Y}^{2}$ are defined in (13a) and (13b),
respectively. In (32), $\mu_{u}$ and $\sigma_{u}^{2}$ are defined as

$\displaystyle\mu_{u}$ $\displaystyle=$
$\displaystyle\frac{\Gamma\left(m_{u}+1/2\right)}{\Gamma\left(m_{u}\right)}\left({\frac{\xi_{u}}{m_{u}}}\right)^{1/2},$
(33a) $\displaystyle\sigma_{u}^{2}$ $\displaystyle=$
$\displaystyle\xi_{u}\left(1-\frac{1}{m_{u}}\left(\frac{\Gamma\left(m_{u}+1/2\right)}{\Gamma\left(m_{u}\right)}\right)^{2}\right).$
(33b)

Next, the expectation term in (31a) can be computed as

$\displaystyle\mathbb{E}\\!\left[{1/\gamma^{*}}\right]\approx\mathbb{E}\\!\left[{1/\tilde{\gamma}^{*}}\right]={1}\big{/}{\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]}+{\mathbb{V}\mathrm{ar}\\!\left[{\tilde{\gamma}^{*}}\right]}\big{/}{\left(\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right)^{3}},$
(34)

where $\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]$ is defined in (32).
Here, $\mathbb{V}\mathrm{ar}\\!\left[{\tilde{\gamma}^{*}}\right]$ is derived
as (see Appendix D-B)

$\displaystyle\mathbb{V}\mathrm{ar}\\!\left[{\tilde{\gamma}^{*}}\right]=\bar{\gamma}^{2}\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]-\left(\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\right)^{2},$
(35)

where $\mathbb{E}{[\tilde{R}^{4}]}$ is given by

$\displaystyle\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]$ $\displaystyle=$
$\displaystyle\sum_{n=0}^{4}\binom{4}{n}\mathbb{E}\\!\left[{\alpha_{u}^{(4-n)}}\right]\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right].$
(36)

Here, $\mathbb{E}\\!\left[{\alpha_{u}^{n}}\right]$ for $n\in\\{1,2,3,4\\}$ can
be derived as

$\displaystyle\mathbb{E}\\!\left[{\alpha_{u}^{n}}\right]=\frac{\Gamma\left(m_{u}+n/2\right)}{\Gamma\left(m_{u}\right)}\left(\frac{\xi_{u}}{m_{u}}\right)^{n/2}.$
(37)

Moreover, $\mathbb{E}{[\tilde{Y}^{n}]}$ for $n\in\\{1,2,3,4\\}$ is given by

$\displaystyle\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]=\frac{\psi}{2\sqrt{\pi}}\sum_{i=0}^{n}\binom{n}{i}\left(\sqrt{2\sigma_{Y}^{2}}\right)^{n-i}\mu_{Y}^{i}I\left(n-i,\frac{-\mu_{Y}}{2\sigma_{Y}^{2}}\right),$
(38)

where $I(\cdot,\cdot)$ is given as

$\displaystyle
I\left(m,t\right)=\begin{cases}(-1)^{m}\gamma\left(\frac{m+1}{2},t^{2}\right)+\Gamma\left(\frac{m+1}{2}\right),&\text{for}\,\,t\leq
0,\\\ \Gamma\left(\frac{m+1}{2},t^{2}\right),&\text{otherwise},\end{cases}$
(39)

where $\gamma(\alpha,x)=\int_{0}^{x}\mathrm{e}^{-t}t^{\alpha-1}dt$ is the
lower incomplete Gamma function [25, Eq. (8.350.1)]. Finally,
$\mathcal{R}_{lb}$ and $\mathcal{R}_{ub}$ can be derived as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mathcal{R}_{lb}=\mathrm{log}_{2}\left(\\!\\!1\\!+\\!\frac{\bar{\gamma}\left(\sigma_{u}^{2}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2}\right)^{3}}{\sum\limits_{n=0}^{4}\binom{4}{n}\frac{\Gamma\left(m_{u}+\frac{(4-n)}{2}\right)}{\Gamma\left(m_{u}\right)}\left(\frac{\xi_{u}}{m_{u}}\right)^{(4-n)/2}\\!\\!\\!\frac{\psi}{2\sqrt{\pi}}\sum\limits_{i=0}^{n}\binom{n}{i}\left({2\sigma_{Y}^{2}}\right)^{(n-i)/2}\mu_{Y}^{i}I\left(n-i,\frac{-\mu_{Y}}{2\sigma_{Y}^{2}}\right)}\right)\\!\\!,$
(40)

$\displaystyle\mathcal{R}_{ub}=\mathrm{log}_{2}\left(1+\bar{\gamma}\left(\sigma_{u}^{2}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2}\right)\right),$
(41)

where $\sigma_{u}^{2}$, $\sigma_{Y}^{2}$, $\mu_{u}$, and $\mu_{Y}$ are defined
in (33b), (13b), (33a), and (13a), respectively.

####  IV-B1 Asymptotic achievable rate as the number of IRS elements grows
large $(L\rightarrow\infty)$

To obtain further insights, the asymptotic achievable rate can be derived as
the number of IRS elements grows without bound $(L\rightarrow\infty)$. In this
operating regime, we observe that the transmit power can be scaled inversely
proportional to the square of the number of IRS elements:
$\lim_{L\rightarrow\infty}P=P_{E}/L^{2}$. By using this transmit power scaling
law, the lower and upper rate bounds in (40) and (41), respectively, can be
shown to approach an asymptotic limit as follows (see Appendix D-C):

$\displaystyle\lim_{L\rightarrow\infty}\mathcal{R}_{lb}\longrightarrow\mathcal{R}^{\infty},\qquad\text{and}\qquad\lim_{L\rightarrow\infty}\mathcal{R}_{ub}\longrightarrow\mathcal{R}^{\infty},$
(42)

where the asymptotic achievable rate $\mathcal{R}^{\infty}$ in large $L$
regime is given by

$\displaystyle\lim_{L\rightarrow\infty}\mathcal{R}\longrightarrow\mathcal{R}^{\infty}=\mathrm{log}_{2}\left(1+\bar{\gamma}_{E}\left(\sum_{n=1}^{N}\eta_{n}\sqrt{\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\frac{\Gamma\left(m_{h}+1/2\right)\Gamma\left(m_{g}+1/2\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)}\right)^{2}\right),$
(43)

where $\bar{\gamma}_{E}=P_{E}/\sigma^{2}_{w}$ and $\eta_{n}=\eta_{nl}$,
$\forall n$.

Remark 3: Our asymptotic analysis reveals that our lower (40) and upper (41)
rate bounds converge to a common limit, and hence, these rate bounds are
asymptotically exact in large $L$ regime. Thus, a finite rate at $D$ is
achievable once the transmit power at $S$ can be scaled inversely proportional
to $L^{2}$.

###  IV-C Average Symbol Error Rate (SER)

The average SER is the expectation of conditional error probability
($P_{e|\gamma^{*}}$) over the probability distribution of $\gamma^{*}$ [30].
For a wide-range of coherent modulation schemes, $P_{e|\gamma^{*}}$ is given
by [30]

$\displaystyle
P_{e|\gamma^{*}}=\omega\mathcal{Q}\left(\sqrt{\vartheta\gamma^{*}}\right),$
(44)

where $\omega$ and $\vartheta$ depend on the modulation scheme [30]. For
instance, the pair of parameters ($\omega,\vartheta$) for binary phase-shift
keying (BPSK) and quadrature phase-shift keying (QPSK) are given by
($\omega=1,\vartheta=2$) and ($\omega=2,\vartheta=1$), respectively [30].
Then, the average SER can be defined as

$\displaystyle\bar{P}_{e}=\mathbb{E}\\!\left[{\omega\mathcal{Q}\left(\sqrt{\vartheta\gamma^{*}}\right)}\right].$
(45)

Then, by using $\tilde{\gamma}^{*}$ in (16), a tight approximation for the
average SER can be derived as

$\displaystyle\bar{P}_{e}\approx\tilde{P}_{e}=\omega\mathbb{E}_{\alpha_{u}}\\!\left[{\mathbb{E}_{\tilde{Y}}\\!\left[{\mathcal{Q}\left(\vartheta_{b}(\alpha_{u}+\tilde{Y})\right)}\right]}\right],$
(46)

where $\vartheta_{b}\\!=\\!\sqrt{\vartheta\bar{\gamma}}$. A tight
approximation to the average SER $\tilde{P}_{e}$ is derived as (see Appendix
E)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\tilde{P}_{e}$ $\displaystyle\approx$
$\displaystyle\omega
B\sum_{k=0}^{2m_{u}-1}\\!\\!\left(\frac{-u_{1}}{2v_{1}}\right)^{2m_{u}-1-k}\\!\\!\\!\\!\\!\Gamma(\gamma_{b})\sum_{i=0}^{\gamma_{b}-1}\\!\\!\frac{v_{1}^{i-\gamma_{b}}}{i!}\left(\\!\mathcal{Q}(\sqrt{2}r)\left(\frac{u_{1}}{2v_{1}}\right)^{2i}\\!\\!\\!-\\!\frac{s\mathrm{e}^{\frac{u_{2}^{2}}{4v_{2}}-\frac{u_{1}^{2}}{4v_{1}}-r}}{\sqrt{\pi}}\sum_{j=0}^{2i}q_{b}^{2i-j}I_{b}\right)\\!,$
(47)

where $\gamma_{b}=(k+1)/2$, $q_{b}=u_{1}/2v_{1}-u_{2}/2v_{2}$, and $I_{b}$ is
defined as

$\displaystyle
I_{b}=\frac{v_{2}^{(j+1)/2}}{2}\begin{cases}\Gamma\left(\frac{j+1}{2},\frac{u_{2}^{2}}{4v_{2}}\right),&\text{for}\,\,\frac{u_{2}}{2v_{2}}>0,\\\
(-1)^{j}\gamma\left(\frac{j+1}{2},\frac{u_{2}^{2}}{4v_{2}}\right)+\Gamma\left(\frac{j+1}{2}\right),&\text{otherwise}.\end{cases}$
(48)

Moreover, $B$, $s$, $r$, and $\\{u_{i},v_{i}\\}$ for $i\in\\{1,2\\}$ are
defined as follows:

$\displaystyle B$ $\displaystyle=$
$\displaystyle\frac{m_{u}^{m_{u}}\psi\mathrm{exp}\left(\frac{u_{1}^{2}}{4v_{1}}-\frac{\mu_{Y}^{2}}{2\sigma_{Y}^{2}}-\frac{1}{4a_{1}}\left(d\vartheta_{b}-\frac{\mu_{Y}}{\sigma_{Y}}\right)^{2}\right)}{2\Gamma(m_{u})\xi_{u}^{m_{u}}\sqrt{2\pi\sigma_{Y}^{2}a_{1}}},\quad\text{and}\quad
s=\frac{c\vartheta_{b}^{2}}{2\sqrt{a_{1}}},$ (49a) $\displaystyle r$
$\displaystyle=$
$\displaystyle\frac{1}{2\sqrt{a_{1}}}\left(d\vartheta_{b}-\frac{\mu_{Y}}{\sigma_{Y}^{2}}\right),\quad\text{and}\quad
u_{1}=d\vartheta_{b}-\frac{c\vartheta_{b}^{2}}{a1}\left(d\vartheta_{b}-\frac{\mu_{Y}}{\sigma_{Y}^{2}}\right),$
(49b) $\displaystyle u_{2}$ $\displaystyle=$ $\displaystyle
u_{1}+2sr,\quad\text{and}\quad
v_{1}=c\vartheta_{b}-\frac{c^{2}\vartheta_{b}^{2}}{a_{1}}+\frac{m_{u}}{\xi_{u}},\quad\text{and}\quad
v_{2}=v_{1}+s^{2},$ (49c)

where $a_{1}=c\vartheta_{b}^{2}+1/2\sigma_{Y}^{2}$, $c=0.374$, and $d=0.777$
[31].

###  IV-D Achievable Diversity Order

The diversity order is the negative slope of the average SER or outage
probability versus the average SNR curve in a log-log scale in the high SNR
regime [32], and it can be defined as

$\displaystyle
G_{d}=-\lim_{\gamma\rightarrow\infty}\frac{\mathrm{log}\left(\bar{P}_{e}\right)}{\mathrm{log}\left(\bar{\gamma}\right)}=-\lim_{\gamma\rightarrow\infty}\frac{\mathrm{log}\left(P_{out}\right)}{\mathrm{log}\left(\bar{\gamma}\right)}.$
(50)

The diversity order provides insights on how the outage probability or average
SER decays with the average SNR. Thus, we derive the asymptotic outage
probability and the average SER in the high SNR regime. The achievable
diversity order and array/coding gain are also quantified.

####  IV-D1 Asymptotic outage probability

The outage probability can be asymptotically approximated in the high SNR
regime as [32]

$\displaystyle\lim_{\bar{\gamma}\rightarrow\infty}{P}_{out}={P}^{\infty}_{out}\approx
O_{c}\left(\frac{\gamma_{th}}{\bar{\gamma}}\right)^{G_{d}}+\mathcal{O}\left(\bar{\gamma}^{-(G_{d}+1)}\right),$
(51)

where $O_{c}$ is a measure of the array/coding gain and $G_{d}$ is the
diversity order [32]. Then, an asymptotic approximation for $P_{out}$ can be
derived as (see Appendix F)

$\displaystyle{P}^{\infty}_{out}=\Omega\Phi(N,L)\left(\frac{\gamma_{th}}{\bar{\gamma}}\right)^{G_{d}}+\mathcal{O}\left(\bar{\gamma}^{-(G_{d}+1)}\right),$
(52)

where the diversity order ($G_{d}$) is given by

$\displaystyle G_{d}=NL\min{(m_{h},m_{g})}+m_{u}.$ (53)

Moreover, in (52), $\Omega$ and $\Phi(N,L)$ are defined as

$\displaystyle\Omega=\frac{2m_{u}^{m_{u}}\Gamma(2m_{u})}{2G_{d}\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2G_{d})},\qquad\text{and}\qquad\Phi(N,L)=\prod_{n=1}^{N}\prod_{l=1}^{L}\eta_{nl}\theta_{n},$
(54)

where $m_{l}=\max{(m_{h},m_{g})}$, $m_{s}=\min{(m_{h},m_{g})}$, and

$\displaystyle\theta_{n}$ $\displaystyle=$
$\displaystyle\alpha^{\prime}\pi\left(4\sqrt{\frac{m_{s}m_{l}}{\xi_{s_{n}}\xi_{l_{n}}}}\right)^{(m_{s}-m_{l})}\frac{\Gamma(2m_{s})\Gamma(2m_{s}-2m_{l})}{\Gamma(m_{l}-m_{s}+1/2)},$
(55a) $\displaystyle\alpha^{\prime}$ $\displaystyle=$
$\displaystyle\frac{4}{\Gamma(m_{s})\Gamma(m_{l})}\left(\frac{m_{s}m_{l}}{\xi_{s_{n}}\xi_{l_{n}}}\right)^{(m_{s}+m_{l})/2}.$
(55b)

Furthermore, $\xi_{s_{n}}$ and $\xi_{l_{n}}$ are the scaling parameters
related to $m_{s}$ and $m_{l}$, respectively, in the Nakagami-$m$ fading
channels. Finally, $O_{c}$ in (51) can be defined as

$\displaystyle O_{c}=\Omega\Phi(N,L).$ (56)

####  IV-D2 Asymptotic average SER

Similarly, the average SER can be approximated in the high SNR regime as [32]

$\displaystyle\lim_{\bar{\gamma}\rightarrow\infty}\bar{P}_{e}=\bar{P}^{\infty}_{e}\approx\left(G_{a}\bar{\gamma}\right)^{-G_{d}}+\mathcal{O}\left(\bar{\gamma}^{-(G_{d}+1)}\right),$
(57)

where $G_{d}$ and $G_{a}$ are diversity order and array gain, respectively. An
asymptotic approximation for the average SER (46) in high SNR regime can be
derived as (see Appendix F)

$\displaystyle\bar{P}^{\infty}_{e}$ $\displaystyle=$
$\displaystyle\left[(\Lambda\Phi(N,L))^{-\frac{1}{G_{d}}}\bar{\gamma}\right]^{-G_{d}}+\mathcal{O}\left(\bar{\gamma}^{-(G_{d}+1)}\right),$
(58)

where $G_{d}=m_{s}NL+m_{u}$ and $\Phi(N,L)$ is given in (54). Moreover,
$\Lambda$ is define as

$\displaystyle\Lambda$ $\displaystyle=$ $\displaystyle\frac{\omega
2^{G_{d}-1}m_{u}^{m_{u}}}{\sqrt{\pi}\vartheta^{G_{d}}G_{d}\xi_{u}^{m_{u}}}\frac{\Gamma(G_{d}+1/2)\Gamma(2m_{u})}{\Gamma(2G_{d})\Gamma(m_{u})}.$
(59)

Lastly, the array gain in (57) can be derived as

$\displaystyle G_{a}=(\Lambda\Phi(N,L))^{-{1}/{G_{d}}}.$ (60)

Remark 4: Our asymptotic performance analysis (53) reveals that the diversity
order is a linear function of the number of distributed IRSs, the number of
passive reflective elements in each IRS, and the $m$ parameters of the
Nakagami fading channels. It is worth noting that $S$ and $D$ are equipped
with a single antenna/RF-chain. However, as per analysis, we observe that the
distributed IRSs deployment can be used to recycle the existing EM waves in
the propagation environment without using additional active RF chains to
substantially increase the diversity order, and consequently, the system
reliability can be boosted in terms of lowering the outage probability and
average SER. It is worth noting that the current state-of-the-art wireless
systems, the diversity order improvements are obtained by virtue of increasing
active RF chains and antennas at the transmitter/receiver.

Remark 5: By varying $m$ parameter of Nakagami fading, a myriad of fading
environments can be investigated. For instance, higher the $m$ value, lesser
the fading severity. The diversity order of a single reflected path
corresponding to an IRS depends on the minimum value of Nakagami-$m$ parameter
of $S$-to-IRS and IRS-to-$D$ channels. Thus, more severe channel (or
equivalently, the channel with lesser $m$ value) always determines the
diversity order of a single IRS-aided reflected channel. However, when there
exist distributed IRSs and a direct channel, the overall diversity order is
the sum of the diversity orders of each end-to-end serviceable channel as per
(53). The proposed distributed IRSs-aided system achieves this diversity gain
by intelligently controlling the phase-shifts at the distributed passive
reflective elements such that the signals received via reflected channels are
constructively combined with the direct channel at $D$.

##  V Impact of Imperfect CSI

In this section, the impact of imperfectly estimated CSI is studied. In
practical wireless communication systems, the channel estimates are generally
perturbed by estimation errors [33]. To this end, the true channels in (1) and
(3) can be modeled to capture the detrimental effects of estimation errors as
follows [33]:

$\displaystyle v=\hat{v}+\rho\epsilon_{v},$ (61)

where $v\in\\{u,h_{nl},g_{nl}\\}$, $\hat{v}$ is the estimated channel of $v$,
and $\epsilon_{v}=v-\hat{v}$ is the channel estimation error, which is in turn
independent of $\hat{v}$ [34]. Here, $\rho$ accounts for the channel
estimation quality with $\rho=0$ for perfect CSI case and $\rho=1$ for the
worst-case CSI scenario, respectively. Thereby, the optimal SNR in (10) can be
rewritten for the imperfect CSI case as

$\displaystyle\hat{\gamma}^{*}=\bar{\gamma}\left[\left(\hat{\alpha}_{u}+\rho\varepsilon_{u}\right)+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\left(\hat{\alpha}_{g_{nl}}+\rho\varepsilon_{g_{nl}}\right)\left(\hat{\alpha}_{h_{nl}}+\rho\varepsilon_{h_{nl}}\right)\right]^{2}=\bar{\gamma}\left[\hat{R}+\rho
E\right]^{2},$ (62)

where $\hat{\alpha}_{v}$ and
$\varepsilon_{v}\sim\mathcal{CN}(0,\beta_{v}^{2})$ are the estimated channel
envelop of $v$ and the corresponding estimation error for
$v\in\\{u,h_{nl},g_{nl}\\}$, respectively. Moreover, $\hat{R}$ and $E$ are
defined as

$\displaystyle\hat{R}$ $\displaystyle=$
$\displaystyle\hat{\alpha}_{u}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\hat{\alpha}_{g_{nl}}\hat{\alpha}_{h_{nl}}=\hat{\alpha}_{u}+\hat{Y},$
(63a) $\displaystyle E$ $\displaystyle=$
$\displaystyle\varepsilon_{u}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\left(\rho\varepsilon_{g_{nl}}\varepsilon_{h_{nl}}+\hat{\alpha}_{h_{nl}}\varepsilon_{g_{nl}}+\hat{\alpha}_{g_{nl}}\varepsilon_{h_{nl}}\right),$
(63b)

where $\hat{Y}$ is an one-sided Gaussian distributed RV and the PDF is given
in (12). An achievable rate upper bound defined in (31b) can be derived for
imperfect CSI case as (see Appendix G)

$\displaystyle\mathcal{\hat{R}}_{ub}\approx\mathrm{log}_{2}\left(1+\mathbb{E}\\!\left[{\hat{\gamma}^{*}}\right]\right)=\mathrm{log}_{2}\left(1+\bar{\gamma}\left(\mathbb{E}\\!\left[{\hat{R}^{2}}\right]+\rho^{2}\mathbb{E}\\!\left[{E^{2}}\right]\right)\right),$
(64)

where $\mathbb{E}\\!\left[{\hat{R}^{2}}\right]$ and
$\mathbb{E}\\!\left[{E^{2}}\right]$ are given as

$\displaystyle\mathbb{E}\\!\left[{\hat{R}^{2}}\right]$ $\displaystyle=$
$\displaystyle\sigma_{u}^{2}+\sigma_{Y}^{2}+2\mu_{u}\mu_{Y}+\mu_{u}^{2}+\mu_{Y}^{2},$
(65a) $\displaystyle\mathbb{E}\\!\left[{E^{2}}\right]$ $\displaystyle=$
$\displaystyle\beta_{u}^{2}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}^{2}\left(\rho^{2}\beta_{g_{nl}}^{2}\beta_{h_{nl}}^{2}+{\xi}_{h_{nl}}\beta_{g_{nl}}^{2}+{\xi}_{g_{nl}}\beta_{h_{nl}}^{2}\right),$
(65b)

where ${\xi}_{v}$ is the scaling parameter of the channel $v$ for
$v\in\\{h_{nl},g_{nl}\\}$.

##  VI Numerical Results

In this section, we present our numerical results to investigate the
performance gains and to validate our analysis. In our simulations, the large-
scale fading is modeled as $\zeta_{ab}=\left(d_{0}/d_{ab}\right)^{\nu}\times
10^{\varphi_{ab}/10}$, where $a\in\\{S,n\\}$, $b\in\\{n,D\\}$, $d_{ab}$ is the
distance between nodes $a$ and $b$, $d_{0}=1$ m is a reference distance,
$\nu=2.8$ is the path-loss exponent, and $10^{\varphi_{ab}/10}$ captures log-
normal shadow fading with $\varphi_{ab}\sim(0,8)$ [35]. In the system
topology, $S$ and $D$ are in positioned at fixed locations and $1200\,$m apart
from each other. The IRSs are randomly distributed over an area of $1000\times
2000$ $\mathrm{m}^{2}$. Moreover, the amplitude of reflection coefficients
$\eta_{nl}$ for $n\in\\{1,\cdots,N\\}$ and $l\in\\{1,\cdots,L\\}$ is set to
$0.9$, and unless otherwise specified, the Nakagami-$m$ parameters ($m_{u}$,
$m_{h}$, and $m_{g}$) are set to 3.

Figure 4: The outage probability for different number of IRSs in the
distributed deployment denoted by $N\in\\{1,2,4,6,8\\}$, the number of
reflective elements per IRS is $L=32$, $m_{u}=m_{h}=m_{g}=3$, the threshold
SNR is $\gamma_{th}=0\,$dB, and $\beta_{v}^{2}=0\,$dB for
$v\in\\{u,h_{nl},g_{nl}\\}$.

In Fig. 4, we plot the outage probability as a function of the average
transmit SNR ($\bar{\gamma}$) by varying the number of distributed IRSs
defined by $N\in\\{1,2,4,6,8\\}$. The exact outage curves are generated by
using Monte-Carlo simulation, while the analytical outage curves are plotted
via our closed-form derivation in (28). As per Fig. 4, we observe that the
tightness of our outage probability approximation significantly improves with
$N$. Moreover, we observe that the distributed IRSs deployment outperforms the
single-IRS set-up. For example, the single-IRS deployment requires an average
transmit SNR of $1$ dB to achieve an outage probability of $10^{-4}$, and this
is an $104.2$% increase compared to the same transmit SNR requirement of the
distributed deployment with eight IRSs. Furthermore, the impact of imperfect
CSI on the outage probability is studied via using Monte-Carlo simulations,
and it severely affect the outage performance.

Figure 5: The outage probability for $N\in\\{1,2,4,8\\}$, $L=6$,
$m_{u}=m_{g}=3$, $m_{h}=1.5$, and $\gamma_{th}=10\,$dB.

In Fig. 5, we investigate asymptotic behavior of the outage probability in the
high SNR regime. To this end, our asymptotic outage probability (52) is
plotted as a function of the average transmit SNR, and Monte-Carlo simulations
are also provided for validation purposes. The asymptotic outage probability
curves reveal the achievable diversity orders for different distributed IRSs
deployments; $N\in\\{1,2,4,8\\}$. For instance, the diversity orders for
$(N=1,L=6)$ and $(N=2,L=6)$ cases can be obtained by computing the negative
gradients of asymptotic outage curves in Fig. 5 to be 12 and 21, respectively.
This observation clearly validates our diversity order analysis in (53), i.e.,
$G_{d}=NL\min(m_{h},m_{g})+m_{u}$ with $m_{h}=1.5$, $m_{g}=m_{u}=3$. Thus, the
results in Fig. 5 show that the diversity order increases with the number of
IRS reflective elements, and consequently, the system reliability can be
drastically improved by the proposed distributed IRSs-aided system without
adopting additional RF chains.

Figure 6: The average achievable rate for $N\in\\{1,2,4,6,8\\}$, $L=32$, and
$m_{u}=m_{h}=m_{g}=3$.

In Fig. 6, we plot the achievable average rate for different number of
distributed IRSs as per $N\in\\{1,2,4,6,8\\}$. The achievable rate lower/upper
bounds are plotted via (31a) and (31b), respectively. The accuracy of our rate
analysis is also validated via Monte-Carlo simulation of the exact achievable
rate. Fig. 6 clearly depicts that our lower/upper bounds are tight even for
smaller $N$. Fig. 6 reveals that higher the number of distributed IRSs, higher
the achievable rate. Moreover, Fig. 6 can be used to quantify the rate gains
of the distributed IRSs deployment compared to the direct transmission. For
instance, the single-IRS ($N=1$) case provides a rate gain of $80.7$% compared
to the direct transmission at an average transmit SNR of $-10$ dB. This rate
gain increases to $182.5$%, $266.2$%, and $405.9$% for dual-, quadruple-, and
octuple-IRS cases, respectively, with respect to the direct channel.

Figure 7: An average achievable rate comparison for different system
configurations. For all cases $L=32$ and $m_{u}=m_{h}=m_{g}=m$. The cases are
defined as follows: Case-1: direct channel only, Case-2/Case-3: IRSs are
located in closer with/without direct channel, Case-4/Case-6/Case-7: IRSs are
located far-away and with direct channel, and Case-5: IRSs are located far-
away and without direct channel, respectively.

In Fig. 7, the effects of the IRS positions and the shape parameter ($m$) of
the Nakagami fading are investigated. To this end, $S$ and $D$ are positioned
at fixed locations at $500$ m apart. The Case-1 plots the achievable rate of
the direct channel. For Case-2 and Case-3, the distances between $S$-IRS and
$D$-IRS are $750$ m. Moreover, the distances between $S$-IRS and $D$-IRS for
Case-4 to Case-7 are greater than $1600$ m. In Case-3 and Case-5, the
achievable rates are plotted by assuming that the direct channel is not
serviceable, and the end-to-end communication is established via the reflected
channels. Fig. 7 shows that the rate of the direct channel is higher than that
of the reflected channels when the IRSs are positioned far away from $S$ and
$D$ (see Case-1 and Case-5). When comparing Case-2 and Case-3, we can observe
that having the direct channel is beneficial in boosting the achievable rate.
Although the IRSs are placed far away in Case-4 than that in Case-3, the
achievable rate of the former is higher than that of the latter. This is
because Case-4 has two additional IRSs than the Case-3. Thus, the distributed
IRSs deployment can circumvent the rate losses incurred due to larger
transmission distances in the reflected channels. Moreover, Case-6 and Case-7
show the impact of different $m$ parameter of Nakagami fading. In Case-6,
$m=1$ leads to have Rayleigh fading channel, whereas in Case-7, $m=2$ depicts
the achievable rate of Rician fading with $1.5$ dB Rician factor.

Figure 8: The average achievable rate versus $\rho$ for $N\in\\{2,4,6,8\\}$,
$L=32$, $m_{u}=m_{h}=m_{g}=3$, and $\beta_{v}^{2}=0\,$dB for
$v\in\\{u,h_{nl},g_{nl}\\}$ .

In Fig. 8, the impact of imperfect CSI is studied by plotting the average
achievable rate versus $\rho$, which captures the amount of CSI imperfection.
Several rate curves are plotted for different number of distributed IRSs as
$N\in\\{2,4,6,8\\}$. The exact upper bounds of the average achievable rate
curves are plotted by using Monte-Carlo simulations, while the analytical
curves are plotted via our analysis in (64). As expected, the perfect CSI case
$(\rho=0)$ achieves the best rate performance, while the worst-case CSI
$(\rho=1)$ renders the lowest rates. For instance, the perfect CSI case
achieves rate gains of $18.2\%$ and $7.6\%$ compared to the worst-case CSI
case for $N=2$ and $N=8$, respectively.

Figure 9: The average BER of BPSK versus average transmit SNR ($\bar{\gamma}$)
for $N\in\\{1,2,4,6,8\\}$, $L=32$, $m_{u}=m_{h}=m_{g}=3$, $\omega=1$, and
$\vartheta=2$.

In Fig. 9, the average bit error rate (BER) of BPSK is plotted against the
average transmit SNR ($\bar{\gamma}$) for different number of IRSs;
$N\in\\{1,2,4,6,8\\}$. The exact average BER curves are generated via Monte-
Carlo simulations, where as the analytical curves are plotted by using (47).
From Fig. 9, we observe that our average BER analysis in (47) is tight for
moderately larger $NL$ regime. The average BER for a single-IRS set-up is also
plotted for the purpose of comparison. Fig. 9 reveals that the distributed
IRSs deployment outperforms the single-IRS set-up. For instance, to achieve an
average BER of $10^{-4}$, the single-IRS set-up requires an average transmit
SNR of $0$ dB, which is an $22$ dB increase compared to the transmit SNR
requirement of the distributed IRSs set-up with eight IRSs. Moreover, the
achievable average BER significantly decreases when the number of distributed
IRSs grows large. For example, at an average SNR of $-20$ dB, the average BER
of $N=4$ case is lowered by $99.8$% by doubling the number of IRSs to $N=8$.

Figure 10: A comparison of the average BER with/without IRSs for $N=2$,
$L\in\\{16,32,64\\}$, $m_{u}=m_{h}=m_{g}=3$, $\omega=1$, and $\vartheta=2$.

In Fig. 10, we investigate the effects of increasing number of IRS reflective
elements by plotting the average BER of BPSK versus the average SNR. A dual-
IRS system is considered and the number of reflective elements is varied as
$L\in\\{16,32,64\\}$. Analytical average BER is plotted via (47), and it is
validated through Monte-Carlo simulations. Fig. 10 shows that the dual-IRS
system outperforms the direct channel in terms of the average BER. For
example, the direct transmission requires an average SNR of $28$ dB to achieve
an average BER of $10^{-5}$, while this SNR requirement can be lowered by $36$
dB with a deployment of a dual-IRS set-up with $L=32$. Moreover, at an average
SNR of $-10$ dB, the average BER can be lowered from $1.82\times 10^{-2}$ to
$1.72\times 10^{-4}$ by doubling the number of reflective elements per IRS
from $L=16$ case to $L=32$ case. Thus, the proposed distributed IRSs
deployment can be used to significantly boost the reliability of wireless
communication systems.

##  VII Conclusion

In this paper, the fundamental performance metrics of a distributed IRSs-aided
communication system have been investigated. The optimal received SNR has been
derived by controlling the phase-shifts of the distributed IRS elements. To
facilitate a mathematically tractable performance analysis, this SNR has been
statistically characterized by deriving the tight approximations to the PDF
and CDF for Nakagami-$m$ fading. The KL divergence has been used to
analytically verify the accuracy of our statistical characterization of the
received SNR. Then, the outage probability, average SER, and achievable rate
approximations/bounds have been derived. Useful insights on the diversity
order and array gain have been drawn by deriving the asymptotic outage
probability and the average SER in the high SNR regime. By virtue of an
asymptotic rate analysis, it has been shown that the transmit power can be
scaled inversely proportional to the square of the number of reflective
elements, and in operating regime, our lower and upper rate bounds have been
shown to be asymptotically exact. A rigorous set of Monte-Carlo simulations
has been presented to validate the accuracy of our analysis. The performance
of the proposed distributed IRSs-assisted set-up has been compared with a
conventional/direct and a singe IRS-aided counterparts. Our analysis and
numerical results reveal that distributed IRSs deployment can be used to
significantly enhance the performance of future wireless communications by
recycling the existing EM waves via passive IRS reflective elements without
generating new EM waves with active RF chains. Thus, the concept of
distributed IRSs-aided communications can be an energy and cost effective
paradigm-shift in future wireless designs.

##  Appendix A The derivation of $\mu_{Y}$ in (13a), $\sigma_{Y}^{2}$ in
(13b), and PDF of $Y$ in (12)

First, we define ${Y}=\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}z_{nl}$, where
$z_{nl}=\alpha_{h_{nl}}\alpha_{g_{nl}}$. Since $\alpha_{h_{nl}}$ and
$\alpha_{g_{nl}}$ are independent Nakagami-$m$ variables, the PDF of $z_{nl}$
can be evaluated as

$\displaystyle f_{z_{nl}}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}f_{\alpha_{h_{nl}}}(u)f_{\alpha_{g_{nl}}}(x/u)\times\frac{1}{|u|}du=\alpha
x^{2m_{g}-1}\int_{0}^{\infty}u^{2(m_{h}-m_{g})-1}\mathrm{e}^{-\frac{m_{h}u^{2}}{\xi_{h_{n}}}-\frac{m_{g}x^{2}}{\xi_{g_{n}}u^{2}}}du$
(66) $\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}$
$\displaystyle\frac{\alpha}{2}x^{2m_{g}-1}\int_{0}^{\infty}t^{m_{h}-m_{g}-1}\mathrm{e}^{-\frac{m_{h}t}{\xi_{h_{n}}}-\frac{m_{g}x^{2}}{\xi_{g_{n}}t}}dt\stackrel{{\scriptstyle(b)}}{{=}}\alpha^{\prime}x^{m_{h}+m_{g}-1}\mathcal{K}_{m_{h}-m_{g}}\left(2x\sqrt{\frac{m_{h}m_{g}}{\xi_{h_{n}}\xi_{g_{n}}}}\right)\\!,$

where
$\alpha\\!=\\!4m_{h}^{m_{h}}m_{g}^{m_{g}}\\!/\\!\\!\left(\Gamma\\!(m_{h})\Gamma\\!(m_{g})\xi_{h_{n}}^{m_{h}}\xi_{g_{n}}^{m_{g}}\\!\right)$,
$\alpha^{\prime}\\!\\!=\\!\alpha\\!\left(\\!m_{g}\xi_{h_{n}}\\!/\\!m_{h}\xi_{g_{n}}\\!\right)^{\\!\frac{m_{h}\\!-\\!m_{g}}{2}}$,
and $\mathcal{K}_{v}(\cdot)$ is the Modified Bessel function of second kind
[25, Eq. (8.407)]. The step ($a$) is obtained by letting $t\\!=\\!u^{2}$ and
the step ($b$) is computed via [25, Eq. (3.471.9)]. Then, the $n$th moment of
$z_{nl}$ is derived as

$\displaystyle\mathbb{E}\\!\left[{z_{nl}^{n}}\right]$ $\displaystyle=$
$\displaystyle\alpha^{\prime}\int_{0}^{\infty}x^{m_{h}+m_{g}+n-1}\mathcal{K}_{m_{h}-m_{g}}\left(2x\sqrt{\frac{m_{h}m_{g}}{\xi_{h_{n}}\xi_{g_{n}}}}\right)dx$
(67) $\displaystyle\stackrel{{\scriptstyle(c)}}{{=}}$
$\displaystyle\frac{\alpha^{\prime}}{4}\left(\frac{m_{h}m_{g}}{\xi_{h_{n}}\xi_{h_{n}}}\right)^{\frac{-m_{h}-m_{g}-n+1}{2}}\Gamma\left(\frac{2m_{h}+n}{2}\right)\Gamma\left(\frac{2m_{g}+n}{2}\right),$

where the step ($c$) is written via [25, Eq. (6.561.16)]. The mean and
variance of $Y$ are given as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mu_{Y}$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{Y}\right]=\sum\limits_{n=1}^{N}\sum\limits_{l=1}^{L}\eta_{nl}\mathbb{E}\\!\left[{z_{nl}}\right]\quad\text{and}\quad\sigma_{Y}^{2}=\mathbb{V}\mathrm{ar}\\!\left[{Y}\right]=\sum\limits_{n=1}^{N}\sum\limits_{l=1}^{L}\eta_{nl}^{2}\left(\mathbb{E}\\!\left[{z_{nl}^{2}}\right]-\mathbb{E}\\!\left[{z_{nl}}\right]^{2}\right).$
(68)

By substituting (67) into (68), $\mu_{Y}$ and $\sigma_{Y}^{2}$ can be computed
as (13a) and (13b), respectively. Then, the PDF of $Y$ can be approximated by
an one-sided Gaussian distribution as given in (12) by invoking the central
limit theorem [24].

##  Appendix B The derivation of PDF of $\tilde{R}$ in (14)

Since $\alpha_{u}$ and $Y$ are independent RVs, the PDF of
$\tilde{R}=\alpha_{u}+\tilde{Y}$ can be derived as follows:

$\displaystyle f_{\tilde{R}}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}f_{u}(u)f_{\tilde{Y}}(x-u)du=2a^{m_{u}}\lambda\mathrm{e}^{-\frac{(x-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}\int_{0}^{\infty}u^{2m_{u}-1}\mathrm{e}^{-au^{2}+bu}du$
(69) $\displaystyle=$ $\displaystyle
2a^{m_{u}}\lambda\mathrm{e}^{-\frac{(x-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}\mathrm{e}^{\frac{b^{2}}{4a}}\underbrace{\int_{0}^{\infty}u^{2m_{u}-1}\mathrm{e}^{-a\left(u-\frac{b}{2\sqrt{a}}\right)^{2}}du}_{I_{\tilde{R}}},$

where $b\\!=\\!(x\\!-\\!\mu_{Y})\\!/\\!\sigma_{Y}^{2}$. Here, $a$ and
$\lambda$ are defined in (15). The solution to $I_{\tilde{R}}$ in (69) is
given by

$\displaystyle I_{\tilde{R}}$ $\displaystyle=$
$\displaystyle\frac{1}{a^{m_{u}}}\int_{\frac{-b}{2\sqrt{a}}}^{\infty}\left(t+\frac{b}{2\sqrt{a}}\right)^{2m_{u}-1}\\!\\!\\!\\!\\!\mathrm{e}^{-t^{2}}dt=\frac{1}{a^{m_{u}}}\sum_{k=0}^{2m_{u}-1}\binom{2m_{u}-1}{k}\left(\frac{b}{2\sqrt{a}}\right)^{2m_{u}-1-k}\int_{\frac{-b}{2\sqrt{a}}}^{\infty}t^{k}\mathrm{e}^{-t^{2}}dt$
(70) $\displaystyle\stackrel{{\scriptstyle(d)}}{{=}}$
$\displaystyle\frac{1}{2a^{m_{u}}}\sum_{k=0}^{2m_{u}-1}\binom{2m_{u}-1}{k}\left(\frac{b}{2\sqrt{a}}\right)^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},\frac{b^{2}}{4a}\right),$

where the step ($d$) is due to [25, Eq. (2.33.10)]. By substituting $(b)$ and
(70) into (69), the PDF of $\tilde{R}$ is derived as (14).

##  Appendix C The derivation of CDF of $\tilde{R}$ in (18)

By substituting (14) into (18), $I_{k}$ in (18) can be given as

$\displaystyle I_{k}$ $\displaystyle=$
$\displaystyle\int_{x}^{\infty}\mathrm{e}^{-\Delta\left(\frac{u-\mu_{Y}}{q}\right)^{2}}\left(\frac{u-\mu_{Y}}{q}\right)^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},\left(\frac{u-\mu_{Y}}{q}\right)^{2}\right)du$
(71) $\displaystyle\stackrel{{\scriptstyle(e)}}{{=}}$ $\displaystyle
q\int_{l_{min}}^{\infty}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt,$

where $q=2\sigma_{Y}^{2}\sqrt{a}$, $l_{min}=(x-\mu_{Y})/q$, and the step ($e$)
is written by changing the variable of integration. Then, we divide $I_{k}$
into two integrals: (i) $I_{k}^{o}$ for odd values of $k$ and (ii) $I_{k}^{e}$
for even values of $k$. Thus, for odd values of $k$ and when $x-\mu_{Y}<0$,
this integral can be written as

$\displaystyle I_{k}^{o}$ $\displaystyle=$ $\displaystyle
q\underbrace{\int_{l_{min}}^{0}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt}_{I_{k1}^{o}}+q\underbrace{\int_{0}^{\infty}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt}_{I_{k2}^{o}}.$
(72)

Then, the first integral $I_{k1}^{o}$ in (72) can be computed as

$\displaystyle I_{k1}^{o}$ $\displaystyle\stackrel{{\scriptstyle(f)}}{{=}}$
$\displaystyle\int_{0}^{-l_{min}}\\!\\!\\!\mathrm{e}^{-\Delta
y^{2}}y^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},y^{2}\right)dy\stackrel{{\scriptstyle(g)}}{{=}}(\gamma_{o}-1)!\sum_{i=0}^{\gamma_{o}-1}\frac{1}{i!}\int_{0}^{-l_{min}}\mathrm{e}^{-(\Delta+1)y^{2}}y^{2m_{u}-1-k}dy$
(73) $\displaystyle\stackrel{{\scriptstyle(h)}}{{=}}$
$\displaystyle\frac{(\gamma_{o}-1)!}{2}\sum_{i=0}^{\gamma_{o}-1}\frac{(\Delta+1)^{k/2-m_{u}-i}}{i!}\int_{0}^{l_{min}^{2}(\Delta+1)}v^{m_{u}+i-k/2}\mathrm{e}^{-v}dv$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}$
$\displaystyle\frac{(\gamma_{o}-1)!}{2}\sum_{i=0}^{\gamma_{o}-1}\frac{(\Delta+1)^{k/2-m_{u}-i}}{i!}\gamma\left(m_{u}+i-\frac{k}{2},(\Delta+1)l_{min}^{2}\right),$

where the step ($f$) and step ($h$) are obtained by substituting $t\\!=\\!-y$
and $v\\!=\\!(\Delta+1)y^{2}$, respectively. The step ($g$) is computed via
the fact that
$\Gamma(n,x)=(n-1)!\mathrm{e}^{-x}\sum\nolimits_{m=0}^{n-1}x^{m}/m!$ [25, Eq.
(8.352.7)]. The step ($i$) is written via [25, Eq. (8.350.1)]. The second
integral $I_{k2}^{o}$ in (72) is computed as

$\displaystyle I_{k2}^{o}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt\stackrel{{\scriptstyle(j)}}{{=}}(\gamma_{o}-1)!\sum_{i=0}^{\gamma_{o}-1}\frac{1}{i!}\int_{0}^{\infty}t^{2m_{u}+2i-1-k}\mathrm{e}^{-(\Delta+1)t^{2}}dt$
(74) $\displaystyle\stackrel{{\scriptstyle(k)}}{{=}}$
$\displaystyle\frac{(\gamma_{o}-1)!}{2}\sum_{i=0}^{\gamma_{o}-1}\frac{(\Delta+1)^{k/2-m_{u}-i}}{i!}\Gamma\left(m_{u}+i-\frac{k}{2}\right),$

where the step ($j$) is derived via a similar technique used in the step
($g$), and the step ($k$) is evaluated by using [25, Eq. (2.33.10)]. Then, for
$x-\mu_{Y}>0$, $I_{k}^{o}$ can be calculated as

$\displaystyle I_{k}^{o}$ $\displaystyle=$
$\displaystyle\int_{l_{min}}^{\infty}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt=(\gamma_{o}-1)!\sum_{i=0}^{\gamma_{o}-1}\frac{1}{i!}\int_{l_{min}}^{\infty}t^{2m_{u}+2i-1-k}\mathrm{e}^{-(\Delta+1)t^{2}}dt$
(75) $\displaystyle=$
$\displaystyle\frac{(\gamma_{o}-1)!}{2}\sum_{i=0}^{\gamma_{o}-1}\frac{(\Delta+1)^{k/2-m_{u}-i}}{i!}\Gamma\left(m_{u}+i-\frac{k}{2},(\Delta+1)l_{min}^{2}\right).$

By combining (72) and (75), $I_{k}^{o}$ can be derived as (20). Then, for even
values of $k$, the integral in (71) can be evaluated as

$\displaystyle I_{k}^{e}$ $\displaystyle=$ $\displaystyle
q\int_{l_{min}}^{\infty}\mathrm{e}^{-\Delta
t^{2}}t^{2m_{u}-1-k}\Gamma\left(\frac{k+1}{2},t^{2}\right)dt$
$\displaystyle\stackrel{{\scriptstyle(l)}}{{=}}$
$\displaystyle\frac{(\gamma_{e}-1)!}{2}\sum_{j=0}^{\gamma_{e}-1}\left[-\frac{\mathrm{e}^{-\Delta
t^{2}}t^{2j}}{j!\Delta^{\gamma_{e}-1}}\Gamma\left(\frac{k+1}{2},t^{2}\right)\right]_{l_{min}}^{\infty}\\!\\!-(\gamma_{e}-1)!\sum_{j=0}^{\gamma_{e}-1}\frac{\Delta^{\gamma_{e}-1}}{j!}\int_{l_{min}}^{\infty}\\!\\!t^{k+2j}\mathrm{e}^{-(\Delta-1)t^{2}}dt$
$\displaystyle=$
$\displaystyle\frac{(\gamma_{e}\\!-\\!1)!}{2}\sum_{j=0}^{\gamma_{e}-1}\\!\frac{\mathrm{e}^{-\Delta
l_{min}^{2}}l_{min}^{2j}}{j!\Delta^{\gamma_{e}-1}}\Gamma\\!\left(\frac{k\\!+\\!1}{2},l_{min}^{2}\right)\\!-\\!\frac{(\gamma_{e}\\!-\\!1)!}{2}\sum_{j=0}^{\gamma_{e}-1}\\!\frac{\Delta^{\gamma_{e}-1}}{j!}\frac{\Gamma\left(j+\frac{k}{2}+\frac{1}{2},(\Delta+1)l_{min}^{2}\right)}{(\Delta+1)^{j+k/2+1/2}},$

where the step ($l$) is obtained by using the part-by-part integration
technique.

## Appendix D

###  D-A The Derivation of Expectation of $\tilde{\gamma}^{*}$ in (32)

The expectation term in (32) can be computed as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\tilde{\gamma}^{*}}\right]\\!$
$\displaystyle=$
$\displaystyle\\!\mathbb{E}\\!\left[{\bar{\gamma}\tilde{R}^{2}}\right]\\!\\!=\\!\bar{\gamma}\mathbb{E}\\!\left[{\\!(\alpha_{u}\\!+\\!\tilde{Y})^{2}\\!}\right]\\!\\!=\\!\bar{\gamma}\\!\sum_{n=0}^{2}\\!\\!\binom{2}{n}\mathbb{E}\\!\left[{\alpha_{u}^{(2-n)}}\right]\\!\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]\\!\\!=\\!\bar{\gamma}\\!\left(\\!\sigma_{u}^{2}\\!+\\!\mu_{u}^{2}\\!+\\!\sigma_{Y}^{2}\\!+\\!\mu_{Y}^{2}\\!+\\!2\mu_{u}\mu_{Y}\\!\right)\\!\\!,$
(77)

where $\mu_{Y}$, $\sigma_{Y}^{2}$, $\mu_{u}$, and $\sigma_{u}^{2}$ are given
in (13a), (13b), (33a), and (33b), receptively.

###  D-B The Derivation of Variance of $\tilde{\gamma}^{*}$ in (35)

The first expectation term in (35) can be evaluated as

$\displaystyle\mathbb{E}\\!\left[{\tilde{R}^{4}}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{(\alpha_{u}+\tilde{Y})^{4}}\right]=\sum_{n=0}^{4}\binom{4}{n}\mathbb{E}\\!\left[{\alpha_{u}^{(4-n)}}\right]\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right].$
(78)

Thus, the $n$th moment of $\alpha_{u}^{n}$ denoted by
$\mathbb{E}\\!\left[{\alpha_{u}^{n}}\right]$ can be computed as

$\displaystyle\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\alpha_{u}^{n}}\right]$
$\displaystyle=$
$\displaystyle\\!\int_{0}^{\infty}\\!\\!\\!x^{n}f_{u}(x)dx=\frac{2m_{u}^{m_{u}}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}\int_{0}^{\infty}\\!\\!\\!x^{2m_{u}+n-1}\mathrm{e}^{-\frac{m_{u}x^{2}}{\xi_{u}}}dx\stackrel{{\scriptstyle(m)}}{{=}}\frac{\xi_{u}^{n/2}}{\Gamma(m_{u})m_{u}^{n/2}}\Gamma\left(m_{u}+n/2\right)\\!,$
(79)

where the step ($m$) is evaluated by using [25, Eq. (2.33.10)]. Then,
$\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]$ can be derived as

$\displaystyle\mathbb{E}\\!\left[{\tilde{Y}^{n}}\right]$ $\displaystyle=$
$\displaystyle\frac{\psi}{\sqrt{2\pi\sigma_{Y}^{2}}}\int_{0}^{\infty}y^{n}\mathrm{e}^{-\frac{(y-\mu_{Y})^{2}}{2\sigma_{Y}^{2}}}dy\stackrel{{\scriptstyle(n)}}{{=}}\frac{\psi}{\sqrt{\pi}}\int_{\frac{-\mu_{Y}}{\sqrt{2\sigma_{Y}^{2}}}}^{\infty}\left(\sqrt{2\sigma_{Y}^{2}}t+\mu_{Y}\right)^{n}\mathrm{e}^{-t^{2}}dt$
(80) $\displaystyle\stackrel{{\scriptstyle(o)}}{{=}}$
$\displaystyle\frac{\psi}{2\sqrt{\pi}}\sum_{i=0}^{n}\binom{n}{i}\left(\sqrt{2\sigma_{Y}^{2}}\right)^{n-i}\mu_{Y}^{i}I\\!\left(n-i,\frac{-\mu_{Y}}{2\sigma_{Y}^{2}}\right),$

where the step ($n$) is due to a changing of variable and the step ($o$) is
obtained by expanding $\left(\sqrt{2\sigma_{Y}^{2}}t+\mu_{Y}\right)^{n}$ based
on $n$ value. Moreover, $I(\cdot,\cdot)$ is given in (39).

###  D-C The Derivation of the Asymptotic Achievable Rate
$(\mathcal{R}^{\infty})$ as $L\rightarrow\infty$ in (43)

When the number of IRS elements grows without bound ($L\rightarrow\infty$),
$\mu_{Y}$ in (13a) and $\sigma_{Y}^{2}$ in (13b) can be approximated as

$\displaystyle\mu_{Y}$ $\displaystyle\approx$ $\displaystyle
L\bar{\mu}_{Y}=L\left(\sum_{n=1}^{N}\eta_{n}\sqrt{\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\frac{\Gamma\left(m_{h}+1/2\right)\Gamma\left(m_{g}+1/2\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)}\right),$
(81) $\displaystyle\sigma_{Y}^{2}$ $\displaystyle\approx$ $\displaystyle
L\bar{\sigma}_{Y}^{2}=L\left(\sum_{n=1}^{N}\eta_{n}^{2}\left({\frac{\xi_{h_{n}}\xi_{g_{n}}}{m_{h}m_{g}}}\right)\frac{\Gamma\left(m_{h}+1\right)\Gamma\left(m_{g}+1\right)}{\Gamma\left(m_{h}\right)\Gamma\left(m_{g}\right)}-\bar{\mu}_{Y}^{2}\right),$
(82)

where $\eta_{n}=\eta_{nl}$. Thus, the SNR term in $\mathcal{R}_{ub}$ (41) can
be asymptotically evaluated as

$\displaystyle\lim_{L\rightarrow\infty}\gamma_{ub}$ $\displaystyle=$
$\displaystyle\lim_{L\rightarrow\infty}\bar{\gamma}L^{2}\left(\frac{\sigma_{u}^{2}}{L^{2}}+\frac{\bar{\sigma}_{Y}^{2}}{L}+\frac{2\mu_{u}\bar{\mu}_{Y}}{L}+\frac{\mu_{u}^{2}}{L^{2}}+\bar{\mu}_{Y}^{2}\right)=\bar{\gamma}_{E}\bar{\mu}_{Y}^{2},$
(83)

where $\bar{\gamma}_{E}=P_{E}/\sigma^{2}_{w}$ and
$\lim_{L\rightarrow\infty}P=P_{E}/L^{2}$. Similarly, in the asymptotic regime,
the SNR term in $\mathcal{R}_{lb}$ (40) can be derived as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\lim_{L\rightarrow\infty}\\!\gamma_{lb}$
$\displaystyle=$
$\displaystyle\\!\lim_{L\rightarrow\infty}\bar{\gamma}\\!\\!\left(\\!\\!\frac{L^{6}\left(\frac{\sigma_{u}^{2}}{L^{2}}+\frac{\bar{\sigma}_{Y}^{2}}{L}+\frac{2\mu_{u}\bar{\mu}_{Y}}{L}+\frac{\mu_{u}^{2}}{L^{2}}+\bar{\mu}_{Y}^{2}\right)^{\\!\\!3}}{L^{4}\left(\frac{\psi}{2\sqrt{\pi}}\bar{\mu}_{Y}^{4}I\left(0,\frac{-\sqrt{L}\bar{\mu}_{Y}}{\sqrt{2\bar{\sigma}_{Y}^{2}}}\right)+\frac{\Upsilon(L^{3})}{L^{4}}\right)}\\!\\!\right)\\!\\!=\\!\\!\lim_{L\rightarrow\infty}\bar{\gamma}L^{2}\\!\\!\left(\\!\frac{\bar{\mu}_{Y}^{2}}{\frac{\psi}{2\sqrt{\pi}}I\left(0,-\infty\right)}\\!\\!\right)\\!\\!\stackrel{{\scriptstyle(p)}}{{=}}\\!\\!\bar{\gamma}_{E}\bar{\mu}_{Y}^{2},$
(84)

where $\Upsilon(L^{3})$ is lower order ($L<4$) terms of the expansion of the
SNR in (40). Moreover, the step ($p$) is written by using the fact that
$\lim_{L\rightarrow\infty}\psi=\lim_{L\rightarrow\infty}1/\mathcal{Q}\left(-\sqrt{L}\bar{\mu}_{Y}/\bar{\sigma}_{Y}\right)=\lim_{L\rightarrow\infty}1/\mathcal{Q}\left(-\infty\right)=1$
and $I\left(0,-\infty\right)=2\Gamma(1/2)=2\sqrt{\pi}$. Then, the asymptotic
achievable rate $\mathcal{R}^{\infty}$ can be derived as given in (43).

##  Appendix E The derivation of the average SER in (47)

The inner expectation term with respect to $\tilde{Y}$ in (46) can be
evaluated via a tight approximation [31] for the Gaussian $Q$-function:
$\mathcal{Q}(x)\\!\approx\\!\mathrm{exp}\left(-cx^{2}-dx\right)/2$, where
$c\\!=\\!0.374$ and $d\\!=\\!0.777$ as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}_{\tilde{Y}}\\!\left[{\mathcal{Q}\left(\vartheta_{b}(\alpha_{u}+\tilde{Y})\right)}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{0}^{\infty}\mathrm{e}^{-cx^{2}-dx}f_{\tilde{Y}}(y)dy=\frac{\psi\mathrm{e}^{-c\vartheta_{b}^{2}\alpha_{u}^{2}-d\vartheta_{b}\alpha_{u}-\frac{\mu_{Y}^{2}}{2\sigma_{Y}^{2}}}}{2\sqrt{2\pi\sigma_{Y}^{2}}}\underbrace{\int_{0}^{\infty}\mathrm{e}^{-a_{1}y^{2}-a_{2}y}dy}_{I_{\tilde{Y}}},$
(85)

where $a_{1}=c\vartheta_{b}^{2}+1/2\sigma_{Y}^{2}$ and
$a_{2}=2c\vartheta_{b}^{2}\alpha_{u}+d\vartheta_{b}-\mu_{Y}/\sigma_{Y}^{2}$.
Then, the integral $I_{\tilde{Y}}$ in (85) can be evaluated by first using a
substitution $t=y+a_{2}/(2a_{1})$ and then invoking [25, Eq. (2.33.16)]

$\displaystyle I_{\tilde{Y}}$ $\displaystyle=$
$\displaystyle\mathrm{exp}\left(\frac{a_{2}^{2}}{4a_{1}^{2}}\right)\int_{0}^{\infty}\mathrm{exp}\left(-a_{1}\left(y+\frac{a_{2}}{2a_{1}}\right)^{2}\right)dy=\mathrm{exp}\left(\frac{a_{2}^{2}}{4a_{1}^{2}}\right)\int_{a_{2}/2a_{1}}^{\infty}\mathrm{exp}\left(-a_{1}t^{2}\right)dt$
(86) $\displaystyle=$
$\displaystyle\frac{1}{2}\sqrt{\frac{\pi}{a_{1}}}\mathrm{exp}\left(\frac{a_{2}^{2}}{4a_{1}^{2}}\right)\mathrm{erfc}\left(\frac{a_{2}}{2\sqrt{a_{1}}}\right).$

Then, by substituting $a_{2}$ and (86) into (85), the inner expectation in
(46) can be rewritten as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}_{\tilde{Y}}\\!\left[{\mathcal{Q}\left(\vartheta_{b}(\alpha_{u}\\!+\\!\tilde{Y})\right)}\right]\\!$
$\displaystyle=$
$\displaystyle\\!\frac{\psi\mathrm{exp}\left(\\!-c\vartheta_{b}^{2}\alpha_{u}^{2}-d\vartheta_{b}\alpha_{u}-\frac{\mu_{Y}^{2}}{2\sigma_{Y}^{2}}\right)}{4\sqrt{2\sigma_{Y}^{2}a_{1}}}\mathrm{exp}\left(-v^{\prime}\alpha_{u}^{2}-u_{1}\alpha_{u}\right)\mathrm{erfc}\left(s\alpha_{u}\\!+r\right)\\!,$
(87)

where $v^{\prime}=v_{1}-m_{u}/P_{u}$. Moreover, $s$, $r$, $u_{1}$, and $v_{1}$
are defined in (49). By substituting (87) into (46), the outer expectation
with respect to $\alpha_{u}$ can be written as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}_{\alpha_{u}}\\!\left[{\mathbb{E}_{\tilde{Y}}\\!\left[{\mathcal{Q}\left(\vartheta_{b}(\alpha_{u}+\tilde{Y})\right)}\right]}\right]$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}\mathbb{E}_{\tilde{Y}}\\!\left[{\mathcal{Q}\left(\vartheta_{b}(\alpha_{u}+\tilde{Y})\right)}\right]f_{u}(\alpha_{u})d\alpha_{u}$
(88) $\displaystyle=$ $\displaystyle
2\sqrt{2\pi}B\underbrace{\int_{0}^{\infty}\alpha_{u}^{2m_{u}-1}\mathrm{e}^{-v_{1}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2}}\mathcal{Q}\left(\sqrt{2}(s\alpha_{u}+r)\right)d\alpha_{u}}_{I_{\alpha_{u}}},$

where $B$ is defined in (49). Then, the integral $I_{\alpha_{u}}$ in (88) can
be evaluated as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!I_{\alpha_{u}}=\int_{0}^{\infty}\alpha_{u}^{2m_{u}-1}\mathrm{e}^{-v_{1}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2}}\mathcal{Q}\left(\sqrt{2}(s\alpha_{u}+r)\right)d\alpha_{u}$
$\displaystyle\stackrel{{\scriptstyle(q)}}{{=}}$
$\displaystyle\frac{-1}{2}\sum_{k=0}^{2m_{u}-1}\\!\\!\left(\frac{-u_{1}}{2v_{1}}\right)^{\\!\\!2m_{u}-1-k}\Gamma(\gamma_{b})\\!\left[\mathcal{Q}\left(\sqrt{2}(s\alpha_{u}+r)\right)\mathrm{e}^{-v_{1}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2}}\sum_{i=0}^{\gamma_{b}-1}\frac{\left(\alpha_{u}+u_{1}/2v_{1}\right)^{2i}}{i!v_{1}^{\gamma_{b}-i}}\right]_{0}^{\infty}$
$\displaystyle-\frac{s}{2\sqrt{\pi}}\\!\\!\sum_{k=0}^{2m_{u}-1}\\!\\!\left(\frac{-u_{1}}{2v_{1}}\right)^{\\!2m_{u}-1-k}\\!\\!\Gamma(\gamma_{b})\sum_{i=0}^{\gamma_{b}-1}\frac{1}{i!v_{1}^{\gamma_{b}-i}}\int_{0}^{\infty}\\!\\!\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{\\!\\!2i}\\!\\!\mathrm{e}^{-v_{1}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2}-(s\alpha_{u}+r)^{2}}d\alpha_{u}$
$\displaystyle=$
$\displaystyle\frac{\mathcal{Q}\left(\sqrt{2}r\right)}{2}\sum_{k=0}^{2m_{u}-1}\left(\frac{-u_{1}}{2v_{1}}\right)^{2m_{u}-1-k}\Gamma(\gamma_{b})\sum_{i=0}^{\gamma_{b}-1}\frac{\left(u_{1}/2v_{1}\right)^{2i}}{i!v_{1}^{\gamma_{b}-i}}$
$\displaystyle-\frac{s}{2\sqrt{\pi}}\sum_{k=0}^{2m_{u}-1}\left(\frac{-u_{1}}{2v_{1}}\right)^{\\!\\!2m_{u}-1-k}\\!\\!\Gamma(\gamma_{b})\\!\sum_{i=0}^{\gamma_{b}-1}\frac{1}{i!v_{1}^{\gamma_{b}-i}}\underbrace{\int_{0}^{\infty}\\!\\!\left(\\!\alpha_{u}\\!+\\!\frac{u_{1}}{2v_{1}}\\!\right)^{\\!2i}\mathrm{e}^{-v_{1}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2}-(s\alpha_{u}+r)^{2}}d\alpha_{u}}_{I_{\alpha_{u}^{\prime}}},$

where $\gamma_{b}=(k+1)/2$. The step ($q$) is written by using the part-by-
part integration technique. Then, the integral $I_{\alpha_{u}^{\prime}}$ can
be rearranged as follows:

$\displaystyle I_{\alpha_{u}^{\prime}}$ $\displaystyle=$
$\displaystyle\mathrm{e}^{\frac{u_{2}^{2}}{4v_{2}}-\frac{u_{1}^{2}}{4v_{1}}-r}\int_{0}^{\infty}\left(\alpha_{u}+\frac{u_{1}}{2v_{1}}\right)^{2i}\mathrm{e}^{-v_{2}\left(\alpha_{u}+\frac{u_{2}}{2v_{2}}\right)^{2}}d\alpha_{u}$
(90) $\displaystyle\stackrel{{\scriptstyle(r)}}{{=}}$
$\displaystyle\mathrm{e}^{\frac{u_{2}^{2}}{4v_{2}}-\frac{u_{1}^{2}}{4v_{1}}-r}\sum_{j=0}^{2i}q_{b}^{2i-j}\int_{\frac{u_{2}}{2v_{2}}}^{\infty}t^{j}\mathrm{e}^{-v_{2}t^{2}}dt=\mathrm{e}^{\frac{u_{2}^{2}}{4v_{2}}-\frac{u_{1}^{2}}{4v_{1}}-r}\sum_{j=0}^{2i}q_{b}^{2i-j}I_{b},$

where $q_{b}=u_{1}/2v_{1}-u_{2}/2v_{2}$ and $I_{b}$ is defined in (48).
Further, the step ($r$) is obtained via the substitution
$t=\alpha_{u}+u_{2}/(2v_{2})$. Finally, $u_{2}$ and $v_{2}$ are defined in
(49).

##  Appendix F The derivation of ${P}^{\infty}_{out}$ in (52) and
$\bar{P}_{e}^{\infty}$ in (58)

From (66) in Appendix A, the PDF of $z_{nl}=\alpha_{h_{nl}}\alpha_{g_{nl}}$
can be written as

$\displaystyle f_{z_{nl}}(x)$ $\displaystyle=$
$\displaystyle\alpha^{\prime}x^{m_{s}+m_{l}-1}\mathcal{K}_{m_{s}-m_{l}}\left(2x\sqrt{{m_{s}m_{l}}/{(\xi_{s_{n}}\xi_{l_{n}})}}\right),$
(91)

where $\alpha^{\prime}$ is given in (55b), $m_{s}=\min{(m_{h},m_{g})}$, and
$m_{l}=\max{(m_{h},m_{g})}$. Moreover, $\xi_{s_{n}}$ and $\xi_{l_{n}}$ are the
scaling parameters of the respective channel models (3). The moment generating
function (MGF) of $z_{nl}$ can be derived by evaluating the Laplace transform
of $f_{z_{nl}}(x)$ as [24]

$\displaystyle\mathcal{M}_{z_{nl}}(s)$ $\displaystyle=$
$\displaystyle\\!\int_{0}^{\infty}\\!\\!\\!f_{z_{nl}}(x)\mathrm{exp}(-sx)dx=\alpha^{\prime}\\!\\!\int_{0}^{\infty}\\!\\!\\!x^{m_{s}+m_{l}-1}\mathrm{e}^{-sx}\mathcal{K}_{m_{s}-m_{l}}\left(2x\sqrt{{m_{s}m_{l}}/{(\xi_{s_{n}}\xi_{l_{n}})}}\right)dx$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\stackrel{{\scriptstyle(s)}}{{=}}\\!\frac{\alpha^{\prime}\pi\left(2b_{n}\right)^{(m_{s}-m_{l})}\Gamma(2m_{s})\Gamma(2m_{l})}{\Gamma(m_{s}+m_{l}+1/2)}\left(s\\!+\\!b_{n}\right)^{2m_{s}}F\left(\\!2m_{s},m_{s}\\!-\\!m_{l}\\!+\\!\frac{1}{2};m_{s}\\!+\\!m_{2}\\!+\\!\frac{1}{2};\frac{s\\!-\\!b_{n}}{s\\!+\\!b_{n}}\right),$

where $b_{n}=2\sqrt{m_{s}m_{l}/(\xi_{h_{n}}\xi_{g_{n}})}$ and
$F(\cdot,\cdot;\cdot;\cdot)$ is the Gauss Hypergeometric function [25, Eq.
(9.100)]. In (F), the step $(s)$ is derived by using [25, Eq. (6.621.3)]. The
behavior of the PDF of $z_{nl}$ at the origin is governed by the asymptotic
value of $\mathcal{M}_{z_{nl}}(s)$ as $s\rightarrow\infty$ [32]. To this end,
by using the fact that $\lim_{s\rightarrow\infty}$
$({s-b_{n}})/({s+b_{n}})\rightarrow 1$, and since $m_{s}<m_{l}$, we invoke
[25, Eq. (9.122.1)] to evaluate $\mathcal{M}_{z_{nl}}(s)$ when
$s\rightarrow\infty$ as follows:

$\displaystyle\mathcal{M}_{z_{nl}}^{\infty}(s)$ $\displaystyle=$
$\displaystyle\theta_{n}s^{-2m_{s}},$ (93)

where $\theta_{n}$ is defined in (55a). Thus, the MGF of
${Y}=\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}z_{nl}$ as $s\rightarrow\infty$ is
given by

$\displaystyle\mathcal{M}_{Y}^{\infty}(s)$ $\displaystyle=$
$\displaystyle\prod_{n=1}^{N}\prod_{l=1}^{L}\eta_{nl}\theta_{n}s^{-2m_{s}}=\Phi(N,L)s^{-2m_{s}NL},$
(94)

where $\Phi(N,L)$ is given in (54). The PDF of the direct channel envelope
$(u)$ is expanded at the origin via the Maclaurin series expansion of the
exponential function [25, Eq. (0.318.2)] as

$\displaystyle f_{\alpha_{u}}(x)$ $\displaystyle=$
$\displaystyle\frac{2m_{u}^{m_{u}}x^{2m_{u}-1}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}\mathrm{exp}\left(\frac{-x^{2}}{\xi_{u}}\right)=\frac{2m_{u}^{m_{u}}x^{2m_{u}-1}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}\sum_{k=0}^{\infty}\frac{1}{k!}\left(\frac{-x^{2}}{\xi_{u}}\right)^{k}$
(95) $\displaystyle=$
$\displaystyle\frac{2m_{u}^{m_{u}}x^{2m_{u}-1}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}\left(1-\frac{x^{2}}{\xi_{u}}+\frac{x^{4}}{2\xi_{u}^{2}}-\frac{x^{6}}{6\xi_{u}^{3}}+\cdots\right).$

Then, the PDF of the direct channel envelope $(u)$ for $x\rightarrow 0^{+}$
can be approximated as

$\displaystyle f_{\alpha_{u}}^{0^{+}}(x)$ $\displaystyle=$
$\displaystyle\frac{2m_{u}^{m_{u}}x^{2m_{u}-1}}{\Gamma(m_{u})\xi_{u}^{m_{u}}}+\mathcal{O}\left(x^{2m_{u}}\right).$
(96)

Via the Laplace transform of (96), we obtain the asymptotic MGF of
$\alpha_{u}$ as [25, Eq. (17.13.2)]

$\displaystyle\mathcal{M}_{\alpha_{u}}^{\infty}(s)$ $\displaystyle=$
$\displaystyle\frac{2m_{u}^{m_{u}}\Gamma(2m_{u})}{\Gamma(m_{u})\xi_{u}^{m_{u}}}s^{-2m_{u}}.$
(97)

Then, the asymptotic MGF of $R=\alpha_{u}+Y$ for $s\rightarrow\infty$ can be
derived as

$\displaystyle\mathcal{M}_{R}^{\infty}(s)$ $\displaystyle=$
$\displaystyle\mathcal{M}_{z_{nl}}(s)\mathcal{M}_{\alpha_{u}}^{\infty}(s)=\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{\Gamma(m_{u})\xi_{u}^{m_{u}}}s^{-2m_{s}NL-2m_{u}},$
(98)

where $\Phi(N,L)$ is defined in (54). By taking inverse Laplace transform of
(98) [25], the PDF of $R$ for $x\rightarrow 0^{+}$ can be computed as

$\displaystyle f_{R}^{0^{+}}(x)$
$\displaystyle\stackrel{{\scriptstyle(t)}}{{=}}$
$\displaystyle\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2m_{s}NL+2m_{u})}x^{2m_{s}NL+2m_{u}-1}+\mathcal{O}\left(x^{2m_{s}NL+2m_{u}}\right),$
(99)

where the step $(t)$ is written via [25, Eq. (17.13.3)]. Next, the CDF of $R$
for $x\\!\rightarrow\\!0^{+}$ is given as

$\displaystyle\\!\\!F_{R}^{0^{+}}(x)$ $\displaystyle=$
$\displaystyle\int_{0}^{x}\\!\\!\\!f_{R}^{0+}(u)du\\!=\\!\int_{0}^{x}\\!\\!\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2m_{s}NL+2m_{u})}u^{2m_{s}NL+2m_{u}-1}du+\mathcal{O}\left(u^{2m_{s}NL+2m_{u}}\right)$
(100) $\displaystyle=$
$\displaystyle\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{2G_{d}\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2G_{d})}x^{2G_{d}}+\mathcal{O}\left(x^{2G_{d}+1}\right),$

where $G_{d}=m_{s}NL+m_{u}$. Then, the CDF of $\gamma^{*}=\bar{\gamma}R^{2}$
for $x\rightarrow 0^{+}$ is approximated as [24]

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!F_{\gamma^{*}}^{0^{+}}(x)$
$\displaystyle=$ $\displaystyle\mathrm{Pr}\left(\gamma^{*}\leq x\right)\approx
F_{\tilde{R}}^{0^{+}}\left(\sqrt{x/\bar{\gamma}}\right)=\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{2G_{d}\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2G_{d})}\left(\frac{x}{\bar{\gamma}}\right)^{G_{d}}\\!\\!\\!+\mathcal{O}\\!\left(\\!\left(\frac{x}{\bar{\gamma}}\right)^{\\!\\!G_{d}+1}\right)\\!.$
(101)

Thereby, the asymptotic outage probability can be computed as

$\displaystyle{P}^{\infty}_{out}\approx
F_{\gamma^{*}}^{0^{+}}(\gamma_{th})=\Omega\Phi(N,L)\left(\frac{\gamma_{th}}{\bar{\gamma}}\right)^{G_{d}}+\mathcal{O}\left(\left(\frac{\gamma_{th}}{\bar{\gamma}}\right)^{G_{d}+1}\right),$
(102)

where $\Omega$ is given in (54). Finally, the asymptotic average BER can be
derived as follows:

$\displaystyle P_{e}^{\infty}$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}\omega\mathcal{Q}\left(\sqrt{\vartheta
x}\right)f_{\gamma^{*}}^{0^{+}}(x)dx=\frac{\omega}{2}\sqrt{\frac{\vartheta}{2\pi}}\int_{0}^{\infty}x^{-1/2}\mathrm{exp}\left(\frac{-\vartheta
x}{2}\right)F_{\gamma^{*}}^{0^{+}}(x)dx$ (103) $\displaystyle=$
$\displaystyle\frac{\omega}{2}\sqrt{\frac{\vartheta}{2\pi}}\frac{\Phi(N,L)2m_{u}^{m_{u}}\Gamma(2m_{u})}{2G_{d}\Gamma(m_{u})\xi_{u}^{m_{u}}\Gamma(2G_{d})}\bar{\gamma}^{-G_{d}}\int_{0}^{\infty}x^{G_{d}-1/2}\mathrm{exp}\left(\frac{-\vartheta
x}{2}\right)dx+\mathcal{O}\\!\left(\\!\left(\frac{x}{\bar{\gamma}}\right)^{G_{d}+1}\right)$
$\displaystyle\stackrel{{\scriptstyle(u)}}{{=}}$
$\displaystyle\Lambda\Phi(N,L)\bar{\gamma}^{-G_{d}}+\mathcal{O}\left(\bar{\gamma}^{-(G_{d}+1)}\right),$

where $\Lambda$ is given in (59) and the step $(u)$ is computed via [25, Eq.
(3.326.2)].

##  Appendix G The derivation of $\mathcal{\hat{R}}_{ub}$ in (64)

The expectation term in (64) can be simplified by substituting (62) into (64)
as follows:

$\displaystyle\mathbb{E}\\!\left[{\hat{\gamma}^{*}}\right]=\bar{\gamma}\mathbb{E}\\!\left[{\left(\hat{R}+\rho
E\right)^{2}}\right]=\bar{\gamma}\left(\mathbb{E}\\!\left[{\hat{R}^{2}}\right]+\rho^{2}\mathbb{E}\\!\left[{E^{2}}\right]+2\mathbb{E}\\!\left[{\hat{R}E}\right]\right).$
(104)

The first expectation term $\mathbb{E}\\!\left[{\hat{R}^{2}}\right]$ in (104)
can be calculated by using (63a) as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\\!\left[{\hat{R}^{2}}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{(\hat{\alpha}_{u}+\hat{Y})^{2}}\right]=\sum_{n=0}^{2}\binom{2}{n}\mathbb{E}\\!\left[{\hat{\alpha}_{u}^{(2-n)}}\right]\mathbb{E}\\!\left[{\hat{Y}^{n}}\right]=\left(\sigma_{u}^{2}+\mu_{u}^{2}+\sigma_{Y}^{2}+\mu_{Y}^{2}+2\mu_{u}\mu_{Y}\right).$
(105)

Then, the second expectation term can be evaluated via substituting (63b) as

$\displaystyle\mathbb{E}\\!\left[{{E}^{2}}\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}\\!\left[{\left(\varepsilon_{u}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}\left(\rho\varepsilon_{g_{nl}}\varepsilon_{h_{nl}}+\hat{\alpha}_{h_{nl}}\varepsilon_{g_{nl}}+\hat{\alpha}_{g_{nl}}\varepsilon_{h_{nl}}\right)\right)^{2}}\right]$
(106) $\displaystyle\stackrel{{\scriptstyle(v)}}{{=}}$
$\displaystyle\mathbb{E}\\!\left[{\varepsilon_{u}^{2}}\right]+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}^{2}\left(\rho^{2}\mathbb{E}\\!\left[{\varepsilon_{g_{nl}}^{2}}\right]\mathbb{E}\\!\left[{\varepsilon_{h_{nl}}^{2}}\right]+\mathbb{E}\\!\left[{\hat{\alpha}_{h_{nl}}^{2}}\right]\mathbb{E}\\!\left[{\varepsilon_{g_{nl}}^{2}}\right]+\mathbb{E}\\!\left[{\hat{\alpha}_{g_{nl}}^{2}}\right]\mathbb{E}\\!\left[{\varepsilon_{h_{nl}}^{2}}\right]\right)$
$\displaystyle=$
$\displaystyle\beta_{u}^{2}+\sum_{n=1}^{N}\sum_{l=1}^{L}\eta_{nl}^{2}\left(\rho^{2}\beta_{g_{nl}}^{2}\beta_{h_{nl}}^{2}+{\xi}_{h_{nl}}\beta_{g_{nl}}^{2}+{\xi}_{g_{nl}}\beta_{h_{nl}}^{2}\right),$

where the step $(v)$ is written via the fact that $\hat{\alpha}_{v}$ and
$\varepsilon_{v}$ are independently distributed for
$v\in\\{u,h_{nl},g_{nl}\\}$. The final expectation term in (104) becomes zero
due to independence of the channel estimates and the estimation errors.

## References

  * [1] D. L. Galappaththige, D. Kudathanthirige, and G. Amarasuriya Aruma Baduge, “Performance Analysis of Distributed Intelligent Reflective Surface Aided Communications,” in _IEEE Glob. Commun. Conf. (GLOBECOM)_ , May 2020, pp. 1–6, (Accepted).
  * [2] C. Liaskos _et al._ , “A New Wireless Communication Paradigm Through Software-Controlled Metasurfaces,” _IEEE Commun. Mag._ , vol. 56, no. 9, pp. 162–169, 2018.
  * [3] M. D. Renzo _et al._ , “Smart Radio Environments Empowered by Reconfigurable AI Meta-Surfaces: An Idea Whose Time Has Come,” _EURASIP J. Wireless Commun. Net._ , May 2019.
  * [4] T. Sekitani _et al._ , “Stretchable Active-Matrix Organic Light-Emitting Diode Display Using Printable Elastic Conductors,” _Nature Materials_ , vol. 11, no. 11, p. 494?99, 2009.
  * [5] S. Zhang, Q. Wu, S. Xu, and G. Y. Li, “Fundamental Green Tradeoffs: Progresses, Challenges, and Impacts on 5G Networks,” _IEEE Commun. Surveys Tuts._ , vol. 19, no. 1, pp. 33–56, 2017.
  * [6] E. Basar _et al._ , “Wireless Communications Through Reconfigurable Intelligent Surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019\.
  * [7] W. Tang _et al._ , “Wireless Communications with Reconfigurable Intelligent Surface: Path Loss Modeling and Experimental Measurement,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2020.
  * [8] J. Chen, Y. Liang, Y. Pei, and H. Guo, “Intelligent Reflecting Surface: A Programmable Wireless Environment for Physical Layer Security,” _IEEE Access_ , vol. 7, pp. 82 599–82 612, 2019.
  * [9] D. Kudathanthirige, D. Gunasinghe, and G. Amarasuriya, “Performance Analysis of Intelligent Reflective Surfaces for Wireless Communication,” in _IEEE Int. Conf. on Commun. (ICC)_ , June 2020, pp. 1–6.
  * [10] E. Basar, “Transmission Through Large Intelligent Surfaces: A New Frontier in Wireless Communications,” _European Conf. on Netw. and Commun. (EuCNC)_ , pp. 112–117, 2019.
  * [11] Y. Han, W. Tang, S. Jin, C. Wen, and X. Ma, “Large Intelligent Surface-Assisted Wireless Communication Exploiting Statistical CSI,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 8, pp. 8238–8242, Aug 2019.
  * [12] M. Jung, W. Saad, Y. Jang, G. Kong, and S. Choi, “Reliability Analysis of Large Intelligent Surfaces (LISs): Rate Distribution and Outage Probability,” _IEEE Wireless Commun. Lett._ , vol. 8, no. 6, pp. 1662–1666, 2019.
  * [13] M. Jung, W. Saad, M. Debbah, and C. S. Hong, “On the Optimality of Reconfigurable Intelligent Surfaces (RISs): Passive Beamforming, Modulation, and Resource Allocation,” _arXiv:1910.00968_ , 2019.
  * [14] S. Hu and F. Rusek, “Spherical Large Intelligent Surfaces,” _IEEE Int. Conf. on Acoust., Speech and Signal Process._ , pp. 8673–8677, 2020.
  * [15] C. Psomas, I. Chrysovergis, and I. Krikidis, “Random Rotation-based Low-Complexity Schemes for Intelligent Reflecting Surfaces,” in _2020 IEEE 31st Ann. Int. Symp. on Personal, Indoor and Mobile Radio Commun._ , 2020, pp. 1–6.
  * [16] Q. Wu and R. Zhang, “Intelligent Reflecting Surface Enhanced Wireless Network via Joint Active and Passive Beamforming,” _IEEE Trans. Wireless Commun._ , pp. 1–1, 2019.
  * [17] ——, “Beamforming Optimization for Wireless Network Aided by Intelligent Reflecting Surface with Discrete Phase Shifts,” _IEEE Trans. Commun._ , vol. 68, no. 3, pp. 1838–1851, 2020.
  * [18] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable Intelligent Surfaces for Energy Efficiency in Wireless Communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [19] S. Abeywickrama, R. Zhang, and C. Yuen, “Intelligent Reflecting Surface: Practical Phase Shift Model and Beamforming Optimization,” in _IEEE Int. Conf. on Commun. (ICC)_ , 2020, pp. 1–6.
  * [20] S. Zhang and R. Zhang, “Intelligent Reflecting Surface Aided Multiple Access: Capacity Region and Deployment Strategy,” in _2020 IEEE 21st Int. Workshop on Signal Process. Adv. in Wireless Commun. (SPAWC)_ , 2020, pp. 1–5.
  * [21] C. Guo, Y. Cui, F. Yang, and L. Ding, “Outage Probability Analysis and Minimization in Intelligent Reflecting Surface-Assisted MISO Systems,” _IEEE Commun. Lett._ , pp. 1–1, 2020.
  * [22] M. Cui, G. Zhang, and R. Zhang, “Secure Wireless Communication via Intelligent Reflecting Surface,” _IEEE Wireless Commun. Lett._ , vol. 8, no. 5, pp. 1410–1414, 2019.
  * [23] X. Guan, Q. Wu, and R. Zhang, “Intelligent Reflecting Surface Assisted Secrecy Communication: Is Artificial Noise Helpful or Not?” _IEEE Wireless Commun. Lett._ , vol. 9, no. 6, pp. 778–782, 2020\.
  * [24] A. Papoulis and S. U. Pillai, _Probability, Random Variables, and Stochastic Processes_ , 4th ed. McGraw Hill, 2002.
  * [25] I. Gradshteyn and I. Ryzhik, _Table of Integrals, Series, and Products_ , 7th ed. Academic Press, 2007\.
  * [26] Q. Wu and R. Zhang, “Towards Smart and Reconfigurable Environment: Intelligent Reflecting Surface Aided Wireless Network,” _IEEE Commun. Mag._ , vol. 58, no. 1, pp. 106–112, 2020.
  * [27] T. M. Cover and J. A. Thomas, _Elements of Information Theory_. USA: Wiley-Interscience, 1991.
  * [28] S. Atapattu, C. Tellambura, and H. Jiang, “A Mixture Gamma Distribution to Model the SNR of Wireless Channels,” _IEEE Trans. Wireless Commun._ , vol. 10, no. 12, pp. 4193–4203, 2011.
  * [29] Q. Zhang, S. Jin, K.-K. Wong, H. Zhu, and M. Matthaiou, “Power Scaling of Uplink Massive MIMO Systems with Arbitrary-Rank Channel Means,” _IEEE J. Sel. Areas Signal Process._ , vol. 8, no. 5, pp. 966–981, Oct. 2014.
  * [30] J. Proakis, _Digital Communications_ , 4th ed. New York:McGraw-Hill, Inc., 2001.
  * [31] A. Mastin and P. Jaillet, “Log-Quadratic Bounds for the Gaussian Q-function,” _arXiv e-prints_ , 2013.
  * [32] Z. Wang and G. B. Giannakis, “A Simple and General Parameterization Quantifying Performance in Fading Channels,” _IEEE Trans. Commun._ , vol. 51, no. 8, pp. 1389–1398, Aug. 2003.
  * [33] G. Amarasuriya, C. Tellambura, and M. Ardakani, “Performance Analysis of Hop-by-Hop Beamforming for Dual-Hop MIMO AF Relay Networks,” _IEEE Trans. Commun._ , vol. 60, no. 7, pp. 1823–1837, 2012.
  * [34] S. M. Kay, _Fundamentals of Statistical Signal Processing: Estimation Theory_. Englewood Cliffs, NJ: Prentice Hall, 1993.
  * [35] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, _Fundamentals of Massive MIMO_. Cambridge University Press, Cambridge, UK, 2016.