As a distinguished interdisciplinary frontier field in contemporary research, flexible electronic technology spans electronics and photonics. Characterized by their ultra-thin, flexible, and lightweight design, flexible electronic devices are gradually replacing traditional rigid devices and present substantial potential in applications such as electronic skin, real-time health monitoring, and industrial surface monitoring [1–6]. Recently, the rapid development of terahertz (THz) technology has expanded the operating frequency range of flexible electronic devices, bridging the gap between microwave and infrared light [7]. THz waves are highly valued in various technological fields for their ability to penetrate non-conductive materials and provide high-resolution imaging, particularly in wireless communications, medical imaging, security screening, and radar detection [8–12]. However, integrating THz technology into flexible devices requires innovative materials and structures to maintain electromagnetic stability under mechanical stress.

Metamaterials (MMs) are composite structures composed of periodically arranged subwavelength units, designed to effectively manipulate electromagnetic waves for various applications, such as negative refraction, super lenses, polarizers, holographic imaging, and perfect absorption [13–20]. As an engineered artificially material, MMs can precisely control the phase and amplitude of THz waves by adjusting the structure parameters. This ability to customize electromagnetic responsiveness positions MMs as a crucial technology for a wide array of THz applications. Traditionally, MMs relied on rigid silicon-based substrates [21–23]; however, their inherent brittleness and rigidity are not ideal for applications requiring flexibility. In response to these limitations, researchers are increasingly exploring flexible polymeric materials like PDMS, PI, and PET as viable alternatives to traditional silicon substrates [24–27]. These polymers not only offer superior mechanical flexibility and environmental durability but also high transparency in the THz frequency range. Recent studies have shown that the spacing between resonant elements can be adjusted through mechanical deformation of the substrate, allowing for tuning capabilities [28–31]. MMs utilizing flexible substrates are primarily used in applications that involve measuring material deformation or in tunable detection tasks [32–35]. However, there is a lack of structural designs that maintain stable resonant peaks in terms of frequency and amplitude under mechanical strain, which are essential for passive non-destructive evaluation or sensing in curved environments. Furthermore, flexible substrates also exhibit dielectric losses that attenuate electromagnetic wave transmission through the medium, thereby reducing the overall sensitivity and resolution.

Fortunately, metal-aperture MMs (MA-MMs) offer an ingenious and effective solution to the aforementioned limitations. Utilizing the “cold processing” capabilities of femtosecond laser technology with short pulses and high energy, MA-MMs are manufactured by ablating subwavelength apertures on micrometer-scale metal films, efficiently avoiding the dielectric loss of silicon substrates [36]. In addition, MA-MMs preserve the inherent stability, flexibility, and metallic properties of the thin films, enabling them to be bent without damaging the resonant structures. Despite being well suitable for the development of THz flexible devices, MA-MMs are limited by the inherent Ohmic losses, which leads to resonance peaks with low quality factor (Q-factor). Consequently, the strong coupling effects among multiple surface plasmon modes have been proposed to compensate for the performance defect of metallic MMs [37,38]. Unlike single resonance modes, the hybridization behavior between resonant modes provides a pathway for energy exchange, enhancing the interaction with light and matter. By designing multiple resonant structures within the unit cell to induce multimode coupling, it is possible to generate tunable Rabi oscillations with high quality factors [39–41].

In this paper, an MA-MM with two analogous resonant structures is proposed, which can stimulate a strong coupling effect to generate a Rabi splitting. To replicate actual bending conditions, we innovatively conducted full-wave simulations in the time domain module, observing the variations of frequency and amplitude for each resonant peak at various bending angles. The simulated results were verified by experimental measurements, which were used to calculate the Q-factor of the Rabi splitting and estimate the stability of the response performance. Meanwhile, the transmission responses under various bending configurations were also investigated through both simulations and experiments. The results aim to verify the suitability of flexible MMs for diverse dynamically curved environments.

