1. Introduction

Optical metasurfaces, subwavelength structures described by an effective sheet impedance [1–10], 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 [11–14]. This is a good approximation for surfaces (or unit cells) that are varying slowly, a fact that is closely connected to the “adiabatic theorem” [15,16] for waves propagating through slowly varying media. However, although there are countless papers and books on modeling propagation through slowly varying media [15,17,18], 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 [1–10] in subsuming the microscale into effective surface impedances.

 

Fig. 1 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 and magnetic permeability μ₀ > 0.

Download Full Size | PDF

In particular, we use the technical machinery of surface-integral equations (SIEs) [19,20] 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 ε². 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 [6,21] 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 [22], 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 [23]) 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 [10,24–30]. (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 [10,31–33].) 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 [34]. Yet another closely related approximation is provided by scalar diffraction theory [35].) 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” [36] eigenfunctions of each cross-section [17] or period [16], 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 intermodal scattering [17]. 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 [18, 37] (also known as a volume-current method [38], Kirchhoff approximation [34], 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 ε [19,20]. 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 [39], the corrections are expressed in terms of integrals involving G_p. These integrals must be computed numerically by a “quadrature” technique [40] (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 a certain extent the surface variation. In the far field, G_p simplifies to an expression $G_{\text{p}}^{\text{ff}}$ 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 ε², and more generally as ε^2N+2 if we include Nth-order corrections (Fig. 5). In the presence of guided modes, which correspond to poles that appear in G_p at certain wavevectors [37], 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).

2. 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 and μ₀ > 0, respectively (Fig. 1). The metasurface is characterized by the so-called generalized sheet transition conditions [1]:

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₁, H₂, 0), and can be obtained from the z-component of the electric field by means of the relations:

Relations in Eqs. (2) and (3), on the other hand, yield that the transition conditions in Eq. (1) can be equivalently expressed as

In order to treat both E_z⁻ and H_z-polarization cases, we define the metasurface parameters α and β as

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 Eqs. (4), (5) and (6) that the transition conditions can be equivalently expressed as

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 ∇²u^inc + k²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 in Eq. (7) we obtain that u^scat satisfies

3. 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 the problem in Eq. (8). Such formulae involve the incident surface currents and the Green’s function of the boundary value problem in Eq. (8) and are given in terms of integrals on the metasurface only,

3.1. Exact integral representation

The Green’s function G of the boundary problem in Eq. (8) can be physically interpreted as the total field produced by a point source excitation placed off of the metasurface [37, 41]. 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:

As is shown in Appendix A, Green’s third identity together with Eqs. (8) and (9) can be combined to show that the scattered field u^scat admits the integral representation

