Bound states in the continuum (BICs) are bound states that lie inside the continuum. Under such states the energy can be perfectly confined within the system and without any radiation. Therefore, a BIC can be considered as a resonance with zero linewidth or an infinite Q factor [1,2]. In theory, BICs are dark modes with infinite radiative lifetime, which are “invisible” for far-field measurement, but they can be observed at near-field experiments [3]. However, due to the finite size of structures, material absorption, and other external perturbations, BICs tend to collapse to a Fano resonance with a limited radiative Q factor; such resonances are known as quasi-BICs [4–6]. There are two main types of BICs in metamaterials: symmetry protected BIC [1] and Friedrich–Wintgen (F–W) BIC [7]. The coupling coefficient could vanish due to the symmetry reason when the spatial symmetry of the mode is incompatible with the symmetry of the outgoing radiating waves. This kind of BIC is called symmetry protected BIC. The ideal symmetry protected BIC can be transformed into quasi-BIC by constructing a radiation channel through oblique incidence or breaking the structural symmetry [8–11], while the F–W BIC is generated by the destructive interference of the two eigenmodes, which requires precise adjustment of parameters of the structure. If the structural parameters deviate from the ideal BIC point, it becomes a quasi-BIC with a finite Q value [12–18]. In addition, BICs can be found and originated from guided mode resonances (GMRs) or surface lattice resonance in dimer or multimer configurations [19–22]. They become GMRs or bright modes (electric and magnetic dipole modes) with a finite Q factor when the lattice symmetry is broken by changing the position or size of its building blocks. A BIC provides a new way for the realization of high-Q metasurfaces and has been applied in high-sensitivity sensing [23–26], low-threshold micro-nano lasers [27,28], and nonlinear harmonic generations [29,30].

In metasurfaces, high-Q Fano resonance was mainly obtained by breaking the symmetry of the structure [31–34], and its Q value was determined by the asymmetry parameter of the structure. Theoretically, the smaller the asymmetry degree, the higher Q value of the metasurface. Hence, metasurfaces with arbitrary high Q value can be easily designed. However, it is very difficult to be achieved in experiments. So far, the reported highest Q value in the terahertz (THz) band is 227 [32], and 1011 in the optical regime [35–37]. In addition, high-Q resonances can also be achieved by exploiting surface lattice resonances in the symmetrical plasmonic metasurface, and a Q value of 2340 was obtained experimentally at the optical wavelengths [38–40]. Recent research on BICs has greatly promoted the development of high-Q dielectric metasurface with both symmetric and asymmetric structures, especially at optical frequencies [12,41–44]. For example, in a silicon cubic metasurface supporting symmetry protected BIC, a Q value of a quasi-BIC obtained reached up to 18,511 [43], and meanwhile, optical second and third nonlinear harmonics were observed simultaneously. However, it remains challenging for long years to achieve high-Q resonances in terahertz metasurfaces because of the absorption loss of dielectric materials, the design and fabrication of nanostructures, and high-signal-to-noise ratio high-resolution spectral measurements. To date, the highest Q value reported experimentally is only 250 in a silicon BIC terahertz metasurface [14].

In this work, we investigate and experimentally demonstrate an ultra-high-Q quasi-BIC resonance in an all-silicon metasurface with perforated air hole tetramers. The metasurface supports two BICs originated from GMRs in a frequency range of 0.4–0.9 THz, which collapses into two quasi-BIC Fano resonances with finite Q values through changing the radius of air holes in diagonal configuration. Thanks to the designed metasurface with strong coupling of the unit cells and without substrate avoiding material absorption loss, we successfully achieved a polarization-insensitive quasi-BIC with a record Q factor as high as 1049 in experiments. Further, we demonstrate a terahertz intensity modulation with very low optical pump power by exploring photon-induced carrier absorption of semiconductors. The perforated air-hole all-silicon terahertz metasurfaces have been rarely reported experimentally [45,46]. The work in this paper could pave the way not only for the research and application of high-Q terahertz metasurfaces, but also for other functional devices [47,48].

