Part 2: Spatially Dispersive Metasurfaces -- IE-GSTC-SD Field Solver with Extended GSTCs

An Integral Equation (IE) based field solver to compute the scattered fields from spatially dispersive metasurfaces is proposed and numerically confirmed using various examples involving physical unit cells. The work is a continuation of [1], which proposed the basic methodology of representing spatially dispersive metasurface structure in the spatial frequency domain, $\boldsymbol{k}$. By representing the angular dependence of the surface susceptibilities in $\boldsymbol{k}$ as a ratio of two polynomials, the standard Generalized Sheet Transition Conditions (GSTCs) have been extended to include the spatial derivatives of both the difference and average fields around the metasurface. These extended boundary conditions are successfully integrated here into a standard IE-GSTC solver, which leads to the new IE-GSTC-SD simulation framework presented here. The proposed IE-GSTC-SD platform is applied to various uniform metasurfaces, including a practical short conducting wire unit cell, as a representative practical example, for various cases of finite-sized flat and curvilinear surfaces. In all cases, computed field distributions are successfully validated, either against the semi-analytical Fourier decomposition method or the brute-force full-wave simulation of volumetric metasurfaces in the commercial Ansys FEM-HFSS simulator.


I. INTRODUCTION
Electromagnetic Metasurfaces are constructed using subwavelength resonators of diverse geometrical shapes and material characteristics, which give rise to their tailored macroscopic responses. These microscopic resonators act as engineered electric and magnetic scatterers operating on the incoming incident fields to shape the electromagnetic fields in either space or time or both [2]. This has led to a powerful wave engineering paradigm, leading to a variety of applications across the electromagnetic spectrum, ranging from cloaking and illusions to holograms, and real-time reconfiguration of wireless environments, to name a few [3]- [10].
Since the metasurfaces are constructed from sub-wavelength resonating particles but are typically electrically large, they are naturally multi-scale in nature. As a result, the determination of electromagnetic fields scattered off them for a given incident field is typically a computationally expensive task if performed using standard brute force commercial fullwave simulators. Consequently, surface susceptibilities,χ(ω) have recently been used as compact full-wave simulation Tom J. Smy, João G. Nizer Rahmeier, Jordan Dugan, and Shulabh Gupta are with Carleton University, Ottawa, Canada (e-mail: tjs@doe.carleton.ca). models of wide-variety of practical metasurfaces with good success [11], [12]. They represent zero thickness sheet models involving electric and magnetic surface polarization densities, which when coupled with the Generalized Sheet Transition Conditions (GSTCs) [13], [14], can describe the average macroscopic fields around the metasurface and their angular scattering properties for arbitrary incident fields [15]- [18].
The surface susceptibilities relate the average fields with the induced surface polarizations, and so far in the literature, they have been typically restricted to modeling local interactions with the surface only. Such a point-by-point interaction consequently results in angular independent surface susceptibilities, and the metasurfaces can be referred to as spatially nondispersive. This has been found to be sufficiently adequate for modeling deeply sub-wavelength resonant structures. However, typical practical metasurfaces are not always deeply subwavelength, and their unit cell periodicities may easily reach close to a wavelength in size, as, for instance, in all-dielectric Huygen's structures [19]- [22]. In such cases, the field interaction with the metasurfaces is typically non-local, and the metasurfaces become spatially dispersive. Such structures, in general, cannot be modeled using constant dipolar surface susceptibility models (even if the normal surface polarizations are included) to describe their angular scattering behavior.
Very little work has been done to model spatially dispersive metasurfaces and has been typically limited to weak spatial dispersion (SD) with limited success [23]- [26]. To address this problem, in Part-1 of this work [1], we have proposed a simple method to model a general spatially dispersive metasurface, whose angle dependence of the surface susceptibilities have been expressed as a ratio of two polynomials of the spatial frequency, k || (i.e., the transverse wave vector). This, when transformed into the spatial domain via inverse spatial Fourier transform, leads to an extended form of GSTCs, which involves spatial derivatives of both the difference and the average fields around the metasurface. To the best of our knowledge, this has been the most general treatment of spatially dispersive metasurfaces in the literature.
As mentioned above, surface susceptibilities are compact zero thickness models of otherwise finite thickness volumetric metasurfaces, which are ideal alternatives to rigorously and efficiently compute scattered macroscopic fields. While several works have been reported in the literature integrating standard GSTCs into a variety of numerical platforms such as Finite Difference (FD) [27]- [29] and Integral Equation (IE) based methods [7], [10], [11], [30], for instance, they naturally have been limited to modeling spatially non-dispersive metasurfaces with angle independent surface susceptibilities. This work (Part-2) continues the general treatment of spatially nondispersive metasurfaces from Part-1 and integrates the general spatial boundary conditions via the extended GSTCs into the IE simulation framework. We develop the integrated matrix formulation of extended GSTCs and the IE-based field propagation, rigorously and efficiently describing field scattering in response to arbitrarily specified incident fields.
This work is structured as follows. Sec. II presents the extended GSTCs and develops the fields equations that describe a general spatially dispersive metasurface and present an important specialization based on a physically motivated Lorentz oscillator. Sec. III presents the integrated IE-GSTC formulation equipped with spatial dispersion, i.e., IE-GSTC-SD. This section transforms the desired fields equations into a matrix form and integrates with the IE field propagators, resulting in a system matrix formulation that must be selfconsistently solved to compute the unknown surface currents and the subsequent scattered fields. Sec. IV presents the verification of IE-GSTC-SD by comparing it with a semianalytical Fourier decomposition method applied to uniform structures excited with 2D Gaussian beams. Next, the method is used to model an example practical structure composed of an electric dipole formed using a finite-length conducting wire, which exhibits an elementary Lorentzian form of spatial dispersion (Sec. V). Fields scattered from finite-sized flat and curvilinear structures are presented and compared with bruteforce 3D full-wave simulations performed in Ansys FEM-HFSS to successfully retrieve intricate interference field patterns that are otherwise missed using spatially non-dispersive models. Finally, conclusions are presented in Sec. VI.

