Excitation of Tunable Dual Quasi-Bound States in the Continuum in Graphene Metasurface and Terahertz Sensing Application

This study proposes an all-graphene metasurface supporting dual band symmetric bound state in the continuum (BIC) in the terahertz (THz) range for the first time. The structure consists of a unit cell containing a double-gap split ring resonator periodically on a dielectric substrate. By introducing symmetry breaking, two plasmonic quasi-BIC (Q-BIC) transmission dips are observed with finite Q-factors and high modulation depth. Simulation and analysis results simultaneously exhibit that the Q-factors of Q-BICs follow an inverse square dependency with the asymmetry degree. Via changing the graphene’s chemical potential, the Q-BICs’ operating frequency range can be actively expanded. The dual Q-BICs are immune to the variation of incident angle and have a significant slow light effect with time delay up to 23.8 ps. In addition, the sensing performance in THz region is investigated. A maximum sensitivity of 267.5 GHz/RIU is obtained with a FOM of 24.08 RIU−1. Our work shows an alternative way to design the class of tunable Q-BIC metasurface, which will provide a valuable reference for future dynamic sensor and other fields.


Introduction
In 1929, the idea of bound states in the continuum (BIC) was firstly reported by Neumann and Wigner in quantum mechanics during the process of unique Schrödinger equation calculation [1].The BIC, however, is often brittle, necessitates uniquely tuned potentials, and is hence challenging to observe in the quantum experiments.Because the Maxwell equation and the Schrödinger equation are formal equals, the emergence of these exotic non-radiating states in optics is strongly influenced by the quantum mechanical system [2].Similar to quantum systems, the concept of BIC in optics is the result of two or more waves interfering destructively to reduce radiative loss, where radiation resonances can generally carry power from or to infinity.At the operating frequency of BIC, the embedded eigenvalues vanish with the leakage power and cause eigenstates having an infinitely high-quality factor (Q-factor) [3].Generally, there are two ways to excite BIC resonance.One is tuning the parameters of optic systems to make outgoing waves interfering destructively, and this type of BIC is named as accidental BIC [4].The other way is breaking the in-plane symmetry to generate mismatch between the symmetry of outgoing wave and the resonance mode, and this type is symmetry-protected BIC [5].Due to its novel ability to have the high Q-factor, BIC has been widely studied in various optical systems of quantum dot lasers [6], filters [7], sensors [8], and nonlinear optics [9].
Recently, BIC has been introduced into the metasurfaces [10] owing to the advantage of electromagnetic wave confinement in sub-wavelength, precise wavefront control, and the simplicity of fabrication in two-dimensional planes [11].The periodic arrangements of artificial sub-wavelength unit cells in metasurfaces result in some diffraction open channels.The symmetry-protected BIC is easier to be existed in the absence of coupling to these open channels by engineering the periodic sub-wavelength structures [12].Interestingly, true optical BIC (ideal BIC) is a mathematically defined concept with the infinite value of Q-factor, and its resonance linewidth infinitely tends towards 0. This is only possible in an ideal lossless infinite optical system with extreme structure parameters.In practice, suffering from material absorption, finite lateral size of structures, and other leakage, the observable BIC often refers to quasi-BIC (Q-BIC) with high, yet finite Q-factor [13].Although Q-BIC has received extensive attention in metasurfaces, most of the devices reported in the literatures are lack of dynamic capability.Once the manufacturing parameters are fixed, the operating frequency of Q-BIC is synchronous fixed.So the ability to expand the Q-BIC's operating frequency range without changing the geometrical parameters of the structure is very important and attractive.
Currently, graphene has exhibited powerful tuning ability for active devices [14,15].By controlling the chemical potential E f of the graphene, the target operating frequency region can be dynamically shifted [16].However, the active tuning of Q-BIC reported in the literatures is mainly focused on the switching characteristics controlled by graphene absorption.For example, the authors in ref. [17] proposed dual Q-BICs in a silicon stair-like metasurface, where the transmittance can be controlled with respect to a change in the chemical potential of graphene layer.Similarly, the authors in ref. [18] have reported graphene-metal hybrid terahertz metasurfaces with a similar kind of control.So far, to the best knowledge of us, only two groups have reported dynamic tuning the operating frequency of Q-BIC via the changing of the chemical potential E f of the graphene.Wang et al. have investigated the symmetry-protected Q-BIC by using the structure of two contacting graphene nanoribbons [19], and Roy et al. have proposed rectangular dimer-shaped graphene metasurfaces that supported asymmetric Q-BIC [20].Therefore, the tunable aspects of graphene in Q-BIC are still lack of research.Besides that, tunable Q-BIC with only one operating band also faces time delay introduced by modulation, especially in sensing applications in detecting a variety of gas and chemical solvents.
In this article, the exciting of symmetry-protected dual Q-BICs in an all-graphene-based metasurface is theoretically demonstrated for the first time.The metasurface is made of a square lattice of which unit cell consists of a double-gap split ring resonator (DSRR).By displacing the two gaps from the center, the in-plane C 2 rotational symmetry is broken, as shown in Fig. 1.This symmetry breaking provides an efficient way to transform BICs into Q-BICs, manifested as two new sharp dips in the transmission spectrum.
The Q-factor of dual Q-BICs can be well evaluated by the theoretical normalized Fano line shape equation and follow the inverse square dependence on the asymmetry degree.The dual Q-BICs are further analyzed by the near-field electric field distributions and surface current flows, indicating that the leakage of an ideal bound state can be caused by a slight perturbation on the effective electric dipole moments.We show that the operating frequency of the dual Q-BICs can be actively tuned by varying the chemical potential E f of the graphene.Our Q-BIC metasurface is immune to the variation of incident angle and has a significant slow light effect with time delay up to 23.8 ps.Moreover, our Q-BIC metasurface can operate as a refractive index sensor in the THz region.A maximum sensitivity of 267.5 GHz/RIU is obtained with a FOM of 24.08 RIU −1 .Our work shows an alternative way to design the class of tunable Q-BIC metasurface, which will provide a valuable reference for future dynamic sensor and other fields.

