An Equivalence Principle-Based Hybrid Method for Propagation Modeling in Radio Environments With Reconfigurable Intelligent Surfaces

In this article, we present an equivalence principle-based hybrid method to model wave propagation in wireless communication channels with reconfigurable intelligent surfaces (RISs). Our method computes received signal strength (RSS) in RIS communication channels, considering all practical factors that influence them, such as mutual coupling of RIS cells, edge effects, and multipath propagation toward and from the RIS. We show that the accuracy of our method is comparable with that of full-wave analysis, both near and far from the RIS, through simple examples that are manageable by the finite-element method. Also, we experimentally validate it by comparing simulated and measured data for an indoor radio environment with an anomalous reflection metasurface.


I. INTRODUCTION
R ECONFIGURABLE intelligent surfaces (RISs) can dynamically redirect electromagnetic (EM) waves in a radio environment to selected directions, improving the received signal strength (RSS) in non-line-of-sight (NLoS) environments [1], [2], [3].Assisted by RISs, the performance of wireless communication networks can be improved in a controllable, programmable, and cost-effective manner [4], [5].Thus, RISs are considered a promising technology for current, emerging, and future generations of communication systems [6], [7].
In recent years, considerable research has been dedicated to RIS channel modeling.This has progressed from ideal singlelink scenarios [8], [9], [10], [11], [12], [13], [14], where only Fig. 1.RIS communication channels.(a) Ideal RIS channel, where the link from the Tx to Rx is fully obstructed and only the link from a Tx to Rx via an RIS is considered and (b) indoor RIS channel, where the multipath propagation influences the RIS performance.Selected ray paths are displayed, with green highlighting the edges of diffraction and yellow highlighting the planes that intersect with rays.
the link from a transmitter (Tx) to a receiver (Rx) via an RIS is considered [Fig.1(a)], to modeling actual 3-D radio environments [Fig.1(b)] [15], [16], [17], [18], [19], [20], [21].In real-world environments, the Tx radiates fields toward the RIS and other minor lobe directions.These fields interact with the environment, as illustrated by the red dotted rays in Fig. 1(b).The associated reflected, transmitted, and diffracted fields may impinge upon the RIS with nonnegligible amplitude compared with the direct wave, contributing to the scattered fields of the RIS.In fact, many wireless communication scenarios occur in highly reverberant environments, such as dense urban areas and indoor geometries [e.g., the office shown in Fig. 1(b)] [22].In these cases, rich scattering may cause the deviation of an RIS from its intended behavior [23].
This can be examined by considering multiple incident waves on the RIS that produce multipath effects, as shown by the blue dotted rays in Fig. 1(b).
This work considers the modeling of wave propagation in specific radio environments with RISs.We propose a hybrid ray-tracing/full-wave method based on the equivalence principle.First, we determine multipath incident fields at the RIS.Then, we derive the equivalent sources that correspond to the scattered field of the RIS, computed by full-wave analysis.Hence, mutual coupling between RIS cells and edge effects are fully accounted for.By incorporating these equivalent sources into ray-tracing and applying near-and far-field computational algorithms, our method accurately computes the electric field at Rx points both near and far from the RIS.We show that the accuracy of this method is comparable with the full-wave analysis.
The principle of this method has been briefly introduced in [24], in preliminary form.Here, we provide a comprehensive presentation, including detailed derivations, in-depth discussions, and experimental validation.This is the first method that accurately and efficiently models wave propagation in realistic and 3-D radio environments with RISs, including near and far fields of the RIS.
The rest of this article is organized as follows.Section II reviews existing methods for analyzing RIS channels, to establish the challenges addressed by our work.Section III introduces the modeling of an RIS via equivalent sources.Section IV explains how these are combined with ray tracing to derive a site-specific propagation model.Section V presents the two case studies to validate the proposed method by comparing the simulation results with HFSS finite-element analysis (FEA).Section VI numerically demonstrates the crucial role of multipath propagation in RIS channels.Experimental validation is presented in Section VII.Section VIII summarizes our contributions.

II. REVIEW OF RIS CHANNEL MODELING: FROM A SINGLE LINK TO MULTIPATH PROPAGATION
This section reviews existing methods for RIS channel modeling.We include methods that only consider a single link [Fig.1(a)] and those that consider multipath propagation in realistic radio environments [Fig.1(b)], to establish the research gap addressed by our work.

