Holographic RIS Empowered THz Communications with Hardware Imperfections under Adverse Weather Conditions

Alexandros-Apostolos A. Boulogeorgos, Senior Member, IEEE, Stylianos Trevlakis, Member, IEEE, and Theodoros A. Tsiftsis, Senior Member, IEEE

A.-A. A. Boulogeorgos is with the Department of Electrical and Computer Engineering, University of Western Macedonia, 50100 Kozani, Greece (e-mail:

S. Trevlakis is with the Research & Development Department, InnoCube P.C., Thessaloniki 55535, Greece. (e-mail: trevlakis@innocube.org).

Theodoros A. Tsiftsis is with the Department of Informatics & Telecommunications, University of Thessaly, Lamia 35100, Greece (e-mail:tsiftsis@uth.gr).

This paper focuses on providing a theoretical framework for the evaluation of the performance of holographic reconfigurable intelligent surface (HRIS) empowered terahertz (THz) wireless systems under fog conditions. In more detail, we present a comprehensive methodology for evaluating the geometric losses of the end-to-end channel. Moreover, the stochastic nature of the end-to-end channel is characterized by novel closed-form probability density and cumulative distribution functions. Building upon them, the outage probability and the throughput of the system are extracted in closed form. These formulas account for the impact of transceiver hardware imperfections and are expected to become useful tools for the design of HRIS empowered THz wireless systems due to its remarkable engineering insights.

Fog, hardware imperfections, holographic reconfigurable intelligent surface (HRIS),  outage probability, stochastic characterization, throughput.

§ INTRODUCTION

§.§ State of the Art & Motivation

Bringing the fiber quality of experience into the wireless world through the development of wireless fiber extenders capable of supporting data rates in the Tb/s region is a fundamental promise of terahertz (THz) technologies <cit.>. THz wireless fiber extenders aim to establish low-cost medium to long range high-bandwidth links for mid- and back-hauling without demanding the time-consuming processes of optical fiber installation <cit.>. As a consequence, from a business point of view, flexible almost-instance high-data-rate connectivity is expected to bring important economic benefits, while allowing internet democratization <cit.>.

Nevertheless, THz wireless communication systems are susceptible to a number of drawbacks. First of all, they suffer from high pathloss due to the high-transmission frequencies as well as molecular resonances that lead to molecular adsorptions <cit.>. Motivated by this, a great amount of research effort was put on modeling the deterministic path loss of THz wireless links <cit.>. In more detail, in <cit.>, the authors employed the radiation theory in order to extract the channel characteristics of nano-scale THz wireless links. In <cit.>, an electromagnetic analysis was conducted for the propagation characterization of in-body THz links, whilst in <cit.>, the authors presented a path loss model for in-body nano-scale THz wireless systems. A simplified path loss model for THz systems operating in the $275-400\,\rm{GHz}$ region was presented in <cit.>.

Another important drawback of THz wireless systems is the signal fluctuation due to atmospheric conditions. The source of this phenomenon lies in the variation of the reflection index, as a result of  inhomogeneities in atmospheric pressure as well as temperature and water molecular density across the propagation path <cit.>.  Inspired by this, the authors of <cit.> used a weather emulating chamber in order to quantify in laboratory conditions the impact of rain in THz and infrared links. In <cit.>, the authors quantified the joint impact of rain and misalignment in THz wireless systems in terms of outage probability. The detrimental impact of turbulence in THz wireless systems was discussed in <cit.>. The joint impact of snow and misalignment in THz wireless systems was evaluated in terms of bit error rate and channel capacity in <cit.>. The authors of <cit.> quantified the impact of fog in THz wireless systems in terms of bit error rate in laboratory environment. Finally, in <cit.>, the authors experimentally demonstrated the distortion to the THz signal that was caused due to fog.

To counterbalance the significant path loss as well as the signal fluctuation due to the atmospheric conditions and establish long-range connections, high-directional antennas, such as Cassegrain antennas, are employed at both the transmitter (TX) and the receiver (RX) ends, establishing high-directional links <cit.>. However, these links are sensitive to blockage, due to the high penetration loss of the THz transmission. Therefore, the need to create alternative paths between the TX and RX arises.  On another front, with the explosive growth of data and bandwidth requirements in beyond fifth generation (5G) and sixth generation (6G) wireless communication networks, there is a unprecedented demand for improving key performance indicators (KPIs) such as energy efficiency, spectral efficiency, coverage and serviceability at very high levels <cit.>.

To this end, reconfigurable intelligent surfaces (RIS) are emerging as the revolutionary technology based on programmable metamaterial-based surfaces and seem to be the fundamental component of 6G communication infrastructures in improving the above mentioned KPIs <cit.>. In particular, conventional RIS plates comprises an array of nearly passive reflecting meta-elements, which are able to control impinging signals through a field programmable gate array (FPGA) microcontroller. As a result, RIS adapts both the amplitude and phase of incident signals to the specific wireless environment conditions in order to either increase the coverage or cancel interference, and highlight its effectiveness in various applications such as drone-based communications, integrated sensing and communications, full-duplex communications, etc. <cit.>.

Above and beyond all other considerations, due to the adaptability of RIS in creating alternative wireless paths in highly penetrated path loss environments i.e., THz communication systems, RIS appears to be a reliable technology to to combat communication blockage <cit.>. Specifically, the authors in <cit.> presented for the first time key performance results (path loss and error rate) of RIS-enabled wireless communication systems, and conducted theoretical comparisons with well-established technologies such as such as relaying, multiple-input multiple-output (MIMO)
beamforming, and backscatter communications. In <cit.>, the path loss of RIS-assisted wireless
systems was analytically investigated based on the vector generalization of Green’s theorem. In <cit.>, an improved path loss model for RIS-enabled wireless communications at mmWave bands was proposed. In particular, practical electromagnetic phenomena (e.g., gain patterns, received power, phase errors, and specular reflection loss per meta-element), were investigated. A similar approach is followed in <cit.> where an angle-dependent  path loss model for RISs was proposed based on the radiation patterns of all involved parts of the communication system. However, the proposed prototype in <cit.> are based on PIN diodes, which is impractical for sub-THz bands <cit.>. Furthermore, RIS is highly affected by the frequency range of the incident signal to the RIS meta-elements. Therefore, RIS's reflection and phase adjustment capability employing either varactors or PIN diodes is limited over RF frequency range and, specifically, between 100 MHz - 10 GHz. Therefore, for much higher frequencies at the sub-THz bands, i.e., 0.1-0.3 THz, micro-electromechanical system (MEMS), mechanical approach, liquid crystal, and microﬂuidics appear to be the most suitable tuning RIS technologies <cit.>. It is worth mentioning that the latter two tuning technologies cover a wide range of THz bandwidths compared to the two first mentioned tuning ones   <cit.>. The operation of RIS at sub-THz is also feasible with tuning technologies involving complementary metal-oxide-semiconductor (CMOS) transistors, Schottky diodes,
or high electron mobility transistors (HEMTs) <cit.>. Very recently, the authors in <cit.>, have demonstrated for the first time a THz point-to-point RIS-aided transmission based on real-time beam tracking. In this experiment, the THz signals  near  $0.34$ THz are controlled via a RIS plate consisted of GaN HEMTs

Despite the fact that RIS is potentially considered as the key enabling technology to drive the evolution of beyond 5G and 6G era, it is true that confronts critical challenges. For example, its nearly passive nature and the absence of any active components makes conventional RIS to perform channel estimation and beam tracking as one of the most challenging issues. Consequently, in real life setups, conventional RISs constraint the bandwidth used for transmission and the achievable data rate.
To address the above issues in conventional RIS, the holographic RIS  (HRIS) (a.k.a. as reconfigurable holographic surface), has been emerged recently as an alternate technology <cit.>. HRIS is considered as a limited surface area that integrates a virtually inﬁnite number of tiny metamaterial-based radiation elements in order to form a spatially continuous transceiver aperture that can achieve the holographic beamforming. HRIS is able to support channel estimation and act like continuous surface for higher frequencies up to the THz bandwidth compared to conventional RIS <cit.>.

In one of the most representative works on HRIS <cit.>, the authors  proposed the concept of HRIS  and studied its application to massive MIMO in THz frequency range. More, in this pioneer work, the authors have addressed the transmission design of nearly-passive HRISs with
spatially continuous apertures by designing the beam pattern of a RIS with discrete meta-elements and generalize their beamforming framework to a surface with closely spaced (i.e. continuous) meta-elements yielding the HRIS design. In <cit.>, the authors studied the problem of beam pattern design of RIS-assisted THz communication systems with two-dimensional ﬁnite impulse response ﬁlter design, and the above method proved to facilitate applications such as wireless positioning and over-the-air computation.
In <cit.>, the joint phase and delay wideband precoding was investigated to tackle the severe problem of array gain loss of RIS-aided THz communications due to beam split at RIS. Additionally, the authors in <cit.> have investigated the performance of HRIS non-orthogonal multiple access (NOMA) networks, whereas the same authors in <cit.> have very recently evaluated the outage of HRISs-enabled THz NOMA systems in the presence of misalignment errors and for either perfect or imperfect successive interference cancellation.

Motivated by the fact that HRIS supports the same functionalities as conventional RIS, but most importantly due to its holographic pattern is an attractive candidate for RIS-aided wireless communication systems working at the sub-THz band, in this paper we study  the performance of HRIS-empowered THz wireless systems with hardware imperfections at both ends and under foggy weather conditions.

§.§ Our Contributions

In contrast that most of the published works study the application of conventional RIS in either sub-6GHz bands or mmwave, in this paper, the communication theoretical framework of HRIS in sub-THz bands over severe weather conditions (fog) is analytically studied under the presence of transceiver hardware imperfections. The contributions of our work can be summarized as follows:

* The end-to-end geometric losses of the proposed system are calculated by considering the  following channel coefficients: free-space geometric channel gain, molecular absorption gain, and fog gain. Our analysis shows that the geometric losses are quite affected by the following parameters: the relative position of HRIS in respect to the TX, HRIS size and its orientation, TX/RX antenna gains, the signal frequency, as well as the weather conditions.

* The statistics of the end-to-end channel coefficients of HRIS-enabled THz systems is analytically evaluated in closed form. Based on the latter derivations, new closed-form expressions for the outage probability and data throughput are extracted when hardware imperfections take place at both ends of the considered system in foggy weather conditions.

* New engineering insights are extracted such as how the geometric loss is affected by various frequencies in the sub-THZ band, and what is the optimal placement of HRIS between transmitter and receiver. Additionally, our paper provides a comprehensive study on the impact of the placement and orientation of HRIS to both outage and throughput performance metrics under fixed values of frequency range, spectral efficiency, and transceiver hardware imperfection parameters.

§.§ Organization of the paper

This paper is structured as follows. In the sequel, a few basic mathematical symbols used in this paper are given. The system and signal model studied in this paper are introduced in Section II. In Section III, the geometric losses and the corresponding statistical analysis of the end-to-end channels are analytically presented, while both the outage probability and throughput are obtained in closed form. In Section IV, the analytical expressions are corroborated and cross-compared with Monte Carlo simulation results. In the same Section, some useful engineering insights are revealed and discussed in much detail. Finally, some concluding remarks are given in Section V.

§.§ Notations

In this paper, $\sqrt{x}$ returns the square root of $x$. The cosine and exponential functions are denoted by $\cos\left(\cdot\right)$ and $\exp\left(\cdot\right)$, respectively. The natural logarithm is represented as $\ln(\cdot)$, while $\log_{a}(\cdot)$ is the logarithm to base $a$. The operator $\min\left(a, b\right)$ returns the minimum between $a$ and $b$. The gamma and the upper incomplete gamma functions are respectively represented as $\Gamma(\cdot)$ <cit.> and $\Gamma\left(\cdot,\cdot\right)$ <cit.>. The Kummer confluent hypergeometric function is denoted by $\,_1F_1\left(\cdot,\cdot,\cdot\right)$ <cit.>. Moreover, $H_A\left(a,b;c;x_1, x_2\right)$ stands for the Lauricella hypergeometric function <cit.>. Finally, the Pochhammer operator is denoted by $(n)_m$.

§ SYSTEM AND SIGNAL MODEL

Typical model of an HRIS empowered THz communication system.

As demonstrated in Fig. <ref>, we consider a long-range THz wireless fiber extender that is established through a HRIS. By assuming that the transmitter (TX)-HRIS and HRIS-receiver (RX) links experience fog, the baseband equivalent received signal can be expressed as

\begin{align}
r = h_g\,h_1 \, h_2 \, (s + n_t) + n_r + n,
\label{Eq:r}
\end{align}

where $h_g$ stands for the end-to-end geometric gain channel coefficient, while $h_1$ and $h_2$ are independent random processes that model the impact of fog in the TX-HRIS and HRIS-RX links, respectively. The transmitted signal is denoted by $s$, whereas $n_t$ and $n_r$ stands for the distortion noises due to hardware imperfections at the TX and RX, respectively. The additive white Gaussian noise is modeled by the zero-mean complex Gaussian random process $n$ of variance $\sigma_n^2$.

The end-to-end geometric gain channel coefficient can be analyzed as

\begin{align}
h_g = h_{g,f}\,h_{g,m}\,h_{f},
\label{Eq:hg}
\end{align}

where $h_{g,f}$ is the free space geometric gain channel coefficient that can be expressed as <cit.>

\begin{align}
h_{g,f} = \frac{c\,\sqrt{G_t\,G_r}\,l_h\,l_v}{4\pi\,f\,d_1\,d_2}\cos\left(\psi\right),
\label{Eq:h_gf}
\end{align}

with $c$ being the speed of light, $G_t$ and $G_r$ standing for the transmission and reception antenna gains, respectively, while $l_h$ and $l_v$ representing the length and width of the HRIS. Also, $f$ is the transmission frequency, whereas $d_1$ and $d_2$ are the TX-HRIS and HRIS-RX distances, respectively. Finally, $\psi$ denotes the beam incident angle at the HRIS.

Additionally, $h_{g,m}$ is the molecular absorption gain coefficient that can be obtained as

\begin{align}
h_{g,m} = \exp\left(\frac{1}{2}\,\kappa_m\left(f, T, \phi, p\right)\,\left(d_1+d_2\right)\right),
\label{Eq:h_gm}
\end{align}

where $\kappa_m$ is the molecular absorption coefficient, $T$ is the atmospheric temperature, $\phi$ stands for the relative humidity, and $p$ represents the atmospheric pressure.

Based on <cit.>, the molecular absorption coefficient can be written as

\begin{align}
\kappa_m\left(f, T, \phi, p\right)=\sum_{k,l} \kappa^{g}_{k,l}\left(f, T, \phi, p\right),
\end{align}

where $ \kappa^{g}_{k,l}\left(f, T, \phi, p\right)$ represents the molecular absorption coefficient for the isotoplogue $k$ of the $l-$th gas and can be expressed as

\begin{align}
\kappa^{g}_{k,l}\left(f, T, \phi, p\right) = \frac{p}{p_o} \, \frac{T_o}{T} Q^{g}_{k,l}\,\sigma_{k,l}^{g}\left(f, T, \phi, p\right),
\end{align}

with $p_o$ and $T_o$ denoting the standard pressure and temperature respectively, while $Q^{g}_{k,l}$ and $\sigma_{k,l}^{g}\left(f, T, \phi, p\right)$ respectively being the molecular volumetric density and the absorption cross section for the isotopologue $k$ of the $l-$th gas. The molecular volumetric density can be evaluated as

\begin{align}
Q^{g}_{k,l} = \frac{p}{R\,T}\,q_{k,l}\,N_{A},
\end{align}

where $R$ is the gas constant, $N_A$ stands for the Avogadro constant, and $q_{k,l}$ represents the mixing ratio for the isotopologue $k$ of the $l-$th gas. Note that in the high resolution transmission (HITRAN) database, each isologue contribution is scaled based on its natural abundance in the medium. As a consequence, instead of $q_{k,l}$, the mixing ratio of the specific gas is employed.

The absorption cross section for the isotopologue $k$ of the $l-$th gas can be calculated as

\begin{align}
\sigma_{k,l}^{g}\left(f, T, \phi, p\right) = S_{k,l}\,G_{k,l}\left(f, T, \phi, p\right),
\end{align}

where $S_{k,l}$ is the absorption strength of a specific type of molecule and can be directly extracted from the HITRAN database, while $G_{k,l}\left(f, T, \phi, p\right)$ is the spectral line shape for the isotopologue $k$ of the $l-$th gas that can be evaluated as in <cit.>

\begin{align}
G_{k,l}\left(f, T, \phi, p\right) = \frac{f}{f_{k,l}} \, \frac{\tanh\left( \frac{h\,c\,f}{2\,k_b\,T}\right)}{\tanh\left( \frac{h\,c\,f_{k,l}}{2\,k_b\,T}\right)}\,F_{k,l}\left(f\right),
\end{align}

with $h$ and $k_b$ respectively being the Plank and Bolzmann constants, $f_{k,l}$ representing the resonant frequency for the isotopologue $k$ of the $l-$th gas, and $F_{k,l}\left(f\right)$ denoting the Van Vleck-Weisskopf asymmetric line shape that can be calculated as in <cit.>

\begin{align}
& F_{k,l}\left(f\right) = 100\,c\frac{\alpha_{k,l}}{\pi}\frac{f}{f_{k,l}} \nonumber \\ & \times \left( \frac{1}{\left(f-f_{k,l}\right)^2+\alpha_{k,l}^2} + \frac{1}{\left(f+f_{k,l}\right)^2+\alpha_{k,l}^2} \right),
\end{align}

where $\alpha_{k,l}$ is the Lorentz half-width that can be expressed as

\begin{align}
\alpha_{k,l} = \left( \left(1-q_{k,l}\right)\alpha_o^{a} + q_{k,l}\,a_{k,l}^{o}\right) \frac{p}{p_o} \left(\frac{T_o}{T}\right)^{t},
\end{align}

with $t$ being the temperature broadening coefficient, while $\alpha_o^{a}$ and $a_{k,l}^{o}$ respectively denoting Lorentz half-width of air and the reference value of the Lorentz half-width for the isotopologue $k$ of the $l-$th gas. Note that $t$, $\alpha_o^{a}$, and $a_{k,l}^{o}$ can be directly extracted from HITRAN.

Note that except from the aforementioned model for the molecular absorption coefficient, a number of approximations including <cit.> have been published. Although these approximations are very accurate and have been widely used, since their application is in the area of $100-500\,\rm{GHz}$, in this paper, we decided to use a more general model that can cover all the THz band.

Finally, $h_{f}$ denotes the fog gain, which based on <cit.>, can be analyzed as

\begin{align}
h_f = 10^{-\frac{1}{2}\kappa_f\,M\,\frac{d_1+d_2}{1000}},
\label{Eq:h_f}
\end{align}

where $M$ is the liquid water density in the fog and $\kappa_f$ is the fog attenuation coefficient, which can be obtained as

\begin{align}
\kappa_f = \frac{0.819\,f}{\epsilon_i\left(1+\eta^2\right)},
\end{align}

\begin{align}
\eta = \frac{2+\epsilon_i}{\epsilon_r}.
\label{Eq:eta}
\end{align}

In (<ref>), $\epsilon_r$ and $\epsilon_i$ are the real and imaginary part of the dielectric permittivity of water, which can be expressed as

\begin{align}
\epsilon_r = \frac{\epsilon_0-\epsilon_1}{1+\left(\frac{f}{f_p}\right)^2} + \frac{\epsilon_1-\epsilon_2}{1+\left(\frac{f}{f_s}\right)^2}+\epsilon_2
\end{align}

\begin{align}
\epsilon_i = \frac{f}{f_p} \frac{\epsilon_0-\epsilon_1}{1+\left(\frac{f}{f_p}\right)^2} + \frac{f}{f_s} \frac{\epsilon_1-\epsilon_2}{1+\left(\frac{f}{f_s}\right)^2},
\end{align}

respectively. The parameters $f_p$ and $f_s$ denote the principal and secondary relaxation frequencies, which can be evaluated as

\begin{align}
f_p =& 20.2\times 10^{-9} - 146\times 10^{-9}\left(\theta-1\right)
\nonumber \\ &
+ 316\times 10^{-9}\left(\theta-1\right)^2
\end{align}

\begin{align}
f_s = 39.9\,f_p,
\end{align}

\begin{align}
\epsilon_0 &= 77.66 + 103.3\,\left(\theta-1\right),\\
\epsilon_1 &= 0.0671\,\epsilon_0, \\
\epsilon_2 &= 3.52
\end{align}

\begin{align}
\theta = \frac{300}{T}.
\end{align}

As discussed in <cit.>, $h_i$ with $i=1,2$, follows a logarithmic distribution with probability density function (PDF) and cumulative distribution function (CDF) that can be, respectively, obtained as

\begin{align}
f_{h_i}(x) = \frac{\zeta_i^{k_i}}{\Gamma(k_i)} x^{k_i-1} \left(\ln\left(\frac{1}{x}\right)\right)^{k_i-1},\, \text{with } 0<x\leq 1
\label{Eq:f_h_i}
\end{align}

\begin{align}
F_{h_i}(x) = \frac{\Gamma\left(k_i,\zeta_i\,\ln\left(\frac{1}{x}\right)\right)}{\Gamma(k_i)},\quad\quad\quad\quad \text{with } 0<x\leq 1,
\label{Eq:F_h_i}
\end{align}

\begin{align}
\zeta_i = \frac{4.343}{\beta_i\,d_i}.
\label{Eq:zeta}
\end{align}

In (<ref>)–(<ref>),  $k_1$ and $\beta_1$ stands for the foggy conditions of the TX-HRIS link, while $k_2$ and $\beta_2$ stands for the foggy conditions of the HRIS-RX link.  For instance, $(k_i, \beta_i)$, with $i\in\{1, 2\}$, refers to light fog, if $(2.32, 13.12)$, to moderate fog, if $(5.49, 12.06)$, to thick fog, if $(6, 23)$, and to dense fog, if $(36.06, 11.91)$.

According to <cit.>, the distortion noise due to the TX hardware imperfections can be modeled as a zero-mean complex Gaussian random process with variance that can be obtained as

\begin{align}
\sigma_{n_t}^2 = \kappa_t^2 \, P_s,
\end{align}

where $\kappa_t$ is the TX error vector magnitude (EVM) and $P_s$ is the transmission power.
Similarly, for a given channel realization, the distortion noise due to the RX hardware imperfections can be modeled as a zero-mean complex Gaussian random process with variance that can be expressed as

\begin{align}
\sigma_{n_r}^2 = \kappa_r^2 \, h_g^2\,h_1^2\,h_2^2\, P_s,
\end{align}

where $\kappa_r$ is the RX EVM. Of note, the EVM in THz wireless systems is in the range of $[0.07, 0.4]$ <cit.>.

§ PERFORMANCE ANALYSIS

This section is devoted to analyze the performance of HRIS-empowered THz wireless systems. In this direction, Section <ref> provides a closed-form expression for the HRIS-empowered THz wireless system geometric losses. Section <ref> presents the statistical characterization of the end-to-end channel in terms of PDF and CDF. The outage probability of the system is presented in Section <ref>, while the throughput is given in Section <ref>.

§.§ Geometric Losses

From (<ref>), it becomes evident that the geometric losses can be obtained as

\begin{align}
P_L = \frac{1}{h_g^2},
\end{align}

or equivalently

\begin{align}
P_L = \frac{1}{h_{g,f}^2}\,\frac{1}{h_{g,m}^2}\,\frac{1}{h_f^2}.
\end{align}

Thus, from (<ref>), (<ref>), and (<ref>), we observe that the end-to-end geometric loss depends on the relative position of the HRIS in respect to the TX, the HRIS size, the TX and RX antenna gains, the transmission frequency, as well as the atmospheric conditions, i.e. temperature and liquid water density.

§.§ Statistical Characterization of the End-to-End Channel

\begin{align}
A = h_1\,h_2,
\end{align}

the following theorem returns the PDF of $A$.

The PDF and CDF of $A$ can be, respectively, expressed as

\begin{align}
f_{A}(x) &= \frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{\Gamma\left(k_1+k_2\right)}\,x^{\zeta_2-1}\,\left(-\ln\left(x\right)\right)^{k_1+k_2-1}
\nonumber \\ & \times\,_1F_1\left(k_1;k_1+k_2;(\zeta_1-\zeta_2)\ln(x)\right)
\label{Eq:f_A}
\end{align}

\begin{align}
\frac{\left(\ln\left(x\right)\right)^{k_1+k_2}}{\Gamma\left(k_1\right) \Gamma\left(k_2\right)}
\nonumber \\ & \times
\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1; \left(\zeta_1-\zeta_2\right)\ln(x),\zeta_2\ln(x)\right),
\label{Eq:F_A}
\end{align}

