Sideways adiabaticity: Beyond ray optics for slowly varying metasurfaces

Optical metasurfaces (subwavelength-patterned surfaces typically described by variable effective surface impedances) are typically modeled by an approximation akin to ray optics: the reflection or transmission of an incident wave at each point of the surface is computed as if the surface were"locally uniform", and then the total field is obtained by summing all of these local scattered fields via a Huygens principle. (Similar approximations are found in scalar diffraction theory and in ray optics for curved surfaces.) In this paper, we develop a precise theory of such approximations for variable-impedance surfaces. Not only do we obtain a type of adiabatic theorem showing that the"zeroth-order"locally uniform approximation converges in the limit as the surface varies more and more slowly, including a way to quantify the rate of convergence, but we also obtain an infinite series of higher-order corrections. These corrections, which can be computed to any desired order by performing integral operations on the surface fields, allow rapidly varying surfaces to be modeled with arbitrary accuracy, and also allow one to validate designs based on the zeroth-order approximation (which is often surprisingly accurate) without resorting to expensive brute-force Maxwell solvers. We show that our formulation works arbitrarily close to the surface, and can even compute coupling to guided modes, whereas in the far-field limit our zeroth-order result simplifies to an expression similar to what has been used by other authors.


Introduction
Optical metasurfaces, subwavelength structures described by an effective sheet impedance [25,20,18,19,43,42,17,16,38,2], are now being designed for large-area optical devices using models in which the far-field reflection/transmission coefficients are computed at each point assuming a uniform (or periodic) surface-as explained below, we refer to these as "ray-optics" models.This is a good approximation for surfaces (or unit cells) that are varying slowly, a fact that is closely connected to the "adiabatic theorem" [23,21] for waves propagating through slowly varying media.However, although there are countless papers and books on modeling propagation through slowly varying media [26,41,23], exploiting the rate of change as a small parameter ε, the "sideways" problem of scattering off a slowly varying surface (Fig. 1) is relatively unstudied.In this paper, we address the following key questions: how quickly does the ray-optics approximation converge as ε → 0, can we quickly compute the low-order corrections (both to improve accuracy and to validate ray optics), and how do we compute both far-field and near-field scattering (e.g.coupling to guided modes)?A typical metasurface has two scales: the subwavelength scale of the microstructure and the macroscale of the nonuniformity.In this paper, we address corrections due to the macroscale nonuniformity, which we allow to be completely arbitrary, while we follow other authors [25,20,18,19,43,42,17,16,38,2] in subsuming the microscale into effective surface impedances.
In particular, we use the technical machinery of surface-integral equations (SIEs) [13,29] and a "locally uniform approximation" of the metasurface to show that the ray-optics approximation is the far-field zero-th order term in a convergent series (Section 4 and Appendix C), that each successive correction can be computed simply by performing integrals (not by solving any PDE or other system of equations), and that the next-order correction scales as ε 2 .Moreover, our series allows us to compute the full Green's function of the surface: the fields in response to arbitrary sources or incident fields, including the near-field terms (fields and/or sources close to the surface).We show that these near fields allow us to compute the coupling of an incident wave to guided modes on the surface [27,42] and that they also appear in the zeroth-order locally uniform approximation.For rapidly varying metasurfaces, such as those designed to reflect light at a very oblique angle [46], we show that even the far-field accuracy is substantially improved by including the next-order correction.Perhaps more importantly, the ability to compute the next-order correction provides a way to validate a ray-optics design for very large-area metasurfaces, where brute-force Maxwell simulations are impractical and there was previously no way to evaluate the ray-optics accuracy short of a laboratory experiment.
Since typical metasurface designs lead to large computational domains (often hundreds of wavelengths [24]) that are intractable by standard simulation techniques, e.g.finite-difference and finite-element methods, previous work on metasurfaces has made extensive use of numerical simulations based on ray-optics approximations.In particular, authors typically compute reflection/transmission coefficients for periodic surfaces with a variety of unit cells, they assume that these coefficients remain accurate even for an aperiodic surface where each unit cell is different, and then they select the unit cell at each point on the surface to achieve a desired optical functionality [44,37,2,3,48,49,47,36]. (For subwavelength unit cells where there is only a single "specular" reflected/transmitted wave, the reflection/transmission coefficients can also be fitted to an effective sheet impedance, giving a "homogenized" effective medium at each point [30,40,39,2].)Because these works described the surfaces by a single far-field (planewave) reflection and/or transmission coefficient at each point, they can be thought of as "ray-optics" approximations even if they were expressed in the language of wave optics.(A closely related approximation-curved surfaces treated as locally flat-is called a "tangent-plane" or "Kirchhoff" approximation [45].Yet another closely related approximation is provided by scalar diffraction theory [32].)Here, we assume that the metasurface is subwavelength enough to be described as an effective sheet impedance at each point, but we do not only compute the scattering assuming that the impedance is locally uniform: our goal is to take the macro-scale spatial variation (the aperiodicity) explicitly into account by computing correction to the locally uniform approximation.(Potential extensions to slowly varying periodic structures where the micro-scale is treated explicitly, perhaps to include additional diffraction orders for large-period structures, are discussed in Section 5.) Wave propagation through slowly varying media is usually treated by coupled-mode theory: one expands the wave in the basis of "instantaneous" [11] eigenfunctions of each cross-section [26] or period [21], and then obtains a set of coupled differential equations in the mode coefficients (typically truncated to only a few guided modes).As the medium varies more slowly, the coefficients tend to constants, corresponding to "adiabatic" transport of modes without inter-modal scattering [26].Unfortunately, this approach appears awkward to apply to the problem of scattering off of a slowly varying medium, both because there is a continuum of radiating modes and because one wants to describe the basis of incoming/outgoing planewaves independently of the varying surface in order to connect to a ray-optics (zero-th order) approximation.The fact that an incident planewave is not an eigenfunction of the cross-section at each point means that one cannot simply quote the standard adiabatic theorem to justify the metasurface ray-optics approximation, for example.Another type of technique for approximating the scattering from a weakly perturbed surfaces is a Born approximation [10,41] (also known as a volume-current method [22], Kirchhoff approximation [45], etc.), which handles perturbations like surface roughness that are small in amplitude but not necessarily slowly varying, whereas our goal is to handle variations that are slow but not small.
A surface impedance Z(εx) (defined precisely in Section 2) varies more and more slowly as ε → 0, and our goal is an expansion with terms proportional (in a certain norm) to powers of ε [13,29].This expansion is achieved through an SIE.In particular, since the media above and below the surface are homogeneous, we express the problem in terms of an SIE in which the unknowns reside only on the surface.Our derivations start with an approximate Green's function G p that is a building-block of the locally uniform approach (Section 4), and then we insert this into an exact SIE, obtained by enforcing to transition conditions on the metasurface, to derive a series of corrections (Appendix C)-like a Born-Dyson series [5], the corrections are expressed in terms of integrals involving G p .These integrals must be computed numerically by a "quadrature" technique [14] (Appendix G), but such computations are simple summations on a computer that are far easier than solving the large systems of equations arising in brute-force computational methods, and also have the advantage of parallelizing perfectly (fields at different points can be computed completely independently).Truncating to the zeroth-order term in the series does not correspond to setting ε = 0 (a uniform surface), so even the lowest order locally uniform approximation captures to certain extent the surface variation.In the far field, G p simplifies to an expression G ff p that can be written in closed form (eliminating an integral), recovering the usual ray-optics approximation at zero-th order (Appendix F).Using this approach, we demonstrate through numerical experiments that the ray-optics approximation (i.e., the far-field of the locally uniform approximation) produces far-field errors that vanish as ε 2 , and more generally as ε 2N +2 if we include N th-order corrections (Fig. 5).In the presence of guided modes, which correspond to poles that appear in G p at certain wavevectors [10], we show that this also simplifies the integrals in our perturbative expansion (via a steepest-descent approximation) if one is mainly interested in coupling to guided modes (Appendix F and Figs.6 and 10).

