Artificially engineered materials, usually known as metamaterials, have exhibited novel electromagnetic phenomena which do not exist in natural materials, such as negative refractive index [1, 2], electromagnetic cloaking [3–6], sensing [7, 8], and near-perfect absorption [9, 10]. The elimination of absorption via quantum interference in an atomic medium is termed as electromagnetically induced transparency (EIT) [11–15] but now can be mimicked by non-quantum approaches. Recently, analogs of the EIT-like behavior, plasmon-induced transparency (PIT), and metamaterial-induced transparency (MIT) have attracted much attention in coupled resonators [16, 17], electric circuits [18], and plasmonic structures in the terahertz (THz) range [19–23]. These effects significantly modify the dispersive properties of an otherwise opaque medium, which can lead to slow light phenomena and enhance nonlinear effect in the light–matter interactions. The most frequently used theoretical methods to mimic those effects are the RLC equivalent circuit model and the Lorentz oscillator model. The equivalent circuit models, simplifying the metal structures to lumped RLC elements, efficiently represent their electrical performance for circuit simulation [24]. In our previous studies, we have calculated the induced currents within the split ring resonator (SRR) loops, electric dipoles at the split gaps, and the corresponding radiation spectra for both co- and cross-polarizations using the RLC circuit model [25]. On the other hand, the effect of external field on the metamaterials can also be equivalent to that on artificial atoms to drive the vibration of effective charge based on the Lorentz model. The coupling in this model can be treated as bright–dark mode coupling (only one bright mode resonator couples to the incident wave directly) or bright–bright mode coupling (all resonators interact with the external electric field). For the bright–dark mode, the PIT effect on opposite-directed U-shaped SRRs was calculated under different photoexcitation powers to modulate the amplitude of the broadband transparency window [26]. Yahiaouia et al. proposed a hybrid EIT metamaterial and implemented numerically to tune the EIT window by incorporating photosensitive silicon pads in the split gap region of the resonators [27]. The near-field coupling between two bright modes was reconfigured to realize the deep modulations of the EIT window in frequency position and transmission amplitude in the metamaterial consisting of a graphene cut wire resonator and two graphene closed ring resonators that enable the dynamic controlling transparency window [28]. Additionally, the Lorentz oscillator model has been applied to the three coupled resonators to study the interference property of the EIT metamaterial with a cut wire and two pairs of split-ring resonators [29], as well as the coupling tailored with photoexcitation between the square rings [30]. It is shown that the Lorentz model is more generalized and independent of the specific microstructure, which is especially suitable for the structures that cannot form a circuit, such as metallic bar array resonators. Moreover, few studies have been carried out to find out the relationship between the changed fitting parameters and the underlying physical coupling mechanism. The parameters have great impact on tuning and manipulating spectral line shapes, which are closely related to the sensing performance about the quality factors of the resonances in the coupled metamaterials. Metamaterial-based THz sensor can provide the sensitive detection of the analyte owing to their localized field enhancement properties; exhibiting the variation of the dielectric environment will induce the obvious change of sensing spectra involved with the resonant frequency and linewidth [31, 32]. Recent excellent studies show the apparent shifts of both the resonant frequency and depth can be observed in THz toroidal metasurface for sensitive distinction of lung cancer cells [33]. The EIT-like resonance can experience frequency and magnitude variations when the properties of analytes change, showing that the glioma cells can be distinguished directly based on this EIT resonance [34]. It is shown that the influence of geometrical and theoretical parameters on the spectral sensing ability needs to be clarified to further optimize the metamaterial sensor in the application.

