Analysis and design of two-dimensional compound metallic metagratings using an analytical method

The recently proposed concept of metagrating enables wavefront manipulation of electromagnetic (EM) waves with unitary efficiency and relatively simple fabrication requirements. Herein, two-dimensional (2D) metagratings composed of a 2D periodic array of rectangular holes in a metallic medium are proposed for diffraction pattern control. We first present an analytical method for diffraction analysis of 2D compound metallic metagrating (a periodic metallic structure with more than one rectangular hole in each period). Closed-form and analytical expressions are presented for the reflection coefficients of diffracted orders for the first time. Next, we verify the proposed method's results against full-wave simulations and demonstrate their excellent agreement. As a proof of principle, two applications are presented using the proposed analytical method. The first application is a perfect out-of-plane reflector that transfers a normal transverse-magnetic (TM) polarized plane wave to an oblique transverse-electric (TE) polarized plane wave in the $y-z$ plane. The second one is a five-channel beam splitter with an arbitrary power distribution between channels. Using the proposed analytical method, we designed these metagratings without requiring even a single optimization in a full-wave solver. The performance of the designed metagratings is better than previously reported structures in terms of power efficiency and relative distribution error. Our analytical results reveals that 2D metagratings can be used for manipulating EM waves in the plane and out of the plane of incidence with very high efficiency, thereby leading to extensive applications in a wide range of frequencies from microwave to terahertz (THz) regimes.


I. INTRODUCTION
Wavefront shape manipulation has always been an interesting topic in electromagnetism due to its fundamental role in several applications such as radars, imaging, and communication systems [1]. In the last two decades, artificial structures have reached the peak of attention for controlling electromagnetic waves [2][3][4][5][6]. Metasurfaces, in particular, have been used to devise novel devices with significant practical and scientific applications thanks to their high potential [7][8][9][10][11][12][13][14]. Metasurfaces are two-dimensional thin planar patterned structures formed by spatially arranged building blocks called meta-atoms. They can be designed to control the amplitude, phase, and polarization of EM waves. A major class of metasurfaces is gradient metasurfaces (GMS) which can manipulate EM waves by imparting local momentum to the incoming EM waves through the gradient phase profile of the structure [15]. GMS can realize a wide range of electromagnetic and optical functionalities from beam focusing to holographic imaging [16][17][18][19][20]. However, passive and local GMS have been shown to suffer from low power efficiency and require precise and high-resolution fabrication processes [21][22][23]. These problems restrict many of the applications mentioned above. To address these problems, the concept of metagrating was proposed by Ra'di et al [24].
Metagratings, a sparse periodic array of subwavelength scatterers (meta-atoms), have attracted considerable interest in the last few years because they allow the realization of diverse phenomena such as anomalous reflection, * khavasi@sharif.edu beam splitting, beam steering, and beam focusing [25][26][27][28][29][30][31][32][33]. Their power efficiency is not restricted by any fundamental physical bounds, and they require much less fabrication complexity than metasurfaces [24]. The working principle of metagratings can be understood using Floquet-Bloch (FB) theory, according to which when a plane wave impinges to a periodic structure, it will be diffracted into several discrete waves in a certain direction. Meta-atom properties have a significant effect on diffracted waves (FB mode); hence, by engineering meta-atoms, we can tailor the desired diffraction patterns. Different geometries for meta-atoms, such as loaded thin wires [34], one-dimensional (1D) grooves [35], and graphene ribbons [36] have been used to design metagratings to realize various functionalities.
Most of the designed metagratings are periodic in 1D and are sensitive to incident wave polarization. Therefore 2D-metagratings (periodic array in two directions) are proposed for realizing polarization-independent anomalous reflection with high diffraction efficiency [37][38][39]. However, 2D metagratings have not been designed based on an analytical method and thus have a time-consuming design procedure. For example, in [38], an all-metallic metagrating is proposed for polarization-independent perfect anomalous reflection. It is analyzed by the rigorous coupled-wave analysis (RCWA) technique and anomalously reroutes an obliquely incident wave with high efficiency for both polarization. However, the RCWA technique and other numerical methods can not present a closed-form expression for the reflection coefficients; therefore, an analytical method for the analysis of 2D metagratings is in demand for accelerating design procedures. Furthermore, all of the existing metagratings are used for in-plane control of EM waves (The wave-vector of the diffracted wave lies in the plane of incidence). Although [40] has attempted to transfer the incident power to some directions out of the plane of the incidence, due to its sophisticated method design, the designed metagratings have a low power efficiency. Outof-plane manipulation of EM waves (The wave-vector of the diffracted wave does not lie in the plane of incidence) has interesting potential applications such as 2D planar lenses with high numerical apertures, flat polarization converters, and radar cross-section reductions [41][42][43][44]. Metasurfaces are common elements for realizing this phenomenon, which, as noted earlier, have low power efficiency. To the best of the authors' knowledge, out-ofplane manipulation of EM waves has not been realized by metagratings.
In this paper, we present an analytical method for analyzing two-dimensional compound metallic metagratings (2D-CMGs) and show these metagratings enable in-plane and out-of-plane EM wave control. The proposed metagrating consists of the 2D-periodic repetition of a finite number of rectangular holes carved out of a thick metal slab. To derive our method, we first expand the electromagnetic field by FB theory and extract the reflection coefficients of the zeroth diffracted order and higher diffracted orders by using appropriate boundary conditions in the conjunction mode-matching technique. The accuracy of the proposed method is verified through numerical examples. Using this analytical method, we designed an out-of-plane anomalous reflector with unitary efficiency at normal incidence. Next, we propose a fivechannel beam splitter using 2D-CMG. Two of these channels are in the plane perpendicular to the plane of incidence. The power distribution between these channels are arbitrary. The performance of the designed beam splitters was better than that of the previously reported metasurface-and metagratings-based beam splitters. Finally, some practical aspects of the experimental realization of the designed devices are discussed.

