Terahertz (THz) radiation, electromagnetic radiation with frequencies spanning from 0.1 to 10 THz, holds a wealth of spectral information ripe for exploitation in molecular fingerprint spectroscopy [1,2]. Recently, the intersection of terahertz technology and metasurface technology has been a focus of research, yielding significant practical implementations in communication [3], security screening [4], and biosensing [5]. Metamaterials, artificial composite materials first proposed by Walser in 2001 [6], are often constituted by periodically arranged, sub-wavelength structural units of diverse shapes, enabling unique and custom properties. As the field evolved, the creation of negative refractive index materials through metamaterials was not the terminal achievement. Scientists have manipulated the geometric shape and constituent materials of structures to regulate properties such as polarization, phase, and amplitude of electromagnetic waves, catalyzing the rapid development of novel devices for polarization control, beam shaping, and holography [7–10]. In the present state of affairs, terahertz metasurfaces represent a burgeoning research area, pointing to a promising future in diverse applications. Research on metasurfaces can be categorized into two types: metallic metasurfaces [11–13] and dielectric metasurfaces [14–16]. While the inherent Ohmic losses of metals restrict their applications [17], dielectric metasurfaces, boasting superior thermal properties and compatibility with complementary metal oxide semiconductor (CMOS) processes, hold promise for low-loss photonics [18,19]. As a result, an increasing number of scholars are turning their attention to the study of dielectric metasurfaces. Metasurfaces can support various resonance modes, including Mie resonance [20], Fano resonance [21], and SPR resonance [22]. These resonance modes can produce localized field confinement, which enhances the interaction between light and matter, leading to improved performance parameters such as nonlinear conversion efficiency and sensing sensitivity [16]. The quality factor (Q) has become an important indicator in resonance mode analysis [23]. For sensing applications, a high Q resonance mode, indicating a narrow linewidth, signifies the capacity to detect slight frequency shifts.

In recent times, the phenomenon of bound states in the continuum (BICs) has garnered considerable interest in photonics, owing to their potential to attain low-loss, high Q resonances [24]. Originating from quantum mechanics, bound states in the continuum were first introduced by Von Neumann and Wigner in 1929 [25]. However, their relevance in the realm of optics was only recognized in 2008 [26]. In theoretical terms, BICs are entirely decoupled from incident radiation, lacking any outlets for radiation escape, resulting in infinite Q factors [27,28]. These are also identified as trapped modes, embedded eigenvalues, or dark states [29]. Metamaterial applications involve the classification of BICs into two categories: symmetry-protected BICs and accidental BICs, characterized by their design and excitation principles. Symmetry-protected BICs are achieved by tilting the incident angle or breaking the structural symmetry [30], while accidental BICs are unexpectedly generated by adjusting one or more structural parameters, allowing for continuous tuning of the system parameters while being decoupled from the radiation continuum [31]. In current research on BIC-based metasurfaces, silicon and relatively low-loss materials such as quartz are frequently used [32–34]. As research advances, other low-loss dielectric materials like lithium tantalate [35,36] and lithium niobate are being explored more extensively.

Lithium niobate, often referred to as “optical silicon,” is a significant material in photonics due to its remarkable acousto-optic, piezoelectric, and second-order optical nonlinearity properties [37]. It possesses stable properties and low losses, which have attracted extensive research in fields such as optical fiber communication [38], quantum communication [39], and microwave photonics [40]. There are two primary applications of LiNbO3-based metamaterials based on BICs: nonlinear applications [41–43] and electro-optic modulation [44–46]. Hence, while we acknowledge and discuss the already established uses of lithium niobate metasurfaces, our study seeks to shed light on their untapped potential in the field of sensing, pushing the boundaries of what is currently achievable in this area. Current terahertz tunable metasurfaces predominantly rely on graphene [47–49] and phase-change materials such as vanadium dioxide [50]. However, lithium niobate materials provide a wider electro-optic modulation range, faster response times (enabling nanosecond-level modulation), and additional benefits, including stable properties and low losses.