Here, we have investigated the characteristics of THz spectra of three structures with simple bar arrays, especially the impact of the spacing between bars on the spectral evolution for double-bar and 3-bar structures. It is found in the experiment that the coupling intensity between bars increases initially and then decreases with the bar spacing, which indicates that not only intracoupling but also intercoupling together play important roles in electromagnetic response. We have calculated and analyzed systematically using the Lorentz oscillator models from a single resonator to three resonators, and obtained the agreeable results compared with experimental data. Further calculation was performed to study the influence of the fitting parameters, damping rate, and coupling coefficient on the transmission spectra. By simulating the spectral evolution with the overlayers of different refractive indices, the relation between the sensing ability and the geometrical and theoretical parameters is discussed to explore the optimized metamaterial sensor. Our studies could help to understand the coupling mechanisms and provide prediction and guidance in THz modulation and sensing devices.

We designed three single-bar resonators as shown in Figure 1A, termed as M1, M2, and M3, with length L of 36, 46, and 56 μm, respectively. Then we combined configurations between two or three of them to form new structures of M12 (combination of M1 and M2) and M123 (combination of M1, M2, and M3), as shown in Figures 1B,C. The spacing parameter S ₁ of M12 between 2 bars is a variable and set to be 2, 12, 24, and 32 μm, respectively. The spacing parameter S ₂ of M123 between every 2 bars is set to be 2, 5, 14, and 18 μm, respectively. Other structure parameters in the unit cell are illustrated in Figure 1. Numerical simulations of spectral responses, surface currents, and electric field distributions of those samples were performed using the finite integration method.

Periodic golden arrays with a thickness of 200 nm were fabricated on the 620-μm-thick semi-insulating GaAs substrate using the photolithography method. The laser source is a Ti:sapphire regenerative amplifier delivering ultrashort optical pulses with duration of 50 fs and a central wavelength of 800 nm at a repetition rate of 1 kHz. The THz radiation is normally incident to the sample plane, and the wave polarization is along the bar direction. The transmitted terahertz wave is detected by free-space electro-optic sampling in a 1-mm-thick <110> ZnTe crystal with the probe pulse [35]. The experiment was carried out at room temperature, and the GaAs substrate was served as a reference.

To elucidate the underlying mechanism of the spectral response and coupling effect, and to calculate resonant dips and the MIT transmission spectra, a coupled Lorentz model combined with the effective medium theory is adopted. We first consider an oscillating motion in a single oscillator. If a particle moves under an external driving force from incident field E(ω) = E ₀e ^iωt , then the motion can be given by [36]. x ¨ + γ   x ˙ + ω 0 2 x = gE ( ω ) , (1) where x is the displacement of the oscillator, ω ₀ is the natural frequency of the oscillator, γ is the damping rate, and g is the geometric parameter indicating the strength of each oscillator coupled to the incident field and related to effective charges and effective masses.

If another oscillator is added in this system, we can obtain the dynamics of a pair of oscillators coupled under the incident field termed as the bright–bright mode. The coupled equations can be analytically described by [27, 28]. x ¨ 1 + γ 1 x ˙ 1 + ω 1 2 x 1 + κ x 2 = g 1 E ( ω ), (2) x ¨ 2 + γ 2 x ˙ 2 + ω 2 2 x 2 + κ x 1 = g 2 E ( ω ), (3) where x ₁, x ₂, γ ₁, γ ₂, ω ₁, and ω ₂ are the displacements, the damping rates, and resonance angular frequencies of the resonance modes of each oscillator, respectively. κ is the coupling coefficient between two oscillators, while g ₁ and g ₂ are the geometric parameters. When g ₂ = 0, the coupling behavior is usually known as the bright–dark mode, showing that only one oscillator interacts with the external field.

