Surface Susceptibility Synthesis of Spatially Dispersive Metasurfaces for Space Compression and Spatial Signal Processing

An analytical method is proposed to synthesize the angle-dependent surface susceptibilities, <inline-formula> <tex-math notation="LaTeX">$\chi $ </tex-math></inline-formula>, of spatially dispersive (SD) or nonlocal zero-thickness metasurfaces. The proposed method is based on the extended generalized sheet transition conditions (GSTCs), whereby spatially dispersive metasurfaces are modeled using angle-dependent surface susceptibilities that take the form of rational polynomial functions of the transverse wave vector, <inline-formula> <tex-math notation="LaTeX">$k_{\parallel } $ </tex-math></inline-formula>. The suggested method derives the rational polynomial form of <inline-formula> <tex-math notation="LaTeX">$\chi (k_{\parallel })$ </tex-math></inline-formula>, which can then be expressed in the space-domain using spatial derivatives of the fields, resulting in a corresponding higher order spatial boundary condition to achieve the desired field operation. The proposed synthesis method is illustrated using variety of examples such as a space plate, spatial filters, and field absorbers, which are then validated using an integral equation (IE) solver, in which the corresponding higher order boundary conditions are integrated to predict the scattered fields. The proposed method thus not only represents a simple way to synthesize ideal zero-thickness metasurfaces but also helps establishes a way to define fundamental operational limits of spatially dispersive metasurfaces. This is illustrated by considering the space plate example and deriving the fundamental tradeoff between operation bandwidth and the achievable space compression.


I. INTRODUCTION
M ETASURFACES are 2-D arrays of subwavelength unit cells, where the macroscopic electromagnetic properties of the surface are determined by the microscopic properties of the unit cell.By tuning the unit cell [1], a variety of interesting wave transformations can be achieved.Because of Jordan Dugan, Tom J. Smy, and Shulabh Gupta are with the Department of Electronics, Carleton University, Ottawa, ON K1S 5B6, Canada (e-mail: JordanDugan@cmail.carleton.ca).
Francesco Monticone is with the School of Electrical and Computer Engineering, Cornell University, Ithaca, NY 14850 USA (e-mail: francesco.monticone@cornell.edu).
Color versions of one or more figures in this article are available at https://doi.org/10.1109/TAP.2024.3426072.
Digital Object Identifier 10.1109/TAP.2024.3426072 this, metasurfaces have found many applications from the RF spectrum to optics [2], [3], [4], [5], [6], [7], [8], [9].Generally, metasurfaces are designed by assuming the unit cells are deeply subwavelength.In this case, the induced polarization on the surface can typically be related locally to the fields at the surface.However, a general metasurface may be nonlocal, with induced polarizations at a given point related to the fields in an extended region along the surface.Such metasurfaces are also referred to as spatially dispersive.
Recently, there has been much interest in spatially dispersive metasurfaces due to their ability to manipulate signals in the spatial frequency domain.This allows for the implementation of signal processing operations in situ, at the speed of light and with lower power consumption than digital methods.For these reasons, spatially dispersive metasurfaces have the potential to find applications in high-speed image processing tasks such as edge detection [10], [11], [12], [13], [14], [15], [16], [17].It is possible to perform these operations using local metasurfaces and lenses in a 4f optical correlator configuration [18] but this requires a significantly larger length than that required of a single spatially dispersive metasurface.Another interesting application of spatially dispersive metasurfaces is the realization of space plates [19], [20], [21], [22].These surfaces emulate the effects of propagation through a much larger volume of free space, making them very useful in shrinking optical systems where lenses are required to focus light at a certain distance, such as imaging systems.Space plates may also find applications in radio frequency applications, such as in shrinking the size of lens-based antennas.
Spatially dispersive metasurfaces have, so far, been typically designed using first principles, using specific desired transfer functions expressed in the spatial frequency domain, k.This is limited to only handful of metasurface operations, and to date, there exist no method to systematically synthesize these surfaces for general wave transformations.Moreover, a lack of rigorous synthesis methods limits the understanding of spatially dispersive metasurfaces, their inherent tradeoffs, and practical limitations.Therefore, there is a clear need to devise rigorous synthesis methods to help design metasurfaces taking spatial dispersion (SD) into account.This problem is addressed in this work.
To this end, the metasurface synthesis problem can be very efficiently addressed using the effective medium description of spatially dispersive metasurfaces.Metasurfaces with subwavelength unit cells can be modeled as a zero-thickness sheet characterized by dipolar surface susceptibilities, χ, as their constitutive parameters [23], [24], where the fields around the surface must satisfy the generalized sheet transition conditions (GSTCs) [25], [26].The GSTCs combined with the surface susceptibilities form a boundary condition, suitable for either their integration in numerical methods or analytical synthesis, which can accurately predict the macroscopic fields away from the surface [23], [27], [28].
For deeply subwavelength unit cells, the surface is nonspatially dispersive and can be characterized using angleindependent1 susceptibilities (at least approximately, as some degree of angle dependency is always present, which, however, is often neglected or treated as an issue to mitigate).As cells become electrically larger (but still subwavelength with an upper limit of λ/2), spatial dispersion becomes dominant, and consequently, such a surface can no longer be accurately modeled with angle-independent surface susceptibilities.Recently, the extended GSTC method has been proposed to model such spatially dispersive metasurfaces with angle-dependent dipolar susceptibilities that take the form of rational polynomial functions in the spatial frequency domain [29].Furthermore, this method was integrated into an integral equation (IE) solver and a periodic Floquet solver yielding very accurate results for both uniform and nonuniform surfaces [29], [30], [31].
Therefore, the problem of synthesizing spatially dispersive metasurfaces may be rephrased as follows: What should be the form of effective surface susceptibility tensors and their dependence on the incidence angle, to achieve a specified wave transformation (i.e., generating desired transmission and reflected scattered fields in response to a given incident field)?This problem is thus addressed here in this article, where we propose a method to synthesize a set of angle-dependent surface susceptibilities to realize a desired reflection and transmission response, in addition to integrating those desired surface susceptibilities into an IE solver.A method to synthesize metasurfaces for a given response in the spatial frequency domain was proposed in [32]; however, this method synthesizes nonspatially dispersive surfaces with susceptibilities that are independent of the incident angle, thereby limiting the spatial-frequency-domain responses that can be achieved.
This article is structured as follows, Section II provides an overview of the extended GSTC model of spatially dispersive metasurfaces, as well as details of the proposed synthesis procedure to determine the required angle dependence to achieve a given wave transformation.Section III presents an application of the proposed synthesis procedure to a variety of examples, such as angular filters and wave absorbers.Section IV further demonstrates the usefulness of the proposed method, by applying this to the problem of space-plate design, by establishing some inherent fundamental tradeoffs in space plates, along with potential solutions.Finally, conclusions are provided in Section V.

