Diffractive Interface Theory: Nonlocal polarizability approach to the optics of metasurfaces

We present a formalism for understanding the elecromagnetism of metasurfaces, optically thin composite films with engineered diffraction. The technique, diffractive interface theory (DIT), takes explicit advantage of the small optical thickness of a metasurface, eliminating the need for solving for light propagation inside the film and providing a direct link between the spatial profile of a metasurface and its diffractive properties. Predictions of DIT are compared with full-wave numerical solutions of Maxwell's equations, demonstrating DIT's validity and computational advantages for optically thin structures. Applications of the DIT range from understanding of fundamentals of light-matter interaction in metasurfaces to efficient analysis of generalized refraction to metasurface optimization.


I. INTRODUCTION
Metasurfaces, optically thin structures with engineered diffraction, in the past few years have gained attention as a new platform for controlling the flow of light for ultra-compact photonics [1][2][3][4][5][6][7][8][9][10][11][12]. However, despite significant progress in fabrication and experimental characterization of metasurfaces, there is lack of analytical and numerical tools aimed at understanding and optimizing their optical behavior. As a result of the quasi-twodimensional geometry of metasurfaces, the majority of light-matter interaction occurs in the near-field proximity of the interfaces. At the same time, conventional techniques for numerical calculation of light interaction with composite systems (metamaterials) have been developed with volumetric materials in mind, [13][14][15] and as result, typically are inefficient in understanding the optics of metasurfaces. Likewise, traditional effective medium theories [16,17], which assign effective permittivities and permeabilities to volumetric structures cannot be used to describe the diffractive optics of metasurfaces that, in some sense, are more akin to the twodimensional graphene than to the traditional volumetricmetamaterials. Here we present a formalism that takes advantage of the small optical thickness of a metasurface and can successfully be used to predict its optical response. This technique can be used for development of analytical description of light-metasurface interaction and for rapid design and prototyping of metasurfaces for a wide range of optical applications.
Metasurfaces are typically based on arrays of planar optical resonators, fabricated at the interface between two dielectrics, arranged to produce diffracted waves of pre-designed amplitude and direction. Previous studies have focused on the induced phase discontinuity at a metasurface. Generalized Snell's Law [4,7] as well as standard diffraction theory [18] can be used to predict the * christopher roberts@student.uml.edu † viktor podolskiy@uml.edu direction of the diffracted beam. At the same time, calculation of the phase gradient or calculations of the amplitudes of the diffracted beam require full-wave solutions of Maxwell's equations. The vast majority of studies rely on finite element method (FEM) [14] or finite-difference time domain (FDTD) [15] techniques to solve for light interaction with resonant grating. However, while possible, these solutions require meshing the relatively large (multiple wavelength) regions of space with deep subwavelength (order of grating element) resolution, leading to resource-intensive and time-consuming calculations. An alternative numerical technique, rigorous coupled wave analysis (RCWA) [13], requires solving a large eigenvalue problem, once again making it impractical for design and optimization of real-life structures. More importantly, even when possible, calculation of field distribution is rarely sufficient for understanding the underlying physics which yields this distribution. Successful attempts to gain insight into basics of light interaction with elements of metasurfaces are limited to a few simple resonator shapes [1,4,11,12,19].
In this work, we present a general formalism for understanding the link between the spatial distribution of polarizability of a metasurface and the resulting diffraction performance of the structure. In contrast to the majority of existing analytical and numerical techniques, the proposed approach approximates the metasurface as FIG. 1. Schematic geometry of metasurface structure; DIT approximation is valid in regime h ≪ λ0 a thin polarizable sheet, and thus it takes explicit advantage of quasi-two-dimensional geometry. The spatial inhomogeneity of the polarization distribution results in nonlocal (wavevector-dependent) polarizability, which in turn results in multiple reflected and refracted (diffracted) beams [20,21]. Amplitudes of these beams can be related to the amplitude of the incident beam via generalized boundary conditions. Our technique mostly resembles the work that characterizes metasurfaces using generalized sheet transition currents [22]. However, this previously developed method relies on a priori knowledge of the optical response of 0 th order diffraction, a major limitation for the design of complex metasurfaces, is addressed by our proposed formalism. In some sense, our formalism merges RCWA, Transfer Matrix Method (TMM) [23,24], and Effective Medium Theories (EMT), resulting in a general approach to optics of optically thin diffractive systems: diffractive interface theory (DIT).

