# Performance and Optimization of Reconfigurable Intelligent Surface Aided THz
Communications

Hongyang Du, Jiayi Zhang, Ke Guan, Dusit Niyato, Huiying Jiao, Zhiqin Wang,
and Thomas Kürner H. Du is with the School of Computer Science and
Engineering, the Energy Research Institute @ NTU, Interdisciplinary Graduate
Program, Nanyang Technological University, Singapore (e-mail:
hongyang001@e.ntu.edu.sg).J. Zhang is with the School of Electronic and
Information Engineering, Beijing Jiaotong University, Beijing 100044, China.
(e-mail: jiayizhang@bjtu.edu.cn).K. Guan is with the State Key Laboratory of
Rail Traffic Control and Safety, Beijing Jiaotong University, Beijing 100044,
China. (e-mail: kguan@bjtu.edu.cn).D. Niyato is with the School of Computer
Science and Engineering, Nanyang Technological University, Singapore (e-mail:
dniyato@ntu.edu.sg)H. Jiao and Z. Wang are with China Academy of Information
and Communications Technology, Beijing 100191, P. R. China.T. Kürner is with
Institute for Communications Technology (IfN), Technische Universität
Braunschweig, 38106 Brunswick, Germany. (e-mail<EMAIL_ADDRESS>
###### Abstract

TeraHertz (THz) communications can satisfy the high data rate demand with
massive bandwidth. However, severe path attenuation and hardware imperfection
greatly alleviate its performance. Therefore, we utilize the reconfigurable
intelligent surface (RIS) technology and investigate the RIS-aided THz
communications. We first prove that the small-scale amplitude fading of THz
signals can be accurately modeled by the fluctuating two-ray distribution
based on two THz signal measurement experiments conducted in a variety of
different scenarios. To optimize the phase-shifts at the RIS elements, we
propose a novel swarm intelligence-based method that does not require full
channel estimation. We then derive exact statistical characterizations of end-
to-end signal-to-noise plus distortion ratio (SNDR) and signal-to-noise ratio
(SNR). Moreover, we present asymptotic analysis to obtain more insights when
the SNDR or the number of RIS’s elements is high. Finally, we derive
analytical expressions for the outage probability and ergodic capacity. The
tight upper bounds of ergodic capacity for both ideal and non-ideal radio
frequency chains are obtained. It is interesting to find that increasing the
number of RIS’s elements can significantly improve the THz communications
system performance. For example, the ergodic capacity can increase up to
$25\%$ when the number of elements increases from $40$ to $80$, which incurs
only insignificant costs to the system.

###### Index Terms:

Fluctuating two-ray fading, channel fitting, reconfigurable intelligent
surface, swarm intelligence optimization, THz communications.

##  I Introduction

In the past few years, wireless data traffic has been witnessing unprecedented
expansion. Several technological advances, such as artificial intelligence and
massive multiple-input multiple-output (MIMO) systems, have been presented [1]
while the wireless world moves towards the sixth generation (6G). Developing
new spectrum resources for future 6G communications systems is urgently
demanded.

Because of the high spectrum availability, TeraHertz (THz) communications
emerge as a key candidate in future generations [2, 3, 4, 5]. However, THz
signals have considerably high sensitivity to blockage, i.e., in urban
environments with high obstacle density. Thus, the direct usage of THz-band in
wireless access scenarios is challenging. Recently, the reconfigurable
intelligent surface (RIS) technology has drawn increasing attention [6, 7] due
to its ability to improve both spectrum and energy efficiency. RIS is based on
massive integration of low-cost tunable passive elements that can be easily
coated on existing infrastructures such as walls of buildings, facilitating
low-cost and low-complexity implementation [8, 9, 10, 11]. By smartly
adjusting the scattering properties of its reflecting elements, the RIS can
form equally strong beams in any direction for THz signals, which leads to a
more controllable wireless environment. Compared to conventional relays, RISs
can enhance the transmission performance without being impaired by self-
interference and thermal noise. More importantly, RIS can be deployed to
leverage the line-of-sight (LoS) components with respect to both the base
station (BS) and the users to assist their end-to-end THz communications,
which potentially reduce the path attenuation that originates from the high
frequencies.

However, how to model the small-scale amplitude fading of THz signals in the
RIS-aided system is still an open question. Methods of modeling and evaluating
the performance of THz communications have been extensively studied in the
literature [12, 13, 14, 15, 16, 17, 18, 19, 20]. For example, the authors in
[19] delivered a simplified molecular absorption loss model that was used in
[20] for the quantification of the THz link capacity. However, the THz system
performance is overestimated by neglecting small-scale fading. Although the
small-scale fading is taken into consideration in [17, 18], the channel
modelings were not based on real THz signal measurements. Furthermore, in
[21], the envelope of the fading coefficient was proved by the authors
experimentally to follow the Nakagami-$m$ distribution under both non-line-of-
sight (NLoS) and LoS conditions. Similarly, experimental verification of the
existence of shadowing in the $300$ GHz band were exploited in [22]. Very
recently, a THz measurement campaign was conducted in a train test center with
trains, tracks, and lampposts [23]. We fit the measurement data to Gaussian,
Rician, Nakagami-$m$, and fluctuating two-ray (FTR) distributions [24],
respectively, and prove that the FTR distribution, which is a generalization
of the two-wave with diffuse power fading model and allows the constant
amplitude specular waves of LoS propagation to fluctuate randomly, performs a
much better fit than others.

Therefore, with the help of the FTR fading model, we investigate the RIS-aided
THz communications. We propose an algorithm to obtain the optimal design of
RIS’s reflector array. By designing the phase shift, we can then theoretically
derive the maximum of SINR and corresponding probability density function
(PDF) and cumulative distribution function (CDF) expressions. The main
contributions of this paper are summarized as follows:

  * •

To model the small-scale amplitude fading of THz signals, we prove that the
FTR distribution performs a much better fit than those of the other fading
models with the help of the THz measurement data [23, 25]. Moreover, we derive
exact PDF and CDF expressions of end-to-end SNR and SNDR in terms of the
multivariate Fox’s $H$-function, respectively. The impact of transceiver
antenna gains, the level of transceiver hardware imperfections, the
transmission distance, the misalignment between the RIS and the destination,
the operating frequency, and the environmental conditions including
temperature, humidity, and pressure, are all comprehensively considered in our
analysis.

  * •

For the practical RIS with a finite number of phase-shifters at each element,
we propose a swarm intelligence-based optimization (SIO) algorithm to obtain
the optimal phase-shifts. Furthermore, the impact of the number of discrete
phase-shift levels on the transmission performance is revealed, and we find
that RIS with discrete phase-shifts can achieve performance close to that of
the continuous phase-shifts.

  * •

We derive exact closed-form expressions of outage probability for the RIS-
aided THz communications system with both ideal and non-ideal radio frequency
(RF) chains. Moreover, asymptotic analysis is presented when the SNDR or the
number of RIS’s elements is high, which is easier to calculate and obtain
insights. For the ergodic capacity, the exact closed-form expression for ideal
RF chains is obtained and the tight upper bound expressions for both ideal and
non-ideal RF chains are derived. Important physical and engineering insights
have been revealed. The derived results show that the propagation loss is more
detrimental than the molecular absorption loss, and the performance is
improved significantly by increasing the number of reflecting elements
equipped in an RIS.

The remainder of the paper is structured as follows. We briefly introduce the
RIS-aided THz communications system and signal models in Section II. Section
III presents the optimal phase-shifts of the RIS’s reflector array under
discrete level constraints. In Section IV, we derive exact PDF and CDF
expressions of end-to-end SNR and SNDR of the considered system. Section V
presents the asymptotic analysis to obtain several insights. Important
performance metrics, such as outage probability and ergodic capacity, are
obtained in Section VI. Numerical results accompanied with Monte-Carlo
simulations are presented in Section VII. Finally, the conclusions are drawn
in Section VIII.

##  II System and Signal Models

###  II-A System Description

Figure 1: System model of the RIS-aided THz communications.

We focus on a single cell of the THz communications system. As shown in Fig.
1, the Source is equipped with $M$ antennas, and the Destination has one
antenna. The interference among THz links can be significantly small because
the signal in THz communications is highly directional. Thus, the interference
from other cell-users is ignored. Because THz signals suffer from high path
and penetration losses, an RIS is set up to assist the communications. The RIS
is a uniform planar array (UPA) containing $L$ nearly passive reconfigurable
reflecting units, where each element is capable of electronically controlling
the phase of the incident electromagnetic waves. To establish a practical RIS-
aided THz communications model, we consider the following coefficients and
give the corresponding reasons:

  * •

Small-Scale Fading Coefficient ${\mathbf{h}}_{F}$: Although some of the
existing literature on RIS-aided THz communications considers the utilization
of small-scale fading models, classical models such as Rayleigh and Rician
fadings are typically used. However, recent real measurements of THz signals
show that these classical models cannot fit well the THz signals. Although the
authors of [25] claim that $\alpha\\!-\\!\mu$ fading obtains good results and
passes the Kolmogorov-Smirnov (K-S) test, we can observe that the
$\alpha\\!-\\!\mu$ fading is still not accurate enough in fitting the bimodal
characteristics of THz signal’s amplitude fading. Therefore, a more accurate
channel model is needed.

  * •

Misalignment Fading Coefficient $h_{P}$: Random building sways, which will be
caused by wind loads and thermal expansions, will result in pointing errors at
Destination. Because the THz signals are highly directional, we need to take
the misalignment into consideration.

  * •

Path Gain Coefficient $h_{L}$: Because of the THz signals’ inability to
penetrate solid objects, the path gain in the THz communications system
concludes the propagation gain and the molecular absorption gain.

  * •

Hardware Impairments $n_{S}$ and $n_{D}$: Physical transceivers have hardware
impairments that create distortions that degrade the performance of
communications systems. The vast majority of technical contributions in the
area of relaying neglect hardware impairments and, thus, assume ideal
hardware. Such approximations make sense in low-rate systems, but can lead to
very misleading results when analyzing future high-rate systems. This paper
quantifies the impact of hardware impairments on the RIS-aided THz system.

By considering the above coefficients and the characteristics of RIS
carefully, we can derive the expression of the received signal, which allows
us to further analyze the statistical characteristics of SNR and performance
metrics. A more detailed analysis of the above coefficients are given as
follows:

####  II-A1 Small-Scale Fading Measurement and Modeling

We consider the RIS is electrically small and the correlation between the
amplitudes of the signals is weak enough to be ignored [26, 27]. Let
${{\mathbf{g}}_{1}}\in{\mathbb{C}^{M\times L}}$ and
${{\mathbf{g}}_{2}}\in{\mathbb{C}^{1\times L}}$ denote the channels from the
THz Source to the RIS, and from the RIS to the Destination, respectively,
${\bf{\Phi}}=\beta{\rm{diag}}\left({\bf{\varphi}}\right)$ denote the diagonal
reflection coefficients matrix of the RIS, with $\beta\in(0,1]$ is the phase-
shifts matrix111The value of $\beta$ does not affect the algorithm design.
Without loss of generality, we set $\beta=1$ for simplicity.. The small-scale
fading between BS and RIS is defined as
${\mathbf{h}_{F}}\triangleq{{\mathbf{g}}_{2}}{\bf{{\Phi}}}{{\mathbf{g}}_{1}^{\rm
H}}$, where
${\bf{\varphi}}={\left({{e^{i{\varphi_{1}}}},\ldots,{e^{i{\varphi_{L}}}}}\right)^{T}}$,
$\varphi_{\ell}$ is the phase-shift of the $\ell_{\rm th}$ element of the RIS,
$i={\sqrt{-1}}$,
${{\bf{g}}_{2}}=\left({g_{2}^{\left(1\right)},\ldots,g_{2}^{\left(L\right)}}\right)=\left({\left|{g_{2}^{\left(1\right)}}\right|{e^{-i\theta_{2}^{\left(1\right)}}},\ldots,\left|{g_{2}^{\left(L\right)}}\right|{e^{-i\theta_{2}^{\left(L\right)}}}}\right)$,
and

$\displaystyle{{\bf{g}}_{1}}$
$\displaystyle=\left({\begin{array}[]{*{20}{c}}{g_{1}^{\left({1,1}\right)}}&\cdots&{g_{1}^{\left({1,L}\right)}}\\\
\vdots&\ddots&\vdots\\\
{g_{1}^{\left({M,1}\right)}}&\cdots&{g_{1}^{\left({M,L}\right)}}\end{array}}\right)$
(4)
$\displaystyle=\left({\begin{array}[]{*{20}{c}}{\left|{g_{1}^{\left({1,1}\right)}}\right|{e^{-i\theta_{1}^{\left({1,1}\right)}}}}&\cdots&{\left|{g_{1}^{\left({1,L}\right)}}\right|{e^{-i\theta_{1}^{\left({1,L}\right)}}}}\\\
\vdots&\ddots&\vdots\\\
{\left|{g_{1}^{\left({M,1}\right)}}\right|{e^{-i\theta_{1}^{\left({M,1}\right)}}}}&\cdots&{\left|{g_{1}^{\left({M,L}\right)}}\right|{e^{-i\theta_{1}^{\left({M,L}\right)}}}}\end{array}}\right).$
(8)

The elements in ${{\mathbf{g}}_{1}}$ and ${{\mathbf{g}}_{2}}$ represent the
small-scale fading of the THz signals, for which we have to find the
appropriate channel fading model. Recently, the FTR fading [24] has received a
great amount of research interest, because the FTR model yields a better fit
than those of the other conventional channel models with experimental data in
high-frequency communications. The PDF of the squared FTR RV $\gamma$ can be
expressed as [28]

${f_{\gamma}}\left(\gamma\right)=\frac{{{m^{m}}}}{{\Gamma(m)}}\sum\limits_{j=0}^{N_{T}}{\frac{{{K^{j}}{d_{j}}}}{{j!}}}\frac{{{\gamma^{j}}}}{{\Gamma(j+1){{\left({2{\sigma^{2}}}\right)}^{j+1}}}}\exp\left({-\frac{\gamma}{{2{\sigma^{2}}}}}\right),$
(9)

where

$\displaystyle{d_{n}}\buildrel\Delta\over{=}\sum\limits_{k=0}^{n}{\left(\\!{\begin{array}[]{*{20}{l}}n\\\
k\end{array}}\\!\right)}{\left(\\!{\frac{\Delta}{2}}\\!\right)^{k}}\sum\limits_{l=0}^{k}{\left({\begin{array}[]{*{20}{c}}k\\\
l\end{array}}\right)}\Gamma\\!(n+m+2l-k){e^{\frac{{\pi(2l-k)j}}{2}}}$ (14)
$\displaystyle\times\\!\\!{\left(\\!{{{(m\\!+\\!K)}^{2}}\\!-\\!{{(K\Delta)}^{2}}}\\!\right)^{\frac{{-(n+m)}}{2}}}\\!P_{n+m-1}^{k-2l}\left(\\!\\!{\frac{{m+K}}{{\sqrt{{{(m\\!+\\!K)}^{2}}\\!-\\!{{(K\Delta)}^{2}}}}}}\\!\\!\right)\\!,$
(15)

where $P_{\cdot}^{\cdot}\left(\cdot\right)$ denotes the Legendre function of
the first kind [29, eq. (8.702)], $m$ is the fading severity parameter, $K$ is
the average power ratio of the dominant wave to the scattering multipath, and
$\Delta$ is a parameter varying from $0$ to $1$ representing the similarity of
two dominant waves. The received average SNR is given by
$\upsilon=2\sigma^{2}(1+K)$. Note that $N_{T}$ is the truncation term. Using
the property of PDF, we have
$\int_{0}^{\infty}{{f_{\gamma}}\left(\gamma\right)}=1$. With the help of [29,
eq. (3.351.3)] and [29, eq. (8.339.1)], the $N_{T}$ should be set to satisfy
$\sum\limits_{j=0}^{{N_{T}}}{{d_{j}}}\to 1$. Typically, we can use $40$ terms
to make the series converge for all cases considered in this paper.
Furthermore, the truncation error, which is defined as
$1-\sum\limits_{j=0}^{{N_{T}}}{{d_{j}}}$, is less than $10^{-6}$.

To prove that the FTR model can fit the THz signal’s amplitude fading best, we
consider the real measurement data from two campaigns [23, 25], and evaluate
the fitting of those theoretical fading distributions to the empirical data
with the help of K-S goodness-of-fit statistical test [30].

  * •

Very recently, a measurement campaign was conducted in a train test center
[23]. The authors use a novel $M$-sequence correlation-based ultra-wideband
channel sounder to measure the THz channel. In the measurement, the channel
sounder was operated at $304.2$ GHz with $8$ GHz bandwidth. The small-scale
fading was extracted for each channel impulse response based on the data. We
have fitted the measurement data to a Gaussian, Rician, Nakagami-$m$, and FTR
distribution, respectively, which are also compared with the CDF of the
measurement data. As shown in Figs. 2, we find that the FTR distribution
performs the best fit.

  * •

In [25], three sets of experimental measurements, which have been conducted in
a shopping mall, an airport check-in area, and an entrance hall of a
university towards different time periods, are used to accurately model the
THz signal fading distribution. Analysis shows that conventional
distributions, e.g., Rayleigh, Rice, and Nakagami-$m$, lack fitting accuracy,
whereas the $\alpha\\!-\\!\mu$ distribution has an almost-excellent fit.
However, we can observe that the $\alpha\\!-\\!\mu$ distribution is still
deficient to show the characterizations of the experimental results: 1) The
experimental PDF appears bimodality in many test environments, e.g.,
Transmitter 4 (${\rm{TX}}_{4}$) in the entrance hall or the shopping mall (see
Fig. 2 in [25]), which cannot be fitted by the PDF of $\alpha\\!-\\!\mu$
distribution. 2) As the variable $x$ increases, the slope of CDF gradually
increases and then decreases to $0$, but the CDF of $\alpha\\!-\\!\mu$
distribution cannot fit the process of the increasing slope of CDF when x is
small. As shown in Fig. 3, we can observe that the FTR fading is better than
the $\alpha\\!-\\!\mu$ fading. Furthermore, although both the
$\alpha\\!-\\!\mu$ and FTR distribution can pass the K-S test, the values for
the FTR distribution is always smaller.