A. Modeling of Ideal RIS Channels
The analysis of an ideal single-link RIS channel requires the computation of field patterns of the RIS due to the incident wave from the Tx so that the path loss and link budget can be evaluated.The models in [8], [9], [10], and [11] treat RISs as arrays of independently scattering elements.They determined the scattered field of a single unit cell first and applied field superposition to compute the total scattered field of the RIS for multiple configurations, i.e., different combinations of the phase shifts of unit cells.However, unit cells are strongly coupled due to their subwavelength spacing [25].This leads to significant discrepancies between calculated and actual scattered electric fields [26].Consequently, these models have limited accuracy [17].This limitation is addressed by methods that capture mutual coupling effects [12], [17].
In [17], we introduced the concept of the effective complex radar cross section (CRCS) to compute the scattered fields of an RIS for multiple configurations.The mutual coupling was accounted for by the effective CRCS, which described the far-field scattering of a unit cell in an array.This method produced accurate scattered fields along the mainlobes of an RIS.In addition, Gradoni and Di Renzo [12] presented a circuit-based method, where RIS unit cells were modeled as perfectly conducting thin wires, whose mutual coupling effects were captured by the mutual impedance of two thin wires.
Moreover, Degli-Esposti et al. [13] and Vitucci et al. [14] considered other factors that contribute to the scattering properties of an RIS, such as edge diffractions.In addition to modeling single-link scenarios, attempts for modeling RIS channels with multiple Txs and Rxs, referred to as multipleinput multiple-output (MIMO) systems [27], [28], [29], have also been made.These models consider the links from all Txs to all Rxs via an RIS, while still disregarding multipath propagation.

B. Modeling of Actual RIS Channels
To fully understand how an RIS behaves in an actual radio environment, site-specific simulations that account for multipath propagation are needed.Full-wave methods [30] can model environmental effects.Yet, full-wave analysis of realistic 3-D geometries requires excessive computational resources.This motivated us to use ray tracing to model wave propagation of RIS channels in [15] and [16].We computed the CRCS of the RIS by full-wave analysis, thus determining its scattered fields.These fields were ray-traced to intended Rx points, in addition to the fields that existed in the environment without the RIS.In [17], this method was combined with the presented model for computing scattered electric fields from an RIS, to achieve efficient modeling of RIS channels for different RIS configurations.However, authors [15], [16], [17] did not include multipath components in the incident field at the RIS.Furthermore, the propagation modeling method cannot produce accurate results near the RIS, since it relies on the CRCS, a far-field quantity.
Subsequent research by other groups [18], [19], [20], [21] also adopted ray tracing for analyzing RISs in realistic environments.The method of [18] integrated the model proposed in [13] into a ray tracer to account for the contribution of the direct field from an RIS to the received power at an Rx.It relied on the reciprocity of radio links to account for multiple bounce rays from an RIS.Thus, the multiple bounce rays from the Tx, and the RIS could not be simultaneously considered.Also, Pyhtilae et al. [19] directly employed the RIS model proposed in [11], which used analytical expressions for conducting planes based on physical optics to calculate the scattered fields of unit cells.These analytical expressions are good approximations for electrically large planes [31], but the unit cells of an RIS have subwavelength dimensions.
Multipath propagation due to the far-field scattering of an RIS was also analyzed by ray tracing in [21] and [20].In [20], the RIS was modeled in an idealized way, assuming that all its input power is radiated in the mainlobe.Moreover, Huang et al. [21] applied geometric optics (GOs) to model the anomalous reflection of an RIS, by manually adjusting the normal vector of each unit cell.This changed the computed specular direction of a surface and captured the anomalous reflection of an RIS to some extent.However, this method broke the underlying assumption of GO that the plane of reflection should be electrically large, leading to considerable errors.The methods in [19], [20], and [21] also ignored cellto-cell coupling, which has been demonstrated to play an important role in the scattering properties of RISs [17], [25].
Compared with the methods discussed above for modeling practical RIS channels, our work: 1) uses full-wave analysis to compute the scattered fields of the RIS including mutual coupling and edge effects in a straightforward manner; 2) accounts for all multipath interactions toward and from the RIS; and 3) achieves accuracy comparable with the fullwave computations of field maps in RIS channel geometries, in regions near and far from the RIS, while at a fraction of the full-wave execution time.
The intended use of RISs is, indeed, as the last node of a communication link to an obstructed Rx via a single reflected beam.However, a modeling tool should not assume that this design goal is met by the RIS a priori.Including all multipath components toward and from the RIS, our model is a robust, channel-aware analysis and design tool for RIS geometries.

III. RIS MODELING IN RAY TRACING
VIA EQUIVALENT SOURCES Let us consider a scenario with a single Tx and a single RIS, modeled as a secondary source in an environment.The RIS is assumed to be deployed in the far-field region of the Tx, and all incident waves are approximated as plane waves.This is a practical assumption of RIS channels, where the physical dimensions of environments are considerably larger than the wavelength.
To trace rays from the Tx and the RIS, we can use the image method [32] or the shooting and bouncing rays (SBRs) method [33], [34].The image method determines exact ray paths but suffers from high complexity when modeling complicated geometries.The SBR method, in comparison, is more efficient but less accurate, because it computes approximate ray paths.However, the SBR method with an image-based ray path correction process combines the efficiency of SBR and the accuracy of image method [35].Hence, it is employed for ray tracing in this work.This ray path determination algorithm is briefly explained in Section IV-A.
To model multipath propagation, we launch rays from the Tx and the RIS and determine the exact ray paths that exist in the environment.The electric fields associated with the rays from the Tx are computed by the conventional field computational algorithm [32], with the known radiation pattern of the Tx.This section details the determination of the scattered fields from the RIS, assuming that the incident rays on the RIS and all the ray paths from the RIS to the Rx points of interest have been determined beforehand, by ray tracing.

