We start with the investigation of transmission spectra from Q-BIC metasurface with in plane *C**2* rotational symmetry, which two split gaps of the graphene DSRR are located at the center of the left and right arms and corresponding to asymmetric parameter *δ* = 0. One can see from the black line in Fig. 2(a) that the graphene DSRRs show an obvious plasmonic resonance at 0.59 THz when excited by the incident THz wave along *x*-direction. The transmittance at this frequency point is 0.1. In this simulation, the chemical potential of the graphene is *E**f* = 0.4 eV. The form of the electric field distribution at the frequency of the minimum value of transmittance demonstrates the trapped mode nature of the resonance, as shown in Fig. 3(c). The electric field of plasmonic resonance mode is mostly concentrated in the left and right branches of the graphene DSRR. There are the in-phased surface current flows parallel to the *x-* and *y*-axis of the DSRR elements (white arrows in Fig. 3(c)). The amplitude of the surface current flow parallel to the *x*-axis is greater than that parallel to the *y*-axis. This amplitudes distribution of currents provides a strong coupling between the trapped mode and the incident THz wave [23].

We then plot the transmission spectrum of the periodic graphene DSRRs with distinct asymmetric parameters *δ* in Fig. 2(a). Besides the plasmonic resonance at 0.59 THz, there are two new transmission dips appear with the introduction of asymmetric parameter *δ* = 0.5 µm. The exhibited two new narrow resonances are strongly altered with increasing *δ* from 0 to 4 µm. The two new transmission dips become prominent with the increase of *δ*, and the linewidth gradually increase. Meantime, the two new transmission dips are slightly red shifted. As the non-zero value of asymmetric parameters *δ* breaks the *C**2* rotational symmetry of the graphene DSRR, the leaky energy channel between the resonances and the free space continuum can be built. At *δ* = 0, the two new transmission dips vanish in the transmission spectrum, which implies that the energy leakage channel is closed by the in plane *C**2* rotational symmetry. This is the ideal BIC point as marked by red dashed circles in Fig. 2(a). The two new resonances have finite linewidth induced by the leakage of ideal BICs are two quasi-BICs, and they are named as Q-BIC-I in the low frequency and Q-BIC-II in the high frequency, respectively. To further confirm, the relationship between Q-factors and asymmetry degree *α* are discussed in Fig. 2(b). The dimensionless asymmetry degree is defined to be, *α* = 2*δ* / (*l**1* + *l**2*) × 100%. Here, *l**1* and *l**2* are the total length of upper and lower branches of the graphene DSRR, respectively. When *δ* = 0.5 µm, and *l**2* = *l**1* − 2*δ =* 14 µm, *α* equals to 3.45%. The value of Q-factors of Q-BIC-I and Q-BIC-II are 349.75 and 1456.46, respectively. With the decrease in the asymmetry degree *α*, the increasing Q-factors of modes Q-BIC-I and Q-BIC-II are observed. The fitting solid lines in Fig. 2(b) indicate the Q-factors follows the asymmetry degree as Q∝*α*−2, and diverge at the ideal BIC points.

All the Q-factor values of modes Q-BIC-I and Q-BIC-II discussed above are analyzed by using the following normalized Fano’s equation [24]:

$${T_{Fano}}=a+\frac{{{{(b+c)}^2}}}{{({b^2}+1)(1+{c^2})}}$$

2

,

