Efficient Propagation Modeling for Communication Channels With Reconfigurable Intelligent Surfaces

We present a hybrid numerical method that combines full-wave analysis with ray-tracing to efficiently model propagation in reconfigurable intelligent surface (RIS)-enabled communication channels. Full-wave simulation is used to obtain the radiation pattern of the transmitter and the complex radar cross section (RCS) of the RIS. Then, the RIS is imported into a ray-tracer as a secondary transmitter. That enables the modeling of the interaction of the RIS with complex environments. We present the formulation of the proposed technique and its validation against full-wave (finite element) results in representative scenarios of indoor propagation in the presence of an RIS.


I. INTRODUCTION
R ECONFIGURABLE intelligent surfaces (RISs) can redirect incident waves to selected directions in a wireless communication channel, to improve coverage bypassing obstructing objects [1]. This is a promising technology for future wireless networks, particularly those operating at millimeter wave bands [2]. In the past two years, analytical work on the derivation of path loss models of RIS-enabled links [3]- [5] and investigations on their operating principles [6] have been presented.
In practice, several factors affect the performance of an RIS-enabled communication channel. For example, both the transmitter and the RIS radiate waves according to a pattern that includes a main lobe along with sidelobes, leading to enhanced multipath propagation. This effect can compromise the intended operation of an RIS and deserves careful analysis.
In [7], reflections from an RIS were computed by assigning a reflection coefficient to each unit cell of the RIS, in the direction of the receiver, and computing the associated link budget. In practice though, the RIS unit cells interact collectively rather than individually with incident wavefronts. The mutual coupling between unit cells of metasurfaces is a well-known problem that has been thoroughly studied in the literature [8]. This effect has a significant impact on RIS channel propagation models [9], because it influences the scattering properties of the RIS. This motivated the use of full-wave analysis in [10] in a free-space communication link including an RIS. Yet, the full-wave analysis of realistic indoor and outdoor environments requires excessive computational resources.
To address the need for accurate analysis of the RIS-channel interaction in a computationally feasible manner, we propose a hybrid approach that employs full-wave analysis to determine the complex radar cross section (RCS) of the RIS only. Then, the complex RCS is used to implement the RIS as a secondary radiation source in a ray-tracer [11]. Our goal is to compute the received signal strength (RSS) in complex indoor environments. This approach is associated with the standard method for modeling installed antennas, where the free-space or in situ antenna pattern is ray-traced to capture the interaction of the antenna with the environment [12]. The goal of such simulations is the computation of radiated fields by installed antennas (e.g., on ships or vehicles), using the radiation pattern of the antenna that acts as a primary source. On the other hand, we compute direct and scattered fields from the transmitter and the RIS (modeled as a secondary source through its complex RCS), along with their multiple interactions with the environment (reflections, transmissions, diffractions), to determine the total RSS over large areas.
The rest of this letter is organized as follows. Section II introduces the proposed hybrid method and explains its implementation. Section III demonstrates the accuracy of the method against finite element results in geometries where finite element analysis is feasible. Fig. 1 shows a typical scenario of RIS-enabled communications: an RIS scatters fields from the transmitter toward non-lineof-sight (NLOS) receiver points, bypassing an obstructing wall. The RIS effectively acts as a secondary transmitter for these NLOS points. This secondary transmitter can be embedded in ray-tracing, if we know its pattern and transmit power. To this end, we express the power density of the scattered field of the RIS at a receiver point, as [13]

A. RIS as a Secondary Transmitter
where S inc is the incident power density at the RIS, σ = σ(θ, ϕ) is its scalar RCS, and R 2 is the center-to-center distance from the RIS to the receiver. We approximate the incident power on the RIS as where A ph is the physical area of RIS.
1536-1225 © 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. Then, (1) becomes By inspection of (3), we conclude that the RIS creates the same power density at the receiver, as an antenna of gain and transmit power The far-field of such an antenna at a distance R can be expressed as [14] E =ê where η 0 = 120π Ω is the free-space impedance. Based on this observation, we can embed the RIS into a ray-tracer as a secondary transmitter and compute its far-field from (6), with g t and P t from (4) and (5). The scattered field from the RIS can be superimposed to the transmitted field from the transmitter, in the absence of the RIS, to determine the total electric field throughout a given environment. In the following, we discuss how to obtain: (a) g t and the phase ψ in (6) from the complex RCS [15] of the RIS and (b) P t = P RIS,inc , so that we can use (6) to compute the scattered field of the RIS.

