Surface Susceptibilities as Compact Full-Wave Simulation Models of Fully-Reflective Volumetric Metasurfaces

While metasurfaces (MSs) are constructed from deeply-subwavelength unit cells, they are generally electrically large and full-wave simulations of the complete structure are computationally expensive. Thus, to reduce this high computational cost, non-uniform MSs can be modelled as zero-thickness boundaries, with sheets of electric and magnetic polarizations related to the fields by surface susceptibilities and the generalized sheet transition conditions (GSTCs). While these two-sided boundary conditions have been extensively studied for single sheets of resonant particles, it has not been shown if they can correctly model structures where the two sides are electrically isolated, such as a fully-reflective surface. In particular, we consider in this work whether the fields scattered from a fully reflective metasurface can be correctly predicted for arbitrary field illuminations, with the source placed on either side of the surface. In the process, we also show the mapping of a PEC sheet with a dielectric cover layer to bi-anisotropic susceptibilities. Finally, we demonstrate the use of the susceptibilities as compact models for use in various simulation techniques, with an illustrative example of a parabolic reflector, for which the scattered fields are correctly computed using a integral equation (IE) based solver.


I. INTRODUCTION
In the past 20 years, a variety of approaches have been applied to model the behaviour of electromagnetic metamaterials, each with trade-offs in complexity, physical insight, and computational burden.MSs present an inherently "multiscale" modelling problem: on one hand, they are composed of sub-wavelength scattering elements, which produce strong variations of the fields at the microscopic level, while on the other hand the complete MSs are generally electrically large.Thus, the most rigorous method-full-wave numerical simulations-is computationally expensive.For a uniform metasurface, the application of periodic boundaries reduces the computation region to a single unit cell, but this is not possible for a non-uniform surface where the complete structure must be modelled.Doing so provides the fields at the microscopic level-which may provide physical insight-but it is not efficient for an iterative design flow due to the large simulation model.
For this reason, equivalent models have generally been used for design, which approximate the structure as a zero-thickness boundary [1].One possibility is the impedance boundary conditions (IBCs), which model the metasurface as a sheet of electric and magnetic currents, related to the tangential fields [2]- [5].These can be useful for multi-layer structures such as stacks of metallic patterned layers between dielectric layers, since they can be cascaded akin to elements on a transmission line [6].However, a disadvantage of the IBCs is that the impedances depend on the angle of incidence.Thus, there is no unique set of impedances which truly characterizes the MS, independent of the angle of incidence.Furthermore, IBCs do not take into account normal polarization currents that may be induced; e.g., in planar split ring resonators.
Alternatively, surface susceptibilities (χ) can be used to represent the surface, in conjunction to the generalized sheet transition conditions (GSTCs) which provide boundary conditions on the tangential (and normal) components of the fields on either side of the surface [7]- [11].The susceptibilities relate the acting fields to the induced electric and magnetic polarization densities, where the acting fields have generally been defined one of the two following ways.One possibility is the Tretyakov-Simovski (TS) model, where the acting fields are the incident fields.The second is the Holloway-Kuester (HK) model, which defines the acting fields as the average of the total fields on either side.The latter has an advantage in that resonances are easier to identify in the constitutive parameters (susceptibilities) [12] and so we use the HK approach in this work.These provide true constitutive parameters that can predict the scattered fields, regardless of the incident field 1 , provided that the susceptibilities are correctly selected and extracted.Consequently, the GSTCs have been implemented in a number of simulation techniques, such as finite difference methods [13], [14], the finite element method [15], and integral-equation methods [16], which provide efficient simulations of electrically large-and possibly curvilinear-MSs and their coupling to other scattering objects.
However, the GSTCs as they are generally written as in [7]- [11], have only been shown to rigorously model structures which are composed of a single layer.Generally, these are metasurfaces composed of an array of resonators [2], [7], while the mapping for a single dielectric layer has also been shown [17]- [19].On the other hand, one may wonder: is it possible to apply the GSTCs and the HK model to model a cascaded structure, such as dielectric layer, backed by a PEC ground-plane?While this has been considered using singlesided boundary conditions, which involve the fields which interact with the slab and reflect from the ground-plane [20]- [22], in this paper, we consider whether the two-sided GSTCs can be used to model such a structure.The model should have asymmetrical behaviour, behaving as a PEC with a cover layer on one side, and a PEC on the other side.The two sides are independent (electrically isolated), so it is not immediately clear that the HK model-which uses the average of the total fields as the acting fields-should work.Surprisingly, we will show that it is indeed possible with some limitations.
The paper is structured as follows.In Section II, we provide the GSTCs and simplify the susceptibility tensors for the problem at hand; these are used to derive expressions for the S-parameters of a uniform surface.Next, Section III gradually builds susceptibility models in increasing complexity: a PEC sheet, dielectric slab, a dielectric slab with a ground-plane, and finally the aforementioned reflective metasurface.In doing so, we use an accessible approach with notation for the GSTCs and susceptibilities that have been used in recent literature, and show the limitations of the models.In Section IV, a reflective unit cell is used to design and simulate a parabolic reflector, and a comparison is made to a full-wave simulation to show the accuracy of the susceptibility model for nonuniform metasurfaces.Finally, we conclude in Section V.

