Double-stacked hyperbolic metamaterial waveguide arrays for efficient and broadband terahertz quarter-wave plates

We demonstrate how it is possible to achieve weak dispersion in the phase delay between two orthogonal polarization states by using double-stacked hyperbolic metamaterial (HMM) waveguide arrays. The weak dispersion in the phase delay originates from the different signs of phase delay from the two different HMM waveguide arrays. The condition of dispersion-free phase delay for the transmitted waves has been theoretically derived from the transmission matrix as the propagation characteristic of the HMM waveguide is involved. We further reveal that the designed double-stacked HMM waveguide array can function as an efficient quarter-wave plate that enables the conversion of linearly polarized light to circularly polarized light within a broad frequency band. In addition, the bandwidth over which the degree of linear polarization is nearly unity and over which the angle of linear polarization is kept at approximately 45° is basically consistent with the phase bandwidth. This offers a promising approach for developing a practical polarization converter in the terahertz domain.

metamaterials (HMMs) due to their outstanding electromagnetic properties, including the strong enhancement of spontaneous emission, negative refraction and superlensing effects [28][29][30][31] . It has also been demonstrated that an HMM waveguide can be constructed for wideband photon harvesting 32 and light guiding at the ultra-deep subwavelength scale 33 . Quite recently, we have demonstrated that an HMM waveguide array with a rectangular waveguide cross-section exhibits a giant modal birefringence index 34 , which is dozens of times higher than that of conventional quartz birefringent crystals for THz waves 35 . The birefringence index of the HMM waveguide array is comparable to or even higher than that based on various metamaterials in the THz domain [36][37][38] . More interestingly, the designed polarization manipulation devices with such an HMM waveguide array show the capability of converting linearly polarized light waves to circularly polarized light waves with high transmission. However, their operation bandwidth is very limited, since the HMM polarization manipulation components show strong dispersion in the phase delay between two orthogonal directions of light. Once the device operates away from the optimal wavelength, the phase delay is no longer kept at a constant value, which seriously restricts the HMM polarization manipulation components for broadband applications.
In this article, we have explored double-stacked HMM waveguide arrays to engineer the phase delay for the broadband manipulation of light polarizations. By incorporating the propagation characteristic of the HMM waveguide into the general transmission matrix, we have derived the conditions for achieving dispersion-free phase delay for transmitted waves. In cases where the cross-section parameters of the HMM waveguides are given, the phase delay can be kept constant within a broad frequency range of interest by properly arranging the heights of the HMM waveguides. Numerical simulation results demonstrate that the designed quarter-wave plates enable the conversion of linearly polarized light to circularly polarized light within a wide spectral range. In addition, the double-stacked HMM polarization manipulation components show an excellent figure of merit in terms of the degree of linear polarization (DoLP) and the angle of linear polarization (AoLP) due to the weak dispersion in the amplitude transmission for different polarizations.

