{'fraction_non_alphanumeric': 0.03712664209645691, 'fraction_numerical': 0.011033983765573066, 'mean_word_length': 4.491432828408245, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Robust Optimization of RIS in Terahertz under Extreme Molecular Re-radiation Manifestations

Terahertz (THz) communication signals are susceptible to severe degradation because of the molecular interaction with the atmosphere in the form of subsequent absorption and re-radiation. Recently, reconfigurable intelligent surface (RIS) has emerged as a potential technology to assist in THz communications by boosting signal power or providing virtual line-of-sight paths. However, the re-radiated energy has either been modeled as a non-line-of-sight scattering component or as additive Gaussian noise in the literature. Since the precise characterization is still a work in progress, this paper presents the first comparative investigation of the performance of an RIS-aided THz system under these two extreme re-radiation models. In particular, we first develop a novel parametric channel model that encompasses both models of the re-radiation through a simple parameter change, and then utilize that to design a robust block-coordinate descent (BCD) algorithmic framework which maximizes a lower bound on channel capacity considering imperfect channel state information. In this framework, the original problem is split into two sub-problems: a) receive beamformer optimization, and b) RIS phase-shift optimization. We also analytically demonstrate the limited interference suppression capability of a passive RIS. Our numerical results also demonstrate that slightly better throughput is achieved when the re-radiation manifests as scattering.


I. INTRODUCTION
With the standardization of 5G new radio (NR), it is now well-accepted that the traditional sub-6 GHz spectrum by itself is not sufficient to meet the ever-expanding network demands in the near future [2], [3]. This has led to the pursuit of utilizing higher frequency bands, which ultimately resulted in the recent commercialization of mmWave communication. However, with the advent of extended reality (xR) technologies, even higher data rates -up to 1 Tbit/s -for which mmWave bandwidths are not sufficient anymore, are required [4]. The reason is that the xR ecosystem imposes very stringent requirements on throughput of the wireless communication technologies sustaining it. Once realized, the xR applications are expected to revolutionize many industry sectors, including, but not limited to, healthcare, entertainment, and eCommerce. To support such applications, there has been a recent interest in exploring the possibility of utilizing the THz (0.1-10 THz) spectrum, which lies above the mmWave band [2]. Recent breakthroughs in the research of high-power THz sources [5], [6] have further increased the viability of utilizing this spectrum.
However, THz communication links are highly susceptible to blockages, both by static objects, and by dynamic objects including the users operating the VR [7]. Static blockages consistently prevent suitable quality of experience (QoE), while dynamic blockages result in a sudden decrease in throughput and are detrimental to the immersion of xR. A further impairment for THz signals is the molecular re-radiation that can manifest as either noise or NLOS component of the signal.
Inspired by the recent standardization efforts by various organizations, a potential solution is to deploy the emerging RIS technology that can create virtual LOS links to enhance throughput in situations where direct LOS links are blocked. Yet, the integration of RIS with THz communication links presents the following challenges: a) accurate characterization of the molecular re-radiation, b) channel estimation issues due to RIS, and c) subsequent operation of RIS under imperfect CSI. To overcome these challenges, we first develop a parametric THz channel model in this paper that captures both manifestations of re-radiation, and then use that model to present a novel alternating RIS optimization framework, where the RIS phase shift and receive beamformer are jointly optimized. [14] recently characterized the THz channel as a Rician channel, with the Rician factor computed from the molecule absorption coefficient. In the literature, both manifestations (i.e., as noise and as scattering) have been employed separately, and determining the prevalence of each effect is difficult. To be more specific, the exact effect will most likely exist as a combination of these two extreme situations. There is no way to accurately define the exact effect without comprehensive measurement studies. Given this, one of the reasonable things to do with the current information is to explore the two extreme circumstances and quantify their influence on the RIS performance which is what we do in this paper.
Due to the peculiarities of the THz links, the techniques developed in the beamforming literature on RIS [15]- [20] in the sub-6 GHz spectrum (that deal with the joint optimization of RIS phase-shifts and receiver beamformer) cannot be trivially extended. The interplay between RIS and THz band has been recently studied in [7], [21]- [26]. The authors of [21] proposed a sub-optimal search method to optimize RIS discrete phase-shifts while the authors of [22] jointly optimized the RIS location, phase-shift and THz sub-bands to improve system performance. A deep reinforcement learning-based algorithm has been used to optimize the reliability and rate for RIS-operated virtual reality systems in the THz band [7]. A physically consistent near-field channel model for RIS-THz systems was developed in [23] while the secrecy rate for an RISaided THz system was optimized in [24]. A particle swarm optimization-based method with limited channel estimation was used to optimize the RIS phase shifts in a THz band [25]. The error performance of an RIS-assisted low earth orbit satellite network has been analyzed in [26]. However, the above prior works studying RIS in the THz band neglected either the two possible manifestations of the molecular re-radiation, or the cumulative effect of this re-radiation along with the RIS configuration on the receiver noise [7], [21]- [26]. Further, they did not account for the natural challenge of imperfect CSI resulting from the passive nature of RIS elements, and non-cooperation from the interfering nodes.
The robust optimization of RIS-aided THz systems (against the imperfect CSI) has only been considered in a handful of recent works [27]- [29]. These works used semidefinite programming (SDP) techniques while ignoring the peculiar characteristics of RIS-THz integration.
Such techniques suffer from high computational cost that decreases the energy efficiency of the network, and hence defeat the purpose of low-cost RISs [30]. We bridge this gap by developing a parametric THz channel model that accounts for both assumptions of re-radiation, as well as three BCD-based joint optimization approaches of varying complexity for the proposed channel model under imperfect CSI. We use a lower bound on the channel capacity as our objective as the exact channel capacity is unknown for considering interference in our system model [31].
Different from [27], [28], the achievable throughput expression in our objective function assumes that the receiver only has access to imperfect CSI, which reflects the reality more precisely. Our objective function is also consistent with the discussion of the uplink spectral efficiency under imperfect CSI in [31, eq. (4.1)]. To the best of our knowledge, no comparative study exists that analyzes a jointly-optimized multi-antenna system in a THz environment with two extreme assumptions regarding molecule re-radiation.