Similarly, considering the three oscillators, the resonances and their mutual coupling behavior can be described through the following equations [30]: x ¨ 1 + γ 1 x ˙ 1 + ω 1 2 x 1 + κ 12 x 2 = g 1 E ( ω ), (4) x ¨ 2 + γ 2 x ˙ 2 + ω 2 2 x 2 + κ 12 x 1 + κ 23 x 3 = g 2 E ( ω ), (5) x ¨ 3 + γ 3 x ˙ 3 + ω 3 2 x 3 + κ 23 x 2 = g 3 E ( ω ), (6) where x ₁, x ₂, x ₃, γ ₁, γ ₂, γ ₃, ω ₁, ω ₂, and ω ₃ are the displacements, damping rates, and resonance angular frequencies of M1, M2, and M3, respectively. κ ₁₂ and κ ₂₃ are the coupling coefficients of M1–M2 and M2–M3, respectively. g ₁, g ₂, and g ₃ are the geometric parameters. By solving Eqs 4–6, we can obtain x 1 ( ω ) = κ 23 2 g 1 − κ 12 κ 23 g 3 − C 2 C 3 g 1 + C 3 κ 12 g 2 C 3 κ 12 2 + C 1 κ 23 2 − C 1 C 2 C 3 E ( ω ), (7) x 2 ( ω ) = C 3 κ 12 g 1 − C 1 C 3 g 2 + C 1 κ 23 g 3 C 3 κ 12 2 + C 1 κ 23 2 − C 1 C 2 C 3 E ( ω ), (8) x 3 ( ω ) = κ 12 2 g 3 − κ 12 κ 23 g 1 − C 1 C 2 g 3 + C 1 κ 12 g 2 C 3 κ 12 2 + C 1 κ 23 2 − C 1 C 2 C 3 E ( ω ), (9) where C _j = −ω ²+iωγ _j +ω _j ² (j = 1, 2, 3). In addition, the electromagnetic susceptibility, which relates the intensity of polarization P(ω) of the oscillator to the strength of incoming electric field, is expressed as χ _e(ω) = P(ω)/ε ₀ E(ω), where ε ₀ is the permittivity of vacuum. Then the susceptibility of the samples can be expressed using γ _i , ω _i ,   g i (i = 1, 2, 3), κ ₁₂, and κ ₂₃. Due to it being thin enough, the metamaterial layer can be considered an effective medium layer around air and substrate. We can use the Fabry–Perot interference transmission equation and the Fresnel formula to calculate the transmittance [30]: lim d → 0 | t ˜ ( ω ) | = lim d → 0 | t ˜ T M t ˜ a i r − s a p | = | c ( 1 + n s u b ) c ( 1 + n s u b ) − i ω χ ˜ e | , (10) where t ˜ T M is the transmission of the effective medium layer, t ˜ a i r − s a p is the transmission at the air–substrate interface, n _sub = 3.57 is the refractive index of the GaAs substrate from our experiment, and c is the light velocity in vacuum. Frequency-dependent THz transmission spectra of our designed structures can be obtained from the aforementioned formula.

Figure 2A shows the calculated results of samples M1, M2, and M3 with the same value of γ = 0.3 rad/ps. The transmission dips appear at 1.322, 1.134, and 0.932 THz, respectively. The increase of bar length is accompanied with the decrease of the resonant frequency, broadening of bandwidth, and reduction in amplitude, which are consistent with the simulation results. Figure 2B shows the changed line shape of resonant transmission of M1 with γ while keeping other parameters constant. The frequency position of the resonant dip does not move, but the quality factor decreases with the increased γ from 0.3 to 0.7 rad/ps. The relationships between γ and the full width half maximum (FWHM) of the resonant dips for M1, M2, and M3 are plotted in Figure 2C, showing the FWHM increases with γ and will saturate by further increasing the damping rate. In addition, the longer bar corresponds to a higher value of the FWHM, which is in accordance with the spectral line shape in Figure 2A. Additionally, the curve slope slightly rises with the bar length, suggesting that γ has more obvious influence on the metastructure with a longer bar length.