with $x\in(0,1]$.

For brevity, the proof of Theorem 1 is given in Appendix A.

The PDF of the end-to-end channel coefficient

\begin{align}
A_t = h_g\,A,
\end{align}

can be evaluated as

\begin{align}
f_{A_{t}}(x) = \frac{1}{h_{g}}\, f_{A}\left(\frac{x}{h_g}\right).
\label{Eq:f_A_t}
\end{align}

By applying (<ref>) to (<ref>), we obtain

\begin{align}
f_{A_t}(x) &= \frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{h_g\Gamma\left(k_1+k_2\right)}\,\left(\frac{x}{h_g}\right)^{\zeta_2-1}\,\left(-\ln\left(\frac{x}{h_g}\right)\right)^{k_1+k_2-1}
\nonumber \\ & \times\,_1F_1\left(k_1;k_1+k_2;(\zeta_1-\zeta_2)\ln(\frac{x}{h_g})\right).
\label{Eq:f_At}
\end{align}

Additionally, the CDF of $A_t$ can be expressed as

\begin{align}
F_{A_t}(x) = F_{A}\left(\frac{x}{h_g}\right),
\end{align}

which, by applying (<ref>), yields

\begin{align}
\frac{\left(\ln\left(\frac{x}{h_g}\right)\right)^{k_1+k_2}}{\Gamma\left(k_1\right) \Gamma\left(k_2\right)}
\nonumber \\ & \times
\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1; \left(\zeta_1-\zeta_2\right)\ln\left(\frac{x}{h_g}\right),
\right. \nonumber \\ & \left.
\zeta_2\ln\left(\frac{x}{h_g})\right)\right).
\end{align}