B. Contributions
We study an RIS-aided THz system setting that consists of a single-stream transmitter (Tx) communicating with an RIS-aided multi-antenna receiver (Rx) in the THz band in the presence of potentially multiple single-stream interferers. For this setup, our objective is to jointly optimize the RIS's phase shift and receive beamformer while assuming imperfect CSI knowledge. Our key contributions in this paper are listed next.
A novel parametric THz channel. We propose a new parametric THz channel model that accounts for the following two extreme manifestations of re-radiation in the THz spectrum through a single parameter change: a) re-radiation is assumed as Gaussian noise, and b) reradiation is assumed as an NLOS component of the signal. We also characterize the cumulative effect of molecular re-radiation and the RIS configuration utilizing this parameter.
Three robust BCD algorithms. We formulate an optimization problem in which we jointly optimize the RIS's phase shift vector and receive beamformer vector with the objective of maximizing a lower bound on the channel capacity. Due to the coupling between the two sets of optimization variables (i.e., the RIS's phase shift vector and receive beamformer vector) in its objective function, the formulated problem turns out to be non-convex, and hence its global optimal solution cannot be obtained using standard convex optimization techniques. Because of that, we aim to obtain an efficient solution through the BCD algorithm. In this algorithm, we split the original problem of two sets of optimization variables into the following two sub-problems of one set of variables each: a) receive beamforming vector optimization problem, and b) RIS's phase shift optimization problem. These sub-problems are then solved in an alternative manner until they converge to an efficient solution of the original problem. As the latter sub-problem does not have a closed-form solution unlike the former, we propose three algorithms of varying complexity to solve the RIS sub-problem. First, we propose a conventional SDR approach as a baseline. Due to the high time complexity of the SDR approach, we then present the SA approach for its speed, where the expected receive signal strength is maximized rather than the original objective function. This approach provides a good sub-optimal solution when the interference power in the network is low. However, in a network with a moderate amount of interference, we can achieve better performance without sacrificing any speed by utilizing the gradient descent algorithm, which is our third proposed approach. These approaches consider the direct links of both users and interferers under imperfect CSI. Our objective function also caters to the non-robust counterpart by simply assuming no error. Finally, the time complexities associated with these approaches are explicitly derived.
System design insights. We analytically characterize the performance loss associated with the SA solution by deriving the stationary points of a one-element RIS-aided system with one interferer. This allows us to demonstrate that the passive RIS has limited capability to suppress interference when the direct link of the interferer is much stronger than the reflected link. Multiple system design insights can also be drawn from the numerical results. For example, our numerical results reveal that when re-radiation manifests as scattering, the corresponding throughput of the optimized system is slightly higher than when it manifests as noise. They also show that the gap in performance of the two cases depends on the visibility of the interferer direct links and frequency. Under perfect CSI, throughput is shown to increase linearly and logarithmically with the increasing number of RIS elements and Rx antennas, respectively. Further, we do not observe much penalty in performance by assuming the nature of molecular re-radiation in the optimization method incorrectly whenever perfect CSI is available. The results also show that the proposed robust algorithms perform better than the non-robust counterparts under imperfect CSI. In particular, our results demonstrate that the BCD-GD algorithm is superior in terms of the runtime and SER performance.
Notations: The scalar, vector and matrix are denoted by x, x and X, respectively. All the vectors are column vectors unless stated explicitly. For a matrix X, X T , X H , Tr (X), [X] i,j , Re (X) and X 0 denote its transpose, conjugate transpose, trace, (i, j)-th element, real part and positive semidefiniteness, respectively. The operation vec(X) results in a vector with every element of X. For a vector x, diag (x) denotes the diagonal matrix with the elements of x as its diagonal elements. The element-wise product is denoted by . The distribution of a standard complex normal random variable is denoted by CN (0, 1).