II. PROPOSED SYNTHESIS METHOD A. Extended GSTCs
When modeling a metasurface, it is common to treat the surface as a zero-thickness sheet discontinuity in the electromagnetic field where the fields around the sheet must satisfy the GSTCs [25], [26] where ψ = ψ 2 −ψ 1 represents the difference in fields across the surface, ∇ T represents the component of the gradient operation that is tangential to the surface, and P and M are the electric and magnetic induced polarizations on the surface, respectively, For a nonspatially dispersive surface, the polarizations can be related to the local fields through surface susceptibilities as follows [23], [24]: where each of the tensors χα,β is (3 × 3) matrices containing both tangential and normal susceptibility components.Equations (1) and (2) combine to form the standard zero-thickness model of a nonspatially dispersive metasurface.For a spatially dispersive metasurface, the polarizations are now related to the fields in an extended region along the surface.This nonlocal dependence can be written as a convolution as follows [29]: where * represents spatial convolution.Equation (3) can be written in the spatial frequency domain through a Fourier transform, converting the convolutions into simple multiplication Comparing ( 4) with (2), we see that a nonlocal metasurface response leads to angle (k y )-dependent susceptibilities as opposed to constant susceptibilities in the case of no spatial dispersion.It was shown in [29] that a rational polynomial form appropriately represents each susceptibility component in the spatial frequency domain.This formulation essentially takes a Pade approximation of the metasurfaces' response in the spatial frequency domain, which is very useful as it allows us to capture poles and resonances in the spatial frequency domain, and it also allows us to conveniently represent the metasurface in the spatial domain with a set of differential equations.
Consider next a metasurface in the yz plane containing only tangential susceptibilities, illuminated with TE polarization.Let us assume for simplicity that such a surface can be modeled by χ zz ee and χ yy mm susceptibilities, which would take the rational polynomial form as follows: For this simplified case, (4) reduces to Substituting the angle-dependent susceptibilities into (6) yields on which a spatial inverse Fourier transform can be applied, resulting in the following differential equation: This equation can now be combined with the standard GSTCs (1) to form the extended GSTCs Here, we see that the extended GSTCs take the form of a set of differential equations relating the discontinuity of the fields across the surface to the average fields at the surface and involves spatial derivatives of the fields, essentially modeling the metasurface using higher order boundary conditions [33], [34].These equations can then be integrated into numerical methods, such as an IE solver, as was done in [30], to efficiently simulate electromagnetic scattering from these structures.