Considering that the unit-cell module fails to simulate the bending deformation of the whole structure, the time domain solver was employed to calculate the transmission response of the proposed MMs. We first constructed a flat MM array and subsequently employed the “bend tool” function within the “Model” module to introduce curvature to the array. In order to keep the simulation results only affected by bending deformation, the position of the emitting port was set 100 μm above the highest point of the curved MMs, and the THz wave was always perpendicular to the plane where the highest point was located. The receiving port was set 100 μm below the lowest point of the curved MMs. Full-wave simulations were conducted using two different port configurations [Fig. 1(d)]. The structure was composed of a 10×10 periodic unit array (3000μm×3000μm), where the first port covers an area equivalent to the full structure, and the second port spans an area of 500μm×500μm. Simulation comparisons from Fig. 1(b) indicate that the result for the full port aligns with that of the unit-cell module. However, the result of using the localized port exhibits a reduction in amplitude due to a decrease in the area available to excite surface plasmons, as evidenced by the electric field distribution of the Rabi splitting in Fig. 1(e). The red outline in the figure indicates the boundaries of the port. Under the full port configuration, the electric field nearly spans the entire structure. As for the localized port configuration, the excited electric field exhibits higher energy but is primarily concentrated within the port and its immediate surroundings. In actual detection, THz waves typically shine on the surface in a focused form, which do not cover the entire MMs. To ensure that future simulations accurately reflect real measurement conditions, all simulations will be conducted using the local port approach.

Figure 1(c) shows physical and optical micrographs of MA-MMs, respectively. MA-MMs are fabricated by using femtosecond laser direct writing technique to ablate periodic structures on thin metal films. Femtosecond laser ablation harnesses the high energy density and ultrashort pulse duration to rapidly deliver energy to the surface of the material, inducing localized melting or evaporation [45]. This technology facilitates the formation of high-precision micrometer-scale machining areas while capably reducing thermal effects. Utilizing position control via a displacement platform, femtosecond lasers can ablate metal films to produce surface structures with intricate patterns. Compared to traditional lithography technology, this one-molding processing method has fewer steps, lower cost, and shorter processing time. Independent MA-MMs not only effectively avoid the dielectric loss caused by the substrate effect, but also retain the original physical and chemical stability and flexibility of metal films, which is especially suitable for the application of flexible THz devices. As shown in the red solid line in Fig. 1(b), the experimentally measured transmission peaks and Rabi splitting are located at 0.6 THz, 0.66 THz, and 0.76 THz, respectively, which are consistent with the simulation results.

The transmission spectra were calculated by the time-domain solver as the proposed MM was bent from 0° to 180° [Fig. 2(a)]. Comparing the simulation results of each bending angle in Fig. 2(b), there is no significant frequency shift for the resonance peaks and the Rabi splitting, indicating that bending deformation does not change the excitation mode of surface plasmons. Figures 2(d) and 2(e) depict the frequency fitting curves and frequency shift rates of three peaks for thoroughly analyzing the variations in the transmission spectra. The frequency of the second transmission peak exhibits a slight fluctuation within 0.05 THz, while the frequencies of the first peak and the Rabi splitting keep virtually unchanged. Moreover, the transmission amplitude of each peak was selected for linear fitting, as shown in Fig. 2(c). It is evident that the transmission amplitude of resonance peaks decreases significantly as the bending angle increases, dropping from 0.86 and 0.91 to 0.53 and 0.59, respectively.

For further observing the influence of bending deformation on transmission amplitude, the electric field intensity distribution of the cross section under different angles is monitored. As can be seen from Fig. 2(f), bending deformation enlarges the excitation range of surface plasmons, while simultaneously reducing the electric field intensity within the port. At a bending angle of 180°, THz waves traversing the MMs generate extensive diffraction during propagation, which exceeds the reception capabilities of the output port. THz waves that exceed the port range fail to be captured, resulting in diminished transmission efficiency. The electric field distribution on the incident plane also shows that the excitation range enlarges with the increase of the bending angle (Fig. 8 in Appendix B). These observations confirm that the transmission amplitude is sensitive to the changes in bending deformation, as is evidenced by the excitation range of surface plasmons. Unlike the resonance peaks, the transmittance of the Rabi splitting exhibits relative consistency, staying below 0.03. It suggests that bending deformation does not impede the strong coupling effect between the two LSPP modes. The amplitude of the Rabi splitting is determined by the coupling strength rather than the excitation intensity of any single mode. Consequently, MA-MMs with multiple resonance modes can preserve the original response of the Rabi splitting under any bending conditions, which is crucial for the practical application of THz flexible devices.