§.§ Outage Probability

From (<ref>), the signal-to-distortion-plus-noise-ratio (SDNR) can be obtained as

\begin{align}
\gamma = \frac{h_g^2\,A^2\,P_s}{h_g^2\,A^{2}\left(\kappa_t^2+\kappa_r^2\right)\,P_s + \sigma_n^2},
\end{align}

or equivalently

\begin{align}
\gamma = \frac{A^2}{A^2\left(\kappa_t^2+\kappa_r^2\right)+\frac{1}{\rho}},
\label{Eq:gamma}
\end{align}

\begin{align}
\rho = \frac{h_g^2\,P_s}{\sigma_n^2}.
\end{align}

The outage probability is defined as

\begin{align}
P_o\left(\gamma_{\rm{th}}\right) = \Pr\left(\gamma\leq\gamma_{\rm{th}}\right),
\label{Eq:Po_def}
\end{align}

where $\gamma_{\rm{th}}$ stands for the SNR threshold.
The following proposition returns a closed-form expression for the outage probability.

The outage probability can be evaluated as in (<ref>), given at the top of the next page.

\begin{align}
&P_o\left(\gamma_{\rm{th}}\right) = (-1)^{k_1+k_2=1}\frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{k_1+k_2}\,
\frac{\left(\ln\left(\sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}\right)\right)^{k_1+k_2}}{\Gamma\left(k_1\right) \Gamma\left(k_2\right)}
\nonumber \\ & \times
\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1; \left(\zeta_1-\zeta_2\right)\ln\left(\sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}\right),\zeta_2\ln\left(\sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}\right)\right)
\label{Eq:P_o_final}
\end{align}

In (<ref>), the following condition should hold:

\begin{align}
\gamma_{\rm{th}}\leq \min\left(\frac{1}{\kappa_t^{2}+\kappa_r^2}, \frac{\rho}{\left(\kappa_t^2+\kappa_r^2\right)\rho+1}\right).
\label{Eq:gamma_th_inequality}
\end{align}

Otherwise, $P_o\left(\gamma_{\rm{th}}\right)=1$.

For brevity, the proof of Theorem 1 is given in Appendix B.

§.§ Throughput

The throughput is defined as

\begin{align}
D = W\, r_{t}\,\left(1-P_o\left(\gamma_{\rm{th}}\right)\right),
\label{Eq:D}
\end{align}

where $W$ and $r_{t}$ are respectively the bandwidth and the spectral efficiency of the modulation and coding scheme that is used. Note that $r_t$ and $\gamma_{\rm{th}}$ are connected through

\begin{align}
r_t = \log_2\left(1+\gamma_{\rm{th}}\right).
\end{align}

As a result, (<ref>) can be rewritten as

\begin{align}
D = W\, r_{t}\,\left(1-P_o\left(2^{r_t}-1\right)\right).
\label{Eq:D1}
\end{align}

By applying (<ref>) to (<ref>), we obtain (<ref>), given at the top of the next page.

\begin{align}
D = W\, r_{t}\,&\left(1-
\frac{\left(\ln\left(\sqrt{\frac{2^{r_t}-1}{\rho}}\,\frac{1}{\sqrt{1-\left(2^{r_t}-1\right)\left(\kappa_t^2+\kappa_r^2\right)}}\right)\right)^{k_1+k_2}}{\Gamma\left(k_1\right) \Gamma\left(k_2\right)} \right.
\nonumber \\ & \left. \times
\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1;
\left(\zeta_1-\zeta_2\right)\ln(\sqrt{\frac{\left(2^{r_t}-1\right)}{\rho}}\,\frac{1}{\sqrt{1-\left(2^{r_t}-1\right)\left(\kappa_t^2+\kappa_r^2\right)}}),
\right.\right. \nonumber \\ & \left.\left.
\zeta_2\ln(\sqrt{\frac{\left(2^{r_t}-1\right)}{\rho}}\,\frac{1}{\sqrt{1-\left(2^{r_t}-1\right)\left(\kappa_t^2+\kappa_r^2\right)}})\right) \right)
\label{Eq:D2}
\end{align}

Notice that (<ref>) is valid, if and only if (<ref>) is satisfied. Or equivalently, if

\begin{align}
r_{\rm{th}}\leq & \min\left(\log_2\left(\frac{1}{\kappa_t^{2}+\kappa_r^2}+1\right), \right.
\nonumber \\ & \left.\log_2\left(1+\frac{\rho}{\left(\kappa_t^2+\kappa_r^2\right)\rho+1}\right)\right).
\label{Eq:r_th_inequality}
\end{align}

Otherwise, the throughput is equal to $0$.

§ RESULTS & DISCUSSION

This section focuses on  verifying the  theoretical framework by means of  Monte Carlo simulations and extracting engineering insights. In what follows, lines stand for analytical results, while markers are used for simulations. Unless otherwise stated, black, red, orange, and blue colors are used for light, moderate, thick, and dense fog, respectively.

Geometric loss as a function of $f$ for different values of $l_h=l_v$.