Figure 2: (a). PDF of the measurement data [23] and fitted distributions. The parameters of fitted distributions are obtained by ${\rm fitdist(\cdot)}$ function in Matlab [31] and are shown in Table I (b). CDF of the measurement data and fitted distributions with the same parameters as Fig. 2 (a). Figure 3: (a). PDF of the measurement data [25] and fitted distributions. The parameters of $\alpha\\!-\\!\mu$ distribution are given in [25, Table 5, ${\rm{TX}}_{4}$], i.e., $\alpha=8.54867$, $\mu=0.23377$, and $\beta=8.49874$. The parameters of FTR distribution are $\sigma=0.8991$, $m=30$, $K=30$, and $\Delta=0.8$ (b). CDF of the measurement data [25] and fitted distributions. The parameters of $\alpha\\!-\\!\mu$ distribution are given in [25, Table 1, ${\rm{TX}}_{1}$ and ${\rm{TX}}_{15}$]. The parameters of FTR distribution are shown in Table I. TABLE I: Kolmogorov–Smirnov goodness of fit tests. Measurement |  Channel Model |  Parameters |  K-S Test Values
---|---|---|---
[23] |  FTR Fading |  $\sigma=4.84\times 10^{-3}$, $m=0.63$, $K=2.9$, $\Delta=0.3$ |  $0.0383$ ($\sqrt{}$)
[23] |  Gaussian Fading |  ${\mu_{G}}=1.83\times 10^{-4}$, ${\sigma_{G}}=2.01\times 10^{-4}$ |  $0.1870$ ($\times$)
[23] |  Rician Fading |  ${s_{R}}=7.78\times 10^{-6}$, ${\sigma_{R}}=1.92\times 10^{-4}$ |  $0.3687$ ($\times$)
[23] |  Nakagami-$m$ Fading |  ${\mu_{N}}=0.323$, ${\omega_{N}}=7.35\times 10^{-8}$ |  $0.1058$ ($\times$)
[25] ${\rm{TX}}_{1}$ |  FTR Fading |  $\sigma=2.7729$, $m=6$, $K=7$, $\Delta=0.2$ |  $0.0175$ ($\sqrt{}$)
[25] ${\rm{TX}}_{1}$ |  $\alpha\\!-\\!\mu$ Fading |  $\alpha=3.35467$, $\mu=1.11365$, and $\beta=11.58104$ |  $0.0183$ ($\sqrt{}$)
[25] ${\rm{TX}}_{15}$ |  FTR Fading |  $\sigma=2.382$, $m=3$, $K=6$, $\Delta=0.04$ |  $0.0203$ ($\sqrt{}$)
[25] ${\rm{TX}}_{15}$ |  $\alpha\\!-\\!\mu$ Fading |  $\alpha=3.11872$, $\mu=0.91962$, and $\beta=9.8334$ |  $0.0566$ ($\sqrt{}$)

Moreover, Table I shows that The FTR distribution achieves the smallest K-S
test values in the fit of the experimental results of each THz signal
measurement. Thus, in the following, we use the FTR fading model to illustrate
the small-scale amplitude fading of THz communications.

####  II-A2 Misalignment Fading Coefficient

We consider that Destination has a circular detection beam, and the radius of
the beam is $\alpha$. Furthermore, the RIS has a circular beam of radius
$\rho$, where ${\rm{0}}\leq\rho\leq{\omega_{d}}$, and $w_{d}$ denotes the
maximum radius of the beam. Moreover, we assume that $r$ is the pointing
error. Due to the symmetry of the beam shapes, the misalignment fading
coefficient that represents the fraction of the power collected by
Destination, $h_{P}$, depends only on the radial distance $r$. Thus, $h_{P}$
can be approximated as [32, eq. (11)]

$\displaystyle{f_{{h_{P}}}}(x)={\gamma^{2}}A_{o}^{-{\gamma^{2}}}{x^{{\gamma^{2}}-1}},\>\>\>\>\>0\leqslant
x\leqslant{A_{o}},$ (16)

where ${A_{o}}={\left[{{\mathop{\rm
erf}\nolimits}\left(u\right)}\right]^{2}}$,
$u\triangleq{\sqrt{\pi}a}/({\sqrt{2}w_{d}})$,
$\gamma=\frac{w_{eq}}{2\sigma_{S}}$, $\sigma_{S}^{2}$ denotes the variance of
the pointing error displacement at Destination, ${\omega_{eq}}$ is the
equivalent beam-width, and
$w_{eq}^{2}=w_{d}^{2}\sqrt{\pi}\operatorname{erf}(u){\left({2u\exp\left({-{u^{2}}}\right)}\right)^{-1}}$.

####  II-A3 Path Gain Coefficient

The path gain coefficient can be evaluated as ${h_{L}}={h_{FL}}{h_{AL}}$,
where ${h_{FL}}$ is the propagation gain, and ${h_{AL}}$ is the molecular
absorption gain. Using Friis equation, we can express $h_{FL}$ as

$h_{FL}=\frac{c^{2}\sqrt{G_{t}G_{r}}}{\left(4\pi
f\right)^{2}\left(d_{1}d_{2}\right)},$ (17)

where $G_{t}$ and $G_{r}$ represent the antenna orientation dependent
transmission and reception gains, respectively, $c$ is the speed of light, $f$
stands for the operating frequency, $d_{1}$ is the distance between Source and
RIS, and $d_{2}$ is the distance between RIS and Destination. Furthermore,
${h_{AL}}$ can be evaluated as [19]

$h_{AL}=\exp\left(-\frac{1}{2}\kappa_{\alpha}(f)\left(d_{1}+d_{2}\right)\right),$
(18)

where $\kappa_{\alpha}(f)$ denotes the absorption coefficient which can be
expressed as [19]

$\displaystyle{\kappa_{\alpha}}(f)=$
$\displaystyle\frac{{A(v)}}{{B(v)+{{\left({\frac{f}{{100c}}-{c_{1}}}\right)}^{2}}}}+\frac{{C(v)}}{{D(v)+{{\left({\frac{f}{{100c}}-{c_{2}}}\right)}^{2}}}}$
$\displaystyle+{p_{1}}{f^{3}}+{p_{2}}{f^{2}}+{p_{3}}f+{p_{4}},$ (19)

where $c$ denotes the speed of light, $c_{1}=10.835$ $\mathrm{cm}^{-1}$,
$c_{2}=12.664$ $\mathrm{cm}^{-1}$, $p_{1}=5.54\times 10^{-37}$
$\mathrm{Hz}^{-3}$, $p_{2}=-3.94\times 10^{-25}$ $\mathrm{Hz}^{-2}$,
$p_{3}=9.06\times 10^{-14}$ $\mathrm{Hz}^{-1}$, $p_{4}=-6.36\times 10^{-3}$,
$A(v)=0.2205v(0.1303v+0.0294)$, $B(v)={(0.4093v+0.0925)^{2}}$,
$C(v)=2.014v(0.1702v+0.0303)$, and $D(v)={(0.537v+0.0956)^{2}}$ [19].
Moreover, $v$ can be evaluated as

$v\\!=\\!\phi\\!\left({0.06116{p^{-1}}\\!+\\!2.1148\\!\times\\!{{10}^{-7}}}\right)\exp\\!\left(\\!{\frac{{17.502T}}{{240.97\\!+\\!T}}}\\!\right),$
(20)

where $T$ is the temperature, $\phi$ and $p$ denote the relative humidity and
the atmospheric pressure, respectively. Note that the molecular absorption
loss model that we considered here provides an accurate model for the $275$ –
$400$ GHz frequency band. Meanwhile, for frequencies higher than $400$ GHz, we
can obtain the value of $\kappa_{\alpha}(f)$ from the high resolution
transmission (HITRAN) database [33, 34]. Using (17)-(20), we finally obtain
the path gain coefficient of the considered RIS-aided THz communications
system.

###  II-B Received THz Signal Expression Analysis

The signal received from Source through RIS for Destination can be expressed
as222Because the THz signals have high susceptibility to blockages [33], we
assume that the direct link between the source and destination is negligible
due to the severe attenuation.

$y=P{h_{L}}{h_{P}}{{\mathbf{h}}_{F}}{\bf{w}^{T}}\left({x+{n_{S}}}\right)+{n_{D}}+n,$
(21)

where $n\sim{\mathcal{CN}}\left(0,N_{0}\right)$,
${\mathcal{CN}}\left({t_{1},t_{2}}\right)$ denotes the white Gaussian process
with mean $t_{1}$ and variance $t_{2}$, $x$ is the transmitted signal
satisfying ${\mathbb{E}}\left[{{{\left|x\right|}^{2}}}\right]=1$,
${\bf{w}}\in{\mathbb{C}^{M\times 1}}$ is the beamforming vector,
${\mathbb{E}}\left[{{\left\|{{{\bf{w}}}}\right\|^{2}}}\right]=1$,
$\left\|\cdot\right\|$ denotes the Frobenius norm, $P$ is the transmit power,
the parameters ${n_{S}}\sim{\mathcal{CN}}\left({0,\kappa_{S}^{2}P}\right)$ and
${n_{D}}\sim{\mathcal{CN}}\left(0,\kappa_{D}^{2}P{{\left|{h_{L}}{h_{P}}{h_{F}}\right|}^{2}}\right)$
are distortion noise from hardware impairments in both Source and Destination,
respectively [35], and ${\kappa_{S}}$ and ${\kappa_{D}}$ are non-negative
parameters that characterize the level of hardware impairments [36] in Source
and Destination, respectively. In the RIS-aided communications, the $P$ and
${\bf{\Phi}}$ are adjustable. Considering the optimal beamforming vector
design, e.g., using the maximum ratio transmitting [37, 38], we can define
${\bf{w}}$ as
${\bf{w}}={{{\mathbf{h}}_{F}}}/\left\|{{{\mathbf{h}}_{F}}}\right\|$, and the
corresponding small-scale fading coefficient can be expressed as

$\displaystyle{{{\mathbf{h}}_{F}}}$
$\displaystyle=\\!\left(\\!{{h_{F,1}},\\!\ldots\\!,{h_{F,M}}}\\!\right)\\!$
$\displaystyle=\left(\\!{\left|{\sum\limits_{\ell=1}^{L}{g_{2}^{\left(\ell\right)}{e^{j{\varphi_{\ell}}}}g_{1}^{\left({1,\ell}\right)}}}\right|,\\!\ldots\\!,\left|{\sum\limits_{\ell=1}^{L}{g_{2}^{\left(\ell\right)}{e^{j{\varphi_{\ell}}}}g_{1}^{\left({M,\ell}\right)}}}\right|}\\!\right)\\!\\!.$
(22)

The square of the Frobenius norm of ${{{\mathbf{h}}_{F}}}$ is upper-bounded as

${\left\|{{{\mathbf{h}}_{F}}}\right\|^{2}}\\!=\\!\\!{\left(\\!{\sum\limits_{m=1}^{M}{\left|{\sum\limits_{\ell=1}^{L}{g_{2}^{\left(\ell\right)}{e^{j{\varphi_{\ell}}}}g_{1}^{\left({1,\ell}\right)}}}\right|}}\\!\right)^{2}}\\!\\!\leqslant\\!\\!{\left(\\!{\sum\limits_{m=1}^{M}{\sum\limits_{\ell=1}^{L}{\left|{g_{2}^{\left(\ell\right)}}\right|\left|{g_{1}^{\left({m,\ell}\right)}}\right|}}}\\!\right)^{2}}\\!,$
(23)

where the error between the upper-bound and the exact value of
${\left\|{{{\mathbf{h}}_{F}}}\right\|^{2}}$ decreases as RIS’s phase shift
design becomes more accurate. However, the right-hand side of (23) is still
hard to be used in the performance analysis. The reason is that
$\left|{g_{2}^{\left(\ell\right)}}\right|$ and
$\left|{g_{1}^{\left({m,\ell}\right)}}\right|$ follow i.n.i.d. FTR fading
distribution. Deriving the PDF of ${\left\|{{{\mathbf{h}}_{F}}}\right\|^{2}}$
is equivalent to deriving the PDF of square of two-fold sum of product of FTR
RVs, which is extremely difficult, if not impossible, to be obtained.
Fortunately, the distribution of the sum of FTR random variables (RVs) can be
accurately approximated by the single FTR distribution by following the method
proposed in [39], i.e.,
$\sum\limits_{m=1}^{M}{\left|{g_{1}^{\left({m,\ell}\right)}}\right|}\approx\left|{g_{1}^{\left(\ell\right)}}\right|$.
Thus, with the optimal phase shift design, we have

$\displaystyle{\left\|{{{\mathbf{h}}_{F}}}\right\|^{2}}$
$\displaystyle\\!\approx\\!{\left({\sum\limits_{m=1}^{M}{\sum\limits_{\ell=1}^{L}{\left|{g_{2}^{\left(\ell\right)}}\right|\left|{g_{1}^{\left({m,\ell}\right)}}\right|}}}\right)^{2}}\\!=\\!{\left(\\!{\sum\limits_{\ell=1}^{L}{\left|{g_{2}^{\left(\ell\right)}}\right|}\sum\limits_{m=1}^{M}{\left|{g_{1}^{\left({m,\ell}\right)}}\right|}}\\!\right)^{2}}$
$\displaystyle\approx\\!{\left({\sum\limits_{\ell=1}^{L}{\left|{g_{2}^{\left(\ell\right)}}\right|}\left|{g_{1}^{\left(\ell\right)}}\right|}\right)^{2}},$
(24)

where the two-fold summation is reduced to a one-fold summation, thus giving
the possibility of further performance analysis. Note that (II-B) is also the
exact small-scale fading expression when considering a single antenna
transmitter. In the following, we analyze the SNR expressions in two cases.

####  II-B1 Non-Ideal RF Chains

With the non-ideal RF chains, the instantaneous SNDR can be expressed as

${\gamma_{N}}=\frac{{{{\left\|{{h_{F}}}\right\|}^{2}}{{\left|{{h_{P}}}\right|}^{2}}{{\left|{{h_{L}}}\right|}^{2}}P}}{{{\kappa^{2}}{{\left\|{{h_{F}}}\right\|}^{2}}{{\left|{{h_{P}}}\right|}^{2}}{{\left|{{h_{L}}}\right|}^{2}}P+{N_{o}}}},$
(25)