Design and Simulation of Dual Q-BIC Metasurface
The schematic of the proposed Q-BIC metasurface is shown in Fig. 1a.On a thick SiO 2 substrate, a monolayer graphene is patterned into periodic DSRRs along x and y directions.The relative permittivity of the THz transparent SiO 2 substrate is taken as 3.9.Figure 1b shows the unit cell structure.Here, the periods along x and y directions are P x and P y , and P x = P y = 50 μm.The arm length of the DSSR is l = 31 μm, and the arm width is w = 4 μm.The gap length g = 2 μm.In order to break the in-plane C 2 rotational symmetry of the DSSR, we displace both two gaps from the center by a distance "δ."In this work, we use open source software "East FDTD" [21,22] based on the finite-difference time-domain method to design our structure and calculate the transmission spectrum numerically.Along the x-and y-axis, we choose periodic boundary conditions to simulate a system of infinite size.Along the z-axis direction, the metasurface is excited by a normally incident broadband terahertz plane wave polarizing along x-direction, and two perfectly matched layers (PMLs) behind Floquet periodic ports are used to absorb the scattered light.For the numerical calculations, graphene is treated as a surface conductive layer with the thickness of 1 nm, and its surface conductivity is defined by analytical Kubo's formula [21,22]: where ω represents the frequency of incident terahertz wave, T represents the environmental temperature, e represents an electron charge, k B represents the Boltzmann's constant, ħ represents the reduced Planck constant, and E f represents the chemical potential of graphene DSSRs.In our numerical calculations, room temperature is taken as T = 300 K, and the relaxation time of graphene DSSRs is taken as τ = 2 ps.

