Issue
EPJ Appl. Metamat.
Volume 4, 2017
Artificial materials for advanced applications in electromagnetics and mechanics
Article Number 3
Number of page(s) 13
DOI https://doi.org/10.1051/epjam/2016014
Published online 13 January 2017

© N.L. Tsitsas, Published by EDP Sciences, 2017

Licence Creative CommonsThis 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 beam-to-surface-wave couplers, distributed feedback amplifiers and lasers, multiplexers, and demultiplexers [2], and resonant filters utilized in laser-cavity mode selectors [3] and security and anti-counterfeit 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 [59]. 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 near-field optical enhancement [12]. More recently, it has been elaborated that dielectric materials with low-loss electromagnetic responses, realized by using completely transparent and high-refractive-index dielectric building blocks, may offer the possibility of subdiffraction confinement and guiding of light without metals [13]. To this end, all-dielectric 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 differential-equation techniques and are the rigorous coupled-wave analysis [17] and the modal method [18]. These methods are in certain cases computationally intensive, require time-consuming root searches, large numerical matrices, and multiple time-intensive 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 semi-analytical 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 non-grating 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 narrow-band 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 epsilon-near-zero (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 sn and distances an from z = 0 (n = 1, …, N). The fact that the number N of the rectangles per unit cell and their respective characteristics an and sn may be arbitrarily chosen offers additional degrees of freedom in designing the metasurface with respect to desired operational characteristics. The semi-infinite 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 $ \widehat{\mathbf{y}}$.

thumbnail Figure 1.

The cross-section of the gradient dielectric metasurface under consideration, composed of a rectangular periodic grating with period Λ, thickness w and permittivity εg. A unit-amplitude plane wave impinges on the metasurface at an angle θ.

A TE-polarized 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) E inc ( r ) = Ψ inc ( r ) y ̂ = exp ( - j k inc r ) y ̂ , $$ {\mathbf{E}}^{\mathrm{inc}}\left(\mathbf{r}\right)={\mathrm{\Psi }}^{\mathrm{inc}}\left(\mathbf{r}\right)\widehat{\mathbf{y}}=\mathrm{exp}\left(-j{\mathbf{k}}^{\mathrm{inc}}\cdot \mathbf{r}\right)\widehat{\mathbf{y}}, $$(1)where kinc is the incident wave vector k inc = - k ε c sin θ   x ̂ + k ε c cos θ   z ̂ , $$ {\mathbf{k}}^{\mathrm{inc}}=-k\sqrt{{\epsilon }_c}\mathrm{sin}\theta \enspace \widehat{\mathbf{x}}+k\sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace \widehat{\mathbf{z}}, $$and k = ω/c the free-space wave number.

The electric field induced in every region is of the form E = Ψ ( x , z ) y ̂ , $$ \mathbf{E}=\mathrm{\Psi }\left(x,z\right)\widehat{\mathbf{y}}, $$(2)and the problem is reduced to the determination of the unknown scalar electric field factor Ψ. In the absence of external sources, time-harmonic 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] Ψ ( x , z ) = Ψ 0 ( x , z ) + k 2 ( ε g - ε d ) S g G ( x , z ; x , z ) Ψ ( x , z ) d x d z ,   ( x , z ) R 2 , $$ \mathrm{\Psi }\left(x,z\right)={\mathrm{\Psi }}_0\left(x,z\right)+{k}^2\left({\epsilon }_g-{\epsilon }_d\right)\underset{{S}_g}{\iint }G\left(x,z;x\mathrm{\prime},z\mathrm{\prime}\right)\mathrm{\Psi }\left(x\mathrm{\prime},z\mathrm{\prime}\right)\mathrm{d}x\mathrm{\prime}\mathrm{d}z\mathrm{\prime},\enspace \hspace{1em}\left(x,z\right)\in {\mathbb{R}}^2, $$(3)where Ψ0 is the electric field induced on the non-grating structure due to the plane incident wave (1), G the Green’s function of the non-grating structure with line-source excitation inside the dielectric slab (see Appendix A), and Sg is the total transverse cross-section of the rectangles.

According to the Floquet-Bloch theorem, the electric field factor satisfies the Bloch property [24] (pseudoperiodicity condition) Ψ ( x , z + m Λ ) = exp ( - jk   ε c cos θ   m Λ ) Ψ ( x , z ) ,   m Z , $$ \mathrm{\Psi }\left(x,z+m\mathrm{\Lambda }\right)=\mathrm{exp}\left(-{jk}\enspace \sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace m\mathrm{\Lambda }\right)\mathrm{\Psi }\left(x,z\right),\enspace \hspace{1em}m\in \mathbb{Z}, $$which implies that Ψ ( x , z ) = exp ( - jk   ε c cos θ   z ) u ( x , z ) , $$ \mathrm{\Psi }\left(x,z\right)=\mathrm{exp}\left(-{jk}\enspace \sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace z\right)u\left(x,z\right), $$(4)with u(x, z) a Λ-periodic function of z.