II. SPATIALLY DISPERSIVE METASURFACES A. Generalized Sheet Transition Conditions (GSTCs)
The electromagnetic field interaction with an equivalent zero thickness metasurface model is governed by the GSTCs as given by whereP andM are the temporal frequency domain electric and magnetic surface polarizations. These surface polarizations, for a spatially non-dispersive metasurfaces, following local field interactions are given bỹ where the susceptibility tensors are 3 × 3 in size and defined by: where ab ∈ [ee, mm, me, em]. For later simplicity, but without loss of generality, let us ignore the bi-anisotropic terms and the normal polarization, so that we have, Each of the susceptibility terms present in these equations creates a contribution to the polarizations in (2); for example, χ zz ee creates a component (e.g., for a metasurface lying along y−axis and located at x = 0 excited with a TE mode, with E z , H y and H x components), where the superscript on P z z indicates this is the portion of P z generated by E z,av . The total P z would thus be, This equation represents that the electric (and analogously magnetic) surface polarization at any position of the metasurface is given by the spatial product of the susceptibility and the average fields at that position only, i.e., a local response. Alternatively, it can be seen as assuming a set of angle independent susceptibilities that only captures the field scattering behavior of spatially non-dispersive metasurface unit cells.

B. Spatial Dispersion and the Extended GSTCs
For a spatially dispersive metasurface, the induced electric (and magnetic) surface polarizations on the surface not only depend on the local average fields but also the fields across the metasurface [1]. Assuming a vertical surface at x = 0 we can write,P z z (y) = 0 χ zz ee (y) * Ẽ z,av (y) where * represents a spatial convolution. This relationship can be equivalently expressed in the spatial frequency domain, k y (each k y represents an obliquely propagating plane-wave), using the Fourier transform F y {·}, so that 1 which is a simple product of the average fields and the surface susceptibilities, as opposed to convolution in space, making it a suitable choice for compact unit cell description and later numerical computation. In order to describe a general unit cell, we can represent χ zz ee (k y ) as a ratio of two polynomials to capture the angular dependence of the metasurface. Such a form incorporates possible zeros and poles of the surface susceptibilities, which effectively captures the metasurface response for the sweeping angle of incident plane-waves [1]. Such a form reads, where a m and b n are known complex coefficients, which can be extracted from unit cell simulations of a given metasurface structure.
To express this relationship in the spatial domain and represent a spatial boundary condition across the metasurface, we utilize the GSTCs of (1). For a vertical surface with surface normaln = {1 0 0}, we note that ∆H y = jωP z , and thus we can associate a portion of ∆H y with each portion ofP z and transform back to the spatial domain. This would give for thẽ P z z (k y ) portion, n b n j n ∂ n ∆H z y ∂y n = jω m a m j n ∂ mẼ z,av ∂y m (9) resulting in general form of (4) that incorporates spatial derivatives of both ∆H z y andẼ z,av . The total field difference is then given by, with an equation of the form of (9) being defined for all three terms. Likewise, we would define similar relationships for the other components (∆Ẽ z ), i.e.
The two boundary equations (9) and (11), are now referred to as the extended GSTCs, which can now be applied to the general fields scattering problem as illustrated in Fig. 1.

