Introduction

Wave plates are common devices for adjusting the polarization state of an incident wave and are often used for phase delay1, polarization rotation2,3, holographic imaging4, sensing5,6and other functions7,8. The most commonly used wave plates are the half wave plate (HWP) and the quarter wave plate (QWP). The HWP is capable of converting the direction of polarization of an incident wave into an orthogonal direction, whereas the QWP is able to transform a linearly polarized wave into a circularly polarized one. Conventional wave plates are typically composed of natural birefringent materials, such as calcite, which exhibit different refractive indices in orthogonal directions and induce disparate phase changes when light is transmitted9. Terahertz (THz) electromagnetic wave has a frequency within the range of 0.1 to 10 THz. It is associated with millimeter waves in the long wavelength band and infrared light in the short wavelength band. The distinctive spectral characteristics of the THz wave, including its low energy, high penetration, and fingerprint spectrum, have made it a valuable tool in a multitude of fields, including radar, remote sensing, non-destructive detection, biomedicine and aerospace10,11,12,13,14. The key to promoting the application of THz waves is to better control their polarization. However, the lack of natural materials with a strong birefringence effect in the THz band has resulted in the traditional wave plates having the disadvantages of a large size and low efficiency. These factors have constrained the further development of THz technology.

Metasurfaces are two-dimensional simplified metamaterials that are able to achieve special electromagnetic responses by artificially arranging periodic structures. With the increasing research on metasurfaces, a variety of metasurface wave plates have been designed to achieve the modulation of THz waves15,16,17. Metal is the most common material for making metasurface wave plates, Grady et al.18 periodically arranged arrays of metal cut wires on the dielectric layer can trigger electromagnetic resonance that can make the phase difference of the excited scattered light, and achieve the function of HWP in the range of 0.8–1.36 THz. Ako et al.19 designed a metal grating-dielectric-“H” shaped metal resonator-dielectric-metal grating structure of a metasurface wave plate, which is able to form a Fabry-Perot resonance cavity in two dielectric layers when THz waves are incident, achieving efficient linear polarization conversion with ultrabroadband in the range of 0.2–1.0 THz. However, metallic materials inevitably have the limitation of high ohmic losses, and in order to ameliorate this problem, all-dielectric metasurface wave plates have also begun to receive attention. Zi et al.20 created an all-silicon metasurface wave plate by gouging an elliptical cylindrical air slot into a silicon cube, which resulted in a strong birefringence. The metasurface wave plate can be used as a perfect QWP at 1.76 THz. Wang et al.21 designed an all-silicon metasurface wave plate in which they induced multiple Mie resonances through silicon blocks, and by tuning the multilevel interference between these resonances, the wave plate is able to be used as a QWP at 0.79 THz and as a HWP at 0.91 THz, with insertion losses of 0.54 and 1.25 dB, respectively. With the research on tunable materials, such as \(\hbox {VO}_2\), graphene, GST, Dirac semimetal, etc., multifunctional metasurface wave plates based on these materials have been proposed22,23,24,25. Zhang et al.26 fabricated a multifunctional metasurface wave plate using \(\hbox {VO}_2\) and Dirac semimetal. The electrical conductivity of \(\hbox {VO}_2\) can be altered by adjusting the temperature, and when the \(\hbox {VO}_2\) is in the dielectric state, the wave plate can be used as a perfect QWP in the asymmetric transmission in the range of 2.023–5.971 THz, and a reflective QWP in the range of 2.058–3.423 THz, 4.753–5.6 THz and 6.496 THz when \(\hbox {VO}_2\) presents a metallic state. In each of these modes of operation, a tunable spectral response was achieved by modifying the Fermi energy of the DSM material. Yuan et al.27 designed and fabricated a graphene metasurface wave plate that enables broadband linear polarization transitions and linear-to-circular polarization transitions by tuning the Fermi energy levels of graphene. The above THz metasurface wave plates are able to achieve the modulation of THz waves, but still have lower efficiency, narrower bandwidth, larger dispersion loss or other problems, which are still the basic problems and challenges in the design of THz metasurface wave plates.