A. Equivalent Surface Current Sources
We apply Huygens' principle [31] to find equivalent surface current sources that generate the scattered field of the RIS for incident fields (E inc , H inc ).These equivalent sources will be integrated with ray tracing to compute the scattered electric field of the RIS across the geometry.We select a bounding surface [shown by the dashed line in Fig. 2(a) and (b)] that is at an infinitesimal distance above the RIS, at z = 0.The electric and magnetic fields below the bounding surface at z < 0 are outside the region of interest, and are set to 0 [see Fig. 2(b)] [31].According to the equivalence principle, we introduce the equivalent electric and magnetic current densities on the bounding surface that satisfy the boundary conditions [36] where E tot and H tot are the total electric and magnetic fields on the surface and can be computed using the full-wave analysis with known incident waves.The total fields can be expressed as follows: where E inc and H inc are the incident electric and magnetic fields, and E sc and H sc are the scattered electric and magnetic fields.
We restrict the computation to the portion of the bounding surface that is over the RIS and coincides in size with the RIS [designated as surface S in Fig. 2(c)].Now, the RIS for the configuration of interest is replaced by the surface S with current densities J s and M s .In this work, we employ the Ansys HFSS finite-element solver with perfectly matched layer (PML) boundary conditions to simulate RISs and compute E tot and H tot , as shown in Fig. 2(d).
To integrate these equivalent sources on S with ray tracing, we need to compute the scattered electric field due to J s and M s .This topic is addressed next.

B. Near-and Far-Scattered Fields of the Equivalent Sources
By inspection of (1), both J s and M s have only xand y-components (in the local coordinate system of the RIS) Fig. 3 shows the surface S, with a differential element ds ′ .The x-, y-, and z-components of the scattered electric field due to J s at the observation point (x, y, z) are [31] where η 0 is the impedance of the free space, k 0 is the wavenumber in free space; x ′ , y ′ , and z ′ denote the coordinates of the differential element, R is the distance from each differential element to the observation point, and ds ′ = d x ′ dy ′ , and The three components of the electric field due to M s are Then, the total scattered electric field is In its far-field region, the scattered electric field of S in its local spherical coordinate system can be approximated as follows [31]: where x cos θ cos φ + J y cos θ sin φ e jk 0 r ′ cos ψ ds ′ (9a) The scattered electric field expressions ( 7) and ( 8) are utilized for the field computation in ray tracing, introduced in Section III-C.The fields in (7) and ( 8) are computed numerically, with discrete, equal-size differential area elements.Over each element, the fields are considered constant; they are computed at the center of the element by full-wave analysis of the RIS scattering problem [Fig.2(d)].Note that we determine the size, ds ′ , of the differential elements by studying the convergence of the computed scattered electric field, as we reduce ds ′ .
Large elements ds ′ can be used if the bounding surface is chosen at a larger distance from the RIS.In that case, the dimensions of the surface need to increase as well.This would necessitate larger dimensions for the surface S to accurately capture the current sources, thereby requiring a larger domain and increased computational resources for the full-wave analysis.