Moreover, by representing the total grating’s cross section Sg as countable union of the cross sections of the grating’s rth unit cell S r = n = 1 N [ d - w , d ] × [ a n + r Λ , a n + s n + r Λ ] ,   r Z , $$ {S}_r=\bigcup_{n=1}^N[d-w,d]\times [{a}_n+r\mathrm{\Lambda },{a}_n+{s}_n+r\mathrm{\Lambda }],\hspace{1em}\enspace r\in \mathbb{Z}, $$and employing the derived Fourier integral (A9) of the Green’s function as well as (4), we reformulate the integral representation (3) as u ( x , z ) = Ψ 0 ( x , z ) exp ( jk ε c cos θ   z ) + k 2 ( ε g - ε d ) 4 π r = - + S r { - + exp [ - ( z - z ) ] μ ( λ , x , x ) d λ } u ( x , z ) exp [ jk ε c cos θ ( z - z ) ]   d x d z , $$ \begin{array}{l}u(x,z)={\mathrm{\Psi }}_0(x,z)\mathrm{exp}({jk}\sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace z)\\ +\frac{{k}^2({\epsilon }_g-{\epsilon }_d)}{4\pi }\sum_{r=-\mathrm{\infty }}^{+\mathrm{\infty }}\underset{{S}_r}{\iint }\left\{\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\mathrm{exp}[-{j\lambda }(z-z\mathrm{\prime})]\mu (\lambda,x,x\mathrm{\prime})\mathrm{d}\lambda \right\}u(x\mathrm{\prime},z\mathrm{\prime})\mathrm{exp}[{jk}\sqrt{{\epsilon }_c}\mathrm{cos}\theta (z-z\mathrm{\prime})]\enspace \mathrm{d}x\mathrm{\prime}\mathrm{d}z\mathrm{\prime}\end{array}, $$(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 S0 (the grating’s cross section on [0, Λ]). Then, by using the Poisson’s summation formula [25] for the Dirac function r = - + exp [ jr Λ ( λ - k ε c cos θ ) ] = 2 π Λ p = - + δ ( λ - k ε c cos θ - 2 π p Λ ) , $$ \sum_{r=-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{exp}\left[{jr}\mathrm{\Lambda }(\lambda -k\sqrt{{\epsilon }_c}\mathrm{cos}\theta )\right]=\frac{2\pi }{\mathrm{\Lambda }}\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}\delta \left(\lambda -k\sqrt{{\epsilon }_c}\mathrm{cos}\theta -\frac{2{\pi p}}{\mathrm{\Lambda }}\right), $$and taking into account the basic property-definition of the Dirac function, we rewrite (5) as u ( x , z ) = Ψ 0 ( x , z ) exp ( jk ε c cos θ   z ) + k 2 ( ε g - ε d ) 2 Λ p = - + exp ( - j 2 π p Λ z ) [ S 0 u ( x , ζ ) exp ( j 2 π p Λ ζ ) μ ( k ε c cos θ + 2 π p Λ , x , x ) d x d ζ ] . $$ \begin{array}{l}u(x,z)={\mathrm{\Psi }}_0(x,z)\mathrm{exp}({jk}\sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace z)+\frac{{k}^2({\epsilon }_g-{\epsilon }_d)}{2\mathrm{\Lambda }}\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{exp}\left(-j\frac{2{\pi p}}{\mathrm{\Lambda }}z\right)\left[\underset{{S}_0}{\iint }u(x\mathrm{\prime},\zeta \mathrm{\prime})\mathrm{exp}\left(j\frac{2{\pi p}}{\mathrm{\Lambda }}\zeta \mathrm{\prime}\right)\mu \left(k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\frac{2{\pi p}}{\mathrm{\Lambda }},x,x\mathrm{\prime}\right)\mathrm{d}x\mathrm{\prime}\mathrm{d}\zeta \mathrm{\prime}\right]\\ \end{array}. $$(6)

3 Solution of the integral equation by an entire-domain 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 S0 is expanded in the Fourier series with respect to z u ( x , z ) = n = - + φ n ( x ) exp [ - j ( 2 π n Λ ) z ] ,   ( x , z ) S 0 , $$ u\left(x,z\right)=\sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{\phi }_n(x)\mathrm{exp}\left[-j\left(\frac{2{\pi n}}{\mathrm{\Lambda }}\right)z\right],\enspace \hspace{1em}\left(x,z\right)\in {S}_0, $$(7)where the Fourier coefficients φn are the transverse components of the spatial harmonics [26], satisfying here the ordinary differential equation φ n ( x ) - g 3 , n 2 φ n ( x ) = 0 ,   d - w x d , $$ {\phi }_n^{\Prime }(x)-{g}_{3,n}^2{\phi }_n(x)=0,\enspace \hspace{1em}d-w\le x\le d, $$(8)where g 3 , n = g 3 ( k ε c cos θ + 2 π n Λ ) = [ ( k ε c cos θ + 2 π n Λ ) 2 - k 2 ε g ] 1 / 2 . $$ {g}_{3,n}={g}_3\left(k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\frac{2{\pi n}}{\mathrm{\Lambda }}\right)={\left[{\left(k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\frac{2{\pi n}}{\mathrm{\Lambda }}\right)}^2-{k}^2{\epsilon }_g\right]}^{1/2}. $$

The functions φn are then expressed as linear combinations of the fundamental solutions of (8) (representing waves traveling parallel to the $ \mp \widehat{\mathbf{x}}$ directions) φ n ( x ) = c n + exp { g 3 , n [ x - d + ( w 2 ) ] } + c n - exp { - g 3 , n [ x - d + ( w 2 ) ] } , $$ {\phi }_n(x)={c}_n^{+}\mathrm{exp}\left\{{g}_{3,n}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}+{c}_n^{-}\mathrm{exp}\left\{-{g}_{3,n}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}, $$(9)with $ {c}_n^{\pm }$ 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 u ( x , z ) = Ψ 0 ( x , z ) exp ( jk ε c cos θ   z ) + k 2 ( ε g - ε d ) 2 p = - + n = - + [ J p - n exp ( - j 2 π p Λ z ) ( c n + Q np + ( x ) + c n - Q np - ( x ) ) ] ,   ( x , z ) R 2   , $$ \begin{array}{l}u\left(x,z\right)={\mathrm{\Psi }}_0\left(x,z\right)\mathrm{exp}({jk}\sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace z)+\frac{{k}^2({\epsilon }_g-{\epsilon }_d)}{2}\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}\sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\left[{J}_{p-n}\mathrm{exp}\left(-j\frac{2{\pi p}}{\mathrm{\Lambda }}z\right)\left({c}_n^{+}{Q}_{{np}}^{+}(x)+{c}_n^{-}{Q}_{{np}}^{-}(x)\right)\right],\enspace (x,z)\in {\mathbb{R}}^2\\ \enspace \end{array}, $$(10)where the functions Jpn and $ {Q}_{{np}}^{\pm }$ 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) ∈ Sg, 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)) exp { ± g 3 , m [ x - d + ( w 2 ) ] } exp [ j ( 2 π m Λ ) z ] ,   ( x , z ) S g ,   m Z , $$ \mathrm{exp}\left\{\pm {g}_{3,m}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}\mathrm{exp}\left[j\left(\frac{2{\pi m}}{\mathrm{\Lambda }}\right)z\right],\enspace \left(x,z\right)\in {S}_g,\enspace m\in \mathbb{Z}, $$and carrying out the resulting integrations, we formulate the following algebraic infinite non-homogeneous square linear system of equations with respect to the unknown coefficients $ {c}_n^{\pm }\enspace (n\in \mathbb{Z})$ n = - + J m - n ( c n + K mn ± + + c n - K mn ± - ) = J m   V m ± + k 2 ( ε g - ε d ) n = - + p = - + J p - n J m - p ( c n + Q mnp ± + + c n - Q mnp ± - ) , $$ \sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}{J}_{m-n}\left({c}_n^{+}{K}_{{mn}}^{\pm +}+{c}_n^{-}{K}_{{mn}}^{\pm -}\right)={J}_m\enspace {V}_m^{\pm }+{k}^2\left({\epsilon }_g-{\epsilon }_d\right)\sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}{J}_{p-n}{J}_{m-p}\left({c}_n^{+}{Q}_{{mnp}}^{\pm +}+{c}_n^{-}{Q}_{{mnp}}^{\pm -}\right), $$(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 Nt, in order to reduce the infinite system to the (4Nt + 2) × (4Nt + 2) linear system [ A + + A + - A - + A - - ] [ c + c - ] = [ b + b - ] , $$ \left[\begin{array}{ll}{\mathbf{A}}^{++}& {\mathbf{A}}^{+-}\\ {\mathbf{A}}^{-+}& {\mathbf{A}}^{--}\end{array}\right]\left[\begin{array}{l}{\mathbf{c}}^{+}\\ {\mathbf{c}}^{-}\end{array}\right]=\left[\begin{array}{l}{\mathbf{b}}^{+}\\ {\mathbf{b}}^{-}\end{array}\right], $$(12)where the (2Nt + 1) × (2Nt + 1) matrices $ {\mathbf{A}}^{\pm \pm }$ are given by ( A ± ± ) mn = - [ J m - n + p = - + k 2 ( ε g - ε d ) g 3 , n 2 - g 1 , p 2 J p - n J m - p ] K mn ± ± + p = - + k 2 ( ε g - ε d ) g 3 , n 2 - g 1 , p 2 J p - n J m - p R mnp ± ± $$ ({\mathbf{A}}^{\pm \pm }{)}_{{mn}}=-\left[{J}_{m-n}+\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{{k}^2\left({\epsilon }_g-{\epsilon }_d\right)}{{g}_{3,n}^2-{g}_{1,p}^2}{J}_{p-n}{J}_{m-p}\right]{K}_{{mn}}^{\pm \pm }+\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{{k}^2\left({\epsilon }_g-{\epsilon }_d\right)}{{g}_{3,n}^2-{g}_{1,p}^2}{J}_{p-n}{J}_{m-p}{R}_{{mnp}}^{\pm \pm } $$(13)

($ {R}_{{mnp}}^{\pm \pm }$ are defined in the Appendix B), c± are 2Nt + 1 column vectors of the coefficients $ {c}_n^{\pm }$, b± the 2Nt + 1 column vectors ( b ± ) m = - J m exp [ ± g 3 , m ( w 2 ) ] { B exp ( jk ε d sin θ 1 d ) 1 - exp [ - ( ± g 3 , m + jk ε d sin θ 1 ) w ] ( ± g 3 , m + jk ε d sin θ 1 ) + C exp ( - jk ε d sin θ 1 d ) 1 - exp [ - ( ± g 3 , m - jk ε d sin θ 1 ) w ] ( ± g 3 , m - jk ε d sin θ 1 ) } , $$ \begin{array}{c}({\mathbf{b}}_{}^{\pm }{)}_m=-{J}_m\mathrm{exp}\left[\pm {g}_{3,m}\left(\frac{w}{2}\right)\right]\left\{B\mathrm{exp}({jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\frac{1-\mathrm{exp}[-(\pm {g}_{3,m}+{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1)w]}{(\pm {g}_{3,m}+{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1)}+C\mathrm{exp}(-{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\frac{1-\mathrm{exp}[-(\pm {g}_{3,m}-{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1)w]}{(\pm {g}_{3,m}-{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1)}\right\}\end{array}, $$(14)while θ1, B, and C are defined by (A4), (A7), and (A8), and g i , n = g i ( k ε c cos θ + 2 π n Λ ) ( i = 0,1 , 2 ) . $$ {g}_{i,n}={g}_i\left(k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\frac{2{\pi n}}{\mathrm{\Lambda }}\right)\hspace{1em}\left(i=\mathrm{0,1},2\right). $$

The required truncation order Nt is determined by applying a convergence control to the solutions for increasing Nt. It has been found that small values of Nt 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 2Nt + 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 non-grating structure.

4 Reflected and transmitted fields

After having determined the coefficients c± by solving the non-homogeneous linear system (12), we express the reflected Ψr and transmitted Ψt electric fields by means of the basic representation (10), as follows Ψ r ( x , z ) = p = - + r p exp [ - j ( k x , p r x + k z , p z ) ] ,   x > d ,   z R , $$ {\mathrm{\Psi }}^r(x,z)=\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}{r}_p\mathrm{exp}[-j({k}_{x,p}^rx+{k}_{z,p}z)],\enspace x>d,\enspace z\in \mathbb{R}, $$(15) Ψ t ( x , z ) = p = - + t p exp [ - j ( k x , p t x + k z , p z ) ] ,   x < - d ,   z R , $$ {\mathrm{\Psi }}^t\left(x,z\right)=\sum_{p=-\mathrm{\infty }}^{+\mathrm{\infty }}{t}_p\mathrm{exp}\left[-j\left({k}_{x,p}^tx+{k}_{z,p}z\right)\right],\enspace x<-d,\enspace z\in \mathbb{R}, $$(16)where the wavevectors’ components are expressed by k z , p = k ε c cos θ + ( 2 π p Λ ) , $$ {k}_{z,p}=k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\left(\frac{2{\pi p}}{\mathrm{\Lambda }}\right), $$ k x , p r = - j { [ k ε c cos θ + ( 2 π p Λ ) ] 2 - k 2 ε c } 1 / 2 , $$ {k}_{x,p}^r=-j{\left\{{\left[k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\left(\frac{2{\pi p}}{\mathrm{\Lambda }}\right)\right]}^2-{k}^2{\epsilon }_c\right\}}^{1/2}, $$(17) k x , p t = j { [ k ε c cos θ + ( 2 π p Λ ) ] 2 - k 2 ε s } 1 / 2 , $$ {k}_{x,p}^t=j{\left\{{\left[k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\left(\frac{2{\pi p}}{\mathrm{\Lambda }}\right)\right]}^2-{k}^2{\epsilon }_s\right\}}^{1/2}, $$while the complex reflection and transmission coefficients are, respectively, given by r p = δ p 0 R + [ k 2 ( ε g - ε d ) 2 ] exp ( g 0 , p d ) n = - + [ J p - n ( c n + ρ np + + c n - ρ np - ) ] ,   p Z , $$ {r}_p={\delta }_{p0}R+\left[\frac{{k}^2\left({\epsilon }_g-{\epsilon }_d\right)}{2}\right]\mathrm{exp}\left({g}_{0,p}d\right)\sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\left[{J}_{p-n}\left({c}_n^{+}{\rho }_{{np}}^{+}+{c}_n^{-}{\rho }_{{np}}^{-}\right)\right],\hspace{1em}\enspace p\in \mathbb{Z}, $$(18) t p = δ p 0 T + [ k 2 ( ε g - ε d ) 2 ] exp ( g 2 , p d ) n = - + [ J p - n ( c n + τ np + + c n - τ np - ) ] ,   p Z . $$ {t}_p={\delta }_{p0}T+\left[\frac{{k}^2\left({\epsilon }_g-{\epsilon }_d\right)}{2}\right]\mathrm{exp}\left({g}_{2,p}d\right)\sum_{n=-\mathrm{\infty }}^{+\mathrm{\infty }}\left[{J}_{p-n}\left({c}_n^{+}{\tau }_{{np}}^{+}+{c}_n^{-}{\tau }_{{np}}^{-}\right)\right],\enspace \hspace{1em}p\in \mathbb{Z}. $$(19)

Coefficients R and T are given by (A5) and (A6), δp0 is the Kronecker symbol (:δ00 = 1 and δp0 = 0|p ≠ 0), while ρ np ± = [ ( g 1 , p 2 - g 3 , n 2 ) Λ ( g 0 , p , g 2 , p , g 1 , p , 2 d ) ] - 1 { exp ( ± g 3 , n w 2 ) [ P o ( g 2 , p , g 1 , p , 2 d ) P e ( g 3 , n , g 1 , p , 2 d ) + P e ( g 2 , p , g 1 , p , 2 d ) P o ( g 3 , n , g 1 , p , 2 d ) ] - exp ( g 3 , n w 2 ) [ P o ( g 2 , p , g 1 , p , 2 d ) P e ( g 3 , n , g 1 , p , 2 d - 2 w ) + P e ( g 2 , p , g 1 , p , 2 d ) P o ( g 3 , n , g 1 , p , 2 d - 2 w ) ] } , $$ \begin{array}{l}{\rho }_{{np}}^{\pm }=[({g}_{1,p}^2-{g}_{3,n}^2)\mathrm{\Lambda }({g}_{0,p},{g}_{2,p},{g}_{1,p},2d){]}^{-1}\\ \left\{\mathrm{exp}\left(\pm \frac{{g}_{3,n}w}{2}\right)\right.\left[{P}_o\left({g}_{2,p},{g}_{1,p},2d\right){P}_e\left(\mp {g}_{3,n},{g}_{1,p},2d\right)+{P}_e\left({g}_{2,p},{g}_{1,p},2d\right){P}_o\left(\mp {g}_{3,n},{g}_{1,p},2d\right)\right]\\ -\mathrm{exp}\left(\mp \frac{{g}_{3,n}w}{2}\right)[{P}_o\left({g}_{2,p},{g}_{1,p},2d\right){P}_e\left(\mp {g}_{3,n},{g}_{1,p},2d-2w\right)\left.+{P}_e({g}_{2,p},{g}_{1,p},2d){P}_o(\mp {g}_{3,n},{g}_{1,p},2d-2w)]\right\},\end{array} $$(20) τ np ± = [ ( g 1 , p 2 - g 3 , n 2 ) Λ ( g 0 , p , g 2 , p , g 1 , p , 2 d ) ] - 1 { exp ( ± g 3 , n w 2 ) [ P o ( g 0 , p , g 1 , p , 2 d ) P e ( g 3 , n , g 1 , p , 2 d ) - P e ( g 0 , p , g 1 , p , 2 d ) P o ( g 3 , n , g 1 , p , 2 d ) ] - exp ( g 3 , n w 2 ) [ P o ( g 0 , p , g 1 , p , 2 d ) P e ( g 3 , n , g 1 , p , 2 d - 2 w ) - P e ( g 0 , p , g 1 , p , 2 d ) P o ( g 3 , n , g 1 , p , 2 d - 2 w ) ] } . $$ \begin{array}{l}{\tau }_{{np}}^{\pm }=[({g}_{1,p}^2-{g}_{3,n}^2)\mathrm{\Lambda }({g}_{0,p},{g}_{2,p},{g}_{1,p},2d){]}^{-1}\\ \left\{\mathrm{exp}\left(\pm \frac{{g}_{3,n}w}{2}\right)\right.[{P}_o({g}_{0,p},{g}_{1,p},2d){P}_e(\mp {g}_{3,n},{g}_{1,p},2d)-{P}_e({g}_{0,p},{g}_{1,p},2d){P}_o(\mp {g}_{3,n},{g}_{1,p},2d)]\\ -\mathrm{exp}\left(\mp \frac{{g}_{3,n}w}{2}\right)[{P}_o\left({g}_{0,p},{g}_{1,p},2d\right){P}_e\left(\mp {g}_{3,n},{g}_{1,p},2d-2w\right)\left.-{P}_e({g}_{0,p},{g}_{1,p},2d){P}_o(\mp {g}_{3,n},{g}_{1,p},2d-2w)]\right\}.\end{array} $$(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 p-reflected or p-transmitted diffracted order. The diffracted orders, which propagate along the x-axis, 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 0-order fields are always propagating; notice also the contribution of coefficients R and T in (18) and (19) in the 0-reflected and transmitted orders stemming from the fields induced on the non-grating structure. For orders p ≠ 0, we define the thresholds p r ± = ± ( Λ λ ) ε c ( 1 cos θ ) , $$ {p}_r^{\pm }=\pm \left(\frac{\mathrm{\Lambda }}{\lambda }\right)\sqrt{{\epsilon }_c}\left(1\mp \mathrm{cos}\theta \right), $$ p t ± = ± ( Λ λ ) ( ε s ε c cos θ ) , $$ {p}_t^{\pm }=\pm \left(\frac{\mathrm{\Lambda }}{\lambda }\right)\left(\sqrt{{\epsilon }_s}\mp \sqrt{{\epsilon }_c}\mathrm{cos}\theta \right), $$(22)for which holds:

  1. for $ p<{p}_r^{+}$ or $ p>{p}_r^{-}$ and $ p<{p}_t^{+}$ or $ p>{p}_t^{-}$ the wavevectors’ components $ {k}_{x,p}^r$ and $ {k}_{x,p}^t$ are real and the p-order reflected and transmitted fields are propagating along the x-axis.

  2. for $ p>{p}_r^{+}$ or $ p<{p}_r^{-}$ and $ p>{p}_t^{+}$ or $ p<{p}_t^{-}$ the wavevectors’ components $ {k}_{x,p}^r$ and $ {k}_{x,p}^t$ become purely imaginary and the p-order fields are evanescent.

Hence, $ {p}_r^{\pm }$ and $ {p}_t^{\pm }$ 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 Nt 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 108) and (ii) convergence to the solution with increasing Nt for all the grating and the incident wave parameters [17]. Some representative convergence patterns of the 0-order 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 a1 = 0 and s1 = 0.5 Λ, and (b) N = 2 grooves with a1 = 0, a2 = 0.5 Λ, and s1 = 0.1 Λ, s2 = 0.3 Λ. It is observed that Nt = 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 Nt = 2 and Nt = 3), while Nt = 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).

thumbnail Figure 2.

Convergence patterns of the 0-order reflection coefficient |r0| with respect to the required number Nt 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 a1 = 0 and s1 = 0.5 Λ, and (b) N = 2 grooves with a1 = 0, a2 = 0.5 Λ, and s1 = 0.1 Λ, s2 = 0.3 Λ.

For all the results presented hereafter, it has been checked that the choice Nt = 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, an, and sn.

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 |r0| of the 0-reflected 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 a1 = 0, and s1 = 0.5 Λ, while those of Figure 3b to a grating with three grooves per unit cell, i.e. N = 3, with a1 = 0, s1 = 0.2 Λ, a2 = 0.3 Λ, s2 = 0.1 Λ, a3 = 0.6 Λ, and s3 = 0.2 Λ.

thumbnail Figure 3.

Amplitude |r0| of the 0-order 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 a1 = 0 and s1 = 0.5 Λ, and (b) N = 3 grooves with a1 = 0, a2 = 0.3 Λ, a3 = 0.6 Λ, and s1 = 0.2 Λ, s2 = 0.1 Λ, s3 = 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 |r0| 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 (|r0| < 0.2 in the interval from 17 to 20 GHz). Metamaterials with permittivities values like εg = 0.1 belong to the class of epsilon-near-zero (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 low-loss 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 |r0| 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 an and sn 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 so-called 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 $ {p}_r^{+}<1,\enspace {p}_t^{+}<1$ and $ -2<{p}_r^{-}<-1,\enspace -2<{p}_t^{-}<-1$, which when εs = εc and thus $ {p}_r^{\pm }={p}_t^{\pm }$, lead to max { ε c ( 1 - cos θ ) , ε c 1 + cos θ 2 } < λ Λ < ε c ( 1 + cos θ ) . $$ \mathrm{max}\left\{\sqrt{{\epsilon }_c}(1-\mathrm{cos}\theta ),\sqrt{{\epsilon }_c}\frac{1+\mathrm{cos}\theta }{2}\right\}<\frac{\lambda }{\mathrm{\Lambda }}<\sqrt{{\epsilon }_c}\left(1+\mathrm{cos}\theta \right). $$(23)

It has been recently shown experimentally that anomalous reflection effects can be generated by high-contrast metastructures (HCM) which consist of periodic structures of a single high-refractive-index layer surrounded by a low-index 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 0-order of reflection (following Snell’s law) and simultaneously enhance the −1-order 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 low-index material. Figure 4b depicts the amplitudes of the 0- and −1-order 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, a1 = 0, and s1 = 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 |r0| in the range of θ which is mainly of interest (excluding the near-grazing incidence angles). Particularly, the locally maximum distances between |r–1| and |r0| 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.

thumbnail 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, a1 = 0, and s1 = 0.5 Λ. (b) Amplitudes |r0| and |r–1| of the 0- and −1-order 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 |r0| (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 |r0|. 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).

thumbnail Figure 5.

(a) Metasurface with εc = εs = 1, d = w = 180 nm, εd = 17.14 (Silicon), λ = Λ = 532 nm, N = 1, a1 = 0, s1 = 0.5 Λ, and εg = 2.05 (Teflon). (b) Amplitudes |r0| and |r−1| of the 0- and −1-order 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 −1-reflected order compared with the one of the 0-order 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 6b6e the amplitudes |r0|, |r−1|, |t0|, and |t−1| of the 0- and −1-order 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 −1-transmission order compared to the 0-transmission order is attained, while simultaneously both the 0- and −1-reflection 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 −1-transmission order is still higher than the 0-transmission order but their difference starts to reduce (see Figure 6c). The situation is reversed in Figure 6d, corresponding to εg = 0.5, where the 0-transmission order becomes larger than the −1-order. 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 −1-orders, suppression of the transmission of the −1-order and nearly one transmission coefficient of the 0-order (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.

thumbnail Figure 6.

(a) Geometry of the gradient dielectric metasurface (εc = εs = εd = 1, Λ = 800 nm, d = Λ/2, w = 2d, N = 1, a1 = 0, s1 = 0.5 Λ, and θ = 60°). (b)–(e) Reflection and transmission coefficients |r0|, |r−1|, |t0|, 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 semi-analytical 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 all-dielectric 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

  1. S. Tibuleac, R. Magnusson, T.A. Maldonado, P.P. Young, T.R. Holzheimer, Dielectric frequency-selective structures incorporating waveguide gratings, IEEE Trans. Microwave Theory Tech. 48 (2000) 553–561. [CrossRef] [Google Scholar]
  2. S.S. Wang, R. Magnusson, Design of waveguide grating filters with symmetrical line shapes and low sidebands, Opt. Lett. 19 (1994) 919–921. [CrossRef] [Google Scholar]
  3. S.M. Norton, T. Erdogan, G.M. Morris, Coupled mode theory of resonant-grating filters, J. Opt. Soc. Am. A 14 (1997) 629–639. [CrossRef] [Google Scholar]
  4. J. Turunen, F. Wyrowski (Eds.), Diffractive optics for industrial and commercial applications, Chap. 12, Akademie, Berlin, 1997. [Google Scholar]
  5. 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] [Google Scholar]
  6. 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] [Google Scholar]
  7. 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] [Google Scholar]
  8. F. Monticone, N.M. Estakhri, A. Alù, Full control of nanoscale optical transmission with a composite metascreen, Phys. Rev. Lett. 110 (2013) 203903. [CrossRef] [Google Scholar]
  9. 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. [Google Scholar]
  10. 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] [Google Scholar]
  11. Y. Cui, et al., Plasmonic and metamaterial structures as electromagnetic absorbers, Laser Photon. Rev. 8 (2014) 495–520. [CrossRef] [Google Scholar]
  12. A.A. Darweesh, S.J. Bauman, J.B. Herzog, Improved optical enhancement using double-width plasmonic gratings with nanogaps, Photon. Research 4 (2016) 173–180. [CrossRef] [Google Scholar]
  13. S. Jahani, Z. Jacob, All-dielectric metamaterials, Nature Nanotechnol. 11 (2016) 23–36. [Google Scholar]
  14. D. Ohana, U. Levy, Mode conversion based on dielectric metamaterial in silicon, Optics Express 22 (2014) 27617–27631. [CrossRef] [Google Scholar]
  15. L. Zhu, J. Kapraun, J. Ferrara, C.J. Chang-Hasnain, Flexible photonic metastructures for tunable coloration, Optica 2 (2015) 255–258. [CrossRef] [Google Scholar]
  16. S. Miyanaga, T. Akasura, Intensity profiles of outgoing beams from tapered grating couplers, Radio Sci. 17 (1982) 135–143. [CrossRef] [Google Scholar]
  17. M.G. Moharam, E.B. Grann, D.A. Pommet, T.K. Gaylord, Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings, J. Opt. Soc. Am. A 12 (1995) 1068–1076. [CrossRef] [Google Scholar]
  18. R.S. Chu, J.A. Kong, Modal theory of spatially periodic media, IEEE Trans. Microwave Theory Tech. MTT-25 (1977) 18–24. [CrossRef] [Google Scholar]
  19. 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] [Google Scholar]
  20. Y. Nakata, M. Koshiba, Boundary-element analysis of plane-wave diffraction from groove-type dielectric and metallic gratings, J. Opt. Soc. Am. A 7 (1990) 1494–1502. [CrossRef] [Google Scholar]
  21. G. Athanasoulias, N.K. Uzunoglu, An accurate and efficient entire-domain basis Galerkin’s method for the integral equation analysis of integrated rectangular dielectric waveguides, IEEE Trans. Microwave Theory Tech. 43 (1995) 2794–2804. [CrossRef] [Google Scholar]
  22. 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] [Google Scholar]
  23. 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] [Google Scholar]
  24. 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] [Google Scholar]
  25. P.M. Morse, H. Feshbach, Methods of Theoretical Physics, Part I, McGraw-Hill, Tokyo, 1953. [Google Scholar]
  26. R.E. Collin, Field theory of guided waves, IEEE Press, New York, 1991. [Google Scholar]
  27. 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] [Google Scholar]
  28. 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] [Google Scholar]
  29. H.L. Bertoni, L.-H.S. Cheo, T. Tamir, Frequency-selective reflection and transmission by a periodic dielectric layer, IEEE Trans. Antennas Propagat. 37 (1989) 78–83. [CrossRef] [Google Scholar]
  30. M.G. Silveirinha, N. Engheta, Tunneling of electromagnetic energy through subwavelength channels and bends using ε-near-zero materials, Phys. Rev. Lett. 97 (2006) 157403. [CrossRef] [PubMed] [Google Scholar]
  31. A. Alù, M.G. Silveirinha, A. Salandrino, N. Engheta, Epsilon near-zero metamaterials and electromagnetic sources: Tailoring the radiation phase pattern, Phys. Rev. B 75 (2007) 155410. [CrossRef] [Google Scholar]
  32. 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 epsilon-near-zero metamaterial at microwave frequencies, Phys. Rev. Lett. 100 (2008) 023903. [CrossRef] [PubMed] [Google Scholar]
  33. A. Coves, B. Gimeno, J. Gil, M.V. Andrés, A.A. San Blas, V.E. Boria, Full-wave analysis of dielectric frequency-selective surfaces using a vectorial modal method, IEEE Trans. Antennas Propagat. 52 (2004) 2091–2099. [CrossRef] [Google Scholar]
  34. A. Sommerfeld, Partial differential equations in physics, Academic Press, New York, 1949. [Google Scholar]
  35. C.T. Tai, Dyadic Green’s functions in electromagnetic theory. 2nd ed., IEEE Press, New York, 1994. [Google Scholar]

