Accurate and Fast Analysis of Reflective Surfaces and Metasurface Antennas with Sheet Impedance Boundary Conditions

IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 1 Accurate and Fast Analysis of Reflective Surfaces and Metasurface Antennas with Sheet Impedance Boundary Conditions Jean Cavillot, Member, IEEE, Modeste Bodehou, Member, IEEE, Adam Abazi, Student Member, IEEE, and Christophe Craeye, Senior Member, IEEE Abstract—Simulating the fine geometry of Metasurfaces (MTS) structures in a conventional way is a difficult task requiring huge computational resources. On the other hand, the metallization can usually be modeled as an impedance sheet with modulation scale of the order of the operating wavelength. Even then, direct solution of the system of equations is usually not possible due to memory saturation. As a consequence, one has to resort to iterative methods. However, the analysis of an impedance sheet lying on a grounded slab based on an iterative solution of the Method of Moments (MoM) may lead to ill-conditioning and a large number of iterations. This is certainly the case when the range of impedance spans both the capacitive and inductive domains. Such impedances range is in practice required for Reflective Intelligent Surfaces (RIS), and for some MTS antennas. This paper proposes a preconditioner aiming to solve bad convergence issues caused by the wide range of the surface impedance. The preconditioner involves a multiplication by the conjugate of the MoM matrix followed by a block diagonal preconditioner. The block diagonal preconditioner has memory and multiplication complexity N3/2. It is precalculated in an accelerated scheme relying on FFTs. Besides, multiplication with the MoM matrix is carried out with complexity Nlog(N)thanks to the use of FFTs. Index Terms—Reflective Intelligent Surfaces (RISs), metasurface antennas, Method of Moments (MoM), impedance boundary condition (IBC), fast Fourier transform (FFT), Toeplitz matrices. I. INTRODUCTION REflective Intelligent Surfaces (RISs) [1] and metasurface (MTS) antennas [2] are currently at the center of attention of the antenna and propagation community. RISs are often described as a key factor for the success of the sixth generation of wireless communications due to their ability to engineer the propagation channel through manipulation of reflected fields [3]. Low-profile MTS antennas on the other hand are designed to generate radiating fields leaking from a surface wave (SW) [2], [4]. The feeder being incorporated in the MTS substrate, MTS antennas are particularly suited This work was supported in part by the Belgium Fonds National de la Recherche Scientifique (F.R.S.-FNRS); and in part by the HORIZON EUROPE EIC Pathfinder Open programme under grant agreement No. 101098996 ”Flexible intelligent near-field sensing skins (FITNESS)”. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them. (Corresponding author [email protected]) Jean Cavillot, Adam Abazi and Christophe Craeye are with University of UCLouvain. Modeste Bodehou is with with University d’Abomey-Calavi. for applications requiring high gain with low-profile and lowcost platforms. Both technologies have recently experienced a rapid development. For example, anomalous reflection has been demonstrated in [5], and multibeam reflection has been shown in [6], [7] using RISs. MTS antennas design for pencil beam [8], multi-beaming [9], [10], and shaped beams [4] have also been reported in the literature. All these meta-structures share a feature: the MTS and the RIS are usually realized with a dense texture of sub-wavelength metallization printed on a grounded slab. The metallization layer implements a smooth impedance sheet distribution over the surface. While being numerically easier to analyze as compared to the simulation of the physical structure of the MTS, several papers have shown that the sheet impedance model of the printed metallization provides an accurate prediction of the fields because it properly accounts for the spatial dispersion of the grounded slab [11], [12]. Despite the drastic reduction of the number of unknowns in the sheet impedance formulation, direct solving by Gaussian elimination becomes rapidly impractical due to memory saturation when the size of the MTS exceeds a few wavelengths. Therefore, fast and memory-efficient simulation methods have been developed using the sheet impedance model in an electric field integral equation (EFIE) formalism [13]–[16] based on the Method of Moments (MoM) [17]. Although the resulting system of equations matrix is particularly well-conditioned for strongly capacitive sheet impedances (as is usually the case for large aperture MTS antennas), conditioning issues rapidly appear when dealing with an impedance distribution spanning from capacitive to inductive (or even low reactance capacitive) [12], [16]. Such impedance ranges are in practice useful for RISs [11], [18]. Low reactance capacitive impedances can be adopted for improving the feeder efficiency of MTS antennas [19]. In this paper, we designate by ”low reactance capacitive impedance”, a reactive impedance with negative imaginary part which does not exceed a few hundreds of ohms. The conditioning problems associated with such impedances prevent the use of iterative solvers because the convergence becomes very slow. It is well known that the conditioning of the EFIE MoM system of equations matrix deteriorates when approaching low frequency or with growing discretization density [20]. Over the past few decades, several approaches have been developed to solve this issue. Considering antenna arrays of identical elements, enhanced versions of the block diagonal preconditioner can be used to improve the condition number This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 2 of the MoM matrix [21], [22]. Conditioning can also be improved using the generalized shielded block preconditioner, by dividing the whole structure into overlapping subdomains [23], [24]. These preconditioners are however very demanding in terms of memory allocation and are not fit for dense and connected structures in layered media [25]. Another class of preconditioners have been developed using the Calderon identities [26], [27]. Calderon preconditioners exploit the selfregularizing properties of the EFIE by invoking the square of the EFIE operator. An extension of the Calderon preconditioner including uniform surface impedance is provided in [28]. However, to the best of the authors’ knowledge, this method has not been extended to MTS antennas and RIS modeled as tensor sheet impedances. A paper proposing a different approach has been recently published in [29]. In the latter, the authors tackle the ill-conditioning for the specific case of scalar and spatially constant sheet impedance by using a regularization scheme based on entire-domain basis functions. The present paper proposes a preconditioner adapted to the analysis of a wide range of modulated sheet impedance lying on a grounded slab, as required for RIS and MTS antennas. Slow convergence has been observed in the past when the impedance range approaches the inductive domain. This includes large perfect electrical conductors (PEC) printed on a grounded substrate, which corresponds to the case where the sheet impedance is zero. As explained in [30], the convergence of an iterative solver such as the Generalized Minimal RESiduad (GMRES) solver can be dramatically impacted by the clustering of the eigenvalues. For example, faster convergence is usually observed when the eigenvalues are located on the positive real axis. Accordingly, the proposed preconditioner significantly improves the convergence of the iterative solver by first redistributing the eigenvalues on the positive real axis after left-multiplying with the conjugate of the original matrix using Fast Fourier Transforms (FFTs), hence with a O(Nlog N)complexity. In a second step, a block diagonal preconditioner, which keeps the eigenvalues on the positive real axis, is applied. Besides, FFTs are also used to precalculate the block diagonal preconditioner. Considering a squareshape domain, this matrix has memory and matrix-vector multiplication complexities of O(N3/2). The fast convergence along with the fast matrix-vector multiplications enables the analysis of finite MTSs. Consequently, a MTS of size larger than ten wavelengths, described with hundreds of thousands of basis functions, can be analyzed in less than one hour. Besides, accurate results are obtained with respect to the ones obtained by direct inversion of the initial system of equations. The presented method can also be easily extended to MTSs with arbitrary shapes as shown in [31]. The paper is structured as follows. Section II briefly recalls a fast MoM formulation adapted to high reactance capacitive IBCs. Section III presents the proposed preconditioner, and Section V provides numerical results. II. METHOD OF MOMENTS FORMULATION Let us consider a RIS or a MTS modeled as a sheet impedance boundary condition (IBC) lying on a grounded substrate, as depicted in Fig. 1. A Cartesian system of coordinates Fig. 1. Schematic illustration of a MTS modeled as a sheet IBC laying on a grounded slab. is adopted with unit vectors ˆxand ˆyon the surface. Using the surface equivalence theorem, the sheet IBC can be replaced by an equivalent surface current distribution expanded into a set of basis functions as: J(r′) = N X n=1 injn(r′),(1) where jn(r′)corresponds to the nth basis function, r′= (x′, y′)is a vector providing the local coordinates on the surface, and inis the unknown weight associated with the nth basis function. The sheet IBC Zsimposes a direct relation between the average electric field Eav on the surface and the equivalent surface current. ˆn ×Eav =ˆn ×Zs·J,(2) where ˆn is the normal to the surface. The electric field to be averaged corresponds to the sum of the incident field Einc, which depends on the excitation, and the fields scattered by the MTS, Escat. The latter can be obtained from the convolution between the equivalent surface current and the medium Green’s function GEJ (r,r′): ˆn ×Eav =ˆn ×(Einc +Escat)(3) =ˆn ×Einc +ˆn ×ZZS′ GEJ (r,r′)J(r′) dS′, (4) where it has been assumed that the MTS sheet is thin enough for the electric field to be considered continuous across the sheet [12]. The electric field is tested with a set of testing functions corresponding to the chosen set of basis functions (Galerkin testing). From there, the following equations are obtained: ZZS jm(r)·Einc dS + N X n=1 inZZS jm(r)ZZS′ GEJ (r,r′)jn(r′) dS′dS = N X n=1 inZZS jm(r)Zs(r)jn(r′) dS;m= 1,2, ..N (5) This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 3 which can be written in the following matrix form: (ZG−ZIBC)I=V(6) where ZGis the substrate matrix whose entries are given by ZG(n, m) = ZZ jn(r′)ZZ GEJ (|r−r′|)jm(r)dS′dS, (7) and where the entries of the IBC matrix ZIBC are given by: ZIBC(n, m) = ZZ jn(r)ZS(r)jm(r)dS. (8) Finally the entries of the excitation vector Vare obtained after testing the excitation field V(m) = −ZZ Einc ·jm(r)dS. (9) In this paper, x-oriented and y-oriented rooftop basis functions [32] are used: jn(r′) =    (r′−rp)·ˆnj hp ˆnj·,if x∈Sp 0,otherwise (10) where the subscript pdenotes the rising (R) or the decaying (D) part of the basis function as shown in Fig. 2; hpis the height of the half basis function. Bold fonts are used to denote vectors. The outer normal to the surface ˆnjcorresponds to ˆx for x-oriented basis functions and to ˆyfor y-oriented basis functions. Sjis the surface of the half basis function. The Fig. 2. x-oriented and y-oriented rooftops basis functions. length of the basis function is set to λ0/10 where λ0is the free-space wavelength. A regular grid of x-oriented and y-oriented basis functions is used as in [16]. A schematic representation of the mesh is given in Fig. 3 where the relative size of the rooftop basis functions is exaggerated for clarity. The regularity of the mesh allows one to exploit the structure of the substrate matrix ZGfor faster solution. Let us consider the mesh as Nxshifted rows of nxx-oriented basis functions and Nyshifted rows of nyy-oriented basis functions; hence, a total of Nx×nxx-oriented basis functions and Ny×nyyoriented basis functions. In this case, the substrate matrix can be divided into four Toeplitz-Block Block-Toeplitz (TBBT) sub-matrices [16]: ZG="Txx Tyx  Txy Tyy #(11) Fig. 3. Illustration of a sheet IBC meshed with rooftop basis functions. where Txx and Tyy are the matrices of self interaction of xoriented and y-oriented basis functions, respectively. The offdiagonal matrices account for the coupling between x-oriented and y-oriented basis functions. Each of these submatrices has a TBBT structure: Tij =          Tij 11 Tij 12 ··· ··· Tij 1Nj Tij 21 Tij 11 Tij 12 ··· Tij 2Nj . . .Tij 21 Tij 11 .... . . . . ..........Tij 12 Tij Ni1Tij Ni2··· Tij 21 Tij 11          (12) where i,jare either xor y. All the submatrices Tij are Toeplitz. Thanks to the structural properties of matrix ZG, the filling of this matrix is fast and grows with linear complexity with respect to the number of unknowns [16]. Indeed, knowing the first column (or row) of each of the four submatrices of ZGis sufficient to know the entire matrix. Regarding the IBC matrix, non-zero entries only concern overlapping basis functions. Hence, the matrix ZIBC is highly sparse. As explained in [16], the structure of the matrices can be exploited to perform fast matrix-vector multiplications in an iterative search of the solution. Multiplying by a TBBT matrix can be seen as a two dimensional convolution [33] and hence can be performed with 2-D FFTs. Besides, multiplying by the IBC matrix is fast given its sparse structure and can be achieved with linear complexity. Overall, a matrix-vector multiplication with the system matrix given in (6) can be achieved with N logN complexity. Provided that the system of equations is well conditioned, this feature can be exploited in iterative solvers such as GMRES [34], as done in [16]. III. PRECONDITIONING AND SOLUTION The linear system of equations (6) is well conditioned when considering MTSs with highly capacitive sheet IBCs [12]. However, the system becomes ill-posed when the sheet impedance gets closer or enters the inductive range [12], [16], [35]. In this paper, we propose the usage of a preconditioner Pdefined as: P=iD (ZG−ZIBC)H,(13) This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 4 where Hdenotes the conjugate transpose operator, and iD is a block diagonal preconditioner defined as: iD ="iDxx 0 0 iDyy #.(14) The matrix iDxx is composed of Nxsubmatrices of size nx×nxand the matrix iDyy is composed of Nysubmatrices of size ny×ny: iDii =         iDii 10··· ··· 0 0 iDii 20··· 0 . . .0....... . . . . .. . ........ . . 0 0 ··· 0 iDii Ni         (15) where iis either xor y. Considering a square domain, this matrix has memory complexity O(N3/2),Nbeing the total number of basis functions. Each block diagonal submatrix in iDii is obtained after inverting the corresponding block diagonal of the matrix Z= (ZG−ZIBC )H(ZG−ZIBC ). In other words, if we define the operator D{·} which selects the first Nxblock diagonals of size nx×nxand the following Nyblock diagonals of size ny×nyof a given matrix, iD is given by: iD =D{Z}−1(16) where the superscript −1denotes the matrix inverse. That means, rather than solving (6), we are actually solving: iD(ZG−ZIBC)H(ZG−ZIBC)I=iD(ZG−ZIBC)HV (17) In an iterative solving of (17), one needs to compute at each iteration a product between the matrix iD ×(ZG− ZIBC )H(ZG−ZIBC ) = iD ×Zand a given vector U. This product can be computed extremely fast without explicitly computing Z. Indeed, first, the product ZU is computed after expanding Zas Z=ZH GZG−ZH GZIBC −ZH IBC ZG+ZH IBC ZIBC (18) while noting that multiplying by ZGHcan also be carried out using 2D FFTs, and multiplication with ZH IBC is extremely fast given its sparsity. The block diagonal preconditioner iD can be obtained by selecting, summing and then inverting the block diagonals of the four terms in (18): iD =D{ZH GZG}−D{ZH GZIBC }(19) −D{ZH IBC ZG}+D{ZH IBC ZIBC }−1 (20) =DGG −DGI −DGI,H +DII −1 (21) =DGG +DIBC −1 (22) where all the IBC dependent matrices have been grouped in DIBC . The first block diagonal of size nx×nxis obtained by inverting the sum of four block diagonals each of size nx×nx corresponding to the first block diagonal of the four matrices in (18): iDxx 1=DGG,xx 1−DGI,xx 1−DGI,xx,H 1+DII,xx 1−1 (23) where the third term is the conjugate transpose of the second one. Let us focus on the first term, which does not depend on the sheet impedance. This term can be obtained by performing a block by block multiplication of matrices ZGand ZGH. Thus, the first block diagonal of DGG can be obtained as follows: DGG,xx 1= Nx X i=1 Txx,H 1iTxx i1+ Ny X i=1 Txy,H 1iTyx i1(24) where each matrix-matrix multiplication can be accelerated with the FFT. It is important to note that, thanks to the block Toeplitz matrix properties, the jth block DGG,xx jcan be obtained from minor additional operations on the previous block DGG,xx j−1, for example: DGG,xx 2= Nx X i=1 Txx,H 2iTxx i2+ Ny X i=1 Txy,H 2iTyx i2(25) =DGG,xx 1−Txx,H 1NxTxx Nx1+Txx,H 12 Txx 21 (26) −Txy,H 1NxTyx Nx1+Txy,H 12 Tyx 21 (27) The terms related to the IBC matrix can also be obtained using FFTs and exploiting the sparse structure of ZIBC. Let us consider an anisotropic sheet impedance, in this case DGI,xx 11 can be obtained as DGI,xx 11 =Txx,H 11 Zxx IBC,11 +Txy,H 11 Zxy IBC,11 (28) where we used the same block designation as in Eq. (11) for the ZIBC matrix and where FFTs can be used to compute the matrix-matrix multiplication. Finally, the block diagonal DII,xx 11 can be obtained by multiplying sparse blocks of the IBC matrix and of its conjugate: DII,xx 11 =Zxx,H IBC,11Zxx IBC,11 +Zxy,H IBC,11Zxy IBC,11 (29) The same strategy can be applied to rapidly obtain all the remaining diagonal blocks of the preconditioner. Once the preconditioner has been calculated, the following new system of equations can be solved using GMRES: iD Z I =iD (ZG−ZIBC)HV,(30) where Zand ZGare never explicitly calculated. IV. EIGENVALUES,PRECONDITIONER BLOCK SIZE AND DISCRETIZATION The proposed preconditioner significantly accelerates the convergence of the iterative solver. One of the reasons is that the eigenvalues of iD Z are clustered on the positive real axis. In this section, we consider three isotropic homogeneous sheets of impedance: (i) highly capacitive (−500j Ω)), (ii) zero impedance (perfect electric conductor (PEC)) and (iii) highly inductive (500j Ω)). We show the distribution of the eigenvalues of both the initial and the proposed system of equations in Fig.4. The three sheets of impedance lie on a This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 5 grounded substrate of relative permittivity ϵr= 3.66 and thickness d= 1.524 mm. The sheet impedance is imposed on a3.2λ0×3.2λ0domain and a frequency of 24 GHz is considered. Considering first the capacitive case, the eigenvalues cluster away from zero and the initial system of equations is usually well conditioned. In the PEC and the inductive cases, the eigenvalues of the MoM system span both the positive and negative parts of the imaginary plane and approach zero, thus leading to ill-conditioning. However, the eigenvalues always cluster on the positive real axis when using the proposed preconditioner. The impact of finer discretization and of the block diagonal preconditioner size on the convergence are shown in Table I. The number of iterations required for GMRES to reach a residual norm of 1e-6 for the three aforementioned sheet impedances. Two discretizations densities are considered: λ0/10 and λ0/5. The required number of iterations is also given when block size of the block diagonal matrix iD is doubled in (15). In this case, the preconditioner is noted P2. In the capacitive case, finer discretization does not seem to impact the convergence while it clearly does for the PEC and inductive cases. Finally, it can be observed that using a preconditioner with larger block size than the one given in (15) does not necessarily reduce the required number of iterations (while it does increase the required memory storage). For the rest of the paper, the discretization is set to λ0/10 and the block size of (15) indicated in Section III is used. TABLE I REQUIRED NUMBER OF ITERATIONS TO REACH A RESIDUAL NORM OF 1E-6. Discretization λ0/10 λ0/5 Preconditioner None P P2None P P2 Capacitive 31 29 29 27 28 27 Metal 526 219 197 404 179 194 Inductive 1001 457 333 73 65 72 V. NUMERICAL EXAMPLES This section provides analysis results of RIS and MTS antennas with the proposed preconditioner. A. RIS analysis Let us consider a RIS operating at 24 GHz. The RIS lies on a grounded slab with relative permittivity ϵr= 3.5and thickness d= 1.524 mm. The dimensions of the MTS is fixed to 7λ0 in the x direction and 4λ0in the y-direction. The isotropic sheet impedance distribution derived from the desired local reflection coefficient, is given by: Zs(x) = −jZ0 tan kx sin θr 2+√ϵrcot (k√ϵrd),(31) with Z0= 376.73 Ω being the free-space impedance and kis the free-space wavenumber at the operating frequency. This sheet impedance distribution provides the linear phase evolution (exp (−jkx sin θr)[36]) required to transform an incident TM polarized beam from broadside to a reflected TM polarized beam in a direction θr, in the plane y= 0. -150 -100 -50 0 Real part 0 100 200 300 Imaginary part Initial system - capacitive IBC New system - capacitive IBC (a) -150 -100 -50 0 Real part -100 -50 0 50 100 Imaginary part Initial system - PEC New system - PEC (b) -150 -100 -50 0 Real part -300 -200 -100 0 Imaginary part Initial system - inductive IBC New system - inductive IBC (c) Fig. 4. Distribution of the eigenvalues for the initial EFIE-IBC system of equations (blue) and for the preconditioned system (red) considering (a) a capacitive IBC, (b) a zero IBC (PEC) and (c) an inductive IBC. A Gaussian incident field Eifrom broadside is considered, with expression at the MTS level given by: Ei(x, y)=e−(x2/(2s2 x)+y2/(2s2 y))ˆx,(32) with sx= 1.4λ0and sy= 0.8λ0. A reflection angle θr= 20◦ is targeted. The convergence of the initial system of equations (6) with and without diagonal preconditioner is compared to the convergence of the new system of equations (30) in Fig 5. As observed, the proposed preconditioner significantly improves the convergence while a conventional diagonal preconditioner, which can be expressed as D{ZG−ZIBC}−1, applied to the initial system of equations deteriorates the convergence. This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 6 0 100 200 300 400 500 Number of iterations 10-6 10-4 10-2 100 Norm of residue Initial system - no preconditioner Initial system - diagonal preconditioner New system - diagonal preconditioner Fig. 5. GMRES convergence: the proposed preconditioned system of equations (29) is compared with the initial system (6) with and without diagonal preconditioner. The radiation pattern and the equivalent currents computed with the proposed algorithm are now compared with the ones obtained using an in house MoM code (based on a direct solution of the system of equations) in Figs. 6 and 7, respectively. Besides, a good agreement is obtained between currents. It can be observed that the broadside specular reflection has been canceled and a pencil beam has been formed at an elevation angle θr= 20◦from broadside. The in-house MoM code has already been validated in [37], [38] where good comparisons were obtained with commercial software. However, the proposed preconditioned iterative algorithm is Fig. 6. Radiation pattern of the 7λ0×4λ0RIS generating anomalous reflection at 20◦. The result of the presented iterative method is compared to direct solving. The error in dB between both results is provided. The inset provides the imaginary part of the sheet impedance. significantly faster and requires much less memory resources. Indeed, on our single conventional computer, the proposed iterative method can deal with MTSs of size 18λ0, meshed with more than 250000 basis functions, in less than 45 minutes. The brute force code is limited to about 30000 basis functions for memory reasons and even then already requires multiple hours of calculation. Next, a broadside impinging plane wave excitation is considered and the system of equations is again solved. In this case, the specular reflection has been obtained Fig. 7. Equivalent currents obtained on the RIS surface using (a) the iterative method and (b) direct solving. The relative error between both currents in dB is shown in (c). -50 -20 0 20 50 -20 -10 0 10 20 30 Directivity [dBi] Directivity iterative Directivity direct error 200 600 1000 10-5 10-2 Fig. 8. Radiation pattern of the 7λ0×4λ0RIS generating anomalous reflection at 20◦when excited by a broadside impinging plane wave. The inset shows the convergence of GMRES with the proposed system of equations. with the PO approximation as explained in [39]. The resulting radiation pattern is shown in Fig. 8 along with an inset showing the convergence of the proposed system of equations. The required number of iterations is three times larger than the one associated with the Gaussian beam excitation showing that a different excitation vector can impact the number of iterations. However, the system of equations still converges very well after preconditioning. RISs implemented with the sheet impedance given by (31) usually produce undesired sidelobes which can be related to additional Floquet modes [40]. The RIS shown in the inset of Fig. 6 supports approximately two periods of the sheet impedance. Consequently, the sidelobes corresponding to the Floquet modes do not clearly appear due to the small size of the RIS. However, the modes start to appear when considering larger RISs. For example, the radiation pattern of Fig. 9 is obtained with a RIS supporting 10 periods. In this case the sidelobes corresponding to Floquet modes clearly appear. This RIS is analyzed with 63920 basis functions. Given the high number of unknowns, direct solution of the system of equations is impractical on a conventional computer. However, using the presented method, the iterative solver reaches a norm of residue of 10−6in about 5 minutes. The performance of the algorithm is now analyzed with This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 7 TABLE II EVALUATION TIME OF THE PROPOSED METHOD ON A SINGLE STANDARD COMPUTER Size of the RIS Number of unknowns Green’s function tabulation [min] ZGfilling time [min] DGG filling time [min] ZIBC filling time [min] DIBC filling time [min] Perform block inverse in (19) [min] GMRES solving time [min] 10 λ0×10 λ079600 0.938 0.33 0.044 0.015 0.03 0.022 3.58 12 λ0×12 λ0114720 1.12 0.48 0.094 0.021 0.06 0.04 5.58 14 λ0×14 λ0156240 1.4 0.7 0.16 0.033 0.1 0.053 12.8 16 λ0×16 λ0204160 1.45 0.9 0.25 0.038 0.15 0.078 24.9 18 λ0×18 λ0258480 1.76 1.075 0.284 0.045 0.2 0.1 42.6 Fig. 9. Radiation pattern of a 10 periods RIS generating anomalous reflection at 30◦, considering a Gaussian field incident from broadside. The inset provides the imaginary part of the sheet impedance. respect to the size of the RIS. To do so, the sheet impedance of the previous example is selected and the domain of the RIS is gradually increased while observing the simulation time and the convergence of the method. We consider a square-shape RIS and we solve the system of equations (30). It is shown in Fig. 10 that the method still converges well when the size of the structure and hence the number of unknowns increases. The main steps and their respective computational times using a single standard computer (Intel Core i5-7500 3.4 GHz and 24 GB of RAM) are given in Table II. B. Leaky-wave MTS antenna analysis This section is devoted to the analysis of leaky-wave MTS antennas [8]. This class of antenna usually implements a modulated capacitive sheet impedance. As shown in [16], a low reactance capacitive MTS exhibits a slow convergence due to the poor conditioning of the system of equations matrix. A low reactance capacitive sheet impedance may be required to increase the feed efficiency [19] or to facilitate the surface impedance implementation with an adopted set of patches. In this section, the convergence of the proposed method is analyzed considering squared-shape MTSs operating at 19 GHz with a low reactance capacitive sheet impedance laying on a substrate with relative permittivity ϵr= 6 and thickness d= 1.5mm. The antenna is excited with a vertical elementary dipole placed in the middle of the substrate and in the center 0 200 400 600 800 1000 Number of iterations 10-6 10-4 10-2 100 Norm of residue 8 0 8 0 10 0 10 0 12 0 12 0 16 0 16 0 18 0 18 0 Fig. 10. Convergence of the method with respect to the structure size. of the aperture. An anisotropic sheet impedance profile is considered [41]: Zρρ =jX0[1 + Mcos(kswρ−ϕ)] Zρϕ =Zϕρ =jX0Msin(kswρ−ϕ)(33) Zϕϕ =jX0[1 −Mcos(kswρ−ϕ)] , with the following parameters: X0=−222.42 Ω,M= 0.5, ksw = 1.8k0, where k0is the free-space wavenumber. (ρ, ϕ) corresponds to the polar system of coordinates centered on the MTS aperture. The Zρρ component of this impedance profile is depicted in Fig.13. The GMRES solver is applied with two different sizes for the square-shaped MTS antenna: 5λ0and 7λ0, meshed with 19800 and 38920 rooftop basis functions, respectively. The convergence of the preconditioned system of equations (30) is compared with the one of the initial system of equations (6) in Fig. 11, the latter with and without diagonal preconditioner. As observed, a faster convergence is obtained using the proposed preconditioner. Again, the application of a simple diagonal preconditioner to the initial system of equations (6) significantly deteriorates the convergence. The radiation pattern obtained after analyzing the MTS of size 5λ0is now compared with that resulting from a direct inversion of the system of equations matrix (6) in Fig. 12. Indeed, this MTS is described with 19800 basis functions; a system of equations of that size can still be directly solved on a standard computer. An excellent agreement between the This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 8 0 200 400 600 800 1000 Number of iterations 10-8 10-6 10-4 10-2 100 Norm of residue New system, 5 5 Initial system, 5 5 Initial system, 7 7 New system, 7 7 5 5 - diagonal prec 7 7 - diagonal prec Fig. 11. Convergence of the proposed method compared to the convergence of the initial system of equations with and without diagonal preconditioner for the low reactance capacitive square-shaped MTS antenna. -50 0 50 -60 -40 -20 0 20 Directivity [dBi] co-pol iterative x-pol iterative co-pol direct x-pol direct error co-pol error x-pol Fig. 12. Radiation pattern of the low reactance capacitive square-shaped MTS antenna of size 5λ0×5λ0, obtained with the iterative solver and with the direct method. The absolute error between patterns (in dB) is also provided. Fig. 13. Imaginary part of the Zρρ component of the impedance profile given in (33). patterns is observed, thus validating the proposed iterative solution. C. Sparse Array Metasurface In this section, we use the proposed preconditioned to analyse a 5λ×14λMTS for which the initial MoM system of equations is ill-conditioned and we validate the result at the printed patch level. Therefore, we consider the design proposed in [42] with the scalar surface impedance given by [42, Equ. (8)]. This impedance profile spans both capacitive and inductive ranges and leads to an ill-conditioned IBC-EFIE system of equations. The proposed preconditioner is used to solve the latter. The impedance profile and the corresponding surface currents are shown in Figs. 14 (a) and (b), respectively. The resulting radiation pattern is shown in Fig. 15 where it is (a) (b) Fig. 14. Sparse array MTS proposed in [42]: (a) imaginary part of the impedance profile, (b) normalized equivalent currents. -50 -20 0 20 50 -20 -10 0 10 20 Directivity [dBi] IBC model Printed patches 200 600 900 10-5 10-2 Fig. 15. Radiation pattern of the sparse array MTS proposed in [42] obtained with the proposed preconditioner at the IBC level and compared with the one obtained at the printed patch level in [42]. The inset shows the convergence (norm of residue vs number of iterations) of the present method (blue line) compared to the convergence of the initial system of equations (red line). compared to the one obtained considering a full-wave analysis of the metasurface described at the patch level, which is shown in Fig. 16. More details regarding the patch implementation of the sparse array MTS can be found in [42]. In Fig. 15, an inset is provided with the convergence of the proposed method. This example shows that the method can be used to solve an ill-conditioned IBC-EFIE system of equations for This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ IEEE TRANSACTIONS ON ANTENNAS AND PROPAGATION, VOL. X, NO. XX, 20XX 9 Fig. 16. Implementation of the sparse array MTS proposed in [42] at the patch level. which the solution agrees well with the one of the EFIE when considering the printed patch implementation. VI. CONCLUSION Analyzing structures at the sheet impedance level has often led to ill-conditioning when the range of impedance spans both the capacitive and inductive domains. Therefore, a well conditioned Method of Moments (MoM) technique is proposed for the fast analysis of metasurfaces modeled as a modulated sheet impedance boundary condition. The solution of the system of equations matrix is obtained iteratively with the Generalized Minimal RESidual (GMRES) method. At each iteration, FFTs are used to accelerate the matrix-vector multiplication with the EFIE-IBC matrix and with its conjugate (O(Nlog N complexity). An additional matrix-vector multiplication with a block diagonal matrix is performed (O(N3/2)complexity considering a square-shaped domain). This matrix is precalculated using FFTs. The proposed preconditioner ensures a rapid convergence of the algorithm whatever the range of the sheet impedance profile.As shown in the paper, this method can be used to analyze large MTSs or RISs described with more than 250000 basis functions with a conventional computer. REFERENCES [1] X. Cao, Q. Chen, T. Tanaka, M. Kozai, and H. Minami, “A 1-bit timemodulated reflectarray for reconfigurable-intelligent-surface applications,” IEEE Trans. Antennas Propag., vol. 71, pp. 2396-2408, March 2023. [2] M. Faenzi, G. Minatti, D. Gonz´ alez-Ovejero, G. Caminita, E. Martini, C. Della Giovampaola, and S. Maci “Metasurface antennas: new Models, applications and realizations,” Sci. Rep., vol. 9, July 2019. [3] Q. Bu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: a tutorial,” IEEE Com. Surveys and Tutorials, vol. 69, pp. 3313-3351, May 2021. [4] M. Bodehou, K. Alkhalifeh, S. N. Jha, and C. Craeye, “Direct numerical inversion methods for the design of surface-wave based metasurface antennas: fundamentals, realization, and perspectives,” IEEE Antennas Propag. Mag., vol. 64, pp. 24-36, Aug. 2022. [5] Y. Xie, C. Yang, Y. Wang, Y. Shen, X. Deng, B. Zhou, and J. Cao, “Anomalous refraction and reflection characteristics of bend V-shaped antenna metasurfaces,” Scientific Reports, vol. 9, April 2019. [6] H. F. Ma, Y. Q. Liu, K. Luan, and T. J. Cui,“Multi-beam reflections with flexible control of polarizations by using anisotropic metasurfaces,” Sci. Rep., vol. 6, Nov. 2016. [7] J. L. Keyrouz, V. G. Ataloglou, and G. V. Eleftheriades “Experimental demonstration of a dual-input/dual-output reflective impedance metasurface,” Appl. Phys. Lett., vol. 123, Oct. 2023. [8] G. Minatti, M. Faenzi, E. Martini, F. Caminita, P. De Vita, D. Gonz´ alezOvejero, M. Sabbadini and S. Maci, “Modulated metasurface antennas for space: synthesis, analysis and realizations,” IEEE Trans. Antennas Propag., vol. 63, pp. 1288-1300, Apr. 2015. [9] D. Gonz´ alez-Ovejero, G. Minatti, G. Chattopadhyay, and S. Maci, “Multibeam by metasurface antennas,” IEEE Trans. Antennas Propag., vol. 65, pp. 2923-2930, June 2017. [10] M. Bodehou, E. Martini, S. Maci, I. Huynen, and C. Craeye, “Multibeam and beam scanning with modulated metasurfaces,” IEEE Trans. Antennas Propag., vol. 68, pp. 1273-1281, 2020. [11] C. Yepes, S. Maci, S. A. Tretyako, and, E. Martini, “On the role of spatial dispersion in boundary conditions for perfect non-specular reflection,” EPJ Applied Metamaterials, vol. 9, p. 17, Jul. 2022. [12] M. A. Francavilla, E. Martini, S. Maci, and G. Vecchi, “On the numerical simulation of metasurfaces with impedance boundary condition integral equations,” IEEE Trans. Antennas Propag., vol. 63, pp. 2153-2161, May. 2015. [13] F. Verni, M. Righero, and G. Vecchi, “On the use of entire-domain basis functions and fast factorizations for the design of modulated metasurface,” IEEE Trans. Antennas Propag., vol. 68, pp. 3824-3833, Jan. 2020. [14] D. Gonz´ alez-Ovejero and S. Maci, “Gaussian ring basis functions for the analysis of modulated metasurface antennas,” IEEE Trans. Antennas Propag., vol. 63, pp. 3982-3993, Sep. 2015. [15] M. Bodehou, D. Gonz´ alez-Ovejero, C. Craeye and I. Huynen, “Method of Moments simulation of modulated metasurface antennas with a set of orthogonal entire-domain basis functions,” IEEE Trans. Antennas Propag., vol. 67, pp. 1119-1130, 2019. [16] J. Cavillot, M. Bodehou, S. Hubert, and C. Craeye, “Efficient analysis of planar, arbitrarily shaped and (bi)-anisotropic metasurface antennas,” IEEE Trans. Antennas and Propag., vol. 70, no. 1, pp. 536-546, Jan. 2022. [17] R. F. Harrington, Field Computation by Moment Methods. Piscataway, NJ, USA: IEEE Press, 1993. [18] C. Yepes, M. Faenzi, S. Maci, and E. Martini, “Perfect non-specular reflection with polarization control by using a locally passive metasurface sheet on a grounded dielectric slab,” Appl. Phys. Lett., vol. 118, no. 23, Jun. 2021. [19] G. Minatti, E. Martini, and S. Maci, “Efficiency of metasurface antennas,” IEEE Trans. Antennas Propag., vol. 65, pp. 1532-1541, Apr. 2017. [20] S. B. Adrian, A. D´ ely, D. Consoli, A. Merlini and F. P. Andriulli, “Electromagnetic integral equations: insights in conditioning and preconditioning,” IEEE Open J. Antennas Propag., vol. 2, pp. 1143-1174, Dec. 2021. [21] R. W. Kindt, K. Sertel, E. Topsakal and J. L. Volakis, “Array decomposition method for the accurate analysis of finite arrays,” IEEE Trans. Antennas and Propag., vol. 51, no. 6, pp. 1364-1372, Jun. 2003. [22] P. Janpugdee, P. H. Pathak, P. Mahachoklertwattana, and R. J. Burkholder, “An accelerated DFT-MoM for the analysis of large finite periodic antenna arrays,” IEEE Trans. Antennas Propag., vol. 54, no. 1, pp. 279-283, Jan. 2006. [23] D. Pissoort, E. Michielssen, D.V. Ginste, and F. Olyslager, “A rankrevealing preconditioner for the fast integral-equation based characterization of electromagnetic crystal devices,” Microw. Opt. Technol. Lett., vol.48, no.4, pp. 783-789, Apr. 2006. [24] D. Tihon, S. Hubert, C. Craeye and N. A. Ozdemir, “SVD postcompression combined with shielded-block preconditioner,” in Proc. Int. Conf. Electromag. in Ad. App. (ICEAA), Palm Beach, Aruba, 2014, pp. 503-506. [25] S. Hubert, S. N. Jha, and C. Craeye, “Analysis of large nonregular printed scatterers using the Contour-FFT,” IEEE Trans. Antennas Propag., vol. 66, no. 11, pp. 6115-6127, Nov. 2018. [26] G. C. Hsiao and R. E. Kleinman, “Mathematical foundations for error estimation in numerical solutions of integral equations in electromagnetics,” IEEE Trans. Antennas Propag., vol. 45, no. 3, pp. 316-328, Mar. 1997. [27] F. P. Andriulli, K. Cools, H. Ba˘ gci, F. Olyslager, A. Buffa, S. Christiansen, and E. Michielssen, “A multiplicative Calderon preconditioner for the electric field integral equation,” IEEE Trans. Antennas Propag., vol. 56, no. 8, pp. 2398-2412, Aug. 2008. This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2025.3617109 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/