II. GENERAL FORMALISM
As mentioned above, we aim to derive the formalism for understanding reflection, refraction, and diffraction at an interface between two homogeneous materials covered with an optically thin inhomogeneous layer (Fig.1). We assume that all elements of the structure are nonmagnetic and have a linear responses to electric fields. Specifically, this linear response yields a polarization field inside the material with P = χE with E being the excitation field and χ being polarizability. In bulk materials, polarizability is related to permittivity via ǫ = 1 + 4πχ. (1) Permittivity, in turn, can be used to calculate the dispersion of the waves propagating inside such a material, as well as amplitudes and directions of these waves when the material is illuminated by incident light etc. In 2D films, however, the introduction of bulk permittivity is not meaningful. Therefore, we represent distribution of polarizability across the structure as Here, χ (b) represents polarizability of bulk components of the system, while χ (g) describes polarizability of material constructing the diffractive (grating) part of the structure, and h is the physical height of the grating layer. Note that both χ (b) and χ (g) correspond to volumetric polarizability of the materials comprising metasurface and, for implementation purposes, can be related to permittivity via Eq. (1). Assuming a harmonic timedependence E, H ∝ e −iωt , and the absence of free charges (∇ · E = −4π∇ · P), Maxwell's equations can be represented as resulting in As shown, for example, in Ref. [25], the set of waves propagating inside the system (modes of the structure) are completely determined by bulk components of the system through bulk [χ (b) ] component of polarizability. The polarization sheet χ (g) affects only boundary conditions.
It can be explicitly shown that the tangential components of both magnetic and electric fields are often discontinuous across a metasurface of finite, but small, thickness. The discontinuities of the magnetic field can be straightforwardly calculated by integrating Eq.(4) over "Amperian loops" in the xy and yz-planes, leading to where Similarly, integration of Eq.(5) yields discontinuities for the tangential components of the electric field: An alternate but equivalent derivation of these boundary conditions can be found in Ref. [26] III. DIT FOR PERIODIC METASURFACES We now apply the generalized boundary conditions presented above to derive diffractive interface theory (DIT), the formalism relating the refractive and diffractive properties of the periodic metasurface to its geometry and composition. We assume that the metasurface is periodic in the xy plane ( Fig.1) and is surrounded by two homogeneous materials. For simplicity, we will characterize the quantities referring to the materials above, below, and inside the grating by indices (b) = (1), (2), and (g), respectively. We assume that the structure is excited by a plane-wave with in-plane component of wavevector As with any periodic system, the Bloch theorem implies that electromagnetic field in the regions (1) and (2) can be represented as linear combination of (plane) waves with in-plane components of the wavevector separated by the multiple of the inverse lattice vector. Explicitly, the α Cartesian components of the field in the region (b) are written as Here, we utilize a common diagonalization technique used in implementations of RCWA and PWEM [27] in which the index j spans the domain of integers for onedimensional gratings and the domain of integer ordered pairs for two-dimensional gratings. The in-plane components of the wavevector of the waves are independent in the bulk regions (b) and are given by with the remaining component controlled by the dispersion of the bulk layers k It can be shown that out of six Cartesian components of the electromagnetic fields, only four components are independent; here we chose in-plane (xy) components of electric and magnetic fields as the independent fields which are used to enforce the boundary conditions [Eqs. (6,7)]. We further relate these four fields to the four sets of amplitudes of plane waves propagating in +z and −z direction with TE or TM polarization (see Appendix).
where (E, H, C ± ) are column-vectors with entries spanning the set of j indices.
Therefore, the problem of solving Maxwell's equations in the space surrounding the metasurface is reduced to the problem of finding the amplitudes a plane-wave expansion of electromagnetic fields in the bulk layers surrounding the metasurface.
The polarizability of the periodic metasurface can be expressed as The parameterχ (g) j ≡χ (g) (k j,x , k j,y ) describes the effective polarizability of the metasurface. Its nonlocal nature [20,21] (dependence on the wavevector) reflects long-range correlation and leads to diffraction, and allows the coupling of modes with different wavevectors during the interaction of light with the metasurface.
After a direct substitution of Eq.(10) into the Fourier transformations of Eqs.(6,7) and some straightforward mathematical transformations (see Appendix), we obtain the following block-matrix relationship: where I represent identity matrices, O represent matrices filled with zeros, X xy = − 2πiωh cχ (g) and X z = − 2πich ωχ (g) (I + 4πχ (g) ) −1 , with the mn element of matrixχ (g) given byχ gmn =χ gm−n , and K α being diagonal matrices with elements K αmm = k m,α .
Eq.(12) represents the main result of this work. It essentially represents an extension of transfer matrix formalism [23,24] to the domain of diffractive systems. As such, it can be used to • calculate the dispersion of the guided modes of the system by solving for dependence of k 0 (ω) that produces the solution with C (1)+ = C (2)− = 0; C (1)− , C (2)+ = 0 • calculate the diffraction (generalized refraction) of metasurfaces by computing amplitudes of diffracted fields C (1)− and C (2)+ as a function of incident fields C (1)+ and C (2)− • design, optimize, or retrieve the parameters of the metasurface by calculating the matrixχ g , and in the end, distribution χ (g) (x, y) that realizes the desired or observed optical performance We now illustrate these applications of the proposed formalism on several representative examples, aiming to test, and at the same time, demonstrate the benefits and limitations of diffractive interface theory. We begin by calculating the dispersion of the guided modes supported by a thin homogeneous metallic film. The solutions of Maxwell's equation in this simple, yet illustrative, system are well known. These represent symmetric and anti-symmetric combinations of surface plasmon polaritons. As the film becomes thinner, the mode corresponding to symmetric distribution of the magnetic field approaches light line; while the mode that has an anti-symmetric distribution of magnetic field becomes increasingly confined to the surface of the film, with its modal index n eff = ck x /ω ∝ λ 0 /h [28]. The behavior of these full-wave solutions of Maxwell's equations is illustrated in Fig.(2) with solid lines.
From the metasurface perspective, χ (g) (x, y) = (ǫ m − 1)/4π. As can be explicitly verified, the matrices K x , K y , X xy , and X z become diagonal, eliminating mixing between waves corresponding to different in-plane wavevectors. The terms corresponding to k 0 can thus be isolated, and the solution to DIT can be found analytically. The dispersion of the two waves corresponding to the two guided waves supported by the film is shown in Fig.(2) with dashed lines. It is clearly seen that in the limit as h ≪ λ 0 DIT predictions converge to exact solutions of Maxwell's equations. For the parameters chosen here (λ 0 = 1[µm], ǫ (1) = ǫ (2) = 1, ǫ (g) = −10 + 1i), DIT adequately describes the optical response of these modes when h λ 0 /20.