C. Embedding Equivalent Sources in Ray-Tracing
To compute the scattered electric field due to the equivalent sources in a given environment with ray tracing, the ray paths from the RIS to the Rx points are determined first.
Let us consider a typical indoor RIS channel as shown in Fig. 4, where the role of the RIS is to redirect the wave from the Tx to enhance the signal strength in room 2. As previously mentioned, the RIS is treated as a secondary source in the environment.Thus, the computation of ray paths from the RIS is done in the same way as for a primary Tx.We consider two Rxs of interest: Rx-1 and Rx-2, which are located in the far-and near-field regions of the RIS, respectively.The far-field distance of an RIS can be determined by comparing Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.the near-field pattern obtained by (7) with the far-field pattern computed using (8), and examining their difference.
With the ray path determination algorithm that will be explained in Section IV, rays that originate from the RIS toward the two Rx points are determined.For simplicity, we display a selection of determined rays in Fig. 4, namely, the direct ray r s 11 , the single-reflected ray r 12 , and the diffracted ray r 13 that reach Rx-1, as well as the direct ray r s 21 , the single-reflected ray r 22 , and the diffracted ray r 23 that reach Rx-2.The superscript s denotes that the corresponding ray is scattered from the RIS.The rays r 12 , r 13 , r 22 , and r 23 are, respectively, linked to the rays r s 12 , r s 13 , r s 22 , and r s 23 , as shown in the figure .The electric field at the Rx points is found by associating these fields with the scattered electric fields of the RIS in (7) and (8).
The far-zone scattered electric field at a distance d 0 , E s sc (r = d 0 , θ, φ), where (r, θ, φ) are the local spherical coordinates of the RIS, can be computed by (8).Then, the scattered electric field of the rays from the RIS at (r, θ, φ) can be calculated using [32] With (11), the electric field of a direct ray that reaches a far Rx point (e.g., ray r s 11 ) can be directly computed.Note that the term E s sc (r = d 0 , θ, φ)d 0 e jk 0 d 0 describes the far-zone scattering characteristics of the RIS, which is combined with the reflection coefficients, the transmission coefficients, and the uniform geometric theory of diffraction (UTD) coefficients [32], [37], [38], to compute the electric fields of nondirect rays.For example, the electric field of ray r 12 is [32] In (12), is a diagonal 2 × 2 matrix consisting of TE and TM reflection coefficients in the facet-fixed coordinate system of the intersection facet (where the reflected ray is generated), θ 12 and φ 12 are angles in the local spherical coordinate system of the RIS for the ray r s 12 , and d 12 is the unfolded ray length from the RIS to Rx-1 (i.e., the sum of the ray lengths of r s 12 and r 12 ).Note that E s sc is expressed as a 2 × 1 matrix that contains TE and TM components.
In the near-field region of the RIS, the electric field of the direct ray from the RIS to a Rx point (e.g., ray r s 21 ) is computed Fig. 5. Two-dimensional view of a realistic scenario that multiple incident waves may exist.The ray path correction process is also displayed: the single-reflected red ray is captured by the reception sphere, and its ray path is corrected based on the image method.
using (7).The fields associated with nondirect rays (e.g., rays r s 22 and r s 23 ) are computed with the same field computation algorithm as the Rx points in the far field.

IV. SITE-SPECIFIC PROPAGATION MODEL
FOR RIS CHANNELS In this section, we present our end-to-end modeling method.We first discuss how the ray tracing is used to determine the incident waves on an RIS, along with our ray path computation algorithm.Then, we summarize our methodology in four steps.

A. Incident Field at the RIS
To determine incident waves on an RIS, we place a virtual Rx at the RIS.Then, all incident waves (including their directions, polarization, and field levels) are determined by ray tracing the Tx, as shown in Fig. 5.Note that combinations of different types of rays (e.g., reflected-transmitted and reflected-diffracted) are accounted for.For simplicity, only four incident rays (the direct ray r i 1 , two single-reflected rays r i 2 and r i 3 , and a diffracted ray r i 4 ) at the RIS are displayed in Fig. 5.All incident waves determined by ray tracing can be included in the full-wave simulation round as wave excitations to compute E tot and H tot .
Fig. 5 also shows the process of computing ray paths.We apply the widely used (nearly) uniform ray launching method by tessellating an icosahedron [34].Depending on the angular separation of the launched rays and the ray lengths, reception spheres are used to capture rays that contribute to the total field at an Rx point.Then, the image method is applied to correct the ray paths.For ray r i 3 , the red and green dashed ray paths come from two adjacent launched rays (i.e., the red and green dotted rays), while only the red one is captured by the reception sphere.Then, the image source of the Tx with respect to the intersection facet is determined, and exact ray paths and their lengths from the Tx to Rx are obtained (the black dotted and solid ones).Note that the angular separation of the launched rays can affect the accuracy of the results.A large angular separation can result in some ray paths being lost during the SBR process, while a smaller angular separation can lead to increased computational cost due to more launched rays.Therefore, a convergence test is applied to determine the optimal number of launched rays.
When a complex environmental results in an excessive number of incident waves, it is possible to reduce the number of incident waves considered in the computation of the equivalent sources, by omitting rays with negligible amplitudes.For example, let the ray tracing of the Tx return a total of Q incident rays.The qth ray is included in the full-wave analysis of RIS scattering as an incident wave if its magnitude at the RIS, where e is a set threshold.

B. Summary
The modeling steps are summarized as follows.1) We first use ray tracing to determine the incident waves on the RIS.A virtual Rx point is placed at the RIS to compute all rays from the Tx that reach the RIS.The number of rays can be optionally reduced by applying (13).2) We set the determined incident waves as excitations in the full-wave simulation, as shown in Fig. 2(d).E tot and H tot of the RIS operating for a specific configuration, illuminated by the determined waves, are obtained.3) We ray trace the Tx.The rays that impinge upon the RIS do not generate any reflected rays.Instead, new rays are launched from the RIS and traced.This step results in M rays from the Tx and N rays from the RIS that reach an Rx. 4) At the Rx, we compute the fields associated with rays from the Tx with the conventional field computation algorithm [32].The fields of the rays from the RIS are computed by the field computation algorithm presented in Section III-C.Finally, we compute the total electric field at the Rx point as follows: where E Tx m and E RIS n are the electric fields of the mth ray from the Tx and the nth ray from the RIS, respectively.Note that the fields are expressed and added in the global coordinate system.This process returns the electric field at one Rx point.Steps 3 and 4 can be repeated to compute the electric field at multiple Rx points.Our methodology is applicable to RISs configured for diverse functionalities, such as anomalous reflectors, focusing lenses [13], and reflective polarization Fig. 6.Configuration of the used unit cell.ON/OFF states of the switch alter the current distribution, resulting in phase variation of the reflected field.
converters [39].The scattered electric fields obtained from J s and M s capture the intended functionality of an RIS.

