1. Introduction

Terahertz (THz) waves with frequencies ranging from 0.1 to 10 THz have emerged as a new analytical tool and shown tremendous potential in a variety of fields [1,2], such as communications [3], defense [4], waveguide [5], chemistry [6], and biomedicine [7]. In recent years, metasurfaces that consist of periodically arranged resonators have received considerable attention owing to the novel ways they provide for controlling electromagnetic waves and enhancing light-matter interactions [8] in THz band. The sensitivity of THz technology in the field of detection can be greatly improved by using metasurfaces with high quality (Q) factor resonances. This is due to the fact that a resonance having a high Q factor can enhance the interaction between the THz waves and the analyte [9,10]. However, there still exist fundamental obstacles in achieving high Q factor resonances because of the existence of the radiative losses and the nonradiative losses (ohmic losses) in metasurfaces [11,12]. For metallic metasurfaces, the Drude metals (gold, copper, aluminum, et al.) have inevitable ohmic losses. For all-dielectric metasurfaces, the ohmic losses can be suppressed effectively, which makes it possible to achieve ultrahigh Q factor resonances. Although the ohmic losses of all-dielectric metasurfaces are weak, their radiative losses play the dominant role in hindering the further improvement of the Q factor of the involved resonances [13–15].

Recently, reducing radiative losses for a high Q resonance based on the concept of nonradiative states including anapole mode [16,17] and bound states in the continuum (BIC) [18,19] mode has been attracting extensive attention. Anapole mode can be viewed as the completely destructive interference effect of far-field radiation between electric dipole moments and toroidal dipole moments [20]. The resonances with high Q factors of about $3.8\times 10^{6}$ for metallic metasurfaces [21] and of $2.5\times 10^{6}$ for all-dielectric metasurfaces [22] were successfully achieved by exploiting the excitation of anapole mode. The BIC mode occurs in the country as a result of destructive interference, which is a perfectly localized state with infinite lifetimes in the continuum [23]. In practice, an ideal BIC mode is transformed into a quasi-BIC mode that possesses a finite yet sharply high Q factor and narrow line width because of the existence of materials losses and other perturbations. The high Q factor of resonances governed by BIC mode could be achieved for both metallic metasurfaces [24] and all-dielectric metasurfaces [25,26]. It has been demonstrated that exciting Fano resonances can minimize the radiative losses due to their weak free space coupling enables long decay time [27], and multiple Fano resonances have become a hot research topic owing to their enormous potential for multichannel biosensors and multiwavelength surface enhancement [28,29]. However, the high Q factor of Fano resonance is obtained at the expense of its modulation depth, which makes the resonances difficult to be measured with low resolution and low signal-to-noise systems [15,30]. Therefore, it is of great significance to excite multiple Fano resonances with high Q factor as well as large modulation depth.

Here, we engineer an asymmetric metasurfaces structure composed of silicon double D-shaped resonator (DDR) arrays to realize multiple Fano resonances having not only ultrahigh Q factor but also near-unity modulation depth driven by nonradiative modes at THz range. An enhancing anapole mode with stronger electromagnetic field enhancement and higher Q factor can be observed by increasing the gap of DDR. In addition, two BIC modes are generated, and they can be transformed to quasi-BIC modes with an ultrahigh Q factor of $5\times 10^{7}$ by adjusting the asymmetry parameter of DDR. As far as we know, this is the first study to excite multiple Fano resonances having not only ultrahigh Q factor but also near-unity modulation depth governed by nonradiative modes. Our study will provide a new insight on realizing highly efficient THz devices.

2. Metasurfaces design

