Part 1: Spatially Dispersive Metasurfaces: Zero Thickness Surface Susceptibilities&Extended GSTCs

A simple method to describe spatially dispersive metasurfaces is proposed where the angle-dependent surface susceptibilities are explicitly used to formulate the zero thickness sheet model of practical metasurface structures. It is shown that if the surface susceptibilities of a given metasurface are expressed as a ratio of two polynomials of tangential spatial frequencies, $\boldsymbol{k_{||}}$ with complex coefficients, they can be conveniently expressed as spatial derivatives of the difference and average fields around the metasurface in the space domain, leading to extended forms of the standard Generalized Sheet Transition Conditions (GSTCs) accounting for the spatial dispersion. Using two simple examples of a short electric dipole and an all-dielectric cylindrical puck unit cells, which exhibit purely tangential surface susceptibilities and reciprocal/symmetric transmission and reflection characteristics, the proposed concept is numerically confirmed in 2D. A single Lorentzian has been found to describe the spatio-temporal frequency behavior of a short dipole unit cell, while a multi-Lorentzian description is developed to capture the complex multiple angular resonances of the dielectric puck. For both cases, the appropriate spatial boundary conditions are derived.


I. INTRODUCTION
Electromagnetic Metasurfaces are 2D arrays of subwavelength resonating particles that derive their macroscopic field response from their geometry and material arrangements at the microscopic scale [1]. By engineering these resonating particles, a wide variety of macroscopic fields transformations may be achieved, which has led to a myriad of exotic applications across the electromagnetic spectrum, ranging from cloaking, illusions to holograms extended to real-time reconfiguration of wireless environments [2]- [8].
From their very nature, metasurfaces have a multi-scale architecture with sub-wavelength resonators arranged in an electrically large array. Moreover, computation of the scattered fields from these metasurfaces inherently needs to resolve the sub-wavelength features, which is generally a significant computational task -affecting both design synthesis and subsequent field analysis of electrically large metasurfaces. As a result, exploiting their electrically thin characteristics, they have been modeled as zero thickness surfaces, as spatial João G. Nizer Rahmeier, Tom J. Smy, Jordan Dugan, and Shulabh Gupta are with Carleton University, Ottawa, Canada (e-mail: JoaoNizer@cmail.carleton.ca). discontinuities, described in terms of dipolar tensorial electric and magnetic surface susceptibilities,χ(ω, r) [9], [10]. The Generalized Sheet Transition Conditions govern the resulting macroscopic fields (GSTCs) [11], [12] involving electric and magnetic surface polarization densities.
An essential step in building equivalent zero thickness models of practical metasurfaces is mapping geometrical/material characteristics to tensorial surface susceptibilities at specified design frequencies, accounting for a complete angular scattering of the surface beyond paraxial wave propagation. Most of the work in the literature has been focussed on one class of these structures: spatially non-dispersive metasurfaces, where the surface polarizations are induced due to local field interaction only. In such cases, the surface susceptibilities, being characteristic properties of the metasurface, are entirely independent of the angle of incidence of the incoming waves. As a result, the angular scattering behavior of the metasurfaces can thus be described in terms of the normal components of the surface polarizations, assuming they are physically supported by the structure [9], [13], [14].
Owing to a virtually unlimited number of metasurface structures proposed and demonstrated in the general area of electromagnetic metasurfaces, the foundational resonators structures range from deeply sub-wavelength sizes to close to wavelength dimensions. Of which a typical example, out of many, are all-dielectric resonator-based Huygens' structures, for instance [15]- [18]. Thus, a priori, metasurface structures are not necessarily spatially non-dispersive and may exhibit a non-local response, i.e., spatial dispersion [19]. Such structures, in general, cannot be modeled using standard dipolar surface susceptibility models and, in particular, the normal surface polarization alone.
Some recent works have explored surface susceptibility models of spatially dispersive metasurfaces, particularly in terms of multipolar description of the resonators and hypersurface susceptibilities [20]- [22]. However, these techniques have been typically demonstrated for weak spatial dispersion and have not replicated the complete angular scattering response of an arbitrary metasurface, which may even exhibit multiple resonances across the angular spectrum. Therefore, a general zero thickness modeling of spatially dispersive structures, particularly for strong dispersion, has been an open problem in the literature. In the present work, we seek to describe general metasurface structures in terms of a compact surface susceptibility model compatible with the GSTCs, which represents general spatial boundary conditions. This paper represents Part 1 of this work, where we propose a simple angle-dependent surface susceptibility modeling technique that forms an equivalent zero thickness sheet model of an arbitrary spatially dispersive metasurface structure. We specifically exploit angle-dependent surface susceptibilities and describe them as a ratio of polynomials of the transverse wave-vector k || . We further describe them in terms of physically motivated Lorentz oscillator model with angle-dependent resonator parameters, which leads to a very compact way to describe the complete angular scattering of the metasurfaces. Moreover, they feature convenient spatial domain forms with spatial derivatives of both the difference and average fields at the metasurface, resulting in extended GSTCs, which can be straightforwardly integrated into electromagnetic fields solvers as generalized boundary conditions, e.g., an Integral Equation (IE) solver, as presented in Part II of this work [23].
The paper is structured as follows. Sec. II reviews the conventional GSTCs and proposes the principle of modeling an arbitrary spatially dispersive metasurface structure based on angle-dependent surface susceptibilities. A general representation of angle-dependent surface susceptibilities as a ratio of two polynomials is proposed in Sec. III, which leads to the extended form of the GSTCs involving spatial derivatives of difference and average fields at the metasurface. Next, a Lorentz oscillator model is proposed, which forms a compact model to represent arbitrarily complex structures with multiple angular resonances. The proposed technique is then demonstrated numerically in Sec. IV for two structures: short dipoles and an all-dielectric resonator. Sec V further discusses a connection between spatial dispersion and the normal surface susceptibility components. Finally, conclusions are provided in Sec. VI, summarizing the work and describing the future steps.