B. Synthesis of Spatially Dispersive Metasurfaces
Now we will explore how we can synthesize a set of susceptibilities to produce a desired reflection R(k y ) and transmission responses T (k y ) in the spatial frequency domain.Again, we will consider the case of TE incidence and tangential-only susceptibilities, for sake of illustration and simplicity, allowing us to only consider χ zz ee and χ yy mm susceptibility components.While this may seem limiting, a wide variety of desired reflection and transmission responses in the spatial frequency domain can be achieved with just these two susceptibilities.Moreover, many practical structures, such as all-dielectric Huygens' metasurfaces, can be accurately modeled with these two susceptibility components [29].Here, we extend the method proposed in [24] by writing the following relationship between the spatial frequency components of the fields and susceptibilities: If we normalize the incident field illuminating the surface at the spatial frequency of interest to unity, we can expand that difference and average fields in terms of the desired R(k y ) and T (k y ) as follows: where η 0 is the impedance of free space.Substituting (11) into (10) yields an expression for the spatial-frequency-domain surface susceptibilities in terms of the desired reflection and transmission profile Examining ( 12), we see that the susceptibilities will not generally be in rational polynomial form of ( 5) and as a result do not have a convenient representation in the space domain which cannot be integrated directly into the IE solver as in [30].The use of the rational polynomial form of the susceptibilities, in the prior work of [30], was physically motivated by modeling the metasurface unit cells' temporal response by a sum of Lorentzian resonators, with the resonance properties being a function of the transverse wave vector, as further discussed in [29].Therefore, we approximate the susceptibilities as rational polynomial functions so that they can be incorporated into the IE solver as a higher order boundary condition, and second, to express them in a form that lends itself to realizable physical unit cells incorporating Lorentzianlike resonators, ensuring the stability and causality of the realized metasurface.We therefore express the relationships in (12) using where the coefficients a n -d n can be determined by taking a Taylor series or Pade approximation of any nonpolynomial terms that may appear in the numerator or denominator of ( 12).This synthesis procedure can yield susceptibilities that are active or lossy.If one wishes to realize a metasurface using passive and lossless unit cells, the phase difference between R and T response must be 90 • as shown in the Appendix.It is further shown in the Appendix, that the Pade approximant susceptibilities will also be passive and lossless.Once the susceptibilities are synthesized, we can validate the synthesis by calculating the synthesized reflection and transmission coefficients as follows [27]: In principle, the synthesized spatial-frequency-domain susceptibilities can be rational polynomials of any order which can lead to higher order derivatives in the space domain.However, for numerical and physical considerations, it is useful to express higher orders using a summation of second-order functions (without any loss of generality).We thus decide to express the induced polarization as follows [29]: where the general response has been decomposed into a sum of Lorentzians, and a similar expression can be derived for χ yy mm .With this form, the extended GSTCs can now be implemented into a numerical solver with derivatives of no higher order than two.Any given rational polynomial susceptibility in the spatial frequency domain can be written in this form using a partial fraction expansion.It should be noted that while we are modeling these SD surfaces as zero-thickness sheets, realizing a general nonlocal response will require some thickness as pointed out in [22], [35], and [36].However, even if the realization of the susceptibilities requires some significant thickness, the extended GSTCs will still correctly predict the reflection and transmission response of the surface.Furthermore, a general nonlocal response could also be realized by cascading near-zero-thickness surfaces.In Sections III and IV, we will use this synthesis method to design and analyze some interesting and practically relevant spatially dispersive metasurface structures.