The integral representation formula in Eq. (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 in Eq. (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 Eq. (10) in scattering simulations. In the next section we derive an approximation for G based on ray-optics principles.

3.2. Ray-optics approximation

From the viewpoint of Huygens’ principle (formalized by the principal of equivalence), equation in Eq. (10) represents the scattered field by a source term corresponding to each point along the wavefront incident upon the surface ($f_{\pm }^{\text{inc}}$) [39]. 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 the integral expression Eq. (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

In order to construct $G_{\text{p}}^{\text{ff}}$, 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 $G_{\text{p}}^{\text{ff}}=G^{\text{inc}}+{\tilde{G}}^r$ of the (primary) incident field $G^{\text{inc}}(\mathit{r}|{\mathit{r}}^{\prime })=\frac{i}{4}{H}_0^{(1)}\left(k\left|\mathit{r}-{\mathit{r}}^{\prime }\right|\right)$ 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_0^{(1)}$ is the Hankel function of first kind and order zero [44].) 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

 

Fig. 2 Variables used in the derivation of the approximation in Eq. (16) of the Green’s function of the problem under consideration defined in Eq. (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 Eq. (13) assuming specular reflection.

Download Full Size | PDF

Similarly, below the metasurface the total field at a point r ∈ Ω₋ corresponds to the transmitted field which can be approximated as

For a point source placed below the metasurface at a point r′ in Ω₋, the approximate Green’s function $G_{\text{p}}^{\text{ff}}$ can be derived from symmetry arguments. A rigorous derivation of $G_{\text{p}}^{\text{ff}}$ based on asymptotic analysis is presented in Appendix F.

Note that $G_{\text{p}}^{\text{ff}}$ 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:

We will show in Section 4 that the expression in Eq. (17) is a simplified version of the one in Eq. (12) in the limit of sources and fields far from the surface and, furthermore, that the expression in Eq. (17) is 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. Figure 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 [8,28,29]. As is shown in Appendix D, for a given incident planewave u^inc in the direction k̂^inc = (cos θ^inc, sin θ^inc), −π/2 ≤ θ^inc ≤ 0, metasurface parameters of the form

 

Fig. 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°.

Download Full Size | PDF

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 far-field 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.

4. 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 [39]. We here follow a standard indirect integral equation formulation procedure [19,20] 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_{\text{p}}^{\text{ff}}$ 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 [37]). 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_{\text{p}}^{\text{ff}}$ 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

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 ray-optics 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

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:

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₁ and T₂, vanish as ε → 0 and, therefore, in the limit when ε = 0 the exact SIE solutions are simply ${\mu }_j=g_j^{\text{inc}}$, j = 1, 2. For small but nonzero values of ε, in turn, convergent Neumann-series solutions

The Nth-order approximations of the density functions φ and ψ can thus be defined as

Finally, it follows from the definitions above that the ray-optics approximation (17) is simply the zeroth-order approximate far-field $u_0^{\text{tot},\text{ff}}$, i.e., the N = 0 instance of the formula (28b). In order to see this it suffices to note that ${\phi }_0=f_{+}^{\text{inc}}$ and ${\psi }_0=f_{-}^{\text{inc}}$, which followed directly the definition of $f_{\pm }^{\text{inc}}$ in (11) and the fact that

In summary, the Nth-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:

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.

 

Fig. 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 $\left|u-u_N^{\text{tot},\text{ff}}\right|$, 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_0^{\text{tot},\text{ff}}$. The color scales were adjusted according to the maximum error displayed in each one of the figures.

Download Full Size | PDF

Our next example concerns the dependence of the rate of convergence of the Nth-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 real part of the total field parameters α(x) = a(εx) and β = b(εx)—which depend on a small parameter ε > 0—that tend to constants α₀ = a(0) and β₀ = b(0) as ε → 0. Figure 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(ε²), O(ε⁴) and O(ε⁶), 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 surface-wave 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_{\text{p}}^{\text{ff}}$ is a far-field approximation (which is valid at a certain distance from the metasurface), $u_0^{\text{tot},\text{ff}}$ does not capture at all the aforementioned surface-wave modes.

 

Fig. 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₀{1 − e^iεx} and β(x) = b(εx) = c₀{1 − e^−iεx}. The parameter ε > 0 controls the smoothness of metasurface parameters. The color curves display the maximum of the absolute error $\left|u^{\text{tot}}(\mathit{r})-u_N^{\text{tot}}(\mathit{r})\right|$ in the near-field approximations (28a) evaluated at the spatial points r ∈ {−1, 0, 1} × {−10, 10}.

Download Full Size | PDF

 

Fig. 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.

Download Full Size | PDF

 

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

Download Full Size | PDF

5. 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 [10, 25–30].

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).

Appendices

A. Exact integral representation

This appendix is devoted to the derivation of the integral representation formula (10). In order to achieve that, we show first the symmetry of the exact Green’s function (9). Consider then the functions w(r) = G(r|r₁) and v(r) = G(r|r₂) for the r₁, r₂ ∈ Ω₊. By Green’s formula and the radiation condition we have

(30)$$v({\mathit{r}}_1)-w({\mathit{r}}_2)={\int }_{{\mathrm{\Omega }}_{+}}\left\{w{\nabla }^{2}v-v{\nabla }^{2}w\right\}\text{d}\mathit{r}=-\underset{-\infty }{\overset{\infty }{\int }}\left\{w^{+}{v}_y^{+}-v^{+}{w}_y^{+}\right\}\text{d}s,$$

and, similarly

(31)$$0={\int }_{{\mathrm{\Omega }}_{-}}\left\{w{\nabla }^{2}v-v{\nabla }^{2}w\right\}\hspace{0.17em}\hspace{0.17em}\text{d}\mathit{r}=\underset{-\infty }{\overset{\infty }{\int }}\left\{w^{-}{v}_y^{-}-v^{-}{w}_y^{-}\right\}\hspace{0.17em}\hspace{0.17em}\text{d}s.$$

From the transition conditions (7), on the other hand, it follows that

(32)$$w^{+}{v}_y^{+}-v^{+}{w}_y^{+}=w^{-}{v}_y^{-}-v^{-}{w}_y^{-}.$$

where v^±(s) = v(s, 0^±) = lim_δ→0^± v(s, δ) and $v_y^{\pm }(s)=v_y(s,0^{\pm })={\text{lim}}_{\delta \to 0^{\pm }}$ v_y(s, δ) and similarly for w. Combining (30), (31) and (32) it is obtained that v(r₁)− w(r₂) = G(r₁|r₂)−G(r₂|r₁) = 0, and thus G(r|r′) = G(r′|r) for all r, r′ ∈ Ω₊. The identity G(r|r′) = G(r′|r) for all r, r′ ∈ Ω₋ can be derived in a completely analogous way. Consider now the functions w(r) = G(r|r₁) and v(r) = G(r|r₂) but with r₁ ∈ Ω₊ and r₂ ∈ Ω₋. Integration by parts yields the identity

(33)$$v({\mathit{r}}_1)=-\underset{-\infty }{\overset{\infty }{\int }}\left\{w^{+}{v}_y^{+}-v^{+}{w}_y^{+}\right\}\hspace{0.17em}\hspace{0.17em}\text{d}s=w({\mathit{r}}_2)$$

in this case, which clearly implies that G(r₁|r₂) = G(r₂|r₁). The same result can be easily obtained in the case r₁ ∈ Ω₋ and r₂ ∈ Ω₊.

Finally, the identities

(34)$$\left[\left[G_{y^{\prime }}(\mathit{r}|{\mathit{r}}^{\prime })\right]\right]=-k\alpha (x^{\prime })\left\{\left\{G(\mathit{r}|{\mathit{r}}^{\prime })\right\}\right\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\left\{G_{y^{\prime }}(\mathit{r}|{\mathit{r}}^{\prime })\right\}\right\}=-k\beta (x^{\prime })\left[\left[G(\mathit{r}|{\mathit{r}}^{\prime })\right]\right]$$

involving the normal derivatives of the Green’s function on the metasurface Γ for r ∈ Ω₊ ∪ Ω₋, follow straightforwardly from the reciprocity condition G(r|r′) = G(r′|r) established above.

Having established the symmetry of the Green’s function, i.e., G(r|r′) = G(r′|r) for all r, r′ ∈ Ω₋ ∪ Ω₊, we can now use it to prove the identity (10). Indeed, it follows from (8), the symmetry of G, and Green’s theorem that

(35)$$\underset{-\infty }{\overset{\infty }{\int }}\left\{G_{y^{\prime }}(\mathit{r}|s,0^{+})u^{\text{scat},+}(s)-G(\mathit{r}|s,0^{+})u_y^{\text{scat},+}(s)\right\}\text{d}s=\{\begin{array}{ll}{u}^{\text{scat}}(\mathit{r}),\hfill & \mathit{r}\in {\mathrm{\Omega }}_{+},\hfill \\ 0,\hfill & \mathit{r}\in {\mathrm{\Omega }}_{-},\hfill \end{array}$$

and

(36)$$-\underset{-\infty }{\overset{\infty }{\int }}\left\{G_{y^{\prime }}(\mathit{r}|s,0^{-})u^{\text{scat},-}(s)-G(\mathit{r}|s,0^{-})u_y^{\text{scat},-}(s)\right\}\text{d}s=\{\begin{array}{ll}0,\hfill & \mathit{r}\in {\mathrm{\Omega }}_{+},\hfill \\ u^{\text{scat}}(\mathit{r}),\hfill & \mathit{r}\in {\mathrm{\Omega }}_{-}.\hfill \end{array}$$

Writing the the normal derivatives of G and u^scat in terms of their limit values from above and below Γ, it follows that

(37)$$\begin{array}{l}{G}_{y^{\prime }}(\mathit{r}|s,0^{-})=A(s)G(\mathit{r}|s,0^{-})+B(s)G(\mathit{r}|s,0^{+}),\\ G_{y^{\prime }}(\mathit{r}|s,0^{+})=-B(s)G(\mathit{r}|s,0^{-})-A(s)G(\mathit{r}|s,0^{+}),\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}_y^{\text{scat},-}(s)=A(s)u^{\text{scat},-}(s)+B(s)u^{\text{scat},+}(s)-f_{-}^{\text{inc}}(s)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{u}_y^{\text{scat},+}(s)=-B(s)u^{\text{scat},-}(s)-A(s)u^{\text{scat},+}(s)-f_{+}^{\text{inc}}(s),\end{array}$$

where

(38)$$A(s)=ik\left\{\frac{\alpha (s)+\beta (s)}{2}\right\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}B=ik\left\{\frac{\alpha (s)-\beta (s)}{2}\right\}$$

and $f_{\pm }^{\text{inc}}$ defined in (11). Appropriately combining these expressions we arrive at

(39)$$\begin{array}{ll}{G}_{y^{\prime }}(\mathit{r}|s,0^{-})u^{\text{scat},-}(s)-G(\mathit{r}|s,0^{-})u_y^{\text{scat},-}(s)\hfill & =G(\mathit{r}|s,0^{-})f_{-}^{\text{inc}}(s)+\hfill \\ \hfill & B(s)(G(\mathit{r}|s,0^{+})u^{\text{scat},-}(s)-G(\mathit{r}|s,0^{-})u^{\text{scat},+}(s))\hfill \end{array}$$

(40)$$\begin{array}{ll}{G}_{y^{\prime }}(\mathit{r}|s,0^{+})u^{\text{scat},+}(s)-G(\mathit{r}|s,0^{+})u_y^{\text{scat},+}(s)\hfill & =G(\mathit{r}|s,0^{+})f_{+}^{\text{inc}}(s)+\hfill \\ \hfill & B(s)\left(G(\mathit{r}|s,0^{+})u^{\text{scat},-}(s)-G(\mathit{r}|s,0^{-})u^{\text{scat},+}(s)\right)\hfill \end{array}$$

Finally, from the identities above, and adding (35) and (36), the integral representation formula (10) for the field u^scat solution of the boundary value problem (8) is obtained.

B. Sommerfeld-integral Green’s function approximation

This appendix in devoted to the derivation of a Sommerfeld-integral [37] 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

(41)$$G^{\text{inc}}(\mathit{r}|{\mathit{r}}^{\prime }):=\frac{i}{4}{H}_0^{(1)}\left(k\left|(\mathit{r}-{\mathit{r}}^{\prime }\right|\right)=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\int }}\frac{{\text{e}}^{i{k}_x\left(x-x^{\prime }\right)+i{k}_y\left|y-y^{\prime }\right|}}{k_y}\text{d}{k}_x.$$

The square root $k_y=k_y(k_x):=\sqrt{k^2-k_x^2}$ is defined in the complex k_x plane as the product $\sqrt{k-k_x}\sqrt{k+k_x}$ 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. Figure 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.

 

Fig. 8 Signs of the real and imaginary parts of square root $k_y(k_x)=\sqrt{k_x^2-k^2}$ in the complex k_x plane.

Download Full Size | PDF

Using G^inc as incident field, the total field—which corresponds to G_p—can be expressed as

(42a)$$G_{\text{p}}(\mathit{r}|{\mathit{r}}^{\prime })=\{\begin{array}{ll}{G}^{\text{inc}}(\mathit{r}|{\mathit{r}}^{\prime })+G^r(\mathit{r}|{\mathit{r}}^{\prime }),\hfill & \mathit{r}\in {\mathrm{\Omega }}_{+},{\mathit{r}}^{\prime }\in {\mathrm{\Omega }}_{+},\hfill \\ G^t(\mathit{r}|{\mathit{r}}^{\prime }),\hfill & \mathit{r}\in {\mathrm{\Omega }}_{-},{\mathit{r}}^{\prime }\in {\mathrm{\Omega }}_{+},\hfill \end{array}$$

where the reflected and transmitted fields admit the integral representations

(42b)$$G^r(\mathit{r}|{\mathit{r}}^{\prime }):=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\oint }}R(k_x,x^{\prime })\frac{{\text{e}}^{i{k}_x(x-x^{\prime })+i{k}_y\left|y+y^{\prime }\right|}}{k_y}\text{d}{k}_x,$$