To validate the characteristics of the designed MA-MMs, a home-made flexible stretching fixture was utilized for experimental measurements. The sample was securely clamped at both ends of the fixture using bolts. One end of the fixture was held in place while the opposite end was adjusted by a displacement platform to induce stretching or compression for bending. As shown in Fig. 3, the platform was displayed by a displacement (D) ranging from 0 to 3 mm, leading to the MA-MMs displaying four different bending states. Transmission spectral responses of the flexible MMs with a bandwidth of 0.4–1.0 THz were measured by the THz time-domain spectroscopy system, as presented in Fig. 3(b). The transformation law of transmittance observed in four bending states is identical to the simulation results. The strong coupling effect consistently induces a pronounced splitting peak even at the largest bending angles. As shown in Fig. 3(c), the amplitude of the transmission peaks decreases from 0.4 and 0.44 to 0.06 as the bending angle increases, while the amplitude of the Rabi splitting always remains below 0.05. Meanwhile, the frequency shift and offset rates of three characteristic peaks are also depicted in Figs. 3(d) and 3(e). The frequency of the Rabi splitting changes from an initial 0.66 THz to 0.63 THz, reaching its maximum offset of 4.2% at S=3mm. The second resonance peak exhibits the most noticeable frequency shift, changing from 0.75 THz to 0.81 THz, with a maximum offset rate of 8%. The differences in frequency shift arise from the interference effect caused by bending MMs. In practical measurements, interference effects often cause fluctuations in the transmission amplitude curves, leading to inaccuracies in the peak frequencies due to the limited resolution of the system. Additionally, to evaluate whether the MM maintains stable response performance under bending conditions, the Q-factors for four displacement states were calculated in Fig. 3(f). The Q-factor was 12.28 without bending but slightly declined to 10.33 at S=2mm, where it reached its lowest point. When the MM is further bent to S=3mm, the Q-factor increases to 10.46. This indicates that bending has minimal impact on the Q-factor, which is crucial for applying the designed MMs in flexible devices.

In order to further verify the stable response performance for the designed MMs, the electromagnetic responses under four bending states were analyzed. Figure 4(b) shows optical microscopic images of MA-MMs with different L1 and L2 after femtosecond laser ablation. Similarly, we moved the displacement platform from 0 to 3 mm causing the MA-MM to bend to varying degrees. Measured transmission spectra are shown in Figs. 4(c)–4(f), which show that the transmittance of the four samples all evidently decrease with the increase of the bending angle. In this case, the Rabi splitting induced by the strong coupling effect always exists even when the amplitude of the resonant peak drops below 0.1. As depicted in Fig. 4(h), the frequencies for each sample remain stable across four different bending angles, with a maximum coefficient variation of 4.2% (0.661 THz to 0.633 THz) observed during deformation. Specifically, the coupling strength between the two resonant modes increases as L1 decreases, leading to a blue shift in Rabi splitting and a significant rise in Q-factor. The Q-factors for each Rabi splitting were also calculated and depicted in Fig. 4(g) to further explore the variation in Q-factor following bending deformation. Notably, when L1 is reduced to match the length of L2, a higher coupling strength is observed, resulting in a larger Q-factor of 67.45. Analysis of the samples post-bending deformation reveals a relatively stable Q-factor. Additionally, investigations were conducted on samples with different L2 but the same L1 showing a consistent trend in Q-factor variation across different bending angles, highlighting the stable electromagnetic response of MA-MMs with varying coupling strengths. Unlike samples with different L1, the frequency of Rabi splitting decreases as L2 increases, while the Q-factor continues to rise. These phenomena regarding the influence of structural dimensions are crucial for the design of MA-MMS.

Due to differences in geometric dimensions, the bending extent cannot be simply determined by the bending angle. Therefore, the curvature is used to represent the deformation extent under the same bending angle, which can be calculated by the bending radius. The relationship between curvature (C) and bending radius (R) can be described as C=1/R. As shown in Fig. 5(a), the unit number was increased from N=5 to N=15, resulting in an increase in the bending radius from R=0.75mm to R=2.25mm. The localized port setting is 500μm×500μm, the bending angle is R=180°, and the periodicity is P=300μm. Figure 5(b) shows the inverse relationship between the curvature and the number of units. The curvature decreases as the number of units increases, indicating a reduction in bending deformation. Figure 5(c) depicts the simulated transmission spectra for different array sizes, with dashed lines and highlights marking the frequency ranges of three characteristic peaks. The frequencies of the two transmission peaks exhibit irregular shifts with accompanying amplitude fluctuations, particularly at higher curvatures. This phenomenon arises from bending deformation, which disrupts the uniformity of the resonance structure array, leading to interference effects during the propagation of THz waves. The superposition or cancellation of these waves enhances or attenuates the original transmission amplitude, thereby causing frequency shifts in the resonance peaks. The frequency fitting curve in Fig. 5(d) indicates frequency fluctuation in the second transmission peak at C=0.95 (N=11), while a larger offset is observed for the first peak at C=1.49 (N=7). Besides, it is evident that the frequency of the splitting peak shows little variation. Owing to the strong coupling effect between the two surface plasmon modes, the induced low-transmittance trough is resistant to frequency shifts caused by the interference effects.