For the unequal double-bar structures, the impact of the spacing S ₁ on the spectral properties has been investigated, as shown in Figure 3. In the experiment, when S ₁ increases from 2 to 24 μm, the low-frequency resonance dip1 presents a red-shift and its slope becomes sharper, while the high-frequency resonance dip2 shows a blue-shift with a broader bandwidth. The transparent window between dip1 and dip2 gradually emerges and exhibits a red-shift. However, if we further increase S ₁ to 32 μm, the coupling becomes weaker, and the aforementioned changes begin to recover. Figure 3B shows the calculated results of M12 using Eqs 2, 3 based on the bright–bright mode coupling. The fitting parameters under different spacing are listed in Table 1. It is observed that damping rate γ ₁ decreases significantly but γ ₂ increases when S ₁ varies from 2 to 24 μm. Combined with the experimental spectrum, we find that the decreased damping rate corresponds to the increased value of quality factor. The trend of the coupling coefficient κ supports that the coupling strength between two bright modes is essentially enhanced with the increased S ₁. Notably, S ₁ of 24 μm results in a uniform distance of all bars in the sample. After further increasing S _1, γ ₁ increases and γ ₂ decreases, which is consistent with recovery phenomena when S ₁ is 32 μm. Meanwhile, κ decreases in an obvious manner, which is in good agreement with the weakened coupling, as mentioned previously. It is evident that the strongest coupling occurs at S ₁ of 24 μm, and it could be attributed to the interaction between the bars in one unit cell and adjacent unit cells.

To further study the influence of the parameter on the spectral evolution in the coupled metamaterials, we then vary each parameter while keeping other parameters unchanged. We have modulated the parameters γ ₁, γ ₂, and κ, as shown in Figure 4, based on the fitting parameters for S ₁ of 12 μm. The damping rate represents the radiation loss from the metallic bar and affects the quality factor of the resonance and the MIT magnitude. The coupling coefficient indicates the strength of near-field interaction between the substructures in the metastructure. It is found that, for sample M12, the damping rate mainly determines the FWHM of its resonant dip and the MIT window. As shown in Figure 4A, the FWHM of dip1 increases almost linearly, while the value of MIT window shows an exponential decrease with γ ₁ from 0.4 to 0.9 rad/ps. Figure 4B also exhibits the similar influence of γ ₂ on the FWHM of dip2 and the MIT window. It can be seen that κ determines the FWHM of both resonant dips, the MIT window, and the difference in frequency between the two dips (Δf = f _dip2-f _dip1), as shown in Figures 4C,D. The FWHM of dip1 remains lower than that of dip2, which is consistent with the obtained results in Figure 3A. The FWHM of dip1 shows a linearly decreasing trend and that of dip2 increases exponentially, showing that κ has the opposite influence on the FWHM of two dips. It implies that the stronger coupling leads to the higher (lower) value of quality factor of dip1 (dip2) and more obvious asymmetric spectral line shape. The influence of κ on Δf and the value of MIT window is presented in Figure 4D, with both exponentially increasing behaviors which are in good agreement with the red-shift of dip1 and the blue-shift of dip2 in our experimental data.

Similarly, each bar in the M123 structure can individually serve as a bright mode with electric dipole oscillation excited directly by the incident THz wave. By varying the spacing S ₂ between the bars, the modulations of the transmission properties for the 3-bar sample were observed experimentally, as shown in Figure 5A, with three resonant dips and two MIT peaks. When S ₂ increases from 2 to 14 μm, dip1 and dip2 present the red-shift, but the former turns into a sharper slope in line shape. Simultaneously, dip3 exhibits an asymmetric line shape with much broader bandwidth. During this process, the values of MIT vary correspondingly, showing the higher value in MIT2 at small S ₂ and almost the same values in MIT1 and MIT2 at S ₂ of 14 μm. It can be easily seen that the distance of every 2 bars is uniform in the entire sample for S ₂ of 14 μm. However, when S ₂ is further increased to 18 μm, there is no obvious recovery trend, which is totally different from the behaviors in the double-bar structures. Here, dip1 and dip3 have no remarkable difference compared with those for S ₂ of 14 μm, but the significant changes occur with a blue-shift in dip2 and the lower value of MIT2. This could be attributed to the fact that the couplings between the middle bar M2 and M1 (κ ₁₂) or M2 and M3 (κ ₂₃) in one unit cell are decreased with the increased bar spacing, while the coupling between M1 and M3 in the adjacent unit cells might affect the spectra. The analytical fitting to the measured amplitude transmission is shown in Figure 5B, which is consistent with the experimental data. The corresponding fitting parameters are listed in Table 2. It is found that for sample M123, damping rate γ ₁ is nearly invariable with the increasing S ₂ from 2 to 14 μm, accompanied with a little enhancement in γ ₂. However, γ ₃ increases in an obvious manner, suggesting the decreasing quality factor. κ ₁₂ and κ ₂₃ first increase with S ₂ and then decrease after further increasing S ₂ to 18 μm. The coupling weakens if the equidistantly arranged structure is broken. Therefore, the coupling strength can be expected to be modulated via modification of the structural geometry. The simulated results are presented in Figure 5C, showing the consistent spectral evolution behaviors compared with experimental data and theoretical fitting.