II. SYSTEM MODEL
We consider an RIS-aided THz system setup inspired by an indoor xR scenario with multiple Tx-Rx pairs communicating simultaneously in the same THz frequency band. The Tx of interest (Tx 0 ) is assumed to be a single-stream device without any active beamforming capabilities (like a VR/AR user equipment with a small form factor) [32], while the Rx of interest (Rx 0 ) is considered to have multiple antennas (similar to a mobile edge computing server) [7] and is assisted by a passive RIS. The multi-antenna Rx 0 has N R receive antennas, while the RIS has N elements. Additionally, there are N I co-channel single-stream users that interfere at the direction of Rx 0 . It is assumed that each of the communicating devices can have two paths to Rx 0 , one link coming directly from the Tx, and another link reflected from the RIS. The system model is illustrated in Fig. 1. Note that the N I interfering users are communicating with their own Rxs. These Rxs do not affect our analysis and are therefore not included in Fig. 1. Both RIS and Rx 0 are assumed to be uniform rectangular arrays (URAs) with half-wavelength spacing to sufficiently decrease mutual coupling [33]. The array factor of a general URA with N 0 elements is defined as follows [31]: where k(ϕ, θ) = 2π λ [cos(θ) cos(ϕ), cos(θ) sin(ϕ), sin(θ)] T is the wave vector, λ denotes the wavelength, u i denotes the vector of Cartesian co-ordinates of the i-th URA element, and U denotes the matrix [u 1 , u 2 , . . . , u N 0 ]. The azimuth angle ϕ ∈ [−π, π) is measured from the positive x-axis and the elevation angle θ ∈ [− π 2 , π 2 ) is measured from the x-y plane. The array is assumed to be on the positive y-z plane with the origin as the reference. Note that the array factor is defined as a row vector.

A. Terahertz Channel Model
Owing to the molecular absorption phenomena in the THz band, a fraction 1 − τ (f, d) of the propagating signal is absorbed and re-radiated by the media. The remaining fraction τ (f, d) = e −k(f )d is termed as the transmittance of the channel, where k(f ), f and d denote the molecular absorption coefficient, the operating frequency, and the link distance, respectively. Further, we use a simple LOS channel model that is valid for 200-450 GHz [34] for the molecular absorption coefficient calculation. The equation dictating the value of k(f ) is: where µ denotes the volume mixing ratio of water vapor, y i (f, µ) denotes the absorption coefficient for the i-th absorption line, and the polynomial function g(f, µ) is an equalization factor.
The volume mixing ratio µ can be calculated through the following equation: where φ and p denote the relative humidity and pressure, respectively. The function p w (T, p) can be described by the Buck equation [35]: p w (T, p) = 6.1121 1.0007 + 3.46 × 10 −6 p exp 17.502T 240.97 + T , where the pressure p is in hectopascals and the temperature T is in Celsius. Next, we express {y i (f, µ)} and g(f, µ) as described in [34]: where the frequency f is in Hertz and the various inner entities are described below: Recall from Section I that the absorbed power is emitted again in the same band, and there are two extreme modeling assumptions regarding how this molecular re-radiation manifests. Unlike the forward scattering we witness in visible light owing to tiny particles, this re-radiation occurs virtually isotropically, and the energy is spread in all directions for the THz band [36]. However, it is commonly assumed that the whole absorbed power is available at the Rx node through reradiation [11], [14]. This result follows from the assumptions that a single absorption/re-radiation event occurs throughout the entire propagation route, and that all re-radiated power is directed in the direction of the Rx node, despite the fact that this is never explicitly stated. In our paper, we also apply these assumptions. The assumptions of the two extremes of molecular re-radiation are provided next for a more systematic exposition.
Assumption 1. Molecular re-radiation is modeled as additive Gaussian noise.
Assumption 2. Molecular re-radiation is modeled as a scattering event where the affected channel response includes an NLOS component.
The re-radiated signal, which has 1−τ (f, d) of the total signal power appears as either additive noise or as the NLOS component. Because of the nature of their modeling, the two assumptions are at opposite extremes of the spectrum, with reality falling somewhere in the middle. The Rician factor K d in Assumption 2 for distance d is .
The two assumptions can also be interpreted as resulting from different possible receiver structures. Firstly, we have to recognize that radiation trapping generally leads to both delay and frequency dispersion, with the former arising from the finite lifetime of the excited molecular states, and the latter arising from (partial or complete) frequency redistribution occurring about relaxation of the excited state into the ground state [10]. If the frequency redistribution is over a band that is larger than the transmission bandwidth, part of the absorbed energy will be completely lost, while the remainder stays in the considered band. Depending on the amount of dispersion and its time variance, as well as the considered modulation and coding scheme, the scattered energy can be exploited for signal detection, or act as interference. For this reason, this paper considers the two limiting cases of "fully useful" and "fully noise".
By including a variable ζ in the standard Rician channel model [14], we unify both Assumptions 1 and 2 concerning the re-radiation process for analytical convenience. In the first scenario, the unified model should only include the LOS component since the re-radiation signal appears as Gaussian noise, whereas the channel model in the second case should include both the LOS and NLOS components because re-radiation appears as scattering. In other words, the models corresponding to Assumptions 1 and 2 can be recovered from the unified model by setting ζ = 1 and ζ = 0, respectively. Further details on the modified Rician channel with the ζ parameter are provided next.
The channel response between X and Y is denoted by h XY or H XY (depending on the number of antenna elements at X and Y), and can be expressed as: where denotes a random phase uniformly distributed in [−π, π), X = S or R represents RIS or the Rx 0 , and Y = R or T i represents the the Rx 0 or the i-th Tx. The distance between X and Y is denoted by d while f is the transmission frequency. The geometric channel related to the LOS path is denoted by f LOS (or F LOS ) and the NLOS counterpart is denoted byh XY (orH XY ) which is a complex vector (or matrix) with each entry being an independent and identically distributed (i.i.d.) complex normal random variable with zero mean and unit variance. Note that, when ζ = 1, the NLOS component vanishes and K d K d +1 = τ (f, d) corresponds to the transmittance, as expected. The distances and fixed channels associated with the different combinations of X and Y in the system model are compiled in Table I. We define the following stacked channels for notational ease: where I i is an indicator function that either takes 0 or 1 depending on the visibility of the direct link with probability P L i .