C. Lorentz Oscillators
Although the derivatives of the difference and average fields present in (9) and (11) can be of use and are amenable to the methods presented in this paper, we shall consider the analysis of an important subset which can be physically interpreted and understood, i.e., a Lorentzian oscillator. Let us consider that a collection of Lorentz resonators can describe the resonant response of a metasurface unit cell in the temporal-and spatialfrequency domain as, again usingP z z as an example, where {ω 0 , ω p , γ} are the resonant frequency, plasma frequency, and the damping coefficient, respectively, which all depend on the geometrical and electrical characteristics of the unit cell. Moreover, they also depend on the angle of incidence of the incoming plane-waves [1]. We can do a Taylor expansion, for instance, around the normal incidence (i.e. k 0 sin θ = k y = 0), so that Using (13) in (12), as shown in [1], we can relate the average fields with the polarization in the following form in the spatial frequency domain: which upon an inverse spatial Fourier transform leads to [1], Equation (15) is a simplification of (9) with only 2 nd order terms on the left-hand side. A more complex form considers 2 nd order derivatives of the average fields as well, and is considered in a general form in Sec. II-D. It is worth mentioning that for a symmetric spatial dispersion condition, the term ξ 1 in (14) is zero. That is the case for the practical metasurface cells we will analyze later in this paper for simplicity and illustration. This form thus represents a very simple physically motivated case of a spatial dispersion, which is exhibited by a simple array of short conducting wires, as shown in [1] and to be used later in the example section of this work. For general and more complex unit cell architectures (such as Huygens' structures, for example), we can assume that the surface polarization components are described using N L Lorentz spatial resonances, so that Eq. (14) assumes the form, Taking the inverse spatial Fourier transform, As before, we use ∆H z y,i = jωP z z,i , then using superposition and the independence of the various polarization components, we get, defining, where the 0 th order term is Which describes the surface relationship between a component of ∆H y andẼ z,av . 2 To complete this formulation, we can write each Lorentzian equation associated with the tensorial form of the susceptibilities and introduce a generalized differential operator L, where, with α ∈ {x, y, z}, β ∈ {x, y, z} and γ is defined byn × {x, y, z} with respect to value of α.

D. Generalized Surface Equations
The previous section shows the surface formulation for the case of (16) with a second-order denominator and constant numerator (no zeros). A more general (but still 2 nd order equation) can be obtained by using a 2 nd order polynomial in k y for the numerator of (16), This form incorporates spatial derivatives of the average fields, when transformed to the space domain, and can be related to the normal components of the surface polarizations in some special cases. This extension will later allow for a very simple test to validate the IE-GSTC-SD formulation using a well characterized unit cell with constant normal surface susceptibility. By following the same procedure as used for the field differences on the left-hand side of (18) to handle the spatial derivatives of the average field, we define two more operators associated with each component, we can generalize (23) to This generalized form can capture more complex angular dependence arising from the structure of the unit cell [1]. We should note that we have assumed a 2 nd order for both operators, but higher orders could be used. However, odd orders will produce reflection asymmetry and may be omitted for physical reasons.

III. IE-GSTC-SD FORMULATION
Our next task is to integrate the zero-thickness boundary condition of (27) into the bulk Maxwell's equations, to develop a general-purpose field scattering solver accounting for spatial dispersion in an IE approach with one or more surfaces. Firstly, we formulate the field equations at the surface for a discretized surface by transforming equations such as (22) into matrix equations suitable for incorporation into a system-level solution. Secondly, we form the propagation equations for the scattered fields generated by an implicit incident field, which requires propagation matrices for a self-consistent scattering solution. Finally, the system matrix and source vector must be formed to provide an equation to be solved. We will now take each of these tasks in turn.

A. Surface Discretization
The IE approach, in general, requires discretized surfaces and follows the Boundary Element Method (BEM) techniques [31]- [33]. These surfaces connect regions of homogeneous material properties, coupling the regions together through transmissive and reflective properties. Although a defined surface must surround a region, a portion of the surface may be at infinity (where the fields are assumed to be zero) or implicitly present and not actually modeled. The characteristics of the implicit surface will depend on the assumptions about the implicit excitation.
The problems addressed in this paper will be 2D scattering cases and the surfaces are curvilinear line elements, as illustrated in Fig. 1. For surfaces which are to be physically modeled, we impose a discretization using a uniform segmentation of each surface. Each discrete element of the surface is characterized by a center position r p,i = [x p , y p , 0] i , a length i and a normaln i which are collected into vectors such as r p = r p,1 . . . r p,m . The field quantities are stored in vectors of the same form (with m surface elements): The surface field vector S F holds the fields present on both sides of a surface (+ or −) and will be useful in the formulation of the IE problem [11]. The fields are stored sequentially along the surface and can be thought of as triplets of fields (E i and H i ). Therefore, facilitating the operators' use to implement the spatial derivatives in the general formulation given above.

B. Surface Equations
The GSTCs shown in (1) when formulated using matrix operators and omitting the bi-anisotropic terms for simplicity are where both N T and R T extract the tangential components of a field, but R T incorporates the rotation induced by the cross-product. The addition of the GSTCs with constant angular independent susceptibilities (i.e., spatially non-dispersive structures) to the IE formulation has been reported extensively elsewhere, and we will refer to it as the standard IE-GSTC formulation [8], [17], [33].
We will now develop a spatially dispersive implementation of the GSTCs within the IE framework (IE-GSTC-SD). If we assume a vertical surface at x = 0 running along the y axis, then we have, which extracts the tangential components of the field at the surface, for example E T = E y E z T and We now simplify the following derivations (without losing any generality) by only retaining the primary tangential terms and have, To generalize the relationship between ∆H andẼ av two generalized Lorentzian Operators are introduced [X ee /L ee , X mm /L mm see (27)] and we use the relationshipP = jω∆H for each component, where L ab ee/mm,i and X ab ee/mm,i are defined by (24) and (26) The primary complication in the implementation of these equations is the presence of the spatial derivatives in the L and X operators defined by (24) and (26). However, a similar issue was solved in [17] to allow for the incorporation of the terms involving the gradient of the normal component of the polarizations in (1), and we will follow a similar approach here.
We define central difference operators for the n th segment triplet: where Ψ would be a field difference associated with a particular resonance, such as ∆Ẽ z,i , and obtain, which allows us to write, and where U = [0, 1, 0]. By applying this equation to each segment on the surface we can arrive at the surface operators, which are diagonal matrices created from the triplet operators and would operate respectively, on all of the field differences defined by, and the surface scattered fields S F and applied incident fields S i to obtain, which is the surface level equivalent of (33), a generalization of the GSTCs, and completely defines the tangential field relationships at the surface.

C. IE-GSTC-SD Field Equations
As with the surface equations, the integral expressions derived from Maxwell's equations that capture the propagation need to be put into a discretized form. The EM fields radiated into free-space from electric and magnetic current sources, {J, K}, can be generally expressed using an IE formulation as [17], [31], [32]: with r being the point of interest, r the position of the source current; and E s and H s the radiated (scattered) fields from the surface. 3 The field operators are given by: withC ∈ {J,K}. G(r, r ) represents the Green's function which, for a 2D case, is given by the Hankel function of the 2 nd kind and represents outwardly propagating radial waves.
Using this discretization, (39) is transformed into a set of algebraic equations relating surface currents C = J K T to the scattered fields at r S , If r p = r S then this equation determines the selfpropagation from every element to every other element and we can define a surface propagation operator, and obtain,

D. System Formulation
The final task is to assemble the surface (38) and propagation equations (42) to solve for the unknowns in the scattering problem. The unknowns in the scattering problem can be identified as the element currents C, the surface fields S F , and the field differences associated with the Lorentz resonators S ∆ . To create a complete system matrix, we have to introduce equations to force the unknown currents to be tangential to the surface, to sum the field differences, and link them to the propagation equations.
To force the surface currents to be tangential to the surface, we introduce an operator, this would impose the conditionn·C = 0 for all elements. For our simple case of a vertical surface, this implies thatJ x and K y are zero for all elements. To link the propagation equations to surface equations, we define two new operators such that, where L S is a matrix that sums the field differences associated with each resonance to form the total field differences ∆Ẽ z and ∆H y at every element -D is the matrix that finds the difference of the fields on the two side of the surface. This equation thus links the propagating fields S F to the differences associated with the resonances present S ∆ . These equations are finally assembled into a system matrix equation, resulting in which is square and can be solved directly to determine the unknowns. Once C is known, the propagation matrix defined in (40) can be used to compute the scattered fields in the entire or desired simulation region.
IV. NUMERICAL VERIFICATION To numerically verify the IE-GSTC-SD formulation, we will use two methods using 2D Gaussian beam illumination of uniform metasurfaces. The first approach will compare a semi-analytical technique using Fourier Decomposition (FD) of the incident waves into plane wave components and then determine reflected and transmitted fields. The second method will take a specific form of SD (angular dependence) that can model an otherwise spatially non-dispersive metasurface with two constant susceptibilities -one of which operates on a normal component of the magnetic field. This second approach allows for a direct comparison between the standard IE-GSTC [11] and the proposed IE-GSTC-SD formulations for the same surface.

A. Fourier Decomposition Method
Let us consider a case where an incident field is specified, and we wish to determine the scattered fields from a uniform spatially dispersive metasurface. Since it is a linear problem, we can use the principle of superposition by expressing the incident field as a sum of uniform plane waves [34]. The two GSTCs for the case of a vertical surface with simple tangential susceptibilities and an assumed TE polarization (w.r.t the normal x) given by {E z , H x , H y } can be expressed in the spatial frequency domain as, Introducing T =Ẽ t /Ẽ 0 and R =Ẽ r /Ẽ 0 , we can show that the reflection of single plane wave at an incident angle is given by, [1] R(θ) = 2jk 0 {cos θ 2 χ yy mm − χ zz ee } {jk 0 χ zz ee + 2 cos θ}{jk 0 cos θχ yy mm + 2} where χ yy mm and χ zz ee are functions of the incidence angles. The above equation gives the transmitted and reflected plane waves for a particular spatial component of a general incident field. These can then be integrated over all the k y components (propagating terms) of the incident field to construct the complete reflected and transmitted fields [34]. Specifically, if we represent the incident field at the surface asẼ i z (y), then we can transform this to the spatial Fourier domain using, where we have decomposed the incident field into plane wave components. We can then form the reflected and transmitted waves, with sin θ = k y /k 0 . Using an inverse spatial Fourier transform, we can obtain the scattered fields at the surface, This methodology can be easily implemented using discrete Fourier transforms (DFTs) and can now be used to validate the IE-GSTC-SD implementation.

B. 2D Gaussian Beam Propagation
Consider a uniform metasurface described using tangential susceptibilities only and with unit symmetry about both x− and y− axis. Assume that a single Lorentz resonator can describe it for both χ zz ee and χ xx mm , as The Lorentzian parameters are next synthesized so that the surface reflection is symmetrical with respect to the incident angle (ξ ee/mm,1 = χ ee/mm,1 = 0) and that at normal incidence we have T = 0.9j and R = j 1 − |T | 2 (arbitrarily chosen).
To determine the Lorentzian parameters of (50), we can invert (46) for θ = 0 • and find the two parameters, χ zz ee (0) and χ yy mm (0). We then chose to set χ zz ee,0 = |χ zz ee (0)| and χ zz ee0 = χ zz ee (0)−χ zz ee,0 . The same procedure was used for χ yy mm . Finally, we set ξ zz ee,2 and ξ yy mm,2 , which produces resonances at specified incidence angles of 20 • and 30 • , respectively (again arbitrarily chosen). These synthesized susceptibilities of the form (50), are shown as a function of incidence angle (or alternatively vs spatial frequency k y = k 0 sin θ) in Fig. 2. As specified, one angular resonance is placed in each susceptibility component with symmetric response about θ = 0 • . We now wish to determine the scattered fields when it is illuminated with a 2D Gaussian beam, i.e. spatially broadband in k y .
A 2D Gaussian beam (GB) is a solution of the Paraxial Helmholtz equation, with ∂/∂z = 0, so that there is no spread along the z−direction and the propagation is confined within the x − y plane only. The various field components of a 2D Gaussian beam propagating normally to the surface along x−axis can be easily shown to be: where w is the width of the beam at the waist x = 0, k 0 = ω/c is the free-space wave-number. The magnitude of the Gaussian beam with waists of w = 2λ and w = 0.75λ, respectively are also shown in Fig. 3(a). Moreover, their angular spectra at the metasurface location are also shown as a function of incidence angle θ in Fig. 2. 4 The angular content of the two beams can be compared to the position of the resonances in the susceptibilities, and it can be seen that the wider spatial beam content lies within the two resonances, and little spatial dispersion is expected. On the other hand, the narrow spatial 4 The incident field given by (51b) was transformed to the spatial domain and then sin θ = ky/k 0 was used to translate this into the angular domain. beam has a substantial amount of energy at angles at or above the resonances, and we expect to see a substantial distortion of this beam as it interacts with the metasurface. The first case presented in Fig. 3 (top four plots) is for the moderately wide Gaussian beam (w = 2λ) with, as we have seen, little angular content. As would be expected, the reflected fields for all three cases are quite similar, as only in the presence of significant values of k y (high angular content) are the methods expected to differ. For the FD and the IE-GSTC-SD, we do see some marginal angular dispersion of the reflected and transmitted wave and an excellent match between the two. In Fig. 3(b) (top) with no SD with ξ zz ee,2 = ξ yy mm,2 = 0, we essentially see a simple reflection/transmission of the Gaussian beam recreating the shape of original incident fields. Some distortion in the FD and IE-GSTC-SD fields is present, indicating a small amount of spatial dispersion. Nevertheless, the match between the two methods is excellent.
The second case (shown in the bottom four plots) presents  5. Comparison of the scattered field generated from a loop resonator unit cell-based metasurface when excited with a 2D Gaussian beam, using two independent methods of spatially non-dispersive tangential and normal surface susceptibilities (Case 1) [11], [35], and an equivalent spatially dispersive model with angle-dependent tangential surface susceptibilities only (Case 2). a) The unit cell configuration. b) Scattered fields obtained using standard IE-GSTC for Case 1 and proposed IE-GSTC-SD for Case 2, for normal incidence and oblique incidence. The operating frequency is 60 GHz, and various susceptibility values are tabulated in Tab. I. a much clearer spatial dispersion effect. The Gaussian beam now has a small waist of 0.75λ, and angular content is very obvious in the incident field shown in Fig. 3 bottom. A very strong distortion of the Gaussian beam is observed in both reflection and transmission, indicating a major re-arrangement of various spatial frequencies resulting from spatial dispersion. Again, an excellent match is seen with the FD method, indicating the correct implementation of various spatial derivatives in the extended GSTCs. In contrast, if we were to model this surface by ignoring the spatial dispersion, with constant χ zz ee and χ mm yy and equal to their nominal values at normal incidence, we observe a simple reflection and transmission of the incident Gaussian beam, which naturally is not correct. To further confirm the match between FD and IE-GSTC-SD, the transmitted and reflected fields for the case of a broadband 0.75λ beam are compared along the observations lines (dotted lines on the plot) 2λ away from the surface and are shown in Fig. 4. They are practically identical. This result thus indicates the agreement of two independent methods, IE-GSTC-SD and FD, for capturing the spatial dispersion due to the angular resonance and is an initial validation of the proposed IE-GSTC-SD methodology.