In order to further clarify the underlying coupling mechanism of the 3-bar structures, we have simulated surface current and electric field distributions at resonant dip1, dip2, and dip3, respectively. As shown in Figure 6, the surface currents exhibit the classical dipole resonance with linear current distributed in bars and their edges. The electric fields are mainly localized around bar ends. For instance, at the frequency of dip1 with small bar spacing, we can observe that only M2 and M3 are strongly excited by the incident THz wave with the surface currents out of phase and the overlap of electric field distributions between those 2 bars, suggesting that dip1 is the result of coupling between M2 and M3. Similarly, it is evident that dip2 originates the coupling among 3 bars, and dip3 reveals the coupling between M1 and M2. Gradually, the currents and electric fields present the localized distributions at the top and bottom edge of the unit cell with the increased S ₂, suggesting the coupling appears between two adjacent unit cells. Here, the distance to another bar in the adjacent unit cell decreases with S ₂. Moreover, it is interesting that the upper–lower and left–right adjacent couplings can occur in the longest bar M3 by observing the electric field distribution.

Further fitting results are given in Figure 7 to understand the influence of the fitting parameters on the theoretical curve of transmission. We modulate the parameters γ _i (i = 1,2,3), κ ₁₂, and κ ₂₃ based on the fitting parameters for S ₂ of 14 μm. Figure 7A shows the FWHM of dip1 increases, and the value of the MIT1 window shows a distinct decrease when γ ₁ changes from 0.4 to 0.7 rad/ps. However, γ ₂ determines the FWHM of dip2, MIT1, and MIT2 windows simultaneously, as shown in Figure 7B. The FWHM of dip2 increases accompanied with an obvious decrease of the value of the MIT1 window and a slight increase of the value of the MIT2 window. It is obvious that γ ₂ has the greater impact on the value of the MIT1 window than that on the value of the MIT2 window in our studied range. Figure 7C indicates that γ ₃ mainly affects the FWHM of dip3 and the value of the MIT2 window. The coupling coefficients also affect the FWHM. Because it represents the interaction between two oscillators, the coupling coefficient will simultaneously alter the frequency positions of two resonant dips interacting with each other and the MIT peak located between those two dips. This can have an effect on the symmetry of spectral line shape and make the FWHM increase or decrease with increasing the coupling at different resonant dips. The influence of coupling coefficient on the THz spectral properties is depicted in Figures 7D,E. It is shown that the value of the MIT1 window is enhanced evidently with the increasing κ ₁₂, accompanied with the attenuation of the FWHM of dip1 and the slight increase of the FWHM of dip2. Similarly, the value of the MIT2 window increases with κ ₂₃, but the FWHMs for both dip2 and dip3 decrease correspondingly. The FWHM of dip3 with a higher value than that of dip2 indicates the lower quality factor, which is in good agreement with the experimental and simulated results. Furthermore, it is noticeable that coupling coefficient also affects the shift of resonant frequency, as shown in Figure 7F, where Δf ₁ = f _dip2-f _dip1 and Δf ₂ = f _dip3-f _dip2. Although both Δf ₁ and Δf ₂ are enhanced with κ ₁₂ and κ ₂₃, respectively, κ ₁₂ has a greater impact on resonant frequency shift than κ ₂₃. It is found that a larger frequency shift presents much stronger coupling. In addition to the aforementioned discussion about the resonant dips and MITs, the influence of parameters on the whole spectral lines is presented in Figure 8. The strips of three dips and two MITs are illustrated clearly. It should be noticed that γ ₁ also affects the MIT2 window slightly, κ ₁₂ has an influence on the MIT2 window, and κ ₂₃ has an impact on the MIT1 window, as shown in the areas indicated with the dashed lines in Figures 8(A,D,E).

