https://doi.org/10.1051/epjam/2016014
Research Article
Efficient integral equation modeling of scattering by a gradient dielectric metasurface
Department of Informatics, Aristotle University of Thessaloniki, 54124, Thessaloniki, Greece
^{*} email: ntsitsas@csd.auth.gr
Received:
10
October
2016
Accepted:
10
November
2016
Published online: 13 January 2017
Planewave scattering by a gradient dielectric metasurface, composed of a dielectric slab on which has been etched an infinitely periodic grating, is analyzed by means of rigorous semianalytical integral equation techniques based on the Method of Moments with entire domain basis functions. The electric field integral equation is considered with unknown function the electric field on the grating and is then solved by applying an entire domain Galerkin’s technique. The proposed methodology is characterized by high numerical stability and controllable accuracy; the obtained solution is analytic with the sole approximation of the final truncation of the expansion functions sets. The operation of the metasurface as a narrowband reflection frequency filter is investigated. Moreover, it is also shown that under certain conditions anomalous reflection and transmission can be generated even for a metasurface composed of conventional materials. The numerical results obtained provide design guidelines, which may pave the way for more systematic investigations of gradient dielectric metasurfaces assisting the control and manipulation of electromagnetic waves.
Key words: Gradient dielectric metasurfaces / Wave scattering / Integral equations / Method of moments / Narrow band filters / Anomalous reflection
© N.L. Tsitsas, Published by EDP Sciences, 2017
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Periodic dielectric grating structures have been traditionally employed in several areas of optics, electromagnetics, acoustics, elasticity, and integrated optoelectronics. Representative applications concern reflection and transmission frequency selective surfaces [1], integrated optics devices such as beamtosurfacewave couplers, distributed feedback amplifiers and lasers, multiplexers, and demultiplexers [2], and resonant filters utilized in lasercavity mode selectors [3] and security and anticounterfeit devices, used in credit and identification cards [4].
The remarkable properties of gradient metasurfaces with respect to controlling and manipulating electromagnetic waves have revived the interest for the investigations concerning the efficient analysis and design of periodic dielectric grating configurations. Efficiently designed gradient metasurfaces can control the propagation direction of an incident wave over subwavelength distances, modify the optical wavefront, bend and focus light, generate anomalous reflection and refraction phenomena as well as prescribed distributions of electromagnetic waves [5–9]. Besides, metamaterial gradient index diffraction gratings, composed of metallic and dielectric parts, may assist in gradient index optics applications [10], in multiband electromagnetic absorbers [11], and in the improvement of the nearfield optical enhancement [12]. More recently, it has been elaborated that dielectric materials with lowloss electromagnetic responses, realized by using completely transparent and highrefractiveindex dielectric building blocks, may offer the possibility of subdiffraction confinement and guiding of light without metals [13]. To this end, alldielectric gradient metamaterials have received particular attention and have been shown to support electromagnetic effects like mode conversion [14] and broadband anomalous reflection [15].
The scattering and diffraction phenomena by dielectric grating structures characterized by a periodic variation of their refractive index have been analyzed by various methodologies in the literature. Initially, the related investigations utilized perturbation techniques which provide physically intuitive but approximate results (valid only for weak grating perturbations, e.g. small grating’s thicknesses) [16]. The most widely used methods are based on differentialequation techniques and are the rigorous coupledwave analysis [17] and the modal method [18]. These methods are in certain cases computationally intensive, require timeconsuming root searches, large numerical matrices, and multiple timeintensive steps to determine the resonance’s width [3]. The finite difference time domain method [19] and the boundary element method [20] have also been applied for the analysis of scattering by dielectric grating configurations. The boundary element method, possesses, as a subdomain technique, drawbacks concerning the creation of fictitious charges and currents on the boundaries of the subdomains due to the discontinuous variation of the electromagnetic field as well as the large number of unknowns, leading to large required memory and CPU time and possible numerical instabilities [21]. A rigorous entire domain integral equation methodology for the modeling of plane wave scattering by dielectric grating waveguides was introduced in [22]. This methodology overcomes the above mentioned instabilities encountered in subdomain integral equation techniques.
In this work, the plane wave scattering problem by an infinitely periodic gradient dielectric metasurface is considered. The metasurface is composed of a dielectric slab on which has been etched a periodic rectangular grating. The applied methodology stems from appropriate proper modifications of the techniques utilized in [22]. The scattering problem is analyzed by a semianalytical integral equation methodology based on a Method of Moments technique using entire domain basis functions. The standard electric field integral equation is employed for the electric field on the grating’s grooves. The Green’s function of the nongrating problem is analytically expressed in the Fourier transform domain. The integral equation is solved by applying an entire domain Galerkin’s technique, utilizing a Fourier series expansion of the electric field on the grooves, leading to a nonhomogeneous linear system with respect to the Fourier series coefficients. The solution of this system provides the field on the grooves, by means of which the reflected and transmitted fields are determined. The methodology provides semianalytic solutions (the sole approximation is the truncation of the expansion functions sets) with high numerical stability, controllable accuracy, and high efficiency, since accurate results are obtained by using only a few expansion terms. In particular, it is shown that the beneficial characteristics of the methodology with respect to the efficiency and numerical stability are retained in the asymptotic region that the permittivities values of the metasurface’s materials are close to zero.
Special emphasis is given in the derived numerical results which reveal the suitable geometrical and physical parameters in order for the metasurface to possess interesting characteristics in its electromagnetic response. First, the metasurface’s efficient operation as a narrowband reflection frequency filter is analyzed. Certain characteristics of reflection resonances are described. It is demonstrated that the choices of the number and widths of the grating’s grooves inside the unit cell can efficiently control the resonances’ characteristics. The resonances’ dependence on the permittivity values of the grating’s grooves is investigated. The case of the grooves composed of an epsilonnearzero (ENZ) material is also considered. Second, the generation of anomalous reflection and transmission phenomena is examined; in such cases the main contribution in the reflected and transmitted fields is not in the directions predicted by Snell’s law but is confined in the same region of space with the incident field. It is shown that under certain conditions anomalous reflection and transmission can be generated even for a metasurface composed of an array of conventional materials.
2 Integral equation formulation of the scattering problem
The gradient dielectric metasurface under investigation is depicted in Figure 1. It is composed of a dielectric slab with relative dielectric permittivity ε_{d} and thickness 2d on which has been etched an infinitely Λperiodic rectangular dielectric grating composed of N rectangles (grooves) per unit cell with relative permittivity ε_{g}, thickness w, widths s_{n} and distances a_{n} from z = 0 (n = 1, …, N). The fact that the number N of the rectangles per unit cell and their respective characteristics a_{n} and s_{n} may be arbitrarily chosen offers additional degrees of freedom in designing the metasurface with respect to desired operational characteristics. The semiinfinite superstrate (cover) and substrate plane regions, above and below the slab, have relative permittivities ε_{c} and ε_{s}, respectively. The entire structure has constant magnetic permeability μ_{0} and is assumed uniform along the direction .
Figure 1.
The crosssection of the gradient dielectric metasurface under consideration, composed of a rectangular periodic grating with period Λ, thickness w and permittivity ε_{g}. A unitamplitude plane wave impinges on the metasurface at an angle θ. 
A TEpolarized incident plane wave of unit amplitude impinges on the metasurface from the cover region at an angle θ (see Figure 1). The incident electric field is given by (under exp(jωt) time dependence)(1)where k^{inc} is the incident wave vectorand k = ω/c the freespace wave number.
The electric field induced in every region is of the form(2)and the problem is reduced to the determination of the unknown scalar electric field factor Ψ. In the absence of external sources, timeharmonic solutions of Maxwell’s equations in the structure of Figure 1 are sought by considering the fields on the grating’s domain as equivalent polarization currents [21]. Thus, the factor Ψ admits the integral representation [22, 23](3)where Ψ_{0} is the electric field induced on the nongrating structure due to the plane incident wave (1), G the Green’s function of the nongrating structure with linesource excitation inside the dielectric slab (see Appendix A), and S_{g} is the total transverse crosssection of the rectangles.
According to the FloquetBloch theorem, the electric field factor satisfies the Bloch property [24] (pseudoperiodicity condition)which implies that(4)with u(x, z) a Λperiodic function of z.
Moreover, by representing the total grating’s cross section S_{g} as countable union of the cross sections of the grating’s rth unit celland employing the derived Fourier integral (A9) of the Green’s function as well as (4), we reformulate the integral representation (3) as(5)where function μ denotes the determined kernel of the Green’s function (see (A10)).
Now, the transformation ζ′ = z′ − rΛ reduces the double integrals of equation (5) to integrals on the basic unit cell S_{0} (the grating’s cross section on [0, Λ]). Then, by using the Poisson’s summation formula [25] for the Dirac functionand taking into account the basic propertydefinition of the Dirac function, we rewrite (5) as(6)
3 Solution of the integral equation by an entiredomain Galerkin technique
The integral equation (6) will be subsequently solved by applying a highly efficient entire domain Galerkin technique. First, the electric field’s factor u(x, z) on the grating’s basic unit cell S_{0} is expanded in the Fourier series with respect to z (7)where the Fourier coefficients φ_{n} are the transverse components of the spatial harmonics [26], satisfying here the ordinary differential equation(8)where
The functions φ_{n} are then expressed as linear combinations of the fundamental solutions of (8) (representing waves traveling parallel to the directions)(9)with under determination coefficients. In the terminology of the Method of Moments, the entire domain functions φ_{n}(x)exp[−j(2πn/Λ)z] are referred to as the electric field expansion (basis) functions.
Next, by substituting (7) into (6), we get(10)where the functions J_{p−n} and are defined in the Appendix B.
Furthermore, the observation vector (x, z) in (10) is restricted on the grating’s domain. Hence, for (x, z) ∈ S_{g}, the function u in the left hand side of (10) is expanded in the Fourier series (7). Then, by considering the inner products of both sides of (10) with the test functions (conjugates of the expansion functions of (7))and carrying out the resulting integrations, we formulate the following algebraic infinite nonhomogeneous square linear system of equations with respect to the unknown coefficients (11)where the involved infinite matrices and vectors are defined in the Appendix B.
This infinite system is solved numerically by truncation. Specifically, we take into account the terms of the expansion in (7) and the test functions in the inner products with maximum absolute order N_{t}, in order to reduce the infinite system to the (4N_{t} + 2) × (4N_{t} + 2) linear system(12)where the (2N_{t} + 1) × (2N_{t} + 1) matrices are given by(13)
( are defined in the Appendix B), c^{±} are 2N_{t} + 1 column vectors of the coefficients , b^{±} the 2N_{t} + 1 column vectors(14)while θ_{1}, B, and C are defined by (A4), (A7), and (A8), and
The required truncation order N_{t} is determined by applying a convergence control to the solutions for increasing N_{t}. It has been found that small values of N_{t} provide sufficient convergence, a fact which constitutes a basic advantage of the method; some representative convergence patterns for the metasurface under consideration are given in Section 5. This issue is elaborated in [22], where a detailed analysis concerning convergence and numerical efficiency aspects is included. Specifically, in Section 6.1 of [22] the method is validated by obtaining results which are in excellent agreement with those of [27] and [28], utilizing, respectively, a multiple reflection series method and a method of eigenfunction expansion in polynomial basis functions. To obtain these coincident results, the present method requires at most 2N_{t} + 1 = 5 expansion functions, while the methods of [27, 28] require NT = 31, N = 21 modes and 59 Legendre basis functions, respectively. This superior numerical efficiency of the present method is mainly justified by the facts that the unknown electric field factor and the entire domain expansion functions satisfy the same physical law, i.e. the Helmholtz equation, and that the layered medium Green’s function, utilized herein, satisfies inherently the boundary conditions in the nongrating structure.
4 Reflected and transmitted fields
After having determined the coefficients c^{±} by solving the nonhomogeneous linear system (12), we express the reflected Ψ^{r} and transmitted Ψ^{t} electric fields by means of the basic representation (10), as follows(15) (16)where the wavevectors’ components are expressed by (17) while the complex reflection and transmission coefficients are, respectively, given by(18) (19)
Coefficients R and T are given by (A5) and (A6), δ_{p0} is the Kronecker symbol (:δ_{00} = 1 and δ_{p0} = 0p ≠ 0), while(20) (21)
Hereafter, a specific field component (a specific term in each of the series (15) and (16)) indexed by p will be referred to as the preflected or ptransmitted diffracted order. The diffracted orders, which propagate along the xaxis, depend (according to (17)) on the period Λ, the wavelength λ, the permittivities ε_{c} and ε_{s} of the cover and substrate regions as well as on the angle of incidence θ. More precisely, the 0order fields are always propagating; notice also the contribution of coefficients R and T in (18) and (19) in the 0reflected and transmitted orders stemming from the fields induced on the nongrating structure. For orders p ≠ 0, we define the thresholds (22)for which holds:

for or and or the wavevectors’ components and are real and the porder reflected and transmitted fields are propagating along the xaxis.

for or and or the wavevectors’ components and become purely imaginary and the porder fields are evanescent.
Hence, and are the thresholds, for which the diffracted fields switch from propagation to cutoff; for a more detailed discussion as well as an illustrative example see [22].
5 Numerical results and discussion
The presented method provides semianalytic solutions with high numerical stability, controllable accuracy, high efficiency and without any restrictions on the problem’s parameters. Moreover, the Green’s function is analytically expressed and all the involved integrations are analytically carried out. Therefore, the method is also very efficient in terms of CPU time and computer memory and the sole approximation in the solution is the final truncation of the expansion functions sets. Besides, the present formulation requires no discretization of the integral equation involved, while the boundary element method requires and depends strongly on the discretization in boundary elements [20].
In the subsequent numerical results, for the determination of the required number N_{t} of expansion functions in (7), the following two criteria, in consistency with literature, have been employed: (i) an energy conservation condition [29] (reflected and transmitted fields must conserve power within 1 part in 10^{8}) and (ii) convergence to the solution with increasing N_{t} for all the grating and the incident wave parameters [17]. Some representative convergence patterns of the 0order reflection coefficient are depicted in Figures 2a and 2b corresponding to ε_{c} = ε_{s} = 1, 2d = w = 0.7 cm, ε_{d} = 2.59, θ = 60°, Λ = 1.65 cm, with ε_{g} = 0.5, 0.3, and 0.1, and (a) N = 1 groove with a_{1} = 0 and s_{1} = 0.5 Λ, and (b) N = 2 grooves with a_{1} = 0, a_{2} = 0.5 Λ, and s_{1} = 0.1 Λ, s_{2} = 0.3 Λ. It is observed that N_{t} = 2 expansion functions are sufficient in the case (a) of N = 1 groove (particularly there seems to be no observable difference in the curves of N_{t} = 2 and N_{t} = 3), while N_{t} = 5 functions are sufficient in the case (b) of N = 2 grooves. Hence, it is evident that the superior numerical efficiency of the applied integral equation methodology carries over to the present case of a metasurface composed of grooves with small positive permittivity (lying in the ENZ regime).
Figure 2.
Convergence patterns of the 0order reflection coefficient r_{0} with respect to the required number N_{t} of expansion functions in (7) for ε_{c} = ε_{s} = 1, 2d = w = 0.7 cm, ε_{d} = 2.59, θ = 60°, Λ = 1.65 cm, with ε_{g} = 0.5, 0.3, and 0.1, and (a) N = 1 groove with a_{1} = 0 and s_{1} = 0.5 Λ, and (b) N = 2 grooves with a_{1} = 0, a_{2} = 0.5 Λ, and s_{1} = 0.1 Λ, s_{2} = 0.3 Λ. 
For all the results presented hereafter, it has been checked that the choice N_{t} = 5 is sufficient so that both criteria are satisfied, a fact that validates the method’s superior numerical efficiency.
5.1 Reflection filters
Plane wave scattering by a metasurface composed of an array of dielectric rectangular cylinders can be treated as a special case of the developed method. More precisely, by selecting ε_{c} = ε_{s} and w = 2d, we investigate the structure of an infinite Λperiodic array composed of adjacent rectangular cylinders with permittivities ε_{d} and ε_{g}, the numbers and lengths of which per unit cell can be adjusted by the parameters N, a_{n}, and s_{n}.
First, we analyze the operation of the gradient dielectric metasurface as a frequency selective structure corresponding in particular to a reflection filter in the microwave regime. The cover and substrate regions are composed of air, i.e. ε_{c} = ε_{s} = 1, while the incidence is normal to the metasurface, i.e. θ = 90°. Figure 3 depicts the dependence of the amplitude r_{0} of the 0reflected order with respect to the operating frequency f for ε_{d} = 2.59 (Plexiglas), w = 2d = 0.7 cm, and Λ = 1.65 cm. Three cases are considered for the grating’s material: (i) ε_{g} = 2.05 (Teflon), (ii) ε_{g} = 0.5 (metamaterial 1: MTM1), and (iii) ε_{g} = 0.1 (metamaterial 2: MTM2). The results of Figure 3a correspond to a grating composed of one groove per unit cell, i.e. N = 1, with a_{1} = 0, and s_{1} = 0.5 Λ, while those of Figure 3b to a grating with three grooves per unit cell, i.e. N = 3, with a_{1} = 0, s_{1} = 0.2 Λ, a_{2} = 0.3 Λ, s_{2} = 0.1 Λ, a_{3} = 0.6 Λ, and s_{3} = 0.2 Λ.
Figure 3.
Amplitude r_{0} of the 0order reflection coefficient versus the operating frequency f for ε_{c} = ε_{s} = 1, 2d = w = 0.7 cm, ε_{d} = 2.59 (Plexiglas), θ = 90°, Λ = 1.65 cm, with ε_{g} = 2.05 (Teflon), 0.5 (Metamaterial 1: MTM1), 0.1 (Metamaterial 2: MTM2), and (a) N = 1 groove with a_{1} = 0 and s_{1} = 0.5 Λ, and (b) N = 3 grooves with a_{1} = 0, a_{2} = 0.3 Λ, a_{3} = 0.6 Λ, and s_{1} = 0.2 Λ, s_{2} = 0.1 Λ, s_{3} = 0.2 Λ. 
From Figure 3a it is observed that a grating of one groove composed of a metamaterial with ε_{g} = 0.5 increases the value of r_{0} at the location of the resonance compared with the respective value when the grating is composed of a conventional material (Teflon). Moreover, when the grating is composed of a metamaterial with ε_{g} = 0.1, then the situation becomes different and a total transmission resonance (zero reflection) is observed at the same location (where the total reflection resonances appeared for ε_{g} = 2.05 and ε_{g} = 0.5). Particularly, for ε_{g} = 0.1, small reflection is generated for a wide range of the examined frequencies (r_{0} < 0.2 in the interval from 17 to 20 GHz). Metamaterials with permittivities values like ε_{g} = 0.1 belong to the class of epsilonnearzero (ENZ) metamaterials; such materials have been shown to possess related remarkable electromagnetic properties, including, among others, controlling the radiation pattern [30] and tunneling of electromagnetic energy [31]. Effects of potential losses in the ENZ materials can also be included in the present integral equation modeling and investigated, although related lowloss experimental realizations have already been proposed [32].
For a metasurface composed of a grating with three grooves, by comparing Figures 3a and 3b, we observe that the locations of the reflection resonances are shifted significantly to the left with the particular values of the resonant frequencies depending on the permittivity of the metamaterial composing the grating. In addition, the respective values of r_{0} become almost 1, and hence these correspond to total reflection resonances. The physical mechanism generating such total reflection (as well as total transmission) resonances is strongly connected to the excitation of the waveguide modes and is discussed thoroughly in [29]. Besides, two additional (total reflection) resonances are generated in the case of N = 3, while the bandwidth of the resonances is broadened when a metamaterial grating is used.
Thus, the modulation inside each grating’s unit cell (number N and geometrical parameters a_{n} and s_{n} of the rectangles) combined with the proper selection of the grating’s material offer flexibility and additional degrees of freedom to the efficient designing of gradient metasurfaces with the desirable resonance characteristics. This offers an alternative route of adjusting the resonances behavior to those traditionally proposed in the literature, concerning mainly the addition of extra dielectric layers in the grating structure [1, 33].
5.2 Anomalous reflection and transmission
Now, we investigate and discuss suitable selections of parameters of the metasurface in order for it to generate the socalled anomalous reflection and transmission phenomena in which cases the main contribution in the reflected and transmitted fields is not in the directions predicted by Snell’s law but is confined in the same region of space with the incident field, hence forcing light to propagate toward the side of incidence. To this end, we require that the grating period Λ, the operating wavelength λ, and the angle of incidence θ are such that the respective structures admit only the 0 and the −1 orders propagating in the cover and substrate regions, which, according to the cases (i) and (ii) following (22), are equivalent to the conditions and , which when ε_{s} = ε_{c} and thus , lead to(23)
It has been recently shown experimentally that anomalous reflection effects can be generated by highcontrast metastructures (HCM) which consist of periodic structures of a single highrefractiveindex layer surrounded by a lowindex material, and a periodicity of nearly one wavelength [13, 15]. Such metastructures are designed in order to suitably manipulate the large index contrast with aim to annihilate the 0order of reflection (following Snell’s law) and simultaneously enhance the −1order which creates anomalous reflection at the metastructure interface. The physical mechanism justifying the anomalous reflection effects is related to the excitation of waveguide array modes each of which possesses a different phase after propagating through the HCM. As the modes return to the input plane, due to the strong mismatch to the incoming plane wave, these modes not only reflect back to themselves but also couple to each other [13, 15].
Here, as a first example of simulating structures, which can produce anomalous reflection phenomena, we consider the basic design parameters of the HCM presented in [15] concerning an array of rectangular cylinders composed of Silicon inside a certain homogeneous background composed of a lowindex material. Figure 4b depicts the amplitudes of the 0 and −1order reflection coefficients (being the only two diffracted orders propagating in the cover region, according to (23)) versus the angle of incidence θ corresponding to plane wave scattering by the array of rectangular cylinders, depicted in Figure 4a, with ε_{g} = 17.14 (Silicon), N = 1, a_{1} = 0, and s_{1} = 0.5 Λ, lying inside a background medium with ε_{c} = ε_{s} = ε_{d} = 1, and d = w = 180 nm. The operating wavelength is λ = 532 nm and the cylinders’ period Λ = 500 nm; a condition of the form Λ ~ λ is usually considered in order for the metastructure to support propagation of two and only two diffraction orders for a wide range of angles of incidence [15]. Actually, the HCM of [15] is fabricated as an array of pixels in a transparent, flexible membrane, while in the present considerations, we considered the array of Silicon cylinders to lie inside a homogeneous vacuum (ε_{c} = ε_{s} = ε_{d} = 1). It is evident that r_{–1} is larger than r_{0} in the range of θ which is mainly of interest (excluding the neargrazing incidence angles). Particularly, the locally maximum distances between r_{–1} and r_{0} are observed at θ = 39° and at θ = 75°; it is not of much interest to examine angles of incidence very near to normal because then anomalous reflection effects cannot be clearly distinguished and exploited.
Figure 4.
(a) Gradient dielectric metasurface (infinite Λperiodic array of dielectric rectangular cylinders) with ε_{c} = ε_{s} = ε_{d} = 1, d = w = 180 nm, ε_{g} = 17.14 (Silicon), λ = 532 nm, Λ = 500 nm, N = 1, a_{1} = 0, and s_{1} = 0.5 Λ. (b) Amplitudes r_{0} and r_{–1} of the 0 and −1order reflection coefficients versus the angle of incidence θ for the metasurface described in (a). 
Next, we investigate alternative ways of enhancing the relative diffraction efficiency between r_{–1} and r_{0} (leading to anomalous reflection) and the range of incident angles for which this effect appears. We consider a Silicon dielectric slab (ε_{d} = 17.14, d = 180 nm) possessing a grating composed of Teflon (ε_{g} = 2.05, d = w) grooves (see Figure 5a). The ratio of the wavelength over the periodicity is taken to be exactly 1 (λ = Λ = 532 nm), and thus, according to (23), only the 0 and −1 diffracted orders propagate. The results, depicted in Figure 5b, indicate that r_{–1} becomes larger than 0.9 and simultaneously is nearly 4.5 times greater than the respective value of r_{0}. It is also worth to note that this very high diffraction efficiency of the −1 reflected order is attained for a metasurface composed of an array of conventional materials (Silicon slab filled with Teflon grooves).
Figure 5.
(a) Metasurface with ε_{c} = ε_{s} = 1, d = w = 180 nm, ε_{d} = 17.14 (Silicon), λ = Λ = 532 nm, N = 1, a_{1} = 0, s_{1} = 0.5 Λ, and ε_{g} = 2.05 (Teflon). (b) Amplitudes r_{0} and r_{−1} of the 0 and −1order reflection coefficients versus the angle θ for the metasurface of (a), (c) same as (b), but for ε_{g} = 0.7. 
Still, the diffraction efficiency of r_{−1}, observed in Figure 5b, occurs for a relatively narrow interval of incident angles. The range of angles for which the anomalous reflection effect appears can be broadened if a different material is utilized in the grating region. In Figure 5c, the case of a Silicon slab filled with grooves of ε_{g} = 0.7 is depicted. Now, the higher diffraction efficiency of the −1reflected order compared with the one of the 0order occurs for a significantly wider interval of incident angles than that of Figure 5b. However, a reduction of the diffraction efficiency of r_{−1} by nearly 25% is observed.
Finally, it is also important to investigate related anomalous transmission effects in the gradient dielectric metasurfaces under consideration. To this end, we plot in Figures 6b–6e the amplitudes r_{0}, r_{−1}, t_{0}, and t_{−1} of the 0 and −1order reflection and transmission coefficients versus the operating wavelength λ for the array of dielectric cylinders, depicted in Figure 6a, with four different permittivities ε_{g} < 1 having the same thickness w = 2d = Λ = 800 nm, and lying in free space (ε_{c} = ε_{s} = ε_{d} = 1). Figure 6b shows that, when ε_{g} is in the ENZ regime, very high diffraction efficiency of the −1transmission order compared to the 0transmission order is attained, while simultaneously both the 0 and −1reflection orders have very small reflection coefficients. So, in this case, anomalous transmission, accompanied by very small reflection, is obtained by using an array of rectangular cylinders with ε_{g} = 0.2. As ε_{g} increases to 0.3, the diffraction efficiency of the −1transmission order is still higher than the 0transmission order but their difference starts to reduce (see Figure 6c). The situation is reversed in Figure 6d, corresponding to ε_{g} = 0.5, where the 0transmission order becomes larger than the −1order. As the permittivity ε_{g} of the cylinders becomes closer to 1 (that of free space), then the entire array starts, as expected, to exhibit the characteristics of a transparent layer, i.e. nearly zero reflection in both the 0 and −1orders, suppression of the transmission of the −1order and nearly one transmission coefficient of the 0order (for ε_{g} = 1 the entire structure is vacuum and hence the incident wave travels unaltered). These conclusions are verified by the results of Figure 6e.
Figure 6.
(a) Geometry of the gradient dielectric metasurface (ε_{c} = ε_{s} = ε_{d} = 1, Λ = 800 nm, d = Λ/2, w = 2d, N = 1, a_{1} = 0, s_{1} = 0.5 Λ, and θ = 60°). (b)–(e) Reflection and transmission coefficients r_{0}, r_{−1}, t_{0}, and t_{−1} versus the wavelength λ for (b) ε_{g} = 0.2, (c) ε_{g} = 0.3, (d) ε_{g} = 0.5, and (e) ε_{g} = 0.7. 
6 Conclusions
Plane wave scattering by a gradient dielectric metasurface was investigated by a semianalytical entire domain integral equation methodology. The scattering problem was formulated by using the standard electric field integral equation with unknown function the electric field on the grating’s grooves. The integral equation was subsequently solved by applying a Galerkin’s technique. The basic advantages of the proposed method concern the high accuracy, the superior numerical efficiency and the analytic calculation of the Green’s function and all the integrals involved. Numerical results demonstrated that the efficient selection of the parameters of the metasurface can provide very useful characteristics in the design of reflection frequency filters and can generate anomalous reflection and transmission effects which can be exploited in a number of potential applications. It has been shown that these anomalous reflection and transmission effects can also occur in metasurfaces realized as arrays of cylinders lying in vacuum and composed of conventional materials. Further analysis of such effects by using the developed integral equation methodology in combination with optimization algorithms for the optimal selection of the physical and geometrical characteristics of the metasurface is expected to increase the applicability of the obtained results in a broader range of applications.
Implications and influences
Dielectric slabs with periodic variation of their refractive index have been traditionally utilized as key elements in fundamental components of various microwave and optical devices. This paper examines a gradient dielectric metasurface realized as a dielectric slab possessing periodic rectangular grooves. Recently, such alldielectric metasurfaces have been experimentally shown to possess remarkable properties concerning the efficient manipulation of light without using metals. The developed and presented integral equation methodology for the analysis of the scattering by the under consideration gradient dielectric metasurface offers a very accurate, fast and efficient numerical modeling tool. Moreover, the numerical results derived show the efficient use of such gradient dielectric metasurfaces in reflection frequency filters as well as in the generation of anomalous reflection and transmission effects. Especially, it is shown that the latter effects can be generated by using arrays of cylinders of conventional materials. Optimizing further the metasurface’s design by using the present integral equation analysis in order to achieve anomalous reflection and transmission effects for a wider range of incidence angles and frequencies is a very promising and challenging direction for future work.
References
 S. Tibuleac, R. Magnusson, T.A. Maldonado, P.P. Young, T.R. Holzheimer, Dielectric frequencyselective structures incorporating waveguide gratings, IEEE Trans. Microwave Theory Tech. 48 (2000) 553–561. [CrossRef] (In the text)
 S.S. Wang, R. Magnusson, Design of waveguide grating filters with symmetrical line shapes and low sidebands, Opt. Lett. 19 (1994) 919–921. [CrossRef] (In the text)
 S.M. Norton, T. Erdogan, G.M. Morris, Coupled mode theory of resonantgrating filters, J. Opt. Soc. Am. A 14 (1997) 629–639. [CrossRef] (In the text)
 J. Turunen, F. Wyrowski (Eds.), Diffractive optics for industrial and commercial applications, Chap. 12, Akademie, Berlin, 1997. (In the text)
 N. Yu, P. Genevet, M.A. Kats, F. Aieta, J.P. Tetienne, F. Capasso, Z. Gaburr, Light propagation with phase discontinuities: Generalized laws of reflection and refraction, Science 334 (2011) 333–337. [CrossRef] [PubMed] (In the text)
 X. Ni, N.K. Emani, A.V. Kildishev, A. Boltasseva, V.M. Shalaev, Broadband light bending with plasmonic nanoantennas, Science 335 (2011) 427. [CrossRef] [PubMed]
 C. Argyropoulos, G. D’Aguanno, N. Mattiucci, N. Akozbek, M.J. Bloemer, A. Alù, Matching and funneling light at the plasmonic Brewster angle, Phys. Rev. B 85 (2012) 024304. [CrossRef]
 F. Monticone, N.M. Estakhri, A. Alù, Full control of nanoscale optical transmission with a composite metascreen, Phys. Rev. Lett. 110 (2013) 203903. [CrossRef]
 N.M. Estakhri, C. Argyropoulos, A. Alù, Graded metascreens to enable a new degree of nanoscale light management, Phil. Trans. R. Soc. A 373 (2015) 20140351. [CrossRef] (In the text)
 Y.J. Tsai, S. Larouche, T. Tyler, G. Lipworth, N.M. Jokerst, D.R. Smith, Design and fabrication of a metamaterial gradient index diffraction grating at infrared wavelengths, Optics Express 19 (2011) 24411–24423. [CrossRef] (In the text)
 Y. Cui, et al., Plasmonic and metamaterial structures as electromagnetic absorbers, Laser Photon. Rev. 8 (2014) 495–520. [CrossRef] (In the text)
 A.A. Darweesh, S.J. Bauman, J.B. Herzog, Improved optical enhancement using doublewidth plasmonic gratings with nanogaps, Photon. Research 4 (2016) 173–180. [CrossRef] (In the text)
 S. Jahani, Z. Jacob, Alldielectric metamaterials, Nature Nanotechnol. 11 (2016) 23–36. [CrossRef] (In the text)
 D. Ohana, U. Levy, Mode conversion based on dielectric metamaterial in silicon, Optics Express 22 (2014) 27617–27631. [CrossRef] (In the text)
 L. Zhu, J. Kapraun, J. Ferrara, C.J. ChangHasnain, Flexible photonic metastructures for tunable coloration, Optica 2 (2015) 255–258. [CrossRef] (In the text)
 S. Miyanaga, T. Akasura, Intensity profiles of outgoing beams from tapered grating couplers, Radio Sci. 17 (1982) 135–143. [CrossRef] (In the text)
 M.G. Moharam, E.B. Grann, D.A. Pommet, T.K. Gaylord, Formulation for stable and efficient implementation of the rigorous coupledwave analysis of binary gratings, J. Opt. Soc. Am. A 12 (1995) 1068–1076. [CrossRef] (In the text)
 R.S. Chu, J.A. Kong, Modal theory of spatially periodic media, IEEE Trans. Microwave Theory Tech. MTT25 (1977) 18–24. [CrossRef] (In the text)
 M.R. Zunoubi, H.A. Kalhor, Diffraction of electromagnetic waves by periodic arrays of rectangular cylinders, J. Opt. Soc. Am. A 23 (2006) 306–313. [CrossRef] (In the text)
 Y. Nakata, M. Koshiba, Boundaryelement analysis of planewave diffraction from groovetype dielectric and metallic gratings, J. Opt. Soc. Am. A 7 (1990) 1494–1502. [CrossRef] (In the text)
 G. Athanasoulias, N.K. Uzunoglu, An accurate and efficient entiredomain basis Galerkin’s method for the integral equation analysis of integrated rectangular dielectric waveguides, IEEE Trans. Microwave Theory Tech. 43 (1995) 2794–2804. [CrossRef] (In the text)
 N.L. Tsitsas, N.K. Uzunoglu, D.I. Kaklamani, Diffraction of plane waves incident on a grated dielectric slab: An entire domain integral equation analysis, Radio Sci. 42 (2007) RS6S22. [CrossRef] (In the text)
 N.L. Tsitsas, D.I. Kaklamani, N.K. Uzunoglu, Rigorous integral equation analysis of nonsymmetric coupled grating slab waveguides, J. Opt. Soc. Am. A 23 (2006) 2888–2905. [CrossRef] (In the text)
 M. Weber, D.L. Mills, Interaction of electromagnetic waves with periodic gratings: Enhanced fields and the reflectivity, Phys. Rev. B 27 (1983) 2698–2709. [CrossRef] (In the text)
 P.M. Morse, H. Feshbach, Methods of Theoretical Physics, Part I, McGrawHill, Tokyo, 1953. (In the text)
 R.E. Collin, Field theory of guided waves, IEEE Press, New York, 1991. (In the text)
 D.M. Pai, K.A. Awada, Analysis of dielectric gratings of arbitrary profiles and thicknesses, J. Opt. Soc. Am. A 8 (1991) 755–762. [CrossRef] (In the text)
 R.H. Morf, Exponentially convergent and numerically efficient solution of Maxwell’s equations for lamellar gratings, J. Opt. Soc. Am. A 12 (1995) 1043–1056. [CrossRef] (In the text)
 H.L. Bertoni, L.H.S. Cheo, T. Tamir, Frequencyselective reflection and transmission by a periodic dielectric layer, IEEE Trans. Antennas Propagat. 37 (1989) 78–83. [CrossRef] (In the text)
 M.G. Silveirinha, N. Engheta, Tunneling of electromagnetic energy through subwavelength channels and bends using εnearzero materials, Phys. Rev. Lett. 97 (2006) 157403. [CrossRef] [PubMed] (In the text)
 A. Alù, M.G. Silveirinha, A. Salandrino, N. Engheta, Epsilon nearzero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern, Phys. Rev. B 75 (2007) 155410. [CrossRef] (In the text)
 R. Liu, Q. Cheng, T. Hand, J.J. Mock, T.J. Cui, S.A. Cummer, D.R. Smith, Experimental demonstration of electromagnetic tunneling through an epsilonnearzero metamaterial at microwave frequencies, Phys. Rev. Lett. 100 (2008) 023903. [CrossRef] [PubMed] (In the text)
 A. Coves, B. Gimeno, J. Gil, M.V. Andrés, A.A. San Blas, V.E. Boria, Fullwave analysis of dielectric frequencyselective surfaces using a vectorial modal method, IEEE Trans. Antennas Propagat. 52 (2004) 2091–2099. [CrossRef] (In the text)
 A. Sommerfeld, Partial differential equations in physics, Academic Press, New York, 1949.
 C.T. Tai, Dyadic Green’s functions in electromagnetic theory. 2nd ed., IEEE Press, New York, 1994.