C. Spatially Non-dispersive Unit Cell with Normal Susceptibilities
For a second verification of the proposed IE-GSTC-SD, method, let us consider a practical unit cell as shown in Fig. 5. It is formed from a resonant loop structure consisting of a Metal-Insulator-Metal (MIM) capacitor printed on a thin dielectric slab, as shown in Fig. 5(a). It has been shown that such a unit cell structure can be modeled using one constant tangential surface susceptibilitiesχ zz ee and a normal susceptibility componentχ xx mm , with no spatial dispersion [11], [35]. Including the gradient of the normal fields in (1), leads to a single field equation governing the field scattering and is given by which incorporates the effect of χ xx mm . As shown in Part-1 [1], we can model this uniform structure using purely tangential surface susceptibility χ zz ee , dependent on the angle of incidence θ, and thus spatially dispersive, i.e.
which is of the form (8), with a second order polynomial in the numerator and zeroth order denominator. Therefore, a standard constant dipolar normal surface susceptibility component accounts for a particular form of angular scattering from a uniform metasurface given by (53). However, it also presents an opportunity to verify the spatial dispersion modeling used in the IE-GSTC-SD, as it can be compared to the standard IE-GSTC formulation using constant surface susceptibilities. Such a demonstration is presented in Fig. 5 for the two cases of Gaussian beam illumination (normal and oblique incidence) on a uniform surface described by the two differently formulated surfaces. The surface parameters for the first case which uses constant susceptibilities, one tangentialχ zz ee,0 and one normalχ xx mm = 0, were reported in [17], [35] and are given in Tab. I. The second case uses only one spatially dispersive tangential component, χ zz ee , described by two parameters χ zz ee,0 and χ zz ee,2 , also tabulated in Tab. I. These two conditions are shown in Fig. 5(b) for the two cases of normal incidence and 45 • , respectively. In both cases, the predicted fields are essentially identical between  the two methods. This result verifies that the SD methodology introduced into the GSTC framework successfully models the angular dependence of a physical cell using SD tangential components of the susceptibility in which the angular dependence is due to a strongly dominant normal component of the susceptibility. Naturally, this demonstration also validates the IE-GSTC-SD implementation further.