B. Imperfect CSI Model
In a practical wireless system, wireless channels need to be estimated before reliable communication links can be established. Unlike the estimation of channels between communication nodes, the estimation of non-cooperative interferer channels is particularly challenging and errorprone. We use an additive error model [37], [38] to model this imperfect reflected channel as follows: where the true channel is Z i , the estimated channel isẐ i , and the elements of the error matrix Similarly, the error vector for the imperfect direct channel is δ i whose every element is i.i.d as CN (0, ρ 2 i ). The imperfect direct channel is expressed as follows: Now, estimated stacked channels are defined below: We consider {ρ i , ρ i } to be parameters in our model that control the extent of the uncertainty in our channel estimates. Note that the channel is perfect when these parameters are set to zero.
Larger values of these parameters denote worse CSI quality.

C. Signal Model
If the signal x i of power E[|x i | 2 ] = P i is transmitted by the i-th Tx, the received signal at Rx 0 is expressed as (12): where θ = [e jϕ 1 . . . e jϕ N ] T is the RIS configuration vector, ϕ n ∈ [0, 2π] is the n-th entry of the vector ϕ for all n ∈ {1, . . . , N } and denotes the n-th element's reflection coefficient, and w denotes the additive Gaussian noise with variance σ 2 w + ζσ 2 m . The variance terms σ 2 w and ζσ 2 m represent the thermal noise and molecular re-radiation noise, respectively. Molecular re-radiation noise under Assumption 1 can be calculated as: is the molecular reradiation noise due to the i-th Tx. Note that, as this ζ = 0 conforms to Assumption 2, the molecular re-radiation noise variance disappears and manifests as fading. Now, we have all the information to make the signal model more compact with the stacked channel structure (8). Using that, the received signal is rewritten as: where θ 0 = [θ T 1] T . Now, we multiply the received signal with the receive beamformer u H from the left and express the resulting received signal as: While the first term in (14) represents our desired signal, the last three terms denote interference, channel estimation errors, and noise, respectively. As the exact channel capacity is unknown for the interference channel, we provide a well-known lower bound on the channel capacity C sys based on the discrete memoryless interference channel [31, Corollary 1.3]: such that the SINR term γ(u, ϕ) can be expressed as follows: where {C e i } are the co-variance matrices for the estimation errors and can be calculated as follows: where (a) follows from the fact that the errors are uncorrelated, and (b) follows from using the Lemma 1. In Assumption I, the re-radiated signal from the Tx i is modeled as additive Gaussian

