Anomalous refraction into free space with all-dielectric binary metagratings

Steering the incoming photonic power toward unusual directions, not predicted by Snell’s laws of reﬂection and refraction, is a fundamental operation behind numerous wavefront engineering applications. In this study, strong anomalous transmission is reported to all-dielectric electrically thin structures comprising only two materials in alternating rectangular posts. These binary metagratings have been thoroughly optimized to support almost perfect negative refraction by suppressing all the rest of diffraction orders; such an effect is found to be very robust with respect to oscillation frequency, incident beam angle, and structural defects. In this way, easy-to-fabricate setups working at different colors of the visible spectrum are determined and may be directly incorporated in various optical integrated systems from lenses and beamformers to ﬁeld enhancers and power splitters.


I. INTRODUCTION
Metasurfaces, the two-dimensional counterparts of volume metamaterials, are fabricated by arranging patterned subwavelength scatterers or holes at an interface in such a way that the desirable effective boundary properties are emulated [1,2].By designing systematically the features of the individual scatterers and properly distributing them in space, one can control amplitude, phase, and polarization of light or even generate arbitrarily shaped optical wavefront within a fraction of the operational wavelength [3].Hence, it makes no wonder that metasurfaces have been employed in several classes of optical devices, including high-quality-factor cavity resonators [4], sharp photonic switches [5], aberration-free lenses [6], and wideband absorbers [7].Importantly, all-dielectric gradient metasurfaces have recently received particular attention due to their potential of steering light without the losses and granularity of metals [8]; therefore, are extensively employed into photonic structures of high transmission efficiency in the visible [9,10].The basic components of such alldielectric structures are alternating blocks of transparent and high-refractive-index media composed of low-loss insulating texture.
Conventional metasurfaces usually require a continuous and fast-varying gradient surface impedance profile to achieve their operational characteristics; such a demand limits inevitably the overall effectiveness in transforming an impinging wavefront since it requires high-resolution fabrication processes.To overcome such a limitation, the concept of metagratings has been introduced as periodic arrays of carefully designed scatterers being simple to fabricate and enabling wavefront engineering with unitary efficiency.By employing such guidelines, highly performing setups were proposed including meta-reflectors for wave channeling [11], combined metagratings for broad-angle scattering [12], and reconfigurable multifunctional metasurface platforms [13].Inverse problems have been also considered and is found that the energy can be channeled into a single diffracted mode [14,15] or that the wavefront can be successfully optimized over a spectrum of oscillation wavelengths and incidence angles [16].
One of the most important effects that metagratings may support concerns the efficient steering of incident illumination toward prescribed directions in the reflection and/or refraction regions [3].In that context, the particular terminology anomalous reflection/refraction is utilized to describe unusual and remarkable channeling of the incoming waves toward directions not predicted by the classic Snell's law.It has been manifested [17] that for a planar layer composed of two periodically alternating dielectric rods, it is possible to obtain significantly enhanced first-order anomalous reflection (−1 order) with respective powers reaching up to 100% of the incoming one for certain wavelengths in the visible range.However, determining the optimal characteristics of simple all-dielectric metasurfaces offering enhanced refraction in the −1-order (anomalous refraction) is also very challenging in view of numerous potential applications, like efficient polarization conversion of electromagnetic fields [18], phase tailoring of the transmissive wave [19], and large splitting of polarized light [20].To this end, systematic analysis and careful engineering design is required in order to determine a realizable metasurface with its effective boundary condition manipulating the fields interplay in such a way that optimal results for anomalous refraction are obtained.
Anomalous refraction has been reported to serve diverse purposes in numerous layouts like cascades of "fishnet" structures [21] or patterned, metallic sheets [22].The same effect occurs by employing suitable arrangements of V-shaped nanoantennas [23,24], resonator arrays on two sides of a thin dielectric sheet [25], or even in a wavefront engineering platform comprising high-contrast dielectric elliptical nanoposts [26].In a different context, anomalous transmission can be achieved by using electrically tunable liquid crystals [27] or gain media [28,29].
In this work, we consider a two-dimensional, planar periodic metagrating composed of two rectangular dielectric rods per unit cell, which is the reason that is called binary, illuminated by a plane wave.Our main purpose is to properly select sizes and materials such that significant anomalous refraction is exhibited, in the sense that the −1-refracted order dominates with respect to all the other reflected and refracted waves.In this way, the major part of the total diffracted power is guided in the refraction region toward the same half-plane (with respect to the normal direction to the structure) that the metasurface is excited, and simultaneously the 0-diffracted orders contributions (dictated by Snell's law) are annihilated.
We use integral-equation formulation to mathematically model the associated diffraction phenomena by the considered metasurface and apply a rigorous and highly accurate entire-domain methodology for the solution; such a technique has been introduced and tested [30][31][32] for the analysis of grating-assisted dielectric structures.After implementing systematic optimization schemes, we detect several optimal designs channeling almost 100% of the incident field's power to the direction of the −1-refracted order.Furthermore, we also examine separately the case of normal illumination on the metagrating.By optimizing the metasurface's thickness and scanning the wavelength spectrum, one detects multiple designs that split the normally incident power and steer it along almost horizontal directions.Numerical solutions computed via the commercial field simulation software validate our results and demonstrate the efficiency of the proposed designs.
Compared to the already presented structures offering anomalous diffraction, we note that our metagrating is unarguably simpler than the setups mentioned above serving the same or similar purposes [21][22][23][24][25][26]; moreover, we do not use media requiring additional pumping equipment, like electrically tunable and/or active components [27][28][29].In addition, the proposed metagrating achieves negative refraction into free space, unlike several competing designs [33,34] where the beam is anomalously transmitted into dielectric multilayers.Importantly, we also emphasize that a major advantage of our approach is the rigorous and fast integralequation methodology we are implementing, without relying to successive costly field evaluations by means of commercial software or resorting to approximate numerical models.Besides, such an efficient solution gives reliable numerical results of controllable accuracy, while makes the optimization of the structure with respect to thickness and texture much easier and faster compared to alternative approaches.

II. RIGOROUS MATHEMATICAL MODELING
The geometrical configuration of the examined twodimensional all-dielectric metagrating is depicted in Fig. 1(a).It comprises a -periodic dielectric slab of thickness w composed per unit cell of two adjacent dielectric materials (making a binary metasurface) with refractive indices n 1 and n 2 , respectively.The duty cycle of the material with refractive index n 2 is denoted by s (for s = 0, we obtain a homogeneous slab of refractive index n 1 ).The metasurface lies in vacuum with refractive index n 0 and is assumed uniform along the direction ŷ.The entire space is considered magnetically inert with constant permeability μ 0 .
A TE z -polarized (or, alternatively, TM y -polarized) incident plane wave impinges on the metagrating at an angle of incidence θ with respect to the x axis.We assume a suppressed e iωt time dependence (with ω = 2π/(λ 0 √ ε 0 μ 0 ) as the angular frequency, t as time, and i = √ −1, where ε 0 , μ 0 , and λ 0 denote, respectively, the permittivity, permeability, and wavelength of free space).The electric field of the incident plane wave is expressed by E i (x, z) = e ik 0 n 0 (cos θ x−sin θ z) ŷ, where k 0 = ω √ ε 0 μ 0 is the free-space wave number.
The total electric field induced in every region of the scattering problem is y polarized, namely, E(x, z) = (x, z)ŷ.The unknown scalar electric-field factor (x, z) admits the integral representation [30][31][32] where S denotes the total transverse cross-section of the dielectric material with refractive index n 2 .Function 0 is the scalar electric-field factor induced on the homogeneous dielectric slab of refractive index n 1 by the plane incident wave.Function G is the Green's function of the homogeneous slab, which can be expressed as a Fourier integral with determined kernel [32].
The integral representation of Eq. ( 1) is reformulated with respect to the function u(x, z) = e ik 0 n 0 sin θ z (x, z), which is a -periodic function of z, as dictated by the Floquet-Bloch theorem [35,36].Then, the double integrals of Eq. ( 1) are reduced to integrals on the cross section S 0 = [− w 2 , w 2 ] × [0, s ], namely, on the part of the basic unit cell filled by the material of refractive index n 2 .The function u(x, z), appearing inside the resulting integrals on S 0 and additionally satisfying the Helmholtz equation, is expanded in Fourier series.By restricting the observation vector (x, z) on S 0 , employing wellestablished analytical techniques [31], and applying a highly efficient entire-domain Galerkin technique, we subsequently determine the associated Fourier coefficients.The integrations related to the projection of the expansion functions series on the set of testing functions are carried out analytically.In this way, a numerically robust linear system (which is initially infinite but finally truncated to size N) with respect to the unknown Fourier coefficients is obtained and accordingly the field into region S 0 is derived.
Such a field directly implies the form of the electric-field factors in the reflection (x > w 2 , z ∈ R) and transmission (x < − w 2 , z ∈ R) regions, being expressed as where each term in the field series of Eq. ( 2) (indexed by p ∈ Z) represents a plane wave characterized by the complex reflection coefficient r p or transmission coefficient t p , and will be referred to as the p-reflected or p-transmitted diffracted order.The wave numbers k x,p , k z,p , which define the direction of the plane waves, are given by The implemented integral-equation methodology provides semianalytical solutions with high numerical stability, controllable accuracy, and high efficiency.The sole approximation of the methodology lies in the final truncation of the expansion and test functions sets to the order of N.However, by applying a convergence control to the solutions for increasing N according to the energy conservation condition [37,38] (reflected and transmitted fields conserving power to the order of 10 −8 ), we conclude that small values of N provide sufficient validity for the obtained solutions; for more details see the numerical analysis and representative convergence patterns [32].Hence, in addition to its superior numerical accuracy, the integral-equation methodology is very efficient in terms of CPU time and computer memory.The methodology is validated by obtaining coincident results with different methodologies requiring significantly larger number of basis functions [31].The aforementioned beneficial characteristics make this methodology very suitable for optimizations of configurations involving multiple input parameters.In other words, successive computations of the reflection and transmission coefficients by means of the rigorous integral-equation methodology can be carried out fast and accurately, and thus brute-force optimization schemes can be properly supported.

A. Input parameters and metrics
It should be stressed that from the infinite terms of the sums in Eq. ( 2), only some of them correspond to propagating plane waves; all the others make vanishing evanescent field distributions.In particular, the reflection and transmission from the structure includes only waves with diffraction orders p, for which holds Under this assumption, the angle of propagation (measured with respect to the x axis) of the wave with order p reads By inspection of Eq. ( 4), it is clear that reflection and transmission for p = 0 is always supported; indeed, they compose waves obeying the classical Snell's law.The sign of θ p plays important role since it determines if the generated waves propagate into the half plane z < 0 that the structure is excited (θ p < 0, anomalous diffraction) or not (θ p > 0).It is also obvious that as the optical period /λ 0 of the metagrating increases, more and more propagating waves are generated due to scattering by the interface.However, since −p − > p + , the first diffraction order that is activated for θ > 0, is the p = −1 one.
The aim of the present study is to investigate the anomalous refraction by a simple gradient-index slab; not to examine the interplay of multiple diffraction orders.Therefore, we confine our research to combinations of (θ, /λ 0 ) yielding only two diffraction orders: p = 0, −1.This is assured by considering [17] max 1 − sin θ, Note that Eq. ( 6) guarantees additionally the anomalous diffraction of the −1 order, namely, that this order is reflected in the second quadrant of the plane (where the incident field impinges on the metasurface) and refracted in the third quadrant; in other words, Eq. ( 6) implies that k z,−1 < 0. The contour plot of the absolute value of the difference |θ + θ −1 | between the angle θ −1 of the −1-refracted order and the angle of incidence θ is depicted in Fig. 1(b) as a function of θ and the optical period /λ 0 ; note that the plus sign has been used in this plot because, as explained above, θ p < 0 in the case of anomalous diffraction.White regions in Fig. 1(b) correspond to pairs of θ and /λ 0 which do not fulfill the conditions of Eq. ( 6), namely, exactly two (0, −1) propagating orders one of which (−1) is anomalous.More specifically, for design/incidence combinations belonging to the lower region, we obtain only the Snell's reflection/transmission, while the upper left parametric corner of the map corresponds to activation of diffraction orders p = 0, ±1; similarly, across the upper right blank area, again three but different (p = 0, −1, −2) waves are propagated [39].We are mainly interested for designs supporting θ −1 ∼ = −θ , which usually [40] send large portion of power (either reflected or transmitted) along the anomalous directions.That is why we additionally represent in the same Fig.1(b), the optical extent (thick gray line) of spatial periods /λ 0 (measured at the right vertical axis) as function of incidence angle θ , for which the aforementioned difference |θ + θ −1 | falls below a specific threshold [taken equal to 15 o in Fig. 1(b)].It is clear that this range of corresponding to the potentially over-performing metagratings is larger for angles close to θ = 30 o (say 20 o < θ < 40 o ); accordingly, we have much better chances to determine significant maxima if optimizing the structure for such incident directions defining much more sizable parametric space.Similar investigations can be performed for other choices of incidence angles θ .The special case of normal incidence, i.e. θ = 0, is also important and will be analyzed separately in Sec.IV below.
The powers carried by the waves of p-reflected and ptransmitted (refracted) order are defined, respectively [41], as According to the definitions in Eq. ( 7), the energy balance criterion takes the form p∈P (P r p + P t p ) = 1, where P is the subset of Z including only the propagating diffracted orders; in our case, P = {−1, 0}.As was discussed in Sec.II, the latter energy balance criterion has been tested to be fulfilled up to the order of 10 −8 for all the numerical results presented hereinafter.

B. Parameters of the unit cell
As implied above, the major target of our work is to select suitable structural and textural parameters so that the power of the −1-refracted order is significantly enhanced, and simultaneously the powers of −1-reflected order and the 0-reflected and refracted orders are significantly annihilated.In this way, we can ultimately propose metagratings steering the main part of an incident wave to the symmetric direction in the refraction region, thus generating negative refraction in free space.Hence, the metric for the degree of negative refraction, is the power P t −1 of the −1-refracted order which will be optimized for six central wavelengths, one for each of the colors of the visible spectrum: λ 0 = 420 nm (violet), λ 0 = 470 nm (blue), λ 0 = 530 nm (green), λ 0 = 580 nm (yellow), λ 0 = 610 nm (orange), and λ 0 = 660 nm (red).
According to the discussion of Sec.III A, we consider incident waves with angles of incidence θ ∼ = 30 o , namely, 20 o < θ < 40 o .The geometrical parameters of the metagrating are determined as follows.The range of variation for the period is different for each angle and dictated by Fig. 1(b).When it comes to the thickness w of the structure, the general requirement is for it to be suitably small, namely in the range of a fraction of the effective wavelength into the employed media; only then the structure can be characterized as a metasurface.Finally, as far as the duty cycle s of the dielectric with refractive index n 2 is concerned, it will cover the entire permitted value range: 0 < s < 1.It is important to consider s as an extra optimization variable (in addition to the other geometrical parameters w and ), because this constitutes an additional degree of freedom (DoF), which can aid in canceling the undesired diffraction orders and hence provide improved results on negative refraction.This will be elaborated in Sec.III C below.
As far as the materials with refractive indexes (n 1 , n 2 ) that will be assumed comprising our free-standing (n 0 = 1) structure, we examine two main alternatives.In the first part of Sec.III C, we consider two media exhibiting substantial contrast: lossless crystalline Silicon (c-Si), where n 1 follows well-known frequency-varying models [42,43] and Teflon AF fluoropolymers [44] with nondispersive n 2 ∼ = 1.3.Notably, metagratings composed of materials with the above-mentioned properties (c-Si with Teflon) were shown to support anomalous reflection effects [10,17].Then, we examined various combinations of actual materials which can be safely considered lossless in the entire visible range and hence provide efficient and realizable designs.From all the examined combinations, we include in the second part of Sec.III C the ones leading to the optimal results with respect to the generation of negative refraction phenomena.These refer to metasurfaces with hafnium dioxide (HfO 2 ), with n 1 ∼ = 2.1, as the dense medium and air, with n 2 ∼ = 1, or magnesium fluoride (MgF 2 ), with n 2 ∼ = 1.4,as the sparse medium.

C. Optimal designs
We performed a greedy optimization based on successive computations of the reflection and transmission coefficients by means of the semianalytical integral-equation methodology described in Sec.II.The superior numerical stability, controllable accuracy, and rapid convergence of the method aid substantially the fast and efficient implementation of the optimizations yielding to designs supporting significant anomalous refraction (maximal transmitted power P t −1 ).One could point out that, since a diffraction order requires canceling both the real and imaginary parts of the field and negative refraction demands suppression of two diffraction orders, it would in general be expected to use four DoFs to achieve the regarded purpose, even though the governing laws are not linear [45].For our binary metagrating, four DoFs, {w, , s, θ} can be considered, and they have been used in the employed optimization schemes.Adding more DoFs, like more layers per unit cell in the configuration of Fig. 1(a), could be expected to allow for achieving other types of challenging targets, like steering an incident beam to an oblique direction.However, considering significantly more First, we examined the case that n 1 follows the c-Si frequency variation and n 2 = 1.3.The results of the brute-force optimization described above are shown in Table I, which includes the maximum P t −1 and the optimal values of w, , s for the six examined colors in the visible range.It is evident that impressively large values of P t −1 (above 98%) are achieved for every color, meaning that almost the entire incident field's power is guided in the third quadrant of the plane: x < 0 and z < 0. It is also important to observe that there are very small differences between incidence θ and anomalous refraction direction θ −1 (up to two degrees).Therefore, the reported optimal metagratings transmit nearly the entire amount of the incident plane wave to the symmetric (with respect to the −z axis) direction of illumination.
To demonstrate the efficiency of the designs presented in Table I, we simulated in COMSOL Multiphysics [46] environment certain structures incorporating them.Instead of a plane wave which needs an infinite metagrating to interact, we used as excitation the so-called complex point source which emulates successively the beam produced by a horn antenna with specific aperture [47].In particular, the incident field can take the form E i (x, z) = H 0 (k 0 |r − r 0 |)ŷ, where H 0 is the Hankel function of zeroth order and second kind and r = x x + zẑ the in-plane position vector.The real part of the complex vector r 0 = (R + ib) cos θ x − (R + ib) sin θ ẑ shows the location of the source (distance R from the origin at direction θ ), while its imaginary part indicates the tilt (angle θ ) of the beam's aperture and its size 2b.Locally, it behaves like a spatially restricted plane wave and as long as the aperture is large enough (2b/ cos θ ) placed quite distantly (R w), it makes a suitable excitation of the examined metagrating.
In Fig. 2(a), we show the variation of the electric field E y (x, z) when the optimal design of Table I is excited at λ 0 = 470 nm by a complex point source which is placed outside the simulation box (at distance R = 100 ∼ = 56 μm from the origin) of a large aperture b ∼ = 20λ 0 ∼ = 9.4 μm locally mimicking a plane wave and illuminating at the optimal angle θ = 25 • .It is apparent that negligible reflections occur and, most importantly, the vast portion of the incident power is channeled toward the lower vacuum space along the symmetric direction with respect to z axis.Such anomalous behavior is feasible by generation of the required phase variation across z = 0 plane through the Silicon and Teflon rods [24].In particular, the denser medium (c-Si) provides the necessary time delay to the signal so that the necessary z-dependence e ik 0 n 0 sin θ z is artificially given to the refractive wave [48].It is expressed by a complete circle of the power flow into each (almost) square rod of c-Si which admits the light to get stored and reemitted in a combined way with the rays passing through the rectangular rod of Teflon to build the required phase variation.
In Fig. 2(b), we represent the magnitude of the electric field |E y (x, z)| when the optimal design of Table I is excited at λ 0 = 660 nm by a complex point source of similar characteristics (R = 100 ∼ = 63 μm, b = 20λ 0 ∼ = 1.3 μm, θ = 32 • ).Contrary to the simulation of Fig. 2(a), we now employ two identical binary metagratings (at specific distance D = 8λ 0 ∼ = 5.3 μm between them) to demonstrate more clearly the negative refraction into free-space which reminds us the operation of a perfect lens [49] with no use of double-negative or gain media [28,29].By inspection of Fig. 2(b), it is again clear that the created phase of the wave along the upper metagrating guarantees both perfect matching and negative refraction; the transmitted wave with these features is used as feeding signal at the lower metagrating, where same effect takes place symmetrically with respect to x axis.
In Fig. 3, we pick four of the proposed designs from Table I and examine their "spectral signatures" with respect to the oscillation wavelength λ 0 and the incidence angle θ around the optimal operational points indicated in Table I (with a range of 30 nm for λ 0 and 20 • for θ ).In Fig. 3(a), we show the transmitted power of the first anomalous order P t −1 for the c-Si/Teflon binary metagrating optimized at λ 0 = 470 nm (blue color).It is remarkable that the represented metric P t −1 exhibits substantial robustness with respect to frequency or direction misalignment; in particular, there are wavelengths of oscillation (λ 0 ∼ = 475 nm) for which the anomalous transmission retains very high values (over 95%) almost for the entire range of incidence directions (15 • < θ < 35 • ).However, the performance of the device becomes poorer and more angleselective when λ 0 gets decreased, namely, for a violet-shift operation.
In Fig. 3(b), we investigate the metagrating from Table I functioning optimally at λ 0 = 580 nm (yellow color) and the power P t −1 takes substantial magnitudes across the considered parametric box [the range of the colorbar is much smaller than that of Fig. 1(a) and covers large values only].Interestingly enough, the best score is not recorded exactly at the operation point (λ 0 = 580 nm) mentioned at Table I but at a slightly smaller wavelength, since we do not optimize our metrics with respect to λ 0 .Naturally, when one moves away from the best oscillation wavelength and direction, P t −1 falls mildly but the design remains highly competitive.In Fig. 3(c), we examine the metagrating from Table I which works with orange color and we notice a significant asymmetry of the response along λ 0 axis.In particular, the anomalous transmission is kept high for λ 0 < 610 nm, while it diminishes rapidly for λ 0 > 610 nm, where the incident power splits in four almost equal parts (Snell and negative reflected/transmitted rays).Finally, in Fig. 3(d), we consider the design working for the red color of visible spectrum and we obtain a similar response to that of Fig. 3(b): remarkable P t −1 across all the regarded spectral combinations which gets slightly deteriorated outside the optimal point.
As indicated above, the Table I, contains the parameters for optimal designs incorporating materials with large contrast: lossless crystalline Silicon n 1 (λ 0 ) as the dense one and Teflon (n 2 ∼ = 1.3) as the sparse one; however, in practice n 1 (λ 0 ) possesses substantial imaginary part for λ 0 < 550 nm, which is quite harmful for the metagrating's effectiveness.For this reason, we now pick a lossless material commonly used in the laboratory for structures fed by visible light, to play the role of dense medium: hafnium dioxide [43] (HfO 2 ) with n 1 ∼ = 2.1.The sparse material can be either air (free-standing HfO 2 rectangular rods, n 2 ∼ = 1) or magnesium fluoride [42] (MgF 2 ) with n 2 ∼ = 1.4.The results of optimization described in Sec.III A and in the first paragraphs of Sec.III C when applied for the HfO 2 -based realistic metagratings are shown in Table II.We notice that, mainly due to the weaker texture contrast, the performances of the devices in Table II are lower than the ones of Table I, but still very high accomplishing, on average, negative refraction of 90% of the incident power.Note that even when MgF 2 is employed (designs for violet, blue and red color) and thus n 1 is only slightly larger than n 2 , significant levels of P t −1 are recorded.It should be stressed that there are are several ways of fabricating designs like that of Fig. 1(a) starting from the set of chemical techniques whose more representative member is the so-called atomic layer deposition (ALD).It is a sequential approach where reactive media are sprayed by other substances (precursors) which enable successive chemical reactions each of which generate a new layer [50].Photonic multilayers of dielectric are a common outcome of ALD [51], while stacks of semiconducting films have been reported by using similar approaches [52].Another fabrication category of planar multilayers consists the physical depositions and especially molecular beam epitaxy (MBE), where a heated substrate gets shot by guns of the successive media in gas form, moleculeby-molecule [53].Again, stacks of dielectrics in the visible and semiconductors [54], can be produced by such a technique.Moreover, lithography may be employed successfully for manufacturing of gradient metasurfaces with rectangular rods which constitute the building block [Fig.1(a)] of the reported designs [55].Finally, similar lamellar structures can be even self-organized [56] by utilizing plasmonic eutectics in combination with suitable micropulling.
However, it is meaningful to test the responses of designs appeared at Table II, under fabrication defects like an error δ in the metagrating's period or an imperfect selection in the duty cycle s for the sparser medium.In Figs. 4, we pick four designs from Table II and we show the variation of P t −1 in the vicinity of the optimal s and for −20 nm < δ < 20 nm; the optimal point is marked by a black × sign.In Fig. 4(a), we examine the HfO 2 /MgF 2 binary metagrating working at violet color (λ 0 ∼ = 420 nm); it is clear that the design does not exhibit substantial endurance with respect to perturbations in since it falls rapidly far from δ = 0 line.As far as the performance variation with changing s is concerned, it keeps its large magnitude for most parametric combinations except for a band of (δ , s) pairs describing a specific size of the MgF 2 gaps at which almost total Snell's reflection occurs.
In Fig. 4(b), we consider that setup from Table II operated under green light (HfO 2 /air) and realize that such a structure is very robust against construction imperfections of metagrating's period .On the contrary, the anomalous refraction drops substantially when the portion of HfO 2 is smaller (where Snell's reflection is enhanced) or larger (where Snell's transmission increases) than the optimal value s ∼ = 0.66 indicated in Table II.In Fig. 4(c), we show again the metric P t −1 for HfO 2 rods with airgaps but optimized for yellow light; similarly to Fig. 4(b), the parametric "plateau" of significant performance is not symmetrically formed around the optimal point despite being quite extended.We also notice a secondary local maximum at higher s which gets separated from the major one by a zone of total reflection; furthermore, the anomalous transmission score gets deteriorated when s is much lower than the optimal, where the structure is almost transparent to the incident wave.Finally, the most robust response is noted in Fig. 4(d), owned by the metagrating working at orange light; only in the vicinity of a line below the optimal point on (δ , s) plane, where Snell's transmission prevails, an exception is made.

IV. NORMAL INCIDENCE
In the previous, we have presented and investigated simple binary metagratings configurations achieving substantial anomalous (negative) refraction under oblique incidence (θ = 0).However, it will be meaningful to examine the case of normal illunimation (θ = 0), with which the solutions are expected symmetric with respect to x axis.Therefore, except for the normal reflection P r 0 and transmission P t 0 , the waves of orders ±1 will travel across symmetric directions and, inevitably, P r 1 = P r −1 and P t 1 = P t −1 .In such a scenario, the (common) refracting angle reads from Eq. ( 5): |θ ±1 | = arcsin(λ 0 / ).If one focuses on the case of almost horizontal reflection and refraction, then one should keep /λ 0 slightly larger than unity, a selection which approximately satisfies the condition of Eq. ( 6) [as also indicated by Fig. 1(b)].
By optimizing the thickness w of our metagrating and scanning the wavelength λ 0 spectrum, one may conclude to multiple designs supporting substantial splitting of the normally incoming power.In Fig. 5, we consider certain  where period is chosen extremely close to the operational wavelength λ 0 , we obtain an anomalous refractive direction almost parallel to z axis (|θ ±1 | ∼ = 87 o ).Two strong resonances are recorded for P t ±1 at w ∼ = 0.56λ 0 and w ∼ = 0.78λ 0 , with simultaneous suppression of the reflected counterpart; however, P r ±1 is also maximized at slightly larger w.In Fig. 5(b), where the same quantities are represented with respect to wavelength, we notice a maximal λ 0 threshold beyond which we do not observe anomalous diffraction effects but only Snell's reflection/transmission (P r 0 /P t 0 ).Remarkably, the maximal nonnormal refraction [P t −1 + P t 1 > 0.6, also indicated by Fig. 5(a)] emerges very close to this threshold, an effect also reported in similar metasurfaces [39].
In Figs.5(c) and 5(d), we regard a metagrating of larger period and accordingly the directions of the waves with nonzero diffraction order are diverging from z axis (|θ ±1 | ∼ = 77 o ).However, the efficiency of the non-Snell response is huge since the maximal value of (P t −1 + P t 1 ) surpasses the 90% of the incident power; note that for the optimal thickness w [shown in Fig. 5(c)], almost the entire incoming energy is channeled quasihorizontally.In Fig. 5(d), where the operating frequency is being swept, we understand that the resonance of the second design is more wideband meaning that the effect is less vulnerable to small changes of the structure, compared to the one examined in Figs.5(a) and 5(b).Again, there is a maximum λ 0 at which anomalous reflection/refraction can occur and it is noteworthy that once all the nonnormal diffraction gets switched-off, the metagrating is almost transparent to the incident wave (P t 0 = 1), a property again met in periodic clusters of infinite cylinders [39].Note that the duty cycles s in Fig. 5 have been implicitly optimized, while the respective graphs are not shown for brevity; the full-scale optimization has been demonstrated in the more complicated case of oblique incidence (Sec.III).
To show the efficiency of similar designs in splitting the normally incident power and channeling it along almost horizontal directions, we pick the metagrating considered by Figs.5(c) and 5(d) (with |θ ±1 | ∼ = 77 o ) and we simulated it with help from COMSOL Multiphysics.In Fig. 6, we represent the spatial distribution of the squared magnitude of the sole electric component |E y (x, z)| 2 for such a HfO 2 /air structure when normally illuminated with green light (λ 0 = 544 nm).One can notice the weak Snell's transmission and the very low (±1)-ordered reflection; indeed, the vast portion of energy is propagated by using as vehicles the transmitted waves of ±1 order.In this way, the operation of such a simple all-dielectric binary metagrating as an effective power divider is indicated.

V. CONCLUSIONS
A rigorous integral-equation formulation is provided to model the wave interaction with an all-dielectric metasurface while the electromagnetic fields are determined by implementing a highly accurate methodology using entire-domain basis functions.The considered metagratings employ only two media and are optimized to suppress all the produced diffracted waves except for one giving strong anomalous transmission.The reported effect is not substantially affected by small changes in the characteristics of the incoming signal (wavelength, direction) or fabrication defects (thicknesses, size proportions).The case of normal incidence is separately analyzed and revealed interesting effects of converting the input signal to effectively surface waves traveling parallel to the considered metagrating.
The same technique can be expanded to treat more complicated configurations incorporating inclusions and scatterers of shapes different from rectangular while stacking them in multilayers.In this way, one may conclude to devices serving more challenging purposes like anomalous refraction along several directions with direct applicability to multipleaccess indoor communications or radio-coverage in urban environments.Finally, the analytical nature of the adopted numerically robust methodology, makes it ideal for inverse designing of similar metasurface configurations [14][15][16], namely, for deciding the structural/textural layout that mimic best the desired response across frequency and angular spectra.

FIG. 1 .
FIG. 1.(a) Geometrical configuration of the examined gradient all-dielectric -periodic metagrating and the incident field.(b) Variations of the difference of anomalous refraction angle θ −1 from incidence angle θ (in degrees) as long as only the orders p = 0, −1 are propagated.Along the horizontal axis we represent the incidence angle θ and along the vertical one, the optical spatial period /λ 0 of the metagrating.The gray line, corresponding to values /λ 0 at the right vertical axis, indicates the electrical extent of for which the aforementioned angle difference |θ + θ −1 | is kept below a threshold (15 o ).

FIG. 2 .
FIG. 2. Spatial distribution of: (a) Electric field E y (x, z) in the presence of a single optimal metagrating of Table I working at blue color (λ 0 = 470 nm).(b) Electric field magnitude |E y (x, z)| in the presence of a couple of optimal metagratings of Table I working at red color (λ 0 = 660 nm).Simulations executed via COMSOL; arrow shows the direction of the incident wave which originates from a complex point source outside of the considered simulation box.

FIG. 3 .
FIG. 3. Contour plots of P t −1 as a function of the operating wavelength λ 0 and the angle of incidence θ for the optimal values of w, and s of Table I.Optimal c-Si/Teflon designs for: (a) blue, (b) yellow, (c) orange, and (d) red color.

FIG. 4 .
FIG. 4. Contour plots of P t−1 as a function of error δ in selecting the optimal period of the metagrating and the duty cycle s for the optimal values of w, λ 0 and θ of TableII.Optimal HfO 2 -based designs, marked by black ×, for: (a) violet, (b) green, (c) yellow, and (d) orange color.

2 FIG. 6 .
FIG. 6.The spatial distribution of the squared magnitude of the sole electric component |E y (x, z)| 2 for the optimized design splitting the power at normal incidence of Figs.5(c) and 5(d).

TABLE I .
Optimal values of the parameters w, , and s of a metagrating with n 1 obeying the dispersion model of c-Si and n 2 = 1.3, yielding maximum P t −1 for each considered wavelength; the angles of incidence θ and of negative refraction θ −1 are also depicted.
complicated structures would inevitably increase exponentially the complexity of the required optimization schemes.In this work, we are presenting a minimal model that can support anomalous diffraction; further structural and textural refinements can certainly constitute interesting future work directions.

TABLE II .
Optimal values of the parameters w, , and s of a metagrating with n 1 obeying the dispersion model of HfO 2 and n 2 = 1 (green, yellow, orange color) or n 2 being MgF 2 (violet, blue, red color), yielding maximum P t −1 for each considered wavelength; the angles of incidence θ and of anomalous refraction θ −1 are also depicted.