Application of cylindrical IE-GSTC to physical metasurfaces

The cylindrical Integral Equation -- Generalized Sheet Transition Condition (IE-GSTC) is validated by applying it to physical metasurfaces, i.e. metasurfaces defined by material properties and dimensions rather than by susceptibility tensor components. The previously reported IE-GSTC, which was formulated for zero thickness GSTC discontinuities, is extended to handle the finite thickness of physical metasurfaces. A simple analytical approach is used to extract the bianisotropic susceptibility tensor of concentric, multilayered, magneto-dielectric shells. Plane wave scattering by a physical metasurface constructed of four segments of multilayered, magneto-dielectric metasurface scatterers is used as an example problem to validate cylindrical IE-GSTC. A second example considers an opening on the cylindrical metasurface, confirming IE-GSTC can handle metasurfaces with openings. Third example is that of plane wave scattering by an eight segment metasurface. Good agreement is obtained between IE-GSTC results and full wave simulation results for all the cases. The paper concludes with a novel technique to handle PEC segments embedded in a cylindrical metasurface.


I. INTRODUCTION
M ETASURFACES, the 2D version of volumetric metamaterials, have provided unprecedented ability in the control and manipulation of electromagnetic waves. In comparison to metamaterials, metasurfaces are easier to fabricate. Metasurfaces are electrically thin surfaces constructed using subwavelength particles called metaatoms. By using engineered particles, metasurfaces can transform the phase, amplitude, polarization and propagation direction of the incident wave. Several review articles and books exist on the design, modeling and applications of planar metasurfaces [1]- [4]. Due to its deeply subwavelength thickness, physical metasurfaces can be abstracted as an infinitesimally thin layer of electric and magnetic surface polarization densities. Metasurfaces achieve their functionality by creating a spatial electromagnetic discontinuity. Mathematically, the discontinuity can be expressed by Generalized Sheet Transition Conditions (GSTCs) which relate the electric and magnetic field discontinuities to the electric and magnetic surface polarization densities [5], [6].
S. Sandeep is currently with the Department of Electronic Systems, NTNU, Trondheim, Norway. e-mail: sandeepsrikumar2013@gmail.com A. J. Gasiewski is with the Department of Electrical and Computer Engineering, Unversity of Colorado, Boulder, CO, 80309 USA.
S. Y. Huang is with the Singapore University of Technology and Design (SUTD), Singapore.
A. F. Peterson is with the School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA, 30332 USA. e-mail: peter-son@ece.gatech.edu The surface polarization densities are related to the average of the fields on either side of the GSTC surface through the bianisotropic susceptibility tensor model [7]. It should be noted that even though physical metasurfaces have a finite subwavelength thickness, GSTCs which assume zero metasurface thickness are a very good approximation to such physical metasurfaces, i.e. the scattered field due to the physical metasurface approximates the field scattered by GSTC interface conditions. GSTCs in combination with the bianisotropic susceptibility tensor model can be effectively used to simulate thin electromagnetic surfaces. Previous works in this approach include the application of a Finite Difference Frequency Domain (FDFD) method [8] and a Finite Element Method (FEM) [9] to the modeling of planar GSTCs. Even though most of the metasurfaces reported to date are planar, other canonical shapes such as cylindrical metasurfaces [10]- [13] and spherical metasurfaces [14], [15] are now being studied. Applications of conformal metasurfaces include illusion transformation [13], cloaking [12], high-gain antennas [16] and so on. GSTC bianisotropic susceptibility tensor model had been previously used for conformal metasurfaces. An Integral Equation -Method of Moment (IE-MoM) approach was applied to circular cylindrical metasurfaces in [10]. This was extended to cylindrical metasurfaces of arbitrary cross-section in [11]. Conformal transformation was used to synthesize surface electromagnetic tensor quantities for noncircular cylindrical metasurface in [17]. In [18], a digital signal processing based method is outlined for analysis and synthesis of cylindrical, bianisotropic metasurfaces. In [19], an inverse source technique was used to synthesize surface susceptibility tensor components of metasurfaces of arbitrary shape. The work addressed metasurfaces capable of generating radiation patterns with specified half-power beamwidth, nulls, polarization etc. Recently, asymptotic methods (GO, GTD) were applied to electrically large metasurfaces [20]. Even though mathematically rigorous, none of these works [10]- [12], [14], [19]- [23] have used physical metasurfaces.
One of the goals of this paper is to move from the bianisotropic susceptibility tensor based conformal GSTC surface models used in [10]- [12], [14], [19]- [23] to physical metasurfaces formed by curved metasurface scattering particles. In this work, validity of the cylindrical GSTC [10] is confirmed by comparing the scattering due to a cylindrical, physical metasurface to scattering by a corresponding GSTC contour. Scattering by a cylindrical GSTC contour was obtained by a slightly modified version of IE-GSTC formulation in [10].
The modification incorporates the finite thickness of physical metasurfaces. This modification is important to get accurate results for IE-GSTCs [10], [11], [22] when compared to physical metasurface scattering.
The organization of the paper is as follows. Section II extends the IE-MoM in [10] to handle a physical metasurface with non-zero thickness. The contents of this section closely follow [10]. Only the significant differences are explained. This is followed in section III by susceptibility tensor extraction of cylindrical, concentric, magneto-dielectric shell using an analytical approach. Section IV provides two examples. The metasurfaces in the two examples are constructed using metasurface particles which are part of the shell structure in section III. The second example has an opening which can be modeled using the IE-GSTC formulation, but with reduced accuracy. The IE-GSTC, using susceptibility tensors extracted using the method in section III, is used to simulate the scattering problem, with the results compared with a full wave simulation of the physical metasurface. The final section provides a detailed formulation for handling PEC segments in cylindrical metasurfaces with one numerical validation example.