The proposed metasurfaces structure consists of a periodically arranged silicon double D-shaped resonator (DDR) arrays on a quartz substrate, see Fig. 1(a). The DDR (unite cell) is actually a split silicon disk realized by cutting a rectangle out of the center of the silicon disk along the y direction. Silicon and quartz was used as the materials of metasurfaces due to their low-loss characteristics and better capturing of the incident electric field in the THz band [31,32]. The optical constants of silicon and quartz are both referred to the Palik Handbook [33]. Figure 1(b) shows the top view (x-y plane) of the DDR, the radius of the silicon disk formed by the DDR is r, the distance between the center of the silicon disk and the left D-shaped resonator is ${d}_1$, and the distance between the center of the silicon disk and the right D-shaped resonator is ${d}_2$. In our study, the geometric parameters of metasurfaces are as follows: $P\rm {_x}={P}\rm {_y}={P}=230$ $\mathrm {\mu }$m, ${t}=103$ $\mathrm {\mu }$m, ${h}=500$ $\mathrm {\mu }$m, ${r}=100$ $\mathrm {\mu }$m and ${d}_1=16$ $\mathrm {\mu }$m, which ${d}_2$ is varied to control the asymmetry degree of the DDR. Here, the asymmetry parameter $\alpha$ is defined as $\alpha =\left |{d}_1-{d}_2 \right |/{d}_1$. Specifically, $\alpha =0$ (${d}_1={d}_2$) means the unit cell is symmetric in the x-y plane, while $\alpha \neq 0$ (${d}_1\neq {d}_2$) indicates the symmetry of the unit cell is broken in the x-y plane. It has been demonstrated that the higher order group $\rm {C}_{4v}$ of the unit cell could degenerate to lower $\rm {C}_{2v}$ and thus the BIC mode could be produced by breaking the in-plane symmetry in a metasurfaces structure [34]. Numerical simulations are carried out by using the finite-element method in our study. Unit cell boundary conditions with Floquet ports are applied in x and y directions. The perfectly matched layers and open boundary conditions are set in z direction of the unit cell. The THz plane waves are incident perpendicularly on the metasurfaces along the z direction with the electric field linearly polarized along the y direction.

 

Fig. 1. (a) Schematic diagram of the proposed all-dielectric metasurfaces consisting of periodically arranged DDR arrays. (b) The top view (x-y plane) of the DDR.

Download Full Size | PDF

The electromagnetic and scattering properties of the proposed metasurfaces are investigated based on the decomposition of the fields in the dielectric elements (DDR arrays) into Cartesian multipole moments. This multipole decomposition method in Cartesian coordinates allows for the calculations of the individual contributions stemming from magnetic dipole (md), magnetic quadrupole (mq), and electric quadrupole (eq) moments to the total scattering cross-section of the DDR arrays (the contribution of orders exceeding the quadrupoles is neglected). More importantly, the separated contributions of electric dipole (p) and toroidal dipole (t) moments can also be identified by employing this method. The crossing point of p and t respective scattering cross-section spectra can provide an indication of the anapole mode [13]. In this method, expressed in a Cartesian basis and assuming the ${e}^{j\omega t}$ convention for the harmonic electromagnetic fields, the induced polarization current density J(r) in the DDR is calculated as $J(r)=j\omega (\varepsilon _{p}-\varepsilon _{b})E(r)$, where $\omega$ is the angular frequency, $\varepsilon _{p}$ and $\varepsilon _{b}$ are the permittivities of the DDR and the background medium, respectively [22]. The detailed calculation of the scattering cross-section of DDR arrays can be referred to the previous works [35,36].

3. Results and discussions

3.1 Optical responses of symmetric silicon DDR arrays