V. COMPARISON WITH FINITE-ELEMENT ANALYSIS
In this section, we validate the proposed method by comparing its results with HFSS FEA.Restricted by the intensive computational cost of FEA, we studied two simple geometries that were manageable in HFSS and contained NLoS points.We investigated one scenario where the direct wave from the Tx to the RIS dominated the incident field on the RIS and another scenario, where the direct field impinged on the RIS along with multipath components of comparable strength.
We used the unit cell structure (Fig. 6) proposed in [40] to construct a 16 × 16 RIS with dimensions of 304 × 304 mm 2 (5.87λ × 5.87λ, where λ is the wavelength) that operated at 5.8 GHz.The unit cell has a diode switch to realize 180 • phase variation of the reflected fields.Thus, it features switch-ON and switch-OFF states.We used the full-wave modeling approach of [40], which applied short and open terminations in place of the diode switches for simplicity.This approach demonstrated good agreement between simulation and measurement results in [40].For more sophisticated modeling of the RIS, equivalent circuit models for both the switch-ON and switch-OFFstates of the diodes can be considered.
All the data presented in this section were collected at 5.8 GHz.In each case study, we applied near-field computation in the region r < 2 m because the near-field pattern at r = 2 m obtained by (7) matched well with the far-field pattern computed using (8).The value of the parameter d 0 in (11) for field computation was set to 10 m.With the test of the convergence of the modeling results, the dimensions of the differential elements were selected as 4 × 4 mm (λ/12.9× λ/12.9).Also, the number of launched rays at both the Tx and the RIS was 26 012 (achieving an angular separation between successive rays of approximately 1.52 • ).Only one diffraction was considered per ray, assuming that multiple diffracted fields are negligible [41].

A. Case Study 1: Direct Wave Dominates the Incident Field
We first studied a simple RIS channel, as shown in Fig. 7.A horizontally polarized rectangular horn antenna with a maximum gain of 18.04 dBi and transmit power of 30 dBm directly illuminated an RIS.The Tx was centered at (0, 3, 0.2 m), Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.which was also selected as the phase reference point (i.e., the phase of the radiated field of the Tx at this point was set to 0, for the determination of fields everywhere else).The RIS was mounted on the center of a wall with a relative permittivity ε r = 2.5, a loss tangent tan δ = 0.15, and dimensions of 2 × 0.4 × 0.01 m 3 .A copper slab with conductivity σ = 5.8 × 10 7 S/m and dimensions 1.5 × 0.2 × 0.01 m 3 was placed at (3, 3.25, 0.2 m) as an obstacle.The RIS was set with a configuration (see Fig. 7) that could redirect an incident wave with an angle of (θ = 45 • and φ = 180 • ) to an angle of (θ = 22 • and φ = 0 • ), bypassing the copper slab.All angles are with respect to the local coordinate system of the RIS in Fig. 6.Due to the high memory requirements of FEA, the simulation domain was restricted to the region 2 m ≤ x ≤ 5 m and 0 m ≤ y ≤ 4 m on the z = 0.2 m plane.We uniformly placed 30 351 Rx points in this region for the ray-tracing simulation.Due to the simplicity of the scenario, both the maximum numbers of reflections and transmissions allowed per ray in ray tracing were set to 2 (in fact, only single-reflected and singletransmitted rays could exist in this geometry).Three incident rays were determined by ray tracing as displayed in Fig. 7: the direct ray from the Tx with a field level of 14.64 V/m, and two diffracted rays with field levels of 9.22 × 10 −5 and 2.60 × 10 −5 V/m, respectively.These three rays were implemented as excitations in the full-wave simulation round.The interface between the near-and far-field computation domains (r = 2 m) is displayed in Fig. 7.
The electric field distribution obtained by the proposed method, as shown in Fig. 7, illustrated that the RIS redirected the incident wave from the Tx to the area blocked by the copper slab.This resulted in enhanced electric field levels at NLoS points (i.e., those in the shadow region of the copper slab).Fig. 8 shows the comparison of the electric fields along several lines within the geometry, obtained by HFSS FEA, the proposed method, and the method presented in [15].Note that only the direct wave was implemented as an excitation to model this scenario with the method of [15].The proposed method and the method of [15] produced very similar electric field levels in the far-field region of the RIS, as shown in Fig. 8(b) and (d).Since the direct wave from the Tx was the dominant component of the incident field, accurate far-field results were obtained by the method of [15].Fig. 8(a) and (c) shows that the proposed method provided significantly closer results to the FEA in the near-field region of the RIS compared to the method of [15].

