Accelerated IE-GSTC Solver for Large-Scale Metasurface Field Scattering Problems Using Fast Multipole Method (FMM)

An accelerated integral equations (IE) field solver for determining scattered fields from electrically large electromagnetic metasurfaces using fast multipole method (FMM) is proposed and demonstrated in 2-D. In the proposed method, practical general metasurfaces are expressed using an equivalent zero thickness sheet model described using surface susceptibilities, and where the total fields around it satisfy the generalized sheet transition conditions (GSTCs). While the standard IE-GSTC offers fast field computation compared with other numerical methods, it is still computationally demanding when solving electrically large problems, with a large number of unknowns. Here, we accelerate the IE-GSTC method using the FMM technique by dividing the current elements on the metasurface into near- and far-groups, where either the rigorous or approximated Green’s function is used, respectively, to reduce the computation time without losing solution accuracy. Using numerical examples, the speed improvement of the FMM IE-GSTC method { $O(N^{3/2})$ } over the standard IE-GSTC { $O(N^{3})$ } method is confirmed. Finally, the usefulness of FMM IE-GSTC is demonstrated by applying it to solve electromagnetic propagation inside an electrically large radio environment with strategically placed metasurfaces to improve signal coverage in blind areas, where a standard IE-GSTC solver would require prohibitively large computational resources and long simulation times.


I. INTRODUCTION
E LECTROMAGNETIC metasurfaces are 2-D arrays of sub-wavelength resonators, which are engineered at the microscopic scale to achieve a desired macroscopic field response. By engineering the shapes, sizes, and structures of these resonators forming the unit cells, metasurfaces can enable a variety of wave transformations when excited with a specified input field, in reflection, transmission, or both. Consequently, they have found a myriad of applications throughout Manuscript received 16  the electromagnetic spectrum, ranging from radio waves to optical frequencies [1]- [7]. While the metasurfaces are constructed using subwavelength resonating particles, they themselves are electrically and physically large enabling a significant interaction with the impinging electromagnetic waves. Optical metasurfaces typically consist of millions of resonators (plasmonic-or dielectric-based), which are engineered at the nanoscale. On the other end of the spectrum, at the radio frequencies (RFs), due to standard and mature printed circuit board fabrication processes, unit cells of the order of ≈λ/10 are quite common, and the metasurfaces are several wavelengths long in typical proof-of-concept demonstrations.
However, several practical system-level applications require large area metasurface sizes beyond laboratory-scale experiments. Two examples are that of creating electromagnetic illusions and camouflage [8], [9] and using metasurfaces to engineer the wireless environment, for instance, [8], [10]- [12]. In both the cases, the metasurfaces are operated in larger environments with a variety of scattering objects of different shapes and sizes present. The full-wave field analysis thus becomes a computationally challenging task, with multiscale features ranging from sub-wavelength resonators of the metasurface to large physical objects placed around them.
To make the computational problem more tractable, its accuracy maybe sacrificed as a trade-off, as is commonly done in high-frequency wave propagations' regimes where ray-tracing [13]- [16] and paraxial-approximation-based methods [17] become more viable. Alternatively, to reduce the large number of computational unknowns and avoid dense structure meshing, physical metasurface structures can be represented by an equivalent zero thickness sheet models, where they are treated as spatial discontinuities. The metasurface field interaction is governed by the so-called generalized sheet transition conditions (GSTCs), in conjunction with the spatially varying electric and magnetic surface susceptibility,χ(r, ω) models of the metasurfaces [18], [19].
While a variety of numerical methods have been proposed to implement GSTCs into standard volumetric Maxwell's equation solvers; including finite difference, finite element, and spectral domain methods, for instance, [20]- [25]; they are still not suitable for electrically large simulations. On the other hand, integral equation (IE)-based methods are particularly attractive as large simulation domains can be handled, since no volumetric meshing is used and only surfaces are discretized [26]- [28]. Consequently, zero thickness sheet models implemented using GSTCs in IE solvers (IE-GTSC) have been shown to very efficiently reproduce the macroscopic field behavior of practical metasurfaces with large number of unit cells, when compared with standard brute-force full-wave computer simulations [8], [9], [29]- [31].
While IE-based numerical solvers are efficient compared with other techniques, they still require large computational resources as simulation domains become complex, increase in size, and have more scattering objects, as is typically expected in applications involving wave propagation in large-scale RF environments, for instance. To accelerate the IE-based solvers further, one attractive and well-known technique is to use the fast multipole method (FMM) [32]- [34]. FMM divides surface elements into localized groups and optimizes the far-group interactions to speed up the entire simulation, without losing accuracy.
In this work, in the context of solving large-scale electromagnetic problems, the IE-GSTC-based solver is coupled with FMM to accelerate the simulation speed and reduce the computational resource requirements. Specially, a 2-D FMM implemented via pulse expansion functions and point matching that can handle arbitrary polarization of electric and magnetic fields/currents is proposed here, which can solve large-scale GSTC/metasurface problems. While an accelerated IE solver has also been applied to large-scale environments before in the literature using a volume IE method [35]- [38], the use of volume integral makes this method unable to incorporate GSTCs. Although FMM is commonly applied to 2-D scattering problems [32], [34], much of what is reported in the literature is developed for either TM z or TE z polarization (in terms of EFIE, MFIE, or CFIE integrals), which are not easily extendable to general metasurface boundaries where both the electric and magnetic fields are tightly coupled through surface susceptibilities in GSTCs. Here, we fill this gap, by combining FMM with GSTCs to describe a general electromagnetic boundary.
The article is structured as follows. Section II starts by building an application case where the concept of using metasurfaces for engineering a wireless channel is used to motivate the requirement of an accelerated IE-GSTC solver. A brief summary of the IE-GSTC method is next described. Section III describes the FMM formulation and develops it further to be integrated with GSTCs. Section IV presents the numerical results by taking an example of a finite-sized uniform surface and compares its performance with standard IE-GSTC implemented using the boundary element method (BEM). It is followed by simulating large room RF environment to illustrate how metasurfaces can be used to increase signal coverage, which otherwise cannot be simulated using standard BEM solver.