A. Wire Dipole Based Unit Cell
To demonstrate the importance of spatial dispersion and capability of the proposed IE-GSTC-SD framework to model angular scattering from physical unit cells, we will consider an example of a simple unit cell based on a short conducting dipole, which exhibits spatial dispersion. The basic unit cell used to form a 2D surface is shown in Fig. 6(a) and consists of a short segment of wire of length and with a finite conductivity σ, placed inside free-space. This is then periodically arranged along y− and z−axis with periods of Λ z and Λ y to form a surface lying in the y − z plane. Due to symmetry considerations and assumed TE mode excitation, this cell can be shown to be modeled using two tangential susceptibilities, χ zz ee and χ yy mm only, which are angle-dependent. This unit cell is very simple yet very insightful, as it only exhibits a single Lorentzian resonance which dominantly depends on the length of the wire, thus acting as a perfect testbed for the IE-GSTC-SD framework.
In Part-1 of this work, this cell was characterized using HFSS to obtain its transmission and reflection characteristics as a function of frequency ω and angle of incidence, θ or spatial frequency k y [1]. Fig. 6(b) presents the extracted susceptibilities for a fixed frequency of f = 60 GHz and the fitted model with a single Lorentzian using (16). These susceptibilities produce an angular dependent reflectivity R(θ) and transmission T (θ) which can be obtained from the susceptibilities using (46), and are also shown in Fig. 6(b). The reflectivity is close to unity and exhibits a small amount of angular dependence varying from about 0.96 to 1. The effect of the two resonances present in the susceptibilities is more pronounced in the reflectivity, R dropping from a maximum of 0.2 at 55 • to a clear minimum at 30 • to nearly zero and also a drop off at high angles.