Results and Discussions
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. 2a 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 (1) of the resonance, as shown in Fig. 3c.The electric field of plasmonic resonance mode is mostly concentrated in the left and right branches of the graphene DSRR.There are the inphased surface current flows parallel to the x-and y-axis of the DSRR elements (white arrows in Fig. 3c).The amplitude of the surface current flow parallel to the x-axis is greater than that parallel to the y-axis.This amplitude 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. 2a.Besides the plasmonic resonance at 0.59 THz, there are two new transmission dips which 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 increases.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. 2a.The two new resonances have finite linewidth induced by the leakage of ideal BICs which 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 α is discussed in Fig. 2b.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 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]: 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] (2) The match of fitting results via Eq. 2 with the numerical spectra is shown in Fig. 3b 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) are shown in Fig. 4a.The spectral contrast ratio is defined as.
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. 3d, 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 (4) 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. 3d.The currents at 0.44 THz 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 observed in Fig. 3e.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 to 0.55 eV.With the increase of chemical potential, there is a blueshift of both two Q-BIC resonances.The spectral contrast ratio of modes Q-BIC-I and Q-BIC-II is 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 that our dual band Q-BIC metasurface is 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 where ф donates the phase change of the transmission spectrum.In Fig. 6a, we show the phase oscillates rapidly from −1.39 to −0.42 rad and from −2.41 to −1.98 rad as close to the Q-BIC-I and Q-BIC-II regions, respectively.The phase ф changes monotonously like usual plane wave while away from the Q-BIC region.Figure 6b shows the largest positive time delay of 23.8 ps and 14.5 ps, which represent the remarkable slow light effect and can be achieved by the modes 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. 7a 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 7b, d show the frequency shifts of the modes Q-BIC-I and Q-BIC-II corresponding to different refractive indexes, respectively.For mode Q-BIC-I in Fig. 7b, 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. 7d, 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] The maximum FOM is calculated to be 5.87 RIU −1 and 24.08 RIU −1 for the modes Q-BIC-I and Q-BIC-II, respectively.
In addition, the sensing performance of the proposed graphene metasurface is 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][28][29][30][31][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.

Conclusions
In summary, an all-graphene metasurface capable of simultaneously exciting dual Q-BICs is presented for the first time.The dual ideal BIC points are established by choosing periodic graphene DSRRs on a terahertz transparent SiO 2 substrate.This happens at the asymmetric parameter δ = 0. Dual Q-BICs are observed by varying the asymmetric parameter δ to break the in-plane C 2 rotational symmetry, which manifested as two new resonance dips in the transmission spectrum.We calculate the finite values of two Q-BICs' Q-factors, which follow the asymmetry degree as Q∝α −2 and diverge at the ideal BIC points.Furthermore, the Q-BICs' operating frequency range can be actively expanded via tuning graphene's chemical potential, and the dual Q-BICs are immune to the variation of incident angle.Meanwhile, a significant slow light effect is reported with time delay up to 23.8 ps.Our proposed Q-BIC metasurface can operate as a THz sensor.A maximum sensitivity of 267.5 GHz/RIU is obtained with a FOM of 24.08 RIU −1 .We envision that the proposed graphene metasurface empowered by Q-BIC resonances will exhibit applications in actively tunable slow light and sensing devices and will promote the development of the class of graphene Q-BIC metasurfaces.
Author Contribution Zhong Huang processed the data and wrote the manuscript.Qing Chen and Ji Ding give the support of numerical calculation in software.Chaojun Tang, Zhendong Yan, and Jing Chen helped with analysis and revised the manuscript.Xiaoxiao Cao revised the manuscript.All authors reviewed the paper and approved the final version to be published.
Institute of AdvancedPhotonics Technology, College of Electronic and Optical Engineering & College of Microelectronics, Nanjing University of Posts and Telecommunications, Nanjing 210023, China

Fig. 1 a
Fig. 1 a Schematic diagram of proposed graphene dual Q-BIC metasurface.b Top view of a unit cell of the graphene metasurface

Fig. 2 a
Fig. 2 a Transmission spectra of the periodic graphene DSRRs with different asymmetric parameters δ. b The relationship between the asymmetry degree α and Q-factors of Q-BIC-I and II

Fig. 3 a
Fig. 3 a Analysis between symmetric and asymmetric graphene DSRR (δ = 0 and δ = 4 μm) transmission spectra.b Fano fitting curves of Q-BIC-I and Q-BIC-II with asymmetric parameter δ = 4 μm.c Electric field distribution and surface current flow distribution (the white arrows) at 0.59 THz with asymmetric parameter δ = 0. d, e Electric field distribution and current flow distribution (the white arrows) of mode Q-BIC-I at 0.44 THz and mode Q-BIC-II at 1.19 THz with asymmetric parameter δ = 4 μm, respectively

Table 1
Comparison of sensing performance of the proposed Q-BIC metasurface