Efficient integral equation modeling of scattering by a gradient dielectric metasurface

Plane-wave 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 semi-analytical 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 narrow-band 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.


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 [5][6][7][8][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 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.

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 e d and thickness 2d on which has been etched an infinitely K-periodic rectangular dielectric grating composed of N rectangles (grooves) per unit cell with relative permittivity e 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 semi-infinite superstrate (cover) and substrate plane regions, above and below the slab, have relative permittivities e c and e s , respectively.The entire structure has constant magnetic permeability l 0 and is assumed uniform along the direction ŷ.
A TE-polarized incident plane wave of unit amplitude impinges on the metasurface from the cover region at an angle h (see Figure 1).The incident electric field is given by (under exp(jxt) time dependence) where k inc is the incident wave vector and k = x/c the free-space wave number.The electric field induced in every region is of the form and the problem is reduced to the determination of the unknown scalar electric field factor W. 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 W admits the integral representation [22,23] W where W 0 is the electric field induced on the nongrating 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 S g 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) which implies that Moreover, by representing the total grating's cross section S g as countable union of the cross sections of the grating's rth unit cell and employing the derived Fourier integral (A9) of the Green's function as well as (4), we reformulate the integral representation (3) as where function l denotes the determined kernel of the Green's function (see (A10)).Now, the transformation f 0 = z 0 À rK reduces the double integrals of equation (5) to integrals on the basic unit cell S 0 (the grating's cross section on [0, K]).Then, by using the Poisson's summation formula [25] for the Dirac function and taking into account the basic property-definition of the Dirac function, we rewrite (5) as 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 S 0 is expanded in the Fourier series with respect to z where the Fourier coefficients u n are the transverse components of the spatial harmonics [26], satisfying here the ordinary differential equation where : The functions u n are then expressed as linear combinations of the fundamental solutions of (8) (representing waves traveling parallel to the Çx directions) with c AE n under determination coefficients.In the terminology of the Method of Moments, the entire domain functions u n (x)exp[Àj(2pn/K)z] are referred to as the electric field expansion (basis) functions.
Next, by substituting ( 7) into (6), we get where the functions J pÀn and Q AE np are defined in the Appendix B N.L. Tsitsas: EPJ Appl.Metamat.2017, 4, 3 Furthermore, the observation vector (x, z) in ( 10) is restricted on the grating's domain.Hence, for (x, z) 2 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 non-homogeneous square linear system of equations with respect to the unknown coefficients 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 where the (2N t + 1) • (2N t + 1) matrices A AEAE are given by while h 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 non-grating structure.

Reflected and transmitted fields
After having determined the coefficients c ± by solving the non-homogeneous linear system (12), we express the reflected W r and transmitted W t electric fields by means of the basic representation (10), as follows where the wavevectors' components are expressed by ; 4 N.L. Tsitsas: EPJ Appl.Metamat.2017, 4, 3 while the complex reflection and transmission coefficients are, respectively, given by Coefficients R and T are given by (A5) and (A6), d p0 is the Kronecker symbol (:d 00 = 1 and d p0 = 0|p 5 0), while 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 K, the wavelength k, the permittivities e c and e s of the cover and substrate regions as well as on the angle of incidence h.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 5 0, we define the thresholds for which holds: x;p and k t x;p become purely imaginary and the p-order fields are evanescent.
Hence, p AE r and p AE t 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].

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 0-order reflection coefficient are depicted in Figures 2a and 2b corresponding to e c = e s = 1, 2d = w = 0.7 cm, e d = 2.59, h = 60°, K = 1.65 cm, with e g = 0.5, 0.3, and 0.1, and (a) N = 1 groove with a 1 = 0 and s 1 = 0.5 K, and (b) N = 2 grooves with a 1 = 0, a 2 = 0.5 K, and s 1 = 0.1 K, s 2 = 0.3 K.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).
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.