II. ANALYTICAL METHOD FOR ANALYSIS OF 2D-CMGS
A. Diffraction of a normal TM polarized plane wave by 2D-compound metallic grating Herein, we first present an analytical and closed-form expression for the reflection coefficients of diffracted orders of 2D-CMG comprising two rectangular holes. Next, we discuss how the proposed method can be generalized to analyze a 2D-CMG containing an arbitrary number of rectangular holes. It should also be noted that the timeharmonic of the form exp(jωt) is assumed throughout this paper.
Consider a 2D-CMG including two rectangular holes made in a metallic medium as shown in Fig. 1. The periods of the structure along the x-and y-axes are P x and P y , respectively. Each hole has a width of a i , length of b i , and height of h i , and is filled with a dielectric medium with a refractive index of n i (region II and III). The whole structure is surrounded by a homogeneous medium with a refractive index of n 1 (region I). We denote the distance between the corner of holes with d x and d y in the x-and y-directions, respectively.
Assume that a normal incident TM polarized plane wave (the magnetic field in the y-direction) propagating along the z-direction illuminates the structure. The tangential electric and magnetic fields in the region z > 0 can be written as [45,46] E 1x = e jkz,00z + R 00 e −jkz,00z + H 1y = −ξ 00 e jkz,00z + ξ 00 R 00 e −jkz,00z + m =0 ξ m0 R T M m0 e −jkx,mx e −jkz,m0z + n =0 ξ T E 0n R T E 0n e −jky,ny e −jkz,0nz Using Maxwell equations, the other components of the tangential fields are obtained as [47] where R T E mn and R T E mn are the reflection coefficients of the TE-and TM-polarized mnth diffracted order, respectively and the subscripts m , and n correspond to the order of the diffracted waves along the x-and y-axes, respectively. Furthermore, k x,m , k y,n , and k z,mn are the wave-vector components of the diffraction order along the x-, y-, and z-directions in region I, respectively, and are given by [46] It should be noted that the branch of the square root for the z-component of the wave-vector is chosen in such a way that either its real part should be positive (propagating wave) or its imaginary part should be negative (evanescent wave). Moreover, ξ T M mn = ωε 0 n 2 1 /k z,mn and ξ T E mn = k z,mn /ωµ 0 are the TM/TE-wave admittance of the mnth diffracted order in region I.
In regions II and III, we assume that the holes are single-mode, and due to the TM polarization of the incident wave, we only take into account the T E 01 mode, which is propagating inside the holes, while assuming that the effects of other order modes are negligible. The validity of this approximation is limited to the operating frequency less than , where the higher modes inside the holes are evanescent. As a result, the magnetic and electric fields in the holes can be written as [47] is the propagation constant of the T E 01 mode supported by a rectangular waveguide. In addition, ξ = β i /ωµ 0 (i = 2, 3) is the wave admittance of the T E 01 mode inside each hole. Now, applying the boundary conditions at z = 0 for the electric fields (the continuity of E x and E y at every point of the unit cell) leads to the following equations , and which are obtained by multiplying the electric fields to e jkx,mx e jky,ny and taking the integral of over one unit cell. Similarly, we apply the continuity of the tangential magnetic fields (H x and H y ) at z = 0. Using (1),(6a) and (7a), and by multiplying the magnetic fields by sin(π(y − d yi )/b i ) and then taking the integral of both sides over each hole, we have . By combining (8), and (10), and after some straightforward mathematical manipulations, reflection coefficients can be derived as ); m, n = 0 (11d) mn ; m, n = 0 (13c) Finally, the diffraction efficiencies (the ratio of diffracted power to the incident wave) can be calculated by the following relation These calculations can be generalized to the case of a 2D-CMG with more than two holes in each period. Fig. 2 depicts a 2D-CMG composed of N (arbitrary number) holes per unit-cell with the lattice constant P x and P y along x− and y−axes, respectively. We denote the corner of the ith hole by d xi and d yi , its height by h i and width and length by a i , and b i , respectively (Fig. 2). The ith hole is filled with a dielectric material with a refractive index of n i . Similarly, the total electric and magnetic field must be expanded in all regions, and appropriate boundary conditions must be applied to derive the reflection coefficients of the diffracted orders. For brevity, the details of these calculations are not presented here.

