Field Scattering Analysis of Cylindrical Spatially Dispersive Metasurfaces

In this letter, a semianalytical technique to simulate spatially dispersive (SD) cylindrical metasurfaces is proposed and demonstrated. This method is based on the well-known eigenfunction expansion (EFE) method and extends the works of Smy et al. (2023) by integrating the extended generalized sheet transition conditions into this preexisting method specific to cylindrical surfaces. This method is then applied to several cylindrical metasurfaces composed of SD wire dipole unit cells, and successfully verified against commercial full-wave simulations. Finally, using a case of a closed metasurface shell, the usefulness of the proposed EFE-SD method is shown to illustrate the effect of metasurface thickness on the scattered fields inside and outside the shell.


Field Scattering Analysis of Cylindrical Spatially Dispersive Metasurfaces
Jordan Dugan , João G. Nizer Rahmeier , Student Member, IEEE, Tom J. Smy , and Shulabh Gupta , Senior Member, IEEE Abstract-In this letter, a semianalytical technique to simulate spatially dispersive (SD) cylindrical metasurfaces is proposed and demonstrated.This method is based on the well-known eigenfunction expansion (EFE) method and extends the works of Smy et al. (2023) by integrating the extended generalized sheet transition conditions into this preexisting method specific to cylindrical surfaces.This method is then applied to several cylindrical metasurfaces composed of SD wire dipole unit cells, and successfully verified against commercial full-wave simulations.Finally, using a case of a closed metasurface shell, the usefulness of the proposed EFE-SD method is shown to illustrate the effect of metasurface thickness on the scattered fields inside and outside the shell.Index Terms-Cylindrical metasurfaces (MS), eigenfunction analysis, electromagnetic metasurfaces, electromagnetic propagation, generalized sheet transition conditions (GSTCs), surface susceptibility tensors.

I. INTRODUCTION
E LECTROMAGNETIC metasurfaces (MS) are 2-D arrays of subwavelength resonators that can be engineered to produce a desired macroscopic field scattering response [3] and have found many applications [4], [5], [6], [7], [8], [9].A full-wave simulation of these structures is very computationally demanding due to the fine meshing needed, and this has motivated the use of zero-thickness sheet models characterized by susceptibility tensors [10], [11].Fields around such a surface must satisfy the generalized sheet transition conditions (GSTCs) [12], [13].This compact model allows for a much more efficient simulation of the surface making it very useful in MS design and analysis, and has been integrated into several numerical techniques [14], [15], [16], [17], [18], [19].
While all these numerical methods may be generally applied to arbitrarily conformal geometries, differences in the scattered fields obtained from the zero-thickness sheet models or the physical MS structures can be attributed to two possible sources: 1) the imperfect convergence of numerical simulations; or 2) the validity of a GSTC-based zero-thickness sheet model representation of curved surfaces itself.To systematically address this, a semianalytical technique was recently proposed and demonstrated in [1] and [2] to simulate cylindrical MS by solving the GSTCs using the well-known eigenfunction expansion (EFE) technique [20], [21], [22], [23], [24].The benefit of the method developed in [1] and [2] over the other methods is the fact that it has been systematically developed to simulate MS with full tensorial susceptibilities, including both normal and tangential components.Scattered fields computed using the EFE method can be considered as the correct answer, due to its inherent analytical nature.This method, thus, provided an important benchmarking tool, in addition to being an efficient field solver for cylindrical geometries with deep subwavelength unit cell geometries, to establish limits and validities of GSTCs-based MS modeling.
While deep subwavelength unit cell periods feature modeling conveniences, practical MS structures are quite often large, reaching up to the λ/2 limit and yet are homogenizable.Induced currents on deeply subwavelength structures are locally related to fields around the surface, and thus can be modeled with angle-independent dipolar susceptibilities (referred to as spatially nondispersive or local).This is, however, not the case for larger unit cells, which exhibit significant spatial dispersion (or nonlocality) whose field scattering characteristics are captured using angle-dependent surface susceptibilities.
Recently, the extended GSTC method of modeling nonlocal MS has been proposed whereby the susceptibilities take the form of rational polynomial functions in the spatial frequency domain [25] and implemented in an integral equation (IE) solver in [26].While this IE-SD method performs very well, it is inherently approximate due to the use of discrete surface meshing, making it difficult to discern between numerical artifacts and physical features of the fields.In this work, we, therefore, develop an EFE-based semianalytical method to solve the extended GSTCs for computing the scattered fields from SD cylindrical MS, whose unit cells are not necessarily deeply subwavelength.This will be referred to as the EFE-SD method.

II. PROPOSED EFE-SD METHOD
The use of the GSTCs to model the cylindrical MS assumes it is a zero-thickness sheet producing a discontinuity in the tangential components of the E and H fields across the surface.Under these conditions, the MS is rigorously described by the following two GSTC equations [11], [12]: 1536-1225 © 2023 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.See https://www.ieee.org/publications/rights/index.html for more information.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
where ΔF = F + − F − , P and M are the electric and magnetic surface polarizations, respectively, n is the direction normal to the surface, and P n and M n are the components of the polarizations normal to the surface.