Appendix A
The total electric field induced on the nongrating structure due to the incident plane wave (1) (A1)is given by (A2)where (A3)with (A4)
The Green’s function G induced on the nongrating structure due to the excitation by a 2D infinite along the yaxis line source, located at an arbitrary point (x′, z′) inside the slab, is analytically expressed as a Fourier integral by applying Sommerfeld’s method [34, 35] and using the techniques of [22] and [23], as follows(A9)where the kernel μ is given by(A10)with (A11)
Appendix B
The following auxiliary functions are utilized for the formulation of the method’s linear systems (11) and (12) (B1)
Cite this article as: Tsitsas NL: Efficient integral equation modeling of scattering by a gradient dielectric metasurface. EPJ Appl. Metamat. 2017, 4, 3.
All Figures
Figure 1.
The crosssection of the gradient dielectric metasurface under consideration, composed of a rectangular periodic grating with period Λ, thickness w and permittivity ε_{g}. A unitamplitude plane wave impinges on the metasurface at an angle θ. 

In the text 
Figure 2.
Convergence patterns of the 0order reflection coefficient r_{0} with respect to the required number N_{t} of expansion functions in (7) for ε_{c} = ε_{s} = 1, 2d = w = 0.7 cm, ε_{d} = 2.59, θ = 60°, Λ = 1.65 cm, with ε_{g} = 0.5, 0.3, and 0.1, and (a) N = 1 groove with a_{1} = 0 and s_{1} = 0.5 Λ, and (b) N = 2 grooves with a_{1} = 0, a_{2} = 0.5 Λ, and s_{1} = 0.1 Λ, s_{2} = 0.3 Λ. 