(42c)$$G^t(\mathit{r}|{\mathit{r}}^{\prime }):=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\oint }}T(k_x,x^{\prime })\frac{{\text{e}}^{i{k}_x(x-x^{\prime })+i{k}_y\left|y+y^{\prime }\right|}}{k_y}\text{d}{k}_x,$$

in terms of the reflection and transmission coefficient R and T defined in (13) and (16), respectively. The special integral sign “∮” 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 ∮ refers to a contour integral along the path SC depicted in Fig. 9.

 

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

Download Full Size | PDF

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

(43)$$\begin{array}{l}i{k}_y(1-R-T)=k\alpha (1+R+T),\\ i{k}_y(1-R+T)=k\beta (1+R-T),\end{array}$$

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

(44)$$G_{\text{p}}(\mathit{r}|{\mathit{r}}^{\prime })=\{\begin{array}{cc}\hfill G^t(\mathit{r}|{\mathit{r}}^{\prime }),\hfill & \hfill \mathit{r}\in {\mathrm{\Omega }}_{+},{\mathit{r}}^{\prime }\in {\mathrm{\Omega }}_{-},\hfill \\ G^{\text{inc}}(\mathit{r}|{\mathit{r}}^{\prime })+G^r(\mathit{r}|{\mathit{r}}^{\prime }),\hfill & \hfill \mathit{r}\in {\mathrm{\Omega }}_{-},{\mathit{r}}^{\prime }\in {\mathrm{\Omega }}_{-}.\hfill \end{array}$$

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_i^e(s,\sigma )=G_{\text{p}}(s,0^e|\sigma ,0^i)$ and ${\partial }_y{g}_i^e(s,\sigma )=G_{\text{p},\text{y}}(s,0^e|\sigma ,0^i)$ for s, σ ∈ ℝ, 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

(45a)$$g_i^e(s,\sigma )=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\oint }}{\widehat{g}}_i^e(k_x,\sigma ){\text{e}}^{i{k}_x(s-\sigma )}\text{d}{k}_x\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}{\partial }_y{g}_i^e(s,\sigma )=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\oint }}{\partial }_y{\widehat{g}}_i^e(k_x,\sigma ){\text{e}}^{i{k}_x(s-\sigma )}\text{d}{k}_x,$$

where

(45b)$$\begin{array}{l}{\widehat{g}}_{+}^{+}(k_x,\sigma )={\widehat{g}}_{-}^{-}(k_x,\sigma )=\frac{1}{k_y}\left\{1+R(k_x,\sigma )\right\}=\frac{1}{k_y}\left\{2-{\widehat{\mathrm{\Gamma }}}_{\alpha }(k_x,\sigma )-{\widehat{\mathrm{\Gamma }}}_{\beta }(k_x,\sigma )\right\},\\ {\widehat{g}}_{+}^{-}(k_x,\sigma )={\widehat{g}}_{-}^{+}(k_x,\sigma )=\frac{1}{k_y}T(k_x,\sigma )=\frac{1}{k_y}\left\{-{\widehat{\mathrm{\Gamma }}}_{\alpha }(k_x,\sigma )+{\widehat{\mathrm{\Gamma }}}_{\beta }(k_x,\sigma )\right\},\\ {\partial }_y{\widehat{g}}_{+}^{+}(k_x,\sigma )=-{\partial }_y{\widehat{g}}_{-}^{-}(k_x,\sigma )=i\left\{1+R(k_x,\sigma )\right\}=2i-i\left\{{\widehat{\mathrm{\Gamma }}}_{\alpha }(k_x,\sigma )+{\widehat{\mathrm{\Gamma }}}_{\beta }(k_x,\sigma )\right\}\\ {\partial }_y{\widehat{g}}_{+}^{-}(k_x,\sigma )=-{\partial }_y{\widehat{g}}_{-}^{-}(k_x,\sigma )=-iT(k_x,\sigma )=i\left\{{\widehat{\mathrm{\Gamma }}}_{\alpha }(k_x,\sigma )-{\widehat{\mathrm{\Gamma }}}_{\beta }(k_x,\sigma )\right\},\end{array}$$

with Γ̂_α and Γ̂_β being defined in (14).

Next we introduce the boundary integral operators

(46)$$(S_i^e\varphi )(s)=\underset{-\infty }{\overset{\infty }{\int }}{g}_i^e(s,\sigma )\varphi (\sigma )\text{d}\sigma \hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left({\partial }_y{S}_i^e(\varphi )\right)(s)=\underset{-\infty }{\overset{\infty }{\int }}{\partial }_y{g}_i^e(s,\sigma )\varphi (\sigma )\text{d}\sigma ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}s\in ℝ,$$

(that must be interpreted in the sense of distributions) which arise when taking limits of (21) and its normal derivative on Γ. From the integral representation (21) it follows that

(47)$$\begin{array}{l}\left[\left[u^{\text{scat}}\right]\right](s)=\left(S_{+}^{+}-S_{+}^{-}\right)(\phi +\psi )(s),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\left\{u^{\text{scat}}\right\}\right\}\left(s\right)=\left(S_{+}^{+}+S_{+}^{-}\right)(\phi -\psi )(s),\\ \left[\left[u_y^{\text{scat}}\right]\right](s)=\left({\partial }_y{S}_{+}^{+}-{\partial }_y{S}_{+}^{-}\right)(\phi -\psi )(s),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\left\{u_y^{\text{scat}}\right\}\right\}\left(s\right)=\left({\partial }_y{S}_{+}^{+}+{\partial }_y{S}_{+}^{-}\right)(\phi +\psi )(s),\end{array}$$

in terms of the integral operators (46), where we have utilized the identities $S_{+}^{+}=S_{-}^{-}$, $S_{+}^{-}=S_{-}^{+}$, ${\partial }_y{S}_{+}^{+}=-{\partial }_y{S}_{-}^{-}$ and ${\partial }_y{S}_{+}^{-}=-{\partial }_y{S}_{-}^{+}$ that result directly from (45). The uncoupled system of SIEs

(48a)$$\left\{{\partial }_y{S}_{+}^{+}-{\partial }_y{S}_{+}^{-}+ik\alpha (s)\left(S_{+}^{+}+S_{+}^{-}\right)\right\}{\mu }_1(s)=-2ik\alpha (s)u^{\text{inc}}(s,0)$$

(48b)$$\left\{{\partial }_y{S}_{+}^{+}+{\partial }_y{S}_{+}^{-}+ik\beta (s)\left(S_{+}^{+}-S_{+}^{-}\right)\right\}{\mu }_2(s)=-2u_y^{\text{inc}}(s,0)$$

for the new density functions μ₁ = φ − ψ and μ₂ = φ + ψ is thus obtained by substituting (47) in the transition conditions (8b)–(8c).

To show that SIE system (48) is in fact of the second-kind, we need to further study the properties of the integral operators on the left-hand-side of (48). Such properties can be determined from the regularity of the integral kernels $g_i^e$ and ${\partial }_y{g}_i^e$ which can in turn be derived from decay estimates for their Fourier transforms in (45). Consequently, utilizing the properties of the Fourier transform, it can be shown that the integral kernels in (48) can be expressed as