B. Finite-Sized Vertical Metasurface
To demonstrate the basic field transformation phenomena of this surface and to allow for a detailed comparison with a commercial full-wave simulator like Ansys FEM-HFSS, a vertical surface along the y−axis, consisting of 60 units cells was formed and illuminated with a line source (f = 60 GHz) located at r s = {x s = −30, y s = 0} mm. The incident field produced by the line source is given by, {0,1} are Hankel functions of the second kind, of orders 1 and 2, and E 0 is the peak field amplitude. The surface is simply placed in free space and the simulation set up, and the source position is shown in Fig. 6(c).
The basic behavior of the surface is presented in Fig. 7 where the total field E z is shown in the log scale to emphasize its low amplitude features. As a reference case, Fig. 7(a) first shows the total fields generated by a Perfect Electric Conductor (PEC) of the same size, which naturally generates zero transmission through the surface and exhibits finite diffraction at the two edges. Fig. 7(b) next shows the fields generates by the short conducting wire surface. Although the reflected fields are similar to those of a PEC, the transmitted fields are significantly transformed by spatial dispersion. The angular filtering is seen due to the drop in transmission at 30 • and the increased transmission at higher angles due to the spatial dispersion effect. A high level of transmission for incident flux at 45-65 • , interferes with the edge diffraction to produce quite a characteristic field pattern. In contrast, Fig. 7(c) shows the behavior of a uniform constant surface where the spatial dispersion is set to zero. We observe that the surface is highly reflective as expected, but there is no spatial filtering of the transmitted field in any significant way. Edge diffraction is present but minimal due to the incident flux at the edges arriving obliquely. These results clearly show that spatial dispersion strongly shapes the transmitted fields and, to some degree, the reflected ones. However, it remains to show that the fields depicted in Fig. 7(b) are indeed correct solutions of the surface phenomena. To assess this, we built the entire volumetric surface in Ansys HFSS and full-wave simulated it. This result is shown in Fig. 8(a), which has a very similar structure to the results obtained using IE-GSTC-SD in Fig. 7(b). The IE-GSTC-SD has captured the structure of the fields and interference of the edge diffraction very well. Hence, confirming that this equivalent zero thickness model is correctly capturing the angular scattering of the metasurface structure 5 .
However, there is an interesting feature in the fields obtained from the HFSS simulations related to the subtle interference pattern imposed on the fields. Although this pattern is at the limit of the HFSS simulation and meshing ability to resolve, it is present and thus investigated next. The simulation in Fig. 7(b) shows none of this fine internal field structure, and spatial dispersion alone does not seem a credible explanation. However, the wire dipole cell is not deeply sub-wavelength and is relatively large with Λ y /λ ≈ 0.6, as compared to the unit cell of Fig. 5(a), which is about λ/10 and thus safely with deep sub-wavelength periodicity. Therefore, the presence of a small dipole resonator in a large cell provides impetus to propose that the metasurface will act as a periodic surface with a finite number of current sources, as opposed to a sheet of 5 It can be noted that HFSS simulation took several hours to complete on a high-end server, where all possible symmetry boundary conditions were exploited, and yet it suffered from limited convergence. In contrast, the IE-GSTC-SD simulation took less than a minute on a desktop workstation once the angular-dependent surface susceptibilities are retrieved from unit cell simulations in HFSS, which were themselves computationally inexpensive. continuous currents. By default, the simulation in Fig. 7(b) used a discretization of 10 divisions per wavelength and models the surface current sources as quasi-continuous over the entire unit cell and the surface -as would be standard in a BEM approach -thus producing the smooth fields presented.
To investigate this phenomenon, we modified the BEM methodology and placed the total current present in the unit cell on a single segment at the position of the wire only. The unit cell had five surface segments across its width; therefore, on 4 of the segments, the surface currents were set to zero and the total cell current placed on the middle segment. The result of this quantization of the surface currents is shown Fig. 8(b). The effect is quite dramatic, unveiling the sought-out subtle interference pattern. Unlike the FEM method of HFSS, the interference pattern produced by the IE-GSTC-SD is very well defined due to the intrinsic nature of the method, which has a pure and simple description of the geometry. Figure 8(c) shows the generated fields when the spatial dispersion is switched off in BEM. The angular filtering features, as expected disappeared, while retaining the fine interference patterns only. Therefore, confirming that this interference pattern is solely due to current quantization and not spatial dispersion.
To further investigate the match between the FEM HFSS simulation and the IE-GSTC-SD results, the fields were plotted for two vertical lines ±0.085 m removed from the surface (arbitrarily chosen but in the macroscopic region). These fields can be seen in Fig. 9. Fields are presented for the three cases: FEM-HFSS, IE-GSTC-SD, and the quantized IE-GSTC-SD. The IE-GSTC-SD fields match the HFSS fields well, but of course, have none of the high-frequency variations. The use of the quantization method introduces the expected highfrequency modulation, and it is of a similar magnitude as those in the HFSS results. The details in the interference pattern do not match, but this is not really to be expected due to the small magnitude of the variation and the limits of HFSS simulation. Although this method of introducing the effect of discrete sources is somewhat ad hoc, it appears to confirm the source of the interference pattern and is quite successful.