An all-dielectric metasurface we designed consists of perforated air-hole tetramers arrays in a high-resistance silicon wafer shown in Fig. 1. The all-silicon metasurface without substrate is compatible with the CMOS technique and ease of fabrication; meanwhile, its Q factor can be improved by avoiding the material absorption loss of the substrate.

In Fig. 1, the radii of the air holes at the two diagonals are represented by r1 and r2, respectively, and the thickness of the silicon wafer is 150 μm. The periods of the unit cell in the x and y directions are Λx=Λy=Λ=300μm, and the distance between the centers of the two circles along the x and y directions is fixed at Λ/2=150μm. Numerical simulations were carried out using the finite element method (COMSOL Multiphysics software), in which periodic boundary conditions are employed in the x and y directions, a perfectly matched layer is employed in the z direction, and the dielectric constant of silicon is set to be 11.67.

First, at fixed r2=55μm, the transmission spectra of metasurfaces with different r1 were calculated and shown in Fig. 2(a). There are three strong Fano resonances in the range of 0.4–0.9 THz, marked as f1, f2, and f3, respectively. The resonance frequencies and Q factors of f1 and f2 are shown in Figs. 2(b) and 2(c), and Q values are extracted from transmissions by Fano fitting as [49,50] T(ω)=T0+A0[q+2(ω−ω0)/γ]21+[2(ω−ω0)/γ]2,(1)where q is the Fano fitting parameter that determines the asymmetry of the resonance curve, ω0 and γ represent the resonance peak angular frequency and the resonance linewidth, respectively, T0 is the transmittance baseline shift, and A0 is the coupling coefficient, therefore, Q=ω0/γ.

As r1 increases from 25 to 75 μm, the three resonances have obvious blueshifts, due to the decrease of the effective refractive index of the metasurface. In particular, for a symmetric metasurface when r1=r2=55μm, both resonances f1 and f2 located at ∼0.577THz and 0.739 THz, respectively, disappear, i.e., the Q factor is infinite. Moreover, on both sides of r1=55μm, when r1 is away from 55 μm, the Q factors of f1 and f2 decrease dramatically, showing typical BIC characteristics. While the resonance f3 is very different from f1 and f2, its Q factor decreases rapidly as r1 increases from 25 to 75 μm and is finite when r1=55μm. In addition, when r1 is in the range of 25–40 μm, the Q value of f1 is much larger than that of f2; while in the range of 60–75 μm, the Q value of f1 is just slightly higher than that of f2. Moreover, the transmittances of the two quasi-BIC resonances are polarization-insensitive according to the geometrical symmetry of the structure.

To interpret the physical origin of the three resonances, the related intrinsic band structures of the TE and TM modes from the metasurfaces at lattice constants of 300 μm and 150 μm were investigated by COMSOL Multiphysics. Here, r1=r2=55μm; thus, the metasurface with a lattice constant of 300 μm degenerates to a metasurface with a lattice constant of 150 μm. As we can see from the TE modes in Fig. 3(a), the band structure of the 300 μm-period unit cell (blue line) is folded from that of the 150 μm-period unit cell (red line). Meanwhile, bound states in the discrete region of the 150 μm-period unit cell would be folded inside the light cone, which results in bound states in the continuum [20–22]. From the enlarged view of our considered modes, there are two eigenmodes TE 1 and TE 2 (related to f2), both with infinite Q factors at Γ point and off Γ point (not shown in the figure), which does not exist at the lattice constant of 150 μm. Considering their magnetic field distributions at Γ point in Fig. 3(c), it can be inferred that the two TE modes are BICs originated from GMRs induced by (1, 1) Rayleigh diffraction from a two-dimensional (2D) metagrating at a 300 μm lattice constant, and propagating along 45° and −45° to the x axis [20,51,52]. And resonance f2 is the result of the co-excitation of these two GMRs.