The contributions of Cartesian md ($\rm {C_{sca}^{md}}$), mq ($\rm {C_{sca}^{mq}}$), eq ($\rm {C_{sca}^{eq}}$), p ($\rm {C_{sca}^{p}}$) and t ($\rm {C_{sca}^{t}}$) in the scattering cross-sections as well as the transmittance spectra for silicon disk arrays (${P}=230$ $\mathrm {\mu }$m, ${t}=103$ $\mathrm {\mu }$m, ${h}=500$ $\mathrm {\mu }$m, ${r}=100$ $\mathrm {\mu }$m) and the corresponding symmetric silicon DDR arrays (${d}_1={d}_2=16$ $\mathrm {\mu }$m) are calculated and shown in Fig. 2(a) and Fig. 2(d), respectively. Their corresponding normalized phases and phase differences between the $\rm {p_y}$ and $\rm {ikt_y}$ are presented in Fig. 2(b) and Fig. 2(e), respectively. As shown in Fig. 2(a), for the silicon disk arrays, a notable transmittance peak analogous to electromagnetically induced transparency (EIT) at 0.594 THz is observed. The frequency of this resonance peak coincides with the frequency of the crossing point of p and t contributions in scattering cross-section where the power radiation of the t moment and p moment are equal but out of phase, see Fig. 2(b). The cancellation of the far-field radiation between p and t occurs at this case, demonstrating that the condition for anapole mode is met [16]. It should be noted that the contributions from mq and md are relatively strong while the contribution from eq is very weak, showing non-zero radiation. The calculated electric field distribution (x-y plane) and magnetic field distribution (x-z plane) at the resonance peak are represented in Fig. 2(c). Two opposite circular surface currents (white arrows) are excited and three hot spots of electric field are observed within the silicon disk, giving rise to t moment oriented parallel to the silicon disk surface [16]. These circular currents induce corresponding magnetic field (orange arrows), which exhibits a particular vortex distribution behavior within the silicon disk. Such currents configuration is a distinctive feature of anapole mode. The calculated field enhancement factors of the electric and the magnetic field at the resonance peak are about 8 and 29, respectively. When a gap with $g=d_1+d_2$ is introduced to form the symmetric silicon DDR arrays, the anapole mode is enhanced, see Fig. 2(d), the EIT-like transmittance peak shows a blue shift and becomes narrower. As with silicon disk arrays, the condition for anapole mode is met, and the total scattering cross-sections exhibits non-zero radiation in this case. According to Fig. 2(e), a $\rm {\pi }$ phase difference of $\rm {p_y}$ and $\rm {ikt_y}$ is observed, indicating destructive interference effect between the p moment and t moment. Two opposite circular surface currents (white arrows) are generated within the DDR, which results in a strong t moment along with the symmetric silicon DDR surface, see Fig. 2(f). In this case, the magnetic field (orange arrows) has a particular vortex distribution within the symmetric silicon DDR in the x-z plane. It is worth to remark that the magnetic field is confined entirely inside the symmetric silicon DDR with a field enhancement factor value reaches about 56, almost double the value in the case of silicon disk arrays. And the electric field is effectively concentrated in the gap of symmetric DDR with a field enhancement factor value of 15 that is almost twice as much as the field enhancement factor in the case of silicon disk arrays, which strengthens the interactions between the THz waves and the analytes.

 

Fig. 2. Transmittance spectra and contributions of multipole modes in the scattering cross sections for (a) silicon disk arrays and (d) symmetric silicon DDR arrays ($g=32$ $\mathrm {\mu }$m). The normalized phase and phase difference between the $\rm {p_y}$ and $\rm {ikt_y}$ for (b) silicon disk arrays and (e) symmetric silicon DDR arrays. The electromagnetic field distributions at the corresponding resonance peak for (c) silicon disk arrays and (f) symmetric silicon DDR arrays.

Download Full Size | PDF

The transmittance spectra of symmetric silicon DDR arrays with gap g ranging from 0 $\mathrm {\mu }$m to 40 $\mathrm {\mu }$m are presented in Fig. 3(a). As g increases from 0 $\mathrm {\mu }$m (disk arrays) to 16 $\mathrm {\mu }$m, the resonance peak experiences a blue shift, and an asymmetric Fano resonance can be observed. When $g=24$ $\mathrm {\mu }$m, an EIT-like resonance can be observed because p and t modes are almost degenerate with each other. By further increasing the gap ($g>24$ $\mathrm {\mu }$m), an asymmetric Fano resonance appears again. These phenomena can be explained by the interference of t the p modes due to the change of gap [22,37]. The 2D map of transmittance spectra as a function of the gap g is shown in Fig. 3(b). With the increase of gap, the line width of resonance gets narrow. The Q factor and electric field enhancement factor of the symmetric DDR arrays with different g are calculated and shown in Fig. 3(c). The Q factor increases exponentially and the corresponding electric field enhancement shows a nearly linear increase with g increases. The Q factor can reach about 1630 and the maximum electric field enhancement factor is about 18.6 when $g=40$ $\mathrm {\mu }$m. These results indicate that the higher Q factor and the stronger electric field are available by using the symmetric silicon DDR arrays.

 