In this paper, we innovatively introduce BICs into the sensing application of lithium niobate metasurfaces. By utilizing a high refractive index, and stable chemical and physical properties of lithium niobate materials, we present two strategies. These encompass disrupting both physical structural symmetry and material refractive index symmetry, which culminate in ultra-high Q factor resonances excited by symmetry-protected BICs within the terahertz range. The presented simulation results analyze the performance of the resonator, with a particular focus on refractive index sensing. This structure exhibits a refractive index sensitivity of up to 947 GHz/RIU and a Q factor of 6.56×104. The findings reveal an effective method for achieving ultra-high sensitivity sensing applications of lithium niobate in the terahertz frequency range. Additionally, manipulating the electro-optic properties of lithium niobate to control the frequency shift of the transmission peak is a promising strategy for the optimization of metasurface sensors. This approach provides new insights into the design and manufacturing of terahertz metasurface sensors.

The schematic structure of the terahertz lithium niobate metasurface is shown in Fig. 1(a). The metasurface is composed of dimer cluster periodic arrays arranged in a square lattice, with cylinders made of LiNbO3 and a substrate made of quartz substrate. The proposed unit structure, including two lithium niobate cylinders of identical radius, is delineated in Fig. 1(b) along with its geometric parameters. The lattice constants are Pz=60μm and Py=45μm, the radius of the cylinders is r=11μm, the height of the cylinders is h=15μm, and the separation between the two cylinders is Pz/2, which is 30 μm. To highlight the structural asymmetry, a top view of the unit cell on the y-z plane is presented in Fig. 1(c). A geometric parameter d is introduced to represent the distance that the right cylinder moves along the z-axis. When d=0, the right cylinder remains stationary, resulting in a unit cell that is symmetric on the y-z plane. When d≠0, it indicates symmetry breaking, where the right cylinder of the dimer cluster is displaced by a distance of d along the z-axis while the position of the left cylinder remains unchanged. The refractive index of LiNbO3 is set to 5.2 within the studied terahertz frequency range [51], while the refractive index of quartz is identified as 2. To analyze the spectral responses and electromagnetic properties, the frequency-domain solver in CST Studio Suite 2022 is used. The unit cell boundary conditions are applied in both the z direction and y direction, while open boundary conditions are used in the x direction during the simulation, and a tetrahedral mesh type is selected. The plane wave with y-polarization and propagation along the x-axis is used as the excitation field.

Lithium niobate can be etched to a scale ranging from tens of micrometers using capacitively coupled plasma reactive ion etching (CCP-RIE) [52] and inductively coupled plasma reactive ion etching (ICP-RIE) [53], to hundreds of micrometers using a femtosecond laser [54,55]. Figures 2(a) and 2(b) depict simplified processing flowcharts. The proposed structure is manufacturable with standard micro-nanofabrication equipment, and its performance can be assessed using terahertz time-domain spectroscopy (THz-TDS) [56]. During practical etching processes, the sidewalls of the lithium niobate structures obtained exhibit certain slope angles. To address this, we implemented a conical dimer metasurface to simulate the impact of this experimental error on sensing performance. The slope angle of the sidewall is denoted as “θ” and used to calculate the Q factor and figure of merit (FOM) of the structure with step sizes of 2° in the range of 74°–80°. Encouragingly, the simulation results of the conical dimer metasurface also reveal quasi-BIC transmission peaks. Furthermore, You et al. [57] similarly underscored the limited influence of nanodisk geometry on BIC within dimer metasurface structures. The results indicate that the Q factor of the conical dimer metasurface structure reaches a minimum of 3.55×104, with the lowest FOM being 2656. These values confirm the structure’s robust and significant sensing performance.

A true BIC is a mathematical entity characterized by an infinite Q factor and vanishing resonant width. Breaking the symmetry of the structure converts the BIC into a quasi-BIC, which is also referred to as a supercavity mode. In this state, both the resonant width and Q factor become finite [58]. We begin our analysis by examining the transmission spectrum of the structure, and the results indicate that the resonant peak is absent in the transmission spectrum when the geometric parameter d is adjusted to 0 μm in Fig. 3(a). This is due to the non-radiative feature of the BIC generated in the symmetric structure. As depicted in Fig. 3(b), setting the geometric parameter d to 2 μm disrupts the symmetry of the crystal cell, enabling the establishment of a radiative channel that connects the BIC to free space. A sharp Fano feature is observed at 4.144 THz in the transmission spectrum, representing a quasi-BIC response. To analyze the transmission spectrum qualitatively, the Fano equation is utilized to fit it: TFano=|a1+ia2+bω−ω0+iγ|2.(1)