(49)$$\begin{array}{l}{\partial }_y{g}_{+}^{+}(s,\sigma )-{\partial }_y{g}_{+}^{-}(s,\sigma )+ik\alpha (s)\left(g_{+}^{+}(s,\sigma )+g_{+}^{-}(s,\sigma )\right)=-{\delta }_s+K_1(s,\sigma ),\\ {\partial }_y{g}_{+}^{+}(s,\sigma )+{\partial }_y{g}_{+}^{-}(s,\sigma )+ik\beta (s)\left(g_{+}^{+}(s,\sigma )-g_{+}^{-}(s,\sigma )\right)=-{\delta }_s+K_2(s,\sigma ),\end{array}$$

where δ_s denotes the Dirac delta distribution supported at s and

(50a)$$K_1(s,\sigma ):=\frac{k}{2\pi }\left\{\alpha (\sigma )-\alpha (s)\right\}\underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_x(s-\sigma )}}{k_y+k\alpha (\sigma )}\text{d}{k}_x,$$

(50b)$$K_2(s,\sigma ):=\frac{k}{2\pi }\left\{\beta (\sigma )-\beta (s)\right\}\underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_x(s-\sigma )}}{k_y+k\beta (\sigma )}\text{d}{k}_x.$$

Using the properties of the Dirac delta distribution and defining the integral operators

(51)$${\text{T}}_j[\mu ](s):=\underset{-\infty }{\overset{\infty }{\int }}{K}_j(s,\sigma )\mu (\sigma )\text{d}\sigma ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}s\in ℝ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}j=1,2,$$

we finally conclude that the relations in Eq. (48) can be equivalently expressed, in abstract form, as the SIEs in Eq. (22) for the μ_j, j = 1, 2.

The key fact about T_j for its use in our series of corrections is that T_j → 0 as ε → 0. We demonstrate this numerically in Appendix E. Analytically it occurs because the coefficients α(σ) − α(s) and β(σ) − β(s) vanish as ε → 0 for continuous functions α(x) = a(εx) and β(x) = b(εx). But, as discussed in Section 5, a technically challenging asymptotic analysis is required to rigorously demarcate the function spaces and norms for which T_j → 0 and to determine the rate of convergence, which we relegate to future work.

D. Generalized laws of reflection and transmission

This appendix is devoted to the derivation of the so-called generalized laws of reflection and transmission [28, 29] for general interface parameters for which standard Bloch–Floquet analysis [8] does not work. As it turns out, these laws can be derived from our zeroth-order approximation. To see this, consider metasurface parameters of the form

(52)$$\alpha (s)={\alpha }_0+\tilde{\alpha }(s)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\beta (s)={\beta }_0+\tilde{\beta }(s),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}s\in ℝ,$$

where α₀ and β₀ are constants and α̃ and β̃ are bounded functions. The total field resulting from the scattering of a planewave u^inc(r) = e^ikinc·r, k^inc = k(cos θ^inc, sin θ^inc), −π/2 ≤ θ^inc ≤ 0, off of the metasurface can be expressed as u = v + u₀ where