III. SYNTHESIS EXAMPLES A. Image Processing Filter
Here, we will consider the synthesis of a spatial frequency filter.Such devices have recently gained much interest in the field of wave-based analog computing, with applications such as high-speed image processing [10], [11], [12], [13], [14], [15], [16], [17].In the absence of spatial dispersion, these filters can be implemented through a 4f optical correlator [18], requiring two lenses in addition to a nondispersive metasurface or phase/amplitude mask, with the overall size of the correlator being equal to four times the focal length of the lens.Spatial dispersion allows for the ability to perform spatial-frequencydomain manipulation of fields without having to take a Fourier transform, negating the need for the lens.As a result, the spatial filtering response can, in principle, be implemented using only a single spatially dispersive metasurface, greatly reducing the size requirement of the system (albeit still subject to certain fundamental limitations when physically implemented, especially at optical frequencies [35]).
To illustrate this, we will consider a spatial filter metasurface that performs basic image-processing tasks.The surface will be synthesized such that the transmitted image will be a blurred version (low-pass-filtered) of the input image, and the reflected image will be an edge-detected (high-pass-filtered) version of the input image.One such simple form of a reflection and transmission characteristic of such a surface is as follows: where θ is the angle of incidence.It is easy to see that the transmission response will give T ≈ 1 for incidence angles near normal, and T ≈ 0 for grazing incidence confirming its low-pass characteristic.Conversely, we have R ≈ 0 for angles near normal and R ≈ 1 for grazing angles.Substituting these reflection and transmission coefficients into (12) yields the following expression for the surface susceptibilities required to generate this response: It is worth noting that this filter as synthesized does not satisfy the passivity requirements outlined in the Appendix, and therefore would require active elements in its unit cells.However, the surface can be made to be passive by adding a Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.• phase shift to either R or T response and resynthesizing the surface susceptibilities.To implement this surface into the IE solver in [30], we need to take a Pade approximant of these susceptibility expressions so they are in rational polynomial form.For this example, we take a fourth-order Pade approximant of (17), assuming the surface will be designed to operate at 60 GHz (arbitrarily chosen), and integrate them into the IE solver.
The simulation results are given in Fig. 1.The surface is illuminated with a Gaussian beam with a beamwidth of λ/2 as shown in Fig. 1(a).The scattered fields (reflected and transmitted) are shown in Fig. 1(b).We see in the reflected fields that the wave vector components near normal incidence (fields close to the center of the beam) have been significantly attenuated.Similarly, in the transmitted field, we see that the wave vector components of the field near normal incidence have been preserved, but the components at high incidence angles (fields close to the edge of the beam) have been attenuated.To further verify this, we can take a Fourier transform of the fields at the surface, at x = 0 ± , to determine the magnitude of the reflection and transmission coefficients of the surface.The results are reported in Fig. 1(c) and (d), where we see that the response of the synthesized surface in the IE solver very closely matches the desired surface response.These results thus confirm the operation of the synthesized filter.We note that while this example takes a very simple form of the spatial filter for the sake of demonstration, more advanced filters can easily be synthesized with a wide variety of reflection and transmission characteristics in both the passbands and stop bands of the filters (such as controlling the filter cut-offs, rolloffs, and ripple controls using Butterworth or Chebyshev filter transfer functions, for instance) using a very rich literature on designing microwave filters [37].

B. Perfect Angle-Independent Absorber
Another interesting application of spatially dispersive metasurfaces is the synthesis of a perfect angle-independent absorber.Such a surface would ideally produce no reflection or transmission over the entire angular range, thus acting as an ideal wave absorber.In contrast, a local metasurface would be able to act as a perfect absorber only at specific incidence angles, limiting its overall angular performance.A spatially dispersive surface owing to its angle-dependent susceptibilities can lift this limitation.To synthesize such a surface, we substitute R = T = 0 into (12) yielding the following form of the susceptibilities: Again, these are not suitable for integration in the IE solver.However, we arrive at very convenient expressions for the approximate susceptibilities using the Taylor expansion of the square root term, for instance using second-order approximation.It is clear that the order of the Taylor expansion will determine the angular operation range of the absorber.These susceptibilities were next integrated into the IE solver and used to analyze a surface at 60 GHz.The surface is illuminated by a Gaussian beam with a waist of 0.2λ.Such a narrow Gaussian will have a very wide angular range, and as a result, any absorber will have to be well-matched at all the angles.To demonstrate this, we attempt to use a nondispersive surface (which is perfectly absorbing only for normal illumination) with χ zz ee = χ yy mm = −2 j/k 0 and the spatially dispersive surface with the susceptibilities given by (19).
The scattered fields generated by the nonspatially dispersive metasurface show a strong angular scattering as seen in Fig. 2(a) and absorption near normal incidence angle as expected.On the other hand, the scattered fields in the spatially dispersive case are strongly attenuated at all the angles as seen in Fig. 2(b), with absorption levels < −50 dB, in contrast to about −20 dB in the local case.From these results, we see that the spatially dispersive absorber clearly outperforms the nondispersive absorber.To further demonstrate this, we plot the magnitude of the reflection and transmission coefficients in the spatial frequency domain for each surface as shown in Fig. 2(c).These results further confirm that the spatially dispersive metasurface vastly outperforms the nondispersive Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.metasurface over the entire angular range with the scattered fields produced by the spatially dispersive metasurface being about 20 dB or more lower than the nondispersive surface across all the angles.