B. Complex Radar Cross Section of the RIS
With known locations of the transmitter and the RIS, the angle of incidence and the polarization vector of the incident plane wave from the transmitter to the RIS can be determined. As a result, the scalar RCS, σ(θ, ϕ) ≡ σ(θ i , ϕ i ; θ, ϕ) of the RIS for the direction (θ i , ϕ i ) of the incident plane wave of the transmitter, can be found via full-wave simulation. To account for the phase of the scattered fields [i.e., ψ in (6)], we use the complex RCS [14], whose phase is defined as where ϕ sc and ϕ inc are the phase of the scattered electric field and the phase of the incident electric field at the center of the RIS, respectively. In (7), ϕ σ can be obtained by full-wave simulation and ϕ inc can be calculated by where ϕ t is the phase of the transmitted field, k = 2π/λ is the wavenumber in free space and R 1 is the center-to-center distance from the transmitter to the RIS, as shown in Fig. 1. Note that (8) is also applicable for NLOS incident rays. In that case, R 1 is the total unfolded path length. For a specific geometry, ϕ σ and ϕ inc can be computed via full-wave analysis and (8), respectively. Then, ϕ sc can be computed from (7). We should, finally, note that this methodology of treating the RIS as a transmitter with a pattern that corresponds to its complex RCS is similar to the technique that was used to account for diffuse scattering in ray-tracing of millimeter-wave communication channels in [16].
The remaining unknown term in (6) is P t ≡ P RIS,inc . Its computation is explained in the following.

C. Scattered Power From the RIS
The incident power at the RIS can be rigorously determined by integrating the power density of the transmitted electromagnetic field on a surface enclosing the RIS. To approximate it by raytracing, we place a mesh of M ×N isotropic receivers, uniformly dividing the RIS area in M ×N cells of area ΔS. We compute, by ray-tracing, the received electric field E i,j (i = 1, . . . M, j = 1, . . . N) at each receiver. In general, this field includes direct and multipath components from the transmitter. Then, the incident power can be approximated as This is an approximation that can be readily deduced by raytracing, at the cost of omitting polarization loss at the RIS.
To determine the number of cells and the area ΔS, we apply the following iterative procedure. 1) We set ΔS to an initial value ΔS 0 and compute the corresponding P 0 RIS,inc , with (9). 2) We increment the number of receivers M , N in each direction, decreasing the cell area ΔS. At the kth iteration, the area is ΔS k and the power computed by (9) is P k RIS,inc . 3) For each iteration, we compute |P k RIS,inc − P k−1 RIS,inc |/ P k−1 RIS,inc , tracking the convergence of the value of the computed power. When this ratio is below a threshold ε, we terminate the iteration and set P RIS,inc = P k RIS,inc .

D. Propagation Modeling
We compute the total signal strength across the communication channel, assuming a single, direct ray-path from the transmitter to the RIS. This assumption is compatible with the operation concept of RISs, whose role is to redirect a narrow beam from the transmitter toward intended receiver points, bypassing obstructions along the line-of-sight [17], [18]. The modeling steps are summarized as follows.
1) We use a standard Shooting-Bouncing Ray (SBR) raytracer, as in [16], to obtain the field E 1 that the transmitter would create in the environment, in the absence of the RIS. In this ray-tracing run, no reflections are allowed from the RIS area, i.e., any reflections from the area that is occupied by the RIS are removed. 2) We compute P RIS,inc , with ray-tracing, using the iterative approach explained above. 3) We compute the complex RCS of the RIS using full-wave analysis. In the results shown below, we use the Ansys HFSS finite-element solver. 4) The complex RCS allows us to model the scattered field from the RIS as a secondary transmitter. We embed this transmitter in the SBR ray-tracer and compute the electric field E 2 it produces throughout the domain. We add the phase ϕ sc from (7) to the phase of E 2 . Notably, this simulation does not include the original transmitter.

5) The total field is
This approach is appealing for its efficiency. However, it relies on the assumption that the interactions of the RIS with its environment can be approximated by ray-tracing its complex RCS. Similar to ray-tracing the far-field pattern of installed antennas in complex environments [12], this approach is approximate, as it disregards near-field interactions of the RIS. In practical indoor geometries that are much larger than typical RISs, this limitation is not severe. In the following, we compare our method to fullwave analysis of a small geometry, with reflecting/diffracting walls in close proximity to the RIS, demonstrating its efficiency and accuracy.