Grating is a common artificial anisotropic structure, which can achieve the birefringence effect by designing the stripe period and duty cycle. However, the traditional grating wave plates have problems such as narrow bandwidth and large size. Gratings can be formed as metasurface by using grating stripes as unit structures and incorporating rotation angles into the grating stripes within each unit structure. In our previous studies, it has been demonstrated that a single-layer grating metasurface wave plate has a larger dispersion bandwidth than conventional grating wave plate28. In this paper, we have investigated the achromatic bandwidth of multilayer all-dielectric grating metasurface QWP. Unlike previous studies, we chose to base our design and analysis on a fixed intercept period due to the possibility of period mismatch in the surface QWP of the multilayer grating element (the difference in period between two layers may be too large), which can lead to waveform leakage. We find that intercepting each layer of gratings according to a fixed intercept length increases the boundary mismatch of each cell structure, which leads to a stronger interferential coupling of THz waves. Meanwhile, the simulation results show that the polarization conversion ability of the trilayer silicon grating metasurface QWP originates from the anisotropy of each layer structure, while the broadband achromatic performance originates from the wave vector coupling between the unit cells and the resonance effect between the grating layers. This paper improves the problems that the current wave plates need more layers to obtain larger achromatic bandwidth and lower transmittance after the combination of layers, points out that a reasonable setting of the intercept period can produce more resonance effects and thus increase the bandwidth, which provides a theoretical basis for the design of ultrabroadband achromatic multilayer wave plates in the future.

Algorithm principle and theoretical analysis

In this paper, we focus on the design of a multilayer THz metasurface grating QWP. In order to achieve broadband achromatic polarization conversion, we need to keep the phase difference between the fast and slow axes of the QWP within \(90^{\circ }\) ± \(3^{\circ }\), which requires considering the input parameters of the grating structure for each layer: (period: P, duty cycle: \(\sigma\) , depth: D, and rotation angle: \(\theta\) ), and at the same time optimizing the two outputs (phase difference and ellipticity).

Fig. 1
figure 1

Schematic diagram of a conventional trilayer grating with silicon as the grating material, when the incident wave is a \(45^{\circ }\) linear polarization wave, the outgoing wave is a circular polarization wave. \(\hbox {P}_j\), \(\hbox {D}_j\) and shown in the figure, are the period, depth and rotation angle of the j layer (j = 1,2,3).

Full size image

Firstly, we have optimized the structural parameters of the conventional grating QWP by theoretical calculations. Silicon is a common dielectric material with low absorption in the THz range, therefore, is a candidate for making low-loss highly birefringent wave plates. The conventional trilayer grating is shown in Fig. 1. It can be found that the conventional multilayer grating just uses each layer as a single-frequency wave plate, and the bandwidth is increased by stacking. When the grating is not rotating and the wavelength of the incident wave is similar to the grating period of each layer or the period \(\hbox {P}_j\) and wavelength \(\lambda\) of each grating layer satisfy \(\frac{P_j}{\lambda }<\frac{1}{n_1}\)29, where \(\hbox {n}_1\) = 3.415, \(\hbox {n}_2\)= 1. The birefringence in the TE and TM directions can be calculated by quasi-static effective medium theory (EMT) theory30:

$$\begin{aligned} \left\{ \begin{array}{l} n_{TE,j}^{2}=\sigma _jn_{1}^{2}+\left( 1-\sigma _j \right) n_{2}^{2}\\ n_{TM,j}^{2}=\left[ \left( \frac{\sigma _j}{n_1} \right) +\left( \frac{1-\sigma _j}{n_2} \right) \right] ^{-1}\\ n_{TE2,j}^{2}=n_{TE}^{2}+\frac{1}{3}\left[ \frac{P_j}{\lambda }\pi \sigma _j\left( 1-\sigma _j \right) \left( n_{1}^{2}-n_{2}^{2} \right) \right] ^2\\ n_{TM2,j}^{2}=n_{TM}^{2}+\frac{1}{3}\left[ \frac{P_j}{\lambda }\pi \sigma _j\left( 1-\sigma _j \right) \left( \frac{1}{n_{1}^{2}}-\frac{1}{n_{2}^{2}} \right) n_{TE}n_{TM}^{2} \right] ^2\\ \end{array} \right. \end{aligned}$$
(1)

From the \(n_{TE2,j}\) and \(n_{TM2,j}\) obtained in the Eq. 1, we can find out the phase difference(PD) obtained without rotation of the wave plate at that frequency:

$$\begin{aligned} \begin{aligned} PD\ =\ 2\pi \frac{\left( n_{TE2,j} -n_{TE2,j} \right) D_j}{\lambda }=\frac{2\pi }{c}\cdot D_j\cdot \left( n_{TE2,j} -n_{TE2,j} \right) \cdot f \end{aligned} \end{aligned}$$
(2)

It is well known that the key condition to satisfy that the transmitted wave is a circular polarization wave is that the TE direction and the TM direction have a phase difference of \(90^{\circ }\). In order to calculate the phase difference after rotation, we introduce the concept of Jones matrix. A THz wave with polarization can be expressed as:

$$\begin{aligned} T=\left[ \begin{array}{c} \cos a\\ e^{ib}\sin a\\ \end{array} \right] \end{aligned}$$
(3)

Where b is the phase difference between the TE, TM components of the wave, and a denotes the angle between the polarization direction and the horizontal direction of the THz wave, when \(b=\left( 2k+1 \right) \frac{\pi }{2}\) and \(a=\left( 2k+1 \right) \frac{\pi }{4}\), the wave is circularly polarized wave. In Fig. 1, the direction along the grating stripes is the fast axis (u-axis) of the grating, while the direction perpendicular to the grating stripes is the slow axis (v-axis) of the grating, and the angle between the fast axis and the x-axis is the angle of rotation \(\theta _j\), and then we can derive the relationship between the fast and slow axes and the x-y axes, as shown in the Eq. 4.

$$\begin{aligned} \left[ \begin{array}{c} u\\ v\\ \end{array} \right] =R\left( \theta _j \right) \left[ \begin{array}{c} x\\ y\\ \end{array} \right] =\left[ \begin{matrix} \cos \theta _j& \sin \theta _j\\ -\sin \theta _j& \cos \theta _j\\ \end{matrix} \right] \left[ \begin{array}{c} x\\ y\\ \end{array} \right] \end{aligned}$$
(4)

R(\(\theta _j\)) in the Eq. 4 denotes the rotation matrix. The transmission matrix of grating metasurface wave plate with subwavelength anisotropy is determined by the structure of the unit cell, which operates in an effective dielectric, resulting in birefringence. According to Eq. 4, it can be seen that the axis of the unit structure can be divided into the u-axis and the v-axis, which will have different phase delays \(e^{i\phi _{u,j}}\) and \(e^{i\phi _{v,j}}\)31. We can then derive the Jones matrix of the wave plate at this frequency32:

$$\begin{aligned} J_{wp,j}=R\left( -\theta _j \right) \left[ \begin{matrix} t_ue^{i\phi _{u,j}}& 0\\ 0& t_ve^{i\phi _{v,j}}\\ \end{matrix} \right] R\left( \theta _j \right) \end{aligned}$$
(5)

Where \(t_u\) and \(t_v\) are denoted as the transmittance in the u-axis and v-axis. We set \(\varphi _j=\frac{\phi _{u,j}+\phi _{v,j}}{2}\), \(\delta _j=\phi _{u,j}-\phi _{v,j}\), \(t_u=t_v=1\), then we can further derive:

$$\begin{aligned} \begin{aligned} &J_{wp,j}\ =e^{i\varphi _j}\left[ \begin{matrix} \cos \frac{\delta _j}{2}+i\sin \frac{\delta _j}{2}\cos 2\theta _j& i\sin \frac{\delta _j}{2}\sin 2\theta _j\\ i\sin \frac{\delta _j}{2}\sin 2\theta _j& \cos \frac{\delta _j}{2}-i\sin \frac{\delta _j}{2}\cos 2\theta _j\\ \end{matrix} \right] \end{aligned} \end{aligned}$$
(6)

where \(\delta _j\) is the phase difference formed between the fast and slow axes of the jth layer grating wave plate, \(\theta _j\) is the rotation angle of the jth layer grating wave plate, and \(\varphi _j\) is the average of the phase delays caused by the fast and slow axes of the jth layer grating wave plate. When the grating wave plate is stacked n layers, we can get the overall expression as:

$$\begin{aligned} J_{sum}=\prod _{j=1}^n{J_{wp,j}} \end{aligned}$$
(7)

Since the incident THz wave is a \(45^{\circ }\) linear polarization wave, which can be expressed as a = \(\frac{\pi }{4}\) and b = 0 in Eq. 3, then the matrix of the incident THz wave can be derived as: \(p_{in}=\left[ \begin{array}{c} \frac{1}{\sqrt{2}}\\ \frac{1}{\sqrt{2}}\\ \end{array} \right]\). And the final transmitted THz wave can be given as: \(p_{out}=J_{sum}\cdot p_{in}\). Then, based on the Jones matrix of the outgoing wave, we can find its phase difference. When designing a QWP, the two most important factors are the transmittance and the phase difference. Since the silicon grating structure used does not cause an excessive transmittance difference in the TE and TM directions, it is required to design the structure to satisfy the phase difference of \(90^{\circ }\) i.e., the phase in the Y-direction lags behind that in the X-direction by about \(90^{\circ }\)as much as possible over a larger bandwidth. The accurate design of the parameters of each grating wave plate (period, duty cycle, depth, rotation angle) is the key to achieving the ultrabroadband, as this process involves parameters proportional to the number of layers, so the normal optimization methods take up a lot of resources. We propose a nested theoretical-simulation optimization process because the theoretical calculations are ideal, the Fabry-Perot effect and the presence of electromagnetic resonances between the layers, moreover, the grating metasurface wave plate we designed and simulated may have a different polarization conversion performance from the theoretically calculated conventional grating wave plate and the resulting structural parameters do not always have the maximum bandwidth. The proposed joint theory-simulation process is 1). Firstly, the theoretical structural parameters are optimized using a common simulated annealing algorithm to derive a theoretical optimum, which, it is worth noting, is usually not the maximum bandwidth. 2). The optimized theoretical structural parameters are then substituted as initial values into the Simulated Annealing-based Particle Swarm Optimization (SA-PSO) algorithm, which is looped through commercial mathematical software and electromagnetic simulation software based on finite element method until an optimized solution with maximum simulation bandwidth is found33. We have used the same loss function in these processes:

$$\begin{aligned} F_{loss}\ =\ \sum _{f=f_a}^{f_b}{\left[ w_{PD}\left[ PD\left( f \right) -90 \right] +w_{EP}\left[ EP \left( f \right) -1 \right] \right] } \end{aligned}$$
(8)

In the Eq. 8, \(\hbox {W}_{PD}\) and \(\hbox {W}_{EP}\) are the weighting coefficients of phase difference and ellipticity, respectively, PD(f) and EP(f) represent the phase difference and ellipticity of the multilayer wave plate at the frequency f, respectively. The ellipticity EP can be calculated by the following equation:

$$\begin{aligned} EP\left( f \right) =\frac{2\left| E_x \right| \left| E_y \right| \sin \left[ PD\left( f \right) \right] }{\left| E_x \right| ^2+\left| E_y \right| ^2} \end{aligned}$$
(9)
Fig. 2
figure 2

Comparison of our designed nested optimization process with the conventional optimization process in electromagnetic simulation software.

Full size image

Where \(\hbox {E}_X\) and \(\hbox {E}_Y\) denote the X-direction component and Y-direction component of the transmitted electric field, respectively. In the optimization process of 1) since there is no way to derive the transmittance electric field, we assume that the x and y components of the transmittance electric field are the same, i.e., only the phase difference is considered, and in the optimization process of 2) the information of the electric field can be derived from the electromagnetic simulation software, and at this time, all the performance indexes of designing a QWP are considered. We then compared the simulation process using the optimization algorithm directly and the simulation with our proposed nested process. We put the interception period as 50um, the simulation frequency is set at 0.1–1.8 THz, the loss function is used as the basis for judging, the achromatic bandwidth is used as the output, and the final optimization result is shown in Fig. 2.

It can be clearly seen that in the case of substituting the same initial values, although all the algorithms used are SA-PSO, our nested process reaches convergence faster than the conventional process, which indicates that the nested process we designed can reduce the time of searching for the optimized value of optimization directly with electromagnetic simulation software, and give a theoretically calculated optimal solution for the phase difference first, and then optimize the other parameters of the grating metasurface wave plate on this basis, which can effectively solve the current problem of the difficulty in designing the structure of the multilayer achromatic wave plates.

Simulation results and discussions