A. Metasurface Application Case
Let us consider an illustrative example of wireless communication inside a closed space, as shown in Fig. 1(a). The space consists of doors, windows, and variety of physical objects, while an access point is assumed outside the rooms, aiming to provide signal coverage inside.
At microwave frequencies (≈<10 GHz), the electromagnetic (EM) wave interaction with the environment is highly diffractive, due to smaller antennas, compared with the physical environment and its operating wavelengths. EM waves bounce off physical objects in the wireless environment and propagate in uncontrolled directions, creating multiple paths and signal copies at the intended receivers. This provides rich channel diversity, which is used to recover useful signal buried inside the noisy non-LOS wireless channel. This is extensively used in typical multiple-input multiple-output (MIMO) systems to enhance the quality of signal transmission across devices [39]. As sources become more directive (typically due to increasing the frequency), waves propagate along prescribed directions, specularly reflect from walls while satisfying Snell's laws, and produce corner diffraction in the presence of sharp edges.
This directive wave propagation characteristic presents a unique opportunity to investigate a paradigm-shifting approach to engineer and design the RF environment itself and route the electromagnetic energy to form the desired and distinct wireless signal paths. Specifically, the RF environment can be modified by introducing controlled local reflection/ diffraction engineered EM objects in strategically chosen locations, instead of relying on passive reflection, and transmission through physical objects. This requirement of controlled reflection/diffraction is possible using engineered metasurfaces with specified wave transformations to improve signal coverage and throughput performance.
Consider an example scenario of Fig. 1(c) for directive communication, where the objective is to obtain the best possible signal at the intended receiver location from the source access point. With a conventional approach that uses the passive reflections off various walls and physical objects, such as in Fig. 1(b), receiving a weak or possibly null signal at the receiver is likely. However, if a network of smart metasurfaces are installed at strategic locations, the directive signal energy can be channeled to form a series of LOS signal paths between the access point and the receiver, as indicated by paths #1 and #2 in Fig. 1(c). Unlike the passive walls and reflectors of Fig. 1(a) and (b), engineered metasurfaces control the output reflection and refraction angles using generalized Snell's laws, which makes it possible to steer beams and avoid physical obstructions and undesired corner diffractions. In principle, several signal paths can be devised between the source and the receiver, from which the best path providing the best signal communication performance can eventually be chosen. In other situations, these metasurfaces can also be used to intentionally transform incoming directive radiation into omnidirectional wavefronts, as shown in Fig. 1(d). This opens two new capabilities.
1) Metasurfaces can channel EM energy in a hybrid directive and omni-directional manner, enabling uniform illumination of selective local regions in the environment. 2) They can locally enhance the channel diversity at mm-Wave frequencies, suitable for MIMO-type systems, which is otherwise generally not available at mm-Wave bands.