Fig. <ref> depicts the geometric loss as a function of $f$ for different values of $l_h=l_v$, assuming $G_{t}=G_r=50\,\rm{dBi}$, $\psi=\pi/4$, $T=20^{o}\rm{C}$, $P=101300.0\,\rm{Pa}$, $d_1=d_2=50\,\rm{m}$ and water vapour density equals to $7.5\,\rm{g/m^3}$. Notice that a water vapour density that is equal to $7.5\,\rm{g/m^2}$ corresponds to moderate fog.  From this figure, we observe that in the range of $100$ to $1000\,\rm{GHz}$, $10$ transmission windows exist. In more detail, the first transmission windows is from $100$ to $180\,\rm{GHz}$; there, the available bandwidth is equal to $80\,\rm{GHz}$. The second transmission window is from $192$ to $322\,\rm{GHz}$. The available bandwidth of the second transmission window is $130\,\rm{GHz}$. Notice that although the available bandwidth of the second transmission window is greater than the one of the first transmission window, the average geometric loss in the second transmission window is also higher than the one of the first  transmission window. The range of the third transmission window is $335-373\,\rm{GHz}$. As a consequence, the available bandwidth of the third transmission window is $38\,\rm{GHz}$. The forth transmission window is from $390$ to approximately $443\,\rm{GHz}$. Thus, the available bandwidth of the forth transmission window is $53\,\rm{GHz}$. The fifth transmission window is from $456$ to $472\,\rm{GHz}$, while the sixth transmission window is from $477$ to $487\,\rm{GHz}$. The seventh transmission window is from $587$ to $616\,\rm{GHz}$. The range of the eighth transmission window is from $630$ to $715\,\rm{GHz}$. The ninth transmission window is from $790$ to $890\,\rm{GHz}$, while the tenth is from $950$ to $970\,\rm{GHz}$. Finally, from a given $f$, as $l_h=l_v$ increases, i.e., as the size of the HRIS increases, the geometric loss decreases.

Geometric loss vs $d_1$ for different values of $f$.

Fig. <ref> demonstrates the geometric loss as a function of $d_1$  for different values of $f$, assuming $d_1+d_2=100\,\rm{m}$, $G_{t}=G_r=50\,\rm{dBi}$, $l_h=l_v=1\,\rm{m}$, $\psi=\pi/4$, $T=20^{o}\rm{C}$, $P=101300.0\,\rm{Pa}$, and water vapour density equals to $7.5\,\rm{g/m^3}$. Note that the selected frequencies are within the transmission windows. From this figure, we observe that for fixed frequency and $d_1<\frac{d_1+d_2}{2}$, as $d_1$ increases, the geometric loss also increases, while for $d_1>\frac{d_1+d_2}{2}$, as $d_1$ increases, the geometric loss decreases. Likewise, for a given $f$, the maximum geometric loss is at $d_1=\frac{d_1+d_2}{2}$. For example, for $f=100\,\rm{GHz}$, as $d_1$ increases from $1$ to $10\,\rm{m}$, the geometric loss increases from $21.72$ to $40.89\,\rm{dB}$, while, for the same $f$, as $d_1$ increases from $90$ to $99\,\rm{m}$, the geometric loss increases from $40.89$ to $21.72\,\rm{dB}$. From this example, it becomes evident that for a given frequency, the same geometric loss is achieved for $d_1$ and $d_t-d_1$, where $d_t$ is the total distance. Moreover, we observe that based on the geometric loss, the optimal position of the HRIS is as near to the TX or the RX as possible.  Finally, for a given $d_1$, as $f$ increases, the geometric loss increases.

Outage probability vs of $f$ for different values of $\frac{P_s}{\sigma_n^2}$.

In Fig. <ref>, the outage probability is plotted as a function of $f$, for different values of $\frac{P_s}{\sigma_n^2}$, assuming  $G_{t}=G_r=50\,\rm{dBi}$, $\psi=\pi/4$, $T=20^{o}\rm{C}$, $P=101300.0\,\rm{Pa}$, $d_1=d_2=50\,\rm{m}$,  $l_h=l_v=1\,\rm{m}$, and ideal RF front-end at both the TX and the RX. Moreover, the water vapour density is set to $7.5\,\rm{g/m^3}$. Notice that these conditions corresponds to moderate fog; thus, $k_1=k_2=5.49$ and $\beta_1=\beta+2=12.06$. As expected, for a given $f$, as $\frac{P_{s}}{\sigma_n^2}$ increases, the outage performance improves. For example, for $f=200\,\rm{GHz}$, the outage probability decreases from $1$ to $0.36$, as $\frac{P_{s}}{\sigma_n^2}$ increases from $60$ to $80\,\rm{dB}$.
Additionally, for a fixed $\frac{P_{s}}{\sigma_n^2}$, a local maximum in the outage probability is observed, at the water resonating frequencies.
Likewise, from this figure, we observe that based on the application reliability requirement, i.e., the maximum allowed outage probability, as well as the required energy consumption, which is translated to $\frac{P_{s}}{\sigma_n^2}$, a different range of frequencies can be employed. For instance, if the maximum allowed outage probability is set to $10^{-6}$, which is a realistic value for backhauling scenarios, and the transmission SNR is set to $100\,\rm{dB}$, the transmission frequency should be in the range of $[100, 151.5\,\rm{GHz}]$. Note that this frequency range is a subset of the first transmission frequency window, which is defined in Fig. <ref>. For the same reliability requirement,  and for a transmission SNR equal to $120\,\rm{dB}$, we observe that the range of the transmission frequency increases to $[100, 309.5\,\rm{GHz}]$. In other words, as the transmission SNR increases, the range  of the transmission frequencies that can be used to support the system increases. As discussed in <cit.>, this highlights the importance of adopting new type of transceivers in THz wireless systems that minimizes their noise figure.

Outage probability vs $\frac{\rho}{\gamma_{\rm{th}}}$ for different values of $d_1=d_2=d$, assuming light (black), moderate (red), thick (orange), and dense (blue) fog.

Fig. <ref> illustrates the outage probability as a function of $\frac{\rho}{\gamma_{\rm{th}}}$ for different values of $d_1=d_2=d$ and fog conditions. Moreover, we assume that the transceivers are equipped with ideal RF front-end, i.e. $\kappa_t=\kappa_r=0$. As expected for given fog conditions and $d$, as $\frac{\rho}{\gamma_{\rm{th}}}$ increases, the outage performance improves. For example for light fog and $d=30\,\rm{m}$, as $\frac{\rho}{\gamma_{\rm{th}}}$ increases from $5$ to $10\,\rm{dB}$, the outage probability decreases for about $2$ orders of magnitude. Moreover, for fixed fog conditions and $\frac{\rho}{\gamma_{\rm{th}}}$, as $d$ increases, the outage probability also increases. For instance, for light fog and  $\frac{\rho}{\gamma_{\rm{th}}}=15\,\rm{dB}$, the outage probability increases from $2.08\times 10^{-5}$ to $7.63\times 10^{-3}$, as $d$ increases from $30$ to $50\,\rm{m}$. Finally, for given $\frac{\rho}{\gamma_{\rm{th}}}$ and $d$, as the fog becomes more severe, the outage performance degrades. For example, for  $\frac{\rho}{\gamma_{\rm{th}}}=15\,\rm{dB}$ and $d=30\,\rm{m}$, as the fog conditions changes from light to moderate, the outage probability increases from $2.08\times 10^{-5}$ to $7.15\times 10^{-3}$. For the same  $\frac{\rho}{\gamma_{\rm{th}}}$ and $d$, the system under thick fog, achieves an outage probability that is equal to $5.9\times 10^{-1}$, while, under dense fog, it is equal to $1$. These examples indicate the importance of accounting for the fog impact on the performance of HRIS-empowered THz wireless systems.

Outage probability vs $d_1$ for different values of $\frac{\rho}{\gamma_{\rm{th}}}$, assuming light (black), moderate (red), thick (orange) fog and $d_1+d_2=100\,\rm{m}$.

Fig. <ref> depicts the outage probability as a function of $d_1$ for different values of $\frac{\rho}{\gamma_{\rm{th}}}$ and fog conditions, assuming $d_1+d_2=100\,\rm{m}$ and ideal RF front-end. From this figure, we observe that for given $\frac{\rho}{\gamma_{\rm{th}}}$ and fog conditions, for $d_1<\frac{d_1+d_2}{2}$, as $d$ increases the outage probability decreases. For example, for $\frac{\rho}{\gamma_{\rm{th}}}=40\,\rm{dB}$ and light fog, as $d_1$ increases from $10$ to $20\,\rm{m}$, the outage probability decreases from $2.17\times 10^{-6}$ to $4.35\times 10^{-7}$. On the other hand, for fixed $\frac{\rho}{\gamma_{\rm{th}}}$ and fog conditions, for $d_1>\frac{d_1+d_2}{2}$, as $d_1$ increases, the outage probability also increases. For instance, $\frac{\rho}{\gamma_{\rm{th}}}=40\,\rm{dB}$ and light fog, as $d_1$ increases from $80$ to $90\,\rm{m}$, the outage probability increases from  $4.35\times 10^{-7}$ to $2.17\times 10^{-6}$. For given fog conditions and $\frac{\rho}{\gamma_{\rm{th}}}$,the outage probability achieves the minimum value for $d_1=\frac{d_1+d_2}{2}$. In other words, the optimal place for the HRIS is at $\frac{d_1+d_2}{2}$. Moreover, we observe that the same outage performance are achieved for $d_1$ and $d_t-d_1$, where $d_t=d_1+d_2$. For example, for $\frac{\rho}{\gamma_{\rm{th}}}=40\,\rm{dB}$ and light fog, for both $d_1=10\,\rm{m}$ and $d_1=90\,\rm{m}$, the outage probability is equal to $2.17\times 10^{-6}$. Additionally, for given $d$ and fog conditions, as $\frac{\rho}{\gamma_{\rm{th}}}$ increases, the outage probability decreases. For instance, for $d=10\,\rm{m}$ and light fog, as $\frac{\rho}{\gamma_{\rm{th}}}$ changes from $20$ to $40\,\rm{dB}$, the outage probability decreases for approximately $5$ orders of magnitude. Finally, for fixed $d$ and $\frac{\rho}{\gamma_{\rm{th}}}$, as the fog density increases, the outage performance of the system degrades.

Outage probability vs $\kappa$ for different values of ${\gamma_{\rm{th}}}$, assuming light (black), moderate (red), thick (orange) fog and $\rho=30\,\rm{dB}$ (a) and $\rho=40\,\rm{dB}$ (b).