Similar to the above two TE modes, two TM modes (related to f1), TM 1 and TM 2, both with infinite Q factor at Γ point and off Γ point, only exist in the metasurface at the lattice constant of 300 μm [shown in Fig. 3(b)]; they are two BICs originated from GMRs propagating along 45° and −45° to the x axis, as shown in Fig. 3(c). In addition, TM 3 and TM 4 are two GMRs induced by (1, 1) Rayleigh diffraction from a 2D metagrating at 300 μm lattice constant and propagating along 45° and −45° to the x axis as well, while TE 3 and TE 4 are two GMRs induced by (1, 0) and (0, 1) diffraction and propagating along the x axis and y axis, respectively. However, TM 5 and TM 6 (related to f3) with a finite Q factor are calculated at the lattice constant of 150 μm [shown in Fig. 3(b)]. According to their electric field distributions shown in Fig. 3(c), they are GMRs induced by first-order diffraction (1, 0) and (0, 1), and propagating along the x and y directions, respectively.

Thus, quasi-BICs or GMRs with a finite Q factor can be experimentally observed when the lattice symmetry is broken [20–22,52]. However, the lattice symmetry of a 2D metasurface has more complexity and diversity than that of a 1D grating [19,20]. For example, the sizes of the air holes can be tuned (by diagonal configuration in Fig. 1, or ipsilateral configuration in Fig. 4, or only changing one air hole) to modulate the TE and/or TM BICs. In addition, the position of the air holes can be shifted either in the x axis or in the y axis to modulate the GMRs as well. But they are maintained at the oblique incidence [20].

Figure 4 shows the GMRs modulated by changing the radius of the air holes on the ipsilateral side. Here, resonance f1 corresponds to a GMR, marked as TE 3 in Fig. 3, while f2 is related to the two GMRs, TM 3 and TM 4. Likewise, as r1 increases, both resonances f1 and f2 are blueshifted. When r1=r2=55μm, the two resonances disappear. Moreover, when r1 is away from 55 μm, the Q factors of the two resonances f1 and f2 decrease rapidly. It can also be seen from the Q characteristics that when r1 increases from 25 to 50 μm, the Q value of the resonance f1 increases from 63 to 1950 and 76 to 1420 for f2. Meanwhile, the Q value of f1 is significantly larger than that of f2 when r1 is in the range of 55–75 μm. However, the Q factors of the two quasi-BICs obtained by changing the ipsilateral air holes are much smaller, as compared with the diagonal configuration; this means that the asymmetry degree of a metasurface in the ipsilateral air holes case is larger than that in a diagonal one.

We fabricated high-Q metasurfaces in a 2 in. high-resistivity silicon wafer (>5000Ωcm) with a thickness of 200 μm. We chose high-resistance silicon as the building material of the metasurface because of its relatively low loss and low dispersion in the terahertz band. Here, we followed a three-step fabrication procedure. First, the silicon wafer was thinned to 150 μm, then patterned by conventional UV photolithography, and finally, etched using a deep reactive ion etching process. Since quasi-BIC resonance with higher Q factor can be calculated from the diagonal configuration, three polarization-insensitive metasurfaces with a size of 15mm×15mm were prepared according to the structure shown in Fig. 1, where r2=55μm, and r1 is chosen as 25, 30, and 35 μm, respectively. It should be noted here that the results in Fig. 2 are calculated by assumption of the cylindrical air holes in the structure; the real sample fabricated by deep etching is of circular truncated cone with an inclination angle about 2°, which can be seen from the sample’s SEM picture in Fig. 1(b). Hence the volume of the air holes is larger than that of the cylindrical air holes, which will influence the quasi-BIC characteristics. Therefore, we recalculated the transmission spectra of the metasurface with circular truncated cone air holes, as shown in Fig. 5(a), and material loss of the high-resistance silicon (loss tangent of 0.0001) is also considered in the calculation.

Compared with the calculations of the cylindrical air hole metasurface in Fig. 2(a), when r1 increases from 25 to 35 μm, the frequency of resonance f1 changes little, while greatly for f2. In addition, it is interestingly found that the Q factor of the resonance f2 of the real metasurfaces is greatly increased compared with that of the corresponding cylindrical air hole metasurfaces. This also shows that the coupling of the unit cells in the real metasurface is greatly enhanced, which is quite beneficial for the realization of high-Q resonance. For example, when r1=25μm, the Q value of f2 increases from 328 to 845, and when r1=30μm, the Q value of f2 increases from 652 to 2088.