B. Standard IE-GSTC Formulation
The problem of modeling electromagnetic wave propagation in wireless radio environments described in Fig. 1 is naturally an electrically large problem which is aggravated by the presence of various metasurfaces making it a multiscale simulation problem. To simplify the computational setup, without losing accuracy, physical metasurface structures can be represented using analytical models. Given that all practical metasurfaces are electrically thin, i.e., δ λ 0 and with sub-wavelength periodicities, λ 0 , they can be represented using a zero thickness sheet model, which captures their effective dipolar electric and magnetic surface polarization responses. The complete electromagnetic properties of the sub-wavelength resonators can thus be described using compact surface susceptibility tensorsχ(r, ω).
The electromagnetic fields around the equivalent zero thickness sheet model are then governed by the GSTCs, which are given byn relating the field difference ({·}) across the metasurface with the average E-and H-fields in terms of its tangential {·} || and normal surface polarizations {·} n , respectively, i.e., If for simplicity we will limit the analysis to tangential terms, the GSTCs can thus be expressed aŝ Each of the surface susceptibilities,χ αβ with αβ ∈ {ee, mm, em, me}, is a 3 × 3 tensor with complex scalar elements, in general. Now, a more general statement of the problem of interest can be formulated, as also illustrated in Fig. 2(a) consisting of a metasurface surrounded by various objects of different shapes, sizes, and material. We wish to model scattering of an electromagnetic field off one or more surfaces (media interfaces, metasurfaces, object surfaces, etc.)-including the coupling between surfaces. Solving the scattering problem is essentially the need to find a self-consistent solution of the excitation, the propagation of the fields, and the interface conditions (e.g., GSTCs for metasurfaces) present for surfaces within the region of interest. If the surfaces present are uncharged, then these interface conditions relate the tangential components of the E and H fields across the interface or the surface. An appropriate formulation of Maxwell's equations to determine the field propagation is an IE approach.
The EM fields radiated into free space from given electric and magnetic current sources, {J, K}, can be generally expressed using an IE formulation as [26], [27] with r being the point of interest, r is the position of the source current; and E s and H s are the radiated (scattered) fields from the surface. 1 The field operators are given by with C ∈ {J, K}. G(r, r ) represents the Green's function which, for a 2-D case, is given by the Hankel function of the 2nd kind where J 0 and Y 0 are the Bessel functions of the 1st and 2nd kind and the function represents outwardly propagating radial waves.
The objective of the overall field scattering problem is thus to find a self-consistent field solution of (3) and (4), 2 for the specified excitation fields, thereby solving for unknown surface currents {J, K} on all surfaces along with the scattered fields in both the reflection and transmission regions. For example, the total field in the reflection region of the metasurface is a sum of the incident fields (ψ inc. ) and all the reflected scattered fields (ψ − s ) generated by unknown currents, {J, K}, on all the surfaces via (4). Similarly, the total fields in transmission are the total scattered fields (ψ + s ) generated by the same unknown currents. The relationship between the average fields (ψ inc. + ψ − s + ψ + s )/2 and their differences across the meta- 3 is given by GSTCs (3). This IE-based approach coupled with GSTCs is referred to here as the IE-GSTC field formulation. For numerical computation, these equations are discretized by segmenting the surfaces following BEM techniques and assembled in a matrix form for further processing [8], [9], [26], [27].