Problem formulation
We consider a metasurface Γ in two spatial dimensions that divides the xy plane into two unbounded half-planes, Ω + = {y > 0} and Ω − = {y < 0}, occupying the regions above and below Γ, respectively.The media Ω + and Ω − surrounding the metasurface are assumed to be homogeneous with electric permittivity and magnetic permeability denoted by 0 > 0 and µ 0 > 0, respectively (Fig. 1).
The metasurface is characterized by the so-called generalized sheet transition conditions [25]: which for the sake of presentation simplicity are assumed to be given in terms of scalar quantities corresponding to surface impedance Z and the surface admittance Y [17].The symbol n in (1) denotes the unit normal vector to Γ = {y = 0} pointing upwards, and E ± and H ± (resp.E ± and H ± ) denote the tangential (resp.entire) fields at y = 0 ± .A time-harmonic incident field u inc with angular frequency ω > 0 impinges on a metasurface with surface impedance Z and surface admittance Y , producing a reflected field propagating in the upper half-plane and a transmitted field propagating in the lower half-plane.The metasurface is surrounded by a homogeneous medium with electric permittivity 0 > 0 and magnetic permeability µ 0 > 0.
It thus follows from Maxwell's equations that in E z polarization the total electromagnetic field (E, H) is given by E = (0, 0, E) and H = (H 1 , H 2 , 0), and can be obtained from the z-component of the electric field by means of the relations: where k = ω √ 0 µ 0 is the wavenumber of the surrounding media, ω > 0 is the angular frequency, and η = µ 0 / 0 is the intrinsic free-space impedance.Similarly, in H z polarization it holds that E = (E 1 , E 2 , 0) and H = (0, 0, H) where Relations ( 2) and (3), on the other hand, yield that the transition conditions (1) can be equivalently expressed as in E z polarization, and in H z polarization, where the notations have been introduced to refer to the jump and the sum of a scalar field u across Γ, where u + (resp.u − ) denotes the limit value of u on Γ from Ω + (resp.Ω − ).
In order to treat both E z -and H z -polarization cases, we define the metasurface parameters α and β as α = η 2Z and β = 2ηY (6a) in E z polarization, and by in H z polarization.Throughout this paper we assume that α and β are continuous complexvalued functions that satisfy Re α ≥ 0 and Re β ≥ 0, which correspond to assuming that both the surface impedance Z and the surface admittance Y are passive but not necessarily lossless.We will eventually consider these quantities to be slowly varying functions of the form α(x) = a(εx) and β(x) = b(εx), where ε > 0 is a small parameter.Letting u tot denote either the total electric field E in E z polarization or the total magnetic field H in H z polarization, it follows from (4), ( 5) and ( 6) that the transition conditions can be equivalently expressed as in terms of the metasurface parameters α and β introduced in ( 6), where we have used the notation u y = ∂u/∂y.
In this paper we consider the problem of scattering that arise when the metasurface is illuminated by a time-harmonic incident field u inc which is assumed to satisfy the Helmholtz equation ∇ 2 u inc + k 2 u inc = 0 in Ω + and Ω − (Fig. 1).In order to properly formulate a scattering problem, we proceed to express the total field as u tot = u scat + u inc , where u scat denotes the scattered field off of Γ. Replacing u tot = u scat + u inc in the Helmholtz equation and the transition conditions (7) we obtain that u scat satisfies In order for (8) to be a well-posed boundary value problem for u scat , the scattered field has to satisfy a certain radiation condition at infinity [13,4,31].Such a radiation condition, which roughly speaking means that u scat corresponds to an up-going wave-field in Ω + and a down-going wave-field in Ω − , can be formally stated in terms of the angular spectral representation by requiring the existence of functions (or more generally, distributions) A + and A − such that where contour C corresponds to the real k x -axis that is suitably dented around the possible poles singularities of A ± [15].