Fig. 3. (a) Transmittance spectra of symmetric silicon DDR arrays with different gaps. (b) The 2D map of transmittance spectra of symmetric silicon DDR arrays for different gaps. (c) The extracted Q factor and the electric field enhancement factor of symmetric silicon DDR arrays as a function of the gap.

Download Full Size | PDF

3.2 Symmetry-protected quasi-BIC of asymmetric silicon DDR arrays

It was demonstrated that symmetry-protected quasi-BIC could be generated by breaking the structural symmetry of the resonator in a metasurfaces system [18,38]. The symmetric breaking is obtained by changing $d_2$ in our present study to investigate the symmetry-protected BIC mode, where $d_1$ is fixed at 16 $\mathrm {\mu }$m. As shown in Fig. 4(a), three modes (marked as mode 1, mode 2 and mode 3) can be observed in the asymmetric DDR arrays. When $d_2$ is getting close to $d_1$, mode 1 and mode 3 show a slight blue shift with an increasingly narrow line width, and disappear at $d_2=d_1$ (marked as the pentagrams), indicating the existence of two symmetry-protected BIC modes. Then they appear again with wider line width as $d_2$ continues to increase. The appearance of mode 1 and mode 3 is attributed to the leakage of energy from the bound state to the free space, and such leakage causes additional radiative losses [39]. However, mode 2 shows a blue shift and its line width gets narrower as the increase of $d_2$. This is due to the fact that the increase of $d_2$ will increase the gap of the asymmetric DDR, which enhances the mode 2. Specifically, the transmittance spectra for the asymmetric DDR arrays with $d_2=11$ $\mathrm {\mu }$m, 13.5 $\mathrm {\mu }$m, 16 $\mathrm {\mu }$m, 18.5 $\mathrm {\mu }$m and 21 $\mathrm {\mu }$m are depicted in Fig. 4(b). Mode 2 can be found for all $d_2$ values and shifts to a higher frequency as $d_2$ increases. Mode 1 and mode 3 can be observed only under the condition of symmetry-broken, indicating that they are driven by the symmetry-protected BIC. When $d_2$ is equal to 16 $\mathrm {\mu }$m ($d_2=d_1$), ideal BIC modes with disappeared line width are excited, which means that there is no leaky energy to the free space. When $d_2\neq d_1$, mode 1 and mode 3 appear with finite line width, demonstrating these two symmetry-protected BIC modes are transformed to the corresponding quasi-BIC modes due to the symmetry-broken. In addition, mode 1 and mode 3 experience a slight blue shift as $d_2$ increases, which can be attributed to the decrease of the effective refractive index of the resonator as $d_2$ increases [40]. It is noted that mode 1 and mode 3 have not only ultrahigh Q factor but near-unity modulation depth. Take the metasurfaces with $d_2=18.5$ $\mathrm {\mu }$m as an example, it is found that mode 1 and mode 3 can be described rigorously by the classical Fano formula [22,41], see Fig. 4(c) and Fig. 4(d). The Q factor of the resonances can be calculated by $Q=\omega _{0}/\gamma$, where $\omega _{0}$ is the resonance frequency and $\gamma$ describes the resonance line width, and the calculated Q factors for mode 1 and mode 3 are $3.9\times 10^4$ and $1.4\times 10^4$, respectively.

 

Fig. 4. (a) The 2D map of transmittance spectra of the asymmetric silicon DDR arrays as a function of $d_2$. The color-coded arrows represent the transmittance for different $d_2$ of (b). (b) The transmittance spectra for the asymmetric silicon DDR arrays with different $d_2$. The results of Fano fitting for (c) mode 1 and (d) mode 3 when $d_2=18.5$ $\mathrm {\mu }$m.

Download Full Size | PDF