B. Numerical results
Here, we present some numerical examples to verify the accuracy of the proposed method. As the first numerical example, in accordance with Fig. 1, we set the parameter of the structure to P y = 0.83P x , .08P x , and n 1 = n 2 = n 3 = 1. For the second example, consider a 2D-CMG with four holes in each period. The parameters of the structure are assumed as P x = 0.77P y , 46P y , n 1 = n 3 = n 4 = n 5 = 1, and n 2 = 1.5. The diffraction efficiencies of the diffracted orders versus the normalized frequency are displayed in Figs. 3 and 4. Here, we define the normalized frequency as ω n = max[P x , P y ]/λ 0 , and λ 0 is the free space wavelength. A full-wave simulation is also carried out to validate the analytical method using the finite integration technique (FIT) in CST Microwave Studio 2019. In CST, periodic boundary conditions are applied in both x-and y-directions, while the perfectly matched layer (PML) boundary condition is applied in the z-direction. Evidently, the results of our proposed analytical method are in excellent agreement with those obtained by using the full-wave simulations.

III. APPLICATIONS AND DISCUSSIONS
In this section, we design anomalous reflectors and beam splitters with near-unitary power efficiency using the proposed analytical method. As mentioned earlier, each term of the summations in (1)-(4) can be interpreted as a plane wave at elevation angle θ mn and azimuth angle ϕ mn which differ from the angles of the incident wave (θ 00 and ϕ 00 ), except for the specular mode (00 order). From (5), θ mn and ϕ mn of the diffracted wave can calculated by following equations cos ϕ mn sin θ mn = k x,m /k 0 (15a) sin ϕ mn sin θ mn = k y,m /k 0 (15b) The five first diffraction orders (i.e. (0, 0), (±1, 0) and (0, ±1)) are depicted in Fig. 5, while the metagrating is illuminated by a normal TM plane wave. In showing how the structure works, we do not depict the higher diffraction orders in Fig. 5 for simplification. In this case, (±1, 0) FB modes lie within the x−z plane at angle θ ±1,0 from the z-axis, and (0, ±1) FB modes lie within the y−z plane at angle θ 0,±1 from the z-axis as shown in Fig. 5. Moreover, the (0,0) mode overlaps at z-axis. Based on Equations (15) and (5), θ 0,±1 , and θ ±1,0 can be expressed and ϕ ±1,0 = 0 • and ϕ 0,±1 = 90 • for the normal incidence. Note that the azimuth angle of a higher diffraction order can take an arbitrary value due to periodicity and operating wavelength. Therefore, we have a multi-channel metagrating, with each channel having a certain elevation and azimuth angle proportional to the period of the structure and the wavelength. In the following, we aim to use a 2D-CMG to manipulate the power distribution between these channels of this metagrating, achieve the desired diffraction pattern, and propose various applications accordingly.