II. GENERALIZED SHEET TRANSITION CONDITIONS (GSTCS) AND SURFACE SUSCEPTIBILITIES
We will use the formulation for surface susceptibilities as presented in [10].Briefly, the boundary conditions at the metasurface, called the generalized surface transfer conditions (GSTCs), are with ∆φ = φ t − (φ i + φ r ) being the difference in fields across the boundary (φ ∈ {E, H}), and we define n = ±ẑ as being the surface normal in the direction of incidence, directed from the side on which the incident field is present to the transmission side.P and M are the electric and magnetic surface polarization densities, respectively, with denoting the projection to the boundary, while z is the normal part.
Meanwhile, the polarization densities are related to the averaged electric fields by the constitutive relations (HK model) with φ av = 1 2 (φ t + φ i + φ r ).There are four sets of tensors, χ, for a total of 36 constitutive parameters.Many of these components can be eliminated or simplified due to reciprocity, symmetry, or energy conservation, depending on the particular surface [11], [23].In this work, we will consider surfaces With the 8 unique terms retained in these selected tensors, no assumptions have been made regarding energy conservation, and there is a possibility for omega-type bianisotropy [24].Since the susceptibilities in (3) do not convert polarization, TE and TM illuminations can be considered separately, as depicted in Figure 1.With the periodicity of the surface being subwavelength, no higher-order diffraction orders are generated [25] and under oblique plane wave illumination at θ, there will be reflected and transmitted plane waves at the same angle, following standard Snell's laws.The "ports" on the left and right sides are denoted 1 and 2, respectively, such that e.g. S {TE, TM} 21 denotes transmission in the forwards direction and S {TE, TM} 12 in the backwards direction.Substituting the expressions for the plane waves along with (3) into ( 1) and ( 2), and solving the resulting system of equations for the Sparameters, we obtain where the top and bottom signs (∓) are taken for 11/21 and 22/12, respectively, and We note in particular that the bianisotropic terms χ {yx,xy} em [and implicitly terms in χ me following (3b)] lead to an asymmetry in reflection when the direction of illumination is flipped, as noted in [11].{12,21} (−θ) implies that the metasurface leads to angular symmetric scattering [11].

III. SURFACE SUSCEPTIBILITY MODELS
We will use (4) to determine the susceptibility models for several simple structures, progressively building up to a reflective metasurface, for which the susceptibilities correctly model the behaviour with illumination from either side.We use a pedagogical approach, starting with simple structures such as a PEC sheet and a dielectric sheet to motivate the susceptibility terms in (3) and understand their role for the final reflective metasurface.= 0.A PEC is the limiting case of a conductor with infinite conductivity (σ → ∞), which in fact, corresponds to a limiting case of χ xx ee and χ yy ee .By eliminating all other susceptibility terms from (4) and taking a limit, we find3 lim χ yy ee →−j∞