Proof: See Appendix A.
We define σ 2 σ 2 m 2 ,i using Lemma 1 for convenient notation. This completes the SINR (16) description.

VECTOR
In this section, we jointly optimize the receive beamforming weights and the RIS configuration vector. We assume that Rx 0 knows the estimation error variances in the robust case. As the exact channel capacity is not computable, we set the objective function as a well-known lower bound of the channel capacity (15) under imperfect CSI. Maximizing a lower bound is useful in its own right as it effectively maximizes the channel capacity. However, maximizing log 2 (1 + γ(u, ϕ)) is equivalent to maximizing γ(u, ϕ) due to the monotonically increasing nature of logarithm function. The new objective function (16) is non-convex due to the coupling between the two sets of variables. The non-convex constraint of unit modulus makes the problem more difficult.
We employ the BCD method to solve this optimization problem by separating it into two subproblems, one for each set of variables, and solving each sub-problem in an alternating manner, to obtain an efficient solution. With a fixed RIS phase vector, the receive beamformer sub-problem can be conveniently posed as a maximization of a Rayleigh quotient problem, and hence has a simple analytical solution. However, the RIS sub-problem does not have an analytical solution.
We propose three methods of different complexities to achieve suboptimal solutions: a) SDP with Gaussian randomization, b) SA, and c) GD approach. Finally, we show the convergence and compare the complexities of the proposed algorithms.

A. Problem Formulation
Our objective is to maximize the effective SINR of Tx 0 given the estimated channel such that the norm of receive beamforming vector is unity, and the RIS elements follow unit modulus constraint. The transmit power P 0 of the user is set as the maximum allowable power to maximize the objective. Now, the optimization problem can be formulated as: 0 ≤ ϕ n < 2π, ∀n = 1, 2, . . . , N.
The general BCD algorithm used for the optimization here is shown in Algorithm 3. Note that under imperfect CSI, initializing the {ρ i , ρ i } to zero results in the non-robust counterpart of this algorithm. The non-robust algorithm treats the estimated CSI as the perfect instantaneous CSI and is used as benchmark for the performance of the robust algorithms.

B. Receive Beamformer Optimization
As noted before, for a given θ, this sub-problem can be expressed as an unconstrained maximization of a Rayleigh quotient: The normalized analytical solution for the given problem is given by [39]: where e 0 =Ĥ 0 θ 0 is the dominant eigenvector of B 0 and the solution is normalized to ensure a unit norm vector.

C. RIS Optimization through SDP
We approach this sub-problem through SDP because optimality can be ensured if we can express it as a convex optimization problem. To that end, we express (16) for a given receive beamformer u as: where Following [40], we introduce an auxiliary variable b ≥ 0 to transform (22) to an epigraph form: [Ψ] l,l = 1, ∀l = 1, 2, . . . , N + 1.
Note that if the above problem is feasible then γ(u, Ψ) * ≥ b, while the opposite condition γ(u, Ψ) * ≤ b holds when the above problem is infeasible. So, using bisection for (24) provides a good solution for Ψ. The problem (24) can be solved by any standard convex optimization package like CVX [41], [42]. However, this solution will generally not be a rank-one solution.
We extract a rank-one solution by Gaussian randomization [40], as discussed next. The unit circle projection (the division of each entry of the vector by its absolute value) of the rank-one solution provides θ * .
Choose the solution with unit circle projectionθ i+1 that provides the highest γ i+1 through (16).
Gaussian Randomization: The Gaussian randomization scheme [43] entails generation of a random vector z ∼ CN (0 N +1 , ϑ) where ϑ is the Cholesky decomposition of Ψ to extract a rankone solution. The rank-one candidate for one such vector is where z N is the vector consisting of the first N elements of z. A number of such candidates are generated to choose the one that is feasible and provides the largest objective value among them. The approximation accuracy of such a scheme in different scenarios is well-investigated by [43]. The worst-case approximation accuracy was proved to be reasonable for a finite number of generations. However, the complexity of this problem is prohibitive for RISs with large number of elements [30]. This encourages us to explore some low-complexity approaches in the subsequent subsections.