A. Angle Dependent Surface Susceptibilities
Consider a field scattering problem from a metasurface which is lying in the y-z plane, and incident with arbitrary fields, thereby generating the scattered fields in the transmission and reflection regions, as shown in Fig. 1. The metasurface is assumed to be a zero thickness sheet, i.e., δ = 0. When excited with an incident plane-wave, ψ 0 (θ, ω), induced electric and magnetic surface currents, {J s , K s } are generated, which then re-radiate to produce the scattered fields in both reflection and transmission regions. For spatially nondispersive metasurfaces, these currents represent the surface response to the local incident fields only (i.e., point-by-point interaction). Thus, their zero thickness sheet model can be described using constant (independent of the angle of incidence, θ) electric and magnetic dipolar surface susceptibility tensors, χ. Hence, the corresponding dipolar surface polarizations are related to the averaged fields around the metasurface [24], as where each of the tensorsχ α,β are (3 × 3) matrices containing both tangential and normal susceptibility components. In addition, the electromagnetic fields at the metasurface follow the Generalized Sheet Transition Conditions (GSTCs) given by 1 , wheren is the surface normal (e.g. along x−axis) and ∆ψ = (ψ + − ψ − ) is the field difference across the metasurface [1], [11], [12]. However, for a spatially dispersive metasurface, the surface susceptibilities are functions of the incoming incidence angle (or spatial frequencies, k y ), as illustrated in Fig. 1. This is analogous to temporal dispersion where the constitutive parameters of a medium depend on the temporal frequencies, i.e. = (ω) and µ = µ(ω). It results in a re-arrangement of the instantaneous temporal frequencies when a broadband signal propagates through the medium, leading to envelope distortion in time. For spatially dispersive surfaces, the induced surface currents depend on θ, so that their constitutive parameters, χ α,β =χ α,β (k y ). This means that the induced currents at a location y on the surface depend on the local, as well as non-local fields across the metasurface so that a pointby-point interaction no longer holds. Moreover, the spatial frequencies of an incoming spatially broadband signal will be re-arranged in space, distorting the spatial field envelope. Therefore, an engineered spatially dispersive surface acts like a spatial frequency filter.
For simplicity, let us consider a symmetric structure with TE mode excitation, so thatχ em, me terms are zero, and we further assume that there are no normal components, so that the surface is completely described by tangential surface susceptibilities, i.e.,χ yy mm andχ zz ee [13], [14]. Using (1) and (2), we can compute the two surface susceptibilities as [25] χ ee (k y = k 0 sin θ) = 2j cos θ k 0 where R and T are the reflection and transmission response of the surface, as a function of the angle of incidence θ of an incoming uniform plane-wave, i.e. δ(k y − k 0 sin θ), at a specified temporal frequency ω. Under these conditions of non-local interaction, (1) does not hold, and a relationship between induced surface currents and the average fields, via the angle-dependent surface susceptibilities, must be modified to account for the non-local field interaction.
{Js, Ks} ky = k0 sin θ angular spectrum Non-dispersive χ α,β (ky) = const Fig. 1. The general field problem of a spatially dispersive metasurface where its surface susceptibilities,χ(θ) are dependent on the angle of incidence, θ of the incoming plane-waves. A uniform surface with purely tangential surface susceptibilities with no bi-anisotropic terms, excited with TE mode incidence fields are assumed throughout this work for simplicity.