Appendix A

The total electric field induced on the non-grating structure due to the incident plane wave (1) E 0 = Ψ 0 ( x , z ) y ̂ , $$ {\mathbf{E}}_0={\mathrm{\Psi }}_0\left(x,z\right)\widehat{\mathbf{y}}, $$(A1)is given by Ψ 0 ( x , z ) = exp ( - jk ε c cos θ   z ) Ψ 0 ( x ) , $$ {\mathrm{\Psi }}_0\left(x,z\right)=\mathrm{exp}\left(-{jk}\sqrt{{\epsilon }_c}\mathrm{cos}\theta \enspace z\right){\mathrm{\Psi }}_0(x), $$(A2)where Ψ 0 ( x ) = { exp ( jk ε c sin θ   x ) + R exp ( - jk ε c sin θ   x ) , x d B exp ( jk ε d sin θ 1 x ) + C exp ( - jk ε d sin θ 1 x ) , - d x d T exp ( jk ε s sin θ 2 x ) x - d , $$ {\mathrm{\Psi }}_0(x)=\left\{\begin{array}{ll}\mathrm{exp}\left({jk}\sqrt{{\epsilon }_c}\mathrm{sin}\theta \enspace x\right)+R\mathrm{exp}\left(-{jk}\sqrt{{\epsilon }_c}\mathrm{sin}\theta \enspace x\right),& x\ge d\\ B\mathrm{exp}\left({jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1x\right)+C\mathrm{exp}\left(-{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1x\right),& -d\le x\le d\\ T\mathrm{exp}\left({jk}\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2x\right)& x\le -d\end{array}\right., $$(A3)with θ 1 = co s - 1 [ ( ε c ε d ) cos θ ] ,   θ 2 = co s - 1 [ ( ε c ε s ) cos θ ] , $$ {\theta }_1=\mathrm{co}{\mathrm{s}}^{-1}\left[\left(\frac{\sqrt{{\epsilon }_c}}{\sqrt{{\epsilon }_d}}\right)\mathrm{cos}\theta \right],\enspace {\theta }_2=\mathrm{co}{\mathrm{s}}^{-1}\left[\left(\frac{\sqrt{{\epsilon }_c}}{\sqrt{{\epsilon }_s}}\right)\mathrm{cos}\theta \right], $$(A4)

R = exp ( 2 jk ε c sin θ   d ) [ exp ( 2 jk ε d sin θ 1 d ) ( 1 + ε s sin θ 2 ε d sin θ 1 ) ( 1 - ε d sin θ 1 ε c sin θ ) + exp ( - 2 jk ε d sin θ 1 d ) ( 1 - ε s sin θ 2 ε d sin θ 1 ) ( 1 + ε d sin θ 1 ε c sin θ ) ]   [ exp ( 2 jk ε d sin θ 1 d ) ( 1 + ε s sin θ 2 ε d sin θ 1 ) ( 1 + ε d sin θ 1 ε c sin θ ) + exp ( - 2 jk ε d sin θ 1 d ) ( 1 - ε s sin θ 2 ε d sin θ 1 ) ( 1 - ε d sin θ 1 ε c sin θ ) ] - 1 , $$ \begin{array}{l}R=\mathrm{exp}(2{jk}\sqrt{{\epsilon }_c}\mathrm{sin}\theta \enspace d)\left[\begin{array}{l}\mathrm{exp}(2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1+\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1-\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\\ +\mathrm{exp}(-2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1-\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1+\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\end{array}\right]\\ \enspace {\left[\begin{array}{l}\mathrm{exp}(2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1+\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1+\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\\ +\mathrm{exp}(-2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1-\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1-\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\end{array}\right]}^{-1}\end{array}, $$(A5)

T = 4   exp [ jk ( ε c sin θ   + ε s sin θ 2 ) d ]   [ exp ( 2 jk ε d sin θ 1 d ) ( 1 + ε s sin θ 2 ε d sin θ 1 ) ( 1 + ε d sin θ 1 ε c sin θ ) + exp ( - 2 jk ε d sin θ 1 d ) ( 1 - ε s sin θ 2 ε d sin θ 1 ) ( 1 - ε d sin θ 1 ε c sin θ ) ] - 1 $$ T=4\enspace \mathrm{exp}[{jk}(\sqrt{{\epsilon }_c}\mathrm{sin}\theta \enspace +\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2)d]\enspace {\left[\begin{array}{l}\mathrm{exp}(2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1+\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1+\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\\ +\mathrm{exp}(-2{jk}\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1d)\left(1-\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right)\left(1-\frac{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}{\sqrt{{\epsilon }_c}\mathrm{sin}\theta }\right)\end{array}\right]}^{-1} $$(A6)

B = ( T 2 ) exp [ jk ( ε d sin θ 1 - ε s sin θ 2 ) d ] { 1 + [ ε s sin θ 2 ε d sin θ 1 ] } , $$ B=\left(\frac{T}{2}\right)\mathrm{exp}\left[{jk}\left(\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1-\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2\right)d\right]\left\{1+\left[\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right]\right\}, $$(A7)

C = ( T 2 ) exp [ - jk ( ε d sin θ 1 + ε s sin θ 2 ) d ] { 1 - [ ε s sin θ 2 ε d sin θ 1 ] } . $$ C=\left(\frac{T}{2}\right)\mathrm{exp}\left[-{jk}\left(\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1+\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2\right)d\right]\left\{1-\left[\frac{\sqrt{{\epsilon }_s}\mathrm{sin}{\theta }_2}{\sqrt{{\epsilon }_d}\mathrm{sin}{\theta }_1}\right]\right\}. $$(A8)

The Green’s function G induced on the non-grating structure due to the excitation by a 2-D infinite along the y-axis 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 G ( x , z ; x , z ) = 1 4 π - + e - ( z - z ) μ ( λ , x , x ) d λ , $$ G\left(x,z;x\mathrm{\prime},z\mathrm{\prime}\right)=\frac{1}{4\pi }\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}{e}^{-{j\lambda }\left(z-{z}^\mathrm{\prime}\right)}\mu \left(\lambda,x,x\mathrm{\prime}\right)\mathrm{d}\lambda, $$(A9)where the kernel μ is given by μ ( λ , x , x ) = [ Λ ( g 0 , g 2 , g 1 , 2 d ) ] - 1 { P o ( g 2 , g 1 , 2 d ) P e ( g 2 , g 1 , 2 d ) [ cosh ( g 1 x ) P e ( g 2 , g 1 , 2 d ) + sinh ( g 1 x ) P o ( g 2 , g 1 , 2 d ) ] exp [ - g 0 ( x - d ) ] , x d g 1 - 1 P e ( g 0 , g 1 , 2 d ) P e ( g 2 , g 1 , 2 d ) P o ( g 2 , g 1 , 2 d ) P o ( g 0 , g 1 , 2 d ) [ cosh ( g 1 x ) P e ( g 2 , g 1 , 2 d ) + sinh ( g 1 x ) P o ( g 2 , g 1 , 2 d ) ] [ cosh ( g 1 x ) P e ( g 0 , g 1 , 2 d ) - sinh ( g 1 x ) P o ( g 0 , g 1 , 2 d ) ] , - d x x d g 1 - 1 P e ( g 0 , g 1 , 2 d ) P e ( g 2 , g 1 , 2 d ) P o ( g 2 , g 1 , 2 d ) P o ( g 0 , g 1 , 2 d ) [ cosh ( g 1 x ) P e ( g 0 , g 1 , 2 d ) - sinh ( g 1 x ) P o ( g 0 , g 1 , 2 d ) ] [ cosh ( g 1 x ) P e ( g 2 , g 1 , 2 d ) + sinh ( g 1 x ) P o ( g 2 , g 1 , 2 d ) ] , - d x x d P o ( g 0 , g 1 , 2 d ) P e ( g 0 , g 1 , 2 d ) [ cosh ( g 1 x ) P e ( g 0 , g 1 , 2 d ) - sinh ( g 1 x ) P o ( g 0 , g 1 , 2 d ) ] exp [ g 2 ( x + d ) ] , x - d $$ \begin{array}{l}\mu (\lambda,x,x\mathrm{\prime})=[\mathrm{\Lambda }({g}_0,{g}_2,{g}_1,2d){]}^{-1}\\ \left\{\begin{array}{ll}\begin{array}{l}{P}_o\left({g}_2,{g}_1,2d\right){P}_e\left({g}_2,{g}_1,2d\right)\left[\frac{\mathrm{cosh}\left({g}_1x\mathrm{\prime}\right)}{{P}_e\left({g}_2,{g}_1,2d\right)}+\frac{\mathrm{sinh}\left({g}_1x\mathrm{\prime}\right)}{{P}_o\left({g}_2,{g}_1,2d\right)}\right]\mathrm{exp}\left[-{g}_0\left(x-d\right)\right],\\ \end{array}& x\ge d\\ \begin{array}{l}{g}_1^{-1}{P}_e\left({g}_0,{g}_1,2d\right){P}_e\left({g}_2,{g}_1,2d\right){P}_o\left({g}_2,{g}_1,2d\right){P}_o\left({g}_0,{g}_1,2d\right)\\ \left[\frac{\mathrm{cosh}\left({g}_1{x}^\mathrm{\prime}\right)}{{P}_e\left({g}_2,{g}_1,2d\right)}+\frac{\mathrm{sinh}\left({g}_1x\mathrm{\prime}\right)}{{P}_o\left({g}_2,{g}_1,2d\right)}\right]\left[\frac{\mathrm{cosh}\left({g}_1x\right)}{{\mathrm{P}}_e\left({g}_0,{g}_1,2d\right)}-\frac{\mathrm{sinh}\left({g}_1x\right)}{{P}_o\left({g}_0,{g}_1,2d\right)}\right],\\ \end{array}& -d\le {x}^\mathrm{\prime}\le x\le d\\ \begin{array}{l}{g}_1^{-1}{P}_e\left({g}_0,{g}_1,2d\right){P}_e\left({g}_2,{g}_1,2d\right){P}_o\left({g}_2,{g}_1,2d\right){P}_o\left({g}_0,{g}_1,2d\right)\\ \left[\frac{\mathrm{cosh}\left({g}_1x\mathrm{\prime}\right)}{{P}_e\left({g}_0,{g}_1,2d\right)}-\frac{\mathrm{sinh}\left({g}_1x\mathrm{\prime}\right)}{{P}_o\left({g}_0,{g}_1,2d\right)}\right]\left[\frac{\mathrm{cosh}\left({g}_1x\right)}{{P}_e\left({g}_2,{g}_1,2d\right)}+\frac{\mathrm{sinh}\left({g}_1x\right)}{{P}_o\left({g}_2,{g}_1,2d\right)}\right],\\ \end{array}& -d\le x\le {x}^\mathrm{\prime}\le d\\ {P}_o\left({g}_0,{g}_1,2d\right){P}_e\left({g}_0,{g}_1,2d\right)\left[\frac{\mathrm{cosh}\left({g}_1x\mathrm{\prime}\right)}{{P}_e\left({g}_0,{g}_1,2d\right)}-\frac{\mathrm{sinh}\left({g}_1x\mathrm{\prime}\right)}{{P}_o\left({g}_0,{g}_1,2d\right)}\right]\mathrm{exp}\left[{g}_2\left(x+d\right)\right],& x\le -d\end{array}\right.\end{array} $$(A10)with P e ( g i , g 1 , 2 d ) = g i cosh ( g 1 d ) + g 1 sinh ( g 1 d ) ,   P o ( g i , g 1 , 2 d ) = g 1 cosh ( g 1 d ) + g i sinh ( g 1 d ) , $$ {{P}_e\left({g}_i,{g}_1,2d\right)={g}_i\mathrm{cosh}\left({g}_1d\right)+{g}_1\mathrm{sinh}\left({g}_1d\right),\enspace P}_o\left({g}_i,{g}_1,2d\right)={g}_1\mathrm{cosh}\left({g}_1d\right)+{g}_i\mathrm{sinh}\left({g}_1d\right), $$(A11)

Λ ( g 0 , g 2 , g 1 , 2 d ) = [ P e ( g 0 , g 1 , 2 d ) P o ( g 2 , g 1 , 2 d ) + P o ( g 0 , g 1 , 2 d ) P e ( g 2 , g 1 , 2 d ) ] 2 , $$ \mathrm{\Lambda }({g}_0,{g}_2,{g}_1,2d)=\frac{[{P}_e({g}_0,{g}_1,2d){P}_o({g}_2,{g}_1,2d)+{P}_o({g}_0,{g}_1,2d){P}_e({g}_2,{g}_1,2d)]}{2}, $$(A12)

g 0 ( λ ) = ( λ 2 - k 2 ε c ) 1 / 2 ,   g 1 ( λ ) = ( λ 2 - k 2 ε d ) 1 / 2 ,   g 2 ( λ ) = ( λ 2 - k 2 ε s ) 1 / 2 . $$ {g}_0(\lambda )=({\lambda }^2-{k}^2{\epsilon }_c{)}^{1/2},\enspace {g}_1(\lambda )=({\lambda }^2-{k}^2{\epsilon }_d{)}^{1/2},\enspace {g}_2(\lambda )=({\lambda }^2-{k}^2{\epsilon }_s{)}^{1/2}. $$(A13)

Appendix B

The following auxiliary functions are utilized for the formulation of the method’s linear systems (11) and (12) J l =   1 Λ n = 1 N a n a n + s n exp [ j ( 2 π l Λ ) ζ ] d ζ ,   l Z , $$ {J}_{\mathcal{l}}=\enspace \frac{1}{\mathrm{\Lambda }}\sum_{n=1}^N\underset{{a}_n}{\overset{{a}_n+{s}_n}{\int }}\mathrm{exp}\left[j\left(\frac{2\pi \mathcal{l}}{\mathrm{\Lambda }}\right)\zeta \mathrm{\prime}\right]\mathrm{d}\zeta \mathrm{\prime},\enspace \hspace{1em}\mathcal{l}\in \mathbb{Z}, $$(B1)

Q np ± ( x ) = d - w d exp { ± g 3 , n [ x - d + ( w 2 ) ] } μ ( k ε c cos θ + 2 π p Λ , x , x ) d x ,   x R , $$ {Q}_{{np}}^{\pm }(x)=\underset{d-w}{\overset{d}{\int }}\mathrm{exp}\left\{\pm {g}_{3,n}\left[x\mathrm{\prime}-d+\left(\frac{w}{2}\right)\right]\right\}\mu \left(k\sqrt{{\epsilon }_c}\mathrm{cos}\theta +\frac{2{\pi p}}{\mathrm{\Lambda }},x,x\mathrm{\prime}\right)\mathrm{d}x\mathrm{\prime},\enspace \hspace{1em}x\in \mathbb{R}, $$(B2)

K mn ± ± = d - w d exp { ± g 3 , m [ x - d + ( w 2 ) ] } exp { ± g 3 , n [ x - d + ( w 2 ) ] } d x ,   m , n Z , $$ {K}_{{mn}}^{\pm \pm }=\underset{d-w}{\overset{d}{\int }}\mathrm{exp}\left\{\pm {g}_{3,m}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}\mathrm{exp}\left\{\pm {g}_{3,n}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}\mathrm{d}x,\enspace \hspace{1em}m,n\in \mathbb{Z}, $$(B3)

Q mnp ± ± = 1 2 d - w d exp { ± g 3 , m [ x - d + ( w 2 ) ] } Q np ± ( x ) d x ,   m ,   n , p Z , $$ {Q}_{{mnp}}^{\pm \pm }=\frac{1}{2}\underset{d-w}{\overset{d}{\int }}\mathrm{exp}\left\{\pm {g}_{3,m}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}{Q}_{{np}}^{\pm }(x)\mathrm{d}x,\hspace{1em}\enspace m,\enspace n,p\in \mathbb{Z}, $$(B4)

V m ± = d - w d Ψ 0 ( x ) exp { ± g 3 , m [ x - d + ( w 2 ) ] } d x ,   m Z , $$ {V}_m^{\pm }=\underset{d-w}{\overset{d}{\int }}{\mathrm{\Psi }}_0(x)\mathrm{exp}\left\{\pm {g}_{3,m}\left[x-d+\left(\frac{w}{2}\right)\right]\right\}\mathrm{d}x,\enspace \hspace{1em}m\in \mathbb{Z}, $$(B5)

R mnp ± ± = P e ( g 0 , p , g 1 , p , 2 d ) P e ( g 2 , p , g 1 , p , 2 d ) P o ( g 2 , p , g 1 , p , 2 d ) P o ( g 0 , p , g 1 , p , 2 d ) 2 ( g 3 , m 2 - g 1 , p 2 ) g 1 , p Λ ( g 0 , p , g 2 , p , g 1 , p , 2 d ) { e ( ± g 3 , m ± g 3 , n ) w / 2 ( P e ( g 3 , m , g 1 , p , 2 d ) P e ( g 2 , p , g 1 , p , 2 d ) + P o ( g 3 , m , g 1 , p , 2 d ) P o ( g 2 , p , g 1 , p , 2 d ) ) ( P e ( g 3 , n , g 1 , p , 2 d ) P e ( g 0 , p , g 1 , p , 2 d ) - P o ( g 3 , n , g 1 , p , 2 d ) P o ( g 0 , p , g 1 , p , 2 d ) ) - e ( ± g 3 , m g 3 , n ) w / 2 ( P e ( g 3 , m , g 1 , p , 2 d ) P e ( g 0 , p , g 1 , p , 2 d ) - P o ( g 3 , m , g 1 , p , 2 d ) P o ( g 0 , p , g 1 , p , 2 d ) ) ( P e ( g 3 , n , g 1 , p , 2 d - 2 w ) P e ( g 2 , p , g 1 , p , 2 d ) + P o ( g 3 , n , g 1 , p , 2 d - 2 w ) P o ( g 2 , p , g 1 , p , 2 d ) ) - e ( g 3 , m ± g 3 , n ) w / 2 ( P e ( g 3 , m , g 1 , p , 2 d - 2 w ) P e ( g 2 , p , g 1 , p , 2 d ) + P o ( g 3 , m , g 1 , p , 2 d - 2 w ) P o ( g 2 , p , g 1 , p , 2 d ) ) ( P e ( g 3 , n , g 1 , p , 2 d ) P e ( g 0 , p , g 1 , p , 2 d ) - P o ( g 3 , n , g 1 , p , 2 d ) P o ( g 0 , p , g 1 , p , 2 d ) ) + e ( g 3 , m g 3 , n ) w / 2 ( P e ( g 3 , m , g 1 , p , 2 d - 2 w ) P e ( g 0 , p , g 1 , p , 2 d ) - P o ( g 3 , m , g 1 , p , 2 d - 2 w ) P o ( g 0 , p , g 1 , p , 2 d ) )   ( P e ( g 3 , n , g 1 , p , 2 d - 2 w ) P e ( g 2 , p , g 1 , p , 2 d ) + P o ( g 3 , n , g 1 , p , 2 d - 2 w ) P o ( g 2 , p , g 1 , p , 2 d ) ) } $$ \begin{array}{l}{R}_{{mnp}}^{\pm \pm }=\frac{{P}_e({g}_{0,p},{g}_{1,p},2d){P}_e({g}_{2,p},{g}_{1,p},2d){P}_o({g}_{2,p},{g}_{1,p},2d){P}_o({g}_{0,p},{g}_{1,p},2d)}{2({g}_{3,m}^2-{g}_{1,p}^2){g}_{1,p}\mathrm{\Lambda }({g}_{0,p},{g}_{2,p},{g}_{1,p},2d)}\\ \left\{{e}^{(\pm {g}_{3,m}\pm {g}_{3,n})w/2}\left(\frac{{P}_e(\mp {g}_{3,m},{g}_{1,p},2d)}{{P}_e({g}_{2,p},{g}_{1,p},2d)}+\frac{{P}_o(\mp {g}_{3,m},{g}_{1,p},2d)}{{P}_o({g}_{2,p},{g}_{1,p},2d)}\right)\right.\left(\frac{{P}_e(\mp {g}_{3,n},{g}_{1,p},2d)}{{P}_e({g}_{0,p},{g}_{1,p},2d)}-\frac{{P}_o(\mp {g}_{3,n},{g}_{1,p},2d)}{{P}_o({g}_{0,p},{g}_{1,\mathrm{p}},2d)}\right)\\ -{e}^{(\pm {g}_{3,m}\mp {g}_{3,n})w/2}\left(\frac{{P}_e(\mp {g}_{3,m},{g}_{1,p},2d)}{{P}_e({g}_{0,p},{g}_{1,p},2d)}-\frac{{P}_o(\mp {g}_{3,m},{g}_{1,p},2d)}{{P}_o({g}_{0,p},{g}_{1,p},2d)}\right)\left(\frac{{P}_e(\mp {g}_{3,n},{g}_{1,p},2d-2w)}{{P}_e({g}_{2,p},{g}_{1,p},2d)}+\frac{{P}_o(\mp {g}_{3,n},{g}_{1,p},2d-2w)}{{P}_o({g}_{2,p},{g}_{1,p},2d)}\right)\\ -{e}^{(\mp {g}_{3,m}\pm {g}_{3,n})w/2}\left(\frac{{P}_e(\mp {g}_{3,m},{g}_{1,p},2d-2w)}{{P}_e({g}_{2,p},{g}_{1,p},2d)}+\frac{{P}_o(\mp {g}_{3,m},{g}_{1,p},2d-2w)}{{P}_o({g}_{2,p},{g}_{1,p},2d)}\right)\left(\frac{{P}_e(\mp {g}_{3,n},{g}_{1,p},2d)}{{P}_e({g}_{0,p},{g}_{1,p},2d)}-\frac{{P}_o(\mp {g}_{3,n},{g}_{1,p},2d)}{{P}_o({g}_{0,p},{g}_{1,p},2d)}\right)\\ +{e}^{(\mp {g}_{3,m}\mp {g}_{3,n})w/2}\left(\frac{{P}_e(\mp {g}_{3,m},{g}_{1,p},2d-2w)}{{P}_{\mathrm{e}}({g}_{0,p},{g}_{1,p},2d)}-\frac{{P}_o(\mp {g}_{3,m},{g}_{1,p},2d-2w)}{{P}_o({g}_{0,p},{g}_{1,p},2d)}\right)\\ \enspace \left.\left(\frac{{P}_e(\mp {g}_{3,n},{g}_{1,p},2d-2w)}{{P}_e({g}_{2,p},{g}_{1,p},2d)}+\frac{{P}_o(\mp {g}_{3,n},{g}_{1,p},2d-2w)}{{P}_o({g}_{2,p},{g}_{1,p},2d)}\right)\right\}\end{array} $$(B6)

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

thumbnail Figure 1.

The cross-section of the gradient dielectric metasurface under consideration, composed of a rectangular periodic grating with period Λ, thickness w and permittivity εg. A unit-amplitude plane wave impinges on the metasurface at an angle θ.

In the text
thumbnail Figure 2.

Convergence patterns of the 0-order reflection coefficient |r0| with respect to the required number Nt 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 a1 = 0 and s1 = 0.5 Λ, and (b) N = 2 grooves with a1 = 0, a2 = 0.5 Λ, and s1 = 0.1 Λ, s2 = 0.3 Λ.

In the text
thumbnail Figure 3.

Amplitude |r0| of the 0-order 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 a1 = 0 and s1 = 0.5 Λ, and (b) N = 3 grooves with a1 = 0, a2 = 0.3 Λ, a3 = 0.6 Λ, and s1 = 0.2 Λ, s2 = 0.1 Λ, s3 = 0.2 Λ.

In the text
thumbnail 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, a1 = 0, and s1 = 0.5 Λ. (b) Amplitudes |r0| and |r–1| of the 0- and −1-order reflection coefficients versus the angle of incidence θ for the metasurface described in (a).

In the text
thumbnail Figure 5.

(a) Metasurface with εc = εs = 1, d = w = 180 nm, εd = 17.14 (Silicon), λ = Λ = 532 nm, N = 1, a1 = 0, s1 = 0.5 Λ, and εg = 2.05 (Teflon). (b) Amplitudes |r0| and |r−1| of the 0- and −1-order reflection coefficients versus the angle θ for the metasurface of (a), (c) same as (b), but for εg = 0.7.

In the text
thumbnail Figure 6.

(a) Geometry of the gradient dielectric metasurface (εc = εs = εd = 1, Λ = 800 nm, d = Λ/2, w = 2d, N = 1, a1 = 0, s1 = 0.5 Λ, and θ = 60°). (b)–(e) Reflection and transmission coefficients |r0|, |r−1|, |t0|, 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

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.