The transmission amplitude fitting curves for the three peaks are depicted in Fig. 5(e), demonstrating that the transmittance of the Rabi splitting consistently maintains a stable level below 0.05. The transmittance of the two transmission peaks increases continuously with curvature; as curvature changes from 0.7 to 2.1, the amplitudes rise by 0.24 and 0.44, respectively. To thoroughly understand the impact of array size on amplitude, the electric field distribution across the cross-section is also plotted, as shown in Fig. 5(f). The excitation range of surface plasmons expands as the number of units increases, while the electric field enhancement concentrates primarily below the input port. When N=5, the port range covers the entire MMs, allowing THz waves passing through the aperture to be fully received by the output port. However, as N increases, the proportion of the port in the array area steadily decreases. THz waves excited outside the port range fail to be captured, resulting in a decline in transmission efficiency.

As illustrated in Fig. 6(a), the fabricated MA-MMs are laid over the bending models produced by 3D printing technology. The processed array will cover the hollow part of the bending model to ensure that the bending radius is the sole influencing factor in the experiment. Figure 6(b) displays four bending models with radii decreasing from 25 mm to 10 mm, and the calculated curvatures are 0.04mm−1, 0.05mm−1, 0.07mm−1, and 0.1mm−1, respectively. With the help of the bending models, the spectral responses at different curvatures are measured, as shown in Fig. 6(c). As predicted by the simulation results, the frequencies of the three characteristic peaks stay extremely steady, and the amplitude of the Rabi splitting maintains a constant value across the four bending models. Figure 6(d) shows that the amplitude increases only slightly from 0.12 to 0.14, while the amplitudes of the two transmission peaks increase with curvature. The amplitude of the first transmission peak rises from 0.56 to 0.7, while that of the second transmission peak increases from 0.64 to 0.77. Figure 6(e) illustrates the frequency shifts of the three peaks, which hold steady at 0.58 THz, 0.66 THz, and 0.76 THz, respectively. The insensitivity of the characteristic frequencies to the array size demonstrates that the designed MA-MMs may function as a flexible fixture for detecting various conditions, such as on the curved surfaces of large machinery or small electronic instruments. Additionally, Fig. 6(f) plots the Q-factor of the Rabi splitting across the four curves, which ranges from 8.5 and 8.73. This minor variation further confirms that the detection performance of the flexible MMs remains stable under bending deformation.

In conclusion, an MA-MM was designed by harnessing the strong coupling effects between two surface plasmon modes to ensure stable response performance under bending deformation. We employed a time-domain solver to compute the transmission spectra of the entire MM, innovatively utilizing localized ports to accurately simulate the changes of the electromagnetic response at different bending deformations. With the increase of the bending angle, the transmission amplitude decreases while the resonance frequency stays constant. The discrepancy is attributed to the bending deformation modifying the excitation range of surface plasmons without altering the resonance mode. Importantly, the experiment results confirm that the response properties of the Rabi splitting generated by the strong coupling effect can remain relatively stable in bending deformation. Despite an 87.6% reduction in the transmittance of the resonant peak, the Q-factor of the Rabi splitting exhibited only a 14.8% decrease. Additionally, the transmission spectra of MA-MMs with varying array sizes were measured through both simulation and experimental methods to evaluate the response performance of each Rabi splitting peak across diverse complex environments. The findings indicate that the frequency and Q-factor of Rabi splitting are also unaffected by the array size, with the maximum deviation in Q-factor being merely 2.63%. It is conceivable that the designed MA-MMs offer a viable solution for flexible THz devices requiring consistent performance across different bending states.