B. Surface Polarizations
In order to avoid mathematical complexity, and focus on the underlying physical concept, consider a uniform metasurface, lying in the y − z plane at x = 0, and excited with an oblique incident uniform plane-wave, which induces surface currents J andK on the surface. The fields radiated by these currents are obtained using the electric and magnetic vector potentials A(r) and F(r) [26]. Since the currents are only on the surface and are fully tangential,J(y) =J s δ(y ), so that [20] The scattered fields radiated due to these currents may be obtained as The total fields on each side of the surface are then given by the sum of the incident and total scattered fields due to these currents. Assuming incident fields only on the reflection side of the metasurface (x < 0), the average fields are given bỹ For local field interaction, the incident fields are related to the surface currents at x = 0, as: where it is assumed thatᾱ α,β are angle-independent, and thus non-dispersive. The currents are thus induced due to local fields only following point-by-point interaction. However, for a general spatially-dispersive metasurface, with non-local field interaction, the surface fields are given by the spatial convolution of the average fields with the polarization response of the surface, i.e.
To avoid convolution operation in space, this can easily be expressed in the spatial frequency domain, i.e., F y {ψ(y)} = ψ(k y ), so that, where k y represents the tangential component of the wavevector k along the surface (i.e., related to the incidence angle of the incoming wave). IsolatingẼ i andH i , substituting them in the spatial Fourier transform of (6), and rearranging the terms, we relate the average fields to the surface currents, i.e.
Finally, expressing surface dipole moments, as surface current densities over the unit cell area S, asJ s = jωP andK s = jωM, and solving for the surface polarizations, we get whereχ αβ are the angle-dependent surface susceptibilities. For a non-spatially dispersive metasurface with local field interaction only,χ(k y ) = const., so that we retrieve the standard surface susceptibility relations of (1) following (7).

III. EXTENDED GENERALIZED SHEET TRANSITION CONDITIONS (GSTCS) A. Susceptibilities as Ratio of Partial Fractions
Let us consider a spatially dispersive metasurface that is characterized using tangential surface susceptibilities only, i.e. χ zz ee and χ yy mm . Using (11), for a TE mode, we get To model such a metasurface, we need to express their angular dependence using a convenient functional representation. We postulate that a general surface susceptibility function can be expressed as a ratio of two polynomials in k y , with known complex coefficients featuring various possible poles and zeros [27] accounting for angular resonances, so that we can writẽ Using (12), this can be expressed as: Taking the inverse spatial Fourier transform, the terms with various polynomial orders turn into spatial derivatives, resulting in Now, we know from the GSTCs that a portion of the field differences can be associated with this polarization, i.e.
Substituting this in the above equation results in an extended form of the GSTC to Therefore, for a metasurface that is explicitly described using tangential surface susceptibilities only, spatial dispersion manifests as spatial derivatives of both the difference fields and the average fields across the metasurface. This is the key result of this work. For a spatially non-dispersive cell, the above extended GSTC naturally reduces to their standard forms as with b 0 = d 0 = 1, a 0 = χ zz ee and c 0 = χ yy mm . It should be noted that for a perfectly symmetric metasurface invariant to the sign of ±k y , only even orders of the polynomial exist in the general representation of (13), i.e., only even derivatives in (17). The extended GSTCs of (17), now represent general boundary conditions which can now be operated on arbitrary incident fields, and be integrated in various standard field solvers such as FDTD [28]- [31] and Integral Equation (IE) methods [6], [7], for instance. The IE-GSTC implementation of (17) is presented in Part 2 of this work [23].