C. Finite-Sized Curvilinear Metasurfaces
The final two examples will involve comparing Ansys FEM-HFSS and the proposed IE-GSTC-SD simulations for more geometrically complex structures created using the simple wire unit cell. We will consider two configurations: 1) an open hexagonal-based structure formed of three sides, and 2) an open semi-circular structure. For both simulations, the incident illumination will be a plane-wave at 60 GHz.
Simulations of the hexagonal-based structure are shown in Fig. 10. The structure consist of three sides of 21 cells and facet length of s = 0.0473 m. In this simulation, the incident plane wave travels left to right at 60 • measured from horizontal, striking the bottom facet of the open surface at normal incidence. Fig. 10(a) shows the fields obtained from HFSS simulations, which exhibit a complex field pattern due to the reflection from two exposed facets, transmission through the facets, and the intrinsic angular filtering of the wire structure. The facet in the shadow of the incident field also has a strong effect as it reflects back the field transmitted through the bottom and side facets. Some strong interference patterns are present due to the multiple reflections present in the interior of the hexagonal. We also see evidence of the individual wires acting as discrete sources in creating a more subtle interference pattern of the field in the interior of the open hexagon.
The second Fig. 10(b) presents a basic IE-GSTC-SD simulation of the structure with continuous current distribution accounting for spatial dispersion. It can be seen that this simulation produces a very close match with the HFSS results. It captures all of the basic features of transmission and reflection, including the interference patterns from the complicated reflection in the interior. Once again, to recreate the interior fine interference features, the currents were quantized, and the generated field patterns are shown in Fig. 10(c). As before, this quantization is speculative and ad hoc; however, as with the simulation of the vertical surface, we see the successful creation of a subtle interference pattern (particularly in the transmitted fields) as observed in the HFSS fields. Finally, to emphasize the importance of spatial dispersion, Fig. 10(d) presents the same simulation but with no spatial dispersion. The internal field structure due to complex interference is simply lost.
The final set of simulations (Fig. 11) are of a semicircular structure with a radius of r = 0.0342 m and consisting of 51 wire unit cells. The incident plane wave travels left to right at 30 • from horizontal striking the front of the open surface. Fig.  11(a) shows the HFSS simulation, which presents a complex transmitted field pattern. Besides, a complex spatial dispersion effect is expected since the plane wave is illuminating a curved structure with the angle of incidence varying along the surface. Due to the angular dependence of transmission, we see significant nulls in the field patterns internal to the curved surface, which also produce interference patterns due to multiple reflections in its interior. There is also, once again, evidence of the discrete nature of the resonators on the surface. Fig. 11(b) presents the IE-GSTC-SD model results. It should be noted that the susceptibility was extracted in a periodically infinite flat surface model [1], implying that some errors would eventually show up when applied to a curved surface. However, we see in Fig. 11(b) that the basic field structure is captured very well -all of the primary features of the interference pattern are predicted accurately. As previously, the introduction of quantization of the surface currents, shown in Fig. 11(c), creates some additional interference features as seen in the HFSS simulation. Finally, a comparison with the Dispersive zero thickness sheet model using IE-GSTC-SD with continuous sheet currents (10 div/λ). c) Dispersive zero thickness sheet model using IE-GSTC-SD but with discrete sources matching the number of wires (i.e., 1 div/cell). d) Fictitious non-dispersive zero thickness sheet with constant tangential surface susceptibilities. Simulation parameters: each side with N = 63 wires of radius d 0 = 0.2 mm, and length = 2.5 mm, separated by Λ = 2.15 mm, operating frequency 60 GHz and angle of incidence 60 • measured from x−axis.
non-spatially dispersive case in Fig. 11(d) shows a significant absence of the fine field structure, establishing the importance of capturing the spatial dispersion property of the unit cell.