where $\kappa^{2}\triangleq\kappa_{S}^{2}+\kappa_{D}^{2}$.

####  II-B2 Ideal RF Chains

Let us consider the case of ideal RF chains. Using (25) and letting
${\kappa=0}$, we can express the instantaneous SNR of the received signal as

${\gamma_{I}}=\frac{{{{{\left\|{{h_{F}}}\right\|}^{2}}{{\left|{{h_{P}}}\right|}^{2}}}{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}.$
(26)

##  III Phase Shift Design Under Discrete Phase Shift Constraints

In this section, we propose an SIO algorithm to obtain the optimal design of
RIS’s reflector array, ${\bf{\Phi}}$. With the optimal phase shift design of
RIS, we can theoretically obtain the maximum of SNDR. In practical systems,
ideal continuous phase-shifts are hard to implement because of hardware
limitation. Thus, we consider the discrete phase-shifts design, which is more
energy efficient [40].

###  III-A Optimal Phase Shift Design of RIS’s Reflecting Elements

####  III-A1 Motivations for proposing the SIO algorithm

There are several challenges for the existing RIS phase shift design
algorithms, which can be summarized as follows:

  * •

When the number of reflecting elements at the RIS, $L$, is large, most of the
existing algorithms [41, 7, 40, 42] require a lot time to converge, which
motivates us to propose an effective algorithm for practical RIS-aided
systems. The running time should not be greatly effected by $L$.

  * •

Eliminating the channel phases is difficult because the RIS is assumed to be
passive. Although a novel and simple algorithm based on the binary search tree
was proposed to obtain ${\bf{\Phi}}_{\rm opt}$ without the channel phases
information [41], the algorithm proposed in [41] is based on an assumption
that the phase-shifts is continuous at each reflecting element, which is
practically difficult to realize due to the hardware limitation of the RIS.

Combining the above discussion, we design an SIO algorithm to find
${\bf{\Phi}}_{\rm opt}$ for the RIS with finite-level phase-shifters without
the channel phases information, inspired from Particle Swarm Optimization
(PSO) and evolution algorithms.

####  III-A2 Introduction of the traditional PSO algorithm

For the traditional PSO algorithm, $n$ particles move in an $L$-dimensional
search space to optimize a fitness function [43]. In a real RIS-aided THz
communications system, the value of fitness function can be obtained by
measuring the expectation of the amplitude of the received signal, $E_{\rm
re}$, by averaging the received signal strength in the time domain. First, the
particles are placed randomly in the search space. Then, for a predefined
number of iterations, each particle evaluates fitness and moves to a new
location based on the history of particles’ own best previous location and
other particles’ best positions. Thus, the PSO algorithm consists of a swarm
of $n$ particles, and the position of each particle represents a possible
solution of the fitness function. The velocity $v_{i,j}$ and position
$x_{i,j}$ of the $d_{\rm th}$ dimension of the $i_{\rm th}$ particle are
updated as follows:

$\displaystyle{v_{id}}\left({k+1}\right)=$
$\displaystyle\omega{v_{id}}\left(k\right)+{c_{1}}{R_{1}}\left({{p_{id}}\left(k\right)-{x_{id}}\left(k\right)}\right)$
$\displaystyle+{c_{2}}{R_{2}}\left({{p_{gd}}\left(k\right)-{x_{id}}\left(k\right)}\right),$
(27)
${x_{id}}\left({k+1}\right)={x_{id}}\left(k\right)+{v_{id}}\left({k+1}\right),$
(28)

where $i=1,2,\ldots,m$, $d=1,2,\ldots,L$, $\omega$ is the inertia weight,
$c_{1}$ and $c_{2}$ are the acceleration constants, $R_{1}$ and $R_{2}$ are
uniformly distributed random numbers in $[0,1]$,
${P_{i}}=\left({{p_{i1}},{p_{i2}}\ldots,{p_{iL}}}\right)$ is the best previous
position of the $i_{\rm th}$ particle, and
${P_{g}}=\left({{p_{g1}},{p_{g2}}\ldots,{p_{gL}}}\right)$ is the best position
among all the particles.

####  III-A3 Design of the SIO algorithm

The main disadvantage of the PSO algorithm is that it is difficult to balance
local and global search and keep the diversity of the population. Thus, by
combining the modified PSO algorithm [44, 45] with evolution algorithms, we
present the SIO Algorithm 1, which maintains the optimal balance between
exploration and exploitation, and reduces computation time, to find
${\bf{\Phi}}_{\rm opt}$. In the following, we introduce the SIO’s working
mechanism and the important parameters:

  * •

Each particle is given a periodic value, ${\rm Pe}$, whose initial value is
$0$. When the particle updates to a better state, the periodic value resets to
$0$, and otherwise, the value increases by $1$. We set the size of swarm
reaches as $n_{\rm begin}$. As the searching time increases, the particle that
has the largest periodic value will be deleted, until the size of the swarm
reaches $n_{\rm end}$. We can set $n_{\rm begin}=30$ and $n_{\rm end}=10$ to
reduce the computational resources and keep the diversity of the population.

  * •

The inertia weight $\omega$ can significantly affect the global and local
search-ability. By linearly decreasing $\omega$ from a relatively larger value
to a smaller value, the algorithm’s global search-ability at the beginning and
the local search-ability near the end of the search procedure can be
efficiently enhanced [46]. In experimental work (i.e., in Figs. 4-6), $\omega$
is kept in between $0.9$ to $0.4$.

  * •

Adding an exploitative component allows the algorithm to utilize local
information to guide the search towards better regions in the search space.
The accelerate constants, i.e., $c_{3}$, $c_{1}$ and $c_{2}$, represent the
particle stochastic acceleration weights toward the best positions of the
local, $i_{\rm th}$ particle and among all the particles. The system designer
can set larger accelerate constants for fast search, or smaller accelerate
constants to reduce oscillations of particles.

In contrast to the continuous phase-shifts, we consider the practical case
where the RIS only has a finite number of discrete phase-shifts that are
equally spaced in $[0,2\pi)$. The set of phase-shifts at each element is given
by $\mathcal{O}=\\{0,\Delta\theta,\ldots,\Delta\theta(K-1)\\}$, where
$K=\left\lceil{\frac{{{\rm{2}}\pi}}{{\Delta\theta}}}\right\rceil$. If we use
$b$ bits to represent each of the levels, we can also express $K$ as
$K=2^{b}$. Denote the best position of the local as
${P_{l}}=\left({{p_{l1}},{p_{l2}},\ldots,{p_{lD}}}\right)$, after finding the
three best values, i.e., ${P_{l}}$, ${P_{i}}$, and ${P_{g}}$, the velocities
and positions are updated with the following formulas:

${v_{id}}\left({k+1}\right)=\left[\left(\omega{v_{id}}\left(k\right)+\Delta
v\left(k\right)\right)/\Delta\theta\right]\times\Delta\theta,$ (29)
$\displaystyle\Delta v\left(k\right)=$
$\displaystyle{c_{1}}{R_{1}}\left({{p_{id}}\left(k\right)-{x_{id}}\left(k\right)}\right)+{c_{2}}{R_{2}}\left({{p_{gd}}\left(k\right)-{x_{id}}\left(k\right)}\right)$
$\displaystyle+{c_{3}}{R_{3}}\left({{p_{ld}}\left(k\right)-{x_{id}}\left(k\right)}\right),$
(30)
${x_{id}}\left({k+1}\right)={x_{id}}\left(k\right)+{v_{id}}\left({k+1}\right),$
(31)

where $[\cdot]$ is the rounding function. From (29)-(31), we can observe that,
in practical communications systems, we only need to measure the $E_{\rm re}$
in the time domain and adjust the RIS’s phase shift with the help of Algorithm
1, and does not need any channel estimation to obtain ${\Phi_{{\rm opt}}}$.

Algorithm 1 The swarm intelligence-based optimization algorithm to find
${\Phi_{{\rm opt}}}$ for the RIS with finite-level phase-shifters.

Input: Input the number of elements on RIS: $L$, the minimum value of phase-
shifts: $\Delta\theta$, maximum iterations: ${\rm Max}$, size of swarm:
$n_{\rm begin}$ and $n_{\rm end}$, inertia weights: $\omega_{\rm begin}$ and
$\omega_{\rm end}$, acceleration constants: $c_{1}$, $c_{2}$ and $c_{3}$,
maximum of velocity: $v_{\rm max}$
Output: The phase-shift variable ${\rm\varphi_{\ell}}$ $(\ell=1,\ldots,L)$ for
each RIS element should be set.

1:Set iteration index $k=0$. Initialize the velocities and positions of $n$
particles, $\omega\leftarrow\omega_{\rm begin}-k/{\rm Max}*\left(\omega_{\rm
begin}-\omega_{\rm end}\right)$.

2:for Every praticle in swarm do

3: Obtain corresponding expectation of the amplitude of the received signal,
$E_{\rm re}$.

4:Update personal best position of each particle, global best position of the
swarm, local best position and corresponding values.

5:while $k<{\rm Max}$ do

6: for Every praticle do

7: Obtain corresponding $E_{\rm re}$.

8: if Current position is personal best position for this particle then

9: Update its personal best position.

10: if Current value of $E_{\rm re}$ is smaller than the last value then

11: Increase the periodic value: ${\rm Pe}\leftarrow{\rm Pe}+1$

12: Update global best position, local best position and corresponding values

13: Calculate the number of particles that will be deleted: $D\leftarrow
n-\left[n_{\rm begin}-k/{\rm Max}*\left(n_{\rm begin}-n_{\rm
end}\right)\right]$

14: if $D>0$ then

15: Delete $D$ particles which have the largest periodic value, ${\rm Pe}$

16: Update the number of particles: $n\leftarrow n-D$

17: for Every praticle do

18: Update $v$ with the help of (29) and (III-A3)

19: if $v>v_{\rm max}$ then

20: Update $v$: $v\leftarrow v_{\rm max}$

21: Update $x$ with the help of (31)

22: Update the inertia weight: $\omega\leftarrow\omega_{\rm begin}-k/{\rm
Max}*\left(\omega_{\rm begin}-\omega_{\rm end}\right)$

23:return In $L$-dimensional search space, the global best position is
${\rm\varphi_{\ell}}$ $(\ell=1,\ldots,L)$

###  III-B Convergence and Effectiveness Analysis

####  III-B1 Convergence Analysis of the SIO

To prove the convergence of our proposed SIO, we first derive the location
iteration function. However, when $\Delta\theta\neq 0$, it is difficult to
substitute (29) into (31) because of the rounding function in (29). To solve
this problem, we rewrite (29) as

${v_{id}}\left({k+1}\right)=\omega{v_{id}}\left(k\right)+\Delta
v\left(k\right)+c\left(k\right),$ (32)

where $c\left(k\right)$ denotes the mismatch between
${\omega{v_{id}}\left(k\right)+\Delta v\left(k\right)}$ and RIS’s practical
discrete phases. We can obtain the value of $c\left(k\right)$ by

$c\left(k\right){\rm{=}}\left[{\left({\omega{v_{id}}\left(k\right)+\Delta
v\left(k\right)}\right)/\Delta\theta}\right]\times\Delta\theta-\left({\omega{v_{id}}\left(k\right)+\Delta
v\left(k\right)}\right).$ (33)

Thus, we have
$-\frac{1}{2}\Delta\theta<c\left(k\right)<\frac{1}{2}\Delta\theta$. The
velocity of each particle is updated as

${v_{id}}\left(k+1\right)={x_{id}}\left(k+1\right)-{x_{id}}\left({k}\right).$
(34)

Let $T_{i}\triangleq c_{i}R_{i}$. Substituting (III-A3) and (32) into (34), we
can derive the location iteration function as

$\displaystyle{x_{id}}\left({k+1}\right)-\left({1-{T_{1}}-{T_{2}}-{T_{3}}+\omega}\right){x_{id}}\left(k\right)+\omega{x_{id}}\left({k-1}\right)$
$\displaystyle={T_{1}}{p_{id}}\left(k\right)+{T_{2}}{p_{gd}}\left(k\right)+{T_{3}}{p_{ld}}\left(k\right)+c\left(k\right).$
(35)

When the RIS’s elements can adjust the phase-shifts continuously, we have
$\Delta\theta=0$ and $c\left(k\right)=0$. For the discrete case, the
$c\left(k\right)=0$ can be regarded as a bounded noise. Note that the
recurrence equation of SIO is linear non-homogeneous. To calculate the
characteristic roots, we consider the secular equation as

${x^{2}}-\left({1-{T_{1}}-{T_{2}}-{T_{3}}+\omega}\right)x+\omega=0.$ (36)

The discrete linear system that has such a secular equation has been well
studied [47]. Because the latent roots of (36) can be derived easily, we can
obtain the expression of $x(k)$ and prove that the SIO algorithm can converge
to the optimal solution. Moreover, SIO suffers from a linear increase in the
search complexity with an increase in the maximum iteration number, ${\rm
Max}$. In addition, $E_{{\rm{re}}}$ is obtained $n$ times in each iteration.
Thus, the complexity of SIO is given as $\mathcal{O}(n{\rm Max})$.

####  III-B2 Effectiveness Analysis of the SIO

Let
$E_{{\rm{re}}}\triangleq{\mathbb{E}}\left[{\left|{{h_{L}}}\right|\left|{{h_{P}}}\right|\left|{{h_{F}}}\right|}\right]$
and
$E_{{\rm{opt}}}\triangleq\left|{{h_{L}}}\right|{\mathbb{E}}\left[{\left|{{h_{P}}}\right|}\right]{\mathbb{E}}\left[{\left|{\sum\limits_{\ell=1}^{L}{{g_{1}^{\left(\ell\right)}}{g_{2}^{\left(\ell\right)}}}}\right|}\right]$,
where $E_{{\rm{re}}}$ is the fitness function in our SIO algorithm which can
be easily measured in the time domain, and $E_{{\rm{opt}}}$ the maximum
achievable $E_{{\rm{re}}}$, ${\mathbb{E}}\left[\cdot\right]$ denotes the
mathematical expectation.

To investigate the effectiveness of the proposed phase optimization method, we
first derive $E_{{\rm{opt}}}$ in the following theorem

######  Theorem 1.

The maximum achievable expectation of the amplitude of the received signal,
$E_{{\rm{opt}}}$, can be derived as

$E_{{\rm{opt}}}\\!=\\!\frac{{\left|{{h_{L}}}\right|{\gamma^{2}}{A_{o}}}}{{{\gamma^{2}}\\!+\\!1}}\sum\limits_{\ell=1}^{L}\prod\limits_{\ell=1}^{2}{\frac{{m_{\ell}^{{m_{\ell}}}}}{{\Gamma\\!\left(\\!{{m_{\ell}}}\\!\right)}}}\\!\sum\limits_{{j_{\ell}}=0}^{{N_{T}}}{\frac{{K_{\ell}^{{j_{\ell}}}{d_{{\ell_{{j_{\ell}}}}}}{{\left({2\sigma_{\ell}^{2}}\right)}^{\frac{1}{2}}}}}{{{j_{\ell}}!\Gamma\\!\left(\\!{{j_{\ell}}\\!+\\!1}\\!\right)}}}\Gamma\\!\left(\\!{\frac{3}{2}\\!+\\!{j_{\ell}}}\\!\right)\\!\\!.$
(37)

###### Proof:

Please refer to Appendix A. ∎

Figure 4: The number of iterations (until $E_{{\rm re}}/E_{{\rm opt}}\geq
90\%$) versus the number of RIS’s elements, with $\Delta\theta=\pi/10$. Figure
5: The expectation error of the proposed SIO algorithm under different minimum
values of phase-shifts versus the iterations, with $L=50$. Figure 6: The
expectation error of the SIO algorithm and the algorithm proposed in [41]
versus the number of measurements of $E_{{\rm{re}}}$.

We use $E_{\rm re}/E_{\rm opt}$ to characterize the expectation accuracy and
verify the convergence and effectiveness of the proposed algorithm. In Figs. 4
and 5, we verify the convergence and effectiveness of the proposed algorithm
for different iterations, by setting $P=40$ ${\rm dBW}$, $m_{\ell,1}=3$,
$m_{\ell,2}=4$, ${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)=20$ ${\rm
dB}$, ${{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=10$ ${\rm dB}$,
$f=300$ ${\rm GHz}$, $d_{1}=30$ ${\rm m}$, $d_{2}=20$ ${\rm m}$, $N_{0}=1$
${\rm dBW}$, $K_{\ell,1}=4$, $K_{\ell,2}=5$, $\Delta_{\ell,1}=0.7$,
$\Delta_{\ell,2}=0.5$, $\kappa_{S}=\kappa_{D}=0.2$, $n_{\rm begin}=30$,
$n_{\rm end}=10$, $c_{1}=1$, $c_{2}=2$, $c_{3}=2$, $\omega_{\rm begin}=0.9$,
and $\omega_{\rm end}=0.4$. In Fig. 4, the proposed SIO is compared with
traditional PSO under the different number of reflecting elements, $L$. As it
can be observed, although the numbers of iterations for both PSO and SIO
increase linearly with the number of RIS elements, the increase rate when SIO
is used is reduced by about $34\%$ compared to that when PSO is used, which
means that the SIO is always superior to the traditional PSO at finding
optimal phase-shifts. Moreover, the red and blue shaded ranges in Fig. 4
indicate the uncertainty in the actual operation of the algorithms. We can see
that the fluctuation range of the number of iterations of the SIO (typically
less than $25$) is smaller than that of the PSO. Figure 5 depicts $E_{\rm
re}/E_{\rm opt}$ versus the amount of reflecting elements on the RIS, with
different choices of $\Delta\theta$. Although the error of expectation
increases as $\Delta\theta$ is increased, we can observe that using the RIS
with large $\Delta\theta$, e.g., $\pi/4$, can achieve the accuracy close to
that when using RIS with continuous phase-shifts. Furthermore, Fig. 6 shows
that the proposed SIO algorithm is better than the binary search tree-based
algorithm proposed in [41]. Although both the proposed SIO algorithm and
algorithm proposed in [41] need to measure the $E_{{\rm{re}}}$, the number of
reflecting elements only affects the number of measurements in the algorithm
proposed in [41]. Thus, especially when the RIS is equipped with a large
number of elements, the SIO algorithm can achieve a higher $E_{\rm re}/E_{\rm
opt}$ than that of the algorithm proposed in [41] with the same number of
measurements.

##  IV Exact Statistics of the End-to-End SNR and SNDR

Substituting $R=\sqrt{\gamma}$ into (9), we can easily obtain the PDF
expression of FTR RVs. Let
$Y=\sum\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{N}{{R_{\iota,\ell}}}}$ be
the sum of product of $N$ FTR RVs, where $R_{\iota,\ell}\sim{\cal F}{\cal
T}{\cal
R}\left({{K_{\iota,\ell}},{m_{\iota,\ell}},{\Delta_{\iota,\ell}},\sigma_{\iota,\ell}^{2}}\right)$
($\iota=1,\ldots,L$ and $\ell=1,\ldots,N$), the exact PDF and CDF of $Y$ have
been respectively derived in multivariate Fox’s $H$-function [48, eq. (A.1)]
as [41, eq. 11] and [41, eq. 12]

In Figs. 4-6, we verify that Algorithm 1 can be used to obtain the phase-
shifts of RIS’s elements which enable the received power tightly close to the
optimal received power. Then, we present the instantaneous optimal SNR and
SNDR for the case of ideal and non-ideal RF chains, respectively.

###  IV-A Non-Ideal RF Chains

To derive a generic analytical expression for the CDF of the instantaneous
SNDR of non-ideal RF chains, we first present the following lemma with the
help of (25).

######  Lemma 1.

Considering hardware impairments, we derive an integral representation for the
PDF and CDF of ${\gamma_{N}}$ in (25) assuming arbitrarily distributed
$\left|{{h_{F}}}\right|$ and $\left|{{h_{P}}}\right|$ as

$\displaystyle{f_{{\gamma_{N}}}}\left(x\right)$
$\displaystyle=\int_{0}^{{A_{o}}}{\frac{1}{y}}{f_{\left|{{h_{F}}}\right|}}\left({\frac{{\rm{1}}}{y}\sqrt{\frac{{x{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}}}\right){f_{\left|{{h_{P}}}\right|}}(y){\rm{d}}y$
$\displaystyle\times\frac{1}{{2\sqrt{x{{\left({1-x{\kappa^{2}}}\right)}^{\rm{3}}}\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}}}},$
(38)
${F_{{\gamma_{N}}}}\left(x\right)=\int_{0}^{{A_{o}}}{{F_{\left|{{h_{F}}}\right|}}}\left({\frac{{\rm{1}}}{y}\sqrt{\frac{{x{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}}}\right){f_{\left|{{h_{P}}}\right|}}(y){\rm{d}}y.$
(39)

###### Proof:

Please refer to Appendix B. ∎

With the help of Lemma 1, for special distributed $\left|{{h_{F}}}\right|$ and
$\left|{{h_{P}}}\right|$ in this paper, the statistics of instantaneous SNDR
can be obtained as the following theorem.

######  Theorem 1.

The PDF and CDF of ${{\gamma_{N}}}$ can be derived as (1) and (1), shown at
the bottom of the next page, where
$\left\\{{\left({{a_{n}}}\right)}\right\\}_{1}^{N}=\left({{a_{1}}}\right),\ldots,\left({{a_{N}}}\right)$,
and
${\upsilon_{{h_{FP}},\ell}}=\left\\{{\left({-{j_{\ell,n}},0.5}\right)}\right\\}_{1}^{\rm{2}},\left({1-{\gamma^{2}},1}\right)$
$(\ell=1,\ldots,L)$.

$\displaystyle{f_{{\gamma_{N}}}}\left(x\right)=\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{{{\gamma^{2}}}}{{\Gamma({j_{\iota,\ell}}+1)}}\frac{1}{{2x\left({1-x{\kappa^{2}}}\right)}}}}}$
$\displaystyle\times\\!\\!H_{1,0:3,{{2}};\ldots;3{{,2}}}^{0,0:1,3;\ldots;1,3}\left(\\!\\!\\!{\left.{\begin{array}[]{*{20}{c}}{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{1,\ell}}}\right)}}\\\
\vdots\\\
{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{L,\ell}}}\right)}}\end{array}}\\!\right|\begin{array}[]{*{20}{c}}{\left({0;1,\ldots,1}\right):{\upsilon_{{h_{FP}},1}};\ldots;{\upsilon_{{h_{FP}},L}}}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots;\left({0,1}\right),\left({-{\gamma^{2}},1}\right)}\end{array}}\\!\\!\\!\\!\right)$
(45)