D. SA Solution to the RIS Sub-problem
A low-complexity approach to the RIS sub-problem is to maximize the received signal strength as opposed to SINR due to the existence of a closed-form solution [15]. As this approach aligns the phases of reflected signal with the phase of the direct signal, we denote this approach as SA and is expressed as follows: This can incur some performance loss from the global optimum solution. Since most RIS phase-shift optimization problems are NP-hard, a global optimum is difficult to obtain [44].
However, it has been demonstrated that the local optimal solutions can improve performance significantly [19], [20]. In that case, the performance loss can be characterized by finding good stationary points. In [19], a stationary point of a general SINR is found with respect to the phase-shift of one element without considering direct links to facilitate the element-wise BCD algorithm. Meanwhile, the authors of [18] obtain a stationary point of an SINR expression in an N -element active RIS-aided network that does not consider interference. This was possible as the active RIS elements do not have the unit modulus constraint of passive RIS. However, a similar analysis has not been done for the passive RIS. To fill this gap, we provide an analysis next for a one-element passive RIS with only one interferer while considering direct links.
One-element RIS Sub-problem: A generic SINR expression for an one-element RIS with only one interferer can be readily written as: where a 1 , b 1 , θ 1 , and c denote the reflected signal channel, reflected interferer channel, RIS phase shifts and noise. The direct signal and interferer channels are denoted by h and g, respectively.
Theorem 1. The SINR for one-element RIS takes one of the two values at the stationary points given by: where C = (L P ) 2 + (M N ) 2 − 2L M N P cos(s − t).
Proof: See Appendix B.

Remark 1.
In the simple setting without interference or P 1 → 0, the SINR γ * 1 → L +M N . This is a result of choosing θ 1 in a way to align the phase of a 1 θ 1 with h, which is the SA method. In absence of interference, the optimal solution is the SA solution.
Proof: The SA solution in this case is x = −s. The resulting SINR is denoted by γ h 2 = L +M N +P cos(s−t) . Both g 1 (k) and g 2 (k) can be calculated easily by evaluating |γ * 1 − γ h 2 |.
Remark 2. Note that g 1 (k) and g 2 (k) both tend to zero when k → ∞. In other words, the stationary points converge to the SA solution when the interferer direct link is too powerful. In a practical scenario, where the interferer direct link is as powerful as the reflected link or k → 1, (4c+c 2 )(2 cos(s−t)+2+c) implying the sub-optimality of the SA solution.
These results confirm the need for a low-complexity algorithm that finds a better solution than the SA, which inspires our next approach.

E. Gradient Descent Approach to the RIS Sub-problem
Achieving a better solution than the SA method while retaining its low-complexity benefit requires a different approach. As noted in the previous subsection, local optimal solutions can provide significant performance improvement. Gradient descent is a natural choice for such a solution because it tries to converge to a local minimum from the initial point in the chosen descent direction. Generally, the chosen descent direction is the steepest one or the negative gradient of the objective function. As shown next, this can be calculated from an alternate formulation of our RIS sub-problem in an unconstrained manner with respect to where The gradient can be calculated as follows: Remark 3. By substituting {C i } as zero vectors, we obtain the same objective function and the gradient as [30].
Armed with the analytical expression of the gradient, a simple GD algorithm works through the simple update rule: where ϕ (t) is the RIS phase vector at t-th iteration, β (t) is the step-size and ∇ ϕ ϕ (t) is the gradient. In practice, convergence time is dependent on choosing a good step-size. A step-size that is either too large or too small can result in slow convergence by either oscillating or moving too slowly in the descent direction. A practical alternative is to implement backtracking line searches. These line searches start with an initial step-size and keep diminishing them in a loop. The largest step-size in that sequence is chosen that ensures a sufficient decrease in the descent direction. Moreover, [45] shows that GD with diminishing step sizes almost always avoids saddle points under random initialization. Saddle points are the critical points that are neither local minima nor local maxima. It is also shown in [46] that saddle points can slow down GD with constant step size considerably to the extent of needing exponential time to escape, even with reasonable random initialization schemes. The objective function is highly non-convex and expected to have multiple saddle points, so being able to avoid saddle points is a desirable property. With these motivating factors, we choose Armijo-Goldstein (AG) line search [47]. This where 0 < ε < 1 is a constant. If the condition (31) is violated, the step-size is decreased by a factor of 0 < < 1 starting from a larger initial step-size. Complete details of the GD approach are shown in Algorithm 2. Note that, our proposed algorithm uses the SA solution as the initial point. With all the sub-problem solutions, the general BCD framework is demonstrated in Algorithm 3.

F. Convergence and Complexity Discussion
The BCD algorithm presented in Algorithm 3 ensures non-decreasing objective values after each iteration. As the objective function is also upper bounded by some value, Algorithm 3
converges regardless of the approach taken to solve the RIS sub-problem. However, all the approaches are of varying complexities. In this subsection, worst-time complexities are derived in the big-O notation for the algorithms when tractable and run-times are compared numerically.
First, for each iteration of the BCD algorithm, we will calculate the time complexity for each sub-problem.  a normal random variable with the covariance matrix Ψ * as per [40]. while the latter has a time complexity of O(4N 2 ). Considering the initial point calculation, the total time complexity is O(I 1 ((N I +1)(N 2 +2N R N )+4N 2 )+N 2 R N ). Now the run-times for one iteration at (N R = 100) and N = 100 are compared in Fig. 2c and compiled in Table II along with the complexities. The algorithm with the SA sub-problem solution is denoted by BCD-SA.