Exact and approximate surface integral representations
In this section we derive an exact and an approximate integral representation formulae for the scattered field u scat solution of (8).Such formulae involve the incident surface currents and the Green's function of the boundary value problem (8) and are given in terms of integrals on the metasurface only,

Exact integral representation
The Green's function G of the boundary problem ( 8) can be physically interpreted as the total field produced by a point source excitation placed off of the metasurface [10,4].In detail, letting r = (x , y ) denote the location of a point source and r = (x, y) denote an observation point, the Green's function G(r|r ) can be found by solving the following boundary value problem: with δ denoting the Dirac's delta distribution and where G is additionally required to satisfy the radiation condition.
As is shown in Appendix A, Green's third identity together with ( 8) and ( 9) can be combined to show that the scattered field u scat admits the integral representation where the current densities f inc ± in (10) are given by f inc ± (s) = u inc y (s, 0) ± ikα(s)u inc (s, 0) (11) in terms of the incident field u inc , and where G(r|s, 0 ± ) = lim δ→0 + G(r|s, ±δ).
The integral representation formula (10) provides an explicit expression for the scattered field which is valid everywhere (in the near and far fields).It has, however, little practical relevance unless an exact or approximate Green's function is available.Unfortunately, a formula for the Green's function ( 9) cannot be easily obtained for "general" spatially varying metasurface parameters α and β, and thus, suitable approximations of G are needed in order to make proper use of (10) in scattering simulations.In the next section we derive an approximation for G based on ray-optics principles.

Ray-optics approximation
From the viewpoint of Huygens' principle (formalized by the principal of equivalence), equation (10) represents the scattered field by a source term corresponding to each point along the wavefront incident upon the surface (f inc ± ) [5].The typical "ray-optics" approximation is to compute the reflection/transmission at each point x as if the surface were uniform in the vicinity of that point.That approach corresponds to approximating (10) by a similar equation, but with the exact Green's function G replaced by an approximate "proto"-Green's function G p defined by the scattering of the source at r = (x , y ) from a uniform surface α(x ) and β(x ).This approximation yields which which turn out to be our zeroth order approximation in Sec. 4. We call G p a proto-Green's function because it is a building-block for our solution, but it is not the Green's function one would get by putting a point source as the incident field in (12).We give an exact Green's function G p for reflection and transmission off a uniform surface in Appendix B. But in the far field (fields far from the surface), as is derived rigorously in Appendix F, this simplifies to a function G ff p that we present in a more elementary fashion here.16) of the Green's function of the problem under consideration (9).The field produced by a point source at r impinges on the metasurface at r .The total field is observed at the point r far away from the metasurface that is denoted by Γ.The total field is decomposed in primary and secondary fields.The latter corresponds to the field reflected off of the metasurface, which is approximated by means of the local reflection coefficient R in (13) assuming specular reflection.
In order to construct G ff p , we consider the scattering configuration depicted in Fig. 2. With reference to that figure, the total wave field observed above the metasurface at a point r = (x, y), y > 0 is given by the superposition 0 (k|r − r |) produced by a point source placed above the metasurface, at r = (x , y ), y > 0, and the (secondary) field G r resulting from the reflection at r = (x , 0) (on Γ) of the ray stemming from r .(The function H (1) 0 is the Hankel function of first kind and order zero [1].)The magnitude and phase of the reflected field are characterized by the local reflection coefficient R(k cos θ, x ) that depends on the reflection angle ϑ = arctan(y/(x − x )) (which is measured with respect to tbe metasurface) of the reflected ray at r (Fig. 2).As is shown in Appendix B, the local reflection coefficient is given by where In order to account for both the magnitude and the direction of the reflected field, we consider the image point source r = (x , −y ) ∈ Ω − which allows the total field field above the metasurface to be approximated as Similarly, below the metasurface the total field at a point r ∈ Ω − corresponds to the transmitted field which can be approximated as in terms of the local transmission coefficient which is shown in Appendix B to be given by in terms of Γ α and Γ β defined in (14).For a point source placed below the metasurface at a point r in Ω − , the approximate Green's function G ff p can be derived from symmetry arguments.A rigorous derivation of G ff p based on asymptotic analysis is presented in Appendix F.
Note that G ff p provides a valid approximation of the exact Green's function G in the far-field zone, and can be thought of as the approximate field scattered from a single point on the surface.In view of the Huygens' principle, to get the total field we must add together all of the surface points resulting in what we refer to as the ray-optics approximation: of the total field u tot is achieved by replacing G by G ff p in the integral representation formula (10), where the functions f inc ± are defined in (11) and where the limits G ff p (r|s, 0 ± ) = lim δ→0 + G ff p (r|s, ±δ) are obtained by setting in formulae (15a) and (15b), respectively.We will show in Section 4 that ( 17) is a simplified version of ( 12) in the limit of sources and fields far from the surface and, furthermore, that ( 17) and ( 12) are the zero-th order terms in a convergent series of corrections for slowly varying surfaces.
In order to examine the accuracy of the ray-optics approximation of the total-field (17), we consider a series of numerical examples.Fig. 3 presents a comparison of the approximate and "exact" total field solution of the problem of scattering of a planewave that impinges at normal incidence upon three different mesaturfaces.The metasurface parameters α and β were selected so that the transmitted fields satisfy the so-called generalized laws of reflection and transmission [49,47,16].As is shown in Appendix D, for a given incident planewave u inc in the direction kinc = (cos θ inc , sin θ inc ), −π/2 ≤ θ inc ≤ 0, metasurface parameters of the form produce (to leading-order asymptotics) a transmitted field corresponding to a single planewave in the direction kt = (cos θ t , sin θ t ).
Although the ray-optics approximation (17) seems to capture qualitatively the main features of the scattered field, quantitatively it exhibits large near-field errors-and also large far-field errors in some cases (Fig. 3(i), for example)-that are just one order of magnitude smaller than the scattered field itself.In what follows of this paper we present a methodology to produce both near-and farfield corrections to the ray-optics approximation (17), which turns out to be just the zeroth-order terms of a series approximation of the exact scattered field in the far-field.