In the previous section, we have optimized the structural parameters of each layer of a multilayer grating metasurface QWP using a nested process. In this section, we will analyze the performance of the optimized multilayer grating metasurface QWP using electromagnetic simulation software, and simulations have been carried out using a time-domain solver and a frequency-domain solver, respectively. Unlike the multilayer wave plate formed by the stacking of gratings, the multilayer wave plate we designed consists of a stack of grating metasurface in each layer, as shown in Fig. 3. The grating metasurface of each layer are arrays of gratings as structural units, which enable THz waves to penetrate the wave plate not only with phase difference due to the birefringence of the unit gratings, but also with stronger interference coupling due to the destructive nature of the structure.

Fig. 3
figure 3

Schematic structure of a trilayer THz metasurface grating QWP after optimization of the nesting process. (a). Overall structure. (b). Structure of the first-layer grating metasurface QWP and the unit cell diagram after interception, where IP is the interception period. (c). Structure of the second-layer grating metasurface QWP and the unit cell diagram after interception, where IP is the interception period. (d). Structure of the third-layer grating metasurface QWP and unit cell diagram after interception, where IP is the interception period.

Full size image

Performance analysis of achromatic polarization conversion

We fix the interception period (IP) to 50um, and after the nested process optimization in the previous section, we can derive the structural parameters of each layer as: \(\hbox {P}_1\)=47um, \(\sigma _1\)=0.8, \(\hbox {D}_1\)=46.85um, \(\theta _1=1.9^{\circ }\); \(\hbox {P}_2\)=34.25um, \(\sigma _2\)=0.71, \(\hbox {D}_2\)=48.4um, \(\theta _2=1.5^{\circ }\); \(\hbox {P}_3\)=34.42um, \(\sigma _3\)=0.64, \(\hbox {D}_3\)=12.65um, \(\theta _3\)=-\(8.9^{\circ }\). It is worth noting that in the optimization process of the simulation, all the materials we have used are silicon with a relative permittivity of 11.9. Firstly, we have carried out the simulation analysis based on the optimized parameters using a frequency domain solver, where in the Z-direction we have chosen an open boundary condition, while in the X-, Y-directions we have chosen a unit cell boundary condition. The incident wave with \(45^{\circ }\) linear polarization angle penetrates the QWP along the direction shown in Fig. 3, and the detector is placed at a distance of 1/4 of the central wavelength in the penetration direction, and finally the transmittance in the X and Y directions can be measured as shown in Fig. 4(a). We can find that the transmittance of our designed multilayer grating metasurface QWP lies in the range of 0.8–1 in both X and Y directions. The upward and downward fluctuations of the transmittance curves originate from the Fabry-Perot (F-P) effect caused by the superposition of the thicknesses of the multilayer wave plates, and at the same time, since THz waves have different losses when passing through silicon grating stripes and air, and the combined effect of the two factors may lead to a large difference in the intensity of the transmitted waves in the X- and Y-directions, the ratios of transmittance in the X-direction and Y-direction have also been calculated. As shown by the red line in Fig. 4(a), the ratio of transmittance fluctuates in the range of 0.9 to 1.2, which indicates that our designed wave plate produces a relatively small difference in the transmittance intensity between the directions in the range of 0.1 to 1.8 THz, and meets the conditions of high-performance QWP.

Fig. 4
figure 4

(a). Transmittance of an incident THz wave in the X and Y directions of multilayer grating metasurface QWP and the ratio of the transmittances with frequency-domain solver. (b). Ellipticity and phase difference of transmitted waves, where sky blue colour indicates regions with ellipticity greater than 0.94, and light red colour indicates regions where the phase difference lies between \(87^{\circ }\) and \(93^{\circ }\).

Full size image

In Fig. 4(b) we plot the phase difference and the ellipticity of the trilayer grating metasurface QWP. According to the relationship between the ellipticity and the axis ratio, it can be seen that when the ellipticity is greater than 0.94 the transmitted wave can be regarded as a circular polarization wave34, which is represented in the figure as a sky blue region with a bandwidth from 0.63 THz to 1.8 THz reaching 1.17 THz. We regard \(90^{\circ }\) ± \(3^{\circ }\) as the acceptable range of phase difference between the X-direction and Y-direction, represented as the light red region in figure. The part where these two regions overlap represents the incident THz wave that can be achromatically converted to a left-handed circularly polarized wave, calculated as 0.8–1.8 THz, and the relative achromatic polarization conversion bandwidth is able to reach 76.9%.