IV. SPACE-PLATE EXAMPLE
The susceptibility synthesis method as presented and demonstrated here for the two example cases not only provides a systematic procedure to design spatially dispersive metasurfaces but also represents a useful tool to gain physical insights and derive associated performance tradeoffs and key characteristics.In this section, using a third demonstration example of a space plate, we elucidate on some of its important features and then derive key characteristics.

A. Space-Plate Synthesis
Space plates have recently emerged as a means of shrinking optical systems by emulating the propagation through some length of free space ℓ using a structure with a thickness δ < ℓ.This is illustrated in Fig. 3, where a free-space length ℓ is ideally compressed to zero, using a zero-thickness sheet metasurface in such a way that the propagating fields are preserved outside the compression region.Metasurfaces, being electrically thin, represent a good candidate for the implementation of a space plate.Several metasurface-based space plates have been reported in the literature [19], [20], many with structures suitable for modeling with the extended GSTC method.
An ideal space plate should have zero reflection and a transmission response with the following form in the spatial Phase response of a free-space body of length ℓ and of the synthesized space plate in the spatial frequency domain, where the space plate is represented as an equivalent zero-thickness sheet model (a physical implementation of this model will necessarily have a nonzero thickness and a finite compression ratio).frequency domain: where the response mimics free space by applying a phase shift to an incoming field that is equal to the phase shift it would experience had it propagated over a distance ℓ.The transfer function of (20) reveals that |T (k y )| = 1, i.e., transmission magnitude is always unity at all the incidence angles, across the angular spectrum, and only the transmission phase is varied with k y .Such form is well-known as an all-pass response in conventional filter theory.This allows us to write (20) in the following form: where This form automatically ensures that |T | = 1 for any ζ .Moreover, ζ can be easily expanded using a Taylor series of a tan(•) function up to any order leading to a polynomial form of T (k y ), as desired for the synthesis.This realization that a space plate is an all-pass spatial filter, thus requiring zero-backscattering and an engineerable transmission phase, further reveals that it can be physically implemented as Huygens' metasurface.Huygens' metasurface is composed of a subwavelength scatterer formed using collocated orthogonal electric and magnetic dipole moments [38], [39], [40].When the two moments are properly balanced, both in phase and magnitude, a zero backscattering is achieved.Consequently, considering a space plate does not rotate the field polarization and can be implemented using Huygen's configuration, only two surface susceptibility components are sufficient to realize a space-plate response: χ zz ee and χ yy mm for the wave propagation along z.Therefore, substituting ( 21) into ( 12) yields the following two expressions of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.required susceptibilities: However, this is not yet suitable for integration with the extended GSTCs.To remedy this, we must approximate these susceptibilities as a rational polynomial.To do this for an arbitrary space plate is a difficult task as the presence of the tangent function and the radical function in (22) will place poles on the k y -axis.Any poles in the neighborhood of |k y | ≤ 1 will have a significant effect on the angular response.As the space-plate effect increases, the number of poles in this range will also increase, resulting in a higher order rational polynomial being needed to achieve a space-plate response over a given angular range.While this is not a fundamental problem per se, it does makes the integration of the corresponding spatial boundary conditions involving higher spatial derivatives into IE solver more difficult.
As an example, we consider a space plate operating at a frequency of 60 GHz, intended to replace a free-space length of say ℓ = 0.01 m.We will synthesize two surfaces, based on second-and fourth-order Pade approximations of (22), respectively.The ideal and synthesized susceptibilities are plotted in Fig. 4(a) and (b).From Fig. 4(a) and (b), we see that the ideal susceptibilities contain four poles along the k y -axis.To determine the quality of the synthesis, we can calculate the synthesized reflection and transmission coefficients using (14).The phases of the transmission coefficient coefficients are plotted in Fig. 4(c), and the magnitudes of the reflection and transmission coefficients are plotted in Fig. 4(d).We see from these results that the second-order Pade approximant reproduces the phase response of the space plate accurately near normal incidence and diverges as the incidence angle is increased.We observe the same trend for the fourth-order Pade approximant but, as expected, the fourth-order approximation matches over a larger angular range.Therefore, as expected a better space-plate response can be obtained with higher order polynomial representation of susceptibilities.
With the susceptibilities synthesized, we can now simulate the space plate with the full-wave IE solver from [30].Fig. 5(a) shows a quasi-Gaussian incident field at 60 GHz with a beam waist at x = 0.03 m.The beam is not an ideal Gaussian as it has a fairly wide angular range, resulting in some distortion due to nonparaxial effects.Fig. 5(a) also shows the total fields if a second-order metasurface space plate is placed at x = −0.03m.Here, we see that the beam waist has been shifted to x = 0.02 m as expected, but the shape of the field has been slightly distorted near the beam waist.This can be attributed to the fact that the angular spectrum of the incident field is greater than the angular range in which the second-order metasurface behaves as a space plate.In Fig. 5(a) further shows the total fields using the fourthorder space plate, where the beam waist has also shifted to x = 0.02 m, but this time the shape of incident field is much better preserved, as the fourth-order metasurface functions as a space plate over a larger angular range.To further verify this phase compression behavior, we can determine the phase response of each metasurface in the spectral domain by taking a Fourier transform of the fields on the surface.This phase is plotted, along with the ideal space-plate phase response in Fig. 5(b) where we see very excellent agreement with the expected spectral response of Fig. 4(c).