Locally uniform approximation and corrections to ray-optics
This section presents an SIE formulation of the problem of scattering (8) from which corrections to the ray-optics approximation (17) can be easily obtained in the form of a Born (or Neumann) series [5].We here follow a standard indirect integral equation formulation procedure [13,29] in which the field is represented by means of a Green's function G p that satisfies the Helmholtz equation in both upper and lower homogeneous media, but does not satisfy the correct sheet transition conditions at the metasurface, and we solve for effective source terms that restore the desired transition conditions.Note that, unfortunately, G ff p cannot be used for this purpose because it does not satisfy the Helmholtz equation, since both reflection and transmission coefficients depend on the location of the source and observation points.
Just as in Section 3.2, we begin by constructing a proto-Green's function G p , which is given in terms of Fourier-like integrals that we deem as Sommerfeld integrals (due to the similarities they share with layered-media Sommerfeld integrals [10]).This Green's function possesses two important features.On one hand G p -as G itself-satisfies the inhomogeneous Helmholtz equation (9a) with a point source excitation and, on the other hand, its far field equals the approximation G ff p of the exact Green's function G, used in the ray-optics approximation (17).The former allows us to properly derive a second-kind SIE, while the latter guarantees that the zeroth-order approximation obtained by truncation of the Born series solution of the second-kind integral equation does indeed correspond to the ray-optics approximation (17) in the far-field zone.As in Section 3.2, we use the term "proto" because G p is just a building block and is not equal to the Green's function of our final zeroth-order approximation.
The key feature of the proto-Green's function is that, instead of satisfying the non-local transition conditions (9b)-(9c), it satisfies the following local transition conditions G p,y (r|r ) = −ikα(x ) G p (r|r ) and G p,y (r|r ) = −ikβ(x ) G p (r|r ) , r ∈ Γ, (20) with metasurface parameters α and β depending on the source point r = (x , y ).This local transition conditions formalize the slowly varying assumption used in Section 3.2 where the rayoptics approximation was derived.Since the metasurface parameters α and β that appear in local transition conditions (20) do not depend on the observation point r = (x, y), they can be treated as constants and, thus, an analytical expression for G p in terms of Sommerfeld integrals can be easily obtained.The idea behind this calculation is to decompose the point-source incident field G inc as a superposition of both propagative and evanescent planewaves.Since specular reflection takes place for each plane wave impinging on the metasurface, the resulting scattered field can be written down as a superposition of reflected and transmitted plane waves weighted by the reflection and transmission coefficients provided in ( 13) and ( 16), respectively.The details of this derivation are presented in Appendix B. Furthermore, it is shown in Appendix F by means of a detailed asymptotic analysis that, to leading asymptotics, G p equals G ff p as |r| → ∞.With an analytical expression for G p in hand (i.e., formulae (42) and ( 44)) we proceed to derive an SIE for the solution of the scattering problem (8) from which corrections to the rayoptics approximation (17) can be computed.Following the exact integral representation (10) of the scattered field u scat we introduce an indirect integral formulation for the scattering problem (8) by setting where ϕ and ψ are (so far) unknown surface density functions.Note that if G p were the exact Green's function G, then u scat in ( 21) would be the exact solution of (8) provided ϕ = f inc + and ψ = f inc − , where f inc + and f inc − are defined in (11).Note further that in virtue of the relationship between G p and G ff 0 established in Appendix F, the substitutions ϕ = f inc + and ψ = f inc − in (21) would produce an approximation of u scat that, in the far-field zone, exhibits the same accuracy of the ray-optics approximation (17).
Continuing with the derivation of the SIE, we observe that in order for u scat (21) to be an exact solution of ( 8) it has to satisfy both the Helmholtz equation (8a) and the transition conditions (8b)-(8c).The problem here is that, although u scat (21) does satisfy the Helmholtz equation (8a) in Ω + and Ω − for any admissible densities ϕ and ψ (since be construction G p (r|r ), r ∈ Γ, satisfies it), it does not necessarily fulfill the correct transition conditions (9b)-(9c) unless (ϕ, ψ) is solution of a certain SIE.Indeed, it is shown in Appendix C that, imposing the transition conditions (9b)-(9c) on u scat in ( 21), an SIE for the unknown density functions (ϕ, ψ) is obtained.The resulting equations correspond to two decoupled second-kind SIEs: for the unknown auxiliary densities µ j , j = 1, 2, which are directly related to the integral densities in ( 21) by The precise definition of the integral operators T j , j = 1, 2, is given in Appendix C and the functions g inc j , j = 1, 2, on right-hand-side of ( 22), are Clearly, the densities µ j , j = 1, 2, can be determined by solving the SIEs (22) and from them, the desired densities ϕ and ψ that make u scat in (21) the exact solution of (8) can be readily obtained.
We now recall that we are here interested in slowly varying interface parameters α and β of the form α(x) = a(εx) and β(x) = b(εx) where ε > 0 is a small parameter.In view of definitions (50) and (51), we observe that both integral operators, T 1 and T 2 , vanish as ε → 0 and, therefore, in the limit when ε = 0 the exact SIE solutions are simply µ j = g inc j , j = 1, 2. For small but nonzero values of ε, in turn, convergent Neumann-series solutions of the SIEs ( 22) can be obtained because the integral operators satisfy T j < 1 in a certain operator norm for sufficiently small ε.
The N th-order approximations of the density functions ϕ and ψ can thus be defined as where µ (N ) j , j = 1, 2, are the truncated Neumann series From (26) we then define the N th-order locally uniform approximation of the total near and far fields: for r ∈ R 2 \ Γ, respectively.Note the zeroth-order term (or any higher order term) in (28) does not correspond to setting ε = 0 (a uniform surface).Finally, it follows from the definitions above that the ray-optics approximation ( 17) is simply the zeroth-order approximate far-field u tot,ff 0 , i.e., the N = 0 instance of the formula (28b).In order to see this it suffices to note that ϕ 0 = f inc + and ψ 0 = f inc − , which followed directly the definition of f inc ± in (11) and the fact that In summary, the N th-order approximation of the total near and far fields resulting from the scattering of an incident field u inc off of a metasurface Γ, can obtained as follows: 1. Evaluate the input data g inc j , j = 1, 2, defined in (24), using the prescribed incident field u inc .
2. Compute the N th-order approximate densities µ (N ) j , j = 1, 2, defined in (27), by repeated application of the integral operators T j to g inc j .
In order to illustrate the accuracy yielded by the higher-order corrections to the ray-optics (zeroth-order) approximation, we present Fig. 4 which concerns the scattering configuration considered above in Figs.3(g)-(i), which corresponds to the scattering of a planewave that impinges at normal incidence on a metasurface that renders a transmitted planewave with wavevector forming an angle of 22.5 • with respect to the metasurface (Fig. 4(a)).Figs.4(b) and 4(c) display the real part of the total fields (incident + reflected, and transmitted) produced by the zeroth and first order approximations of the far field, corresponding to formula (28b) with N = 0, 1.In order to better visualize the convergence of the locally uniform approximations (28b) as N increases, we present Figs.4(d)-(h) that display the absolute value of the zeroth, first, second, third and fourth order far-field errors.The reference "exact" far-field was computed by direct solution of the SIE system (22).These results indicate that the far-field error is roughly reduced by a factor of 0.5 as the order increases.This is explained by the fact that the spectral radii of the discrete versions of the integral operators T j , j = 1, 2, are approximately 0.5.More details on the convergence of the Neumann series approximation are given in Appendix E.
Our next example concerns the dependence of the rate of convergence of the N th-order approximations (28a) and (28b) on the smoothness of the metasurface parameters.As it turns out,  for constant metasurface parameters the zeroth-order near and far field approximations are exact.In this example we thus attempt to quantify how errors depart from zero as the metasurface parameters become non-constant.In order to so we consider slowly-varying metasurface parameters α(x) = a(εx) and β = b(εx)-which depend on a small parameter ε > 0-that tend to constants α 0 = a(0) and β 0 = b(0) as ε → 0. Fig. 5 displays the near-field errors for vanishing values of the smoothness parameter ε > 0. Clearly, the zeroth-, first and third-order approximations exhibit errors of order O(ε 2 ), O(ε 4 ) and O(ε 6 ), respectively, as ε → 0, i.e., as the metasurface parameters tend to constants.
Interestingly, the detailed asymptotic calculations presented in Appendix F also reveal that surface-wave modes appear in the asymptotic expansion of G p for certain constant values of the metasurface parameters (note that for constant α and β, it holds that G = G p ).Such surfacewave modes are also present in the field scattered by metasurfaces with non-constant metasurface parameters.To demonstrate this fact, we present Fig. 6 which displays the total field solution of the problem of scattering of a Gaussian beam by a metasurface for which a surface-wave mode propagates from left to right along.Three difference solution are displayed in that figure: the exact solution, the ray optics approximation (17), and the zeroth-order locally uniform approximation (28a) with N = 0. Since G ff p is a far-field approximation (which is valid at a certain distance from the metasurface), u tot,ff 0 does not capture at all the aforementioned surface-wave modes.