Fig. 5
figure 5

(a). Transmittance of an incident THz wave in the X and Y directions of multilayer grating metasurface QWP and the ratio of the transmittances with time-domain solver. (b). Ellipticity and phase difference of transmitted waves, where sky blue colour indicates regions with ellipticity greater than 0.94, and light red colour indicates regions where the phase difference lies between \(87^{\circ }\) and \(93^{\circ }\).

Full size image

Next, we have analyzed the same trilayer grating metasurface QWP using a time-domain solver. The incident wave was still set to 0.1–1.8 THz with a \(45^{\circ }\) polarization angle. The difference is that the time-domain solver does not have unit cell boundaries, so we chosed periodic boundaries in the X- and Y-directions. We still placed the probe plane at a 1/4 wavelength to receive the transmitted wave, and the final performance results obtained are shown in Fig. 5. In Fig. 5(a), the transmittance in the X- and Y-directions of the trilayer grating metasurface QWP obtained by the time-domain solver and the results of the frequency-domain solver show a similar trend in the high-frequency region, but there are some differences in the low-frequency region, and the ratio of the transmittance is also different in the low-frequency region. Overall, there is also no large difference between the transmittance in the X and Y directions obtained by the time domain solver, and the ratio stays between 0.9 and 1.2, which meets the requirements of the QWP. In Fig. 5(b), the phase difference and the ellipticity based on the time-domain solver are shown, and it can be seen that the ellipticity is greater than 0.94 in the range of 0.637 THz to 1.8 THz, which meets the requirement of the circular polarization wave transmittance and is basically the same as that of the frequency-domain solver. Although the phase difference curve is basically similar to that of the frequency domain solver in terms of curve trend, the smooth band of the curve is beyond the \(90^{\circ }\pm 3^{\circ }\), and the achromatic polarization conversion can be achieved only in the range of 1.45 THz to 1.65 THz.

Fig. 6
figure 6

Top view of the grating metasurface wave plate, the cell structure circled by the red box in the figure does not fully satisfy the floquetent condition under the time-domain solver, resulting in some deviation of the time-domain results from the frequency-domain results.

Full size image

The differences in the results of the solvers are attributed to subtle differences in the boundary conditions. In the frequency-domain solver we use the unit cell boundary condition, which allows us to arrange the metasurface’s cellular structures in a periodic array in both the X- and Y-directions simultaneously, whereas in the time-domain solver we use the periodic boundary condition in the X- and Y-directions separately, which only allows us to form the periodic structures in the X- and Y-directions individually, such as the cellular structures framed in Fig. 6, which are inaccuracy in the time-domain. These structures do not fully satisfy the floquent condition at the time-domain solver, which makes a slight difference between the two results directly. Therefore, we should use the more accurate results of the frequency domain solver.

Fig. 7
figure 7

The physical field simulation results of the trilayer metasurface grating when the frequency of the incident wave is 1.3THz, and for a more detailed representation of the details, we only show the front view of the two unit structures. (a), (b) Electric field distribution and the direction of the electric field. (c), (d) Magnetic field distribution and the direction of the magnetic field.

Full size image

Analysis of the physical mechanism

In this section, we have analyzed the electric field and the magnetic field inside the trilayer grating metasurface QWP through the simulation of electromagnetic software to determine the working principle. In Fig. 7, we draw the electric field diagrams and the magnetic field diagrams of the two cell structures at the central operating frequency (1.3 THz), and it is worth noting that the simulation is carried out in an infinitely extended period even though we have drawn only two unit cells.

According to the electric field diagram and the electric field direction diagram we can find that most of the electric field energy is gathered at the connection between layers and a toroidal current is formed at the connection, which suggests that the grating metasurface QWP can increase the bandwidth of the operation and the stability of the phase variation by increasing the number of layers and forming a toroidal current due to the nano-resonance at the connection between each layer. Based on the magnetic field diagram and the magnetic field direction diagram, it is clear that, unlike the electric field, the magnetic field is mostly clustered at the connection of neighbouring units in the same layer. This suggests that rotating the grating as a unit structure increases the mismatch at the boundaries of the unit structure (breaks in the grating stripes may occur at the boundaries), and that this broken periodicity can increase the interferential coupling of the THz waves, where enhanced magnetic resonances are formed. It should still be noted that the grating cell with birefringence effect is the fundamental reason for the circular polarization conversion of our QWP, while the enhanced magnetic resonance due to the formation of the grating as a metasurface structure and the interlayer loop current due to the stacked wave plates are the physical principles that ensure the achromatic and the stable broadband operation of the QWP.