IV. NUMERICAL RESULTS
In this section, we provide numerical results to quantify the impact of molecular re-radiation, RIS elements, Rx antennas, interferers, and channel estimation error on key performance metrics of an RIS-aided wireless network inspired by a 3D indoor xR setup. We consider a spherical coordinate system of (r, ϕ, θ) where r denotes the distance in meters from the origin, and the angles are defined in (1). In the simulation setup, one corner of the Rx 0 arranged in a square URA pattern with 100 antennas is situated at the origin (0, 0 • , 0 • ) on the positive y-z plane while the location and orientation of the RIS (equipped with 100 elements) are translated 1m along the positive x-axis with respect to the Rx 0 location. The RIS is also arranged in a square URA pattern. The Tx 0 is placed at (1, 60 • , 0 • ) and transmits at 2 W (33 dBm) effective isotropic radiated power (EIRP) (which is within typical operating parameters of THz systems, e.g., see [51]) over a large bandwidth specified later. The coordinates of the interferer are (1.5, 110 • , 0 • ), and it transmits at the same power level. Note that, these distances are chosen to ensure that all the nodes operate in a far-field region. Unless otherwise specified, the system parameters for the simulation configuration are as follows: the transmission carrier frequency of 220 GHz, bandwidth of 10 GHz, relative humidity of 50%, standard atmospheric pressure of 1 atm, and temperature of 27 • C. These parameters, in conjunction, determine the value of k(f ). In addition, the thermal noise variance is assumed to be −174 dBm/Hz. For the general BCD algorithm, the parameter is set to be 10 −6 . The parameters 0 , 1 , and G in the SDR sub-problem are set to be 20, 10 −6 , and 1000, respectively. The parameters ε, th , β 0 , and in the GD sub-problem are set as 0.00005, 10 −6 , 1, and 0.5, respectively. We have used throughput and uncoded symbol error rate (SER) as the performance metrics. Throughput is calculated by the expression BC sys for the perfect CSI case and throughput results are averaged over 2000 iterations. SER is calculated using 4-QAM modulation scheme and SER results are averaged over 10 6 symbols. We also assume that direct signal path is completely blocked to focus on the RIS's capability. This use-case pertains to indoor THz communications, as the visibility of LoS link depends on the transmission frequency and the probability of blockage by the user's own body (self-blockage) or by other user's body (dynamic blockage) [7]. Runtime comparison of different algorithms: In Fig. 2c, we plot run-times of one iteration of the three different algorithms against different numbers of RIS elements. These run-times are averaged from the simulations needed to create Fig. 2a and 2b on a 3.59GHz AMD Ryzen 5 3600 6-Core PC with 16GB RAM. We observe that GD is the fastest algorithm while SDR is the slowest. Note that the run-time for the GD algorithm increases linearly unlike the SA and SDR algorithms. A possible explanation is that the GD sub-problem performs better in the BCD does not provide much information about the assumptions.
Effect of Rx antennas and estimation of ζ on achievable throughput: In Fig. 3a and 3b, we plot throughput versus number of Rx antennas. We investigate the case when we optimize according to Assumption 1 whereas Assumption 2 is the reality and vice-versa. In both the cases, we assume that we have the perfect CSI for the optimization procedure. Our assumption for the optimization procedure only affects the amount of noise in the SINR expression. The cases of perfect and opposite ζ estimation are shown in Fig. 3a and 3b, respectively. We observe from Effect of CSI quality on SER: We plot SER versus CSI quality and number of interferers in Fig 4 where we consider CSI quality of the signal link and the interferer links by defining σ 2 0 = η 2 1 , and σ 2 i = η 2 2 for all i = 1, 2, . . . , N I , respectively. We also conduct the simulations under Assumption 1 and disregard BCD-SDR as it only acts as a computationally expensive baseline. In Fig. 4a, we plot SER against η 2 1 . The SER curves of both robust and non-robust algorithms overlap implying that the robust algorithms have limited benefits when signal link has error. In Fig. 4b, we plot SER with interferer error η 2 2 . When interferer direct links are not present, the SER of the non-robust BCD-SA algorithm is significantly worse than the other Effect of interferers with imperfect CSI on SER: For a more practical scenario, SER is plotted against the number of interferers under no signal error and high interferer channel error (η 2 2 = 10 −11 ) in Fig. 4c. This is a practical use-case when the interferers are non-cooperative and the Rx has limited information about the interferer channels. Further, the interferers are distributed randomly in a uniform manner on a circular ring of radius 2 meter with Rx 0 as the center.
In this figure, SER performance of a robust algorithm is better than the non-robust counterpart in both the scenarios. We observe that the gap between them decreases with the number of interferers when direct links of the interferers are present, but exhibits increasing trend when direct links of the interferers are blocked. This emphasizes the usefulness of BCD-GD algorithm in the presence of non-cooperative interferers without direct links to the Rx.
Effect of transmission frequency on achievable throughput: Fig. 3c shows the variation of achievable throughput with the frequency. As there is not much gap between assumptions in the 'ND' scenario, we focus our discussion to the 'D' scenario. In this case, the performance in Assumption 2 is not affected much by the absorption coefficient peaks due to the NLOS nature of re-radiation unlike the performance in Assumption 1. The performance loss in Assumption 1 results from the high molecular re-radiation noise in the absorption peaks.