B. Lorentz Oscillator Model
To understand the angle dependent unit cell resonances captured by the poles of (13), and their physical origins, consider an arbitrary sub-wavelength unit cell structure excited by an oblique uniform plane wave (i.e. incidence angle θ or k y = k 0 sin θ), which can be described using a standard Lorentz oscillator model: which describes the temporal electric (and magnetic) surface polarization function in response to the average electric (and magnetic) fields around the surface, and where ω 2 p is the plasma frequency, γ is a damping coefficient and ω 0 is the resonant frequency at that specific angle of incidence. In the temporal frequency domain, we can express this using a temporal Fourier transform, F t {·} as: To account for the angular dependence of the unit cell, let us heuristically assume that the damping coefficient and the resonant frequency are polynomial functions of the incoming plane-wave angles such that Consider terms up to the second order for γ, ω 2 p and ω 2 0 . Then, substituting (21) in (20), allows us to relate the average fields with the polarization in the following form, where, .
For the case of spatially symmetric unit cells that is assumed throughout this work, the terms ξ 1 in (22) will be zero. Finally, taking an inverse spatial Fourier transform F −1 y {·} of (22), we get (along with the analogous equation for the magnetic surface polarization) ξ zz ee,2 which appears as a spatial counterpart of (19) and represents the spatial boundary condition across the zero thickness sheet. While the general form of (13) is applicable in general, the Lorentzian form of (22) represents an important special case motivated by physical considerations.

IV. APPLICATION TO PRACTICAL METASURFACE STRUCTURES
To illustrate the proposed method, we will consider two example metasurfaces composed of a 2D array of a) a short electric dipole, and b) a cylindrical dielectric puck, respectively, lying in the y-z plane, with x-y as the plane-ofincidence and TE mode excitation. Both the unit cell structures exhibit symmetry about the origin so that their transmittance/reflectance is an even function of θ and have identical and reciprocal responses for left and right excitations. It can be shown that under these conditions, only χ zz ee , χ yy mm and χ xx mm are the possible non-zero susceptibility tensor components. Furthermore, it can be shown that for both these structures, χ xx mm is also zero, so that these structures are completely described in terms of tangential surface susceptibilities and are thus birefringent [13].

A. Short Metal Dipole
Let us consider the first example of a short electric dipole unit cell formed using a conducting wire, as shown in Fig. 2(a), excited with a TE mode (E z , H y , H x ). It is simulated in Ansys FEM-HFSS using Floquet boundary conditions, where its transmittance and reflectance are computed for a sweeping angle of plane-wave incidence, θ, using (3). At 60 GHz and for different wire lengths, , Fig. 2(b) shows that the spatial resonance is located at normal incidence (θ = 0 • ) when = 2.5 mm, and it moves to higher angles as the wire length decreases. For a fixed length = 2.5 mm, Fig. 2(c) shows a typical surface susceptibility distribution as a function of both temporal frequency, ω and spatial frequency k y . The magnetic susceptibility is found to be negligible and thus not shown. It is clear that the surface susceptibilities are strongly angular dependent, and thus the resulting metasurface is expected to be spatially dispersive. Moreover, we observe a strong resonance migrating towards lower temporal frequencies for increasingly oblique angles.
Next, to capture the angle-dependent surface susceptibility, and in particular, a single angle-dependent resonance of the structure, the Lorentz oscillator model of (22) is used to numerically curve-fit this response. Fig. 2(c-d) shows the reconstructed electric and magnetic surface susceptibility profile across both temporal frequencies, ω and spatial frequencies k y = k 0 sin θ. There is a remarkable agreement between the full-wave simulated susceptibility and the reconstructed one, despite noisy data from HFSS, possibly due to poor convergence specially at higher angles. It thus confirms that the Lorentz model with only six non-zero parameters (see Tab. I) fully describes such a complex response of this unit cell. Consequently, the extended GSTCs for this structure are simply given by (23) with χ yy mm = 0.