METASURFACE
The IE-MoM described in [10] was limited to a cylindrical GSTC contour of zero thickness whose surface polarization densities are related to electric and magnetic fields on either side of the metasurface by the bianisotropic susceptibility tensor model. That formulation needs to be slightly extended to include the non-zero thickness of a physical metasurface. Instead of using a single contour for calculating MoM coefficients and integration of surface fields, two contours (one external and other internal to the metasurface) are used. This is shown in Fig. 1. Throughout this paper, there is no variation of material properties or excitation along the z axis. The inner contour is Γ 2 and the outer contour is Γ 1 , separated by the metasurface thickness δ. In this section, the extension of the formulation in [10] is outlined. A clear understanding of [10] is required to understand this section. The variables and notations follow [10]. The 2D IEs for TM polarization given in equation (6) of [10] are modified to whereρ ′ 1 andρ ′ 2 denote contour integration position variables along Γ 1 and Γ 2 respectively. Such a distinction was not made in [10]. The subscripts 1, 2 on the field variables denote the fields on domain Ω 1 , Ω 2 respectively. G 0 (ρ; ρ is the 2D Green's function. Derivative with respect to n ′ denotes directional derivative along the radial direction (pointing outwards). The 2D IEs for TE polarization given in equation (7) of [10] are modified to Following the same procedure as in [10], four integral equations are obtained using (1), (2) and circular metasurface synthesis equations given in equation (4) of [10]. These four equations are used to solve for the four tangential fields ) and H z1 (ρ ′ 1 ) on the exterior contour Γ 1 . The equations are given as Since the material properties of the metasurface vary azimuthally, the bianisotropic susceptibility tensor components and matrix A(ρ ′ c ) also varies with position. To be precise the matrix A(ρ ′ c ) is defined as The matrices A 1 (ρ ′ c ) and A 2 (ρ ′ c ), which are defined in equations (22), (23) in [10], are repeated in equations (6) and (7) for completeness.
The integration contours in (3a) -(3d) should be carefully noted. They are different from the corresponding equations in [10]. Similar to [10], the IEs (3a) -(3d) are converted to an MoM system of linear equations using pulse basis functions and point matching. The MoM system of linear equations reads     where the unknown vectors are N is the total number of basis functions. The Z matrix coefficients and the excitation vector are as follows where A ij,n is the value of matrix elements A ij at the n th discretization segment. It is these matrix elements which connect IE-GSTC to susceptibility tensor components. The coefficients p 1 mn , q 2 mn , r 1 mn , r 2 mn , s 1 mn , s 2 mn are calculated as s n1 , s n2 denote the n th discretization segments on contours Γ 1 , Γ 2 , respectively.ρ m1 ,ρ m2 are the middle points of the m th segments on contours Γ 1 , Γ 2 . δ mn is the Kronecker delta function. The superscript 1 or 2 in p 1 mn , q 2 mn , r 1 mn , r 2 mn , s 1 mn , s 2 mn denotes whether the integration is performed along the outer ) or inner contour. The excitation vector components are given by These four excitation vectors provide four possibilites: b 1,m is for a TM excitation external to the metasurface, b 2,m is for a TM excitation internal to the metasurface, b 3,m is for a TE excitation external to the metasurface, and b 4,m is for a TE excitation internal to the metasurface.