A. Perfect out-of-plane anomalous reflection
In this subsection, we design a perfect out-ofplane anomalous reflector using the proposed analytical method. The problem geometry is depicted in Fig. 6(a). Our goal is to couple the normal incident TM plane wave, to a TE plane wave in the y − z plane with angle θ 0,−1 from the z-axis. Note that, (0, ±1) FB mode are propagation along a direction in the y − z plane which dose not lie in the plane of incidence (x − z plane). Hence, if the power of incident wave transfer to (0, ±1) FB mode, out-of-plane reflection can be realized.
To simplify the design process, we assume that only the (0, 0) and (0, ±1) are propagating while higherorder diffracted modes are evanescent, which is achieved by choosing P x ,and P y from the range of [0, λ 0 ] and [λ 0 , 2λ 0 ], respectively, and satisfying the following condition : k 2 x,1 + k 2 y,1 < k 2 0 . By eliminating the DE 0,0 and DE T E 0,1 , we achieve unitary efficiency for the (0, −1) mode since the higher-order modes are evanescent and the entire structure is lossless. Based on [34,48], for perfect elimination of N FB modes, N meta-atoms are needed for the structure of 1D metagrating. According to [38], this principle is verified for 2D metagratings. Here, each hole is considered as a meta-atom; therefore, to suppress two FB modes, we use a 2D-CMG with two holes per period.
For a −50 • deflection angle, based on (16b), the periodicity of the structure along the y-axis must be chosen as 1.3λ 0 . To further simplify the fabrication process, we assume that all holes are filled with air (n i = 1). To extract the other parameters, we utilize the genetic algorithm (GA) to minimize the DE 0,0 and DE T E 0,1 of the structure. Using the proposed method, we define the cost function as DE 0,0 + DE T E 0,1 + 1/DE T E 0,−1 in the desired frequency. The optimized parameters of the structure are extracted as P x = 0.59P y , a 2 = 0.083P y , a 3 = 0.073P y , b 2 = 0.637P y , b 3 = 0.4P y , h 2 = 0.669P y , h 3 = 0.427P y , d x = 0.46P y , and d y = 0.288P y . The diffraction efficiencies of the optimized metagratings are plotted in Fig. 6(b), depicting an excellent agreement between the results of full-wave simulation and those predicted by our analytical approach. It can be seen in Fig. 6(b) that almost all the power of the incident wave (99.9%) is transferred to the (0, −1) order in the desired frequency ω n = 1.305. This efficiency is a remarkable achievement compared with previously reported anomalous reflectors [37][38][39][40][41]. The magnitude of the electric field distributions is also depicted in Figs. 6(c) and 6(d). Based on the electric field distributions, the designed metagratings transfer a normal incident TM plane wave to an oblique TE plane wave (with an angle of θ 0,−1 = −50 • ) in the y − z plane.
Moreover, the 2D-CMG can be used for designing outof-plane reflectors with different design angles. Similarly, we repeat the design process for extracting parameters of the anomalous reflector with θ 0,−1 in the range of −35 • to −75 • . The optimized structure parameters and the power efficiency (PE) of the designed metagratings are listed in Table I. In all of the designed anomalous reflectors, we can achieve near-unitary efficiency. It should be noted that according to (5a), the anomalous reflection occurs in the normalized frequency ω n = P y /λ 0 . Note that the 2D-CMG can also be used to realize the inplane anomalous reflectors (coupling the incident power to (±1, 0), (±2, 0), ... FB modes ). Nevertheless, there are some structures with less complication, such as 1D-CMG, that can realize this phenomenon [35]. For brevity, we do not present the results of this application here.