B. Performance Tradeoffs of Space Plates
In Section IV-A, we considered the synthesis of an ideal space plate from a purely mathematical point of view.While this yields a set of susceptibilities capable of reproducing the space-plate response, it does not give us direct insight into how a space plate might be implemented and behave from a physical point of view.Here we consider how one might implement a space plate using a spatially dispersive Huygen's surface.This derivation closely follows the one performed in [19] and [20], adapted here for Huygens' surfaces characterized by surface susceptibilities.
Huygen's surface is a metasurface with a collocated electric and magnetic resonance response allowing for zerobackscattering due to the destructive interference between backscattered waves from the induced electric and magnetic dipole moments.In the temporal frequency domain, the response of a Huygen's cell can often be modeled as a Lorentzian oscillator with spatial dispersion manifesting as an angular dependence on the Lorentzian properties [29].Thus, the susceptibilities can be expressed as follows: where ω 0 is the resonance frequency, ω p is the plasma frequency, γ is the damping coefficient, and A is the amplitude.If we consider a surface where χ zz ee = χ yy mm , then we will have R ≈ 0. This is only strictly true for normal incidence but the backscattering remains small over a wide angular range. 2 We 2 This can be mitigated by assuming the matching condition at all the angles using the spatial dispersion of the surface such that χ zz ee /χ Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.can consequently write the transmission coefficient in terms of the Lorentzian parameters as follows: If we assume a low-loss structure such that ωγ ≪ Ak 0 ω 2 p , then we can write the phase of the transmitted field as Near resonance we have ω ≈ ω 0 and so we can approximate the phase of the transmitted field as follows: Finally, if the unit cell is spatially symmetric, then near the resonance we can model the Lorentzian parameters as second-order polynomials of k y [29] If we design the surface such that β 2 = 0, making ω p constant with respect to spatial frequency, we can express the phase of the transmitted field as follows: Given that the phase of an ideal space plate is we can conclude by comparison that Huygen's surface will operate as a space plate with an effective length of This equation shows that the effective length of the space plate is determined by the strength of nonlocality (ζ 2 ) and the bandwidth of the resonance which is inversely proportional to ω p .This means that there is a tradeoff between the bandwidth within which the surface will operate as a space plate and the effective length of the space plate consistent with [19], [20], and [22].We can further observe a similar tradeoff between the effective length of the space plate and its angular range.This can be obtained by considering the transmission phase φ(k y = 0) ≥ −2π at OFF-resonance frequency ω.Since the value of φ(k y ) is bounded between 0 and −2π , the second term of ( 27) must be bounded between 0 and 2π, i.e., This gives the maximum value of k y and thus the angular range as which using (29) can also be written as This implies that the maximum angular range of operation of a space plate is inversely proportional to the space-compression length, i.e., the stronger the space compression, the smaller the angular range of operation of the space plate.Finally, rearranging (32) yields Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.which is precisely twice the angular range of space plates provided in [19].The reason for the increase in the angular range can be attributed to the fact that Huygen's structure contains two resonances, an electric and magnetic resonance, whereas the derivation in [19] was performed assuming a single resonance.This analysis thus strongly highlights an ability of a Huygens' surface to increase the spatial frequency bandwidth of a space plate, thanks to the deeper insight gained using the proposed surface susceptibility description of the surface.