$\displaystyle{F_{{\gamma_{N}}}}\left(x\right)=\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{1}{{\Gamma({j_{\iota,\ell}}+1)}}}}}{\gamma^{2}}$
$\displaystyle\times\\!\\!H_{1,0:3,{{2}};\ldots;3{{,2}}}^{0,0:1,3;\ldots;1,3}\left(\\!\\!\\!{\left.{\begin{array}[]{*{20}{c}}{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{1,\ell}}}\right)}}\\\
\vdots\\\
{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{L,\ell}}}\right)}}\end{array}}\right|\begin{array}[]{*{20}{c}}{\left({1;1,\ldots,1}\right):{\upsilon_{{h_{FP}},1}};\ldots;{\upsilon_{{h_{FP}},L}}}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots;\left({0,1}\right),\left({-{\gamma^{2}},1}\right)}\end{array}}\\!\\!\right)$
(51)

###### Proof:

Please refer to Appendix C. ∎

######  Remark 1.

The multivariate Fox’s $H$-function cannot be evaluated directly by popular
mathematical packages such as MATLAB or Mathematica. However, many
implementations have been proposed [49, 50]. In the following, we use an
efficient and accurate Python implementation [50] which can provide accurate
results of multivariate Fox’s $H$-function quickly by estimating the
multivariate integrals with simple rectangle quadrature. For example, the
execution time for calculating ternary Fox’s $H$-function is less than a few
seconds.

###  IV-B Ideal RF Chains

######  Theorem 2.

The PDF and CDF of ${{\gamma_{I}}}$ can be derived as (2) and (2), shown at
the bottom of the next page.

$\displaystyle{f_{{\gamma_{I}}}}\left(x\right)=\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{{{\gamma^{2}}}}{{\Gamma({j_{\iota,\ell}}+1)}}\frac{{1}}{2x}}}}$
$\displaystyle\times
H_{1,0:3,{{2}};\ldots;3{\rm{,2}}}^{0,0:1,3;\ldots;1,3}\left({\left.{\begin{array}[]{*{20}{c}}{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{1,\ell}}}\right)}}\\\
\vdots\\\
{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{L,\ell}}}\right)}}\end{array}}\right|\begin{array}[]{*{20}{c}}{\left({0;1,\ldots,1}\right):{\upsilon_{{h_{FP}},1}};\ldots;{\upsilon_{{h_{FP}},L}}}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots;\left({0,1}\right),\left({-{\gamma^{2}},1}\right)}\end{array}}\right)$
(47)

$\displaystyle{F_{{\gamma_{I}}}}\left(x\right)=\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{1}{{\Gamma({j_{\iota,\ell}}+1)}}}}}{\gamma^{2}}$
$\displaystyle\times
H_{1,0:3,{{2}};\ldots;3{{,2}}}^{0,0:1,3;\ldots;1,3}\left({\left.{\begin{array}[]{*{20}{c}}{{A_{\rm{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{1,\ell}}}\right)}}\\\
\vdots\\\
{{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{L,\ell}}}\right)}}\end{array}}\right|\begin{array}[]{*{20}{c}}{\left({1;1,\ldots,1}\right):{\upsilon_{{h_{FP}},1}};\ldots;{\upsilon_{{h_{FP}},L}}}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots;\left({0,1}\right),\left({-{\gamma^{2}},1}\right)}\end{array}}\right)$
(53)

###### Proof:

Substituting $\kappa^{2}=0$ into (1) and (1), we obtain the PDF and CDF of
${{\gamma_{N}}}$ after performing some mathematical manipulations, which
completes the proof. ∎

The statistics of $\gamma_{I}$ and $\gamma_{I}$ that we derived in Theorem 2
and Theorem 3 can be used to obtain exact closed-form expressions for key
performance metrics.

##  V Asymptotic Analysis

Although the previously derived analytical results are obtained in closed-
form, they may not provide many insights as to the factors affecting system
performance. In the following, the asymptotic performance analysis that
becomes tight in high-SNDR and $L$ regimes is presented.

###  V-A Asymptotic analysis in the high SNDR regime

######  Proposition 1.

A simple approximation of the CDF of $\gamma_{N}$ in the high-SNDR regime can
be derived as

$\displaystyle
F_{{\gamma_{N}}}^{HS}\left(x\right)\\!=\\!\\!\\!\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{{N_{T}}}{\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{{N_{T}}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}{\gamma^{2}}}}{{{j_{\iota,\ell}}!\Gamma\\!\left(\\!{{m_{\iota,\ell}}}\\!\right)\\!\Gamma\\!\left(\\!{{j_{\iota,\ell}}\\!+\\!1}\\!\right)}}}}}}$
$\displaystyle\times\\!\prod\limits_{\iota=1}^{L}\frac{{\mathop{\lim}\limits_{{\varsigma_{\iota}}\to-{\varpi_{\iota}}}\left({{\varsigma_{\iota}}+{\varpi_{\iota}}}\right)\prod\limits_{\ell=1}^{2}{\Gamma\\!\left({1\\!+\\!{j_{\iota,\ell}}\\!-\\!\frac{1}{2}{\varpi_{\iota}}}\right)}\Gamma\\!\left({{\varpi_{\iota}}}\right)\Gamma\\!\left({{\gamma^{2}}\\!-\\!{\varpi_{\iota}}}\right)}}{{\Gamma\\!\left({1\\!+\\!\sum\limits_{\iota{\text{
=
1}}}^{L}{{\varpi_{\iota}}}}\right)\\!\Gamma\\!\left({{\gamma^{2}}\\!-\\!{\varpi_{\iota}}\\!+\\!1}\right)}}$
$\displaystyle\times\\!{{\left(\\!\\!{{A_{\text{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1\\!-\\!x{\kappa^{2}}}\right)}}{{x{N_{o}}}}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\\!\right)}^{-{\varpi_{\iota}}}},$
(44)

where
${\varpi_{\iota}}=\min\left\\{{{\gamma^{2}},2+2{j_{\iota,1}},2+2{j_{\iota,2}}}\right\\}$,
and

$\displaystyle\mathop{\lim}\limits_{{\varsigma_{\iota}}\to-{\varpi_{\iota}}}\left({{\varsigma_{\iota}}+{\varpi_{\iota}}}\right)\prod\limits_{\ell=1}^{2}{\Gamma\left({1+{j_{\iota,\ell}}-\frac{1}{2}{\varpi_{\iota}}}\right)}\Gamma\left({{\varpi_{\iota}}}\right)\Gamma\\!\left({{\gamma^{2}}-{\varpi_{\iota}}}\right)$
$\displaystyle=\left\\{\begin{gathered}\prod\limits_{\ell=1}^{2}{\Gamma\\!\left({1+{j_{\iota,\ell}}-\frac{1}{2}{\gamma^{2}}}\right)}\Gamma\\!\left({{\gamma^{2}}}\right),{\varpi_{\iota}}={\gamma^{2}}\hfill\\\
\Gamma\\!\left({{j_{\iota,2}}-{j_{\iota,1}}}\right)\Gamma\\!\left({2+2{j_{\iota,1}}}\right)\Gamma\\!\left({{\gamma^{2}}-2-2{j_{\iota,1}}}\right),{\varpi_{\iota}}=2+2{j_{\iota,1}}\hfill\\\
\Gamma\\!\left({{j_{\iota,1}}-{j_{\iota,2}}}\right)\Gamma\\!\left({2+2{j_{\iota,2}}}\right)\Gamma\\!\left({{\gamma^{2}}-2-2{j_{\iota,2}}}\right),{\varpi_{\iota}}=2+2{j_{\iota,2}}\hfill\\\
\end{gathered}\right.$ (49)

###### Proof:

With the help of the definition of multivariate Fox’s $H$-function [48, eq.
(A-1)], eq. (1) can be expressed in terms of multiple Mellin-Barnes type
contour integrals. When $SNDR\to\infty$, we have $P\to\infty$. The Mellin-
Barnes integrals can be approximated by evaluating the residue at the minimum
pole on the left-hand side of $d{\varsigma_{\iota}}$ as [48, Theorem 1.11].
Therefore, the expression in (1) is derived, which completes the proof. ∎

Note that the CDF of ${{\gamma_{I}}}$ in the high-SNDR regime can be derived
by substituting $\kappa^{2}=0$ into (1).

###  V-B Asymptotic analysis in the high $L$ regime

Considering the RIS typically has hundreds of elements, we next derive closed-
form statistics of the RV
${{h_{F}}}=\sum\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{{R_{\iota,\ell}}}}$
in high-SNDR and $L$ regimes. Employing central limit theorem (CLT) as [51],
we have
${{h_{F}}}\sim\mathcal{N}\left({L{\mathbb{E}}\left[{h_{F}}\right],L{\mathbb{V}}\left[{h_{F}}\right]}\right)$,
where ${\mathbb{V}}(\cdot)$ denotes the variance of an RV. Considering that
${\mathbb{V}}(\cdot)$=${\mathbb{E}}(\cdot^{2})$-${\mathbb{E}}^{2}(\cdot)$, we
can obtain the values of ${\mathbb{E}}\left[{h_{F}}\right]$ and
${\mathbb{V}}\left[{{h_{F}}}\right]$ with the help of [41, eq. (6)].
Therefore, the CDF of ${h_{F}}$ in the high $L$ regime can be expressed as
[52]

$F_{{h_{F}}}^{HL}\\!\left(x\right)\\!=\\!\frac{1}{2}\\!\left(\\!{1\\!+\\!{\rm
erf}\left(\\!{\frac{{x-L{\mathbb{E}}\left[{h_{F}}\right]}}{{\sqrt{2L{\mathbb{V}}\left[{h_{F}}\right]}}}}\\!\right)}\\!\right)\\!=\\!\frac{1}{2}{{\rm
erfc}\\!\left(\\!-{\frac{{x-L{\mathbb{E}}\left[{h_{F}}\right]}}{{\sqrt{2L{\mathbb{V}}\left[{h_{F}}\right]}}}}\\!\right)},$
(50)