B. Five-channel beam splitters
Herein, we design several five-channel beam splitters with an arbitrary power distribution based on the concept of metagratings using analytical expressions derived in the previous section. The new problem geometry is  , and (f) The diffraction efficiencies of the designed beam splitters. The incident power can be arbitrarily distributed between diffraction orders by the 2D-metagratings. The optimum parameters of the structure are mentioned in the text and also listed in Table II, and Table III. Based on the given design angle, beam splitting occurs in the normalized frequency ωn = 1.115, ωn = 1.316, ωn = 1.414 for the first, second, and third beam splitters, respectively. illustrated in Fig. 7. To simplify the design process, we assume that only the first five diffraction orders ((0, 0), (0, ±1), and (±1, 0)) are propagating, and higher orders do not carry any power in region 1. Consequently, we restrict the periodicities of the structure and the operating frequency to the range that satisfies these conditions: k 2 x,1 + k 2 y,1 < k 2 0 , k 2 y,2 < k 2 0 , and k 2 x,2 < k 2 0 . By distributing the incident power between these diffraction orders, a five-channel beam splitter can be realized. Here, we note again that two channels (orders) of the proposed beam splitter lie in the x − z plane(±1, 0), two of them lie in the y − z plane (0, ±1), oriented along the angle θ ±1,0 and θ 0,±1 , respectively, from the z-axis, and these angles can be controlled by changing the periodicities. Furthermore, the final channel (0, 0) overlaps on the z-axis. To attain desired power distribution and complete control over diffraction patterns, four meta-atoms (holes) per period can provide sufficient degrees of freedom, based on what was outlined in [34,48] and due to the passivity condition (Note that the proposed structure is lossless).
As proof of the concept, we designed three devices for five-channel beam splitting using the proposed analytical method. All holes of the designed beam splitters in this subsection are filled with air for more simplicity in the fabrication process (n i = 1). The first beam splitter suppresses the (0, 0) FB mode and divides the incident power uniformly into four directions: two of which lie in the x−z plane, and the others lie in the y − z plane. All directions have an angle of 60 • to the z-axis (|θ ±1,0| =|θ 0,±1 |=60 • ). Since this beam splitter has a symmetric diffraction pattern (DE T E 0,1 = DE T E 0,−1 and DE T M 1,0 = DE T M −1,0 ), for a simpler design, we can use a metallic grating with one rectangular hole in each unit cell as shown in Fig. 7(a). According to the given angle and (16), periodicities must be chosen as P x = P y = 1.155λ 0 . After running an optimization using the proposed method (here the cost function defined as |DE T E 0,1 − 0.25| 2 +|DE T M 1,0 − 0.25| 2 in the normalized frequency ω n = 1.115), the other parameters of the structure (i.e., the width, the length, and the height of the rectangular hole) are extracted as a 2 = 0.65P x , b 2 = 0.479P x , and h 2 = 0.564P x . The diffraction efficiencies of the optimized structure are plotted in Fig. 7(b). Evidently, power uniformly is transferred to (0, ±1), and (±1, 0) in the desired frequency (ω n = 1.115). The relative distribution error (defined here as relative deviation from the desired power distribution) is less than 1%, and the total diffraction efficiencies of orders are more than 99.9%, which is significantly improved in terms of both the power efficiency and relative distribution error compared with the previously published beam splitters [49][50][51][52][53].
In the following, we design two beam splitters with an asymmetric diffraction pattern, unlike the first beam splitter. Therefore, in these cases, we use a 2D-CMG that has four holes in a period. A schematic representation of the second beam splitter is depicted in Fig. 7(c). This beam splitter eliminates (0, 0) and (1, 0) FB modes while transferring 50% of the incident power to (−1, 0) order, 30% to (0, +1) order, and 20% to (0, −1). In this case, (±1, 0) orders are oriented along directions with θ ±1,0 = ±50 • to the z-axis and (0, ±1) orders along directions with θ 0,±1 = ±65 • to the z-axis. Hence, P x and P y must be chosen as 1.3λ 0 and 1.1λ 0 , respectively, according to (16). The third designed metagrating reflect 10% of the incident power to the specular mode and transfer 30% and 25% of the incident power to the channels lying in the x − z plane with angle θ ±1,0 = ±45 • to the z-axis. The rest of the incident power goes to the y − z plane with a splitting ratio of 3 : 4 and angle θ 0,±1 = ±70 • as shown in Fig. 7(e). To design such a beam splitter, periodicity along the x-and y-axes must be chosen as 1.41λ 0 and 1.064λ 0 , respectively, according to (16). Again, we utilize GA to extract other parameters of the second and third beam splitters using the proposed analytical method (Tables II and III). The used cost function for designing these beam splitters respectively for the second and third cases. The diffraction efficiencies of the designed metagratings are plotted in Figs. 7(d) and 7(f), depicting an excellent agreement between the results of the full-wave simulation and our analytical method. The findings demonstrate a near-unitary total efficiency (99.9%) and a relative distribution error of < 1%, which is a remarkable achievement compared with previously reported grating-and metasurface-based beam splitters [28,34,48,52,[54][55][56][57][58]. Therefore, a fivechannel beam splitter with arbitrary power distribution and a near-unitary efficiency can be realized using the proposed method and based on the concept of metagrating.
Next, we consider the effects of the metal ohmic losses on our proposed devices. As noted earlier, we assume that rectangular holes in 2D-CMG are carved on a PEC slab (not a real metallic slab). This approximation is valid for microwave, millimeter-wave, and low THz regimes. Hence, the designed metagratings can be used in a wide range of frequencies. To investigate this effect, we replace PEC with a lossy metal (copper with a conductivity of σ = 5.96 ×10 7 ) in the second designed beam splitter. We consider operating frequency to be 1THz, and other parameters of the structure can obtain from Table II. We perform a full-wave simulation to plot diffraction efficiencies versus frequency ( Fig. 8(a)). The results do not strictly change, and even in this case, the relative distribution error and the total power efficiency are 2.15%, and 99.9%, respectively, which are better than the previously reported beam splitters. The performance of other designed metagratings with lossy metal is similar, and their results are not presented here for brevity.
Finally, we investigate the diffraction pattern of the designed metagrating when truncated to a finite size. Again, we only investigate the second beam splitter made with the lossy metal. To extract the scattering patterns of the 2D-CMG with a finite size, we perform a 3D simulation using CST Microwave Studio 2019. The physical size of CMG is approximately 3.95mm in the x-direction and 3.33mm in the y-direction (10 × 10 unit cell). The truncated metagrating is under a normal TM-plane wave and the far-field patterns are depicted in Fig. 8(b) at 1THz. The results show that almost no power is transferred to the directions with angles θ = ϕ = 0 • ((00) FB mode) and angles θ = 50 • , ϕ = 0 • ((10) FB mode) as we expected from Fig. 8(a). Also, it can be observed that in the operating frequency, the incident power is split into three desired directions with a predesigned ratio.

IV. CONCLUSION
Herein, a 2D-CMG was proposed for manipulating inplane and out-of-plane EM waves based on the concept of metagratings. An analytical method was introduced for diffraction analysis of 2D-CMGs and verified through some numerical examples, indicating excellent agreement with full-wave simulation results. Closed-form and analytical expressions were also presented for the diffraction efficiency of the diffracted orders. By using the proposed method and without needing a single simulation in the full-wave software, we designed out-of-plane reflectors and five-channel beam splitters. The proposed reflectors transferred a normal TM plane wave to an oblique TE plane wave in the y − z plane with angles above −30 • to the z-axis with unitary power efficiency. The designed beam splitter distributed the incident power to five directions with an arbitrary ratio. The total power efficiency of the proposed beam splitters was above 99.9% and their relative distribution error was less than 1%. This proposed method can pave the way for the analytical design of 2D metagratings with various potential applications for microwave and terahertz wavefront manipulation.