A. Spatially Dispersive GSTCs
The polarization present on a zero-thickness shell of an SD MS can be determined through a convolution of the fields with the susceptibilities over the surface where E av and H av are the average fields at the MS.Each of the susceptibility terms present in these equations are tensors and creates a contribution to the polarizations in (2); for example, χ zz ee creates a component for an SD unit cell given by P z z = 0 χ zz ee * E z,av , where the superscript on P z z indicates this is the portion of P z generated by E z,av and the convolution (which captures the spatial dispersion) is along the surface.The total P z would, thus, be (assuming no polarization rotation and bianisotropy) For a flat surface expressed in a Cartesian system, following the work in [25] and [26], we express the susceptibilities as a ratio of polynomials in a spatial frequency domain k y (assuming that the surface normal is x).For clarity, if we consider TE incidence, only χ zz ee is present Setting b 0 = 1 and a 0 = χ zz ee,0 , for simplicity and sake of demonstration, we obtain ( This can be brought back in to the spatial domain using an inverse Fourier transform 1Pz + Each portion of the field differences can be associated with a polarization component generated from a susceptibility term, for example, Δ Hy = jω Pz .Substituting this in the abovementioned equation results in the extended form of the GSTCs (second GSTC is not shown here for brevity) which incorporates nonlocal SD effects.Here, we see that the extended GSTCs lead to a higher order boundary condition model of the MS.Such boundary conditions have been used previously in electromagnetics, for example, in [28] and [29].
The difference between our works on spatial dispersion [25], [26], [27] and these works is that these works develop boundary conditions for thin dielectric layers with specific geometry only, whereas our works on the extended GSTCs allow us to derive the boundary condition for an array of scatterers with any arbitrary geometry.

B. Eigenfunction Scattering Analysis
To perform an EFE analysis of a cylindrical structure, we first write (6) in cylindrical coordinates, r = (ρ, φ, z), using dy = r 0 dφ (see Fig. 1) and for simplicity, limit the sum to second order and, assuming a symmetrical unit cell, dropping odd ordered terms where the total fields (H φ , E z ) are the sum of the scattered (H s φ , E s z ) and incident (H i φ , E i z ) fields, respectively.Exploiting rotational periodicity, the scattered and incident fields are taken to have the form of [30], [31] where R(ρ, n) is a radial function that characterizes the interior and exterior fields, and k 0 is the free space wave vector.The total field solution has to satisfy interface conditions for the tangential fields at the surface.The scattered field must satisfy the radiation condition at ρ = ∞.The interior (F − ) fields have to be finite at ρ = 0, which limits the interior radial function to the Bessel function, J n (ρ), only.The exterior (F + ) scattered fields can be represented by Hankel functions of the second kind, H components of E and H, and we have )e jnφ which are characterized by four sets of unknown constants (a n , b n , c n , d n ) to be determined by the nature of the incident field present and the surface conditions.Other field components can be determined from the z-field components via Maxwell's equations (see [

1, Appendix A]).
Each component of the surface susceptibility on the cylindrical surface can also be described by an angular Fourier series (FS) due to implied angular periodicity of 2π, which is given by, for example As (7) has terms that involve products of implicit sums, we note that representing this equation in the angular domain involves an implied discrete convolution along the surface, i.e., in φ.For example, the product of χ zz ee,0 (φ) and the average E z can be represented as follows: Using this relationship and that for a field prescribed by ( 8), we have We can now equate terms of the order n in (7) to obtain a set of equations given by From ( 1), similar equations for the other tangential field components H z , E φ , and E z can be developed.For example, for an incident polarization of E z , H φ , and H ρ , we have a second tangential equation Equations ( 9) and ( 10) must now be solved self-consistently for the unknown coefficients of Bessel and Hankel functions forming the scattered fields, assuming the FS form of the susceptibilities.Limiting the number of terms by setting a range of n from −N to N in the FS expansion, gives us a consistent set of equations we can put into matrix form (see [1] and [2] for details) and solve for the unknown coefficients (a n , b n , c n , and d n )