Structural parameter analysis

The variation of the structural parameters of the trilayer grating metasurface QWP can influence the birefringence effect and thus achieve different polarization modulation of the incident THz wave, which has already been concluded by numerical analysis in section 2. In this section, the concept of transmission efficiency is also introduced35. This can be defined as \(\rho =\frac{t_{xx}^{2}+t_{yy}^{2}}{2}\), and is used to explore the effect of variations in these parameters on the QWP transmittance. The intercept period (IP) is a value that we have defined to force the period of the unit structure to be the same for each layer of the grating metasurface, and varying this value can lead to a change in the mismatch of the grating structure. We have analyzed this parameter as shown in Fig. 8.

Fig. 8
figure 8

Relationship between the intercept period (IP) and the polarization conversion performance of the trilayer grating metasurface wave plate, where the IP is varied in steps of 10 um starting from 20 um. (a). Variation of transmission efficiency. (b). Variation of phase difference, where the transparent red region indicates the acceptable range of phase difference. (c). Variation of ellipticity, where the transparent red region indicates the interval greater than 0.94.

Full size image

Clearly, when the IP is changed, The transmission efficiency of the QWP changes, the phase difference of the transmitted wave and the ellipticity also change. This is the same as our analysis of the physical mechanism, where a change in IP causes a change in the mismatch of the boundaries, which affects the enhanced magnetic resonance that occurs at the boundaries of each layer. This suggests that we can design the value of IP to achieve modulation of the magnetic resonance, and hence the transmission efficiency and polarization conversion performance. In the actual simulation and fabrication, the interception period should be suitably selected so that more resonances are generated at the appropriate frequency.

Fig. 9
figure 9

Relationship between the thickness of the substrate and the polarization conversion performance of the trilayer grating metasurface wave plate, where the thickness is varied in steps of 2.5 um starting from 2.5 um. (a). Variation of transmission efficiency. (b). Variation of phase difference, where the transparent red region indicates the acceptable range of phase difference. (c). Variation of ellipticity, where the transparent red region indicates the interval greater than 0.94.

Full size image

Next, we have analyzed the thickness of the substrate, in the above, the substrate thickness has been kept at 5um, followed by a certain fluctuation scan. In Fig. 9, it can be seen that by changing the thickness of the substrate, The transmission efficiency of the QWP has not changed much but a change in the polarization of the transmitted THz wave can be achieved. Because the electric field and toroidal current are mostly concentrated in the substrate, when the substrate is thin, there is no way to achieve local accumulation of the electric field, resulting in a weak polarization conversion. And increasing the substrate also increases the overall thickness, which produces a large F-P effect and makes the phase difference not able to be stabilized within the tolerance interval. All these problems affect the polarization conversion ability of multilayer grating metasurface QWP, and it is more important to choose appropriate thickness values of the substrate when the number of layers increases.

Finally, we have investigated the polarization conversion capability of the trilayer grating metasurface QWP under oblique incidence conditions. Previous studies have shown that the birefringence effect of the grating metasurface QWP is the fundamental source of the circular polarization conversion capability, while the modulation of the incident wave can be achieved as long as there is a THz wave with a Z-direction component28. We have scanned the incident angle of the incident wave in steps of \(10^{\circ }\), and the performance variations of the trilayer grating metasurface QWP are shown in Fig. 10.

Fig. 10
figure 10

The relationship between the polarization conversion performance of the trilayer grating metasurface wave plate and the incident angle. The oblique incidence angle is varied in \(10^{\circ }\) steps. (a). Variation of transmission efficiency. (b). Variation of phase difference, where the transparent red area indicates the acceptable range of phase difference. (c). Variation of ellipticity, where transparent red areas indicate intervals greater than 0.94.

Full size image

