Terahertz electromagnetic waves hold unique opportunities for various important applications, including terahertz communications [1], non-destructive measurements [2], and security imaging [3]. Metamaterials, artificial structural materials, are used for arbitrary manipulation of terahertz waves, taking advantage of its unique electromagnetic properties. Essential functionalities achieved include holograms [4–6], generalized Snell’s law [7], metalenses [8, 9], and wavefront controlling [10, 11]. Furthermore, metamaterials have become instrumental for observing non-linear phenomena, providing a new platform for development of non-linear photonics [12]. Compared to the non-linear bulk crystals, metamaterials can not only enhance the localized density of states but also relax the necessity of phase-matching constraints [13, 14], leading to stronger non-linear processes, which are very interesting for practical applications.

Resonance response, associated with strong localized field confinement, has the potential to improve the non-linear responses. Recently, the bound states in the continuum (BICs), which can completely confine the electromagnetic wave without radiation leakage, have been introduced as strong resonances in metamaterials [15]. The most common approach is to realize symmetry-protected BICs in metamaterials with in-plane inversion symmetry [16]. By assigning the bound state and the continuum to different symmetries, the leakage of the bound state is forbidden due to symmetry incompatibility. True BICs are resonances associated with an infinite quality-factor (Q-factor), considered embedded trapped modes. However, the infinite Q-factor BICs are not available under external excitation. In practice, with the structural perturbations, the BICs can transform into quasi-BICs with a finite Q-factor [17]. For metamaterials supporting symmetry-protected BICs, the distinct quasi-BICs are realized by breaking the symmetry [16]. Although BIC metamaterials are widely used to boost the non-linear response in the optical range [18–20], the Q-factor of quasi-BICs at terahertz frequencies is several orders of magnitude lower than that at optical bands, which limits the research on the terahertz non-linear response. Exploring the efficient terahertz non-linear response remains a basic challenge.

In the paper, we proposed a graphene grating-like metamaterial with BICs for enhancing the terahertz THG response. By symmetry breaking of the metamaterial, the symmetry-protected BICs can make the transition to quasi-BICs. It enables high Q-factor and strong localized field enhancement, improving the light–matter interaction in monolayer graphene. By exploiting the large third-order non-linearity of graphene, the THG non-linear process with lower pump power at quasi-BICs is significantly boosted. Our study reveals that the quasi-BICs can significantly boost the THG in the graphene grating-like metamaterial, which offers an approach to achieve practical non-linear devices for on-chip terahertz non-linear conversion applications.

The schematic of the proposed graphene grating-like metamaterial is illustrated in Figure 1A, which contains two monolayer graphene strips, a monolayer graphene nanoribbon, and a metal ground layer separated by the glass spacer layers (where the refractive index of glass is 1.6 [21]). The metamaterial is periodic in the direction of x-axis and considered infinite in the direction of y-axis. As shown in Figure 1B, the structural parameters of the metamaterial include p = 300 nm, w ₁ = 100 nm, w ₂ = 90 nm, d = 50 nm, h ₁ = 12 nm, and h ₂ = 1,350 nm. The thickness of the metal ground layer is greater than the skin depth to reflect the transmitted waves. The Fermi levels of the monolayer graphene strips and the monolayer graphene nanoribbon are E _f1 and E _f2, respectively At terahertz frequencies, the linear conductivity of graphene can be described by the Drude model as [22, 23] σ g = i e 2 E f π ℏ 2 ω + i τ − 1 , (1) where e is the electron charge, E _f is the graphene Fermi level, ℏ is the reduced Planck’s constant, ω is the angular frequency of the incident waves, and τ is the carrier relaxation time, which can conform to the following equation τ = μE _f /(ev _f ²). The Fermi velocity v _f is 10⁶ m/s, and the carrier mobility μ is 1 m²/(V·s), which can be obtained in experiment [24]. The linear conductivity is a complex number, and the imaginary part represents the loss of graphene. In the simulation, we neglected the loss of glass, and the loss of the metamaterial is primarily derived from graphene. The third-order non-linear conductivity of graphene is expressed as [25, 26] σ 3 = i σ 0 ℏ v f e 2 48 π ℏ ω 4 T ℏ ω 2 E f , (2) where σ ₀ = e²/4ℏ, T(x) = 17G(x) - 64G(2x) + 45G(3x), and G(x) = In|(1 + x)/(1 - x)| + iπθ(|x| - 1). θ(z) is the Heaviside step function. The graphene is modeled to be a non-linear surface current, which is calculated by the formula [27]. J = σ g E T H + σ 3 E F F 3 , (3) where E _FF and E _TH are the electric fields of the fundamental frequency and the third-harmonic wave, respectively. The metamaterial is normally illuminated by a transverse magnetic-polarized wave with an input intensity of 100 kW/cm² along the direction of z-axis. The finite element method (FEM) is utilized to compute the linear response and third-order non-linear response in the frequency domain. The perfectly matched layer conditions are employed in the direction of z-axis, and the period boundary conditions are employed in the direction of x-axis.