B. Generalized Refraction
In order to verify the validity and efficiency of the developed technique for a main class of applications of metasurfaces, which involve generalized refraction and polarization control, we compare the results of DIT to the predictions of RCWA [13]. Since both techniques rely on mode expansion, the disagreement between the two essentially reflects the implication of approximating the metasurface as an optically thin layer.
Techniques that rely on mode expansion require a certain number of modes in order to ensure convergence of the solution. In both DIT and RCWA, the number of required modes relates to the number of terms in discrete Fourier transform required to adequately represent the spatial profile of polarizability χ (g) (x, y). Our calculations indicate that DIT, presented in this work, converged as fast or faster than RCWA. The real advantage of the above formalism is computational complexity. The computational bottleneck for RCWA is the solution of an eigensystem [O( 4 3 n 3 )] that is required to calculate the propagation constants of the modes propagating through the periodic structure. At the same time, the limiting step of DIT is matrix division [O( 2 3 n 3 )]. The factor of two advantage can be dramatically improved with graphics processing unit (GPU) computing that strongly favors matrix division over eigenvalue problems [29,30]. In the end, our tests show that DIT runs almost an order of magnitude faster than RCWA for an equal number of modes, providing a compelling case for DIT for the design and optimization platform for optically thin metasurfaces.
To analyze the validity of DIT, we first study the optical response of two classes of metasurfaces, (i) arrays of nano-disks on dielectric substrates and (ii) arrays of nano-antennae. The behavior of the metasurfaces was studied in two different permittivity regimes, corresponding to highly metallic (ǫ = −2683 + 1367i, permittivity of gold in the mid-IR [31]) and plasmonic (ǫ = −10 + 1i, as achievable in the mid-IR with engineered highly-doped semi-conductors [32]) response of the metasurface material. In all calculations, the wavelength of excitation was fixed at λ 0 = 8µm and the period was fixed at Λ x = Λ y = 15.92µm.
We first analyze the utility of DIT to calculate the diffraction efficiency of metasurfaces. To achieve this goal, the diffraction efficiency of a metasurface comprised of nano-disks were studied as a function of the disc radius. The nano-disks are assumed to have a height of h = λ 0 /50, and are deposited at the interface of a dielectric [ǫ (2) = 10.8] and air. The system is excited from the air side at normal incidence. The diffraction efficiencies of reflected and transmitted modes are summarized in Fig.3. Overall, it is seen that DIT adequately represents the optical response of the structure. Minor deviations between DIT and RCWA predictions are seen at a high fill fraction.
Similar to FEM and FDTD calculations, DIT can be used to calculate not only diffraction efficiency, but also the field distribution across the system. However, in contrast to FEM, FDTD and similar techniques that calculate the total field, DIT is capable of analyzing the field distribution of a particular wave in the system. To illustrate this utility we analyze the response of the array of thin (h = λ 0 /100) nano-antennae surrounded by air [ Fig.4] for different lengths and angular orientation of the antennas and use the DIT to calculate the portion of the field corresponding to main reflection order right above the center of an antenna. These results were then compared to RCWA and shown in Fig.4. Once again, it is clearly seen that the DIT formalism adequately describes the field distribution in the system. Our calculations suggest a similar level of agreement in calculations of total field in this system and in the system where antennas are deposited at air-dielectric interface. For a comparable number of {L, θ} combinations, DIT runs almost an order of magnitude faster than RCWA.
It is widely known that efficiency of metasurfaces for generalized refraction and other related applications is rather small. In practice, only 10% of incident light is refracted into target diffraction order [4][5][6]. The efficiency of a metasurface, however, can be substantially enhanced when the metasurface is coupled to a homogeneous reflecting layer. In particular, antenna arrays incorporated into an optical stack with a highly reflective substrate have been recently suggested for applications in polarization conversion [33]. In this structure, the polarization of incident beam gets rotated via multiple refections between metasurface and the perfectly conducting metallic plate below the surface [see Fig.5].
While DIT correctly describes individual interaction of incident beam with a metasurface, the zero-thickness approximation introduces small errors into the optical response of the system. To understand the potential effect of accumulating of these small errors introduced by DIT, we consider the structure where the dielectric space between the metasurface and the PEC layer has variable loss (see Fig.5a). In this geometry, increase of the loss leads to the decrease of the number of successive metasurface-metal reflections. Results of these studies are shown in Fig.5(b-d). It is seen that when the number of successive reflections is relatively small, the pre- dictions of DIT are almost identical to those of RCWA. However, as the number of interactions is increased, the accuracy of DIT is reduced. The number of successive reflections/transmissions through the metasurface serves as a limiting factor for DIT applications C. Optimizing the metasurface Finally, we illustrate the potential of the proposed formalism in design and optimization of metasurfaces. For simplicity, here we optimize the parameters of 1D diffraction grating, maximizing the reflection into 1st diffraction order for a normal incident beam at λ 0 = 8µm. The unit cell of the grating is assumed to be represented as binary mask with fixed minimum width, w = Λ/N (mimicking gratings fabricated with lighography-based approaches). This particular grating design naturally lends itself to a genetic optimization algorithm where the binary profile of the period serves as the chromosome of the member of the population, and the fitness function represents the percent of light diffracted into 1st order. A highly metallic grating (ǫ = −2683 + 1367i), with a height of 10nm, on top of a dielectric substrate (ǫ = 10.8) was optimized for successively finer binary masks. The results of the op-

V. CONCLUSIONS
To conclude, we presented a novel formalism to describe the interaction of light with optically thin diffractive systems (metasurfaces): diffractive interface theory. The formalism, which takes explicit advantage of the quasi-two-dimensional nature of metasurfaces, provides a direct link between the spatial profile of polarizability of metasurface and its diffraction (generalized refraction) properties. In our tests, GPU-based implementations of DIT have been demonstrated to run almost 10 times faster than comparable implementations of RCWA. Applications of DIT for understanding the behavior of guided modes, calculations of generalized reflection of metasurfaces, calculations of field distributions in a system, and metasurface optimization have been demonstrated. The formalism can be straightforwardly extended to incorporate structures with anisotropic components, with strong magnetic response [19], as well as metasurfaces with non-planar geometries. This research is supported by NSF (grants ECCS-#1102183 and DMR-#1209761).