V. CONCLUSIONS
In this work, we investigated the sensitivity of an RIS-aided THz system under two extreme manifestations of molecular re-radiation. In particular, we first developed a parametric THz channel model that accounts for both re-radiation assumptions through a simple parameter change. This channel model was then utilized to jointly optimize the RIS's phase-shift and receive beamformer with the objective of maximizing a lower bound on channel capacity that considers imperfect CSI. Specifically, we proposed an alternating BCD optimization framework that splits the original problem of two sets of variables (i.e., the RIS's phase shift and receive beamformer) into two sub-problems of a single set of variables each. These two sub-problems are then solved iteratively to converge to an efficient solution. In this framework, we approached the RIS sub-problem from three different directions: a) an SDR method that reformulates the sub-problem as a convex optimization problem while relaxing some constraints, and finds a near-optimal solution with high computational complexity, b) a fast SA method that maximizes the numerator of the SINR, and c) GD-based method that converges to a first-order stationary point of the original non-convex sub-problem.
Our analytical results for a one-element RIS-aided system in the presence of a single interferer demonstrated that the SA solution is sub-optimal when the direct link of the interferer has comparable power to the reflected links. Several key system design insights were also obtained from our numerical results. For instance, our results revealed that the throughput of an optimized system is slightly higher in the scattering manifestation of re-radiation, and the exact difference depends on both the LOS probability of the direct links and frequency. On top of that, they showed that the peaks in the absorption coefficient for different frequencies have no additional impact on the throughput performance in the scattering manifestation of the reradiation. Moreover, they demonstrated that the performance loss due to incorrect assumption of the re-radiation model in the optimization procedure is minimal under perfect CSI. Further, the throughput exhibits linear and logarithmic increasing trend against the number of RIS elements and Rx antennas, respectively. The results also demonstrated that the robust algorithms provide better SER when interferer direct links are blocked. Furthermore, they highlighted the efficacy of BCD-GD algorithm through its runtime and superior SER performance in the presence of non-cooperative interferers. To the best of our knowledge, this paper makes the first attempt to investigate the performance sensitivity of an optimized RIS-assisted THz system caused by different assumptions regarding molecular re-radiation. As the numerical and analytical results show that passive RIS has limited capability to combat powerful interference, using an active RIS is left as a promising direction for future work. If the distances between Tx i to RIS, and RIS to Rx 0 are d γ i and d α , respectively, we inspect the signal x of power P i through the m-th element of RIS with reflection coefficient (α m e jθm ) without including path-loss terms for simplicity.
As the reflected power from the RIS element is |α m | 2 P i , n 2 ∼ CN (0, |α m | 2 P i (1 − τ (f, d α ))) is the additive noise for the RIS to Rx 0 path. So, the noise variance due to both the paths is E[|n 1 α m e jθm τ (f, d α ) + n 2 | 2 ] = |α m | 2 τ (f, d α )P i (1 − τ (f, d γ i )) + |α m | 2 P i (1 − τ (f, d α )) Extending this result to an N -element RIS, if the RIS-Rx 0 and Tx i -RIS channels are a h 1 , and a h 2 with their entries as array factors a h 1,m and a h 2,m with ULA assumption for RIS, the received signal for SISO is The molecular absorption noise variance is then ζσ 2 m,i , where σ 2 m,i = σ 2 m 1 ,i + N σ 2 m 2 ,i as this noise will only exist for Assumption 1 or ζ = 1.

B. Proof of Theorem 1
We start with the generic SINR expression: L + M cos(s + x) N + P cos(t + x) .
Differentiating (32) with respect to x and equating it to zero results in the following equations: L + M cos(s + x) N + P cos(t + x) = M sin(s + x) P sin(t + x) , (34) can be derived from (33). However, (33) provides an alternate form of γ 1 at a stationary point. Using (33), the SINR can be expressed as a function of sin(t + x): Note that, (34) can be expressed as a quadratic equation of sin(t+x). The solutions for sin(t+x) plugged in (35) proves the theorem.