Fig. <ref> demonstrates the impact of transceiver hardware imperfections on the outage performance of HRIS-empowered THz wireless systems. In more detail, the outage probability is plotted as a function of $\kappa=\sqrt{\kappa_t^2+\kappa_r^2}$ for different values of $\gamma_{\rm{th}}$, $\rho$, as well as fog conditions, assuming $d_1=d_2=50\,\rm{m}$. As expected, for given  $\gamma_{\rm{th}}$, $\rho$ and fog conditions, as $\kappa$ increases, the impact of transceivers hardware imperfections increases; thus, the outage probability increases. For example, for $\gamma_{\rm{th}}=10\,\rm{dB}$, $\rho=30\,\rm{dB}$, and light fog, the outage probability increases for more than $2$ orders of magnitude, as $\kappa$ increases from $0.1$ to $0.3$. This indicates the importance of taking into account the impact of transceiver hardware imperfections, when designing modulation and coding schemes. From this figure, it become evident that, for given $\kappa$, $\rho$, and fog conditions, as $\gamma_{\rm{th}}$ increases, the impact of hardware imperfections becomes more severe; as a consequence, an outage performance degradation is observed. For instance, for $\kappa=0.2$, $\rho=30\,\rm{dB}$, and light fog, the outage probability increases for more than $1$ order of magnitude, as $\gamma_{\rm{th}}$ increases from $5$ to $10\,\rm{dB}$. Moreover, or given $\kappa$, $\rho$, and $\gamma_{\rm{th}}$, as the density of fog increases, the outage probability also increases. For example, for  $\kappa=0.2$, $\rho=30\,\rm{dB}$, and $\gamma_{\rm{th}}=5\,\rm{dB}$, as the fog conditions change from light to moderate, the outage probability increases for more than $2$ orders of magnitude.  Finally, for given  $\kappa$, $\gamma_{\rm{th}}$, and fog conditions, as $\rho$ increases, the outage probability decreases.

$D/W$ vs $r_t$ for different values of ${\rho}$, assuming light (black), moderate (red), thick (orange) fog and $d_1=d_2=50\,\rm{m}$.

Fig. <ref> presents the throughput to bandwidth ratio as a function of $r_t$ for different values of ${\rho}$ and fog conditions, assuming ideal RF front-end and $d_1=d_2=50\,\rm{m}$. From this figure, we observe that for given $\rho$ and fog condition, an optimal spectral efficiency of the modulation and coding transmission scheme, $r_{t}^{o}$, exists that maximizes the throughput to bandwidth ratio. For $r_t<r_{t}^{o}$, as $r_{t}$ increases, $\frac{D}{W}$ also increases. On the other hand, for  $r_t>r_{t}^{o}$, as $r_{t}$ increases, $\frac{D}{W}$ decreases. Moreover, for fixed $\rho_t$ and fog conditions, as $\rho$ increases, $\frac{D}{W}$ also increases. For example, under light fog, for $r_t=8\,\rm{bits/s/Hz}$, the $D/W$ increases from $4.32$ to $7.96\,\rm{bits/s/Hz}$, as $\rho$ increases from $30$ to $40\,\rm{dB}$. Finally, for given $r_t$ and $\rho$, as the fog density increases, $D/W$ decreases.

$D/W$ vs $r_t$ for different values of ${\kappa}$, assuming light (black), moderate (red), thick (orange) fog, $d_1=d_2=50\,\rm{m}$ and $\rho=50\,\rm{dB}$.

Fig. <ref> presents the throughput to bandwidth ratio as a function of the transmission scheme spectral efficiency for different values of $\kappa$ and under different fog conditions. As a benchmark, the case of ideal RF front-end is plotted. As expected, an optimal transmission spectral efficiency, $r_{t}^{o}$, exists, for which the throughput to bandwidth ratio is maximized. For given $\kappa$, fog density, and $r_t<r_{t}^{o}$, as $r_{t}$ increases, $D/W$ also increases. For example, thick fog and $\kappa=0.07$, the optimal transmission spectral efficiency is equal to $5.5\,\rm{bits/s/Hz}$, where the achievable throughput to bandwidth ratio is $4.125\,\rm{bits/s/Hz}$. Moreover, as the transmission spectral efficiency increases from $4$ to $5\,\rm{bits/s/Hz}$, the throughput to bandwidth ratio increases from $3.58$ to $4.058\,\rm{bits/s/Hz}$. On the other hand, for $r_{t}>r_{t}^{o}$, as the transmission spectral efficiency increases,  the throughput to bandwidth ratio decreases. For instance, for thick fog and $\kappa=0.07$, as the transmission spectral efficiency increases from $6$ to $7\,\rm{bits/s/Hz}$, the throughput to bandwidth ratio decreases from $4.01$ to $2.74\,\rm{bits/s/Hz}$. Likewise, for given fog density and transmission spectral efficiency, as $\kappa$ increases,  the throughput to bandwidth ratio decreases. For example, for thick fog and transmission spectral efficiency equal to $6\,\rm{bits/s/Hz}$,  the throughput to bandwidth ratio decreases from  $4.01$ to $3.26\,\rm{bits/s/Hz}$, as $\kappa$ increases from $0.07$ to $0.1$. Additionally, for fixed transmission spectral efficiency and $\kappa$, as the fog density increases, the throughput to bandwidth ratio decreases. For instance, for $r_t=6\,\rm{bits/s/Hz}$ and $\kappa=0.07$, if the fog condition changes from light to thick, the throughput to bandwidth ratio will decrease from $6.0$ to $3.26\,\rm{bits/s/Hz}$. This indicates the importance of accurately modeling the fog conditions, when estimating the throughput performance of HRIS-empowered wireless THz systems. Finally, from this figure, it becomes apparent that a maximum transmission spectral efficiency exists beyond which the throughput to bandwidth ratio becomes equal to $0$. The maximum transmission spectral efficiency depends on the level of transceiver hardware imperfections. For example, for $\kappa=0.07$, it is equal to $7.68\,\rm{bits/s/Hz}$, while, for $\kappa=0.1$, it is equal to $6.66\,\rm{bits/s/Hz}$. Note that the maximum transmission spectral efficiency is independent of the fog conditions. Also, we observe that the maximum transmission spectral efficiency is lower than the optimal  transmission spectral efficiency for the case of ideal RF front-end. This highlights the importance of accounting for the impact of hardware imperfections, when selecting the transmission scheme.

$D/W$ vs $\psi$ for different values of $\kappa_t=\kappa_r$, assuming $f=100\,\rm{GHz}$ (black) and $f=370\,\rm{GHz}$ (red).

In Fig. <ref>, $D/W$ is illustrated as a function of $\psi$ for different levels of transceiver hardware imperfections and frequencies, assuming $\frac{P_s}{\sigma_n^2}=120\,\rm{dB}$, $r_t=5\,\rm{bits/s/Hz}$, $G_{t}=G_r=50\,\rm{dBi}$, $\psi=\pi/4$, $T=20^{o}\rm{C}$, $P=101300.0\,\rm{Pa}$, $d_1=d_2=50\,\rm{m}$,  $l_h=l_v=1\,\rm{m}$, and water vapor equal to $7.5\,\rm{g/m^3}$, which corresponds to moderate fog conditions. As a benchmark, the ideal RF front-end scenario is also plotted. From this figure, we observe that for given $\kappa_t=\kappa_r$ and $f$, as $\psi$ increases, $D/W$ decreases. For example, for $\kappa_t=\kappa_r=0$ and $f=370\,\rm{GHz}$, as $\psi$ increases from $10^o$ to $45^{o}$, $D/W$ decreases from $4.55$ to $3.84\,\rm{bits/s/Hz}$. This indicates that not only the relative distance from the TX and the RX, but also the orientation of the HRIS plays a pivotal role on the throughput performance of HRIS-empowered wireless THz systems. Additionally, it is become apparent that, for fixed $\psi$ and $\kappa_t=\kappa_r$, as the geometric loss increases through a frequency increase, $D/W$ decreases. For instance, for $\psi=45^{o}$ and $\kappa_t=\kappa_r=0$, $D/W$ decreases from $5.0$ to $3.84\,\rm{bits/s/Hz}$, as $f$ increases from $100$ to $370\,\rm{GHz}$. This example reveals the importance of appropriately allocating the frequencies in order to make the most out of the available bandwidth. Finally, for given $f$ and $\psi$, as the level of transceiver hardware imperfections increase, $D/W$ decreases.

$D/W$ vs $d_1$ for different values of $\psi$ assuming $\kappa_t=\kappa_r=0$ (black), $\kappa_t=\kappa_r=0.1$ (red) and $d_1+d_2=100\,\rm{m}$.

Fig. <ref> depicts $D/W$ as a function of $d_1$ for different values of $\psi$ and levels of hardware imperfections, assuming moderate fog conditions, $d_1+d_2=100\,\rm{m}$, $\frac{P_s}{\sigma_n^2}=80\,\rm{dB}$, $f=100\,\rm{GHz}$, $r_t=5\,\rm{bits/s/Hz}$, $G_{t}=G_r=50\,\rm{dBi}$, $\psi=\pi/4$, $T=20^{o}\rm{C}$, $P=101300.0\,\rm{Pa}$, and  $l_h=l_v=1\,\rm{m}$. The water vapor is set to $7.5\,\rm{g/m^3}$. From this figure, we observe that for given $\psi$, $\kappa_t=\kappa_r$ and $d_{1}<\frac{d_1+d_2}{2}$, as $d_1$ increases, the outage probability also increases; thus, $D/W$ decreases. On the other hand, for fixed $\psi$, $\kappa_t=\kappa_r$ and $d_{1}>\frac{d_1+d_2}{2}$, as $d_1$ increases, the outage probability decreases; hence, $D/W$ increases. For $d_1=\frac{d_1+d_2}{2}$, $D/W$ is minimized. Notice that this observation contradicts Fig. <ref>. This is because in Fig. <ref>, the impact of geometric losses was neglected. This indicates the importance of taking into account both the deterministic and stochastic phenomena  that affect the performance of the HRIS-empowered wireless THz system. Moreover, it becomes evident that the maximum $D/W$ is achieved for $d_1\to 0$ and $d_1\to d_1+d_2$. In other words, the optimal placement of the HRIS is as near either the TX or the RX as possible. Of note, this observation is in line with the results of <cit.>. Likewise, we observe that for fixed $d_1$ and  $\kappa_t=\kappa_r$, as $\psi$ increases, $D/W$ decreases. For example, for $d_1=50\,\rm{m}$ and $\kappa_t=\kappa_r=0.1$, as $\psi$ increases from $45^{o}$ to $75^{o}$, $D/W$ decreases about two orders of magnitude. Finally, the deterministic impact of hardware imperfections is revealed by this figure. In more detail, for given $d_1$ and $\psi$, as the level of transceiver hardware imperfections increases, $D/W$ decreases.

§ CONCLUSION