where *a* donates the real Fano fitting number. *b* = 2(*f* - *f**0*)/*Γ**tot*, where *f* donates the broadband frequency of incident THz wave, *f**0* donates the central plasmonic resonant frequency, and *Γ**tot* is the total loss. *c* is a non-zero Fano parameter for Q-BICs. Since the transmission dips are due to plasmonic Q-BICs, the Q-factor is given by [10]:

$$Q=\frac{{{f_0}}}{{{\Gamma _{tot}}}}$$

3

.

The match of fitting results via Eq. 2 with the numerical spectra is shown in Fig. 3(b) by choosing the asymmetry parameter *δ* = 4 µm. The fitting spectral line (the red solid line) is consistent widely with the simulation curve (the black solid curve). The full transmission spectra of Q-BICs (*δ* = 4 µm) in comparison with ideal BIC (*δ* = 0) is shown in Fig. 4(a). The spectral contrast ratio is defined as:

$$\frac{{{T_{\hbox{max} }} - {T_{\hbox{min} }}}}{{{T_{\hbox{max} }}+{T_{\hbox{min} }}}} \times 100\%$$

4

,

where *T**max* and *T**min* donate the maximum and minimum values of transmittance. Mode Q-BIC-I at 0.44 THz has a spectral contrast ratio larger than 70%, and that of Q-BIC-II at 1.19 THz is close to 50%.

In order to clarify the transition from ideal BICs to Q-BICs more clearly, the electric field distributions and surface current flow distributions at *xoy* sections at each Q-BIC resonance frequency are presented at Fig. 3(d) and (e). Under linearly polarized THz plane wave incidence, for DSRRs with *C**2* rotational symmetry, only collective modes that have net effective electric dipole moments can be excited due to the symmetry incompatibility [25]. Anti-phased and equal magnitude currents parallel to the *x*- and *y*-axis will result in the zero effective electric dipole moments, which cannot cause the coupling of ideal BIC to the incident THz wave. When the gap of the graphene DSRR is displaced (*δ* ≠ 0), the symmetry of surface current flow is broken in terms of magnitude, as shown in Fig. 3(d). The currents at 0.44THz parallel to the *x*- and *y*-axis are anti-phased while the magnitude of current in lower branch is obviously higher than that of upper branch. Therefore, the symmetry of electric dipole moments will be broken not only in terms of magnitude, but also in terms of the orientation, and leads to non-zero net effective electric dipole moments which can cause the coupling of ideal BIC to the incident THz wave. As a result, the ideal BIC transforms into Q-BIC which reveals the characteristic of symmetry protection. Similar trend is observed for mode Q-BIC-II at 1.19 THz. Compared with mode Q-BIC-I which has two branches of current flows, there are four asymmetric branches of current flows are observed in Fig. 3(e). Here, more nodes of asymmetric electric field distribution and more branches of current distributions also explain why the mode Q-BIC-II can be excited at the higher frequency.

The most attractive property of graphene plasmonic resonance is its tunable chemical potential without structure reconfiguration via chemical doping, electrical gating, and ion gel as top gates [26]. To further expand the Q-BIC’s operating frequency range, the transmission spectra of our Q-BIC metasurface with different chemical potentials are calculated in Fig. 4. For the graphene DSRRs with asymmetric parameter *δ* = 4 µm, we change the chemical potential from 0.4 eV to 0.55 eV. With the increase of chemical potential, there is a blueshift of both two Q-BICs resonances. The spectral contrast ratio of mode Q-BIC-I and Q-BIC-II are also increased. At higher chemical potential, graphene becomes more absorbing because of the band filling. Therefore, the total incident wave energy loss due to transmission is less.

To proceed, we perform an angle-dependent calculation for our proposed graphene metasurface. The asymmetric parameter is *δ* = 4 µm and the graphene’s chemical potential is *E**f* = 0.4 eV. The 2D mesh plot in Fig. 5 shows that the mode Q-BIC-I is constantly maintained at 0.44 THz for the incident angle from 0° to 70°, with a spectral contrast ratio larger than 70%. The mode Q-BIC-II is constantly maintained at 0.19 THz as well when changing the incident angle. At higher angles > 60°, its spectral contrast ratio begins to decrease, and the Q-BIC point is slight shift in the spectrum. The operating frequencies of Q-BICs are loosely dependent on incident angle because the geometric dimensions of the graphene DSRR are much smaller than the incident wavelength, thus the phase difference at neighboring graphene DSRR is rather independent of incident angle. Figure 5 shows our dual band Q-BIC metasurface are immune to the variation of incident angle.

We also perform the calculation of slow light effect for our proposed Q-BIC metasurface. In this simulation, the asymmetric parameter is *δ* = 4 µm and the graphene’s chemical potential is *E**f* = 0.4 eV. The slow light effect can be evaluated by group time delay *τ**g*, which is given by:

$${\tau _g}= - \frac{{d\phi }}{{d\omega }}$$

5

,

where *ф* donates the phase change of the transmission spectrum. In Fig. 6(a), we show the phase oscillates rapidly from − 1.39 rad to -0.42 rad and from − 2.41 rad to -1.98 rad as close to the Q-BIC-I and Q-BIC-II region, respectively. The phase *ф* changes monotonously like usual plane wave while away from the Q-BIC region. Figure 6(b) shows the largest positive time delay of 23.8 ps and 14.5 ps, which represent the remarkable slow light effect, can be achieved by the mode Q-BIC-I and Q-BIC-II, respectively. The very slow group velocity when light passing through the metasurface will extremely enhance the light-matter interaction for various applications.

Thereafter, we investigate the sensing characteristics of our proposed Q-BIC metasurface. The asymmetric parameter and the graphene’s chemical potential are fixed to *δ* = 4 µm and *E**f* = 0.4 eV, respectively. We calculate the change of transmission spectrum by varying the refractive index (*n*) of external liquid from 1.33 to 1.37. As shown in Fig. 7(a) and (c), both two Q-BICs’ resonant frequency has a redshift when increasing the refractive index. As we know that, the important performance index of the refractive index sensor is its sensitivity, which is defined as [27]:

.

Figure 7(b) and (d) show the frequency shifts of the mode Q-BIC-I and Q-BIC-II corresponding to different refractive indexes, respectively. For mode Q-BIC-I in Fig. 7(b), the curve shows a relatively linear dependence, and the slope is 96 GHz/RIU. Similar trend is observed for mode Q-BIC-II in Fig. 7(d), while the slope is relative larger and can reach up to 267.5 GHz/RIU. The measurable frequency offset caused by a slight change in the refractive index of external liquid indicates that our proposed Q-BIC metasurface can operate as a THz liquid sensor. The sensing characteristics can also be quantified by calculating the figure of merit (FOM). The FOM is defined as [27]:

$$FOM=\frac{S}{{{\Gamma _{tot}}}}$$

7

.

The maximum FOM is calculated to be 5.87 RIU− 1 and 24.08 RIU− 1 for the mode Q-BIC-I and Q-BIC-II, respectively.

In addition, the sensing performance of the proposed graphene metasurface are compared with the other types of recently reported THz sensors as shown in Table 1. The sensitivities and FOMs of several kinds of sensors are listed. Obviously, our work having high sensitivity is superior to the previous works [27–32] in sensing applications. It is worth noting that the sensitivities of Ref. [30] and [32] are higher than that of our proposed Q-BIC metasurface, however, their FOMs are much lower. When the analyte concentration of external liquid is very low, a sensing platform with high FOM is typically desired. Table 1 shows that our proposed Q-BIC metasurface exhibits improved sensing performance over most of the recently reported THz sensors. The simplicity of our tunable design will provide a valuable reference for future dynamic sensor.

Table 1

Comparison of sensing performance of the proposed Q-BIC metasurface.

Reference | Sensitivity (GHz/RIU) | FOM (RIU− 1) |

[27] | 170.58 | - |

[28] | 187.1 | 285.8 |

[29] | 105 | 7.501 |

[30] | 533.7 | 6.03 |

[31] | 266 | - |

[32] | 597 | 9.48 |

This work | 96 267.5 | 5.87 24.08 |