Concluding remarks
We developed an SIE approach, based on a locally uniform approximation of a metasurface, to establish the accuracy and compute higher-order corrections to a ray-optics approximation commonly used in inverse metasurace design, where metasurfaces are modeled by means of slowly varying surface parameters.
This work opens many research directions that could be pursued in the future.The most important (and straightforward, in principle) is perhaps the extension of the proposed approach to three-spatial dimensions.As a practical matter, however, there are many subtle implementation aspects of this extension, such as the derivation of suitable three-dimensional SIE formulations and the efficient evaluation of the resulting two-dimensional surface integrals, that need to be addressed.
Another future research direction is the extension of the proposed approach to more general classes of metasurfaces that cannot be modeled by means of sheet transition conditions and in particular, to approximate metasurfaces as locally periodic rather than locally uniform.This extension, however, poses new theoretical challenges.Such an extension requires the knowledge of a certain proto-Green's function associated with a periodic transmission problem, which does not admit an expression in terms of Sommerfeld integrals and must be computed numerically.Despite these theoretical challenges, there is both numerical and experimental evidence that a locally periodic approximation is sufficiently accurate for practical metasurface design [37,2,3,48,49,47,36].
Finally, we mention that there remains considerable room for further asymptotic analysis of the integral operators T j , j = 1, 2, and their convergence as ε → 0, to rigorously establish the convergence rate of the Neumann-series solution (25) (i.e. the corrections to ray optics).Although it is clear that T j → 0 as ε → 0, an intricate analysis is required to obtain convergence rates, and to clearly specify for which function spaces convergence is obtained, especially for unbounded surfaces and incident fields where limiting processes are tricky to apply to the surface integrals (51).