B. Case Study 2: Enhanced Multipath Propagation Scenario
The second scenario, used for validation of our method, is shown in Fig. 9.The Tx was a finite-length dipole, with a far electric field [42] where I 0 = 1 A and l = 10 mm.The electric field is expressed in the local spherical coordinate system of the dipole (assuming that the dipo le is along the z-axis), which is different from the global coordinate system defined in Fig. 9.
The dipole was placed parallel to the x y plane at (0, 3, 0.2 m).Two copper slabs with dimensions of 1.5 × 0.2 × 0.01 and 1.2 × 0.2 × 0.01 m 3 were placed at (3, 3.25, 0.2 m) and (0.6, 4, 0.2 m), respectively.We used the same RIS with the same location and orientation as in the geometry of Fig. 7.A different RIS configuration was applied, as shown in Fig. 9.This configuration could redirect an incident wave with an angle of (θ = 45 • and φ = 180 • ) to an angle of (θ = 29 • and φ = 0 • ), in the local coordinate system of the RIS.This simple scenario was designed to examine the operation of an RIS in an enhanced multipath environment (as opposed to the first scenario where the direct ray from the Tx to the RIS was dominant).
Both the maximum number of reflections and transmissions were set to 3. The first ray-tracing round returned ten incident rays, namely the direct ray r 1 , the single-reflected ray r 2 , and the single-diffracted rays r 3 , r 4 , r 5 , r 6 , r 7 , r 8 , r 9 , and r 10 , as displayed in Fig. 9. Particularly, r 7 and r 8 were reflected off copper slab 2 first and were, respectively, diffracted by different edges of copper slab 1; r 5 and r 7 , r 6 and r 8 , and r 9 and r 10 , respectively, had the same ray paths.The electric field   I.All these rays were implemented as excitations in the full-wave simulation round.The finite-element simulation domain and the number of Rx points in this region (for the simulation using the proposed method) were the same as in the case study 1.
Then, we obtained the electric fields across the geometry, as shown in Fig. 10.Fig. 11 shows the comparison of the electric fields along several lines within the geometry, obtained by HFSS FEA, the proposed method, and the method of [15].Only the direct wave was implemented as an excitation to model this scenario with the method of [15].Fig. 11 shows that the proposed method achieved very close field level predictions to the HFSS FEA, as demonstrated by the comparison of the electric field levels along lines in Fig. 11(a)-(d).The method of [15], however, failed to produce accurate results across the geometry as depicted in Fig. 11(a)-(d) because of the multipath incident waves at the RIS.

C. Execution Time Analysis
The simulations were conducted on a desktop computer with 12 cores.Table II lists the comparison of the time required to execute the simulations.Both the proposed method and the method of [15] required a similar amount of time for the full-wave simulation round in each case study, whether for obtaining E tot and H tot or the CRCS of the RIS.The ray-tracing time of the proposed method in the two case studies was 10 and 11 s more than that of the method in [15].This was primarily due to the modified field computation algorithm in the near-field region of the RIS, where the electric field at each Rx point needed to be computed with (7) that involved numerical integration.However, the method of [15] applied far-field computation across the whole region of interest based on the CRCS of the RIS.Its computation of the direct field from the RIS in the near-field region was similar to using (8), which was faster than applying (7).Note that the computational expense of our method is primarily governed by the full-wave simulation part.The number of incident waves  has little impact on the full-wave simulation time.Yet, an RIS with larger electrical dimensions requires a larger full-wave simulation domain, resulting in increased simulation memory and execution time.
The proposed method is orders of magnitude faster than fullwave analysis.Compared with the method of [15], it exhibits higher accuracy in the near-field region of the RIS while consuming almost the same amount of time.Furthermore, it can be easily employed to model scenarios where an RIS is illuminated by multiple waves.