B. Huygens' Metasurface
Next, consider an all-dielectric resonator structure, which consists of a cylindrical dielectric puck made of high permittivity material embedded inside a host medium of lower permittivity (assumed air for simplicity here), as shown in Fig. 3(a). All-dielectric structures are common, especially at optical frequencies as Huygens' structures (co-located orthogonal electric and magnetic dipoles), due to their low-loss characteristics and their zero backscattering property as further shown in Fig. 3(a) [17]. They have also been proposed at millimeter-wave frequencies in both all-dielectric [32]- [34] and standard printed circuit board implementations [18]. One common feature among all these structures is their relatively large unit cell sizes which can approach free-space wavelength, making them weakly sub-wavelength.
The typical angle-dependent electric and magnetic surface susceptibilities of a dielectric puck are shown in Fig. 3(b), simulated in FEM-HFSS, and susceptibilities extracted using (3). Compared to the simpler unit cells of Fig. 2, the dielectric unit cell features a more complicated angular dependence of the susceptibilities. Specifically, it shows multiple angular resonances across the angular spectrum drifting across the temporal frequencies, which suggests that a single Lorentzian oscillator model of (22) is insufficient to model this structure. While one can use a brute force approach to fit this response using rational polynomials of the form of (13), we observe that at each temporal frequency, surface susceptibilities appear as a summation of several angular resonances with Lorentz characteristics. Consequently, we next develop a multi-Lorentz surface description of this structure.
Let us assume that the electric surface polarization components (and analogous development for the magnetic ones) are described using N L Lorentz resonators and that, without loss of generality, each resonator's properties are up to secondorder dependent on k y . Then, (22) assumes the form, where the constant term, χ zz ee0 is added for generality. Taking the inverse spatial Fourier transform, F −1 y {·} on each side, we get Furthermore, from the GSTCs ∆H y,i = jωP z,i and using superposition and the independence of the polarizations, we retrieve the extended form of the GSTCs for this case as where the 0 th order term is, It should be noted that a sum of Lorentzians (24) may equivalently be seen as a ratio of two polynomials of (13) when all combined, which in general may lead to spatial derivatives higher than 2. However, their decomposition in terms of Lorentzians not only provides a physically motivated description, it allows us to write the spatial boundary conditions involving spatial derivatives up to second-order only, which is more suitable for numerical implementation.
To demonstrate the multi-Lorentz description of the alldielectric cell, Fig. 3(c) shows the electric and magnetic surface susceptibilities at couple of arbitrarily chosen frequencies and as a function of angle. An excellent fitting is observed across the entire angular spectrum, confirming the suitability of (24) to model this unit cell. Moreover, a complete description of the surface polarizabilities requires several Lorentz resonators to capture the whole angular response at discrete frequencies. Data in Tab II presents the values of the parameters for each frequency. It further reinforces the advantage of the presented methodology in describing spatial dispersion using a compact set of parameters.