According to the above figure, it can be seen that the phase difference tolerance interval of the trilayer grating metasurface QWP can be maintained at 80% of the positive incidence condition for incidence angles up to \(40^{\circ }\), which indicates that our designed QWP can maintain a large achromatic polarization conversion capability for incidence angles up to \(40^{\circ }\). It is noteworthy that the dispersion bandwidth of our QWP can be maintained at more than 80% of the original bandwidth when the incident angle is \(50^{\circ }\) or less. These results indicate that when the incident wave deviates from normal incidence, the wave vector component in the Z-direction is reduced, although it is still able to form a phase difference at a specific frequency, the achromatic bandwidth may be reduced because of the inability to form a sufficiently enhanced magnetic resonance and the toroidal current.

Comparisons with the QWPs and manufacturing methods

Subsequently, a comparison between the designed trilayer all-silicon grating metasurface QWP and some recently reported devices is presented in Table 1. The Layer in the table indicates the number of layers in the QWPs (each different material is considered a separate layer), the Mode indicates whether the QWPs work in transmission mode (T) or reflection mode (R), and the angular tolerance indicates whether the QWPs are able to maintain their original working efficiency under oblique incidence conditions(NA denotes that the data is not mentioned in the articles). We can see that most of the current QWPs are metal structures, which utilize the resonance effect of metal to achieve the modulation of THz waves, although the ohmic loss of metal materials may lead to the small relative bandwidth of these QWPs. At the same time, most of the metal-structured QWPs require the metal to be attached to the surface of the dielectric layer, which also adds a certain number of layers. Reflective QWPs, although having higher reflectivity and larger bandwidth, suffer from the problem of crosstalk between the incident and outgoing waves, that means the incident wave will affect the outgoing wave to some extent. The trilayer grating metasurface QWP designed in this paper is fabricated from silicon, which reduces the material loss and enables a large achromatic bandwidth while maintaining high transmittance. Meanwhile, the grating metasurface structure with strong birefringence effect is utilized in each layer, which enables our designed QWP to have a small number of layers and also a large incidence angle tolerance. This has implications for the design of future THz achromatic wave plates. However, our all-silicon structure is also more complicated to fabricate than the metal structure.

Table 1 Comparison of the designed trilayer grating metasurface wave plate with other QWPs
Full size table

In fabrication, our current design of processing method is to etch each layer of grating metasurface wave plate in a single layer, and then mechanically combine them. Because of the relative simplicity of the etching of each wave plate, for multilayer grating metasurface wave plate is also only the structure of each layer is processed separately. However, this processing method faces the problem of how to carry out accurate mechanical combination, due to the small size of the wave plate, the mechanical combination process requires optical axis alignment, this process may have a greater impact on the experimental results, in the experiment should always pay attention to the alignment accuracy of each layer of the wave plate. Some newly proposed multilayer processing techniques, such as imprinting, may also reduce the difficulty of fabricating multilayer wave plates.

Conclusion

In this paper, the trilayer grating metasurface QWP is constructed by stacking all-dielectric gratings as the structural unit, determining a good specific value of the interception period, and then arraying to form the metasurface. The trilayer grating metasurface QWP has 12 structural parameters, which will take a lot of time to calculate directly using the optimization algorithm. Therefore, this paper designs a novel nested optimization process, which can quickly and accurately find the optimal parameters. The nested optimized trilayer grating metasurface QWP is able to achieve linear-circular polarization conversion with an ellipticity greater than 0.94 and achromatic (phase difference of \(90^{\circ }\pm 3^{\circ }\)) within 0.8–1.8 THz, and is able to maintain an achromatic bandwidth of more than 80% when the incident angle is less than \(40^{\circ }\). At the same time, the differences that occur between the time-domain solver and the frequency-domain solver are discussed, and it is determined that the frequency-domain solver provides more accurate results in the analysis of metasurface wave plates. We draw these conclusions from the simulations that the polarization conversion capability of the trilayer grating metasurface QWP originates from the anisotropy of the layer structures, while the broadband achromatic performance originates from the enhanced magnetic resonance due to the wavevector coupling between the cell structures and the electric field circulation between the layers. While our designed trilayer grating metasurface QWP has high performance, its design process and mechanism research also have design guidance for future multilayer ultrabroadband achromatic wave plates.