Here, a1,a2, and b denote numerical constants, γ represents the entire damping rate, and ω0 represents the center frequency of the Fano response. The Q factor is calculated using the following formula: Qrad=ω0/(2γ) [59]. The excellent agreement between the numerically computed spectrum and the fitted spectrum (red dashed line) depicted in Fig. 3(c) confirms a high level of accuracy. Thus, the Q factor of the Fano resonance at 4.144 THz is established as 6.56×104. To reveal the existence of BICs in cylindrical dimer metasurface structures, we employed a rigorous computational approach based on band structure analysis. The results were obtained using the commercial finite element software COMSOL Multiphysics, utilizing lateral Floquet periodic boundary conditions. All numerical models are built with 3D structures. In Fig. 3(d), curves composed of solid circles depict dispersion curves along the X′Γ and ΓX directions within the structure. The eigenmode with an infinite Q factor at the Γ point exhibits a frequency of 4.136 THz [60].

To delve deeper into the physical mechanisms driving the resonance, we performed a multipole decomposition of the induced current density, extracted in a Cartesian-coordinate system. This approach is a common method used to elucidate the physical mechanisms of BICs [61]. In the calculation, we disregarded higher-order components with negligible contributions, focusing on five terms: the toroidal dipole (TD), magnetic quadrupole (MQ), electric quadrupole (EQ), magnetic dipole (MD), and electric dipole (ED). Figure 4(a) illustrates that the dominant factor contributing to the resonance is the MQ resonance, followed by the TD resonance. This is due to the fact that the presence of MQ and TD in the multipole decomposition analysis often leads to a pair of MDs oriented in opposite directions [62]. The EQ also plays a role in the resonance response, while the MD and ED are significantly attenuated. Based on the aforementioned findings, it can be inferred that the observed sharp resonance is attributable to the MQ resonance. The prominent existence of TD indicates that this resonance simultaneously exhibits the advantages linked with TD resonances, including high Q characteristics and enhanced nonlinear responses. This result aligns with prior multipole decomposition findings for cylindrical dimer [57,63] and tetramer cluster [35] metasurface structures. At the same time, the prevalence of MQ resonance concentrates the electric field enhancement in the gaps of the structure rather than internally. The symmetric distribution of the electric field in the y-z plane around the center of the cylindrical dimer in Fig. 4(b) further confirms the presence of the MQ resonance at the resonant frequency, and the enhancement of the field inside the cylindrical dimer cluster is greater than that between adjacent structures. This phenomenon facilitates the enhancement of light–matter interactions, proving advantageous for subsequent sensing applications [64]. As depicted in Fig. 4(c), the magnetic field distribution in the y-z plane points to the existence of a strong magnetic vortex current within the cylinder, which generates a robust resonance following the right-hand rule around the center of the dimer.

To gain a deeper understanding of the MQ resonance, we computed the transmission spectra while varying the geometric parameter d within the range of 0–4 μm. The results, as illustrated in Fig. 5(a), reveal that when d=0, the resonance disappears from the transmission spectrum, suggesting an absence of energy dissipation from the bound state into the surrounding free space. As the geometric parameter d increases, the linewidth of the resonance broadens and experiences a blue shift. This indicates that when d≠0, the symmetry of the structure is broken, and the symmetry-protected BIC transitions to a quasi-BIC. The radiation channel of the nanostructure expands, resulting in an increase in radiation energy leakage [65]. We define d>0 to mean that the right cylinder is closer to the y-axis, while d<0 means that the right cylinder is farther away from the y-axis. The correlation between the Q factor of the MQ resonance and the asymmetric parameter α=2d/Pz was also examined, as depicted in Fig. 5(b). The Q factor is determined by Q=f0/FWHM, where f0 signifies the resonant frequency, and FWHM is the full width at half maximum. The Q factor and asymmetric parameter follow the relationship Qrad∝α−2, implying that this Fano resonance, stimulated by the BIC, can have its frequency response’s position and width adjusted through the manipulation of the asymmetric parameter [66]. The excitation of this MQ quasi-BIC resonance with an ultra-high Q factor is associated with a robust near-field enhancement. Berté et al. [67] experimentally validated the attainment of BICs in metasurfaces by breaking the in-plane symmetry in the permittivity of the comprising materials. Their study also demonstrated the equivalence between permittivity-asymmetric qBICs and geometrically asymmetric qBICs. Building on these findings, we further studied the electro-optic modulation properties of LiNbO3 to explore a simpler adjustable mechanism to break the symmetry of the metasurface. By applying an external voltage to change the refractive index of the material, the permittivity symmetry of the structure can be further disrupted, thus achieving high Q factor quasi-BIC resonances. In this structure, we use x-cut LiNbO3. Along the x and y axes, the refractive index no corresponds to the ordinary light (o-light) refractive index, while along the z-axis, the refractive index ne corresponds to the extraordinary light (e-light) refractive index [37]. Notably, the refractive index change is at its maximum when the direction of the electric field aligns with the z-axis of the LiNbO3 crystal. Therefore, in the construction of metal electrodes corresponding to the metasurface, we opt to apply an external voltage along the z-axis to achieve maximum electro-optic modulation. When an electric field Ez exists between the two electrodes, the change in the refractive index along the z direction is Δnz=−ne3γ33Ez/2, where γ33 is a component of the LiNbO3 electro-optic tensor matrix. The electric field intensity Ez=V/G, where V is the driving voltage, and G is the gap between the electrodes; hence Δn can be represented as Δn=ne32γ33VG.(2)