In this paper, we have investigated the outage and throughput performance of HRIS-empowered THz wireless systems that suffer from the joint impact of transceiver hardware imperfections and fog. In more detail, after providing the methodology to evaluate the end-to-end geometric losses and characterizing the stochastic nature of the channel in terms of PDF and CDF, we have extracted novel closed-form expressions for the system's outage probability and throughput. These expressions have revealed a number of engineering insights, such as the existence of transmission windows, the optimum placement of the HRIS, the need to design transceivers with low noise figure, as well as the severe impact of transceiver hardware imperfections and fog. Therefore, they are expected to play a key role in the design of such systems.

§ APPENDICES

§ APPENDIX A

§ PROOF OF THEOREM 1

Since $h_1$ and $h_2$ are independent random variables (RVs), the PDF of $A$ can be expressed as

\begin{align}
f_A(x) = \int_{x}^{1} \frac{1}{y} f_{h_1}(y)\,f_{h_2}\left(\frac{x}{y}\right)\,\rm{d}y.
\label{Eq:f_A_def}
\end{align}

By applying (<ref>) into (<ref>), we obtain

\begin{align}
f_A(x) = \frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{\Gamma(k_1)\,\Gamma(k_2)} x^{\zeta_2-1}\,&\int_{x}^{1} y^{\zeta_1-\zeta_2-1} \left(\ln\left(\frac{1}{y}\right)\right)^{k_1-1}
\nonumber \\ & \times
\left(\ln\left(\frac{y}{x}\right)\right)^{k_2-1}\,\rm{d}y,
\end{align}

or equivalently

\begin{align}
f_A(x) = \frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{\Gamma(k_1)\,\Gamma(k_2)} &x^{\zeta_2-1}\,\int_{x}^{1} y^{\zeta_1-\zeta_2-1} \left(-\ln\left({y}\right)\right)^{k_1-1}
\nonumber \\ & \times
\left(\ln\left({y}\right)-\ln\left(x\right)\right)^{k_2-1}\,\rm{d}y.
\label{Eq:f_A_step_2}
\end{align}

By setting $\ln(y) = t$, (<ref>) can be rewritten as

\begin{align}
f_A(x) = \frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{\Gamma(k_1)\,\Gamma(k_2)} x^{\zeta_2-1}\,&\int_{\ln(x)}^{0}  (-t)^{k_1-1}\,\left(\exp(t)\right)^{\zeta_1-\zeta_2}\nonumber \\ & \times\left(t-\ln(x)\right)^{k_2-1}\,\rm{d}t,
\end{align}

which, by applying <cit.>, yields (<ref>).

The CDF of $A$ can be evaluated as

\begin{align}
F_{A}(x)=\int_{0}^{x} f_{A}(x)\,\rm{d}x,
\end{align}

which, by applying (<ref>), can be rewritten as

\begin{align}
F_{A}(x)= (-1)^{k_1+k_2=1}\frac{\zeta_1^{k_1}\,\zeta_2^{k_2}}{\Gamma\left(k_1\right) \Gamma\left(k_2\right)}\,\mathcal{J},
\label{Eq:F_A_s1}
\end{align}

\begin{align}
\mathcal{J} = \int_{0}^{x} & y^{\zeta_2-1} \left(\ln(y)\right)^{k_1+k_2-1}
\nonumber \\ & \times
\,_1F_1\left(k_1;k_1+k_2;\left(\zeta_1-\zeta_2\right)\,\ln(y)\right)\,\rm{d}y.
\label{Eq:J}
\end{align}

By setting $\ln(y)=t$, (<ref>) can be rewritten as

\begin{align}
\mathcal{J} = - \int_{\ln(y)}^{\infty}& t^{k_1+k_2-1}\,\exp\left(\zeta_2\, t\right)\nonumber \\ & \times
\,_1F_1\left(k_1;k_1+k_2;\left(\zeta_1-\zeta_2\right)\,t\right)\,\rm{d}t.
\label{Eq:J_s2}
\end{align}

With the aid of <cit.>, (<ref>) can be equivalently expressed as

\begin{align}
\mathcal{J} = - \sum_{n=0}^{\infty} \frac{(k_1)_n\left(\zeta_1-\zeta_2\right)^{n}}{(k_1+k_2)_n\, n!} \mathcal{K},
\label{Eq:J_s3}
\end{align}

\begin{align}
\mathcal{K} = \int_{\ln(x)}^{\infty} t^{k_1+k_2+n-1}\,\exp\left(\zeta_2\,t\right)\,\rm{d}t.
\label{Eq:K}
\end{align}

Next, by employing <cit.>, (<ref>) can be written as

\begin{align}
\mathcal{K}=\sum_{m=0}^{\infty} \frac{\zeta_2^m}{m!}\,\int_{\ln(y)}^{\infty} t^{k_1+k_2+n+m-1}\,\rm{d}t,
\end{align}

which, by applying <cit.>, can be rewritten as

\begin{align}
\mathcal{K}=-\sum_{m=0}^{\infty} \frac{\zeta_2^m}{m!}\,\frac{1}{k_1+k_2+m+n}\,\left(\ln(x)\right)^{k_1+k_2+m+n}.
\label{Eq:K_s2}
\end{align}

By applying (<ref>) into (<ref>), we obtain

\begin{align}
\mathcal{J}=\left(\ln\left(x\right)\right)^{k_1+k_2}\,\mathcal{L},
\label{Eq:J_s4}
\end{align}

\begin{align}
\mathcal{L}=\sum_{n=0}^{\infty}\sum_{m=0}^{\infty} & \frac{1}{k_1+k_2+m+n}\,\frac{(k_1)_n}{(k_1+k_2)_n} \frac{1}{n!}\frac{1}{m!}
\nonumber \\ & \times
\left(\left(\zeta_1-\zeta_2\right)\ln(x)\right)^{n}\,\left(\zeta_2\ln(x)\right)^{m}.
\label{Eq:L}
\end{align}

\begin{align}
\frac{1}{k_1+k_2+m+n} = \frac{\left(k_1+k_2+m+n-1\right)!}{\left(k_1+k_2+m+n\right)!},
\end{align}

or equivalently

\begin{align}
\frac{1}{k_1+k_2+m+n} = \frac{\Gamma\left(k_1+k_2+m+n-1\right)}{\Gamma\left(k_1+k_2+m+n\right)},
\end{align}

\begin{align}
\frac{1}{k_1+k_2+m+n} = \frac{(k_1+k_2)_{m+n}\Gamma(m+n)}{(k_1+k_2+1)_{m+n}\Gamma(m+n+1)},
\end{align}

which can be rewritten as

\begin{align}
\frac{1}{k_1+k_2+m+n} = \frac{(k_1+k_2)_{m+n}}{(k_1+k_2+1)_{m+n}(k_1+k_2)_1)}.
\label{Eq:frac_1_k_1_k_2}
\end{align}

Moreover, by accounting for the fact that $(k_1+k_2)_1=k_1+k_2$, (<ref>) can be written as

\begin{align}
\frac{1}{k_1+k_2+m+n} = \frac{1}{k_1+k_2}\frac{(k_1+k_2)_{m+n}}{(k_1+k_2+1)_{m+n}}.
\label{Eq:frac_1_k_1_k_2_s2}
\end{align}

By applying (<ref>) to (<ref>), we get

\begin{align}
\mathcal{L}=\frac{1}{k_1+k_2} \sum_{n=0}^{\infty}&\sum_{m=0}^{\infty} \frac{(k_1+k_2)_{m+n}}{(k_1+k_2+1)_{m+n}} \,\frac{(k_1)_n}{(k_1+k_2)_n} \frac{1}{n!}\frac{1}{m!}
\nonumber \\ & \times
\left(\left(\zeta_1-\zeta_2\right)\ln(x)\right)^{n}\,\left(\zeta_2\ln(x)\right)^{m}.
\end{align}

or equivalently

\begin{align}
\mathcal{L} = \frac{\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1; \left(\zeta_1-\zeta_2\right)\ln(x),\zeta_2\ln(x)\right)}{k_1+k_2}.
\label{Eq:L_final}
\end{align}

By applying (<ref>) to (<ref>), we obtain

\begin{align}
\nonumber \\ & \times
\rm{H}_{A}\left(k_1+k_2, k_1; k_1+k_2+1; \left(\zeta_1-\zeta_2\right)\ln(x),\zeta_2\ln(x)\right),
\label{Eq:J_s5}
\end{align}

Finally, by employing (<ref>) to (<ref>), we obtain (<ref>). This concludes the proof.

§ APPENDIX B

§ PROOF OF PROPOSITION 1

By applying (<ref>) to (<ref>), we obtain

\begin{align}
P_o\left(\gamma_{\rm{th}}\right) = \Pr\left(\frac{A^2}{A^2\left(\kappa_t^2+\kappa_r^2\right)+\frac{1}{\rho}}\leq \gamma_{\rm{th}}\right),
\end{align}

or equivalently

\begin{align}
P_o\left(\gamma_{\rm{th}}\right) = \Pr\left(A^{2}\left(1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)\right)\leq \frac{\gamma_{\rm{th}}}{\rho}\right).
\label{Eq:P_o_s2}
\end{align}

By assuming that $1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)>0$ or equivalently

\begin{align}
\gamma_{\rm{th}} < \frac{1}{\kappa_t^{2}+\kappa_r^2},
\label{Eq:gamma_th_in_1}
\end{align}

(<ref>) can be rewritten as

\begin{align}
P_o\left(\gamma_{\rm{th}}\right) = \Pr\left(A \leq \sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}\right),
\end{align}

\begin{align}
P_o\left(\gamma_{\rm{th}}\right) = F_{A}\left(\sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}\right),
\end{align}

which, by applying (<ref>), returns (<ref>).

From (<ref>), the following inequality should hold:

\begin{align}
0 < \sqrt{\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{\sqrt{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}} \leq 1
\end{align}

or equivalently

\begin{align}
0 < {\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}} \leq 1
\end{align}

From (<ref>),  ${\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}}>0$ is always true. Likewise, in order for  ${\frac{\gamma_{\rm{th}}}{\rho}}\,\frac{1}{{1-\gamma_{\rm{th}}\left(\kappa_t^2+\kappa_r^2\right)}} \leq 1$ to be true, the following condition should be satisfied:

\begin{align}
\gamma_{\rm{th}}\leq \frac{\rho}{\left(\kappa_t^2+\kappa_r^2\right)\rho+1}.
\label{Eq:gamma_th_in_2}
\end{align}

By combining (<ref>) and (<ref>), we obtain (<ref>). This concludes the proof.


[acronym]long-short acronymnohyperfirsttrue

# Lattice study of disordering of inhomogeneous condensates and the Quantum
Pion Liquid in effective $O(N)$ model††thanks: Presented at the Exited QCD
2024, Benasque, Spain, 15-19 January, 2024