In the text 
Figure 3.
Amplitude r_{0} of the 0order reflection coefficient versus the operating frequency f for ε_{c} = ε_{s} = 1, 2d = w = 0.7 cm, ε_{d} = 2.59 (Plexiglas), θ = 90°, Λ = 1.65 cm, with ε_{g} = 2.05 (Teflon), 0.5 (Metamaterial 1: MTM1), 0.1 (Metamaterial 2: MTM2), and (a) N = 1 groove with a_{1} = 0 and s_{1} = 0.5 Λ, and (b) N = 3 grooves with a_{1} = 0, a_{2} = 0.3 Λ, a_{3} = 0.6 Λ, and s_{1} = 0.2 Λ, s_{2} = 0.1 Λ, s_{3} = 0.2 Λ. 

In the text 
Figure 4.
(a) Gradient dielectric metasurface (infinite Λperiodic array of dielectric rectangular cylinders) with ε_{c} = ε_{s} = ε_{d} = 1, d = w = 180 nm, ε_{g} = 17.14 (Silicon), λ = 532 nm, Λ = 500 nm, N = 1, a_{1} = 0, and s_{1} = 0.5 Λ. (b) Amplitudes r_{0} and r_{–1} of the 0 and −1order reflection coefficients versus the angle of incidence θ for the metasurface described in (a). 

In the text 
Figure 5.
(a) Metasurface with ε_{c} = ε_{s} = 1, d = w = 180 nm, ε_{d} = 17.14 (Silicon), λ = Λ = 532 nm, N = 1, a_{1} = 0, s_{1} = 0.5 Λ, and ε_{g} = 2.05 (Teflon). (b) Amplitudes r_{0} and r_{−1} of the 0 and −1order reflection coefficients versus the angle θ for the metasurface of (a), (c) same as (b), but for ε_{g} = 0.7. 

In the text 
Figure 6.
(a) Geometry of the gradient dielectric metasurface (ε_{c} = ε_{s} = ε_{d} = 1, Λ = 800 nm, d = Λ/2, w = 2d, N = 1, a_{1} = 0, s_{1} = 0.5 Λ, and θ = 60°). (b)–(e) Reflection and transmission coefficients r_{0}, r_{−1}, t_{0}, and t_{−1} versus the wavelength λ for (b) ε_{g} = 0.2, (c) ε_{g} = 0.3, (d) ε_{g} = 0.5, and (e) ε_{g} = 0.7. 

In the text 