From this, it is clear that the change in the refractive index of LiNbO3 material is induced by the externally applied voltage. Moreover, for a given material and electrode structure, the magnitude of this change is entirely controlled by the external voltage [68]. For the study, we simulated the transmission spectra of LiNbO3 metasurface structures under various external voltage conditions, while preserving the structural symmetry by keeping the structure parameter d at zero. We applied the external voltage to the outer sidewall of the LiNbO3 cylinder on the right in the z-axis direction, while keeping the ne_left refractive index of the LiNbO3 material on the left cylindrical side constant at 5.2, as illustrated in the inset of Fig. 5(d). The external voltage was varied in steps of 100 V within the range of 0–400 V (with the refractive index ne_right of LiNbO3 changing between 5.2 and 5.24). By directly modulating the refractive index ne_right of the LiNbO3 material, we simulated the effect of the applied external voltage on the structure’s transmission spectrum. As illustrated in Fig. 5(c), when ne_right=ne_left, the resonance disappears from the transmission spectrum. As the external voltage V increases, consequently resulting in a higher ne_right, it leads to an expansion of the linewidth for the transmission peak. This suggests that when ne_right≠ne_left, the symmetry of the structure is disrupted, transitioning the symmetry-protected BIC to quasi-BIC. Contrary to the conventional methodology of disrupting structural symmetry by modifying geometric parameters, a substantial red shift of the transmission peak is observed when the refractive index symmetry of the structure is altered. This phenomenon arises as modifications in the refractive index parameters of the structure lead to changes in the system’s optical path length, thereby adjusting the resonant frequency of the system. Consequently, when changes in the refractive index of the structure occur, the distribution and propagation speed of the light field within the medium are altered, which in turn changes the electromagnetic mode of the system. This has implications for the resonance condition of the system, including the resonant frequency and quality factor. In this instance, we define the asymmetric parameter β=(ne_right−ne_left)/ne_left, and concurrently investigate its correlation with the Q factor, as illustrated in Fig. 5(d). Despite the Q factors being an order of magnitude lower than when the structural parameter d is adjusted, they consistently exceed 103. It follows the relation Qrad∝β−2 with the asymmetry parameter β. The acquired Fano resonances, stimulated by an externally applied voltage altering the refractive index, are also generated by BIC. The positions and linewidths of their frequencies can be modulated by varying the external voltage, thereby bypassing the necessity to modify the physical structural parameters. This innovative approach sets the stage for realizing high Q quasi-BIC resonances.

To examine how the transmission spectrum varies with changes in the structure’s geometric parameters, transmission spectra were computed for the asymmetric metasurface structure using a range of different geometric parameters, while maintaining a constant asymmetric parameter of d=2μm. Except for the variable parameters shown in each figure, the other geometric parameters used in the analysis are identical to those presented in Fig. 1(a). The geometric parameters were adjusted in increments of 1 μm, and the corresponding results are illustrated in Fig. 6. According to Fig. 6(a), by adjusting the periodic constant along the z-axis Pz from 58 to 62 μm, the MQ resonances experience a red shift, and their linewidth broadens with increasing Pz. This is because as Pz increases, the effective refractive index of the structure also increases, resulting in a smaller effective area occupied by the structure within a unit cell. This decrease in the effective area reduces the interaction between adjacent cells, leading to a weakening of the overall influence between the structures. A similar effect is observed when the periodic constant along the y-axis Py is adjusted within the range of 43–47 μm, which also results in a red shift of the resonance peak, as illustrated in Fig. 6(b). As the radius r of the two lithium niobate cylinders increases from 9 to 13 μm, the resonance peak linewidth gradually narrows, and the frequency position experiences a red shift, as shown in Fig. 6(c). Similarly, when the thickness h of the cylindrical dimer increases from 13 to 17 μm, a red shift of the resonance peak is observed, and the linewidth gradually broadens due to the weakening of the coupling strength between the cylinders, as shown in Fig. 6(d). A comparison of the effects of various geometric parameters on the transmission spectra of the asymmetric metasurface structure reveals that the location of the MQ resonance is particularly responsive to changes in the structural parameters r and h of the cylinders, and the resonance peak position can be further fine-tuned by adjusting the structural parameters.