Results and Discussion
A single-sized HMM waveguide array. Figure 1(a) is a schematic representation of a unit cell of an HMM waveguide array with a rectangular cross-section, where a subwavelength Al/polymer benzocyclobutene (BCB) multilayer is placed on the low-loss high-density polyethylene (HDPE) substrate and the surrounding dielectric layer is BCB. The relative permittivity of the Al layers is described by the Drude model ε = − The polymer BCB is widely used in THz applications on account of its ultralow loss, and its relative permittivity can be extracted from experimental data 40 [ Fig. 1(b)]. The HDPE is selected as the substrate due to its low refractive index (1.54), high stability, and small absorption coefficient in the THz   43 indicates that the estimated effective dielectric permittivities parallel and perpendicular to the normal direction of the layered structure are positive and negative [ Fig. 1(c)], respectively. Here, the TM (TE) mode of the HMM waveguide array is defined as the electric field oriented parallel to the x (y) direction. It has been demonstrated that such an HMM waveguide array exhibits a giant modal birefringence for the transmitted wave, which can be dozens of times larger than that of the birefringent crystals for THz waves 34 . Consequently, an HMM waveguide array can manipulate light polarizations with a subwavelength thickness. We first consider an HMM waveguide array with a rectangular cross-section [L x1 > L y1 , Fig. 1(a)], which leads to different cut-off frequencies for different polarizations. In this case, the cut-off frequency for the TM mode (1.48 THz) is much smaller than that for the TE mode [ Fig. 1(e)]. It is well known that slow light occurs if the light frequency approaches the cut-off frequency and hence significantly enhances the photonic density of states 44 , which results in strong absorption of the EM wave 32,[45][46][47][48] . Thus, as light frequency approaches the cut-off frequency for the TM mode, the transmission amplitude of the TM mode is highly suppressed due to the slow-light effect, while that for the TE mode remains high [ Fig. 1(e)]. Meanwhile, the phase delay between the TE and TM modes increases sharply near the cut-off frequency for the TM mode [ Fig. 1(f)], originating from the fact that the propagation constant for the TM mode is significantly increased 34 . However, to operate as an efficient and broadband quarter-wave plate, the transmission of large amplitude and weak dispersion in the phase delay between two orthogonal directions of light should be achieved simultaneously. Despite the fact that the amplitude transmission can be kept at a high level provided light frequency is far from the cut-off frequency [ Fig. 1(e)], such an HMM waveguide array suffers severely from the dispersion in the phase delay between the TE and TM modes, which significantly restricts its broadband applications [ Fig. 1 We next consider another HMM waveguide array with a rectangular cross-section [L y2 > L x2 , Fig. 1(d)]. It is interesting to note that a comparable amplitude transmission is maintained [ Fig. 1 Broadband quarter-wave plate based on double-stacked HMM waveguide arrays. An effective quarter-wave plate should be able to transmit light waves of different polarizations with a high and approximately equivalent output amplitude over a wide spectral band. Further, the phase delay between the TE and TM modes should be kept constant at π/2 in the design frequency range. We therefore determine the conditions for achieving the constant phase delay of π/2 within the design frequency range for a double-stacked HMM waveguide array by incorporating the propagation characteristic of the HMM waveguides into the general transmission matrix (see Methods Section). When Eqs. (6 and 8) are satisfied simultaneously, the phase delay will be held at π/2 in the design frequency range. It has been previously demonstrated that the propagation constants for different polarizations are determined by the side dimension of the HMM waveguide 34 , so k 1 , k 2 and Δk 1 , Δk 2 as a function of frequency f can be uniquely retrieved. Thus, there exists only one pair of (H 1 , H 2 ) that satisfies both Eqs. (6) and (8). If the actual height of the complete device deviates from the designed height, the introduced phase delay is no longer kept at π/2. In addition, it can be easily inferred from Eqs. (6 and 8) that the sign of k 1 , k 2 (Δk 1 , Δk 2 ) should be −1. When the side dimensions for the two HMM waveguide arrays are chosen to satisfy |k 1 | = |k 2 |, there exist no solutions to Eqs. (6 and 8). In other words, |k 1 | should be unequal to |k 2 |. When the side dimensions for the upper and lower HMM waveguides are designed such that the operation frequency approaches both cut-off frequencies for the TE (the lower HMM waveguide) and TM (the upper HMM waveguide) modes, both k 1 and k 2 are relatively large; thus, the resultant heights for the upper and lower HMM waveguide, H 1 and H 2 , can be significantly reduced according to Eq. (6). Figure 3(a,b,c) show the dependence of k 1 and k 2 , phase delay and transmission on the light frequency when the side dimensions are designed to produce relatively large k 1 and k 2 . In the operation frequency range from 1 to 1.2 THz, the estimated Δk 1 and Δk 2 are −0.049 and 0.033 μm −1 , and H 1 and H 2 are approximately 115 and 185 μm, respectively. It can be easily observed from Fig. 3(b) that the phase delay can be made around π/2 in the design frequency region. However, such a scheme results in low transmission [ Fig. 3(c)], which can be attributed to the significant mismatch between the propagation constant of the HMM waveguide arrays and the wave number in air and to an increased absorption loss induced by the slow-light effect near the cut-off frequency. When the side dimensions for the upper and lower HMM waveguides are designed such that the operation frequency is far from each cut-off frequency for the TE (lower HMM waveguide) and TM (upper HMM waveguide) modes, both k 1 and k 2 are relatively small; thus, the resultant height of each HMM waveguide would be significantly increased according to Eq. (6). Figure 3(d,e,f) show the dependence of k 1 and k 2 , phase delay and transmission on the light frequency when the side dimensions are designed to induce relatively small k 1 and k 2 . In the operation frequency range from 1 to 1.2 THz, the estimated Δk 1 and Δk 2 are −0.009 and 0.003 μm −1 , and H 1 and H 2 are approximately 400 and 1125 μm, respectively. Although the device benefits from the nearly dispersion-free phase delay of π/2 [ Fig. 3(e)] and the relatively high transmission within the frequency range of interest [ Fig. 3(f)], the overall device height is nearly five times that of the light wavelength.
To effectively decrease the device height without sacrificing the transmission efficiency, the side dimensions need to be selected appropriately such that the operation frequency is away from the cut-off frequency. In this way, we can effectively avoid the high absorption from the HMM waveguides while still keeping the values of k 1 and k 2 at a relatively high level, which is beneficial for the miniaturization of the HMM wave plate. Figure 4(a,e,i) show the dependence of k 1 (k 2 ) on f for three sets of structural parameters for the upper and lower HMM waveguide arrays. It should be noted here that there exists some oscillation in the dispersion curves, which can be attributed to the Fabry-Perot effect caused by the HMM waveguides. In particular, since the propagation constant of the TM mode of the upper HMM waveguide is much larger than that of the TE mode of the lower HMM waveguide, the oscillation range of the k 1 curves is bigger than that of the k 2 curves. For the three sets of the double-stacked HMM waveguide arrays, −Δk 1 /Δk 2 is approximately 3, 4, and 5 around 1.1 THz, respectively. By incorporating these values into Eqs. (6 and 8), we can obtain the heights of the HMM waveguides that enable phase delay of π/2 with the aim to convert linearly polarized light to left-circularly polarized light. It is interesting to see from Fig. 4(b,f,j) that by utilizing the three sets of structural parameters to form the double-stacked HMM waveguide arrays, the phase delay exhibits a nearly flat response within a wide spectral band around 1.1 THz. If

s E E E E
x y x y 2 , respectively 49 . Here, '*' denotes the complex conjugate. It can be observed from Fig. 4(c,g,k) that the bandwidth over which the DoLP is nearly unity is basically consistent with the phase bandwidth shown in Fig. 4(b,f,j). More interestingly, the AoLP within the phase bandwidth where the DoLP is near unity is kept at approximately 45°. This insensitivity of AoLP to light frequency within the phase bandwidth originates from the weak dispersion in the amplitude transmission for different polarizations [ Fig. 4(d,h,l)]. The weak frequency dependence of the AoLP has the advantage of fixing the fast and slow axes, which is highly desirable when implementing a quarter-wave plate in a practical situation. We have noted in recent studies on quarter-wave plates based on plasmonic metasurfaces that the AoLP changes significantly within the phase bandwidth of operation owing to the strong dispersion in the transmission/reflection coefficients 11,17,24,26 . Figure 4(d,h,l) present the transmission spectra for the three sets of HMM quarter-wave plates, where linearly polarized light with the polarization angle of 45° relative to the x direction illuminates along the -z direction. It gets as high as 0.5 within the operation bandwidth, which is comparable to or even higher than that in Fig. 3(f). The reason that the transmission is limited can be explained as follows. First, there is a mismatch between the propagation constant of the HMM waveguide arrays and the wave number in air; hence, a portion of light will be reflected back into the air. Second, to effectively reduce the device height, an operation frequency band that is near the cut-off frequency has been selected for polarization manipulation, which comes at the cost of significantly increased absorption loss from the HMM waveguide arrays.
For an HMM waveguide, the mode propagation constant is highly dependent on the cross-section parameters. Therefore, the tolerance of the structural parameters of the HMM waveguide cross-section is very important for the design of HMM quarter-wave plates. Figure 5 shows the tolerance of the phase delay and transmission coefficients as a function of the cross-section parameters for the third set of the HMM quarter-wave plates. It can be observed that the phase delay is kept at around π/2 in the design frequency band even if the widths of the HMM waveguides change significantly [ Fig. 5(a,b)]. Meanwhile, the corresponding amplitude transmission for different polarizations changes slightly in the frequency range of interest [ Fig. 5(c-f)], which is beneficial for constructing a quarter-wave with high performance in terms of AoLP and DoLP. It is worth noting here that the variation in L x1 and L y2 has a much stronger influence on the phase delay than that of L y1 and L x2 . This is because the TM (TE) mode for the upper (lower) HMM waveguide is closer to the cut-off frequency and hence more sensitive to L x1 (L y2 ).
We have also numerically studied the influence of misalignment between the upper and lower HMM waveguide arrays on the device performance. The dislocations along the +x and +y directions are denoted as d x and d y [ Fig. 6(a)], respectively. The results show that the phase delay is slightly changed within the operation frequency range of 0.89-1.22 THz even if the misalignment is up to 4 μm [ Fig. 6(b)]. Meanwhile, the amplitude transmission for different polarizations is accordingly kept high compared to the case without misalignment [ Fig. 6(c,d)]. The small influence of misalignment on the device performance may be very suitable for fabricating a practical double-stacked HMM quarter-wave plate since it means that alignment between the upper and lower HMM waveguide arrays is not strictly required. We have also investigated the case of arranging a dielectric substrate for the upper HMM waveguide array so that the upper and lower HMM waveguide arrays are designed separately. The results demonstrate that the device performance can still be maintained (not shown here), indicating that the double-stacked HMM quarter-wave plate might be experimentally implemented by separately fabricating upper  Fig. 4. The phase delay (a,b), amplitude of transmission coefficients for different polarizations, t xx (c,d) and t yy (e,f). In (a,c,e), ΔL x1 and ΔL y1 denote the variation with respect to L x1 (=62 μm) and L y1 (=20 μm), respectively, while L x2 , L y2 are fixed at 23 μm and 47 μm, respectively. In (b,d,f), ΔL x2 and ΔL y2 denote the variation with respect to L x2 (=23 μm) and L y2 (=47 μm), respectively, while L x1 , L y1 are fixed at 62 μm and 20 μm, respectively. and lower HMM waveguide arrays. This may greatly facilitate the design and fabrication process in a practical situation.

Conclusion
In conclusion, we have explored double-stacked HMM waveguide arrays to engineer the phase delay dispersion for broadband manipulation of light polarizations. The conditions for dispersion-free phase delay have been theoretically derived by incorporating the propagation characteristic of the HMM waveguide into the transmission matrix. The designed double-stacked HMM waveguide arrays are shown to be capable of maintaining the phase delay at around π/2 over a wide spectral band, which enables the conversion of linearly polarized light to circularly polarized light. Finally, the double-stacked HMM waveguide quarter-wave plates show excellent performance in terms of DoLP and AoLP, which may facilitate the implementation of a practical polarization converter in the terahertz domain. We emphasize that the presented design approach might be further extended to address the bandwidth issue for various metamaterials-based polarization manipulation devices [36][37][38][50][51][52] .

Methods
For the double-stacked HMM waveguide array [ Fig. 2], the transmission matrix for describing the complex amplitudes of transmitted waves can be retrieved by taking into account the propagation characteristic of the HMM waveguide. Assuming a normal incident plane wave illuminates the double-stacked HMM waveguide arrays from the -z direction with the polarization angle θ relative to the x-axis, the electric field can be expressed as where E x and E y denote the complex amplitudes of electric field along x and y directions, respectively. The input and output electric field for the double-stacked HMM waveguide arrays can be related using the and t xy1 (t xy2 ), t yx1 (t yx2 ) represent the transmission coefficients for x-polarized (y-polarized) electric field with y-polarized (x-polarized) light illumination for the upper (lower) HMM waveguide array, respectively. k x1 (k x2 ), k y1 (k y2 ) denote the real part of the wave vectors of x-and y-polarization for the upper (lower) HMM waveguide array, respectively. Since an HMM waveguide array doesn't convert x-polarized (y-polarized) electric field to y-polarized (x-polarized) electric field, t xy1 , t yx1 , t xy2 and t yx2 should be equal to zero. Consequently, the Jones matrix t for the double-stacked HMM waveguide arrays can be rewritten as Here m = 1 (−1) denotes LCP (RCP). It is worth noting here, an HMM waveguide array has the ability to transmit light waves of different polarizations with a high and approximately equivalent output amplitude over a wide spectral band provided light frequency is far from the cut-off frequency [ Fig. 1(e,g)]. Therefore, according to Eq. (5) it is highly expected that the polarization angle θ can be constantly kept at π/4 over a wide frequency band. The left side of Eq. (6) is the summation of the phase delay caused by the upper and lower HMM waveguide arrays. For simplicity, we rewrite ϕ 1 = k 1 H 1 , ϕ 2 = k 2 H 2 , where k 1 = k y1 − k x1 , and k 2 = k y2 − k x2 .
To ensure that the phase delay shows a flat response, the derivative of left-side of Eq. (6) should be equal to zero in a wide spectral band of interest δ δ where Δk 1 = k 1 (f 2 ) − k 1 (f 1 ), Δk 2 = k 2 (f 2 ) − k 2 (f 1 ). Once the structural parameters of the HMM waveguide cross-sections are given, k 1 and k 2 as a function of f will be determined. As a result, according to Eqs. (6 and 8) H 1 and H 2 can be uniquely retrieved.
In this work, the propagation characteristic of the HMM waveguide array is estimated by numerical simulations with Lumerical finite difference time domain solutions. Periodic boundary condition is employed in the x and y directions, and perfect matched layer absorption condition is applied in the z direction. 5 and 20 mesh grids have been used to represent the thicknesses of Al and BCB layers, respectively, and 160 mesh grids are set in the x and y directions. We have also conducted simulations with even finer grid sizes, and the resultant propagation characteristic of the HMM waveguide array is almost unchanged, indicating that the grid sizes applied are sufficiently enough to accurately retrieve the propagation characteristic of the presented HMM waveguide.