C. Cascaded Huygens' Metasurfaces
Furthermore, following the work in [19], the space-plate limitations elucidated above can be relaxed by implementing the space plate using a serial cascade of metasurfaces.This allows us to use metasurfaces with lower space compression to still achieve a large effective space compression as a result of the cascade, while maintaining the angular range and the bandwidth.Using (33), for a fixed angular range corresponding to k y,max , the compression length ℓ ′ = ℓ/n, so that where n is the total number of metasurfaces used in the cascade.The number of metasurfaces n thus represents another degree of freedom that can thus relax the constraint of (32).Again, (34) is exactly twice the angular range reported in [19].
We note that since we are assuming Huygens' metasurface with near zero backscattering across all angles, such a cascaded configuration can easily be analyzed assuming no electromagnetic couplings between consecutive surfaces.Next we will demonstrate Huygens' space plate using the IE solver.To do this, the Lorentzian parameters need to be mapped to the angle-dependent susceptibilities.This can be done as follows [29]: 0 where these parameters can be substituted into (15) with N L = 1 and χ zz ee 0 = χ zz mm 0 = 0. Again for simplicity, the matching condition is used for normal incidence only.For this example, let us synthesize a surface that operates at 60 GHz with an effective length of say ℓ = 0.02 m.Choosing ω 0 = ω p = 2π 60 × 10 9 rad/s and α = 1000 rad/s, we can determine from (29) that ζ 2 must be equal to 7.11 × 10 16 rad/s/m 2 .This surface was excited with a 60-GHz Gaussian beam with a waist of 1.5λ located at x = 0.03 m as shown in Fig. 6(a).
We next run a series of tests to understand the efficacy of the space plates under various conditions.In Fig. 6(b), we show the total fields as a result of the interaction of the incident field with the Huygens' surface, calculated using the IE method.Here, we see the beam waist has shifted 0.02 m Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.from x = 0.03 m to x = 0.01 m, confirming that the space plate is behaving as expected.Furthermore, we observe no distortion in the shape of the Gaussian beam, suggesting that the incident field has no significant spectral components outside of the allowed angular range of this surface.Now we consider a surface with double the effective length, ℓ = 0.04 m, which is achieved by setting ζ 2 = 1.4212 × 10 17 rad/s/m 2 .The total fields for this surface are given in Fig. 6(c) and we see that while the beam waist is in the expected location of x = 0.01 m, there is slight but visible distortion in the shape of the beam.This is because an increase in the effective length of the surface has caused the usable angular range to decrease as predicted in (33).However, we can implement a space plate with an effective length of ℓ = 0.04 m using two metasurfaces with an effective length of ℓ ′ = 0.02 m.The total fields for this configuration are given in Fig. 6(d) and (e) for the cases where the two metasurfaces are separated by a small and a larger distance, respectively.In both the cases, we see the beam waist is in the expected location of x = 0.01 m and there is no visible distortion to the shape of the beam.This is because the addition of an extra surface has allowed for the angular range of our space plate to double as predicted by (34).Furthermore, we see that the space-plate operation is independent of the separation between the metasurfaces, which can also be attributed to the perfectly matched characteristics of the two Huygens' surfaces used in this demonstrated, effectively decoupling them from each other.Therefore, a designer should seek to make the distance between the metasurfaces as small as practically possible to maximize the compression ratio.
Finally, in Fig. 6(f) we add an additional metasurface to increase the effective length ℓ up to 0.06 m.In this case, the effective length becomes so large that the beam waist is no longer visible.However, a strong space compression effect is clearly visible.Therefore, in principle, any number of metasurfaces can be cascaded to achieve any desired effective length.However, this will increase the thickness of the space plate due to the the free-space separation between consecutive surfaces.The magnitude of this separation is determined by the free-space matching characteristic of the surface, the finite thickness of the surface, and the period of the metasurface (which will determine the evanescent field couplings between the surfaces) and potentially lowers the effective compression ratio achievable in practice.