The linear response for the graphene grating-like metamaterial is first studied. The proposed metamaterial produces symmetry-protected BICs, which can transform into quasi-BICs with high Q-factor and strong localized field enhancement through breaking the symmetry of the metamaterial. The asymmetry property of the metamaterial is controlled by an asymmetric factor defined as α = |w ₁-w ₂|/w ₁, where w ₁ is fixed. For the metamaterial structure with w ₂ = 100 nm, the system can be considered a Fabry–Pérot cavity circular structure, producing a Fabry–Pérot resonance mode, which is explained by [28] ω F P R = c K 2 w n M , (4) where c is the velocity of light in vacuum, w is the width of the graphene strip, and n _M is the effective refractive index of Fabry–Pérot resonance modes.

The absorption spectra as a function of frequency with different widths w ₂ are shown in Figure 2A. When asymmetric factor α = 0 (w ₂ = 100 nm), the proposed metamaterial meets inversion symmetry. Due to symmetry mismatch between the bound state and the continuum, the leakage at the bound state is forbidden at Γ point. Therefore, the left resonance disappears completely, forming symmetry-protected BICs. By varying the width of graphene 2 (w ₂) to break the symmetry, the symmetry-protected BICs are perturbed and leak to free space, leading to quasi-BICs with high Q-factors. Notably, the quasi-BICs are essentially Fabry–Pérot resonance modes. The quasi-BICs move from high to low frequency, with the decreasing width w ₂ explained by Eq. 4. The temporal coupled mode theory is suitable to analyze the coupling in a single-mode optical resonator with multiple ports, and the absorption spectra can be represented by [29]. A = 1 − R = 4 δ γ ω − ω 0 2 + δ + γ 2 , (5) where δ and γ represent the intrinsic loss and external leakage rates, respectively. ω₀ is the resonance angular frequency. Figure 2B shows the Q-factor and the maximum electric field enhancement of the metamaterial versus the asymmetric factor α. The Q-factor and the maximum field enhancement increase gradually as α decreases from 0.2 to 0.04 (corresponding to w ₂ from 80 to 96 nm). The Q-factor is 72 at the metamaterial of α = 0.2 but can be as high as 109 at α = 0.1 and even 135 at α = 0.04 at the quasi-BICs. Generally, the quasi-BICs have a Q-factor of 10⁴ or more at optical bands. At terahertz spectra, the Q-factor of plasmonic metamaterials is significantly reduced by several orders of magnitude compared to optical bands due to inherent ohmic and radiation losses. The order of magnitude of the Q-factor for the quasi-BICs is comparable to that previously reported [17, 30, 31]. Moreover, the quasi-BICs are not only associated with high Q-factors but also with remarkably strong localized field enhancement, which can significantly enhance the non-linear response of the system.

The quasi-BICs are localized by the Fabry–Pérot cavities, which can be regulated by modulating the graphene Fermi level. In the simulation, the relaxation time is dependent on the Fermi level and changing the Fermi level results in a corresponding variation in the relaxation time. Thus, the effects of time and Fermi level on conversion efficiency are discussed together. The Fermi level can be manipulated by the gate voltage, which can be written as follows [32]: E f = ℏ v F π ε 0 ε 1 V e h 1 , (6) where ε ₀ and ε ₁ are the vacuum permittivity and static permittivity of the glass, respectively. Figure 3A displays the absorption of the quasi-BICs with different Fermi levels E _f1. As the Fermi level decreases, E _f1, the resonance frequency of the quasi-BICs, shows a distinct redshift. It can be attributed to the variation in graphene conductivity resulting from the Fermi level. According to Eq. 2, the metamaterial with low Fermi level can achieve strong non-linear response. However, the absorption of the quasi-BICs is below 20% for the Fermi level E _f1 of 0.2 eV, which is negative for boosting the non-linear response. To improve the absorption of the quasi-BICs for the Fermi level E _f1 = 0.2 eV, the effect of the distance h ₁ on the absorption is calculated as shown in Figure 3B. Clearly, the absorption at the quasi-BICs increases as distance h ₁ increases, and the absorption can reach 63% with distance h ₁ of 52 nm. Moreover, the resonance frequency of the quasi-BICs increases as the distance h ₁ increases, and it can be explained by the mechanism of the Fabry–Pérot resonance mode. Graphene plasmons can be accompanied by damping effects due to inelastic scattering with phonons and elastic carrier scattering processes when the plasmon energy exceeds 0.2 eV [33]. The resonance frequency of the quasi-BICs is less than 20 THz, corresponding to the plasmon energy of 0.08 eV. Therefore, damping effects in the proposed system can be neglected.