VI. CONCLUSION
An IE-GSTC field solver to compute the scattered fields from spatially dispersive metasurfaces has been proposed and has been numerically confirmed using variety of examples. The work is a continuation of Part-1 [1], which proposed the basic methodology of representing spatially dispersive metasurface structures in the spatial frequency domain, k y . By representing the angular dependence of the surface susceptibilities in k y as a ratio of two polynomials, standard GSTCs have been extended to include the spatial derivatives of both the difference and average fields around the metasurface. These extended boundary conditions were successfully integrated into the standard IE-GSTC solver, which led to a new IE-GSTC-SD framework. The proposed IE-GSTC-SD platform has been confirmed using the semi-analytical Fourier decomposition method applied to uniform metasurfaces before testing it against a practical short conducting wire unit cell for various cases of finite size flat and curvilinear surfaces. All the results for finite-sized metasurface structures composed of spatially dispersive wire unit cells have confirmed the successful implementation and integration of spatial dispersion into the IE-GSTC simulation framework. Due to its inherent structural symmetry and simplicity, the wire unit cell has been proven to be an excellent example, where the proposed IE-GSTC-SD framework is tested to predict both high and low field amplitude features simultaneously and accurately. Moreover, the framework successfully modeled surfaces with conformal and curvilinear geometries, demonstrating its versatile architecture. Finally, some fine and subtle field interference patterns observed in HFSS have been traced to large unit cell periodicities in practical metasurfaces, which revealed those intricate interference features when accounted for in terms of discrete current sources. This work has so far focused on integrating tangential surface susceptibilities and spatially symmetrical structures to avoid cumbersome mathematical developments and implementation details. Moreover, the focus has been on uniform metasurface structures of relatively small unit cell periodicities for simplicity. A natural extension of this work is to incorporate normal surface susceptibilities and the bi-anisotropic tensor terms, which can model general non-uniform metasurface structures with arbitrary structural symmetries. In addition, a more in-depth analysis must be performed to rigorously explain the current quantization phenomenon observed in this work, which could be of greater importance for electrically large unit cell periodicities. This proposed work thus represents an important developmental step for fast and efficient simulation of practical metasurface structures which are not necessarily deeply sub-wavelength and exhibit fundamental spatial dispersion effects.