We used a high-resolution terahertz frequency-domain spectroscopy system (TeraScan 1550 from Toptica) with a focused beam incidence on the sample to measure the transmissions of three samples, and the diameter of the spot illuminating on the sample is about 6–8 mm. The measurements were done at room temperature and in the dry air condition (humidity ≤1%) to nullify the effect due to water vapor absorption, and the integration constant was set to be 300 ms. The measured transmission spectra are shown in Fig. 5(b), and enlarged views of the resonance f2 and Fano fitting are shown in Fig. 5(c). The results show that a very weak quasi-BIC resonance f1 at 0.565 THz was observed for the sample of r1=25μm, but not clearly seen in the other two samples. This may be attributed to the very high Q values of f1 (≥6000 calculated). The higher the Q factor is, the more it is affected by the material’s absorption, fabrication defects, or device size. And it may also be related to the limited spectral resolution of the measurement instrument (140 MHz). On the contrary, the strong resonance f2 in all three samples was observed. The frequency of f2 measured is about 10 GHz larger than the calculation result, which is caused by the inconsistency of the size between the simulation and the sample. Furthermore, as r1 increases, the intensity of resonance f2 measured becomes weaker, which is in good agreement with the simulation result. In addition, the measured Q value of f2 for the sample of r1=25μm is 787, which is very close to our numerical calculation result of 845. And for the sample of r1=30μm, the calculated Q value of f2 is 2088, while the measured Q value is 1049, which is the highest reported value for THz metasurfaces. However, for the sample of r1=35μm, the measured Q value decreased slightly to 919. This may be related to the limited spectral resolution of the instrument and the fabrication imperfections. Moreover, for the third resonance f3, it is weak for the samples of r1=25μm and 30 μm, but strong resonance is observed when r1=35μm. This agrees with simulations, in which the Q factor of f3 decreases rapidly as r1 increases from 25 to 35 μm.

Moreover, we have also measured the samples using a backward wave oscillator THz spectrometer with spectral resolution better than 10 MHz, but all three resonances f1, f2, and f3 are not observed due to low SNR of the system in the measured frequency band.

Finally, it is well known that high-Q silicon photonic crystals and metasurfaces can be used as efficient THz switches or modulators by exploring the photon-induced carrier absorption of semiconductors [14,53,54]. Here, we demonstrate an optically controllable THz modulation with remarkably reduced optical pump power by using our high-Q metasurface. In the experiment, a continuous-wave laser diode with a wavelength of 808 nm (photo energy 1.53 eV, above the bandgap of silicon 1.1 eV) was used to photoexcite the all-silicon metasurface, and the spot size on the metasurface is about 1mm×2mm. Figure 6 shows transmissions of the metasurface (r1=25μm) when irradiated obliquely by a laser beam with different power. It is clear to see that the high-Q resonance f2 is much more sensitive to the photoexcitation than other frequencies; hence we focus on the resonance f2 shown in Fig. 6(b). The quasi-BIC mode is continuously modulated in terms of resonance intensity as well as Q factor by gradually tuning the photoexcitation pump power. As the pump power increases from 0 to 40 mW, the high-Q resonance of f2 at 0.65 THz gradually becomes weak and has a slight blueshift of 1 GHz as well. And it eventually quenches at a very low pump power of 100 mW. However, a pump power of 600 mW is far from quenching a low-Q-factor mode (0.55–0.8 THz) [see Fig. 6(a)].

To summarize, we have proposed and experimentally demonstrated an ultra-high-Q all-silicon BIC metasurface composed of air-hole tetramer arrays. In the range of 0.4–0.9 THz, the metasurface supports two BICs originated from the GMRs, which collapse into two quasi-BIC resonances with a finite Q factor by breaking the lattice symmetry. The TE and/or TM GMRs can be flexibly modulated by changing the sizes of the air holes in both diagonal configuration and ipsilateral configuration, or by shifting the position of the air holes either in the x-axis or in the y-axis. Furthermore, we demonstrate experimentally a polarization-insensitive GMR with a record Q factor of 1049, which is 3 times higher than the highest reported in terahertz metasurfaces. Furthermore, we demonstrate a terahertz intensity modulation with very low optical pump power by exploring photon-induced carrier absorption of semiconductors. The work in this paper could pave the way not only for the research and application of high-Q terahertz metasurfaces, but also for other functional devices.