where ${\rm erf}(\cdot)$ is the error function [29, eq. (8.250.1)], and ${\rm
erfc}(\cdot)$ is the complementary error function [29, eq. (8.250.4)]. When
$L\to\infty$, the CDF of $\gamma_{N}$ can be derived with the help of (50), by
following the similar method as shown in Appendix C. However, note that ${\rm
erf}\left(\cdot\right)$ or ${\rm erfc}\left(\cdot\right)$ in (50) is a non-
elementary function, which is defined in the form of a definite integral.
Thus, it is hard to obtain new insights. Therefore, we use a five-first-moment
(5FM) based method to obtain the approximate CDF expression.

######  Proposition 2.

An approximation of the CDF of $\gamma_{N}$ in the high-$L$ regime can be
expressed as

$F_{{\gamma_{N}}}^{HL}\left(x\right)={a_{1}}{a_{2}}{\gamma^{2}}G_{3,4}^{3,1}\left({\left.{\frac{\Upsilon}{{{A_{o}}{a_{2}}}}}\right|\begin{array}[]{*{20}{c}}{1,{a_{3}}+1,{\gamma^{2}}+1}\\\
{{a_{4}}+1,{a_{5}}+1,{\gamma^{2}},0}\end{array}}\right),$ (51)

where
$\Upsilon\triangleq\sqrt{\frac{{x{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}}$,
$G_{\cdot,\cdot}^{\cdot,\cdot}\left({\cdot\left|\cdot\right.}\right)$ is the
Meijer’s $G$ function [29, eq. (9.30)], ${\Omega_{i}}$ denotes the $i_{\rm
th}$ moment of $h_{F}$, ${\varphi_{i}}={\Omega_{i}}/{\Omega_{i-1}}$,
${a_{3}}=\left({4{\varphi_{4}}-9{\varphi_{3}}+6{\varphi_{2}}-{\Omega_{1}}}\right){\rm{/}}\left({3{\varphi_{3}}-{\varphi_{4}}-3{\varphi_{2}}+{\Omega_{1}}}\right)$,
${a_{2}}=\frac{{{a_{3}}}}{2}\left({{\varphi_{4}}-2{\varphi_{3}}+{\varphi_{2}}}\right)+2{\varphi_{4}}-3{\varphi_{3}}+{\varphi_{2}}$,
${a_{6}}=\left({{a_{3}}\left({{\varphi_{2}}-{\Omega_{1}}}\right)+2{\varphi_{2}}-{\Omega_{1}}}\right)/{a_{2}}-3$,
${a_{7}}=\sqrt{{{\left({{a_{6}}+2}\right)}^{2}}-4{\Omega_{1}}\left({{a_{3}}+1}\right)/{a_{2}}}$,
${a_{4}}=\left({{a_{6}}+{a_{7}}}\right){\rm{/2}}$,
${a_{5}}=\left({{a_{6}}-{a_{7}}}\right){\rm{/2}}$, and
${a_{1}}={\Gamma\left({{a_{3}}+1}\right)/\left({{a_{2}}\Gamma\left({{a_{4}}+1}\right)\Gamma\left({{a_{5}}+1}\right)}\right)}$.

###### Proof:

Please refer to Appendix D. ∎

By substituting $\kappa^{2}=0$ in (51), after performing some mathematical
manipulations, we can obtain the CDF of ${{\gamma_{I}}}$ in the high $L$
regime.

###  V-C Verification and Insights

Figure 7: Outage probability versus the transmit power with $L=20$,
$A_{o}=0.01$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $m_{\ell,1}=5$, $m_{\ell,2}=7$,
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=10$
${\rm dB}$, $\Delta_{\ell,1}=0.6$, $\Delta_{\ell,2}=0.4$, $N_{0}=1$ ${\rm
dBW}$, $\gamma_{\rm th}=0.5$ ${\rm dB}$, and $\kappa_{S}=\kappa_{D}=0.1$ in
three different cases, where case 1 is $\left|{{h_{L}}}\right|=0.1$ and
$\gamma=0.5$, case 2 is $\left|{{h_{L}}}\right|=1$ and $\gamma=0.5$, and case
3 is $\left|{{h_{L}}}\right|=0.1$ and $\gamma=0.8$.

With the help of the approximate CDF expressions in Propositions 1 and 2, the
difficulty of deriving performance metrics for RIS-aided THz communications
systems can be avoided. For example, the outage probability $P$ can be defined
as the probability that the received SNR per signal falls below a given
threshold $\gamma_{\rm th}$, which can be directly obtained using the derived
CDF expressions. Figure 7 shows the outage probability versus the transmit
power. We can observe that the analytical results ((1) in Theorem 1) match
well with the high-$L$ approximate results obtained by the two methods, i.e.,
CLT and 5FM, in Proposition 2 with $L=20$, which means that the high-$L$
approximate expressions are useful because the RIS can be equipped with
hundreds of reflecting elements in the practical system. Note that when the
number of reflecting elements is small, the CLT method is not accurate, but
the 5FM method can be used because the parameter ${\Omega_{i}}$, which denotes
the $i_{\rm th}$ moment of small-scale fading coefficient, contains the impact
of the number of elements. Furthermore, the accuracy of the approximate CDF
expressions derived in Proposition 1 is verified in Fig. 8 in the following
numerical analysis section. We can observe that at high-transmitting power
regime, (1) approximates well the exact analytical expression.

In addition to simplifying the performance analysis, we can observe several
insights from the approximate expressions. For example, similar to the
diversity order [53], the high-SNR slope is a key performance indicator that
can explicitly capture the impact of channel parameters on OP performance,
which is defined as [54]

${G_{d}}=-\mathop{\lim}\limits_{P\to\infty}\frac{{\log\left({O{P_{{\rm{appr}}}}}\right)}}{{\log\left(P\right)}}.$
(52)

With the help of the accurate 5FM-based approximation (51), we derive the
simple approximate OP in high-SNDR regime as

$\displaystyle{{\rm OP}_{{\rm
appr}}}\left(x\right)\\!=\\!\frac{{\mathop{\lim}\limits_{s\to{a_{\min}}}\\!\left({{a_{\min}}-s}\right){\Gamma^{-1}}\left({{\gamma^{2}}\\!+\\!1\\!-\\!{a_{\min}}}\right)}}{{\Gamma\\!\left({{a_{3}}\\!+\\!1\\!-\\!{a_{\min}}}\right)\Gamma\\!\left({{a_{\min}}\\!+\\!1}\right)}}{\left(\\!{\frac{\Upsilon}{{{A_{o}}{a_{2}}}}}\\!\right)^{{a_{\min}}}}$
$\displaystyle\times{a_{1}}{a_{2}}{\gamma^{2}}\Gamma\left({{a_{4}}+1-s}\right)\Gamma\left({{a_{5}}+1-s}\right)\Gamma\left({{\gamma^{2}}-s}\right)\Gamma\left(s\right),$
(53)

where ${a_{\min}}=\min\left\\{{{a_{5}}+1,{a_{4}}+1,{\gamma^{2}}}\right\\}$ and

$\displaystyle\mathop{\lim}\limits_{s\to{a_{\min}}}\left({{a_{\min}}-s}\right)\Gamma\left({{a_{4}}+1-s}\right)\Gamma\left({{a_{5}}+1-s}\right)\Gamma\left({{\gamma^{2}}-s}\right)\Gamma\left(s\right)$
$\displaystyle=\\!\left\\{\\!\\!\\!{\begin{array}[]{*{20}{l}}{\Gamma\\!\left({{a_{5}}-{a_{4}}}\right)\Gamma\\!\left({{\gamma^{2}}-{a_{4}}-1}\right)\Gamma\\!\left({{a_{4}}\\!+\\!1}\right),}&\\!\\!{{a_{\min}}={a_{4}}+1}\\\
{\Gamma\\!\left({{a_{4}}-{a_{5}}}\right)\Gamma\\!\left({{\gamma^{2}}-{a_{5}}-1}\right)\Gamma\\!\left({{a_{5}}\\!+\\!1}\right),}&\\!\\!{{a_{\min}}={a_{5}}+1}\\\
{\Gamma\\!\left({{a_{4}}\\!+\\!1-{\gamma^{2}}}\right)\Gamma\\!\left({{a_{5}}\\!+\\!1-{\gamma^{2}}}\right)\Gamma\\!\left({{\gamma^{2}}}\right),}&\\!\\!{{a_{\min}}={\gamma^{2}}}\end{array}}\right.$
(57)

###### Proof:

In the high-SNDR regime, e.g., when the transmit power is high, we have
$\Upsilon\to 0$. The Mellin-Barnes integral over $\mathcal{L}_{s}$ can be
approximated by evaluating the residue at the minimum pole on the left-hand
side of $s$ as [48, Theorem 1.11]. With the help of [29, eq. (8.331.1)] and
[29, eq. (9.113)], we can obtain (V-C), which completes the proof. ∎

From (V-C) and (52), we can obtain that
${G_{d}}={a_{\min}}/2=\min\left\\{{{a_{5}}+1,{a_{4}}+1,{\gamma^{2}}}\right\\}/2$.
Recall that ${{a_{4}}}$ and ${{a_{5}}}$ are determined by the moments of the
small-scale fading coefficient, and ${\gamma^{2}}$ is defined as the ratio of
THz signal’s equivalent beam-width and the variance of the pointing error
displacement at Destination. Considering that ${G_{d}}$ indicates how much the
OP decreases with the increase of the transmit power, we can conclude that
when the pointing error is large, the key to improve ${G_{d}}$ lies in
enhancing the beamforming to increase $\gamma$. when the pointing error is
small, ${G_{d}}$ is mainly determined by the small-scale fading. This
observation is verified in Fig. 7: the blue and pink lines, i.e., case 1 and
case 2, respectively, use the same value of ${G_{d}}$, so they have the same
slope in the high SNDR region. However, in the simulation of the red line,
i.e., case 3, the ${G_{d}}$ is increased from $0.5$ to $0.8$, so that the
maximum slope is achieved.

##  VI Performance Analysis

###  VI-A Outage Probability

The outage probability of considered THz communications can be obtained as
${P_{{\rm{out}}}}=P\left({\gamma_{\ell}<{\gamma_{{\rm{th}}}}}\right){{=}}F_{\gamma_{\ell}}\left({{\gamma_{{\rm{th}}}}}\right)$,
where $\ell\in N,I$. For the case of non-ideal RF chains, the outage
probability of the RIS-aided system can be directly obtained by using (1).
Furthermore, we can obtain the outage probability for the case of ideal RF
chains with the help of (2). Moreover, in high SNDR and $L$ regimes, the
outage probability can be derived with the help of (1) and (51), respectively.
Thus, we can observe that an RIS equipped with more reflecting elements will
also achieve a low outage probability.

###  VI-B Ergodic Capacity

The ergodic capacity under optimal rate adaptation with constant transmit
power is given by [55, eq. (29)] as
$C={\mathbb{E}}\left[{{{\log}_{2}}\left({1+{\gamma_{\ell}}}\right)}\right]$,
where $\ell\in N,I$.

####  VI-B1 Non-Ideal RF Chains

For the case of non-ideal RF chains, an exact expression of the ergodic
capacity is hard to be derived due to the complexity of (1). However, the
following proposition provides a closed-form upper bound of the ergodic
capacity, which is more generic, and practically more useful and reliable.
Moreover, in Fig. 10, we show that the derived capacity upper bound is close
to the exact simulation results.

######  Proposition 3.

The effective capacity of the RIS-aid THz communications for the case of non-
ideal RF chains is bounded by

$C_{N}\leq{\log_{2}}\left({1+\frac{{P{{\left|{{h_{L}}}\right|}^{2}}{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]}}{{P{\kappa^{2}}{{\left|{{h_{L}}}\right|}^{2}}{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]+{N_{o}}}}}\right),$
(58)

where ${\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]$ is derived
as (3), shown at the bottom of the next page.

$\displaystyle{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]=\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}{\sum\limits_{{j_{1,{\rm{2}}}},\ldots,{j_{L,{\rm{2}}}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{{{\gamma^{2}}}}{{{e^{2}}\Gamma({j_{\iota,\ell}}+1)}}}}}}$
$\displaystyle\times\\!\\!\\!H_{2,0:3,2;\ldots;3,2}^{0,1:1,3;\ldots;1,3}\left(\\!\\!\\!{\left.{\begin{array}[]{*{20}{c}}{e{A_{{0}}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{1,\ell}}}\right)}}\\\
\vdots\\\
{e{A_{{0}}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{L,\ell}}}\right)}}\end{array}}\\!\\!\right|\\!\\!\begin{array}[]{*{20}{c}}{\left({-1;-1,\ldots,-1}\right),\left({0;-1,\ldots\\!,-1}\right):{\upsilon_{FP,1}};\ldots\\!;{\upsilon_{FP,L}}}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots\\!;\left({0,1}\right),\left({-{\gamma^{2}},1}\right)}\end{array}}\\!\\!\\!\right)$
(64)

###### Proof:

Please refer to Appendix E. ∎

####  VI-B2 Ideal RF Chains

Exact ergodic capacity can be derived for the case of ideal RF chains.

######  Proposition 4.

Effective capacity of the RIS-aid system for the case of ideal RF chains can
be derived as (4), shown at the bottom of the next page.

$\displaystyle C{{=}}$
$\displaystyle\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}{\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{{{\gamma^{2}}}}{{\Gamma({j_{\iota,\ell}}+1)}}\frac{1}{2}}}}}$
$\displaystyle\times\\!\\!\\!H_{2,0:3,{{2}};\ldots;3{{,2,2,2}}}^{0,1:1,3;\ldots;1,3,2,1}\left(\\!\\!{\left.{\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{{\eta_{1}}}\\\
\vdots\\\ {{\eta_{L}}}\end{array}}\\\
{\frac{1}{e}}\end{array}}\\!\\!\right|\\!\\!\begin{array}[]{*{20}{c}}{\left({1;-0.5,\ldots,-0.5,-1}\right),\left({0;1,\ldots,1,0}\right):{\upsilon_{{h_{FP}},1}};\ldots;{\upsilon_{{h_{FP}},L}};\left({0,1}\right)}\\\
{-:\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\ldots;\left({0,1}\right),\left({-{\gamma^{2}},1}\right);\left({0,1}\right),\left({0,1}\right);\left({1,1}\right)}\end{array}}\\!\\!\right)$
(67)

where
${\eta_{\varsigma}}{{=}}{A_{{0}}}\sqrt{e\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{\varsigma,\ell}}}\right)}$.

###### Proof:

Please refer to Appendix F. ∎

######  Proposition 5.

Effective capacity of the RIS-aid system for the case of ideal RF chains are
bounded by

$C_{I}\leq{\log_{2}}\left({1+\frac{{P{{\left|{{h_{L}}}\right|}^{2}}{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]}}{N_{o}}}\right).$
(61)

###### Proof:

With the help of Jensen’s inequality [56], the ergodic capacity can be upper-
bounded as
$C_{I}\leq{\log_{2}}\left({1+{\mathbb{E}}\left[{{\gamma_{I}}}\right]}\right)$.
With the help of (26), we obtain (61) to complete the proof. ∎

From (58), (4), and (61), as expected, we can observe that the multipath
parameters have a positive effect on the ergodic capacity of RIS-aid THZ
communications.

##  VII Numerical Results And Discussions

In this section, we present numerical results to investigate the joint impacts
of misalignment, multipath fading, channel conditions, and transceivers
hardware imperfections on the outage and ergodic capacity performance of the
RIS-aided THz communications. Following the practical parameters in [19, 57,
58, 59, 2], we set the parameters as follows: $T=27$ ${}^{\circ}C$, $p=101325$
${\rm Pa}$, $c=3\times 10^{8}$ ${\rm m/s}$, $w_{d}/a=6$, $G_{t}=G_{r}=40$
${\rm dBi}$, $\sigma_{S}=0.01$ ${\rm m}$, and the number of Monte Carlo
iterations is $10^{6}$.