VI. ROLE OF MULTIPATH PROPAGATION
To demonstrate how multipath propagation can impact channel characteristics, we simulated an L-shaped hallway RIS channel.The channel geometry is displayed in Fig. 12, along with the electric field distribution obtained by the proposed method (by simulating 11 191 uniformly placed Rxs).The width and height of the hallway were 3.5 and 3 m, respectively.All planes were set as concrete with ε r = 5.31 and tan δ = 0.079.A Tx with a gain of 20.42 dBi and transmit power of 30 dBm was located at (−2.5, 3, 1.5 m).The same RIS as in Sections V-A and V-B was placed at (0.5, 0, 1.5 m).It was set to operate with a configuration (see Fig. 12, referred to as configuration A) that could redirect an incident wave with an angle of (θ = 45 • and φ = 180 • ) to an angle of (θ = 3 • and φ = 0 • ).The angles are with respect to the local coordinate system of the RIS in Fig. 6.
Transmission through the walls was omitted.The maximum number of reflections was set to 5 for the Tx.One diffraction per ray was considered, as in previous case studies.The 88 incident rays were obtained by the first round of ray tracing.In this case, we applied (13) to exclude rays of negligible magnitude, setting e = 30 dB.The direct ray from the Tx to the RIS had a field intensity of 19.16 V/m (25.65 dBV/m), while the second strongest incident ray had a field intensity of only 0.56 V/m (−5.04 dBV/m).Therefore, only the direct ray was considered as the incident ray at the RIS to conduct the full-wave analysis.Other simulation parameters were the same as in the previous case studies.
The electric field distributions across the NLoS hallway shown in Fig. 12 demonstrate the significantly improved wireless coverage due to the deployed RIS.In the region of (0 m ≤ x ≤ 3.5 m, 4 m ≤ y ≤ 30 m, and z = 1.5 m), the mean electric field strength increased from −15.95 to 0.21 dBV/m.
To reveal the role of multipath propagation, we set the number of reflections for the rays launched from the RIS to 0 (i.e., disregarding multipath propagation of the RIS) and 5 using the proposed method.The results are displayed in Fig. 13.There are substantial discrepancies between the two curves.These discrepancies show how pronounced multipath effects can develop in RIS channels.
To investigate dynamic wave redirection by the RIS, we modeled the same channel geometry as shown in Fig. 12 for two different configurations.All simulation parameters for propagation modeling remained the same.Configuration B was designed to redirect an incident wave from (θ = 45 • and φ = 180 • ) toward (θ = 14 • and φ = 0 • ), while the RIS worked as a specular reflector for configuration C (all switches off).The results, displayed in Fig. 14, distinctly showcase the impact of the RIS configuration on the electric field distribution within the NLoS hallway.Notably, under configuration C, despite the absence of direct main lobe illumination within the NLoS area, the reflected waves from the main lobe and the presence of sidelobes effectively enhanced wireless coverage.

VII. EXPERIMENTAL VALIDATION
In this section, we present modeling results of an actual indoor channel with an anomalous reflection metasurface, along with a measurement campaign as a means of validation of our method.Note that the metasurface was a non-reconfigurable surface (without switchable components); it was used to mimic an RIS operating with a specific configuration.

A. Channel Setup
In [43], we presented measurement data in an indoor radio environment with an anomalous reflection metasurface to validate the method presented in [15] and [16].The measurement data herein were collected using the same measurement setup as that in [43].The geometry was a hallway junction located on the eighth floor of the Bahen building on the University of Toronto campus, as shown in Fig. 15.A global coordinate system is also defined and displayed in the same figure.
The metasurface consisted of 45 × 33 nonuniform unit cells with metallic rings printed on a full ground plane-backed substrate, enabling anomalous reflection [44], [45], as shown in Fig. 16.The coordinate system is the same as that shown Fig. 15.The substrate was Rogers RO3010 with ε r = 10.2, tan δ = 0.0035, and a thickness of 1.28 mm.The unit cells had dimensions of 12.93 × 12.93 mm 2 (λ/4 × λ/4).The dimensions of the metallic rings were identical along the z-axis, but they varied along the x-axis, as listed in Table III (the unit cell numbers, from left to right, are indicated in Fig. 16).The metasurface was designed to redirect a vertically polarized (i.e., z-polarized) incident wave, with an incident angle of 26.6 • to the direction of 53 • at 5.8 GHz.Its dimensions were 609.6 × 457.2 mm 2 (11.8λ × 8.84λ).The role of the metasurface was to mimic an RIS operating with a specific configuration and to enable coverage within the shadow region of the hallway corner.
The Tx, the metasurface, and the Rx were at the same height of 1.34 m.The Tx was a vertically polarized rectangular horn antenna with a gain of 17.1 dBi and was connected to a signal generator with transmit power of 25 dBm.It was positioned at (−2.5, − 6, 1.34 m) and directly illuminated the metasurface, which was centered at (0.5, 0, 1.34 m).The Rx was a vertically polarized omnidirectional whip antenna with a gain of 2 dBi mounted on a tripod and connected to a spectrum analyzer.All facets including the walls, the floor, and the ceiling, were set as concrete with ε r = 5.31 and tan δ = 0.079 in ray-tracing simulations.The RSS at multiple grid points (white points in Fig. 15) was measured.In this measurement, we focused on the far-field performance of the metasurface.
Due to multipath propagation, small variations of the Rx position can cause large variations of the signal strength [46].In addition, small positioning errors of the devices are inevitable.Therefore, measuring multiple RSS values around each point and taking the average to reduce the impact of Fig. 15.Realistic indoor radio environment with an anomalous reflection metasurface used for validation of our method.The electric field obtained by the proposed method and the measured receiving points are also displayed.small-scale fading is common practice [46].We took 12 RSS measurements at each point by rotating the Rx about the axis of the tripod.The radius of rotation was 2 cm.The mean and the standard derivation of the RSS at each point were used to compare measured and simulated data, taking into account the local fluctuation of RSS due to small-scale fading (represented by the standard deviation).We first calibrated the measurement setup by measuring the RSS versus distance along the mainlobe direction of the Tx in the hallway, to evaluate polarization, impedance mismatch, and other losses.