III. SUSCEPTIBILITY EXTRACTION OF CYLINDRICAL, MULTILAYERED MAGNETO-DIELECTRIC SHELL
In this section, the bianisotropic susceptibility tensor components of a cylindrical, multilayered, magneto-dielectric shell are extracted. The geometry is shown in Fig. 2. Even though a three layered geometry is shown, the method can be used for any number of layers. The medium inside and outside the shell is vacuum. TM z polarization is assumed. Follow-  (4) in [10], the differential field components across the cylindrical metasurface ∆E z , ∆H ϕ , ∆H z , ∆E ϕ are related to the average of fields across the metasurface E z,av , H ϕ,av , E ϕ,av , H z,av , through the following metasurface synthesis expressions: For the case of TM z polarization, they reduce to Since there are four unknown susceptibility components in the above equation, two different excitation problems are solved to obtain the four susceptibility components. The two excitations used are shown in Fig. 3. The first problem is that of an infinite, electric line source at the origin and the second excitation problem is that of a cylindrical surface electric current density of radius a s , where a s > a 4 . The two scattering problems are simple boundary value problems, which can be either solved analytically [24] or simulated using COMSOL. In the case of analytical solution, the vector magnetic potential, electric field and magnetic field in each region can be expressed in terms of Bessel functions: where J 0 (·), Y 0 (·) are Bessel functions of the first and second kind, respectively. k i , µ i are the wavenumber and permeability of region i. For the case of a cylindrical surface current density excitation of the three layered, magneto-dielectric cylinder, the total number of regions is 6 (the surface current density discontinuity splits the external medium into two separate regions). Applying the continuity of electric and magnetic fields across the region boundaries and the magnetic field discontinuity caused by surface current density, the unknowns C i , D i can be obtained. Details of this problem can be found in appendix A. Once the two excitation problems are solved, the susceptibility tensor components are obtained as follows where the superscript (1) denotes line source excitation fields and (2) denotes the cylindrical surface current density excitation. It should be noted that the magneto-dielectric shells and excitations considered are azimuthally invariant, hence the resultant tangential fields and the extractated susceptibilities have no dependance along the ϕ direction.

IV. NUMERICAL VALIDATION
This section presents examples of cylindrical metasurfaces which use IE-GSTC and the susceptibility extraction method described in the the previous sections. COMSOL software is used for full wave simulations.

A. Four segment cylindrical metasurface excited by plane wave
The cylindrical metasurface shown in Fig. 4 consists of four segments (quadrants). Each segment is a multilayered, concentric, magneto-dielectric shell. The susceptibility tensor components of each segment can be extracted using the method in the previous section. Referring to Fig. 4, different annular regions are identified as follows: region 1: 0 ≤ ρ < a 1 , region 2: a 1 ≤ ρ < a 2 , region 3: a 2 ≤ ρ < a 3 , region 4: a 3 ≤ ρ < a 4 and region 5: a 4 ≤ ρ < ∞. The metasurface is surrounded by vacuum, i.e. region 1 and region 5 are vacuum. The material properties of each region (represented by subscript number) in the four quadrants are given in table I. The susceptibility tensors of the quadrants are also given. The dimensions are: a 1 = 5λ, a 2 = 5λ + λ 24 , a 3 = a 2 + λ 24 , a 4 = a 3 + λ 20 . The susceptibility component extraction from the last section is performed four times to get the component values. In the IE-GSTC code, the matrix A is varied depending on the position. For instance, for the first quadrant, A is calculated using tensor components in the first row of table I.
The structure is excited by a plane waveĒ inc (x) = e −jk0xẑ . The total field along a line from (0, 0) to (10λ, 0) is plotted in Fig. 5 and Fig. 6. Good agreement between