V. CONCLUSION
In this article, a method to synthesize the angle-dependent susceptibilities of spatially dispersive, or nonlocal, metasurfaces has been proposed and full-wave-verified by integrating the synthesized susceptibilities into an IE-GSTC solver.The method is based on extended GSTCs which result in the needed higher order boundary conditions to achieve a given wave transformation operation on given incident fields.Using a variety of examples including spatial filtering and spatially broadband wave absorption, the proposed synthesis method has been demonstrated and confirmed using the IE solver.
A particular focus on space plates has further been given to emphasize the usefulness of this method.The proposed method enables a fundamental description of the spatially dispersive metasurfaces in terms of idealized surface susceptibilities, as well as providing deeper physical insights in terms of Lorentz oscillator models of the susceptibilities.This is apparent when the need of a Huygens' metasurface became clear to implement an ideal space plate exhibiting no back reflection.Such an implementation was found to provide more favorable constraints between space-compression ratio and the maximum angular operation range, compared with a singleresonator implementation.Future works will focus on practical implementation of these synthesized susceptibilities, along with engineering practical Huygens' resonators for space-plate applications to achieve desired space compression.A preliminary progress was recently shown in [41] where it was found that an equivalent spatially dispersive response can be realized through a nonspatially dispersive periodic metasurface, establishing a clear conceptual link between angle-dependent susceptibilities and leaky wave antennas/photonic crystal slabs.This may allow for the use of the rich literature on these subjects to aid in the realization of metasurfaces with angle-dependent surface susceptibilitites.Furthermore, while only uniform surfaces have been considered in this work, a natural extension is to introduce nonuniformity of the surface to provide an extensive design freedom to achieve more complex wave transformations using spatial dispersion.

APPENDIX PASSIVE AND LOSSLESS SD METAURFACES
Passivity can be ascertained by considering the metasurface as a two-port reciprocal and symmetric network with S-parameters defined as follows: A network defined by its S-parameters is lossless, and by extension passive, if S is a unitary matrix (I − SS * = 0).It is easy to show that With this, we can see that the S matrix will be unitary if the following two conditions are satisfied: where the second condition is satisfied if either R or T is zero or if R and T have a π/2 phase difference.The question now is will the passivity be preserved when we take the Pade approximant of the susceptibilities?To answer this question, we used the well-known property of passive and lossless susceptibilities, that is, χ T ee = χ * ee χ T mm = χ * mm .For the susceptibility components we are considering, this condition reduces to the following: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
Therefore, if we attempt to synthesize a passive and lossless surface, our desired susceptibilities will be purely real functions of k y .A Pade approximant of a real-valued function will always yield a real-valued function.Therefore, the approximated susceptibilities will be passive as well.

Manuscript received 24
August 2023; revised 9 April 2024; accepted 30 June 2024.Date of publication 16 July 2024; date of current version 9 August 2024.The work of Francesco Monticone was supported by the U.S. Air Force Office of Scientific Research through Dr. Arje Nachman under Grant FA9550-22-1-0204.(Corresponding author: Jordan Dugan.)

Fig. 1 .
Fig. 1.IE simulation of spatial filter metasurface.(a) Incident 60-GHz Gaussian beam with λ/2 waist located at x = 0. (b) Scattered fields produced by the synthesized metasurface.(c) Calculated reflection coefficient magnitude in spatial frequency domain compared with the ideal case.(d) Calculated transmission coefficient magnitude in spatial frequency domain compared with the ideal case. 90

Fig. 2 .
Fig. 2. IE simulation of metasurface perfect absorber assuming an incident 60-GHz Gaussian beam with 0.2λ beam waist located at x = 0. (a) Scattered fields produced by nonspatially dispersive absorber.(b) Scattered fields produced by spatially dispersive absorber.(c) Calculated magnitude of the reflection and transmission coefficients for both the absorber metasurfaces in the spatial frequency domain.

Fig. 3 .
Fig. 3.Phase response of a free-space body of length ℓ and of the synthesized space plate in the spatial frequency domain, where the space plate is represented as an equivalent zero-thickness sheet model (a physical implementation of this model will necessarily have a nonzero thickness and a finite compression ratio).

Fig. 4 .
Fig. 4. Metasurface space plate syntheses.(a) Ideal and synthesized electric susceptibility in the spatial frequency domain.(b) Ideal and synthesized magnetic susceptibility in the spatial frequency domain.(c) Phase response of ideal and synthesized metasurfaces in the spatial frequency domain.(d) Magnitude of the reflection and transmission coefficients of the synthesized space plate in the spatial frequency domain.

Fig. 5 .
Fig.5.IE simulation of synthesized metasurface space plates.(a) Incident quasi-Gaussian beam at 60 GHz with a 1.0λ beam waist at x = 0.03 m and the total fields for simulation with second-order and fourth-order metasurfaces, respectively, with ℓ = 0.01 m.(b) Transmission phase response of the metasurfaces in the spatial frequency domain predicted by the IE method compared with the ideal space-plate response.