Figure 8: Outage probability versus the transmit power with $d_{2}=20$ ${\rm
m}$, $N_{0}=1$ ${\rm dBW}$, $L=40$, $m_{\ell,1}=5$, $m_{\ell,2}=7$,
$K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$, $\Delta_{\ell,2}=0.4$,
$\kappa_{S}=\kappa_{D}=0.1$, and
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=20$
${\rm dB}$. Figure 9: OP versus the threshold, $\gamma_{\rm th}$, with $P=30$
${\rm dBW}$, $m_{\ell,1}=m_{1}$, $m_{\ell,2}=m_{2}$,
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=20$
${\rm dB}$, $f=300$ ${\rm GHz}$, $d_{1}=10$ ${\rm m}$, $d_{2}=10$ ${\rm m}$,
$N_{0}=1$ ${\rm dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$,
$\Delta_{\ell,2}=0.4$, and $\kappa_{S}=\kappa_{D}=0.1$.

Figure 8 illustrates the outage probability performance versus the transmit
power with $d_{2}=20$ ${\rm m}$, $N_{0}=0$ ${\rm dBW}$, $L=40$,
$m_{\ell,1}=5$, $m_{\ell,2}=7$, $K_{\ell,1}=5$, $K_{\ell,2}=6$,
$\Delta_{\ell,1}=0.6$, $\Delta_{\ell,2}=0.4$, $\kappa_{S}=\kappa_{D}=0.1$, and
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=20$
${\rm dB}$ $(\ell=1,\ldots,L)$. It can be observed that the approximate
expressions match the exact results well in the high-SNDR regime. For a given
transmission distance, we can observe that the outage probability increases as
the relative humidity increases. Furthermore, the increase in the transmission
distance of the signal can cause a significant increase in the OP. Comparing
the impact of the propagation loss and the molecular absorption loss, we
observe that the propagation loss has a more adverse impact on the system
performance than that of the molecular absorption loss. The reason is that the
outage probability degradation caused by increasing $d_{1}$ is severer than
the corresponding degradation caused by increasing $\phi_{1}$. Moreover, the
outage probability increases as the frequency increases. When the frequency
increases from $340$ ${\rm GHz}$ to $380$ ${\rm GHz}$, the outage probability
decreases by $58\%$, which is faster than that when the frequency increases
from $300$ ${\rm GHz}$ to $340$ ${\rm GHz}$ ($34\%$). The reason is that a
high molecular absorption exists in $380$ ${\rm GHz}$ because $380$ ${\rm
GHz}$ is one of the water vapor resonance frequencies [60]. Thus, it is
significant for the THz communications system designers to consider the
variation of the environmental conditions and the molecular absorption
phenomenon.

Figure 9 plots the outage probability versus the threshold, $\gamma_{\rm th}$,
with $P=30$ ${\rm dBW}$, $m_{\ell,1}=m_{1}$, $m_{\ell,2}=m_{2}$,
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=20$
${\rm dB}$, $f=300$ ${\rm GHz}$, $d_{1}=10$ ${\rm m}$, $d_{2}=20$ ${\rm m}$,
$N_{0}=1$ ${\rm dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$,
$\Delta_{\ell,2}=0.4$, and $\kappa_{S}=\kappa_{D}=0.1$. We can observe that
the outage probability increases by $72\%$ when $\gamma_{\rm th}$ increases by
$3$ ${\rm dB}$. Moreover, Fig. 9 shows the effects of different channel
conditions. The outage probability decreases as shaping parameters $m_{1}$ and
$m_{2}$ increase because of the better channel conditions. Furthermore, the
number of RIS’s reflecting elements, $L$, has a significant impact on the
performance of the RIS-aided THz communications. For example, when $L$
increases from $40$ to $60$, the outage probability decreases by $97\%$. The
reason is that the RIS can effectively change the phase of reflected signals
to enhance the received signal quality, enabling the desired THz wireless
channel to exhibit a perfect line-of-sight.

Figure 10: Ergodic capacity versus the transmit power with $m_{\ell,1}=5$,
$m_{\ell,2}=7$, $f=360$ ${\rm GHz}$, $d_{1}=8$ ${\rm m}$, $d_{2}=6$ ${\rm m}$,
$N_{0}=0$ ${\rm dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$,
$\Delta_{\ell,2}=0.4$, $L=40$, and
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=25$
${\rm dB}$ $(\ell=1,\ldots,L)$. Figure 11: Ergodic capacity versus the number
of RIS’s reflecting elements with $\kappa_{S}=\kappa_{D}=0$, $m_{\ell,1}=5$,
$m_{\ell,2}=7$, $d_{1}=8$ ${\rm m}$, $d_{2}=6$ ${\rm m}$, $N_{0}=1$ ${\rm
dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$, and
$\Delta_{\ell,2}=0.4$.

Figure 10 depicts the ergodic capacity versus the transmit power with
$m_{\ell,1}=5$, $m_{\ell,2}=7$, $f=360$ ${\rm GHz}$, $d_{1}=8$ ${\rm m}$,
$d_{2}=6$ ${\rm m}$, $N_{0}=1$ ${\rm dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$,
$\Delta_{\ell,1}=0.6$, $\Delta_{\ell,2}=0.4$, $L=40$, and
${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)={{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)=25$
$(\ell=1,\ldots,L)$ for different $\kappa_{S}=\kappa_{D}=\kappa$. We can
observe that the ergodic capacity increases when the transmission power
increases. Thus, increasing the transmission power can be regarded as a direct
measure to reduce the impact of hardware imperfections. Furthermore, Fig. 10
reveals that the hardware imperfections have a small effect on the ergodic
capacity in the low transmit power regime, while having a more adverse impact
in the high transmit power regime. Therefore, the efficient method to increase
the ergodic capacity is to decrease the hardware imperfections when the
transmit power is high. In addition, we observe that the derived capacity
upper bound is close to the Monte Carlo simulation results.

Figure 11 illustrates the ergodic capacity versus the number of RIS’s
reflecting elements with ${{2}}{\sigma_{\ell,1}}\left({1+{K_{\ell,1}}}\right)$
$=$ ${{2}}{\sigma_{\ell,2}}\left({1+{K_{\ell,2}}}\right)$ $=$
${{2}}{\sigma}\left({1+{K}}\right)$, $\kappa_{S}=\kappa_{D}=0$,
$m_{\ell,1}=5$, $m_{\ell,2}=7$, $d_{1}=8$ ${\rm m}$, $d_{2}=6$ ${\rm m}$,
$N_{0}=1$ ${\rm dBW}$, $K_{\ell,1}=5$, $K_{\ell,2}=6$, $\Delta_{\ell,1}=0.6$,
and $\Delta_{\ell,2}=0.4$ for different $f$ and
${{2}}{\sigma}\left({1+{K}}\right)$. From Fig. 11, we can observe that the
ergodic capacity increases by $25\%$ when the number of reflecting elements
increases from $40$ to $80$. This is due to the reason that more elements let
RIS have superior capability in manipulating electromagnetic waves. Thus,
equipping RIS with more reflecting elements can significantly improve the
performance of wireless communications systems. Furthermore, the ergodic
capacity decreases when the operating frequency is changed from $300$ GHz to
$380$ GHz. Thus, it becomes evident that the system’s performance will be
significantly influenced by the selection of operating frequency, which also
verifies the conclusion obtained from Fig. 8.

##  VIII Conclusions

We investigated an RIS-aided THz communications system and presented exact
expressions of the end-to-end SNR and SNDR. Based on two THz signal
measurement experiments conducted in a variety of different scenarios, we
proved that the FTR distribution performs a good fit for modeling the small-
scale amplitude fading of THz signals. Furthermore, an SIO algorithm was
proposed to find the optimal phase-shifts design at RIS under discrete phase-
shift constraints. The impact of the different number of discrete phase-shift
levels was investigated. Moreover, important performance metrics, including
the outage probability and ergodic capacity, are derived in terms of the
multivariate Fox’s $H$-function. To obtain more insights, the asymptotic
analysis was presented when the SNDR or the number of RIS’s elements is large.
The tight upper bound expressions of ergodic capacity for ideal and non-ideal
RF chains are also derived. The numerical results confirmed that the
propagation loss is more detrimental to the performance of the RIS-aided THz
communications system than the molecular absorption loss, and the RIS can
bring a significant performance gain for the proposed THz communications.

##  Appendix A Proof of Theorem 1

With the help of [41, eq. (6)],
${\mathbb{E}}\left[{\left|{\sum\limits_{\ell=1}^{L}{{g_{1}^{\left(\ell\right)}}{g_{2}^{\left(\ell\right)}}}}\right|}\right]$
can be expressed as

${\mathbb{E}}\\!\left[{\left|{\sum\limits_{\ell=1}^{L}{{g_{1}^{\left(\ell\right)}}{g_{2}^{\left(\ell\right)}}}}\right|}\right]\\!\\!=\\!\sum\limits_{\ell=1}^{L}\prod\limits_{\ell=1}^{2}{\frac{{m_{\ell}^{{m_{\ell}}}}}{{\Gamma\\!\left({{m_{\ell}}}\right)}}}\sum\limits_{{j_{\ell}}=0}^{{N_{T}}}{\frac{{K_{\ell}^{{j_{\ell}}}{d_{{\ell_{{j_{\ell}}}}}}{{\left({2\sigma_{\ell}^{2}}\right)}^{\frac{1}{2}}}}}{{{j_{\ell}}!\Gamma\\!\left({{j_{\ell}}\\!+\\!1}\right)}}}\Gamma\\!\left(\\!{\frac{3}{2}\\!+\\!{j_{\ell}}}\\!\right)\\!\\!.$
(A-1)

The $s_{\rm th}$ moment of ${\left|{{h_{P}}}\right|}$ is defined as
$E\left[{{{\left|{{h_{P}}}\right|}^{s}}}\right]{\rm{=}}\int_{0}^{{A_{o}}}{{x^{s}}{f_{{h_{p}}}}(x){\rm
d}x}$. With the help of (16), we obtain

$E\left[{{{\left|{{h_{P}}}\right|}^{s}}}\right]{{=}}\frac{{{\gamma^{2}}}}{{A_{o}^{{\gamma^{2}}}}}\int_{0}^{{A_{o}}}{{x^{{\gamma^{2}}+s-1}}dx}=\frac{{{\gamma^{2}}}}{{{\gamma^{2}}+s}}{A_{o}}^{s}.$
(A-2)

By substituting $s=1$ into (A-2), $E\left[{{{\left|{{h_{P}}}\right|}}}\right]$
is derived. Thus, we obtain (37) which completes the proof.

##  Appendix B Proof of Lemma 1

###  B-A Proof of PDF

With the help of (25), we obtain
$\frac{{{1}}}{{{\gamma_{N}}}}{{=}}{\kappa^{2}}{{+}}\frac{{{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P}}\frac{{\rm{1}}}{{{{\left|{{h_{FP}}}\right|}^{2}}}}$,
where
$\left|{{h_{FP}}}\right|\triangleq\left|{{h_{F}}}\right|\times\left|{{h_{P}}}\right|$.
Let us define that
$V\triangleq\frac{{\rm{1}}}{{{{\left|{{h_{FP}}}\right|}^{2}}}}$,
$W\triangleq\frac{{{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P}}V=\frac{{\rm{1}}}{{{\gamma_{N}}}}-{\kappa^{2}}$.
The PDF of $W$ can be respectively expressed as

$\displaystyle{f_{W}}\left(w\right)=$
$\displaystyle\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}{f_{V}}\left({\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}w}\right)$
$\displaystyle=$
$\displaystyle\frac{{-1}}{{2{w^{2}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}w}}}}}{f_{\left|{{h_{FP}}}\right|}}\left({\sqrt{\frac{{{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}Pw}}}}\right).$
(B-1)

By performing the change of variables
${\gamma_{N}}=\frac{1}{{W+{\kappa^{2}}}}$, ${f_{{\gamma_{N}}}}$ can be
rewritten as

${f_{{\gamma_{N}}}}\\!\left(x\right)\\!=\\!\frac{1}{{2{{\left({1\\!-\\!x{\kappa^{2}}}\right)}^{2}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}\frac{{1-x{\kappa^{2}}}}{x}}}}}}{f_{\left|{{h_{FP}}}\right|}}\\!\left(\\!{\sqrt{\frac{{{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P\frac{{1-x{\kappa^{2}}}}{x}}}}}\\!\right).$
(B-2)

Moreover, with the help of [61, eq. (62)], we obtain the PDF of
$\left|{{h_{FP}}}\right|$ as

${f_{\left|{{h_{FP}}}\right|}}(x)=\int_{0}^{{A_{o}}}{\frac{1}{y}}{f_{\left|{{h_{F}}}\right|}}\left({\frac{x}{y}}\right){f_{\left|{{h_{P}}}\right|}}(y){\rm{d}}y.$
(B-3)

Substituting (B-3) into (B-2), we derive the integral representation for the
PDF of $\left|{{h_{FP}}}\right|$ as (1).

###  B-B Proof of CDF

Following similar procedures as Appendix B-A (i.e., the proof of Lemma 1), we
can derive the CDF of $V$ and $W$ as
${F_{V}}\left(v\right)=1-{F_{\left|{{h_{FP}}}\right|}}\left({\sqrt{\frac{1}{v}}}\right)$
and
${F_{W}}\left(w\right)={F_{V}}\left({\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}w}\right)$.
Making the change of variable ${\gamma_{N}}=\frac{1}{{W+{\kappa^{2}}}}$, we
can derive the CDF of ${\gamma_{N}}$ as
${F_{{\gamma_{N}}}}\left(x\right)={F_{\left|{{h_{FP}}}\right|}}\left({\sqrt{\frac{{x{N_{o}}}}{{{{\left|{{h_{L}}}\right|}^{2}}P\left({1-x{\kappa^{2}}}\right)}}}}\right)$.

With the help of [61, eq. (65)], the CDF of $\left|{{h_{FP}}}\right|$ can be
expressed as

${F_{\left|{{h_{FP}}}\right|}}(x)=\int_{0}^{{A_{o}}}{{F_{\left|{{h_{F}}}\right|}}}\left({\frac{x}{y}}\right){f_{\left|{{h_{P}}}\right|}}(y){{{\rm
d}}}y.$ (B-4)

Substituting (B-4) into ${F_{{\gamma_{N}}}}\left(x\right)$, the integral
representation for the CDF of $\left|{{h_{FP}}}\right|$ as (39) can be derived
to complete the proof.

##  Appendix C Proof of Theorem 1

We first derive the PDF and CDF of
$\left|{{h_{FP}}}\right|\buildrel\Delta\over{=}\left|{{h_{F}}}\right|\times\left|{{h_{P}}}\right|$.

###  C-A PDF of $\left|{{h_{FP}}}\right|$

Using the definition of multivariate Fox’s $H$-function [48, eq. (A-1)], we
can express [41, eq. 11] as