Marc Winstel
Semeon Valgushev  Institut für Theoretische Physik, Goethe-Universität,
Max-von-Laue-Straße 1, D-60438 Frankfurt am Main, Germany.  Department of
Physics and Astronomy, Iowa State University,
Ames, IA, 50011, USA.

###### Abstract

In this talk, we study a scalar $O(N)$ model with a so-called moat regime – a
regime with negative bosonic wave function renormalization – using lattice
field theory. For negative bare wave function renormalization, inhomogeneous
condensates are solutions of the classical equations of motions. Using hybrid
Monte Carlo simulations we demonstrate how bosonic quantum fluctuations
disorder the inhomogeneous condensate. Instead, one finds a so-called Quantum
Pion Liquid, where bosonic correlation functions are spatially oscillating,
but also exponentially decaying.

##  1 Introduction

Scalar $O(N)$ field theories are a useful tool for the study of phase
transitions. This is especially true in the case of Quantum Chromodynamics
(QCD) at finite density where 1) Lorentz symmetry is broken explicitly and
therefore new phases can arise and 2) Monte-Carlo simulations are not
available due to the sign problem. In particular, QCD might exhibit the so-
called inhomogeneous phase (IP)– where in addition to chiral symmetry also
translational invariance is spontaneously broken [1]. In order to study this
phenomenon one can consider effective models such as Nambu-Jona-Lasinio (NJL)
[2] and scalar $O(N)$ models [3]. Instead (or in addition to) of critical end-
point of $O(4)$ universality class one might expect to find so-called Lisfhitz
point where three phases meet: disordered, ordered and the IP.

In previous studies of an NJL-type model in $1+1$ dimensions it was found that
the IP coincides with so-called moat regime [2], where the dispersion relation
of quasi-particles exhibits minimum at non-zero momentum. This unusual feature
manifests itself as periodic spatial oscillations of two-point correlation
functions. However, the connection of the moat regime to IP is not
straightforward. Quantum fluctuations are expected to weaken ordered phases
such as the IP [4, 5], and in fact they can even eliminate them. In
particular, it was argued that in $3+1$ dimensional scalar $O(N)$ model one
can find so-called Quantum Pion Liquid (Q$\pi$L) analogous to Quantum Spin
Liquid (QSL) instead of IP [3] due to transverse fluctuations. In this case,
Q$\pi$L is characterized by small yet always finite mass and is accompanied by
the moat regime. The Q$\pi$L regime can also be generated through mixing
effects of scalar and vector modes [6, 7].

In this work we set ourselves to go beyond limitations of previous works and
address the problem of inhomogeneous phase using ab-initio lattice
calculations of scalar $O(N)$ model at finite $N$.

##  2 Model and algorithm

We consider an $O(N)$ model in $3+1$ dimensions with spatial higher-derivative
terms as defined in [3] and take into account only the static Matsubara mode
$\beta\omega_{0}=0$ where $\beta=1/T$. In this case, an effective
$3$-dimensional Lagrangian of the static mode reads as

$L_{\mathrm{eff}}=\frac{Z}{2}\left(\partial_{j}\vec{\phi}\right)^{2}+\frac{1}{2M^{2}}\left(\sum_{j}\partial_{j}^{2}\vec{\phi}\right)^{2}+\frac{m^{2}}{2}\vec{\phi}^{\,2}+\frac{\lambda
N}{4}(\vec{\phi}^{\,2})^{2}.$ (1)

This Lagrangian can be analyzed in mean-field approximation and in the
large-$N$ limit. In the mean field approach, it features a Lifshitz point at
$Z=0$ and $m^{2}=0$ where three phases meet: disordered, ordered and the IP.
The ground state of the IP is a chiral spiral

$\vec{\phi}=\phi_{0}\left(\cos(k_{0}z),\sin(k_{0}z),\phi_{\perp}=\vec{0}\right)^{T}$
(2)

where $k_{0}$ and $\phi_{0}$ are such that they minimize the free energy. In
contrast, in the large-$N$ limit the Lifshitz point is gone. The IP is
replaced by a Q$\pi$L characterized by dynamically generated mass gap and
oscillating two-point function for certain regions in the
$\left(m^{2},Z\right)$ plane. The Q$\pi$L can exist for both the bare
parameter $Z<0$ but also for positive $Z$. As was argued in Ref. [3], even at
finite values of $N$ transverse fluctuations should be strong enough to
disorder the condensate and transform the IP into a Q$\pi$L.

In order to write down a Lagrangian suitable for numerical calculations we
discretize derivatives:

$[\Delta_{i}]_{x,y}=\frac{1}{12a^{2}}\left(-\delta_{y,x+2e_{i}}+16\delta_{y,x+e_{i}}-30\delta_{y,x}+16\delta_{y,x-e_{i}}-\delta_{y,x-2e_{i}}\right)$
(3)

where $e_{i}$ is a unit lattice step in the $i$-th direction, and perform
standard Hybrid Monte-Carlo (HMC) simulations. We run the calculations for
$N=1,2,4$ and lattice sizes $V=12^{3},16^{3},20^{3}$ and various values of
coupling constants. We generated around $3000-7000$ independent configurations
for each point of the phase diagram and used Jackknife algorithm for error
estimation.

##  3 Results

As discussed above, one expects different alternative scenarios to an IP from
analytical approximations to the partition function of model (1). Both the IP
as well as the Quantum Pion Liquid regime are characterized by a particular
behavior of bosonic correlations functions. Also, direct access of
$\langle\vec{\phi}(x)\rangle$ or $\langle|\vec{\phi}|(x)\rangle$ is not
suitable for detecting IPs due to destructive interference [4]. Thus, a
straightforward choice for an observable characterizing these different
regimes are the spatial correlation functions between the bosonic fields
$\vec{\phi}$

$C^{ij}(\vec{x})=\langle
c^{ij}(\vec{x})\rangle=\tfrac{1}{V}\sum_{\vec{y}}\langle\phi^{i}(\vec{y}+\vec{x})\phi^{j}(\vec{y})\rangle,$
(4)

where the sum over lattice sites $y$ is used to get more statistics. In order
to characterize the different regimes, we use fits of $C^{jj}$ for the
respective regimes:

  * •

Decaying oscillations $C^{ij}(x)\sim\delta_{ij}e^{-m_{r}x}\cos(kx)$ for the
Q$\pi$L, ordinary symmetry-restored phase (OSP) (using $k=0$) and the IP
(using $m_{r}=0$).

  * •

Algebraically decaying oscillations
$C^{ij}(x)\sim\delta_{ij}\frac{\cos(kx)}{x^{\alpha}}$ for possible quasi-long
range order (a.k.a liquid crystal)

As a criterion to compare the fit qualities we use the coefficient of
determination. In the following, we will present results for $M^{2}=1.0$,
$\lambda N=1.0$ and $V=20^{3}$. However, we note that our findings are stable
among different volumes $V=12^{3},16^{3},20^{3}$ except for very small,
negative $Z$, where the decay rates are getting small and larger volumes are
needed.

Figure 1:  Plot of $C(x)=C^{11}\left(\left(x,0,0\right)\right)$ for
$Z=-1.0,m^{2}=0.0$. The preferred fit scenario is determined using the
coefficient of determination. (left) $N=1$. (middle) $N=2$. (right) N = 4.

In fig. 1, we plot $C^{11}$, see eq. 4, as well as the preferred fit scenario
for $Z=-1.0,m^{2}=0.0$ and $N=1,2,4$. As one can see, the Quantum Pion Liquid
fit scenario is preferred independently of $N$. From Ref. [3], one would have
expected to obtain an IP at least for $N=1$ since there is no disordering
through Goldstone modes of $O(N)$ symmery breaking. When further decreasing
$Z$ the obtained exponential decay rate gets smaller in consistency with the
predictions from Ref. [3] from the large-$N$ limit.

In general, the observation of the different regimes in the
$\left(m^{2},Z\right)$ plane is similar to the large-$N$ findings [3] for all
studied $N$. For fig. 2, we depict the found regimes for $N=1$ in comparison
to the large-$N$ boundary lines. Our findings are in agreement with the
large-$N$ prediction, although one has to note that the comparison of
different fit scenarios does not yield conclusive results near any of the
phase boundaries. The phase diagram in the $\left(m^{2},Z\right)$ plane seems
to be identical for all studied values of $Z,m^{2}$ and $N$ up to the current
status of investigation.

Figure 2:  Phase diagram from fits for $N=1$. The green color denotes
Homogeneous Broken Phase (HBP), white — Ordinary Symmetric Phase (OSP) and
blue — Quantum Pion Liquid (Q$\pi$L) regions according to the large-$N$
results from Ref. [3].

###  3.1 Infinite volume limit with the external field

The above findings can be seen as an indication that the IP might not exist in
the $\left(m^{2},Z\right)$ phase diagram of eq. 1. However, one still has to
perform the thermodynamic limit, i.e. $V\rightarrow\infty$, in order to
identify phase transitions. On the lattice we only simulate finite volumes.
Thus, we are in the need of an extrapolation method of our findings to the
infinite volume. While a finite volume scaling of the exponential decay rate
$m_{r}$ in the Quantum Pion Liquid regime is plagued by fitting and
statistical errors, one can go the traditional route of studying phase
transition using an external symmetry breaking parameter. This is introduced
for $N=1$ by modifying eq. 1 using an oscillating external field

$L_{h_{0}}=L_{\mathrm{eff}}-h(x)\phi_{1}(x),\quad
h(x)=\frac{h_{0}}{\sqrt{2\pi}L\,\sigma_{0}}\sum_{n=0}^{L-1}\e^{-\frac{1}{2\sigma}\left(p_{n}-k_{0}\right)^{2}}\cos(p_{n}x_{2})$
(5)

with Gaussian distributed momenta, where $p_{n}=2\pi n/V^{1/3}$ and
$\sigma_{0}=0.1$. In order to determine the peak of the Gaussian in momentum
space, we extract $k_{0}$ from the simulations with $h_{0}=0$. For
$Z=-1.0,m^{2}=0.0$ and $V=20^{3}$ as used in fig. 3 we determine
$k_{0}=0.942\approx 3\times 2\pi/20$. Since the introduction of the external
symmetry breaking term does not only break translational but also rotational
invariance for $h_{0}\neq 0$, we expect that the correlation function still
depends on the relative differences in the $x_{0}$ and $x_{1}$ directions.
Translational symmetry breaking should then have an impact on its dependence
on the $x_{2}$ directions, i.e., one expects

$C(x_{0}-y_{0},x_{1}-y_{1},x_{2},y_{2})=C^{11}(x,y).$ (6)