B Sommerfeld-integral Green's function approximation
This appendix in devoted to the derivation of a Sommerfeld-integral [10] representation of the proto-Green's function G p used in the SIE derivations presented in Section 4 above.
As was mentioned above in Section 4, G p satisfies both the Helmholtz equation (9a) and the radiation condition, but instead of the transition conditions (9b)-(9c), it satisfies the locally uniform transition conditions (20).In order to find an expression for G p we first note that since both metasurface parameters α and β are taken to be functions of the source point r = (x , y ), they are constant as functions of r = (x, y) and thus it is possible to find an exact expression for G p by standard Fourier transform techniques.In fact, for a source point r ∈ Ω + the proto-Green's function G p can be interpreted as the total field produced by the incident field The square root x is defined in the complex k x plane as the product where the first square root has a branch cut along the positive imaginary axis, and the second one has a branch cut along the negative imaginary axis.Fig. 8 depicts the domain of definition of k y along with the curves in the complex plane where the real and imaginary parts of k y change sign.
Using G inc as incident field, the total field-which corresponds to G p -can be expressed as where the reflected and transmitted fields admit the integral representations in terms of the reflection and transmission coefficient R and T defined in ( 13) and ( 16), respectively.The special integral sign "c " introduced in (42) refers to the fact that the path of integration passes below (resp.above) any pole that the integrands may have on the positive (resp.negative) real k x -axis.For the sake of definiteness in what follows of this paper the integral sign c refers to a contour integral along the path SC depicted in Fig. 9.
In order to establish the validity of (42), we note that for G p to satisfy the locally uniform transition conditions (20), R and T have to related by the equations where for notational simplicity we have let R = R(k x , x ), T = T (k x , x ), α = α(x ) and β = β(x ).Solving for R and T from ( 43) we obtain the expressions in ( 13) and ( 16) utilized in the previous section.
Similarly, it follows from the symmetric of G p established in Appendix A that for a point source r ∈ Ω − the total field G p takes the form

C Integral-equation formulation for corrections
This appendix is devoted to the derivations of the SIEs (22).To simplify the notation, we first define the functions g e i (s, σ) = G p (s, 0 e |σ, 0 i ) and ∂ y g e i (s, σ) = G p,y (s, 0 e |σ, 0 i ) for s, σ ∈ R, where the indices i and e correspond to the symbols "+" or "-" that refer to the limit values (from above and below, respectively) on Γ.
From the Sommerfeld-integral representation of G p in ( 42) and ( 44) it thus follows that where with Γ α and Γ β being defined in (14).
Next we introduce the boundary integral operators is the total field resulting from the scattering of the planewave u inc off of a metasurface with constant interface parameters α 0 and β 0 .The field v = u − u 0 , on the other hand, satisfies the Helmholtz equation in R 2 \ Γ, the radiation condition, and the transition conditions Letting then G 0 = G p denote the Green's function in (42) and ( 44) corresponding to constant interface parameters α 0 and β 0 , we obtain from the discussion in Section 4 that the zeroth-order approximation (i.e., (28a) with N = 0) of v is given by where letting the approximate densities ϕ 0 and ψ 0 are given by Replacing G 0 in (53) by its far-field approximation, derived in Appendix E, we find that where, after some algebraic manipulations, the far-field pattern v ∞ (θ|θ inc ) can be expressed as where From (56) it thus follows that the far-field pattern v ∞ in (55) would correspond to a linear combination of planewaves with wavevectors k α = k(cos θ α , sin θ α ) and ∞ were Dirac delta distributions supported at angles θ α and θ β , respectively.Formally, this can be achieved by selecting α(s) = c α e ikdαs and β(s) = c β e ikd β s , with c α , c β ∈ C and d α and d β being such that Finally, the expressions in (19) are obtained by letting α 0 = β 0 = − sin θ inc , c α = ±α 0 , and c β = ∓β 0 , for which a minimal reflection off of the metasurface is achieved.