A. PEC Sheet
Thus, χ {xx,yy} ee → −j∞ precisely models a PEC sheet.4However, this is not suitable for numerical simulation, since it requires a limit.
→ −j∞ would be a PMC.Now, we consider if there is another susceptibility model, having finite susceptibilities, which could be used instead.To this end, we return to considering all 8 terms in (3).Of course, including the bianisotropic terms produces asymmetry with forwards/backwards illumination, and so here we seek a possible compromise that may not rigorously model a PEC but is useful for numerical implementation.We will enforce S {TE,TM} {12,21} (θ) = 0 and S {TE,TM} 11 (θ) = −1 (PEC with forwards illumination); these provide 6 equations from (4) while S {TE,TM} 22 (θ) is not enforced.To reduce the number of unknown terms from 8 to match the number of equations, we consider symmetry.There is rotational symmetry around the z axis, i.e. isotropy, and so the susceptibilities should remain unchanged with this rotation.That is should hold true, where R z (φ) is the transformation matrix which rotates by φ about the z axis [23].This requirement means χ xx {ee,mm} = χ yy {ee,mm} leaving 6 unknown susceptibilities.Solving the resulting system of equations yields the non-zero terms χ yx em = +2j/k and χ xy em = −2j/k.Note that these terms substituted into the bianisotropic tensors (3b) also satisfy (6).
Now, what happens for backwards illumination?Using (4a,c) we find S {TE,TM} 22 (θ) = +1.Thus, the model appears as a PEC with forwards illumination and a PMC with backwards illumination, independent of the angle of incidence.If χ {xy,yx} em are negated, the structure is effectively reflected in the x−y plane.The angular independence with these purely bianisotropic susceptibilities was studied in [26], where it was also noted that such a surface is only possible to fabricate with a physical unit cell at a single frequency, in the limit of zero loss.For this work, however, the inability to physically represent a PEC sheet (with correct behaviour on both sides) does not end up being critical in our our subsequent objective, when we later consider a dielectric cover layer added to the PMC side.

B. Dielectric Sheet
Building towards a resonator on top of a ground-plane with a cover layer, we next consider an isolated sheet of uniform permittivity r = ( r −j r ), having an imaginary part allowing for loss, with a thickness d.This mapping was considered in [19], where only the tangential {ee, mm} susceptibilities were used and earlier in [18] where the normal components (χ zz {ee,mm} ) were included.We will consider the extraction here for completeness, using a simple and accessible approach.
The analytical reflection and transmission through a dielectric slab, are well-known, given by with a ∈ {TE, TM}, φ = kd √ r cos θ, and where ρ a {1,2} (θ) and τ a {1,2} (θ) are the Fresnel coefficients for oblique incidence at the first and second interfaces, dependent on the polarization (TE/TM) [27].The matrix equation can be solved to provide a total of four expressions for S {TE,TM} {11,21} (θ).With θ = 0°, the normal susceptibility terms in (4f) are eliminated, and we set χ {yx,xy} em = 0 so that there is symmetry between forwards and backwards illumination, we can substitute S {TE,TM} {11,21} (0) from ( 7) into (4) to have four equations with the solution Repeating a similar procedure for θ = 0°, the normal susceptibility terms are found to be with γ = r − sin 2 θ.As is, there is an angular dependence in (8c-d), which should not be the case, if χ zz {ee,mm} are to be true constitutive parameters, independent of the incident field [11].However, we can only expect a zero-thickness model to apply for thin slabs (kd 1), and so we proceed by using the first few terms of the Taylor expansion around kd = 0: with the approximation that the third term and higher-order terms are negligible.Within this approximation, χ zz {ee,mm} only depends on the properties of the slab (at a given frequency), and can be considered constitutive parameters.
To validate the assumption and verify ( 8) and ( 9), we consider a numerical example, shown in Figure 2. The reflection and scattering for both TE and TM illuminations is shown, for increasing thickness kd and at three angles, θ ∈ {0°, 30°, 60°}.Firstly, we see that even when the slab is very thin (kd = 0.2, i.e. d ≈ λ/30) the normal components χ zz {ee,mm} are needed to model the scattering at oblique incidence.Secondly, the modelling is very accurate up to about kd = 0.8 (d ≈ λ/8), past which the analytical results diverge.