In the current research, the conventional method for tuning the transmission peak in the transmission spectrum involves altering the parameters of the metamaterial structure [69–71]. However, this approach often introduces significant manufacturing complexities, which consequently limit the widespread application of metasurface sensors [16]. To address this challenge, we have explored a simpler, adjustable mechanism based on the electro-optic modulation properties of LiNbO3, offering enhanced adaptability for sensing needs under different environmental conditions. Our findings reveal that the proposed LiNbO3 metasurface demonstrates remarkable tunability in the realm of optical metasurfaces. For the study, we simulated the transmission spectra of asymmetrical metasurface structures under a variety of external voltages, while maintaining the consistency of other structural parameters as per Fig. 1. The external voltage V was varied in steps of 20 V within the range of −40 to 40 V (with the refractive index ne of LiNbO3 changing between 5.196 and 5.204). We applied the external voltage along the z-axis direction to the outer sidewall of the two LiNbO3 cylinders in the unit structure, as illustrated in the inset of Fig. 7(a). By directly modulating the refractive index ne of the LiNbO3 material, we simulated the effect of the applied external voltage on the structure’s transmission spectrum. Figure 7(a) reveals that an increase in external voltage causes a blue shift in the transmission peak. Given the constancy of the structural asymmetry parameters, the spectral line shape and linewidth remain largely unaffected. This observation validates the robustness and controllability of the LiNbO3 electro-optic effect, allowing for precise tuning of the transmission peak. Further, the transmission spectra at different voltages were fitted with Fano profiles, yielding corresponding Q factors that demonstrate a consistent high-fitting performance. As depicted in Fig. 7(b), the Q factor exhibited an increase from 60,493 to 73,984, indicating that the structures maintained high sensing performance under varying external voltages. The adoption of LiNbO3 substantial electro-optic coefficient for refractive index control via an electric field facilitated the frequency tuning of the transmission peak. This strategy bypasses the complexities associated with adjusting physical structural parameters. This innovative tuning method not only mitigates the manufacturing complexity of terahertz metasurface sensors, but also provides the prospect of achieving a multitude of intricate optical functions.

In summary, we investigated and compiled literature on lithium niobate metasurfaces based on bound states in the continuum that have been published to date. We introduced an ultrasensitive terahertz sensor that utilizes MQ resonance controlled by BICs within an asymmetric lithium niobate metasurface structure. By designing cylinder clusters to form asymmetric clusters, we achieved ultra-high Q factor resonance excited by symmetry-protected BICs in the terahertz band, with a maximum Q factor of 6.56×104. Leveraging the exceptional electro-optical characteristics of lithium niobate, a method was proposed to disrupt its structural symmetry by altering the refractive index. This adjustment led to the realization of ultra-high Q factor resonance, reaching a peak of 3.64×104, which can be effectively regulated by voltage modulation. By adjusting structural parameters such as the spacing and radius of the two cylinders in the unit structure, we also obtained quasi-BICs with controllable Q factors. Furthermore, we studied the physical mechanism of this quasi-BIC using multipole decomposition methods. The electromagnetic field at the resonance is strongly enhanced and confined within the gap between the cylinder clusters. With a sensitivity of up to 947 GHz/RIU, the metasurface exhibits high sensitivity to variations in the refractive index of the surrounding medium, surpassing that of previously reported terahertz all-dielectric metasurface sensors. These findings provide a new approach to the development of ultrasensitive terahertz sensors. Considering the significant potential of lithium niobate, future research endeavors in this area could yield transformative advancements in sensor technology, thereby broadening the horizon of applications in various domains.