III. FAST MULTIPOLE METHOD (FMM) IE-GSTC
The IE part of the IE-GSTC formulation involves the matrix solution of a system-level propagation operator formed from the L{·} and R{·} operators to determine unknown surface currents. In large computational domains, with a large number of current unknowns direct solution of system matrix will become challenging due to very large memory requirements during solution and the computation time will become prohibitively large. An efficient solution to solve this problem and reduce the computation complexity is based on FMM [26]. FMM relies on approximating Green's function and modeling 3 Both electric and magnetic fields, which are coupled by GSTCs, unlike standard EFIE, MFIE, or CFIE formulations. the electromagnetic effect of a collection of surface currents as if they were replaced by a single element. Consequently, the surface elements are divided into near-and far-groups. While the near-group elements are solved using standard L{·} and R{·} operators, the far-group interactions are based on an approximated form of Green's functions. This formulation is very well-suited to the use of an iterative matrix solver as it very efficiently increases the speed of matrix vector multiplications. The use of a such a solver also removes the very large memory requirements of a direct solver [26].

A. Multipole Expansion of Green's Function Between Groups
For far-group computations (i.e., r ≥ r or ρ > ρ ), Green's function of (6), (2-D in our case) can be represented using a series expansion. This allows for the decoupling of the source (r ) and the observation terms (r). Using the addition theorem for the Hankel function of the 2nd kind, the 2-D electromagnetic Green's function can be written as [40] G To further prepare Green's function for far-group calculations, the Bessel functions J n (kr ) can be rewritten using their integral form as follows: where θ is an integration variable and k is a vector with a magnitude equal to the wavenumber, pointing in the θ direction, i.e., k = k{cos θx + sin θŷ}. Substituting (8) into (7) and rearranging yields the following expression for Green's function: (9) where the infinite series has been truncated so that it can be evaluated numerically. Next, two element groups can be formed-source and observation groups-using the decoupled form of Green's function, as illustrated in Fig. 2(b) and thus described in terms of group the geometry. For elements in separated groups, we can write Green's function in terms of the position of the elements within their respective groups and the vector pointing from the center of the source group to the center of the observation group. We begin by noting that [see Fig. 2 As the original form of Green's function is dependent only on the difference between the source and observation vectors, we can write using (10) Therefore, we can rewrite (9), by replacing r and φ with r ab and φ ab , and similarly r with (r s − r 0 ), leading to × e − jn(θ −φ ab + π 2 ) e j k·(r s −r 0 ) dθ (12) which now represents Green's function with respect to the two group centers.

B. FMM Propagation Operators
With the new form of Green's function, (12), for r > r , the two propagation operators, L{·} and R{·} of (4), can now be reformulated. For instance, for the L{·} operator, we can begin by noting for a 2-D case Using (12) (13a), as shown at the bottom of the next page. Similarly, for the R{·}, operator, we note that x .
Using (12) (13b), as shown at the bottom of the next page. These two operators of (13), due to their decoupled source and observation terms, and the form relative to their group centers, may be expressed in a more compact form as in terms of radiation function, T s (θ, r s ), which describes the radiation from the source element to the center of the source group, the transfer function, T ab (θ, r ab ), which describes the radiation between the source and observation centers, and the receive function, R 0 (θ, r 0 ) which describes receiving the radiation from the observation centers to the observation points. The only difference between the two operators is in the way the source element's radiation contributes to the source centers, via the D L and D R matrix operators. This now allows for a very intuitive and convenient interpretation of (14). Consider a single element of L{C} or R{C} where b represents a single source group, and s as all the sources in that group. w p is the weighting applied to the sample θ p , which is determined based on the chosen numerical integration method. In this case, trapezoidal integration was found to give good results when computing (15). All the source element contributions are first aggregated to the source centers via T s (θ, r s ). This forms an effective single source. This is next translated to the center of the observation group via T ab (θ, r ab ), and finally disaggregated to the desired observation points via R 0 (θ, r 0 ). Compared with the brute-force computation of (4), this three-step procedure is significantly more computationally efficient, as transfer, receive, and part of radiation function can be precomputed for all the values of θ and stored for repeated usage. This leads to compute time of the order of O(N 3/2 ) compared with O(N 3 ) when (4) is directly computed [26]. This methodology can now be easily applied to the zero thickness metasurface where source and observation elements are defined on the surfaces, as illustrated in Fig. 2(c). The number of groups N g needs to be defined and is generally chosen as N g = √ N , where N is the total number of meshed surface elements on the metasurface [32]. Field contributions next to the source group are modeled using rigorous L{C} or R{C} operators, while the rest of the observation elements follow (13).

A. Standard IE-GSTC versus FMM IE-GSTC
To validate the implemented FMM algorithm and compare its performance with a standard BEM solver, a simple simulation consisting of a finite-sized uniform metasurface of length is set up, as shown in Fig. 3(a). The surfaces are meshed at a density of ten divisions per wavelength for a total of N = 400 mesh elements. A Gaussian source of same size, , is used to generate a 60 GHz beam to illuminate the surface. The metasurface is designed to achieve a transmission coefficient T = 0.9 j and a reflection coefficient R = 0.1 when illuminated with a plane wave. The surface susceptibilities can then be synthesized using [18] With the surface parameters chosen, now the FMM simulation parameters-number of groups N g , number of terms taken in series expansion of Green's function q in (9), number of points, k, used to numerically integrate with respect to θ in (15)-must be selected so that the FMM algorithm runs at an optimal speed and converges toward the rigorous BEM result. As stated previously, the number of groups N g is typically chosen to be approximately equal to the square root of the number of mesh elements N. This maximizes the speed of the far-field matrix vector multiplication ensuring a compute time which varies as O(N 3/2 ) [32]. In our example, the mesh number of elements is 400, so we set N g = 20. Furthermore, the following empirical formula can be used to determine when to truncate the series expansion of Green's function [32]: where D is the diameter of the groups (D = √ 2). In general, q = 2N e , where N e is the number of elements per group, is enough to ensure convergence and is often smaller than what is given by (17). For this particular setup, we have N = 400 elements divided into N g = 20 groups, resulting in N e = 20 elements per group. Therefore, the series expansion of Green's function is truncated at q = 40. For the numerical integration with respect to θ , we see from (13) that if we were to represent this function using discrete values of θ , we would need q samples to satisfy the Nyquist sampling theorem. For our chosen value of q = 40, we therefore set k = 40.
Next, the iterative solver and its parameters are to be chosen. MATLAB has several iterative matrix solvers built in. In general, GMRES can be applied to most problems, but in some cases, other methods such as the stabilized bi-conjugate gradient method and conjugate gradient squared method may be able to solve the system more efficiently although their convergence is less predictable. For this example, the stabilized bi-conjugate gradient method was used as it was able to converge faster than GMRES [41]. The tolerance of the solver was set to 10 −6 . In addition, a set of preconditioner matrices can be chosen into improve the conditioning of the system and speed up the convergence of the iterative solver. For this example, ILU(0) preconditioning was used with a drop tolerance of 10 −6 .
With the initial FMM parameters set, the convergence of the accelerated FMM IE-GSTC solver can now be tested against a standard BEM-GSTC solver. To do this, the FMM parameters are set to their standard values given above, and we vary q, k and the number of groups, N g . For various cases, the electric current density, J(y), on the metasurface is calculated and compared with the standard BEM result. Fig. 3(b) shows the effect of q and k on J, when the group size N g is fixed. From these results, it is clear that the FMM results converge to the BEM-GSTC results when q = 40 and k = 40, as expected, but if either of these values are lowered, the FMM no longer produces accurate results. Next, the effect of group size on convergence can be examined, and it is shown in Fig. 3(c). We find that when the group size = 2/ √ N = 0.01 m (20 groups in total), the FMM converges to the BEM-GSTC result. However, when we increase the group size to 0.02 m, for instance, the result significantly deviates. This is because as the group size grows, the series expansion truncation, q, and by extension, the number of numerical integration points, k, have to increase accordingly. Since these values were held constant, the far-group representation of Green's function was no longer valid. Now that the FMM results have been verified, the scalability of FMM can be compared with BEM-GSTC. This is done by increasing the length 2 m of the two surfaces (source and the metasurfaces) in Fig. 3(a) and comparing the time it takes both the methods to solve for the current densities under converged conditions. These results are shown in Fig. 3(d). The results show that for smaller lengths, the BEM outperforms the FMM in terms of speed, but as the length of the surfaces is increased, the FMM solves the system of equations significantly faster with a difference of over two orders of magnitude, when the surfaces are 2 m long. This very closely follows the expected simulation time of order O(N 3 ) and O(N 3/2 ) for BEM-GSTC and FMM BEM-GSTC solvers, respectively. From these results, it is clear that to obtain full-wave solutions for metasurfaces placed in electrically large simulations, FMM is a very effective method to accelerate the standard BEM-GSTCs solver.

B. Complex Indoor Environment
Let us apply FMM BMM-GSTC to a practical scenario where the scattered fields need to be computed in a large computational domain, i.e., wireless communication. Fig. 4 shows an example of an arbitrarily chosen two-room scenario for illustration, where each of the room consists of doors and windows, which are connected to each other by a window in the partition wall. While the doors and windows are modeled as transparent surfaces, walls are made up of typical dielectric material (20 cm thick plasterboard with r = 2.02 and tan δ = 0, for example [42]). A fixed source is located along one of the walls in room #1, which is radiating a fixed beam at a chosen operation frequency (e.g., 6 GHz for microwave communication). Our first goal is obtain the field coverage inside this reference environment.
To make some simplifications to this overall large problem, and considering we are interested in fields inside the rooms only, we create an approximate equivalent model of the dielectric walls, by replacing them with a zero thickness surface susceptibility model in terms of χ ee and χ mm . This is achieved using the Fresnel transmission and reflection coefficients of a dielectric slab and using them in (16) to compute the two tangential susceptibilities. Windows and doors are assumed to be perfectly transparent and are modeled by setting all susceptibilities to zero. It should be noted that these equivalent surface susceptibilities are strictly valid for normal plane-wave incidence only, and a more sophisticated tensorial models must be used to properly account for angular scatterings.
Next, the total fields are computed using the FMM BEM-GSTC solver assuming a Gaussian beam generated by the source. Fig. 4 shows the field distribution in the entire computation region. The beam propagates with passive reflections off various walls, windows, and doors, eventually resulting in particularly a weak signal in room #2 and limited signal coverage overall. Let us next consider that we are particularly interested in enhancing this signal coverage in three specific regions (marked regions #1-#3 in Fig. 4). To achieve this goal, we will now design and strategically install metasurfaces in specific locations within the environment.
Among multiple potential solutions to design the wireless environment using metasurfaces, one example plan is shown in Fig. 5(a). We install three metasurfaces (marked MS#1, #2, and #3), one on each of the walls opposite to the source and one on the connecting window between the two rooms. The first metasurface #1 takes the incident source field and reflects a Gaussian beam of smaller width toward the second metasurface, thus acting as a good reflector with steerable beam following a generalized Snell's law. The second metasurface #2 is a reflective beam splitter that generates two Gaussian beams, one illuminating the upper left hand corner of the first room targeting region #1 and the other illuminating the metasurface on the adjoining wall. The third metasurface #3 is an active transmissive beam splitter with gain designed to generate two beams inside Room #2, targeting the desired regions #1 and #2.
Next, to follow this plan, metasurfaces are synthesized assuming the incident fields are Gaussian beams for simplicity and can be described by where x s is the position along the metasurface, r is the distance between the source of the beam and the surface, θ is the incidence angle, σ is the beamwidth at the surface, and σ 0 is Fig. 5. Possible plan to strategically place engineered metasurfaces to route the electromagnetic waves to desired locations, thereby significantly increasing signal coverage. Computed total fields abs{|E|} shown for the two cases of (a) no metasurfaces (base case) and (b) with three specific metasurfaces. The various metasurface design parameters are provided in Table I. MS#1 is 2.5 m long, whereas MS#2 and MS#3 are 2 m long. the initial beamwidth. The increase in the Gaussian beamwidth can be described by standard paraxial wave propagation as [17] σ (z) = σ 0 With these equations, an approximate expression can be found for the incident field on the metasurface based on either the Gaussian source or the specified scattered fields from the previous metasurface. The scattered fields on the metasurface are then specified to achieve the surface's desired functionality. With these approximate expressions for the incident, reflected, and transmitted fields known, the desired surface susceptibilities can be calculated using and are summarized in Table I. With the surface conditions set, the simulation parameters can now be determined. The FMM parameters q, N q , and k are set using the same criteria discussed in the previous section. Now the iterative solver parameters are chosen. Due to the size of this environment, convergence was an issue. For this reason, GMRES was the chosen iterative algorithm as it is the most reliable. As previously, ILU(0) preconditioning was used. This time, the drop tolerance was increased to 10 −4 to speed up the generation of the preconditioner. To further improve the convergence, MATLAB's built-in equilibrate function was used to permute and scale the system of equations to reduce the condition number of the system and lead to a more efficient iterative solution. Finally, the tolerance of the iterative solver was set to 10 −5 for acceptable convergence. The simulation was run on a server with 256 GB of memory. The time it took to solve the system of equations using the accelerated method was 369.5 s. The peak memory used by the simulation was 1.4 GB and the peak memory used to solve the system of equations was 197.4 MB.
The total fields are next computed inside the two rooms in the presence of the three metasurfaces, and the resulting field distribution is shown in Fig. 5(b). A clear and a significant improvement in the signal coverage is observed across the two rooms, in comparison to the reference case of Fig. 4, with the desired functionality of the three metasurfaces successfully achieved, as apparent from their respective scattered fields. Moreover, as desired and expected, the signal in the three targeted regions has been boosted, thanks to the strategic routing through the three metasurfaces. While this example illustrates leveraging the wave transformation capabilities of metasurfaces to control wave propagation inside complex environments, it clearly shows the importance and applicability of the proposed FMM IE-GSTC solver in computationally large problems that a standard IE-GSTC solver may not be able to handle.

V. CONCLUSION
An accelerated IE-GSTC solver for determining the scattered fields from an electrically large electromagnetic metasurface based on FMM has been proposed and numerically demonstrated in 2-D. The proposed method builds upon the standard IE-GSTC simulation framework, where a general metasurface is expressed using an equivalent zero thickness sheet model described using surface susceptibilities and the fields around it satisfy the GSTCs. The field propagation is implemented using an IE formulation. It has been shown that FMM accelerates the computation speed by dividing the unknown current elements on the metasurface into near-and far-groups, where either the rigorous or the approximated Green's function is used, respectively. Using numerical example, the speed improvement of the FMM IE-GSTC method {O(N 3/2 )} over the standard IE-GSTC {O(N 3 )} method has been confirmed. Finally, usefulness of the FMM IE-GSTC has been demonstrated by applying it to solve electromagnetic wave propagation inside an electrically large radio environment equipped with strategically placed metasurfaces to improve signal coverage in blind areas.
While the method described in this article is implemented for 2-D simulations, the development of the matrix operators D L and D R can be extended to 3-D with no issues. In addition, a 3-D full-wave simulation of a large indoor environment is very computationally demanding, reducing the problem to 2-D allowing for these simulations to occur over a reasonable time frame. This is especially important if multiple simulations are to be run in an iterative design procedure. While some accuracy is lost in a 2-D simulation, techniques are being developed to construct an approximate 3-D IE solution from a 2-D solution that shows good agreement between the approximate and exact 3-D solutions [37], [38]. The results generated by the proposed FMM algorithm could potentially be extended using the methods proposed in [37] and [38] to generate approximate simulation results of metasurfaces embedded in a large 3-D environment while maintaining the speed and memory benefits of a 2-D simulation.
Future efforts in this direction may involve subsequent improvements in FMM implementation based on various known techniques such as fast Fourier transforms (FFTs) and multilevel FMM, for instance, [43]- [45] and extension to 3-D. While a single application case of the radio environment has been used here for demonstration, the proposed FMM IE-GSTC is expected to be useful in a variety of other electrically large simulation problems to achieve an optimal use of available computational resources and reduce simulation times.