A. Free-Space RIS Channel
To validate the proposed method, the communication channel shown in Fig. 1 was first modeled without the wall and the obstacle using the proposed method and HFSS. This is a manageable problem for a finite-element solver, purely used for validation of our method.
The channel was simulated at 2.45 GHz. The transmitter, which was a rectangular horn antenna, was centered at (0 m, 2.5 m, 0.3 m). The RIS, which consisted of 24×24 square patches, was centered at (5 m, 0 m, 0.3 m). The RIS area was A ph = 386×386 mm 2 . We computed the field in the region of 5 m ≤ x ≤ 10 m and 2 m ≤ y ≤ 5 m on the z = 0.3 m plane, which was in the far-field of both the transmitter and the RIS. All simulations presented in this letter were performed on a 12-core processor desktop.
The realized gain of the horn antenna and the RCS magnitude of the RIS on the x-y plane are shown in Fig. 2. The obtained patterns were imported into our SBR ray-tracer. The total power of the horn antenna was set to 30 dBm. To compute the incident power on the RIS, P RIS,inc , the threshold ε and the initial dimensions of cells ΔS 0 were set to 0.01% and A ph /(20×20) (with M = N = 20), respectively. After five iterations (ΔS 1 = A ph /(40×40), ΔS 2 = A ph /(60×60), ΔS 3 = A ph /(80×80), ΔS 4 = A ph /(100 × 100)), P RIS,inc converged to 9.56 dBm. Fig. 3 compares the results obtained by the proposed method and HFSS, from which we can observe excellent agreement between them. Fig. 4(a)-(c) shows  |E 1 |, |E 2 |, and |E 1 + E 2 |, respectively, readily obtained by the proposed method in an extended region of 0 m ≤ x ≤ 10 m and 0 m ≤ y ≤ 5 m on the z = 0.3 m plane. We can see the formation of multiple side lobes emanating from the RIS and interference between transmitted and scattered fields over a wide zone that our method can capture. These effects are important for the overall fading profile of the channel and project the importance of the proposed analysis.
It is worth mentioning that the computation time for the full-wave analysis was 1 hr 14 min, while that of the hybrid method was only 15 min 25 sec including full-wave simulations (12 min 17 s) of the radiation pattern of the transmitter and the complex RCS of the RIS.

B. RIS Channel With Non-Line-of-Sight Points
We now compare our method to HFSS for the geometry of Fig. 1, which includes NLOS points, obstructed by an obstacle. The role of the RIS is to create strong propagation paths that bypass the obstacle, effectively serving these NLOS points.
The concrete wall (ε r =2.5, tanδ=0.15) was placed right behind the RIS and a copper slab was centered at (5 m, 3.5 m, 0.3 m) as an obstacle. The wall and the slab had dimensions of length × height × thickness = 4 × 0.6 × 0.02 m 3 and 2×0.3×0.02 m 3 , respectively. The total power of the horn antenna was set to 30 dBm and the total power of the RIS was also obtained as 9.56 dBm with the same set parameters, i.e.,  , ΔS 0 , M , and N , as the free-space modeling. The number of reflections and diffractions were set to 15 and 5, respectively.
Restricted by the large memory usage of the finite element analysis, we computed the electric field in a region 5 m≤ x ≤ 10 m and 2 m ≤ y ≤ 5 m on the z = 0.3 m plane, by simulating five subdomains separately. In this case study, the full-wave analysis took 11 hrs 18 min in total to complete, while the hybrid modeling took only 16 min 45 s (also including 12 min 17 s of full-wave simulation time). By comparing these computation times with those of the free-space case, we note that the full-wave execution time grows dramatically as the environment becomes more complex, while the hybrid technique remains an efficient alternative. Fig. 5(b) shows the magnitude of the electric field obtained by HFSS. We can observe some discontinuities at x = 6, 7, 8, 9 m, at the interface of the simulation subdomains. The magnitude of the electric field obtained by the proposed method is shown in Fig. 5(a). Fig. 6(b)-(g) compares the electric field computed by our method and HFSS, along the lines specified in Fig. 6(a). The far-field distance of the RIS is 2.43 m. Although the electric field presented in Fig. 6(b) and (d) is partially in the near-field (i.e., for (b), y <2.21 m, for (d), x < 6.38 m) of the RIS, no more significant discrepancies compared to other regions are observed. In general, the presented results demonstrate that the hybrid method and finite element analysis result in very similar predictions for the electric field strength across the geometry.

IV. CONCLUSION
We proposed a hybrid technique combining full-wave analysis with ray-tracing to model RIS-enabled communication channels. The proposed method combines the accuracy of full-wave analysis to model the interaction of the RIS with its environment, with the efficiency of ray-tracing that allows us to derive propagation models for realistic geometries in a computationally feasible manner. For validation purposes, we applied our method to a small scenario that was manageable by finite-element analysis, demonstrating excellent agreement between our results and HFSS.