In fig. 3, we plot $C(0,0,x,y)$ as a color map. For the left plot with
$h_{0}=0.04$, one can directly see the signs of translational symmetry
breaking due to external field. On the other hand, translational invariance is
restored for $h_{0}=0.01$ on the right plot of fig. 3 since there clearly is
only a dependence of $C$ on $x-y$. It is important to determine the transition
point $h_{0,c}\in(0.01,0.04)$ between these behaviors as a function of the
volume $V$ in order to clarify the fate of translation symmetry. If an
extrapolation of the computed values of $h_{0,c}$ reveals that
$h_{0,c}(V)=0.0$ for $V\rightarrow\infty$, the phase diagram of model (1)
would in turn consist of an IP for $N=1$ in the thermodynamic limit.
Alternatively, one could rigorously establish that inhomogeneous condensates
can simply be disordered by the inclusion of bosonic quantum fluctuations. It
would also be interesting to study the dependence of $C^{ij}(x,y)$ on
directions which are transverse to $h(x)$ in order to study spontaneous
breaking of the rotational symmetry. This work is planned for a future
publication.

Figure 3:  Color code of $C(0,0,x,y)$ for $Z=-2.0,m^{2}=0.0$. Note that the
scaling of the color bar differs in both plots. (left) $h_{0}=0.04$. (right)
$h_{0}=0.01$.

## References

  * [1] M. Buballa and S. Carignano, _Inhomogeneous chiral condensates_ , _Prog. Part. Nucl. Phys._ 81 (2015) 39 [1406.1367].
  * [2] A. Koenigstein, L. Pannullo, S. Rechenberger, M.J. Steil and M. Winstel, _Detecting inhomogeneous chiral condensation from the bosonic two-point function in the (1 + 1)-dimensional Gross–Neveu model in the mean-field approximation*_ , _J. Phys. A_ 55 (2022) 375402 [2112.07024].
  * [3] R.D. Pisarski, A.M. Tsvelik and S. Valgushev, _How transverse thermal fluctuations disorder a condensate of chiral spirals into a quantum spin liquid_ , _Phys. Rev. D_ 102 (2020) 016015 [2005.10259].
  * [4] J. Lenz, L. Pannullo, M. Wagner, B. Wellegehausen and A. Wipf, _Inhomogeneous phases in the Gross-Neveu model in 1+1 dimensions at finite number of flavors_ , _Phys. Rev. D_ 101 (2020) 094512 [2004.00295].
  * [5] J. Stoll, N. Zorbach, A. Koenigstein, M.J. Steil and S. Rechenberger, _Bosonic fluctuations in the $(1+1)$-dimensional Gross-Neveu(-Yukawa) model at varying $\mu$ and $T$ and finite $N$_, 2108.10616.
  * [6] M. Haensch, F. Rennecke and L. von Smekal, _Medium Induced Mixing and Critical Modes in QCD_ , 2308.16244.
  * [7] M. Winstel, _Spatially oscillating correlation functions in $\left(2+1\right)$-dimensional four-fermion models: The mixing of scalar and vector modes at finite density_, 2403.07430.


# RCSEDv2: Open-source web tools for visualization of imaging and spectral
data

Vladislav Klochkov1,2, Ivan Katkov1,3, Igor Chilingarian1,4, Kirill
Grishin5,1, Anastasia Kasparova1, Vladimir Goradzhanov1,2, Victoria Toptun1,2,
Evgenii Rubtsov1,2 and Sviatoslav Borisov1,6

###### Abstract

We present a set of open-source web tools for visualization of spectral and
imaging data, which we use in the second Reference Catalogue of Spectral
Energy Distributions of galaxies RCSEDv2 (https://rcsed2.voxastro.org/). Using
modern web frameworks Quasar and Vue.js we developed interactive viewers to
visualize spectra and SEDs of galaxies and the diagrams presenting emission
line ratios determined from the analysis of their spectra (BPT diagrams). The
viewers are built in Javascript which puts a minimum load on the server side
while providing full interactivity for the user. The use of modern web
frameworks provides full customization making the viewers easily embeddable
into web-sites of astronomical archives and databases. It also provides
compatibility with popular third-party web-tools such as Aladin Lite.

1Sternberg Astronomical Institute, Moscow State University, Moscow, Russia

2Department of Physics, Moscow State University, Moscow, Russia

3New York University Abu Dhabi, Abu Dhabi, UAE

4Center for Astrophysics – Harvard and Smithsonian, Cambridge, MA, USA

5Université de Paris, CNRS, Astroparticule et Cosmologie, F-75013 Paris,
France

6Department of Astronomy, University of Geneva, Switzerland

##  1 Introduction

The rapid development of Virtual Observatory technologies raises the demand
for customizable and easily embeddable data visualization tools, which can
acquire data from the VO and be used as web components across different
projects. With the exceedingly increasing amount of data, the need for the
interactive visualisation rises accordingly. Such visualisation tools can be
easily implemented with the help of modern web frameworks. In the second
version of Reference Catalog of Spectral Energy Distributions of galaxies
RCSEDv2 we substantially improved the data visualization solutions as compared
to the first release of RCSED (Chilingarian et al. 2017). We have created
tools for interactive data visualisation of 5.7 million spectra from publicly
available sky surveys.

##  2 General structure of the RCSEDv2 web-site

The RCSEDv2 structure is split in two parts, the frontend and the backend.
This allows us to move the data visualisation part on the user side and
eliminates the need for pre-generated templates/images/spectra. The catalog
data is stored in a form of PostgreSQL tables. From there it is processed by
the backend and then sent to the frontend as a JSON file to be read and then
displayed.

### The backend

All catalog data is combined in the PostgreSQL 12.4 RDBMS. Fast cone search is
available through PgSphere (Chilingarian et al. 2004). The database also
supports VO access interfaces via GAVO DaCHS (Demleitner et al. 2014) either
via Simple Spectrum Access Protocol or Table Access Protocol. The data is then
processed with the use of the Django Web Framework. User-submitted queries are
processed with the Django Custom Query. Our search query engine supports
logical operators, parametric search, cone search, SIMBAD-based object name
resolution and queries for the linked table columns.

Figure 1.: Examples of search queries and their results

### The frontend

The frontend part of the website is implemented using the web frameworks
Quasar and Vue.js as well as Plotly.js for the interactive data visualisation.
We implemented the tools to visualize the diagrams presenting emission line
ratios determined from the analysis of their spectra (BPT diagrams, Baldwin et
al. 1981) and interactive viewers for spectra broad-band spectral energy
distributions (SEDs). The spectrum viewer allows one to switch between the
different modes of the full spectrum fitting, i.e. different grids of stellar
population models and star formation histories as implemented in the NBursts
code (Chilingarian et al. 2007a, b), which we use for the data analysis in
RCSEDv2.

Figure 2.: BPT diagrams for J110731.23+134712.8, an AGN powered by an
intermediate-mass black hole in a star-forming galaxy (Chilingarian et al.
2018). Figure 3.: A spectrum with its best-fitting template and a broadband
SED for GMP 4348, a dwarf post-starburst galaxy in the Coma cluster (Grishin
et al. 2021).

The web frameworks we use provide compatibility with third-party web-tools
such as Aladin Lite (Boch & Fernique 2014). We created an Aladin Lite-based
component which allows one to switch between different imaging surveys and
displays the positions of all the spectra available in RCSEDv2 as clickable
markers, which forward the user to web-pages for of the corresponding objects.

Figure 4.: DECaLS images of the two low-surface brightness galaxies with the
spectra available in RCSEDv2, a giant LSB UGC 1382 and a dwarf LSB Ark 18
(Saburova et al. 2021; Egorova et al. 2021) in the embedded Aladin Lite
viewer. Both galaxies have two or more spectra in different parts of their
disks. In the bottom-right region of UGC 1382, there is a very rare compact
elliptical (cE) satellite galaxy (Chilingarian & Zolotukhin 2015) projected
onto its disk, also with the available spectrum. Such an unusual association
sheds the light on the origin of giant LSBs.

##  3 Summary

All the presented visualisation components, which we developed for RCSEDv2 can
be relatively easily embedded into any web-site providing access to
astronomical data for galaxies and/or stars. The flexibility and versatility
of modern web frameworks allows one to develop quite complex interactive
visualization solutions with a minimal investment of time and resources.
Interaction between the web components such the demonstrated integration
between Aladin Lite and the database query engine should definitely be
considered by the developers of next generation astronomical data archives at
large astronomical observatories – this brings the user experience to the next
level and, at the end, increases the science outcome of astronomical archives
and databases.

### Acknowledgments

his project is supported by the RScF Grant 19-12-00281 the Interdisciplinary
Scientific and Educational School of Moscow University “Fundamental and
Applied Space Research”.

## References

  * Baldwin et al. (1981) Baldwin, J. A., Phillips, M. M., & Terlevich, R. 1981, PASP, 93, 5
  * Boch & Fernique (2014) Boch, T., & Fernique, P. 2014, in Astronomical Data Analysis Software and Systems XXIII, edited by N. Manset, & P. Forshay, vol. 485 of Astronomical Society of the Pacific Conference Series, 277
  * Chilingarian et al. (2004) Chilingarian, I., et al. 2004, in Astronomical Data Analysis Software and Systems (ADASS) XIII, edited by F. Ochsenbein, M. G. Allen, & D. Egret, vol. 314 of Astronomical Society of the Pacific Conference Series, 225
  * Chilingarian et al. (2007a) — 2007a, in Stellar Populations as Building Blocks of Galaxies, edited by A. Vazdekis, & R. R. Peletier (Cambridge, UK: Cambridge University Press), vol. 241 of IAU Symposium, 175. arXiv:0709.3047
  * Chilingarian & Zolotukhin (2015) Chilingarian, I., & Zolotukhin, I. 2015, Science, 348, 418. 1504.06990
  * Chilingarian et al. (2018) Chilingarian, I. V., et al. 2018, ApJ, 863, 1. 1805.01467
  * Chilingarian et al. (2007b) — 2007b, MNRAS, 376, 1033. astro-ph/0701842
  * Chilingarian et al. (2017) — 2017, ApJS, 228, 14. 1612.02047
  * Demleitner et al. (2014) Demleitner, M., et al. 2014, Astronomy and Computing, 7, 27. 1408.5733
  * Egorova et al. (2021) Egorova, E. S., et al. 2021, MNRAS, 504, 6179. 2103.00211
  * Grishin et al. (2021) Grishin, K. A., et al. 2021, Nature Astronomy. 2111.01140
  * Saburova et al. (2021) Saburova, A. S., et al. 2021, MNRAS, 503, 830. 2011.01238