V. SPATIAL DISPERSION VS NORMAL SURFACE POLARIZATIONS
So far, we have considered spatially dispersive structures which exhibit angle-dependent tangential surface susceptibilities only and showed how they lead to extended GSTCs, which can describe the angular scattering from the metasurfaces. In a  variety of other structures which are spatially non-dispersive, the angular scattering is described using both the normal and the tangential surface susceptibility components, dictated by their respective physical mechanisms. For example, consider a resonant loop structure consisting of a Metal-Insulator-Metal (MIM) capacitor printed on a thin dielectric slab, as shown in Fig. 4(a). At oblique incidence, the time-varying magnetic flux through the loop induces an electric current around the conducting loop via Faraday's law, which leads to a strong normal magnetic polarization along the x−axis. Consequently, it has been shown that the angular scattering of such a unit cell structure can be accurately modeled using one tangential surface susceptibility, χ zz ee and one normal susceptibility component χ xx mm , with no spatial dispersion [35], [36].
As seen, spatial dispersion and normal surface susceptibility components are related to determining the angular scattering from the surface. Therefore, one may wonder if they are related to each other or represent two independent properties of a given metasurface? For instance, is it possible to model the unit cell structure of Fig. 4(a) using purely tangential surface susceptibilities which are angle-dependent and thus spatially dispersive? To answer this question, we recall that the transmittance and reflectance of a metasurface described using one tangential surface susceptibility,χ zz ee and one normal susceptibility componentχ xx mm , are given by where θ is the angle of incidence of an incoming plane-wave (TE mode), and k 0 is the free-space wavenumber. On the other hand, the transmittance and reflectance of a metasurface described using a single tangential surface susceptibility, χ zz ee , and no normal component, are given by Comparing (25) and (26), it is clear that for R ⊥ = R || and T ⊥ = T || for every angle of incidence θ, we must have, where the effect of normal componentχ xx mm has been absorbed in the new tangential component χ zz ee . To confirm this equivalence, Fig. 4(b) shows the extracted χ zz ee directly from FEM-HFSS using (3), which shows an angle-independent resonance frequency, but whose plasma frequency varies with the angle of incidence (as seen from the line-width broadening of the resonances). Furthermore, Fig. 4(c) presents the comparison between the extracted susceptibilities at 10 GHz and the analytical relation of a fictitious χ zz ee obtained from (27). A very good fit is observed, confirming this equivalence for a uniform plane-wave incidence.
While this demonstration may suggest that a spatially nondispersive surface with normal surface susceptibilities may be represented using a spatially dispersive metasurface with tangential surface susceptibilities only, this equivalence is not generally correct. It is, in fact, applicable for a uniform metasurface only and can be explained by considering the GSTCs for the case of a general non-uniform metasurface. Specifically, for a metasurface described using tangentialχ zz ee and normalχ xx mm , we get, On the other hand, the GSTC for a metasurface with a single tangential χ zz ee described using (27), reads: It is clear that (28) is only equal to (29), if ∂χ xx mm /∂y = 0, i.e. a uniform metasurface. Therefore, we can conclude that for a general nonuniform metasurface, spatial dispersion and the normal surface susceptibilities (if they physically exist) represent two different properties which must be taken into account simultaneously to accurately describe its complete angular scattering.

VI. CONCLUSIONS
A simple method to describe spatially dispersive metasurfaces has been proposed where angle-dependent surface susceptibilities are explicitly used to formulate the zero thickness sheet model of practical metasurface structures. It has been shown that if the surface susceptibilities can be expressed as a ratio of two polynomials of tangential spatial frequencies k || , that captures their zero and pole behaviors, they can be conveniently expressed as spatial derivatives of the difference and average fields around the metasurface in the space domain. They thus represent extended GSTCs accounting for spatial dispersion. Using two simple examples of a short electric dipole and an all-dielectric cylindrical puck unit cells, which exhibit purely tangential surface susceptibilities and symmetric transmission and reflection responses, the proposed concept is numerically confirmed in 2D. A single Lorentzian has been found to describe the spatio-temporal frequency behavior for the short dipole, while a multi-Lorentzian description was necessary to capture the multiple angular resonances of the dielectric puck. In both cases, the appropriate spatial boundary conditions have been provided. We further emphasize that the proposed approach models very complex spatio-temporal responses of the two unit cells considered here, using very few parameters, thereby making them ideal compact simulation models of such structures to be easily integrated in standard electromagnetic field solvers.
As has been shown, the generalized expressions of the surface susceptibilities as a ratio of two polynomials in k y or as a sum of Lorentzian oscillators reveal themselves as spatial derivatives of the fields around the metasurface, resulting in more general boundary conditions than standard GSTCs, that can be incorporated in a variety of numerical methods to compute the scattered fields from spatially dispersive metasurfaces. For example, Part 2 of this work will demonstrate this integration of extended GSTCs in an Integral Equation (IE) based field solver, where the scattered fields may be computed in the temporal frequency domain for an arbitrary incidence wave and a given metasurface configuration.
Moreover, in this work, the analysis has been limited to structures exhibiting purely tangential surface susceptibilities. However, it has been shown that to capture the complete angular scattering properties of a general metasurface, both spatial dispersion, and normal surface susceptibilities must be taken into account if the structure physically supports them. While the proposed method can be extended to include the normal components, developing the method to integrate them will be an important step. Last but not the least, general modeling of spatially dispersive non-uniform metasurfaces will be a natural extension where either the complex coefficients of (13) or the parameters of the Lorentz oscillator become function of space. This work marks the first essential step to achieve this goal.