It is evident that the bar spacing involved with the structure coupling has an impact on the spectral evolution of THz metamaterials with unequal-length bar structures to modulate the resonant dips, the transparency peaks, and the quality factors of the spectra. Hence, those metastructures with different geometrical parameters will possess different spectral sensing properties. We have found that not only the resonant dips but the MIT peaks can also be utilized in analyte sensing. Figure 9 exhibits the simulated frequency shift with the refractive index n of the overlayer on the metamaterial. For metastructure M12, when spacing S ₁ is 2 μm, Δf at MIT is slightly greater than those at the two other dips. With S ₁ of 24 μm, Δf increases significantly for MIT and dip1, but the sensing ability of dip2 decreases because the increased coupling results in the Fano-like asymmetric line shape with sharp dip1 and broadening dip2 accompanied with an obvious MIT peak. For metastructure M123 with S ₂ of 2 μm, three dips and two MIT peaks exhibit the similar frequency shift as the function of n with slightly higher values at MITs. However, at S ₂ of 14 μm, Δf does not increase and its value of dip3 decreases in an obvious manner, indicating the decreased quality factor and poor sensing performance. It is found that for M12 and M123 samples, although the coupling coefficients increase with spacing, the damping rates, which have a great impact on the FWHM, are different with a most remarkable reduction for dip1 at S ₁ of 24 μm compared with the values of other dips, suggesting the good sensing ability of dip1 and its accompanied MIT peak in M12. It can be expected that such metastructures can be selected for sensing experiments in the future. Our obtained results indicate that we could explore and optimize the sensing parameters by investigating the evolution of coupled resonances and their sensing properties with the introduction of dielectric overlayers.

In conclusion, we systematically investigated and analyzed the Lorentz oscillator model from a single resonator to three resonators combined with the experiment and simulation to explore the underlying coupling mechanism, especially the influence of the spacing between bars on the spectral evolution in the double-bar and 3-bar structures. The experimental data exhibit that the resonant dips, the MIT peaks, and the quality factors of the spectra vary with the bar spacing, indicating that not only intracoupling but also intercoupling together play important roles in electromagnetic responses when the spacing between bars is increased. It is noticeable that the spectral evolution with bar spacing in 3-bar structures is different from that in double-bar structures, due to the reason that the middle bar coupling effect is decreased with the other 2 bars’ increased coupling in the M123 sample. Furthermore, the dependences of transmission spectra on the fitting parameters are calculated based on the coupled model. It is found that the damping rate determines the FWHM of its corresponding resonant dip and the MIT peak, which is related to the amplitude of transmission spectra. The coupling coefficient mainly determines the frequency shift and changes the FWHM of its corresponding dip and the MIT window slightly. Further simulated results show that the sensing parameters can be obtained and optimized by investigating the spectral variation of those metastructures after adding the dielectric overlayer. Our research could reveal the mechanism of resonance coupling, clear the role of fitting parameters, help us to design and mimic more complicated coupling mode systems, and establish the foundation to further explore the coupling in THz modulators and sensing devices.