III. EFE-SD METHOD DEMONSTRATION
To demonstrate the EFE-SD formulation, here we will apply our method to a practical MS structure composed of a lossy dipole in free space, as shown in Fig. 2(a).The angle-dependent susceptibilities of this dipole were extracted at a frequency of 60 GHz in [27] and are provided in Fig. 2(a) (note that the susceptibilities for a planar and cylindrical MS are the same).This unit cell was arranged in a half-cylinder formation (i.e., an open configuration) with radius r = 20.53mm (arbitrarily chosen) and simulated using Ansys finite element method high-frequency structure simulator (FEM-HFSS).The total fields predicted by this simulation are shown in the first plot of Fig. 2(b).Microscopic fields close to the wire surface are clearly visible near the structure, while a wave focussing effect with an intricate interference pattern is formed in the macroscopic field region.
The equivalent zero-thickness sheet was next simulated using the EFE method with and without spatial dispersion.The total fields of each simulation are plotted in Fig. 2(b) for comparison with those obtained from FEM-HFSS.There is an excellent match between brute-force HFSS fields and the proposed EFE-SD method, which deteriorates and becomes unacceptable when the SD is removed.To further demonstrate a significant improvement in accuracy when spatial dispersion is considered, the total fields predicted using each method are plotted along an interior circle with a radius of 0.8r at the bottom of Fig. 2(b).This not only shows that extended GSTCs are necessary to model this structure, but also demonstrates the excellent accuracy of the proposed EFE-SD method.Moreover, due to its semianalytical nature, the computation time is orders of magnitude faster than brute-force HFSS simulations (e.g., few tens of seconds versus hours on a typical desktop computer).
As a second example, we consider a complete cylindrical surface of radius r (i.e., a closed configuration).The closed nature of the surface introduces pronounced resonant patterns inside the MS shell, which are quite sensitive to any imperfections in the numerical modeling.We first consider a small sized cylinder of r = 10.25 mm whose scattered fields obtained from FEM-HFSS and EFE-SD are compared in Fig. 3, showing a remarkable agreement in both inside and outside regions.The fields predicted by FEM-HFSS and EFE-SD and EFE with no SD are plotted along the line y = 0 in Fig. 3 confirming very good accuracy in reproducing rapid changes in the field pattern when SD is included and significant disagreement when it is not.Now that we have demonstrated the accuracy of the EFE-SD method with good success, we can explore its usefulness beyond its rapid field modeling capability.We next consider the closed cylindrical surface of twice the size of that in Fig. 3, with an r = mm.This forms an electrically large resonator, which is very sensitive to any variations of its radius.The EFE-SD method uses an equivalent zero-thickness sheet, which now must be placed with care with respect to the physical thickness d 0 of the MS.Fig. 4 shows the EFE-SD computed scattered fields for three different radii (r 0 ) with a relative change of less than 1%.These three radii values correspond to the zero-thickness sheet placed at the inner edge (r 0 = r − d 0 /2), center (r 0 = r), and the outer edge (r 0 = r + d 0 /2) of the dipole array.We observe that such a slight variation in the sheet placement leads to rapid variations in the scattered fields.The MS thickness of d 0 = λ 0 /25 in this resonator configuration, which, otherwise, can be safely considered deeply subwavelength, can no more be ignored.We further notice that for the case of r 0 = r − d 0 /2, the internal fields are much better reproduced compared with the outside fields (exhibiting phase shift errors), with good agreement with HFSS computed fields.This is expected because the susceptibilities were extracted with the reference plane at the edges of the unit cell in the Floquet unit cell simulations in FEM-HFSS.This example, thus, illustrates how the EFE-SD method can reveal subtle aspects and limitations of zero-thickness sheet modeling of practical SD structures, which, otherwise, remain hidden under numerical uncertainties in standard modeling techniques.

IV. CONCLUSION
In this letter, a rigorous semianalytical method for simulating SD cylindrical MS using the EFE technique has been proposed.This letter integrates the extended GSTCs into the preexisting EFE solver of [1] and [2], giving it the capability to analyze nonlocal cylindrical surfaces.This EFE-SD method using various cylindrical surfaces composed of a wire unit cell was successfully validated against brute-force HFSS simulations showing excellent accuracy, in addition to inherent benefit of computing speed.Using the EFE-SD method, the fields within closed electrically large structures were further explored to reveal sensitivities of zero-thickness sheet placements compared with MS thickness.The proposed method is fast, simple, and does not suffer from the nonphysical numerical artifacts that can arise from the discrete meshing required for most numerical techniques, making it an excellent candidate to test the limits of GSTC-based modeling, as well as a benchmarking solution for alternative numerical modeling techniques.

Fig. 1 .
Fig.1.Basic geometry for a single-shell cylindrical MS configuration.The example consists of a partial cylindrical surface characterized by angle-dependent surface susceptibilities that account for spatial dispersion or nonlocality.

( 2 )
n (ρ)-outgoing propagating waves.The scattered field configuration of the two polarizations can be captured by the z Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Fig. 2 .
Fig. 2. Field scattering from a circular array of electric dipole wires exhibiting spatial dispersion.(a) Unit cell configuration with its dispersion parameters for χ zz ee is available in [27, Table I].(b) Half-cylindrical wire array with a TE plane-wave excitation from the left is obtained using full-wave HFSS simulations, and the corresponding EFE simulations with and without spatial dispersion.Also, a plot of the field magnitude along an interior circle at 80% of the radius is shown for comparison.The wire parameters are d 0 = 0.2 mm, = 2.3 mm, Λ y = 2.15 mm, Λ z = 4.3 mm, conductivity σ = 5.8 × 10 4 S/m, and r = 1.

Fig. 3 .
Fig. 3. Comparison of scattered fields of a complete cylindrical surface (r = 10.25 mm) of electric dipole wires computed from a brute-force FEM-HFSS simulation and the EFE-SD zero-thickness sheet model.

Fig. 4 .
Fig.4.Scattered fields, Re{E z (x, y)}, from a complete cylindrical array of wire dipoles computed using FEM-HFSS and EFE-SD sheet models for varying radius r 0 = r ± d 0 /2, where d 0 is the thickness of the wire and r is the radius of the array measured from the center of the dipole unit cell.