Quadrant
Material Susceptibility tensor properties χ ϕz me = j1.45598 × 10 −3 ϵ r2 = 4, µ r2 = 1 χ ϕϕ mm = 9.25367 × 10 −3 ϵ r3 = 3, µ r3 = 1 χ zz ee = 2.80732 × 10 −2 ϵ r4 = 6, µ r4 = 2 χ zϕ em = −j9.74236 × 10 −4 ϵ r5 = 1, µ r5 = 1 Quadrant 4 ϵ r1 = 1, µ r1 = 1 χ ϕz me = j2.83566 × 10 −4 ϵ r2 = 1, µ r2 = 1 χ ϕϕ mm = 6.01626 × 10 −3 ϵ r3 = 6, µ r3 = 2 χ zz ee = 1.44824 × 10 −2 ϵ r4 = 1, µ r4 = 1 χ zϕ em = −j3.78786 × 10 −5 ϵ r5 = 1, µ r5 = 1 IE-GSTC and full wave results obtained using COMSOL software can be seen. For the IE-GSTC result a small gap can be observed around the metasurface radius of x = 5λ. This is because IE results is valid only for regions interior and exterior to the metasurface. Minor discrepancies between the two results can be seen around the metasurface. This is probably due to MoM being inefficient in the near field of the scatterer as well as GSTC being an approximate interface condition. In Fig. 5 and Fig. 6, the results from the zerothickness IE-GSTC [10] are plotted for the same problem to show that the formulation in this paper is required to get more accurate results. COMSOL discretizes and solves the fields inside the thin dielectric layers, while IE-GSTC abstracts it. In table I, it can be seen that even though the metasurface elements are passive, lossless and reciprocal the condition χ ϕz me = −χ zϕ em as stated in [14] is not satisfied. In [14], fig. 4 shows a unit cell characterization system in which a curved, spherical metasurface particle is illuminated by both inward propagating and outward propagating spherical waves. However, [14] provides no details on the numerical implementation of such a system. The susceptibility extraction technique in this paper is essentially an analytical, cylindrical version of the numerical, spherical system shown in fig. 4 in [14]. The condition χ ϕz me = −χ zϕ em is not satisfied because of the curvature, that is the inner contour of the scattering particle is not of the same length as of the outer contour. The inner contour length becomes equal to the outer contour length only when the radius of curvature becomes very large. This can be proved numerically with the susceptibility extraction method described in the last section. In tensor components of the metasurface element in quadrant I are shown as the inner radius a 1 is varied. It can be seen that for large radius of curvature, the condition is satisfied. The results of the susceptibility extraction are robust with respect to different illuminations. For instance, changing the value of outer surface current excitation i.e. a s gives the same results. Furthermore, a one dimensional cylindrical FDTD code [25] was developed to independently obtain the results. A one dimensional cylindrical FDTD with variation along the ρ direction and no variation along ϕ, z directions will automatically take the curvature into account.

B. Cylindrical metasurface with opening excited by a plane wave
This example, as shown in Fig. 7, is a four segment metasurface but with an opening. The opening subtends an angle of π/2 at the center and is equally divided between the first and fourth quadrants. The material properties of the quadrants are same as that of the previous example. GSTC cannot handle segments with material properties same as that of the background medium. In the IE-GSTC approach, an opening is handled by setting A(ρ ′ c ) = diag(1, 1, 1, 1) [10] . This ensures that tangential fields on either side of the opening are equal. This is only an approximation because the opening has a width of fraction of a wavelength. IE-GSTC results and COMSOL results on a circle of radius 8λ are plotted in Fig.  8 and Fig. 9.

C. Eight segment metasurface excited by a line source
In this example, an eight segment cylindrical metasurface with each segment subtending an angle of π/4 is excited by an infinite electric, line source at origin, i.e.J(r) = δ(x)δ(y)ẑ Am −2 . The dimension of the structure is same as that of the example A. The segment angular positions are [0, π/4], [π/4, π/2] and so on. The segments on the top half are exactly the same ones as in table I. The segments on the bottom half is listed in table. III. The magnitude and phase of the electric field on a circle of radius 8λ are shown in Fig. 10 and Fig. 11 respectively. Zero thickness IE-GSTC [10] results are also superimposed. In the COMSOL full wave simulation, the line source is modeled using the external current density option on a circle of radius much smaller than wavelength (∼ λ/100).