The behaviors of the Q factor for mode 1 and mode 3 at different asymmetry parameters $\alpha =\left |{d}_1-{d}_2 \right |/{d}_1$ are investigated. Figure 5(a) and Fig. 5(c) show the dependence of the Q factor on the asymmetry parameter $\alpha$ (log-log scale) for mode 1 and mode 3, respectively. It can be seen that the evolution of Q factor on $\alpha$ follows the universal behavior of BIC mode very well, exhibiting an inverse quadratic trend ($Q \propto {\alpha ^{{\rm {\ }\hbox{-}{\rm \ }}2}}$ ). As the $\alpha$ decreases, the radiation channels become narrower and the radiation loss reduces, resulting in an infinite Q factor resonance when $\alpha =0$ [42]. The Q factors for mode 1 and mode 3 can reach as high as $5\times 10^7$ and $1\times 10^7$ when $\alpha =0.006$, respectively. The scattering cross-sections of multipole moments of mode 1 and mode 3 when $d_2=18.5$ $\mathrm {\mu }$m are illustrated in Fig. 5(b) and Fig. 5(d). For mode 1, see Fig. 5(b), md has considerable scattering intensity and plays a dominant role while the scattering intensities of p and t are strongly suppressed. For mode 3, see Fig. 5(d), mq is the dominant multipole component, and p and t are suppressed.

 

Fig. 5. Dependence of the Q factor on the asymmetry parameter $\alpha$ (log-log scale), the solid black lines are linear fit for (a) mode 1 and (c) mode 3. The calculated scattering cross-sections of multipole moments for (b) mode 1 and (d) mode 3 when $d_2=18.5$ $\mathrm {\mu }$m.

Download Full Size | PDF

To further identify the physics of these two modes, their electromagnetic field distributions are investigated. For mode 1, see the upper panel of Fig. 6(a), the surface currents exhibit a circular behavior (white arrows) and the magnetic field $H _{\rm{z}}$ is mainly localized in the gap of the asymmetric DDR in the x-y plane with a field enhancement factor of 493. In this case, two magnetic field loops (orange arrows) with opposite directions are generated in the x-z plane, which produces md moment oriented along the z direction, see the lower panel of Fig. 6(a). For mode 3, see Fig. 6(b), the electric field and magnetic field are both tightly concentrated within the asymmetric DDR. The electric field distribution with two reversed loops in the y-z plane (while circular symbols) has a field enhancement factor value of 76. The magnetic field distribution in the x-z plane demonstrates two lineal behaviors going in opposite directions along the x axis, and a field enhancement factor value reaches about 316 can be achieved. Such electromagnetic field distributions correspond to mq moment. Therefore, those results indicate that mode 1 and mode 3 can be regarded as md and mq resonances, respectively.

 

Fig. 6. (a) The magnetic field distribution $H _{\rm{z}}$ in the x-y plane (upper panel) and magnetic field distribution in the x-z plane (lower panel) for mode 1. (b) The electric field distribution $H _{\rm{z}}$ in the x-y plane (upper panel) and magnetic field distribution in the x-z plane (lower panel) for mode 3.

Download Full Size | PDF

4. Conclusion

In conclusion, the proposed metasurfaces have been demonstrated the ability to excite the multiple Fano resonances driven by nonradiative modes with not only ultrahigh Q factor but also near-unity modulation depths in the THz region. By increasing the gap of the symmetric DDR, the excited resonance derived from anapole mode can be enhanced with a Q factor of about 1630 and with field enhancement factor as high as 18.6. Two Fano resonances governed by BIC modes can also be generated and they are transformed into quasi-BIC modes with an ultrahigh Q factor of $5\times 10^7$ by elaborately adjusting the asymmetry parameter. The multipole decomposition and electromagnetic field distributions demonstrate that these two BIC modes are dominated by the md and mq, respectively. It is worth noting that the Q factor of the multiple Fano resonance modes is well-controllable by adjusting different geometric parameters, making them ideal for practical applications. Our findings pave the way for the development of high-performance THz devices such as multi-channel sensors and optical modulators.