E On the convergence of the Neumann series (25)
In our next appendix we consider an example to study the convergence of the Neumann series approximation (27) by examining the dependence of the spectral radii of the discretized integral operators T j , j = 1, 2, on the smoothness of the metasurface parameters α and β.We consider here the discretized version of the SIEs (22) which take the form (I − T j )µ j = g j , where I is the identity matrix, T j is the discretized integral operator by the method described in Appendix G, µ j is the unknown vector, and g j is the discretized interface data.It is easy to show that if T j < 1 in some matrix norm, then the relative error in the discretized Neumann series approximation can be bounded by for some constant c > 0 where ρ(T j ) = max n |λ n (T j )| denotes the spectral radius of T j .This bound shows that the spectral radius provides an approximate rate of convergence of the Neumann series approximation as N increases.Moreover, it can be shown that the Neumann series for the discrete linear system converges if and only if the ρ(T j ) < 1.In the following example then, we consider a metasurface parameter α(x) = a(εx) given by the truncated Fourier series where the coefficients c are randomly generated from a uniform distribution and are also adjusted so that the constrain Re α ≥ 0 is satisfied.Clearly, ε > 0 is a parameter that controls the smoothness of α.Fig. 7 displays the spectral radius of the matrix T 1 in log-log scale for a range of values of ε.
There results demonstrate that convergence of the Neumann series is in fact expected for a large range of values of ε, including some of those that give rise to quite rough interface parameters α.
F The far field of G p and its relationship with G ff p This Appendix presents a detailed asymptotic analysis that establishes rigorously the relationship between G p and G ff p , as well as the existence of guided modes.To establish the relationship between G ff p and G p , we derive the far-field asymptotic approximation (as |r| → ∞) of the proto-Green's function G p given in (42) and (44).In order to do so we resort to the method steepest descents for which we follow the analysis of Sommerfeld integrals presented [4,Chapter 8].Similar saddle point calculations can also be found in classical references on layered media scattering, such as [6,10].
Assuming first that r ∈ Ω + and letting r = |r|(cos θ, sin θ), θ ∈ (−π, 0) ∪ (0, π), we have that the resulting reflected and transmitted fields in (42) can be expressed as G r (r|r ) = c SC q r k x , r e |r|φ(kx) dk x and G t (r|r ) = c SC q t k x , r e |r|φ(kx) dk x (59) in terms of the phase and amplitude functions defined as respectively.Note that in this case (r ∈ Ω + ) we are interested in G r for θ ∈ (0, π) and in G t for θ ∈ (−π, 0).Three kinds of critical points have to be taken in account in the steepest descent method approximation of the integrals (59), namely, saddle points of the phase function φ, (possible) poles singularities of the integrands q r and q t , and the branch points of the square root x .We first consider the saddle points of φ, which correspond to solutions k * x ∈ C of the algebraic equation φ (k * x ) = 0.In view of φ (k x ) = i cos θ − ik x | sin θ|/k y follows that there is only one saddle point on SC given by k * x = k cos θ = kx/|r| at which φ(k * x ) = ik.The steepest descent directions from k * x , on the other hand, are given by the angles 3π/4 and −π/4, which were obtained from φ (k * x ) = −i/(k sin 2 θ) = 0. We then conclude that the steepest descent path SD is given implicitly by the equation Im φ(k x ) = k from which it can be shown that SD intersects SC again at k x = k sec θ and that and Im Fig. 9 depicts the steepest descent paths for x > 0 and x < 0. SD coincides with SC when x = 0.
Re k x With this information in hand we then proceed to deform the Sommerfeld contour SC to the steepest descent contour SD that passes through the saddle point k * x .Note that SD does not intersect the branch cuts stemming from ±k and, thus, there is no contribution to the asymptotic expansions from the branch points at ±k.There might be, however, contributions arising from pole singularities of the integrands.In fact, from the expressions for the reflection and transmission coefficients in ( 13) and ( 16), respectively, we have that the poles of both q r (61) and q t (62)-which correspond to the poles of the functions Γ α (x, x ) and Γ β (x, x ) defined in ( 14)-are solutions of the (independent) algebraic equations To find necessary and sufficient conditions on the metasurface parameters α and β for such poles to exist, we first note that the conditions Re α ≥ 0, Re β ≥ 0 and the equations (64) imply that poles of q r or q t could only exist in the regions of the complex k x plane where Re k y (k x ) ≤ 0. It is easy to see that such regions amount to D (Fig. 8).Furthermore, since Im k y = k 2 − k 2 x ≥ 0 in D, poles will exist if and only Im α ≤ 0 or Im β ≤ 0.
By the Cauchy residue theorem we thus have that the poles of q r (resp.q t ), if any, will only contribute to the far-field expansion of G r (resp.G t ) if they lie within the region in the complex plane enclosed by SC and SD.In order to determine whether a pole of q r (resp.q t ) lies inside that region, and in view of the fact that we do not have access to an explicit parametrization of SD, we utilize the asymptotic identities (63).Doing so we conclude that the relevant poles-that are henceforth denoted by k (p) x , p = 1, 2-must meet the conditions k (1)  x = k 1 − α 2 or k (2)  x = k 1 − β 2 , and 0 ≤ Im k (p) x ≤ Re k (p) x | cot θ| − k| csc θ| if x > 0, (65a) and k (1)  x A r p e ik (p) as |r| → ∞, where r is the image point source r = (x , −y ) (Fig. 2).The amplitudes A r 1 and A t 1 of the guided waves in (66) are directly obtained from the residues of q t and q r and read as x in (65) exists, and they equal zero otherwise.Similarly, x in (65) exists, and they equal zero otherwise.The contribution to the asymptotic expansions of the pole singularities corresponds to surface waves that travel away from the point source r = (x , y ) and are confined to a narrow strip containing the metasurface-as they decay exponentially fast toward the upper and lower halfplanes.For example, the contribution of the pole k to the asymptotic expansion of the reflected and transmitted fields equals with + and − corresponding to the reflected and transmitted fields, respectively.In order to establish the relationship between G ff p and G p we recall the asymptotic identity  15) with ( 42), we conclude that Therefore, G ff p (15) can be simply interpreted as an approximation of the proto-Green's function G p for point sources at the metasurface and observations points far away from the metasurface.Note that for such configuration of source and observation points, the surface wave modes do not have any significant contribution as they decay exponentially as |y| → ∞.
Finally, in order to demonstrate the validity of the asymptotic expansions derived in this section we present Fig. 10 that displays the real part of G p (r|r ) for constant interface parameters α = −i and β = −1.2i.Note that for these interface parameters the terms Γ α (k x , s) and Γ β (k x , s) defined in (14) have poles on the real axis which make surface waves mode of the form (69) appear in the asymptotic expansions (66).