V. IE-GSTC MODELING OF CYLINDRICAL METASURFACE -PEC COMPOSITE STRUCTURES
In this section, the IE-GSTC formulation is extended to include PEC segments. GSTC can be used to model dielectric segments, dielectric segments with metallic/non-metallic inclusions, etc. When one of the segments is PEC, the IE-GSTC model needs to be extended. Potential applications include PEC backed metasurface lens antennas, intelligent surface enabled wireless environments [26], and scattering reduction and radiation problems [27].
The problem under consideration is shown in Fig. 12. Similar to the previous examples, there is no variation of material properties or excitation along the z axis. The inner contour is Γ 2 and the outer contour is Γ 1 , separated by metasurface thickness δ. The circular boundary which separates the external domain Ω 1 and internal domain Ω 2 is partly GSTC (mathematical abstraction for physical metasurface) and
where N G , N P are the number of basis functions for the GSTC and PEC parts, respectively, on either the inner or outer surface. Therefore the total number of unknowns are 2(N G + N P ). By using the above expansions and testing functions δ(ρ 1 −ρ 1G,m ), m = 1, 2, · · · , N G for (22a), δ(ρ 1 − ρ 1P,m ), m = 1, 2, · · · , N P for (22b), δ(ρ 2 −ρ 2G,m ), m = 1, 2, · · · , N G for (23a) and δ(ρ 2 −ρ 2P,m ), m = 1, 2, · · · , N P for (23b) an MoM system of equations can be derived. The MoM system of linear equations reads where the unknown vectors are The Z matrix coefficients and the excitation vector are as follows where A ij,n is the value of matrix element A ij at the n th discretization segment. The entries in the matrix elements are calculated as follows: where s n1 , s n2 denote the n th discretization segments on contours Γ 1 , Γ 2 respectively.ρ m1 ,ρ m2 are the middle points of the m th segments on contours Γ 1 and Γ 2 , respectively. δ mn is the Kronecker delta function. The first part of superscript in the above coefficients denotes whether the testing is performed on the inner or outer contour. It also shows whether the testing is done on the GSTC part or PEC part. The second part of the superscript denotes whether the integration is performed in the outer or inner contour's GSTC part or PEC part. For example in t 1P1G mn , the testing is done on the exterior contour's PEC part (i.e. 1P) and the integration is performed on the exterior contour's GSTC part (i.e. 1G).
The excitation vector components are given by These four excitation vectors provide the possibility for TM z excitation internal and external to the metasurface-PEC system. To validate the formulation, consider the metasurface-PEC system example shown in Fig. 13. The four metasurface segments are identical to what is shown in table I. The first quadrant in table I corresponds to a metasurface segment in the range 0 ≤ ϕ ≤ π/4 in Fig. 13 and so on. The dimensions are same as used in previous examples. The excitation is a plane waveĒ inc (x) = e −jk0xẑ . In Fig. 14 and Fig. 15, IE-GSTC-PEC results are compared with COMSOL full wave simulations on a circle of radius 10λ. VI. CONCLUSION This work validates IE-GSTC for circular cylindrical metasurfaces of nonzero thickness. The essence of this work is that physical circular metasurface scattering particles can be abstracted by a few numbers, i.e., susceptibility tensor components. To the author's knowledge, this work is one of the few works in which the GSTC and a bianisotropic susceptibility tensor model is validated for a physical, conformal metasurface. Although the examples and susceptibility extraction are presented for the TM polarization, it is easily extendable to the TE polarization. The combination can be used to analyze and design circular metasurfaces for the case of circularly polarized incident waves. This work is important in the sense that it pushes the previous works on conformal GSTC [10]- [12], [14] to physical implementation and optimization. The IE-GSTC formulation is extended to handle PEC segments which could be useful for   fig. 3. The excitation is a cylindrical surface electric current density of radius a s , represented as J s = δ(ρ − a s )ẑ. Even though fig. 3 shows a three layered shell, this section generalizes it to m − 1 magneto-dielectric layers. The m + 2 regions, i.e. m − 1 magneto-dielectric regions, region internal to the shell, and space outside the shell separated into two regions by the surface current density are identified as follows: Region i : a i−1 ≤ ρ ≤ a i ; i = 1, 2, · · · , m + 2 a 0 = 0 ; a m+1 = a s ; a m+2 = ∞ where a m+1 = a s separates the external space (outside the magneto-dielectric shell) into two regions. There is a field discontinuity between these regions due toJ s = δ(ρ − a s )ẑ. Region i has the constitutive parameters (ϵ i , µ i ) and wave number k i = ω √ ϵ i µ i . Assuming no field or material variation along ϕ and z direction, the vector magnetic potential, electric field and magnetic field in each region for the case of TM z polarization can be expressed in terms of Bessel functions.
where J 0 (·), Y 0 (·) are Bessel functions of the first and second kind, respectively. Hence there are 2(m + 2) unknowns to be determined. The singularity of Y 0 (·) at ρ = 0 and the presence of only an outgoing wave in the outermost region m + 2, i.e. A z,m+2 (ρ) ∝ H 0 (k m+2 ρ) entails the following conditions: The continuity of the tangential electric field across the region interfaces results in: The continuity of the tangential magnetic field across the interfaces and the discontinuity caused by the surface electric current density excitation results in: Y 1 (k i+1 a i ) = δ i,m+1 ; i = 1, 2, · · · , m + 1 The set of equations (35), (36) and (37) constitute 2(m + 2) equations that can be solved for fields everywhere.