$\displaystyle{f_{\left|h_{F}\right|}}(y){{=}}\\!\\!\\!\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1.2}},\ldots,{j_{L,2}}=0}^{{N_{T}}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!\Gamma\\!({m_{\iota,\ell}})\Gamma\\!({j_{\iota,\ell}}+1)}}}}}$
$\displaystyle\times\\!\frac{1}{y}\prod\limits_{\iota=1}^{L}\frac{1}{{2\pi
i}}\int_{{{\cal
L}_{\ell}}}\frac{{\prod\limits_{\ell=1}^{2}{\Gamma\\!\left({1+{j_{\iota,\ell}}+\frac{1}{2}{\varsigma_{\iota}}}\right)}}}{{\Gamma\\!\left({-\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}}\right)\Gamma^{-1}\\!\left({-{\varsigma_{\iota}}}\right)}}{{\left(\\!{{y^{-1}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\\!\\!\right)}^{{\varsigma_{\iota}}}}{\rm
d}{\varsigma_{\iota}},$ (C-1)

where the integration path of ${\cal L_{\ell}}$ $(\ell=1,2,\ldots,L)$ goes
from $a_{\ell}-\infty j$ to $a_{\ell}+\infty j$ and $a_{\ell}\in\mathbb{R}$.
Substituting (C-A) and (16) into (B-3), with the help of Fubini’s theorem and
[29, eq. (8.331.1)], we obtain

$\displaystyle{f_{\left|{{h_{FP}}}\right|}}(x)\\!=\\!\\!\\!\\!\\!\sum\limits_{{j_{1\\!,\\!1}},\ldots,{j_{L\\!,\\!1}}=0}^{N_{T}}\\!\\!\\!\\!\ldots\\!\\!\\!\\!\sum\limits_{{j_{1\\!,\\!N}},\ldots,{j_{L\\!,\\!N}}=0}^{N_{T}}\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!\Gamma\\!({m_{\iota,\ell}})\Gamma\\!({j_{\iota,\ell}}\\!+\\!1)}}}}$
$\displaystyle\times\\!\\!\frac{{{\gamma^{2}}}}{{xA_{o}^{{\gamma^{2}}}}}\prod\limits_{\iota=1}^{L}\frac{1}{{2\pi
i}}\int_{{{\cal
L}_{\ell}}}\frac{{I_{A}}{\prod\limits_{\ell=1}^{2}{\Gamma\\!\left(\\!{1\\!+\\!{j_{\iota,\ell}}\\!+\\!\frac{1}{2}{\varsigma_{\iota}}}\\!\right)}}}{{\Gamma\\!\left(\\!{-\\!\sum\limits_{\iota{\rm{=1}}}^{L}{{\varsigma_{\iota}}}}\\!\right)\Gamma^{-1}\\!\left(\\!{-{\varsigma_{\iota}}}\\!\right)}}{{\\!\left(\\!{{\frac{1}{x}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\\!\\!\right)}^{{\varsigma_{\iota}}}}\\!\\!{\rm
d}{\varsigma_{\iota}},$ (C-2)

where

${I_{A}}{{=\\!\\!}}\int_{0}^{{A_{\rm{0}}}}\\!\\!{{y^{{\varsigma_{\iota}}+{\gamma^{2}}-1}}{\rm{d}}y}{\rm{=}}\frac{{\rm{1}}}{{{\varsigma_{\iota}}\\!+\\!{\gamma^{2}}}}{A_{\rm{0}}}^{{\varsigma_{\iota}}+{\gamma^{2}}}=\frac{{\Gamma\\!\left({{\varsigma_{\iota}}\\!+\\!{\gamma^{2}}}\right)}}{{\Gamma\\!\left({{\varsigma_{\iota}}\\!+\\!{\gamma^{2}}\\!+\\!1}\right)}}{A_{\rm{0}}}^{{\varsigma_{\iota}}\\!+\\!{\gamma^{2}}}.$
(C-3)

Substituting ${I_{A}}$ into (C-A) and using [48, eq. (A-1)], we can obtain
${f_{\left|{{h_{FP}}}\right|}}(x)$.

###  C-B CDF of $\left|{{h_{FP}}}\right|$

Substituting [41, eq. 12] and (16) into (B-4), we can express the CDF of
${\left|{{h_{FP}}}\right|}$ as

$\displaystyle{F_{\left|{{h_{FP}}}\right|}}(x)\\!=\\!\\!\\!\\!\\!\sum\limits_{{j_{1\\!,\\!1}},\ldots,{j_{L\\!,\\!1}}=0}^{N_{T}}\\!\\!\\!\\!\ldots\\!\\!\\!\\!\sum\limits_{{j_{1\\!,\\!N}},\ldots,{j_{L\\!,\\!N}}=0}^{N_{T}}\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!\Gamma\\!({m_{\iota,\ell}})\Gamma\\!({j_{\iota,\ell}}\\!+\\!1)}}}}$
$\displaystyle\times\\!\\!\frac{{{\gamma^{2}}}}{{A_{o}^{{\gamma^{2}}}}}\prod\limits_{\iota=1}^{L}\frac{1}{{2\pi
i}}\int_{{{\cal
L}_{\iota}}}\\!\\!\frac{{I_{A}}{\prod\limits_{\ell=1}^{2}{\Gamma\\!\left({1\\!+\\!{j_{\iota,\ell}}\\!+\\!\frac{1}{2}{\varsigma_{\iota}}}\right)}}}{{\Gamma\\!\left(\\!{1\\!-\\!\sum\limits_{\iota{\rm{=1}}}^{L}{{\varsigma_{\iota}}}}\\!\right)\Gamma^{-1}\\!\left(\\!{-{\varsigma_{\iota}}}\\!\right)}}\\!{{\left(\\!\\!{\frac{{\rm{1}}}{x}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\\!\\!\right)}^{{\varsigma_{\iota}}}}{\rm
d}{\varsigma_{\iota}},$ (C-4)

where ${I_{A}}$ has been solved in (C-3). Substituting ${I_{A}}$ into (C-B),
the CDF of ${\left|{{h_{FP}}}\right|}$ can be expressed in terms of
multivariate Fox’s $H$-function. Therefore, with the help of (1), (39), (C-A),
and ${F_{\left|{{h_{FP}}}\right|}}(x)$, the PDF and CDF of ${{\gamma_{N}}}$
can be derived as (1) and (1), respectively.

##  Appendix D Proof of Proposition 2

To derive the approximate PDF of $\gamma_{N}$, we need first to get the PDF of
$\left\|{{{\bf{h}}_{F}}}\right\|$, where
$\left\|{{{\bf{h}}_{F}}}\right\|=\sum\limits_{\ell=1}^{L}{\left|{g_{2}^{\left(\ell\right)}}\right|\left|{g_{1}^{\left(\ell\right)}}\right|}$.
With the help of a 5FM based method proposed in [62] and the definition of
Meijer’s $G$ function [29, eq. (9.30)], we can obtain the approximate PDF of
$\left\|{{{\bf{h}}_{F}}}\right\|$ as

${F_{\left|{{h_{F}}}\right|}}=\frac{{{a_{1}}}}{{2\pi i}}\int_{\cal
L}{\frac{{x\Gamma\left({{a_{4}}-s}\right)\Gamma\left({{a_{5}}-s}\right)\Gamma\left({s+1}\right)}}{{\Gamma\left({{a_{3}}-s}\right)\Gamma\left({s+2}\right)}}}{\left({\frac{x}{{{a_{2}}}}}\right)^{s}}{\rm
d}s.$ (D-1)

Substituting (D-1) and (16) into (39), we have

$\displaystyle{F_{{\gamma_{N}}}}\left(x\right)$
$\displaystyle={\gamma^{2}}A_{0}^{-{\gamma^{2}}}\frac{{{a_{1}}}}{{2\pi
i}}\int_{\cal
L}{\frac{{\Gamma\left({{a_{4}}-s}\right)\Gamma\left({{a_{5}}-s}\right)\Gamma\left({s+1}\right)}}{{\Gamma\left({{a_{3}}-s}\right)\Gamma\left({s+2}\right)}}}$
$\displaystyle\times\Upsilon{\left({\frac{1}{{{a_{2}}}}\Upsilon}\right)^{s}}{I_{D}}{\rm
d}s,$ (D-2)

where $I_{D}$ can be solved with the help of [29, eq. (8.331.1)] as

${I_{D}}{\rm{=}}\\!\int_{0}^{{A_{o}}}\\!{{y^{-s-2+{\gamma^{2}}}}{\rm{d}}y}=\\!\frac{{{A_{o}}^{-1+{\gamma^{2}}-s}}}{{-1\\!+\\!{\gamma^{2}}\\!-\\!s}}=\frac{{{A_{o}}^{-1+{\gamma^{2}}-s}\Gamma\\!\left({{\gamma^{2}}-s-1}\right)}}{{\Gamma\\!\left({{\gamma^{2}}-s}\right)}}.$
(D-3)

Substituting $I_{D}$ into (D) and using [29, eq. (9.30)], we can obtain (51)
to complete the proof.

##  Appendix E Proof of Proposition 3

With the help of Jensen’s inequality [56], the upper-bound of the ergodic
capacity can be expressed as
$C_{N}\leq{\log_{2}}\left({1+{\mathbb{E}}\left[{{\gamma_{N}}}\right]}\right)$,
where

${\mathbb{E}}\left[{{\gamma_{N}}}\right]\\!{{=}}{\mathbb{E}}\left[\\!{\frac{{{{\left|{{h_{FP}}}\right|}^{2}}{{\left|{{h_{L}}}\right|}^{2}}P}}{{{\kappa^{2}}{{\left|{{h_{FP}}}\right|}^{2}}{{\left|{{h_{L}}}\right|}^{2}}P+{N_{o}}}}}\\!\right]\\!{\rm{=}}\frac{{P{{\left|{{h_{L}}}\right|}^{2}}{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]}}{{P{\kappa^{2}}{{\left|{{h_{L}}}\right|}^{2}}{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{2}}}\right]+{N_{o}}}}.$
(E-1)

The $s_{\rm th}$ moment of ${\left|{{h_{FP}}}\right|}$ is defined as
${\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{s}}}\right]{{=}}\int_{0}^{\infty}{{x^{s}}{f_{{h_{FP}}}}(x){\rm
d}x}$. With the help of (C-A),
${\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{s}}}\right]$ can be expressed
as (E), shown at the bottom of the next page, where
${I_{C}}=\int_{0}^{\infty}{{x^{s-\sum\limits_{\iota=1}^{L}{{\varsigma_{\iota}}}-1}}{\rm
d}x}$.

$\displaystyle{\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{s}}}\right]{{=}}$
$\displaystyle\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\ldots\sum\limits_{{j_{1,N}},\ldots,{j_{L,N}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}}}{{{j_{\iota,\ell}}!}}\frac{{{m_{\iota,\ell}}^{{m_{\iota,\ell}}}}}{{\Gamma({m_{\iota,\ell}})}}\frac{{{\gamma^{2}}}}{{\Gamma({j_{\iota,\ell}}+1)}}}}}$
$\displaystyle\times\prod\limits_{\iota=1}^{L}{\frac{1}{{2\pi i}}\int_{{{\cal
L}_{\ell}}}{\frac{{\prod\limits_{\ell=1}^{2}{\Gamma\left({1+{j_{\iota,\ell}}+\frac{1}{2}{\varsigma_{\iota}}}\right)}\Gamma\left({-{\varsigma_{\iota}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}}\right)}}{{\Gamma\left({-\sum\limits_{\iota{\rm{=1}}}^{L}{{\varsigma_{\iota}}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}+1}\right)}}{{\left({{A_{\rm{0}}}\prod\limits_{\ell=1}^{2}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\right)}^{{\varsigma_{\iota}}}}{I_{C}}}{\rm
d}{\varsigma_{\iota}}}$ (E-2)

Let ${\cal L}\left\\{{p\left(t\right)}\right\\}=P\left(x\right)$. With the
help of the property of Laplace transform, it follows that ${\cal
L}\left\\{{\int_{0}^{t}{p\left(z\right)dz}}\right\\}=\frac{{P\left(x\right)}}{x}$.
According to the final value theorem, we have

$\mathop{\lim}\limits_{t\to{N_{T}}}\left({\int_{0}^{t}{p\left(z\right)dz}}\right)=e\frac{{P\left(e\right)}}{e}=P\left(e\right),$
(E-3)

where $e$ is a number close to zero (e.g., $e=10^{-6}$). Thus, ${I_{C}}$ can
be solved as
${I_{C}}={\left({\frac{1}{e}}\right)^{s-\sum\limits_{\iota=1}^{L}{{\varsigma_{\iota}}}}}\Gamma\left({s-\sum\limits_{\iota=1}^{L}{{\varsigma_{\iota}}}}\right)$.
Substituting ${I_{C}}$ into (E) and using [48, eq. (A-1)], we can rewrite
${\mathbb{E}}\left[{{{\left|{{h_{FP}}}\right|}^{s}}}\right]$ in terms of
multivariate Fox’s $H$-function. With the help of (E-1), letting $s=2$, we
obtain (58) to complete the proof.

##  Appendix F Proof of Proposition 4

For the case of ideal RF chains, the ergodic capacity can be expressed as
${C_{I}}=\int_{0}^{\infty}{{{\log}_{2}}\left({1+x}\right)}{f_{\left|{{\gamma_{I}}}\right|}}\left(x\right){\rm
d}x$. Substituting (2) into ${C_{I}}$ and using [63, eq. (01.04.07.0002.01)],
we obtain

$\displaystyle C_{I}{{=}}$
$\displaystyle\frac{1}{2}\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}{\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{N_{T}}{\prod\limits_{\iota=1}^{L}\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}{\gamma^{2}}}}{{{j_{\iota,\ell}}!\Gamma({m_{\iota,\ell}})\Gamma({j_{\iota,\ell}}+1)}}}}}$
$\displaystyle\times\prod\limits_{\iota=1}^{L}\frac{1}{{2\pi i}}\int_{{{\cal
L}_{\ell}}}\frac{{\prod\limits_{\ell=1}^{2}{\Gamma\left({1+{j_{\iota,\ell}}+\frac{1}{2}{\varsigma_{\iota}}}\right)}\Gamma\left({-{\varsigma_{\iota}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}}\right)}}{{\Gamma\left({-\sum\limits_{\iota{\rm{=1}}}^{L}{{\varsigma_{\iota}}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}+1}\right)}}$
$\displaystyle\times{{\left({{A_{{0}}}\sqrt{\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}}\prod\limits_{\ell=1}^{{2}}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\right)}^{{\varsigma_{\iota}}}}{I_{D_{1}}}{\rm
d}{\varsigma_{\iota}},$ (F-1)

where

$\displaystyle{I_{D_{1}}}=$
$\displaystyle\int_{0}^{\infty}{{{\log}_{2}}\left({1+x}\right)}{x^{-1-\frac{1}{2}\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}}}dx$
$\displaystyle{{=}}\frac{1}{{2\pi i}}\int_{{\cal
L}_{L+1}}{\frac{{\Gamma\left({{\varsigma_{L+1}}+1}\right){\Gamma^{2}}\left({-{\varsigma_{L+1}}}\right)}}{{\Gamma\left({1-{\varsigma_{L+1}}}\right)}}{I_{{D_{{2}}}}}}{\rm
d}{\varsigma_{L+1}},$ (F-2)

and $I_{D_{2}}$ can be solved by following similar procedures as (E-3) as

$\displaystyle{I_{D_{2}}}\\!=$
$\displaystyle\int_{0}^{\infty}{{x^{-1-\frac{1}{2}\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}-{\varsigma_{L+1}}}}}{\rm
d}x$
$\displaystyle\times={\left(\\!{\frac{1}{e}}\\!\right)^{-\frac{1}{2}\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}-{\varsigma_{L+1}}}}\Gamma\\!\left(\\!{-\frac{1}{2}\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}\\!-{\varsigma_{L+1}}}\right).$
(F-3)

Substituting ${I_{D_{2}}}$ and ${I_{D_{1}}}$ into (F), we obtain

$\displaystyle{C_{I}}{{=}}$
$\displaystyle\sum\limits_{{j_{1,1}},\ldots,{j_{L,1}}=0}^{N_{T}}\sum\limits_{{j_{1,2}},\ldots,{j_{L,2}}=0}^{{N_{T}}}{\prod\limits_{\iota=1}^{L}{\prod\limits_{\ell=1}^{2}{\frac{{{K_{\iota,\ell}}^{{j_{\iota,\ell}}}{d_{\iota,\ell}}_{{j_{\iota,\ell}}}{m_{\iota,\ell}}^{{m_{\iota,\ell}}}{\gamma^{2}}}}{{2{j_{\iota,\ell}}!\Gamma({m_{\iota,\ell}})\Gamma({j_{\iota,\ell}}+1)}}}}}$
$\displaystyle\times\prod\limits_{\iota=1}^{L}\frac{1}{{2\pi i}}\int_{{{\cal
L}_{\ell}}}\frac{{\prod\limits_{\ell=1}^{2}{\Gamma\left({1+{j_{\iota,\ell}}+\frac{1}{2}{\varsigma_{\iota}}}\right)}\Gamma\left({-{\varsigma_{\iota}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}}\right)}}{{\Gamma\left({-\sum\limits_{\iota{{=1}}}^{L}{{\varsigma_{\iota}}}}\right)\Gamma\left({{\varsigma_{\iota}}+{\gamma^{2}}+1}\right)}}$
$\displaystyle\times{{\left({{A_{\rm{0}}}\sqrt{e\frac{{{{\left|{{h_{L}}}\right|}^{2}}P}}{{{N_{o}}}}}\prod\limits_{\ell=1}^{\rm{2}}{\left({\sqrt{2}{\sigma_{\iota,\ell}}}\right)}}\right)}^{{\varsigma_{\iota}}}}\frac{1}{{2\pi
i}}$
$\displaystyle\times\int_{{L_{L+1}}}{\frac{{{\Gamma^{2}}\left({-{\varsigma_{L+1}}}\right)\Gamma\left({-\frac{1}{2}\sum\limits_{\iota=1}^{L}{{\varsigma_{\iota}}}-{\varsigma_{L+1}}}\right)}}{{{\Gamma^{-1}}\left({{\varsigma_{L+1}}+1}\right)\Gamma\left({1-{\varsigma_{L+1}}}\right){e^{{\varsigma_{L+1}}}}}}}{\rm{d}}{\varsigma_{L+1}}{\rm{d}}{\varsigma_{\iota}}.$
(F-4)