$$u_0(\mathit{r})=\{\begin{array}{ll}{\text{e}}^{i{\mathbf{k}}^{\text{inc}}\cdot \mathit{r}}+R_0({\theta }^{\text{inc}}){\text{e}}^{i{\mathbf{k}}^r\cdot \mathit{r}},\hfill & \mathit{r}\in {\mathrm{\Omega }}_{+},\hfill \\ T_0({\theta }^{\text{inc}}){\text{e}}^{i{\mathbf{k}}^{\text{inc}}\cdot \mathit{r}},\hfill & \mathit{r}\in {\mathrm{\Omega }}_{-},\hfill \end{array}$$

with k^r = k(cos θ^inc, − sin θ^inc) and

$$R_0(\theta ):=1-\frac{{\alpha }_0}{\left|\text{sin}\theta \right|+{\alpha }_0}-\frac{{\beta }_0}{\left|\text{sin}\theta \right|+{\beta }_0}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{T}_0(\theta ):=-\frac{{\alpha }_0}{\left|\text{sin}\theta \right|+{\alpha }_0}+\frac{{\beta }_0}{\left|\text{sin}\theta \right|+{\beta }_0},$$

is the total field resulting from the scattering of the planewave u^inc off of a metasurface with constant interface parameters α₀ and β₀. The field v = u − u₀, on the other hand, satisfies the Helmholtz equation in ℝ²\Γ, the radiation condition, and the transition conditions

$$\left[\left[v_y\right]\right]=-ik\alpha \left\{\left\{v\right\}\right\}-ik\tilde{\alpha }\left\{\left\{u_0\right\}\right\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left\{\left\{v_y\right\}\right\}=-ik\beta \left[\left[v\right]\right]-ik\tilde{\beta }\left[\left[u_0\right]\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{on}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{\Gamma }.$$

Letting then G₀ = G_p denote the Green’s function in Eqs. (42) and (44) corresponding to constant interface parameters α₀ and β₀, we obtain from the discussion in Section 4 that the zeroth-order approximation (i.e., Eq. (28a) with N = 0) of v is given by

(53)$$v_0(\mathit{r})=\underset{-\infty }{\overset{\infty }{\int }}\left\{G_0(\mathit{r}|\sigma ,0^{+}){\phi }_0(\sigma )-G_0(\mathit{r}|\sigma ,0^{-}){\psi }_0(\sigma )\right\}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{d}\sigma ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathit{r}\in {ℝ}^2\backslash \mathrm{\Gamma },$$

where letting

$$\tilde{A}(s)=ik\left\{\frac{\tilde{\alpha }(s)+\tilde{\beta }(s)}{2}\right\}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\tilde{B}(s)=ik\left\{\frac{\tilde{\alpha }(s)-\tilde{\beta }(s)}{2}\right\},$$

the approximate densities φ₀ and ψ₀ are given by

(54)$$\begin{array}{l}{\phi }_0(s)=-\left\{\tilde{A}(s)T_0({\theta }^{\text{inc}})+\tilde{B}(s)(1+R_0({\theta }^{\text{inc}}))\right\}{\text{e}}^{iks\text{cos}{\theta }^{\text{inc}}},\\ {\psi }_0(s)=\left\{\tilde{B}(s)T_0({\theta }^{\text{inc}})+\tilde{A}(s)(1+R_0({\theta }^{\text{inc}}))\right\}{\text{e}}^{iks\text{cos}{\theta }^{\text{inc}}}.\end{array}$$

Replacing G₀ in Eq. (53) by its far-field approximation, derived in Appendix F, we find that

(55)$$v_0(\mathit{r})~\frac{{\text{e}}^{ik\left|\mathit{r}\right|}}{\sqrt{\left|\mathit{r}\right|}}{v}_{\infty }(\theta ,{\theta }^{\text{inc}})\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|\mathit{r}\right|\to \infty ,$$

where, after some algebraic manipulations, the far-field pattern v_∞(θ|θ^inc) can be expressed as

(56a)$$v_{\infty }(\theta ,{\theta }^{\text{inc}})=\sqrt{\frac{2k}{\pi }}{\text{e}}^{i\frac{3\pi }{4}}\left|\text{sin}\left({\theta }^{\text{inc}}\right)\text{sin}\theta \right|\{\begin{array}{ll}{v}_{\infty }^{\alpha }\left(\theta ,{\theta }^{\text{inc}}\right)+v_{\infty }^{\beta }\left(\theta ,{\theta }^{\text{inc}}\right)\hfill & \theta \in (0,\pi ),\hfill \\ v_{\infty }^{\alpha }\left(\theta ,{\theta }^{\text{inc}}\right)-v_{\infty }^{\beta }\left(\theta ,{\theta }^{\text{inc}}\right)\hfill & \theta \in (-\pi ,0),\hfill \end{array}$$

where

(56b)$$\begin{array}{l}{v}_{\infty }^{\alpha }\left(\theta ,{\theta }^{\text{inc}}\right):={\int }_{-\infty }^{\infty }\frac{\tilde{\alpha }(\sigma ){\text{e}}^{ik\sigma \left(\text{cos}{\theta }^{\text{inc}}-\text{cos}\theta \right)}}{\left(\left|\text{sin}\theta \right|+{\alpha }_0\right)\left(\left|\text{sin}{\theta }^{\text{inc}}\right|+{\alpha }_0\right)}\text{d}\sigma \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\\ v_{\infty }^{\beta }\left(\theta ,{\theta }^{\text{inc}}\right):={\int }_{-\infty }^{\infty }\frac{\tilde{\beta }(\sigma ){\text{e}}^{ik\sigma \left(\text{cos}{\theta }^{\text{inc}}-\text{cos}\theta \right)}}{\left(\left|\text{sin}\theta \right|+{\beta }_0\right)\left(\left|\text{sin}{\theta }^{\text{inc}}\right|+{\beta }_0\right)}\text{d}\sigma .\end{array}$$

From Eq. (56) it thus follows that the far-field pattern v_∞ in Eq. (55) would correspond to a linear combination of planewaves with wavevectors k_α = k(cos θ_α, sin θ_α) and k_β = k(cos θ_β, sin θ_β) if $v_{\infty }^{\alpha }$ and $v_{\infty }^{\beta }$ were Dirac delta distributions supported at angles θ_α and θ_β, respectively. Formally, this can be achieved by selecting

$$\tilde{\alpha }(s)=c_{\alpha }{\text{e}}^{ik{d}_{\alpha }s}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\tilde{\beta }(s)=c_{\beta }{\text{e}}^{ik{d}_{\beta }s},$$

with c_α, c_β ∈ ℂ and d_α and d_β being such that

$$d_{\alpha }=\text{cos}{\theta }_{\alpha }-\text{cos}{\theta }^{\text{inc}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_{\beta }=\text{cos}{\theta }_{\beta }-\text{cos}{\theta }^{\text{inc}}.$$

Finally, the expressions in Eq. (19) are obtained by letting α₀ = β₀ = − sin θ^inc, c_α = ±α₀, and c_β = ∓β₀, for which a minimal reflection off of the metasurface is achieved.

E. On the convergence of the Neumann series in Eq. (25)

In our next appendix we consider an example to study the convergence of the Neumann series approximation in Eq. (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 in Eq. (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

(57)$$\frac{\Vert {\mathit{\mu }}_j-{\mathit{\mu }}_j^{(N)}\Vert }{\Vert {\mathit{\mu }}_j\Vert }\le \Vert T_j^{N+1}\Vert \le c{\left[\rho \left(T_j\right)\right]}^{N+1}$$

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

(58)$$\alpha (x)=a(\epsilon x):=\sum _{\ell =0}^{10}{c}_{\ell }{\text{e}}^{i\ell \epsilon x},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in ℝ,$$

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 α. Figure 7 displays the spectral radius of the matrix T₁ 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_{\text{p}}^{\text{ff}}$

This Appendix presents a detailed asymptotic analysis that establishes rigorously the relationship between G_p and $G_{\text{p}}^{\text{ff}}$, as well as the existence of guided modes.

To establish the relationship between $G_{\text{p}}^{\text{ff}}$ and G_p, we derive the far-field asymptotic approximation (as |r| → ∞) of the proto-Green’s function G_p given in Eqs. (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 [41, Chapter 8]. Similar saddle point calculations can also be found in classical references on layered media scattering, such as [37,45].

Assuming first that r′ ∈ Ω₊ and letting r = |r|(cos θ, sin θ), θ ∈ (−π, 0) ∪ (0, π), we have that the resulting reflected and transmitted fields in Eq. (42) can be expressed as

(59)$$G^r(\mathit{r}|{\mathit{r}}^{\prime })={\oint }_{\mathit{SC}}{q}_r(k_x,{\mathit{r}}^{\prime }){\text{e}}^{\left|\mathit{r}\right|\varphi (k_x)}\hspace{0.17em}\hspace{0.17em}\text{d}{k}_x\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{G}^t(\mathit{r}|{\mathit{r}}^{\prime })={\oint }_{\mathit{SC}}{q}_t(k_x,{\mathit{r}}^{\prime }){\text{e}}^{\left|\mathit{r}\right|\varphi (k_x)}\hspace{0.17em}\hspace{0.17em}\text{d}{k}_x$$

in terms of the phase and amplitude functions defined as

(60)$$\varphi (k_x):=i{k}_x\text{cos}\theta +i{k}_y\left|\text{sin}\theta \right|,$$

(61)$$q_r(k_x,{\mathit{r}}^{\prime }):=\frac{iR(x^{\prime },k_x)}{4\pi }\frac{{\text{e}}^{-i{k}_x{x}^{\prime }+i{k}_y{y}^{\prime }}}{k_y},$$

(62)$$q_t(k_x,{\mathit{r}}^{\prime }):=\frac{iT(x^{\prime },k_x)}{4\pi }\frac{{\text{e}}^{-i{k}_x{x}^{\prime }+i{k}_y{y}^{\prime }}}{k_y},$$

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 in Eq. (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 $k_y=k_y(k_x)=\sqrt{k^2-k_x^2}$.

We first consider the saddle points of ϕ, which correspond to solutions $k_x^{*}\in ℂ$ of the algebraic equation ${\varphi }^{\prime }\left(k_x^{*}\right)=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\text{cos}\theta =kx/\left|\mathit{r}\right|$ at which $\varphi \left(k_x^{*}\right)=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 ${\varphi }^{″}\left(k_x^{*}\right)=-i/\left(k{\text{sin}}^2\theta \right)\ne 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

(63a)$$\text{Im}{k}_x=\text{Re}{k}_x\left|\text{cot}\theta \right|-k\left|\text{csc}\theta \right|\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|k_x\right|\to \infty \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{if}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x>0,$$

and

(63b)$$\text{Im}{k}_x=\text{Re}{k}_x\left|\text{cot}\theta \right|+k\left|\text{csc}\theta \right|\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|k_x\right|\to \infty \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{if}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x<0.$$

Figure 9 depicts the steepest descent paths for x > 0 and x < 0. SD coincides with SC when x = 0.

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 Eqs. (13) and (16), respectively, we have that the poles of both q_r in Eq. (61) and q_t in Eq. (62)—which correspond to the poles of the functions Γ̂_α(x, x′) and Γ̂_β(x, x′) defined in Eq. (14)—are solutions of the (independent) algebraic equations

(64)$$k_y(k_x)=-k\alpha (x^{\prime })\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_y(k_x)=-k\beta (x^{\prime }).$$

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 Eqs. (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 = {Re k_x ≥ k, Im k_x ≥ 0} ∪ {Re k_x ≤ −k, Im k_x ≤ 0} (Fig. 8). Furthermore, since $\text{Im}{k}_y=\sqrt{k^2-k_x^2}\ge 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 in Eq. (63). Doing so we conclude that the relevant poles—that are henceforth denoted by $k_x^{(p)}$, p = 1, 2—must meet the conditions

(65a)$$k_x^{(1)}=k\sqrt{1-{\alpha }^2}\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{k}_x^2=k\sqrt{1-{\beta }^2},\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\hspace{0.17em}0\le \text{Im}{k}_x^{(p)}\le \text{Re}{k}_x^{(p)}\left|\text{cot}\theta \right|-k\left|\text{csc}\theta \right|\hspace{0.17em}\hspace{0.17em}\text{if}\hspace{0.17em}\hspace{0.17em}x>0,$$

and

(65b)$$k_x^{(1)}=-k\sqrt{1-{\alpha }^2}\hspace{0.17em}\hspace{0.17em}\text{or}\hspace{0.17em}\hspace{0.17em}{k}_x^2=-k\sqrt{1-{\beta }^2},\hspace{0.17em}\hspace{0.17em}\text{and}\hspace{0.17em}\text{Re}{k}_x^{(p)}\left|\text{cot}\theta \right|+k\left|\text{csc}\theta \right|\le \text{Im}{k}_x^{(p)}\hspace{0.17em}\le 0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{if}\hspace{0.17em}\hspace{0.17em}x<0.$$

Therefore, accounting for both saddle point and poles contributions, we obtain the following asymptotic expansions

(66a)$$G^r(\mathit{r}|{\mathit{r}}^{\prime })=\frac{{\text{e}}^{ik\left|\mathit{r}\right|-ik\mathit{r}\cdot {\overline{\mathit{r}}}^{\prime }/\left|\mathit{r}\right|+i\frac{\pi }{4}}}{\sqrt{8\pi k\left|\mathit{r}\right|}}R\left(\frac{kx}{\left|\mathit{r}\right|},x^{\prime }\right)+\sum _{p=1,2}{A}_p^r{\text{e}}^{i{k}_x^{(p)}\left|x-x^{\prime }\right|+i{k}_y\left(k_x^{(p)}\right)\left|y+y^{\prime }\right|}+O\left({\left|\mathit{r}\right|}^{-3/2}\right)$$

and

(66b)$$G^t(\mathit{r}|{\mathit{r}}^{\prime })=\frac{{\text{e}}^{ik\left|\mathit{r}\right|-ik\mathit{r}\cdot {\overline{\mathit{r}}}^{\prime }/\left|\mathit{r}\right|+i\frac{\pi }{4}}}{\sqrt{8\pi k\left|\mathit{r}\right|}}T\left(\frac{kx}{\left|\mathit{r}\right|},x^{\prime }\right)+\sum _{p=1,2}{A}_p^t{\text{e}}^{i{k}_x^{(p)}\left|x-x^{\prime }\right|+i{k}_y\left(k_x^{(p)}\right)\left|y-y^{\prime }\right|}+O\left({\left|\mathit{r}\right|}^{-3/2}\right)$$

as |r| → ∞, where r̄′ is the image point source r̄′ = (x′, −y′) (Fig. 2). The amplitudes $A_1^r$ and $A_1^t$ of the guided waves in Eq. (66) are directly obtained from the residues of q_t and q_r and read as

(67)$$A_1^r=A_1^t=-\frac{\alpha }{2\sqrt{1-{\alpha }^2}}$$

if a pole $k_x^{(1)}$ in Eq. (65) exists, and they equal zero otherwise. Similarly,

(68)$$A_2^r=-A_2^t=-\frac{\beta }{2\sqrt{1-{\beta }^2}}$$

if a pole $k_x^{(2)}$ in Eq. (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 half-planes. For example, the contribution of the pole $k_x^{(1)}=k\sqrt{1-{\alpha }^2}$ to the asymptotic expansion of the reflected and transmitted fields equals

(69)$$-\frac{\alpha }{2\sqrt{1-{\alpha }^2}}{\text{e}}^{ik\sqrt{1-{\alpha }^2}\left|x-x^{\prime }\right|-ik\alpha \left|y\pm y^{\prime }\right|},$$

with + and − corresponding to the reflected and transmitted fields, respectively.

 

Fig. 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 in Eq. (69) to propagate along the metasurface. The point source is located at distance 0.25λ above the metasurface. (a) $G_{\text{p}}^{\text{ff}}$ which in this case corresponds to the ray-optics approximation of the exact Green’s function G. (b) Rigorous far-field asymptotic approximation of G_p in Eq. (66). (c) Direct evaluation of the proto-Green’s function G_p defined in Eqs. (42) and (44).

Download Full Size | PDF

In order to establish the relationship between $G_{\text{p}}^{\text{ff}}$ and G_p we recall the asymptotic identity

(70)$$\frac{i}{4}{H}_0^{(1)}\left(k\left|\mathit{r}-{\mathit{r}}^{\prime }\right|\right)=\frac{{\text{e}}^{ik\left|\mathit{r}\right|-ik\mathit{r}\cdot {\overline{\mathit{r}}}^{\prime }/\left|\mathit{r}\right|+i\frac{\pi }{4}}}{\sqrt{8\pi k\left|\mathit{r}\right|}}+O{\left|\mathit{r}\right|}^{-3/2}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|\mathit{r}\right|\to \infty ,$$

that follows from |r − r′| = r − r · r′/|r| + O(|r|⁻¹) as |r| → ∞ and the standard asymptotic expansion of the Hankel function. Using the approximation in Eq. (70), the fact that (x − x′)/|r − r′| = x/|r| + O(|r|⁻¹) as |r| → ∞, and comparing Eq. (15) with Eq. (42), we conclude that

(71)$$G_{\text{p}}^{\text{ff}}\left(\mathit{r}|s,0^{\pm }\right)=G_{\text{p}}\left(\mathit{r}|s,0^{\pm }\right)+O\left({\left|\mathit{r}\right|}^{-3/2}\right)\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\text{as}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\left|y\right|\to \infty .$$

Therefore, $G_{\text{p}}^{\text{ff}}$ in Eq. (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 Γ̂_α(kx, s) and Γ̂_β(kx, s) defined in Eq. (14) have poles on the real axis which make surface wave modes in Eq. (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 Eqs. (42)–(44) and the integral kernels K₁ and K₂ in Eq. (50). This approach, which was originally developed for layered-media scattering problems [46] (see also [47, Section 2.3.5]), is a combination of the contour-integration method described in [48] and the the smooth-windowing approach put forth in [49] for the evaluation of oscillatory integrals.

Proto-Green’s function Consider the Sommerfeld integrals G^r and G^t in Eqs. (42)–(44) which are given by linear combinations of integrals of the form

(72)$$\mathrm{\Phi }(d_1,d_2)=\frac{i}{4\pi }\underset{-\infty }{\overset{\infty }{\oint }}\frac{k\gamma }{k_y+k\gamma }\frac{{\text{e}}^{i\left(k_x{d}_1+k_y{d}_2\right)/k}}{k_y}\text{d}{k}_x,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{d}_1{d}_2\ge 0,$$

where d₁ = k|x − x′|, d₂ = k|y + y′| in the case of G^r or d₂ = k|y − y′| in the case of G^t, and γ = α(x′) or γ = β(x′). Note that d_j, j = 1, 2, are dimensionless variables. To tackle the most challenging integration scenario, in what follows we consider the case when the integrand has a pole at $k\sqrt{1-{\gamma }^2}$ on the real axis.

Making use of the change of variable ξ = k_x/k, using the fact that k_y(k_x) = k_y(−k_x), and letting

(73)$$f(\xi )=\frac{i\gamma }{\pi }\frac{\text{cos}(\xi d_1)}{\sqrt{1-{\xi }^2+\gamma }}\frac{{\text{e}}^{i{d}_2\sqrt{1-{\xi }^2}}}{\sqrt{1-{\xi }^2}},$$

we have that Φ can be expressed as Φ = I₁ + I₂ where $I_j={\int }_{C_j}f(\xi )\text{d}\xi$, j = 1, 2, with C₁ and C₂ being the contours depicted in Fig. 11. The curve C₁ is a simple curve in the fourth quadrant that is parametrized by a smooth complex-valued function ζ : [0, π] ↦ ℂ satisfying ζ(0) = 0 and ζ(π) = L := 1 + ξ_p, where ξ_p is the pole of f at $\sqrt{1-{\gamma }^2}$. For the sake of definiteness, the curve C₁ is here selected as the semi-ellipse

(74)$$\zeta (t):=\left\{\frac{L\left(1+\text{cos}(t+\pi )\right)}{2}+iH\text{sin}(t+\pi ):t\in [0,\pi ]\right\}$$

that passes below all the singularities of the integrand f. The contour C₂, on the other hand, is simply the interval [L, ∞] on the real axis.

 

Fig. 11 Integration contours used in the evaluation of the proto-Green’s function G_p, and the integral kernels K₁ and K₂.

Download Full Size | PDF

Note that on C₁ the function f (ζ(t)) grows exponentially as t increases from 0 to π/2. Indeed, it can be shown [47] that $\left|{\text{e}}^{i{d}_2\sqrt{1-{\xi }^2}}\right|\le {\text{e}}^{d_{2}H}$ and |cos(d₁ξ)| ≤ e^d₁H for ζ ∈ C₁. Thus, in order to control the exponential growth of the integrand on C₁ we select H = (max{10, d₁ + d₂})⁻¹. This simple procedure ensures that the exponential terms of f remain bounded by one along C₁. The resulting expression for the contour integral I₁ is then approximated by means of the Clenshaw–Curtis quadrature rule [40]—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₁ and d₂, the number of quadrature points is chosen to grow linearly with d₁.

In order to evaluate the oscillatory integral I₂, on the other hand, we utilize the windowing method put forth in [49]. Using this procedure I₂ is approximated as

(75)$$I_2\approx {\int }_L^{A+L}f(t)w_A(t-L)\text{d}t,$$

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^∞(ℝ) function

(76)$$\eta (t;t_0,t_1)=\{\begin{array}{ll}1,\hfill & \left|t\right|\le t_0,\hfill \\ \text{exp}\left(\frac{2{\text{e}}^{-1/u}}{u-1}\right),\hfill & t_0<\left|t\right|<t_1,u=\frac{\left|t\right|-t_0}{t_1-t_0},\hfill \\ 0,\hfill & \left|t\right|>t_1,\hfill \end{array}$$

which equals one on [−t₀, t₀] and is supported on the (bounded) interval [−t₁, t₁].

In virtue of the oscillatory behavior of the integrand when d₁ ≠ 0, and the exponential decay of the integrand when d₂ ≠ 0, the integral on the right-hand-side of Eq. (75) converges to I₂ faster than any negative power $\sqrt{d_1^2+d_2^2}A$ as A goes to infinity—as proved in [47, Proposition 2.3.4]. In the special case d₁ = d₂ = 0, however, the f is slowly decaying on C₂ and does not oscillate, thus, it leads to slow (algebraic) convergence of the windowed-integral approximation in Eq. (75) to the integral I₂ as A → ∞. In fact, in that case the error in the windowed-integral approximation decays as O((cA)⁻¹) [47, Proposition 2.3.4]. Therefore, as in the case of I₁, the integral on the right-hand-side of Eq. (75) is here approximated by using Clenshaw–Curtis quadrature. The super-algebraic/exponential convergence of the windowed integral allows I₂ to be approximated with a fix accuracy and a fixed computational cost by choosing A inversely proportional to $\sqrt{d_1^2+d_2^2}$.

Integral kernels In order to numerically evaluate the integral kernels K₁ and K₂ defined in Eq. (50), we resort to the contour-integration procedure described above. The evaluation of these kernels requires the approximation of integrals of the form

(77)$$\mathrm{\Psi }(d)=\underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_{x}d/k}}{k_y+k\gamma }\text{d}{k}_x,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}d\ge 0,$$

where d = k|s − σ|, γ = α(σ) in the case of K₁ and γ = β(σ) in the case of K₂. Unlike the integral in Eq. (72), the Fourier integral in Eq. (77) is only conditionally convergent for d > 0 and diverges at d = 0. In order to separate the singularity of Ψ at d = 0 we use identity

(78)$$H_0^{(1)}(d)=\frac{1}{\pi }\underset{-\infty }{\overset{\infty }{\int }}\frac{{\text{e}}^{i{k}_{x}d/k}}{k_y}\text{d}{k}_x=\frac{1}{\pi }\frac{{\text{e}}^{i\xi d}}{\sqrt{1-{\xi }^2}}\text{d}\xi ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}d\ge 0,$$

which follows from Eq. (41), to obtain

(79)$$\mathrm{\Psi }(d)=\pi H_0^{(1)}(d)-k\gamma \underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_{x}d/k}}{k_y(k_y+k\gamma )}\text{d}{k}_x.$$

Thus, the integral in Eq. (79)—that turns out to be a continuous function of d—can now be evaluated directly by means of the contour integration procedure presented above for all d ≥ 0. We thus have $\mathrm{\Psi }(d)=\pi H_0^{(1)}(d)-\gamma (I_1+I_2)$ where $I_j={\int }_{C_j}f(\xi )\text{d}\xi$ with f now being given by $f(\xi )={\text{e}}^{i\xi d}/\left[\sqrt{1-{\xi }^2}\left(\sqrt{1-{\xi }^2}+\gamma \right)\right]$. Although f is absolutely integrable (it decays as ξ⁻² as |ξ| → ∞) a large value of A > 0 is needed in the windowed approximation of I₂ in Eq. (75) to achieve a prescribed accuracy when d = 0. In order to improve the slow O((cA)⁻¹) convergence rate as A → ∞ when d = 0, we note further that using the identity

(80)$$\underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_{x}d/k}}{k^2-k_x^2}\text{d}{k}_x=-\frac{i\pi }{2k}{\text{e}}^{i{k}_{x}d/k},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}d\ge 0,$$

which follows directly from Jordan’s lemma and Cauchy’s residue theorem, the integral in Eq. (79) can be expressed as

(81)$$\underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_{x}d/k}}{k_y(k_y+k\gamma )}\text{d}{k}_x=-\frac{i\pi }{2k}{\text{e}}^{ikd/k}-k\gamma \underset{-\infty }{\overset{\infty }{\oint }}\frac{{\text{e}}^{i{k}_{x}d/k}}{\left(k^2-k_x^2\right)(k_y+k\gamma )}\text{d}{k}_x.$$

Therefore, since the integrand on the right-hand-side of Eq. (81) decays as ξ⁻³ as |ξ| → ∞, the contour integration procedure described in the previous section yields super-algebraic convergence as A → ∞ for d > 0 and now yields an O((cA)⁻²) convergence rate when d = 0.

Integral operators and potentials Finally, we briefly mention that the (improper) oscillatory integrals in Eqs. (17), (21), (28), and in the definition of the operators T₁ and T₂ in Eq. (51), can be accurately truncated by means of the windowing procedure described above, in a manner similar to the windowed Green function method [50–52]. In order to handle the logarithmic singularity of the integral kernels in Eq. (51), on the other hand, standard singular integrations techniques, such as the spectrally accurate Martensen–Kussmaul quadrature rule [53], can be used.

Funding

Army Research Office and under Cooperative Agreement Number W911NF-18-2-0048.

References and links

1. E. F. Kuester, M. A. Mohamed, M. Piket-May, and C. L. Holloway, “Averaged transition conditions for electromagnetic fields at a metafilm,” IEEE Transactions on Antennas Propag. 51, 2641–2651 (2003). [CrossRef]  

2. C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O’Hara, J. Booth, and D. R. Smith, “An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials,” IEEE Antennas Propag. Mag. 54, 10–35 (2012). [CrossRef]  

3. C. L. Holloway, A. Dienstfrey, E. F. Kuester, J. F. O’Hara, A. K. Azad, and A. J. Taylor, “A discussion on the interpretation and characterization of metafilms/metasurfaces: The two-dimensional equivalent of metamaterials,” Metamaterials. 3, 100–112 (2009). [CrossRef]  

4. C. L. Holloway, E. F. Kuester, and A. Dienstfrey, “Characterizing metasurfaces/metafilms: The connection between surface susceptibilities and effective material properties,” IEEE Antennas Wirel. Propag. Lett. 10, 1507–1511 (2011). [CrossRef]  

5. S. A. Tretyakov, “Metasurfaces for general transformations of electromagnetic fields,” Proc. Royal Soc. Lond. A: Math. Phys. Eng. Sci. 373, 20140362 (2015).

6. S. N. Tcvetkova, D. H. Kwon, A. Díaz-Rubio, and S. A. Tretyakov, “Near-perfect conversion of a propagating plane wave into a surface wave using metasurfaces,” Phys. Rev. B 97, 115447 (2018). [CrossRef]  

7. A. Epstein and G. V. Eleftheriades, “Passive lossless Huygens metasurfaces for conversion of arbitrary source field to directive radiation,” IEEE Transactions on Antennas Propag. 62, 5680–5695 (2014). [CrossRef]  

8. A. Epstein and G. V. Eleftheriades, “Floquet-Bloch analysis of refracting Huygens metasurfaces,” Phys. Rev. B 90, 235127 (2014). [CrossRef]  

9. C. Pfeiffer and A. Grbic, “Metamaterial Huygens’ Surfaces: Tailoring Wave Fronts with Reflectionless Sheets,” Phys. Rev. Lett. 110, 197401 (2013). [CrossRef]  

10. K. Achouri, M. A. Salem, and C. Caloz, “General metasurface synthesis based on susceptibility tensors,” IEEE Transactions on Antennas Propag. 63, 2977–2991 (2015). [CrossRef]  

11. A. Díaz-Rubio, V. S. Asadchy, A. Elsakka, and S. A. Tretyakov, “From the generalized reflection law to the realization of perfect anomalous reflectors,” Sci. advances 3, e1602714 (2017). [CrossRef]  

12. V. S. Asadchy, M. Albooyeh, S. N. Tcvetkova, A. Díaz-Rubio, Y. Ra’di, and S. A. Tretyakov, “Perfect control of reflection and refraction using spatially dispersive metasurfaces,” Phys. Rev. B 94, 075142 (2016). [CrossRef]  

13. N. M. Estakhri and A. Alu, “Wave-front Transformation with Gradient Metasurfaces,” Phys. Rev. X 6, 041008 (2016).

14. B. Groever, W. T. Chen, and F. Capasso, “Meta-Lens Doublet in the Visible Region,” Nano Lett. 17, 4902–4907 (2017). [CrossRef]   [PubMed]  

15. B. Z. Katsenelenbaum, Theory of Nonuniform Waveguides: The Cross-Section Method (IET, 1998). [CrossRef]  

16. S. G. Johnson, P. Bienstman, M. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,” Phys. Rev. E 66, 066608 (2002). [CrossRef]  

17. D. Marcuse, Theory of Dielectric Optical Waveguides (Academic Press, 1974).

18. A. W. Snyder and J. Love, Optical Waveguide Theory (Springer, 2012).

19. D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, vol. 72 (SIAM, 2013). [CrossRef]  

20. J.-C. Nédélec, Acoustic and Electromagnetic Equations: Integral Representations for Harmonic Problems, vol. 144 (Springer, 2001). [CrossRef]  

21. E. Martini and S. Maci, “Metasurface transformation theory,” in “Transformation Electromagnetics and Metamaterials,” (Springer, 2014), pp. 83–116. [CrossRef]  

22. J. P. Wong, A. Epstein, and G. V. Eleftheriades, “Reflectionless wide-angle refracting metasurfaces,” IEEE Antennas Wirel. Propag. Lett. 15, 1293–1296 (2016). [CrossRef]  

23. M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science. 358, 1–8 (2017). [CrossRef]  

24. L. Verslegers, P. B. Catrysse, Z. Yu, W. Shin, Z. Ruan, and S. Fan, “Phase front design with metallic pillar arrays,” Opt. Lett. 35, 844–846 (2010). [CrossRef]   [PubMed]  

25. C. Pfeiffer, N. K. Emani, A. M. Shaltout, A. Boltasseva, V. M. Shalaev, and A. Grbic, “Efficient light bending with isotropic metamaterial Huygens’ surfaces,” Nano Lett. 14, 2491–2497 (2014). [CrossRef]   [PubMed]  

26. F. Aieta, P. Genevet, N. Yu, M. A. Kats, Z. Gaburro, and F. Capasso, “Out-of-plane reflection and refraction of light by anisotropic optical antenna metasurfaces with phase discontinuities,” Nano Lett. 12, 1702–1706 (2012). [CrossRef]   [PubMed]  

27. N. Yu, P. Genevet, F. Aieta, M. A. Kats, R. Blanchard, G. Aoust, J.-P. Tetienne, Z. Gaburro, and F. Capasso, “Flat optics: controlling wavefronts with optical antenna metasurfaces,” IEEE J. Sel. Top. Quantum Electron. 19, 4700423 (2013). [CrossRef]  

28. N. Yu, P. Genevet, M. A. Kats, F. Aieta, J. P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: Generalized laws of reflection and refraction,” Science. 334, 333–337 (2011). [CrossRef]   [PubMed]  

29. N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. 13, 139–150 (2014). [CrossRef]   [PubMed]  

30. R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, “Inverse design of large-area metasurfaces,” arXiv preprint arXiv:1808.04215 (2018).

31. T. Niemi, A. O. Karilainen, and S. A. Tretyakov, “Synthesis of polarization transformers,” IEEE Transactions on Antennas Propag. 61, 3102–3111 (2013). [CrossRef]  

32. M. Selvanayagam and G. V. Eleftheriades, “Polarization control using tensor Huygens surfaces,” IEEE Transactions on Antennas Propag. 62, 6155–6168 (2014). [CrossRef]  

33. C. Pfeiffer and A. Grbic, “Bianisotropic metasurfaces for optimal polarization control: Analysis and synthesis,” Phys. Rev. Appl. 2, 044011 (2014). [CrossRef]  

34. A. G. Voronovich, Wave Scattering From Rough Surfaces, vol. 17 (Springer, 2013).

35. D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather, Diffractive Optics: Design, Fabrication, and Test, vol. 62 (SPIE Press, 2004).

36. C. Cohen-Trannoudji, B. Diu, and F. Laloë, Quantum Mechanics. Volume II (Hermann, Paris, 1973).

37. W. C. Chew, Waves and Fields in Inhomogeneous Media, vol. 522 (IEEE pressNew York, 1995).

38. S. G. Johnson, M. Povinelli, M. Soljačić, A. Karalis, S. Jacobs, and J. Joannopoulos, “Roughness losses and volume-current methods in photonic-crystal waveguides,” Appl. Phys. B 81, 283–293 (2005). [CrossRef]  

39. M. Born and E. Wolf, Principles of Optics (Cambridge Univeristy Press, 1999), 7th ed. [CrossRef]  

40. P. J. Davis and P. Rabinowitz, Methods of Numerical Integration (Courier Corporation, 2007).

41. N. Bleistein, Mathematical Methods for Wave Phenomena (Academic Press, 2012).

42. A. Nosich, “Radiation conditions, limiting absorption principle, and general relations in open waveguide scattering,” J. Electromagn. Waves Appl. 8, 329–353 (1994). [CrossRef]  

43. J. A. DeSanto and P. A. Martin, “On angular-spectrum representations for scattering by infinite rough surfaces,” Wave Motion 24, 421–433 (1996). [CrossRef]  

44. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, vol. 55 (Courier Corporation, 1964).

45. L. Brekhovskikh, Waves in Layered Media, vol. 16 (Elsevier, 2012).

46. C. Pérez-Arancibia and O. P. Bruno, “High-order integral equation methods for problems of scattering by bumps and cavities on half-planes,” J. Opt. Soc. Am. A 31, 1738–1746 (2014). [CrossRef]  

47. C. Pérez-Arancibia, “Windowed integral equation methods for problems of scattering by defects and obstacles in layered media,” Ph.D. thesis, California Institute of Technology (2017).

48. M. Paulus, P. Gay-Balmaz, and O. Martin, “Accurate and efficient computation of the Green’s tensor for stratified media,” Phys. Rev. E 62, 5797 (2000). [CrossRef]  

49. J. A. Monro, “A Super-Algebraically Convergent, Windowing-Based Approach to the Evaluation of Scattering from Periodic Rough Surfaces,” Ph.D. thesis, California Institute of Technology (2008).

50. O. P. Bruno, M. Lyon, C. Pérez-Arancibia, and C. Turc, “Windowed Green function method for layered-media scattering,” SIAM J. on Appl. Math. 76, 1871–1898 (2016). [CrossRef]  

51. O. P. Bruno and C. Pérez-Arancibia, “Windowed Green function method for the Helmholtz equation in the presence of multiply layered media,” Proc. Royal Soc. Lond. A: Math. Phys. Eng. Sci. 473, 20170161 (2017). [CrossRef]  

52. O. P. Bruno, E. Garza, and C. Pérez-Arancibia, “Windowed Green function method for nonuniform open-waveguide problems,” IEEE Transactions on Antennas Propag. 65, 4684–4692 (2017). [CrossRef]  

53. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, vol. 93 (Springer, 2012), 3rd ed.