C. PEC Sheet with Cover Layer
Next, consider a PEC with a dielectric cover layer on one side.In this case, the forwards and backwards illuminations clearly have non-symmetrical behaviour.With the cover layer χ , ( 8) & ( 9) χ with χ zz = 0 Fig. 2. The reflection and transmission through a lossy dielectric sheet ( r = 4 − j0.04) was calculated analytically using (7) to compare with the mapped susceptibilities ( 8), with and without the normal components.on the left side of the PEC sheet in Figure 1, the forwards illumination case is governed by where E a PEC is the amplitude of the forwards-travelling wave at the PEC.Since E a PEC is not of interest, this reduces to two expressions for S after all but the leading angle-independent terms in the Taylor expansion are retained for (11e,f) (the complete expressions for χ zz {ee,mm} are cumbersome and not presented).It is interesting that χ {yx,xy} em are in fact identical to those for a sheet that behaves as a PEC on the back side (S {TE, TM} 22 = −1) and a PMC on the front side (S {TE,TM} 11 = +1), as discussed in Section III-A.However, in this case the other susceptibilities contribute such that the front does not appear as a PMC, but a PEC with a cover layer.
To check that this is is indeed the case, numerical results are shown in Figure 3.For TE polarization, (11) provides a good match for the forwards reflection past kd > 1, in comparison to the analytical expression (10); the resonance at kd ≈ 0.9 is correctly predicted.For TM polarization, there is reasonable agreement if χ zz ee is neglected, but including causes a large disagreement.Thus, it seems that the truncation of the Taylor series in (11e) is a poor approximation, and this term should not be used.Meanwhile, the model predicts S {TE,TM} {21,12} = 0 and S
Finally, note that ( 11) satisfies (6).That is, the susceptibility tensors are isotropic-as should be the case for this rotationally symmetric structure-and thus are valid for any polarization or plane of incidence (or any incident field, in general).

D. Sub-wavelength Resonator on a Dielectric Slab
Next, we extend the structure by including an array of electric dipoles on top of the grounded slab; this can truly be called a metasurface.For simplicity, we will only consider TE polarization, since the unit cell was designed to resonate with TE polarization.Shown in Figure 4, the unit cell is deeply sub-wavelength (λ/10 at 30 GHz).It has a "dogbone"shaped copper dipole loaded with a lumped inductance (L) at the center, placed on a Rogers RO4003C substrate (508 µm thickness; i.e. kd = 0.32), on top of a PEC ground-plane.
While there is no analytical expression for S TE 11 (θ), it can be found through full-wave simulations of the unit cell [2], [10].In particular, two angles of incidence need to be simulated.Then, using (4) to solve the susceptibilities 5 , with all other terms being irrelevant for TE polarization due to the lack of polarization conversion.This extraction was performed using full-wave simulations (Ansys HFSS) of the unit cell with periodic boundaries, using θ ∈ {0°, 60°}, and with 1.2 nH ≤ L ≤ 1.45 nH; the susceptibilities are plotted in Figure 4b.Both χ yy ee and χ zz mm are extracted, while χ yx em only depends on frequency through (12b) (not plotted).Subsequently, the predicted reflection is plotted in Figure 4c and compared to the full-wave simulated one for L = 1.4 nH, at which there is a resonance at 30 GHz.We observe nearly perfect agreement with the HFSS simulations for S TE  11 at the two different angles of incidence, while we also have S TE  22 = −1 corresponding to an ideal PEC at all frequencies (not shown here) thanks to the inclusion of χ yx em .Thus, this combination of χ yy ee , χ yx em , and χ zz mm produces a very accurate model for scattered TE-polarized fields, regardless of the incident field.
IV. NUMERICAL DEMONTRATION: PARABOLIC REFLECTOR Finally, we demonstrate the utility of the two-sided susceptibility model of the reflective unit cell from Section III-D by designing a parabolic reflector; i.e. a metasurface that reflects a plane wave when a cylindrical wave is present from a line source (forward illumination).The incident field produced by the line source is where r s = (x s , z s ) is the location of the source, H {0,1} are Hankel functions of the second kind, of order 1 and 2, and E 0 is the amplitude.For our example, r s = (0 cm, −10 cm), at 30 GHz.Meanwhile, the reflected field should have a uniform phase, while we taper the amplitude so that the surface remains passive.Specifically, E r,y (x, 0) = 0.9 |E i,y (x, 0)| (14a) Using ( 13) and ( 14), we calculate the ideal χ yy ee point-wise along x using (12c), treating χ zz mm as a perturbation. 6We do not aim to optimize the design, but rather wish to show how the two-sided susceptibility model is accurate for a finite-sized and non-uniform surface.This produces the "desired" profile in Figure 5b.However, the unit cell extraction is limited in replicating this profile since there is only one parameter L, so that the real and imaginary parts of χ yy ee cannot be tuned independently.For the sake of this demonstration, the real part is realized, using the look-up plot from Figure 4b.Meanwhile, the (χ yy ee ) (and χ zz mm ) are also included in the susceptibility model.As seen in Figure 5b, the imaginary part deviates from the desired profile.Finally, the profile is discretized into 50 cells, for a length of 50 mm.
To check the accuracy of the susceptibility model for this finite non-uniform MS, we use an integral-equation simulator which implements the complete dipolar GSTCs [16].Meanwhile, the structure is simulated in HFSS using the model in Figure 5a, where we leverage symmetry in the y − z plane to reduce the size of the model by half (25 unit cells are simulated).The magnitude of the total fields are compared in Figure 5c, for a cylindrical wave originating from r s = (0 cm, −10 cm), i.e., forwards illumination.We see a good match between the predicted total fields based on the non-uniform susceptibilities and the full-wave simulation, which takes orders of magnitude longer in time to simulate.Moreover, we see that the susceptibility model correctly models the diffraction around the edges of the finite surface, which is where "backward illumination" behavior is critical.The reflected fields at the z = −10 mm plane are plotted for closer examination in Figure 5d.While an ideal parabolic reflector would produce a uniform amplitude and phase, there are fluctuations, which may be explained by the finite size of the reflector and edge diffraction, and the fact that the unit cell had only a single parameter of control so that only (χ yy ee ) was implemented.However, it is not our purpose to optimize the design, which would require a more complex unit cell with loss, but rather to show the accuracy of the susceptibility model.To this end, we see an excellent match between BEM and HFSS.
Lastly to exemplify the two-sided GSTCs, in Figure 5e we repeat the simulations with a cylindrical wave originating from r s = (0 cm, 10 cm), i.e., backward illumination.As expected, we observe a good match for the total fields, but the interference obscures the difference between the forwards and backwards illumination cases.By plotting the reflected fields alone (Figure 5e), we see that there is a curved phase, which is to be expected for a PEC sheet with an incident cylindrical wave.Thus, the PEC behavior with backwards illumination is correctly modeled and confirmed.