With the help of [48, eq. (A-1)], we can rewrite (F) as (4) that completes the
proof.

## References

  * [1] J. Zhang, E. Björnson, M. Matthaiou, D. W. K. Ng, H. Yang, and D. J. Love, “Prospective multiple antenna technologies for beyond 5G,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1637–1660, Aug. 2020.
  * [2] A.-A. A. Boulogeorgos, A. Alexiou, and T. Merkle, “Terahertz technologies to deliver optical network quality of experience in wireless systems beyond 5G,” _IEEE Commun. Mag._ , vol. 56, no. 6, pp. 144–151, Jun. 2018.
  * [3] C. Chaccour, M. N. Soorki, W. Saad, M. Bennis, and P. Popovski, “Risk-based optimization of virtual reality over terahertz reconfigurable intelligent surfaces,” in _Proc. IEEE Int. Conf. Commun._ , Jun. 2020, pp. 1–6.
  * [4] A. Faisal, H. Sarieddeen, H. Dahrouj, T. Y. Al-Naffouri, and M.-S. Alouini, “Ultramassive MIMO systems at terahertz bands: Prospects and challenges,” _IEEE Trans. Veh. Technol._ , vol. 15, no. 4, pp. 33–42, Apr. 2020.
  * [5] C. Chaccour, M. N. Soorki, W. Saad, M. Bennis, P. Popovski, and M. Debbah, “Seven defining features of terahertz (THz) wireless systems: A fellowship of communication and sensing,” _arXiv preprint arXiv:2102.07668_ , 2021.
  * [6] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, Aug. 2019.
  * [7] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [8] Y. Liu, E. Liu, and R. Wang, “Energy efficiency analysis of intelligent reflecting surface system with hardware impairments,” in _Proc. IEEE Global Commun. Conf. (GLOBECOM)_ , pp. 1–6, Dec. 2020.
  * [9] J. V. Alegria and F. Rusek, “Achievable rate with correlated hardware impairments in large intelligent surfaces,” in _Proc. IEEE Int. Workshop Comput. Adv. Multi-Sensor Adaptive Process. (CAMSAP)_ , pp. 559–563, Dec. 2019.
  * [10] A.-A. A. Boulogeorgos and A. Alexiou, “How much do hardware imperfections affect the performance of reconfigurable intelligent surface-assisted systems?” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 1185–1195, 2020.
  * [11] Z. Xing, R. Wang, J. Wu, and E. Liu, “Achievable rate analysis and phase shift optimization on intelligent reflecting surface with hardware impairments,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 5514–5530, Sept. 2021\.
  * [12] K. Tekbıyık, G. Karabulut Kurt, A. Rıza Ekti, A. Görçin, and H. Yanikomeroglu, “Reconfigurable intelligent surfaces empowered THz communication in LEO satellite networks,” _arXiv preprint arXiv:2007.04281_ , 2021.
  * [13] S. P. Maruthi, T. Panigrahi, and M. Hassan, “Improving the reliability of pulse-based terahertz communication using intelligent reflective surface,” in _Proc. IEEE Int. Conf. Commun. Workshops_ , pp. 1–6, May 2020\.
  * [14] A.-A. A. Boulogeorgos, N. Chatzidiamantis, H. G. Sandalidis, A. Alexiou, and M. Di Renzo, “Cascaded composite turbulence and misalignment: Statistical characterization and applications to reconfigurable intelligent surface-empowered wireless systems,” _IEEE Trans. Veh. Technol._ , to appear, 2022.
  * [15] X. Ma, Z. Chen, W. Chen, Z. Li, Y. Chi, C. Han, and S. Li, “Joint channel estimation and data rate maximization for intelligent reflecting surface assisted terahertz MIMO communication systems,” _IEEE Access_ , vol. 8, pp. 99 565–99 581, 2020.
  * [16] J. Qiao and M.-S. Alouini, “Secure transmission for intelligent reflecting surface-assisted mmWave and terahertz systems,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 10, pp. 1743–1747, Oct. 2020.
  * [17] A.-A. A. Boulogeorgos and A. Alexiou, “Performance analysis of reconfigurable intelligent surface-assisted wireless systems and comparison with relaying,” _IEEE Access_ , vol. 8, pp. 94 463–94 483, 2020.
  * [18] 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, Jun. 2019.
  * [19] J. Kokkoniemi, J. Lehtomaki, and M. Juntti, “Simplified molecular absorption loss model for 275-400 Gigahertz frequency band,” in _Proc. 12th Eur. Conf. Antennas Propag. (EuCAP)_ , pp. 1–5, Apr. 2018.
  * [20] A.-A. A. Boulogeorgos, E. N. Papasotiriou, J. Kokkoniemi, J. Lehtomaeki, A. Alexiou, and M. Juntti, “Performance evaluation of THz wireless systems operating in 275-400 GHz band,” in _Proc. IEEE Veh. Technol. Conf. (VTC)_ , Jun. 2018, pp. 1–5.
  * [21] A. R. Ekti, A. Boyaci, A. Alparslan, İ. Ünal, S. Yarkan, A. Görçin, H. Arslan, and M. Uysal, “Statistical modeling of propagation channels for terahertz band,” in _Proc. IEEE Conf. Standards Commun. Netw._ , pp. 275–280, Sep. 2017.
  * [22] S. Priebe, C. Jastrow, M. Jacob, T. Kleine-Ostmann, T. Schrader, and T. Kurner, “Channel and propagation measurements at 300 GHz,” _IEEE Trans. Antennas Propag._ , vol. 59, no. 5, pp. 1688–1698, May 2011.
  * [23] K. Guan, B. Peng, D. He, J. M. Eckhardt, S. Rey, B. Ai, Z. Zhong, and T. Kürner, “Measurement, simulation, and characterization of train-to-infrastructure inside-station channel at the terahertz band,” _IEEE Trans. THz Sci. Technol._ , vol. 9, no. 3, pp. 291–306, Mar. 2019.
  * [24] J. M. Romero-Jerez, F. J. Lopez-Martinez, J. F. Paris, and A. J. Goldsmith, “The fluctuating two-ray fading model: Statistical characterization and performance analysis,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 7, pp. 4420–4432, Jul. 2017.
  * [25] E. N. Papasotiriou, A.-A. A. Boulogeorgos, K. Haneda, M. F. de Guzman, and A. Alexiou, “An experimentally validated fading model for THz wireless systems,” _Sci. Rep._ , vol. 11, no. 1, pp. 1–14, Jan. 2021.
  * [26] S. Atapattu, R. Fan, P. Dharmawansa, G. Wang, J. Evans, and T. A. Tsiftsis, “Reconfigurable intelligent surface assisted two-way communications: Performance analysis and optimization,” _IEEE Trans. Commun._ , vol. 68, no. 10, pp. 6552–6567, Oct. 2020.
  * [27] Z. Peng, Z. Zhang, C. Pan, L. Li, and A. L. Swindlehurst, “Multiuser full-duplex two-way communications via intelligent reflecting surface,” _IEEE Trans. Signal Process._ , vol. 69, pp. 837–851, Jan. 2021.
  * [28] J. Zhang, W. Zeng, X. Li, Q. Sun, and K. P. Peppas, “New results on the fluctuating two-ray model with arbitrary fading parameters and its applications,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 3, pp. 2766–2770, Mar. 2017.
  * [29] I. S. Gradshteyn and I. M. Ryzhik, _Table of Integrals, Series, and Products_ , 7th ed. Academic Press, 2007\.
  * [30] A. Papoulis and S. U. Pillai, _Probability, random variables, and stochastic processes_. New York, NY, USA: McGraw-Hill, 2002.
  * [31] Matlab, “The matlab help center,” https://www.mathworks.com/help/stats/.
  * [32] A. A. Farid and S. Hranilovic, “Outage capacity optimization for free-space optical links with pointing errors,” _IEEE/OSA J. Lightw. Technol._ , vol. 25, no. 7, pp. 1702–1710, Jul. 2007.
  * [33] C. Chaccour, M. N. Soorki, W. Saad, M. Bennis, and P. Popovski, “Can terahertz provide high-rate reliable low latency communications for wireless VR?” _IEEE Internet Things J._ , to appear, 2022.
  * [34] Y. L. Babikov, I. Gordon, S. Mikhailenko, L. Rothman, and S. Tashkun, “HITRAN on the web–a new tool for HITRAN spectroscopic data manipulation,” in _Proc. 12th Intl. HITRAN Conf._ , 2012, pp. 29–31.
  * [35] E. Björnson, M. Matthaiou, and M. Debbah, “A new look at dual-hop relaying: Performance limits with hardware impairments,” _IEEE Trans. Commun._ , vol. 61, no. 11, pp. 4512–4525, Nov. 2013.
  * [36] T. L. Jensen and T. Larsen, “Robust computation of error vector magnitude for wireless standards,” _IEEE Trans. Commun._ , vol. 61, no. 2, pp. 648–657, Feb. 2013.
  * [37] T. K. Lo, “Maximum ratio transmission,” _IEEE Trans. Commun._ , vol. 47, no. 10, pp. 1458–1461, Oct. 1999.
  * [38] Y. Tokgoz and B. D. Rao, “Performance analysis of maximum ratio transmission based multi-cellular MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 5, no. 1, pp. 83–89, Jan. 2006.
  * [39] J. Zhang, H. Du, P. Zhang, J. Cheng, and L. Yang, “Performance analysis of 5G mobile relay systems for high-speed trains,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 12, pp. 2760–2772, Dec. 2020.
  * [40] Q. Wu and R. Zhang, “Beamforming optimization for wireless network aided by intelligent reflecting surface with discrete phase shifts,” _IEEE Trans. Commun._ , vol. 68, no. 3, pp. 1838–1851, Mar. 2020.
  * [41] H. Du, J. Zhang, J. Cheng, Z. Lu, and B. Ai, “Millimeter wave communications with reconfigurable intelligent surface: Performance analysis and optimization,” _IEEE Trans. Commun._ , vol. 69, no. 4, pp. 2752–2768, Apr. 2021.
  * [42] K. Feng, Q. Wang, X. Li, and C.-K. Wen, “Deep reinforcement learning based intelligent reflecting surface optimization for MISO communication systems,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 5, pp. 745–749, May 2020.
  * [43] S. Cheng, H. Lu, X. Lei, and Y. Shi, “A quarter century of particle swarm optimization,” _Complex Intell. Syst._ , vol. 4, no. 3, pp. 227–239, Mar. 2018.
  * [44] M. Clerc and J. Kennedy, “The particle swarm - explosion, stability, and convergence in a multidimensional complex space,” _IEEE Trans. Evol. Comput._ , vol. 6, no. 1, pp. 58–73, Jan. 2002.
  * [45] H. Yoshida, K. Kawata, Y. Fukuyama, S. Takayama, and Y. Nakanishi, “A particle swarm optimization for reactive power and voltage control considering voltage security assessment,” _IEEE Trans. Power Syst._ , vol. 15, no. 4, pp. 1232–1239, Apr. 2000.
  * [46] M. Juneja and S. K. Nagar, “Particle swarm optimization algorithm and its parameters: A review,” in _Proc. International Conf. Control Comput. Commun. Mater. (ICCCCM)_ , pp. 1–5, 2016.
  * [47] W. Kan and S. Jihong, “The convergence basis of particle swarm optimization,” in _Proc. Int. Conf. Industrial Control Elect. Eng._ , 2012, pp. 63–66.
  * [48] A. M. Mathai, R. K. Saxena, and H. J. Haubold, _The $H$-function: Theory and Applications_. Springer Science & Business Media, 2009.
  * [49] J. D. V. Sánchez, D. P. Moya Osorio, E. E. Benitez Olivo, H. Alves, M. C. Paredes Paredes, and L. Urquiza Aguiar, “On the statistics of the ratio of nonconstrained arbitrary $\alpha$-$\mu$ random variables: A general framework and applications,” _Trans. Emerging Tel. Tech._ , p. e3832, 2019\.
  * [50] H. R. Alhennawi, M. M. El Ayadi, M. H. Ismail, and H.-A. M. Mourad, “Closed-form exact and asymptotic expressions for the symbol error rate and capacity of the $H$-function fading channel,” _IEEE Trans. Veh. Technol._ , vol. 65, no. 4, pp. 1957–1974, Apr. 2015.
  * [51] D. Selimis, K. P. Peppas, G. C. Alexandropoulos, and F. I. Lazarakis, “On the performance analysis of RIS-empowered communications over Nakagami-$m$ fading,” _IEEE Commun. Lett._ , vol. 25, no. 7, pp. 2191–2195, Jul. 2021.
  * [52] P.-S. M. de Laplace, _Théorie analytique des probabilités_ , 1995\.
  * [53] L. Zheng and D. N. C. Tse, “Diversity and multiplexing: A fundamental tradeoff in multiple-antenna channels,” _IEEE Trans. Inf. Theory_ , vol. 49, no. 5, pp. 1073–1096, May 2003.
  * [54] A. Lozano, A. M. Tulino, and S. Verdú, “High-snr power offset in multiantenna communication,” _IEEE Trans. Inf. Theory_ , vol. 51, no. 12, pp. 4134–4151, Dec. 2005.
  * [55] M.-S. Alouini and A. J. Goldsmith, “Capacity of Rayleigh fading channels under different adaptive transmission and diversity-combining techniques,” _IEEE Trans. Veh. Technol._ , vol. 48, no. 4, pp. 1165–1181, Apr. 1999.
  * [56] S. G. Krantz, _Handbook of complex variables_. Springer Science & Business Media, 2012.
  * [57] A.-A. A. Boulogeorgos and A. Alexiou, “Coverage analysis of reconfigurable intelligent surface assisted THz wireless systems,” _IEEE Open J. Veh. Tech._ , vol. 2, pp. 94–110, Jan. 2021.
  * [58] S. Koenig, D. Lopez-Diaz, J. Antes, F. Boes, R. Henneberger, A. Leuther, A. Tessmann, R. Schmogrow, D. Hillerkuss, R. Palmer _et al._ , “Wireless sub-THz communication system with high data rate,” _Nat. Photonics_ , vol. 7, no. 12, pp. 977–981, Dec. 2013.
  * [59] T. Nagatsuma, G. Ducournau, and C. C. Renaud, “Advances in terahertz communications accelerated by photonics,” _Nat. Photonics_ , vol. 10, no. 6, pp. 371–379, Jun. 2016.
  * [60] C. Han, A. O. Bicen, and I. F. Akyildiz, “Multi-wideband waveform design for distance-adaptive wireless communications in the terahertz band,” _IEEE Trans. Signal Process._ , vol. 64, no. 4, pp. 910–922, Feb. 2015.
  * [61] A.-A. A. Boulogeorgos, E. N. Papasotiriou, and A. Alexiou, “Analytical performance assessment of THz wireless systems,” _IEEE Access_ , vol. 7, pp. 11 436–11 453, Feb. 2019.
  * [62] F. El Bouanani and D. B. da Costa, “Accurate closed-form approximations for the sum of correlated Weibull random variables,” _IEEE Wireless Commun. Lett._ , vol. 7, no. 4, pp. 498–501, Apr. 2018.
  * [63] Wolfram, “The wolfram functions site,” http://functions.wolfram.com.