G Numerics
In this appendix we briefly describe a high-order method for the numerical evaluation of the Sommerfeld integrals G r and G t in ( 42)-( 44) and the integral kernels K 1 and K 2 in (50).This approach, which was originally developed for layered-media scattering problems [35] (see also [34,Section 2.3.5]), is a combination of the contour-integration method described in [33] and the the smooth-windowing approach put forth in [28] for the evaluation of oscillatory integrals.
that passes below all the singularities of the integrand f .The contour C 2 , on the other hand, is simply the interval [L, ∞] on the real axis.Note that on C 1 the function f (ζ(t)) grows exponentially as t increases from 0 to π/2.Indeed, it can be shown [34] that e id  .This simple procedure ensures that the exponential terms of f remain bounded by one along C 1 .The resulting expression for the contour integral I 1 is then approximated by means of the Clenshaw-Curtis quadrature rule [14]-which, for the smooth integrand under consideration, yields rapid convergence.In view of the oscillatory behavior of the integrand and in order to maintain the same accuracy for all d 1 and d 2 , the number of quadrature points is chosen to grow linearly with d 1 .
In order to evaluate the oscillatory integral I 2 , on the other hand, we utilize the windowing method put forth in [28].Using this procedure I 2 is approximated as where the window function w A is defined as w A (t) := η(t; cA, A), A > 0, 0 < c < 1, in terms of the C ∞ (R) function η(t; t 0 , t 1 ) = which equals one on [−t 0 , t 0 ] and is supported on the (bounded) interval [−t 1 , t 1 ].
In virtue of the oscillatory behavior of the integrand when d 1 = 0, and the exponential decay of the integrand when d 2 = 0, the integral on the right-hand-side of (75) converges to I 2 faster than any negative power d 2 1 + d 2 2 A as A goes to infinity-as proved in [34,Proposition 2.3.4].In the special case d 1 = d 2 = 0, however, the f is slowly decaying on C 2 and does not oscillate, thus, it leads to slow (algebraic) convergence of the windowed-integral approximation (75) to the integral I 2 as A → ∞.In fact, in that case the error in the windowed-integral approximation decays as O((cA) −1 ) [34,Proposition 2.3.4].Therefore, as in the case of I 1 , the integral on the right-hand-side of (75) is here approximated by using Clenshaw-Curtis quadrature.The superalgebraic/exponential convergence of the windowed integral allows I 2 to be approximated with a fix accuracy and a fixed computational cost by choosing A inversely proportional to d 2 1 + d 2 2 .
Integral kernels In order to numerically evaluate the integral kernels K 1 and K 2 defined in (50), we resort to the contour-integration procedure described above.The evaluation of these kernels

Figure 1 :
Figure1: Scattering problem under consideration.A time-harmonic incident field u inc with angular frequency ω > 0 impinges on a metasurface with surface impedance Z and surface admittance Y , producing a reflected field propagating in the upper half-plane and a transmitted field propagating in the lower half-plane.The metasurface is surrounded by a homogeneous medium with electric permittivity 0 > 0 and magnetic permeability µ 0 > 0.

Figure 2 :
Figure 2: Variables used in the derivation of the approximation (16) of the Green's function of the problem under consideration(9).The field produced by a point source at r impinges on the metasurface at r .The total field is observed at the point r far away from the metasurface that is denoted by Γ.The total field is decomposed in primary and secondary fields.The latter corresponds to the field reflected off of the metasurface, which is approximated by means of the local reflection coefficient R in (13) assuming specular reflection.

Figure 3 :
Figure 3: First column: Ray-optics approximation (17) of the total field u tot = u inc + u scat solution of the problem of scattering of a planewave that impinges at normal incidence on metasurface's that produce transmitted fields that are predominantly planewaves propagating in the directions θ t (a) 45 • , (d) 30 • and (g) 22.5 • with respect to the metasurface.Second column: Full "exact" total field solution of the corresponding for planewave transmitted fields in the directions θ t (b) 45 • , (e) 30 • and (h) 22.2 • .Third column: Absolute errors for planewave transmitted fields in the directions θ t (c) 45 • , (f) 30 • and (i) 22.5 • .
error real part of the total field

Figure 4 :.
Figure 4: High-order corrections to the ray optics approximation (17).(a) Geometrical configuration of the problem under consideration which corresponds to an unit amplitude planewave impinging at normal incidence on a metasurface that renders a transmitted planewave with wavevector forming an angle of 22.5 • with respect to the metasurface.(b) and (c): Real part of the total field 0th and 1st order approximations.(d), (e), (f), (g) and (h): Absolute errors |u − u tot,ff N |, for N = 0, 1, 2, 3 and 4, respectively, in the zeroth (no correction), first, second, third and fourth order corrections to the ray optics approximation u tot,ff 0 .The color scales were adjusted according to the maximum error displayed in each one of the figures.

Figure 5 :
Figure 5: Rate of convergence of the zeroth-, first-, second-, third-and fourth-order near-field approximations (28a) in terms of the smoothness of the metasurface parameters.The plot is in log-log scale.The problem under consideration is the scattering of a unit amplitude plane wave at normal incidence that impinges on a metasurface with metasurface parameters α(x) = a(εx) = c 0 {1 − e iεx } and β(x) = b(εx) = c 0 {1 − e −iεx }.The parameter ε > 0 controls the smoothness of metasurface parameters.The color curves display the maximum of the absolute error |u tot (r) − u tot N (r)| in the near-field approximations (28a) evaluated at the spatial points r ∈ {−1, 0, 1} × {−10, 10}.

Figure 6 :
Figure 6: Total field solution of the problem of scattering of a Gaussian beam by a metasurface with interface parameters α and β engineered to allow for surface wave modes (69) to propagate along the metasurface.

Figure 7 :
Figure 7: Spectral radius of the discretized integral operator T 1 as a function of the smoothness of the interface parameter α in (58).The yellow painted strip indicates the range of parameters ε for which the Neumann series diverges.

Figure 8 :
Figure 8: Signs of the real and imaginary parts of square root k y (k x ) = k 2x − k 2 in the complex k x plane.

Figure 9 :
Figure 9: Sommerfeld contour SC (continuous blue line) and steepest descent contour SD (dashed red line) utilized in the evaluation of the integral (59).Poles of the Γ α and Γ β lying in the shaded regions contribute to the far-field asymptotic expansion of the approximate Green's function.

Figure 10 :
Figure 10: Real part of the total field solution of the problem of scattering of the field produced by a point source by a metasurface with constant interface parameters α = −i and β = −1.2i.This selection of the interface parameters allow a surface wave mode of the form (69) to propagate along the metasurface.The point source is located at distance 0.25λ above the metasurface.(a) G ff p which in this case corresponds to the ray-optics approximation of the exact Green's function G. (b) Rigorous far-field asymptotic approximation (66) of G p .(c) Direct evaluation of the proto-Green's function G p defined in (42) and (44).

Figure 11 :
Figure 11: Integration contours used in the evaluation of the proto-Green's function G p , and the integral kernels K 1 and K 2 .