Reflection filters
Plane wave scattering by a metasurface composed of an array of dielectric rectangular cylinders can be treated as a N.L. Tsitsas: EPJ Appl.Metamat.2017, 4, 3 special case of the developed method.More precisely, by selecting e c = e s and w = 2d, we investigate the structure of an infinite K-periodic array composed of adjacent rectangular cylinders with permittivities e d and e 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. e c = e s = 1, while the incidence is normal to the metasurface, i.e. h = 90°.Figure 3 depicts the dependence of the amplitude |r 0 | of the 0-reflected order with respect to the operating frequency f for e d = 2.59 (Plexiglas), w = 2d = 0.7 cm, and K = 1.65 cm.Three cases are considered for the grating's material: (i) e g = 2.05 (Teflon), (ii) e g = 0.5 (metamaterial 1: MTM1), and (iii) e 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 K, 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 K, a 2 = 0.3 K, s 2 = 0.1 K, a 3 = 0.6 K, and s 3 = 0.2 K.
From Figure 3a it is observed that a grating of one groove composed of a metamaterial with e 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 e 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 e g = 2.05 and e g = 0.5).Particularly, for e 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 e 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 |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].

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 K, the operating wavelength k, and the angle of incidence h 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; p þ t < 1 and À2 < p À r < À1; À2 < p À t < À1, which when e s = e c and thus 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 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 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 h corresponding to plane wave scattering by the array of rectangular cylinders, depicted in Figure 4a, with e g = 17.14 (Silicon), N = 1, a 1 = 0, and s 1 = 0.5 K, lying inside a background medium with e c = e s = e d = 1, and d = w = 180 nm.The operating wavelength is k = 532 nm and the cylinders' period K = 500 nm; a condition of the form K ~k 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 (e c = e s = e d = 1).It is evident that |r -1 | is larger than |r 0 | in the range of h which is mainly of interest (excluding the near-grazing incidence angles).Particularly, the locally maximum distances between |r -1 | and |r 0 | are observed at h = 39°and at h = 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.
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 (e d = 17.14, d = 180 nm) possessing a grating composed of Teflon (e g = 2.05, d = w) grooves (see Figure 5a).The ratio of the wavelength over the periodicity is taken to be exactly 1 (k = K = 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 N.L. Tsitsas: EPJ Appl.Metamat.2017, 4, 3 a metasurface composed of an array of conventional materials (Silicon slab filled with Teflon grooves).
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 e 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 6b-6e the amplitudes |r 0 |, |r À1 |, |t 0 |, and |t À1 | of the 0-and À1-order reflection and transmission coefficients versus the operating wavelength k for the array of dielectric cylinders, depicted in Figure 6a, with four different permittivities e g < 1 having the same thickness w = 2d = K = 800 nm, and lying in free space (e c = e s = e d = 1).Figure 6b shows that, when e 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 e g = 0.2.As e 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 e g = 0.5, where the 0-transmission order becomes larger than the À1-order.As the permittivity e 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 e g = 1 the entire structure is vacuum and hence the incident wave travels unaltered).These conclusions are verified by the results of Figure 6e.

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

Appendix A
The total electric field induced on the non-grating structure due to the incident plane wave (1) is given by where 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 0 , z 0 ) 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 where the kernel l is given by lðk; x; x 0 Þ ¼ ½Kðg 0 ; g 2 ; g 1 ; 2dÞ À1 P o g 2 ; g 1 ; 2d ð Þ P e g 2 ; g 1 ; 2d ð Þ cosh g 1 x 0 ð Þ P e g 2 ; 1 P e g 0 ; g 1 ; 2d ð Þ P e g 2 ; g 1 ; 2d ð Þ P o g 2 ; g 1 ; 2d ð Þ P o g 0 ; g 1 ; 2d ð Þ cosh g 1 x 0 ð Þ P e g 2 ; 1 P e g 0 ; g 1 ; 2d ð Þ P e g 2 ; g 1 ; 2d ð Þ P o g 2 ; g 1 ; 2d ð Þ P o g 0 ; Appendix B The following auxiliary functions are utilized for the formulation of the method's linear systems (11) and ( 12)

Figure 1 .
Figure1.The cross-section of the gradient dielectric metasurface under consideration, composed of a rectangular periodic grating with period K, thickness w and permittivity e g .A unit-amplitude plane wave impinges on the metasurface at an angle h.

r
and p < p þ t or p > p À t the wavevectors' components k r x;p and k t x;p are real and the p-order reflected and transmitted fields are propagating along the x-axis.(ii) for p > p þ r or p < p À r and p > p þ t or p < p À t the wavevectors' components k r