B. Numerical Modeling Results and Experimental Validation
Same as the simulations in Section VI, the maximum numbers of reflections, transmissions, and diffractions were 5, 0, and 1, respectively.The 222 incident rays were obtained by the first round of ray-tracing, and we applied (13) to exclude rays of negligible magnitude, setting e = 30 dB.The strongest component of the incident field corresponded to the direct ray from the Tx, which was 4.65 V/m (13.35 dBV/m).The second highest field level among all the field components was only 0.11 V/m (−19.17 dBV/m), contributed by the black dashed diffracted ray in Fig. 15.Based on (13), only one wave corresponding to the direct ray was included as an excitation in the full-wave simulation round.Other parameters, including Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.d 0 , the number of launched rays at the Tx and the RIS, and the dimensions of the differential elements, were the same as in the previous case studies.The near-field computation was applied in the region r < 2 m.The simulations were performed on the same desktop as the previous case studies.We placed 22 781 uniform Rx in the simulation domain to obtain the electric field distribution displayed in Fig. 15.The time required to execute the proposed method was 78 m, 49 s (including 75 m, 14 s of full-wave simulation time).Fig. 15 shows the electric field distribution across the geometry obtained by our method, which demonstrates that the metasurface was successful in redirecting the incident wave from the Tx to the NLoS locations.Fig. 17 show the comparison of the average RSS obtained by the proposed method and measurements along three lines.For each point, we compute the mean value and the standard deviation of the 12 measured RSS values, RSS and σ RSS , respectively.In the plot, we show the mean value, RSS, along with a bar that ranges from (RSS − σ RSS ) to (RSS + σ RSS ).The simulated data generally fall within the bars, showing good agreement between measurement and simulation results, when the small-scale fluctuation of measured RSS is taken into account.The remaining discrepancies are attributed to the simplified input geometry model of the ray tracer.As shown in Fig. 15, the model included only the walls, omitting complex features such as doors, pillars, and windows.Based on the agreement between the measured and simulated data we obtained, we concluded that this simplified geometry model was sufficient for our validation study.

VIII. CONCLUSION
We presented a new hybrid ray-tracing/full-wave method based on the equivalence principle for rigorously modeling wave propagation in realistic RIS channels.The use of equivalent sources extracted from full-wave analysis ensures that all practical factors that influence the RIS behavior, such as mutual coupling between unit cells and the finiteness of the surface are captured.Our method accurately models electric fields in the near-and far-field regions of an RIS in actual environments, readily accounting for multiple incident waves.The accuracy of our method has been further verified by an indoor channel measurement.Our approach is applicable to reflecting RISs for various applications, such as wave redirection, lensing, and polarization conversion.Moreover, it can be extended to transmitting RISs, by modifying the surface, where J s and M s are computed.Hence, this method is a valuable tool for propagation studies in real-world RIS channels, as well as the site-specific analysis, design, and optimization of RIS structures.

Fig. 2 .
Fig. 2. Modeling problem.(a) Original problem, (b) equivalent problem, (c) problem to be solved, and (d) two-dimensional view of the full-wave setup.

Fig. 4 .
Fig.4.Two-dimensional view of a typical indoor RIS channel with some rays displayed.Rx-1 and Rx-2 are in the far-field and near-field regions of the RIS, respectively.

Fig. 7 .Fig. 8 .
Fig. 7. Simple RIS channel used to validate the proposed method.The electric field obtained by the proposed method, the FEA simulation domain, the sampling points (indicated by red dotted lines) for comparison, the configuration of the RIS (yellow denotes switch-ON and blue denotes switch-OFF), and three incident rays are displayed.

Fig. 9 .
Fig. 9. RIS channel with multiple incident waves impinging upon the RIS used to validate the proposed method.The incident ray diagrams and the configuration of the RIS are displayed.

Fig. 10 .Fig. 11 .
Fig. 10.Electric field distribution obtained by the proposed method.The FEA simulation domain and the sampling points (indicated by red dotted lines) for comparison are displayed.

Fig. 12 .
Fig. 12. Electric field distributions across the geometry for configuration A. Left side: electric field of the Tx only and right side: total electric field.

Fig. 13 .
Fig. 13.Electric field levels in the mainlobe direction of the RIS.

Fig. 14 .
Fig. 14.Electric field distribution across the geometry for (a) configuration B and (b) configuration C.

TABLE I ELECTRIC
FIELD LEVELS OF SELECTED INCIDENT RAYS FOR CASE STUDY 2 (UNIT: V/m)

TABLE II TIME
a REQUIRED TO EXECUTE THE SIMULATIONS levels of these are listed in Table