We next investigate the THG response for the graphene grating-like metamaterial. Utilization of graphene to improve the functionality of absorbers and non-linear effects has been extensively reported in previous works [34–38]. However, the non-linear effects normally require tremendous input radiation intensities to be excited, rendering their practical applications an enormous challenge. The proposed metamaterial with BICs significantly confines the electromagnetic power in subwavelength range, which results in a high Q-factor and strong localized field near the non-linear graphene, enhancing the non-linear effects. Meanwhile, the requirement of complex phase-matching conditions for the proposed metamaterial with ultrathin thickness can be relaxed than for the non-linear body crystals. The fundamental frequency of the incident radiation for the non-linear metamaterial is aligned with the resonance frequency of the quasi-BICs. To represent the power conversion from the fundamental frequency to the third-harmonic wave, the THG conversion efficiency is calculated as η = P _TH(3ω)/P _FF(ω) [39], where P _TH(3ω) is the output THG power for the third-harmonic frequency derived from the integration of the THG power outflow on the top port of the metamaterial and P _FF(ω) is the incident power for the fundamental frequency. The parameters are assumed to be p = 300 nm, w ₁ = 100 nm, w ₂ = 90 nm, h ₁ = 52 nm, h ₂ = 1,350 nm, E _f1 = 0.2 eV, and E _f2 = 1 eV.

Figure 4A depicts the THG conversion efficiency of the graphene grating-like metamaterial versus the fundamental frequency. Obviously, the THG conversion efficiency exhibits a remarkable peak at the resonance frequency of 18.7 THz for the quasi-BICs. The lower or higher fundamental frequencies reduces the THG conversion efficiency, as explained using Eq. 3. With the input intensity around 100 kW/cm², the graphene grating-like metamaterial exhibits a THG conversion efficiency of 2.8% at the quasi-BICs. The electric-field-enhancement distributions for the graphene grating-like metamaterial at the fundamental frequency and the third-harmonic wave are illustrated in the insets of Figure 4A. The strong localized field enhancement is achieved both at the fundamental frequency and at the third-harmonic wave, which is beneficial for improving the THG response. Figure 4B displays the incident intensity of fundamental frequency dependence of THG conversion efficiency for the graphene grating-like metamaterial. The fundamental frequency is fixed at a resonance frequency of 18.7 THz for the quasi-BICs, and the THG conversion efficiency is significantly enhanced with increasing incident intensity, revealing a squared trend. As shown in the inset of Figure 4B, the corresponding log plot displays a linear relationship with a slope of ∼2 between the THG conversion efficiency and incident intensity, as expected for the THG process.

Figure 5A displays the THG conversion efficiency for the graphene grating-like metamaterial at Fermi level E _f1 = 0.2, 0.4, and 0.6 eV. As shown in the figure, the THG conversion efficiency of the quasi-BICs increases as the Fermi level E _f1 decreases. At the quasi-BICs, the THG conversion efficiency under the Fermi level E _f1 decreased from 0.6 to 0.2 eV and improved from 0.6% to 2.8%, which is enhanced by more than four times. This significant enhancement can be explained by the increase in third-order non-linear conductivity as the graphene Fermi level decreases, as given by Eq. 2. Particularly, the variation in the Fermi level can have an effect on the carrier mobility. The maximum conversion efficiency in the proposed system is obtained with a fixed Fermi level of 0.2 eV, so the variation of the Fermi level hardly affects the carrier mobility. The THG conversion efficiency for the graphene grating-like metamaterial with different distances h ₁ is plotted in Figure 5B. It is obvious that the THG conversion efficiency of the quasi-BICs increases as the distance h ₁ increases, which is caused by the increase in the absorption of the quasi-BICs at the fundamental frequency with the increase in the distance h ₁, as shown in Figure 3B.

Based on the aforementioned results, we now investigate how to further enhance the non-linear response of the graphene grating-like metamaterial. As shown in Figure 4B, the THG response has not achieved saturation at the asymmetry factor α = 0.1 (w ₂ = 90 nm). The transition from BICs to quasi-BICs causes the resonance to decay, resulting in leakage of optical radiation to free space, which can reduce the light–matter interaction in non-linear graphene. To solve this problem, we calculate the THG conversion efficiency for the graphene grating-like metamaterial with different widths w ₂, as shown in Figures 6A. Figures 6B displays the variation in the THG conversion efficiency at the quasi-BICs with asymmetry factor α ranging from 0.04 to 0.2 (corresponding width w ₂ from 96 to 80 nm). Remarkably, the THG conversion efficiency can reach 3.1% for the asymmetry factor α of 0.04, which can be largely enhanced by 10% in contrast to the result obtained with the asymmetry factor α of 0.1. The efficiency of the proposed metamaterial is significantly boosted than that of previous ones reported both in simulations [40–46] and experiments [12, 47, 48]. Notably, the THG conversion efficiency at the quasi-BICs can be enhanced at a higher level with the asymmetry factor α close to 0.

In summary, we theoretically demonstrated a graphene grating-like metamaterial with BICs for enhancing the terahertz THG response. The quasi-BICs with accompanying high Q-factor and strong localized field enhancement can be achieved through symmetry-breaking of metamaterials, producing the strong terahertz THG response. By adjusting the graphene Fermi level and structural parameters, the THG conversion efficiency is further improved. Under the input intensity 100 kW/cm², a THG conversion efficiency of 2.8% is attained with the asymmetric factor α of 0.1. Moreover, THG conversion efficiency exceeds 3.1% at α of 0.04. The proposed efficiency non-linear metamaterial indicates great promise for on-chip terahertz non-linear conversion applications.