V. CONCLUSIONS
We have shown that it is possible to use the conventional two-sided GSTCs and the HK model for the surface susceptibilities to predict the fields scattered from a fullyreflective metasurface.The validity of the HK model may not be a priori obvious, since this model defines the acting fields as the average of the fields on either side of the metasurface, which seems to be in contradiction with the considered cases where the two sides are in fact electrically isolated and independent.Nevertheless, this approach proved to be in good agreement with full-wave simulations.Indeed, the retrieved susceptibilities from the model work for both forwards and backwards illumination, behaving as a PEC for the latter, which is important if the surface is finite such that the fields may interact with the reverse side.We have also shown the mapping of the geometrical and electrical properties of a dielectric slab (and a PEC sheet with a dielectric cover layer) to susceptibilities, showing the role of the normal and bi-anisotropic terms.These equivalent surface susceptibilities thus act as compact models for these practically volumetric structures and may be easily integrated into a variety of simulation platforms to enable an efficient computation of scattered field from finite-sized volumetric structure for a faster iterative design flow.

Fig. 1 .
Fig. 1.A depiction of the electric field orientations for TE and TM fields, with forwards and backwards plane wave illumination of a planar structure (e.g. a metasurface).Propagation is in the x − z plane, with the metasurface in the x − y plane.

Fig. 4 .
Fig. 4. A deeply-subwavelength reflective unit cell was designed using an electric dipole, loaded with a lumped inductor (L) on top of a Rogers RO4003C substrate ( r = 3.55, tan δ d = 0.0027) that is on a PEC ground-plane.(a) Unit cell model.(b) χ extraction at 30 GHz using (12), as a function of the loading inductance L. (c) S TE 11 with L = 1.4 nH at two angles of incidence.

Fig. 5 .
Fig. 5.A parabolic reflector was designed using the unit cell from Figure 4, having a focal length of 10 cm with forwards illumination.(a) HFSS simulation model of half of the structure, with symmetry in the y − z plane, and a total of 2 × 25 cells; i.e., total length of 5 cm.(b) Desired and realized susceptibilities for reflecting a plane wave.(c) Comparison of the total field magnitude, using a boundary element method (BEM) which implements the GSTCs and HFSS (full-wave), with an incident cylindrical wave originating from rs = (0 cm, −10 cm).(d) Comparison of reflected fields at z = −10 mm.(e-f) Same as (c-d) but with an incident cylindrical wave originating from rs = (0 cm, +10 cm)