Quadruple Plasmon-Induced Transparency and Dynamic Tuning Based on Bilayer Graphene Terahertz Metamaterial

Abstract

This study proposes a terahertz metamaterial structure composed of a silicon–graphene–silicon sandwich, aiming to achieve quadruple plasmon-induced transparency (PIT). This phenomenon arises from the interaction coupling of bright–dark modes within the structure. The results obtained from the coupled mode theory (CMT) calculations align with the simulations ones using the finite difference time domain (FDTD) method. Based on the electric field distributions at the resonant frequencies of the five bright modes, it is found that the energy localizations of the original five bright modes undergo diffusion and transfer under the influence of the dark mode. Additionally, the impact of the Fermi level of graphene on the transmission spectrum is discussed. The results reveal that the modulation depths (MDs) of 94.0%, 92.48%, 93.54%, 96.54%, 97.51%, 92.86%, 94.82%, and 88.20%, with corresponding insertion losses (ILs) of 0.52 dB, 0.98 dB, 1.37 dB, 0.70 dB, 0.43 dB, 0.63 dB, 0.16 dB, and 0.17 dB at the specific frequencies, are obtained, achieving multiple switching effects. This model holds significant potential for applications in versatile modulators and optical switches in the terahertz range.

1. Introduction

Surface Plasmon Polaritons (SPPs) [1], which are collective resonances produced by the interaction of incident photons and surface electrons [2,3], have been extensively studied in metallic and graphene materials [4,5,6,7,8]. Graphene, a tunable two-dimensional hexagonal lattice material, exhibits similar dielectric constants to noble metals in the terahertz and infrared wavelengths [9]. Recently, researchers have noted the occurrence of the PIT phenomenon in both graphene [10,11,12] and metallic structures [13]. However, those waveguide structures utilizing noble metals only allow for static adjustments and suffer from significant losses [14]. In comparison, graphene offers significant advantages with dynamic modulation through gate voltage tuning, making it highly versatile [15,16]. Because of its low loss [17,18,19], excellent dynamic tunability, and significant local field enhancement effect [18,19,20,21], graphene finds crucial applications in optical sensing [22,23,24], optical switching [25,26], optical storage [27,28], etc.

Graphene-excited surface plasmons show robust localized field enhancement effects in the terahertz and infrared wavelengths. The plasmon excitations on the surface of graphene can generate plasmon-induced transparency (PIT) [29], which is a phenomenon similar to the electromagnetically induced transparency (EIT) [30,31] caused by atomic coherence. However, due to the challenging experimental conditions and high costs associated with EIT, it is not conducive to practical applications [32]. PIT is a transparent phenomenon produced in classical systems by destructive interference. It shares similarities with EIT, but the required experimental conditions are not as stringent as those of EIT [33,34]. Recently, PIT has been gaining significant attention among numerous researchers, owing to its unique characteristics [35,36,37]. Studies have indicated that the interplay of bright and dark modes can generate the PIT effect [38]. Bright modes, which are directly excited by the incident light, induce resonance at specific frequencies, thereby causing a continuous transmission dip within the transmission spectra. Conversely, dark modes, being unable to be directly stimulated by incident light, typically do not impact the spectra. However, in those designed bright–dark mode symbiosis structures, the dark mode is brought into action through an optical field created by the interaction of the incident light and the bright mode, leading to anomalous light transmission and a PIT window.

In this study, we present a dual-layer graphene waveguide structure composed of five graphene ribbons (M1, M2, M3, M4, M5) and double vertical graphene ribbons (DVRG) to achieve multiple PIT windows. The proposed structure is a multilayer terahertz metamaterial, in which a silicon substrate is sandwiched between two graphene layers, each connected to its respective electrodes, respectively. Furthermore, the destructive interference between the structure’s four bright modes and one dark mode results in significant PIT. The theoretical results of quadruple PIT are fitted using the coupled mode theory (CMT), which shows high consistency with the outcomes of the finite difference time domain (FDTD) simulations. Importantly, excellent optical modulation characteristics are achieved for our structure, since the maximum and minimum modulation depths (MDs) reach 97.51% and 88.20%, corresponding to the insertion losses (ILs) of 0.43 dB and 0.17 dB, respectively. Therefore, our proposed structure offers significant prospects in areas such as optical switching and modulation.

2. Simulation Model

The three-dimensional view depicted in Figure 1a showcases a terahertz metamaterial structure consisting of a double-layer graphene. The substrate of the structure is made of silicon, serving as the supporting layer, with a thickness of 0.45μm. The upper layer consists of a single rectangular graphene (M1), while the lower layer is composed of four rectangular graphene sheets (M2, M3, M4, M5) and double vertical rectangular graphene (DVRG) sheets, with each graphene layer having a thickness of 1 nm. Figure 1b provides a distinct depiction of the top view of the proposed metamaterial structure. Figure 1c,d provide the detailed parameters of the upper and lower layers of graphene. The geometric parameters are as follows: ${\mathrm{d}}_1=0.9 \mu \mathrm{m}$, ${\mathrm{d}}_2=0.3 \mu \mathrm{m}$, ${\mathrm{d}}_3={\mathrm{d}}_4=0.8 \mu \mathrm{m}$, ${\mathrm{d}}_5={\mathrm{d}}_6=1.1 \mu \mathrm{m}$, ${\mathrm{d}}_7={\mathrm{d}}_8={\mathrm{d}}_9={\mathrm{d}}_{10}=0.2 \mu \mathrm{m}$, ${\mathrm{d}}_{11}={\mathrm{d}}_{12}=0.6 \mu \mathrm{m}$, ${\mathrm{d}}_{13}=0.5 \mu \mathrm{m}$, ${\mathrm{d}}_{14}=0.7 \mu \mathrm{m}$, ${\mathrm{d}}_{15}={\mathrm{d}}_{16}=0.2 \mu \mathrm{m}$, ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}=2 \mu \mathrm{m}$, ${\mathrm{h}}_1=0.2 \mu \mathrm{m}$, ${\mathrm{h}}_2=0.1 \mu \mathrm{m}$, ${\mathrm{h}}_3=0.15 \mu \mathrm{m}$. The incident plane wave, serving as the excitation source, propagates in the negative direction of the z-axis, perpendicular to the x–y plane. Monitors are positioned at a distance of 0.15 μm above the bottom of the silicon substrate. In the X and Y directions, periodic boundary conditions (PBC) are used and in the Z direction, a perfect matching layer (PML) is applied. PBC is suitable for periodic structures, reducing computational resource requirements and avoiding edge effects. Meanwhile, PML effectively absorbs reflected waves, ensuring the accuracy of simulation results. The mesh step sizes in the x and y directions are set to 0.1 µm, while the mesh step size in the z direction is set to 0.02 µm, with a utilization of 22.5 mesh nodes vertically. The simulation time is set to 15,000 fs. The Fermi level of the graphene ribbon is 0.6 eV.

Within the context of our proposed metamaterial structure, Molecular Beam Epitaxy (MBE) presents a promising approach for fabricating the upper Si regions. MBE is capable of growing crystalline thin films with excellent layer thickness and material quality control. This technique involves the precise deposition of atoms or molecules onto a substrate surface under ultra-high vacuum conditions. By finely adjusting growth parameters such as temperature, deposition rate, and silicon precursor flux, the controlled growth of Si layers with nanoscale thickness can be achieved.

The Kubo model stands as a well-established theoretical framework in the realm of condensed matter physics. It is harnessed to ascertain the electrical conductivity of graphene, which can be formulated as follows [38]:

$$\sigma \left(\omega \right)={\sigma }_{intra}+{\sigma }_{inter}$$

(1)

$${\sigma }_{intra}=\frac{2\mathrm{i}{e}^2{k}_{B}T}{\pi {\hslash }^2\left(\omega +i{r}^{-1}\right)}\ln \left[2\cosh (\frac{E_F}{2k_{B}T})\right]$$

(2)

$${\sigma }_{inter}=\frac{e^2}{4{\hslash }^2}\left[\frac{1}{2}+\frac{1}{\pi }\arctan \left(\frac{\hslash \omega -2E_F}{2k_{B}T}\right)\right]-\frac{e^2}{4{\hslash }^2}\left[\frac{i}{2\pi }\ln \frac{{\left(\hslash \omega +2E_F\right)}^2}{{\left(\hslash \omega +2E_F\right)}^2+4{\left(k_{B}T\right)}^2}\right]$$

(3)

where $e, \hslash , k_B, \omega , E_F$ stand for the electron charge, the reduced Planck constant, the Boltzmann constant, the angular frequency of incident light and the Fermi level of graphene, respectively. $\mu$ is carrier mobility, and $v_F={10}^6 \mathrm{m}{\mathrm{s}}^{-1}$ denotes the Fermi velocity. From $\tau =\mu E_F/e{v}_F^2$ [39], the carrier relaxation time $\tau$ can be deduced. The ambient temperature is configured as $T=300 \mathrm{K}$. By adjusting the bias voltage, the Fermi level of graphene $E_F$ can be modulated from 0.6 eV to 0.72 eV. Under the conditions established in this study, ${\sigma }_{inter}$ can be disregarded, so $\sigma$ can be obtained [40,41]:

$$\sigma \left(\omega \right)=\frac{e^2{E}_F}{\pi {\hslash }^2}\frac{i}{\omega +i{\tau }^{-1}}$$

(4)

In the proposed graphene metamaterial structure, with silicon as the substrate, the propagation constant $\beta$ and effective refractive index $n_{eff}$ are as follows [38,42]:

$$\frac{{\epsilon }_{si}}{\sqrt{{\beta }^2-{\epsilon }_{si}{k}_0{}^2}}+\frac{{\epsilon }_{air}}{\sqrt{{\beta }^2-{\epsilon }_{air}{k}_0{}^2}}=-\frac{i\sigma }{\omega {\epsilon }_0}$$

(5)

$$n_{eff}=\frac{\beta }{k_0}$$

(6)

where ${\epsilon }_{si}, {\epsilon }_{air}, \mathrm{and} {\epsilon }_0$ present the relative dielectric constants of silicon, air and vacuum dielectric, respectively. $k_0$ and ${\eta }_0$ correspond to the wave number in free space and the inherent impedance, respectively.

3. Results of Simulation and Theoretical Analysis

The FDTD simulation process involves numerically solving Maxwell’s equations using a grid-based approach. It discretizes both time and space, allowing us to track the evolution of electromagnetic fields over discrete time steps, enabling analysis of its optical behavior and properties. In the simulation, the graphene layer is discretized into grid cells both in time and space. The electromagnetic field interacts with the graphene layer based on its material properties, and the FDTD algorithm iteratively solves Maxwell’s equations in a timely manner to calculate the propagation of the electromagnetic field and its interaction with the graphene layer. To elucidate the physical mechanism behind the emergence of the quadruple PIT phenomenon, the proposed structure is simulated using the FDTD method, and the transmission spectra are depicted in Figure 2. When an x-polarized beam is vertically incident on the structure, different graphene structures exhibit distinct Lorentz curves. In Figure 2a, the M1, M2, M3, M4, and M5 structures generate five Lorentz modes, considered as the bright modes, which can be readily stimulated by the incident plane wave. The resonant frequencies of these bright modes, denoted as $f_0$, $f_1$, $f_2$, $f_3$, and $f_4$, are 2.75 THz, 3.37 THz, 3.98 THz, 4.63 THz, and 5.44 THz, respectively. When the incident light is launched on the DVRG, similar bright modes are hardly excited, but the DVRG generates a mode with a transmission of approximately 1, which serves as the dark mode. The dark mode is incapable of direct excitation by the incident plane wave, yet it engages with the localized field generated by the bright modes. Therefore, the destructive interferences between the five bright modes and the dark mode lead to a pronounced quadruple PIT phenomenon, as shown by the solid red curve, indicating the formation of the quadruple PIT windows. The four PIT peaks are labeled as $f_a$, $f_b$, $f_c$, and $f_d$ from left to right in Figure 2b with the red curve. At the resonant frequencies of the bright modes ($f_0$, $f_1$, $f_2$, $f_3$, and $f_4$), the originally bright modes transform into the states with transmission rates close to 1, as confirmed in Figure 2a,b.

To further delve into the underlying formation mechanism that governs the creation of the quadruple PIT, the electric field distributions of the bright modes, the dark mode, and the overall structure at the resonant frequencies of the four PIT windows are analyzed, respectively. The SPPs’ excitation of graphene determines the intensity of the electric field, where a higher excitation level results in stronger electric field energy. Figure 3a illustrates the transmission spectra for the cases of only M2, DVRG, and the entire structure. Figure 3b,c represent the electric fields of the bright mode M2 and the dark mode DVRG at the resonant frequency of 3.37 THz, respectively. In Figure 3b, the electric field distribution at the edges of graphene appears in a darker shade of red, indicating higher energy levels. Graphene strongly couples with the incident plane wave for generating the bright mode, and then the energy of the incident wave tightly binds with the graphene. In Figure 3c, when only two vertical graphene strips DVRG are presented, the electric field energy primarily localizes outside the graphene strips, with a very weak energy level. Within the terahertz range, it exhibits no significant response to the incident plane wave, behaving as a dark mode. Figure 3d represents the spatial distribution of the electric field within the entire structure at the first resonant peak $f_a$. The energy of the bright mode M2 experiences significant attenuation and transfers to the dark mode DVRG, induced by the near-field coupling between the bright and dark modes. At the resonance frequency of 3.27 THz, a distinctive and clearly pronounced transparency window emerges as a result of the destructive interference occurring between the two distinct modes. Figure 4a depicts the transmission spectra for the scenarios involving only M3, DVRG, and the entire structure. Figure 4b,c depict the spatial distribution of the electric field within the bright mode M3 and the dark mode DVRG at the resonant frequency of 3.98 THz. It can be observed that when M3 and DVRG exist independently, energy localization occurs independently under the influence of the incident light. In Figure 4b, the energy is mainly contributed by the rectangular graphene strip. In Figure 4c, the weaker energy is localized along the edges of the two vertical graphene strips. Referring to Figure 4d, which depicts the energy distribution of the electric field across the entire structure at the second resonance peak, the enhanced electric field is mostly distributed on the M1, M2 and DVRG, while the energy within the M3 region diminishes. This demonstrates the transfer of electric field energy from the M3 structure to the DVRG and M2, induced by the resonant interaction between the bright and dark modes. Figure 5 depicts the emergence of the second notable PIT transparency window, resulting from the destructive interference between the two modes, which is similar to the previous analysis. Figure 5a illustrates the transmission spectra for the cases of only M4, DVRG, and the entire structure. Figure 5b,c represent the electric fields of the bright mode M4 and the dark mode DVRG at the resonant frequency of 4.63 THz. At the third PIT transparency window, the dark mode is indirectly excited, and the electric field energy primarily propagates from the M4 structure to the DVRG, as depicted in Figure 5d. Figure 6a illustrates the transmission spectra for the cases of only M5, DVRG, and the entire structure. Figure 6b illustrates the direct excitation of M5 by the incident light, resulting in high electric field energy. At the fourth resonance peak $f_{\mathrm{d}}$, the electric field on the M5 structure weakens, and the energy transfers to the M1 and M2 structures and weakly couples with the vertical graphene strips. Furthermore, the bright mode M1 at the resonant frequency of 2.75 THz undergoes destructive interference with other graphene structures, causing the transmittance changed from ~0 to ~1. Figure 7a illustrates the transmission spectra for the cases of only M1, DVRG, and the entire structure. Figure 7b indicates that the energy at 2.75 THz is primarily contributed by M1. Figure 7c demonstrates that the edges of the vertical graphene strips DVRG exhibit a conspicuous presence of extremely weak electric field energy, and a low energy level across the entire structure is revealed in Figure 7d. The energy on M1 is significantly attenuated, and the energy undergoes transfer and diffusion, primarily weakly localizing outside the vertical graphene strips DVRG and graphene structures. Therefore, under the influence of the incident plane waves, the bright mode M1 at the resonant frequency of 2.75 THz experiences destructive interference with the remaining graphene structures. This is investigated in the transmission spectrum with quadruple PIT, where the absorption peak disappears and the transmittance increases from 0.10 to 0.76. In Figure 8, the electric field energy is concentrated at the central local position of M1, and the energy level is relatively low. This phenomenon indicates that at the frequency of approximately 5.64 THz, the geometrical features of structure M1 resonate with the wavelength of the incident electromagnetic wave, leading to the concentration of electric field energy in specific regions. Due to the lower energy level, this resonance behavior is not as pronounced as the main peak. Therefore, a smaller peak appears in the transmission spectrum line represented by the black line in Figure 2a.

Based on the analysis process above, it can be understood that in the absence of the dark mode, the energy of the five bright modes remains independently localized. However, with the inclusion of the dark mode, the destructive interferences between the dark mode and the five bright modes give rise to the formation of quadruple PIT. The originally localized bright modes undergo diffusion and transfer under the influence of the dark mode, resulting in the enhanced energy of the dark mode. This mechanism plays a crucial role in the emergence of the quadruple PIT effect.

4. Math Method

To validate the FDTD simulation results, the CMT is employed to fit the FDTD simulation results. The theoretical coupling diagram derived from the CMT is depicted in Figure 9, where elements A, B, C, D, and E represent the five resonators (corresponding to M2, M3, M4, M5, and DVRG, respectively), and the amplitudes of the resonator modes are labeled as “a, b, c, d, and e”, respectively. The superscripts “in” and “out” of “$A_{\pm }^{\mathrm{in}/\mathrm{out}}, B_{\pm }^{\mathrm{in}/\mathrm{out}}, C_{\pm }^{\mathrm{in}/\mathrm{out}}, D_{\pm }^{\mathrm{in}/\mathrm{out}}, \mathrm{and} E_{\pm }^{\mathrm{in}/\mathrm{out}}$” are respectively the input and output waves of the resonator. The subscripts “+” and “−” of “$A_{\pm }^{\mathrm{in}/\mathrm{out}}, B_{\pm }^{\mathrm{in}/\mathrm{out}}, C_{\pm }^{\mathrm{in}/\mathrm{out}}, D_{\pm }^{\mathrm{in}/\mathrm{out}}, \mathrm{and} E_{\pm }^{\mathrm{in}/\mathrm{out}}$” represent the positive and negative directions of wave propagation, respectively. Drawing upon the concepts of CMT [43,44], the mathematical depiction that elucidates the interconnections among the various resonant modes finds its expression as follows:

$$\left(\begin{array}{ccccc}{\gamma }_1& -i{\mu }_{12}& -i{\mu }_{13}& -i{\mu }_{14}& -i{\mu }_{15}\\ -i{\mu }_{21}& {\gamma }_2& -i{\mu }_{23}& -i{\mu }_{24}& -i{\mu }_{25}\\ -i{\mu }_{31}& -i{\mu }_{32}& {\gamma }_3& -i{\mu }_{34}& -i{\mu }_{45}\\ -i{\mu }_{41}& -i{\mu }_{42}& -i{\mu }_{43}& {\gamma }_4& -i{\mu }_{45}\\ -i{\mu }_{51}& -i{\mu }_{52}& -i{\mu }_{53}& -i{\mu }_{54}& {\gamma }_5\end{array}\right).\left(\begin{array}{c}a\\ b\\ c\\ d\\ e\end{array}\right)=\left(\begin{array}{ccccc}-{\gamma }_{\mathrm{o}1}^{1/2}& 0& 0& 0& 0\\ 0& -{\gamma }_{o2}^{1/2}& 0& 0& 0\\ 0& 0& -{\gamma }_{o3}^{1/2}& 0& 0\\ 0& 0& 0& -{\gamma }_{o4}^{1/2}& 0\\ 0& 0& 0& 0& -{\gamma }_{o5}^{1/2}\end{array}\right).\left(\begin{array}{c}{A}_{+}^{in}+A_{-}^{in}\\ B_{+}^{in}+B_{-}^{in}\\ C_{+}^{in}+C_{-}^{in}\\ D_{+}^{in}+D_{-}^{in}\\ E_{+}^{in}+E_{-}^{in}\end{array}\right)$$

(7)

where ${\mu }_{mn}$ denotes the mutual coupling coefficients that interconnect the five resonators $(m, n=1, 2, 3, 4, 5, m\ne n)$. The variable ${\gamma }_n$ (with $n$ values of 1, 2, 3, 4, 5) is contingent upon the resonant angular frequency ${\omega }_n$ and is explicitly characterized by the equation ${\gamma }_n=i\omega -i{\omega }_n-{\gamma }_{in}-{\gamma }_{on}$ (for $n=1, 2, 3, 4, 5$), where the inter-loss coefficient ${\gamma }_{in}$ is determined by ${\gamma }_{in}={\omega }_n/(2Q_{in})$ and the external loss coefficient ${\gamma }_{on}$ is obtained by ${\gamma }_{on}={\omega }_n/(2Q_{on})$. In addition, the overall quality factor $Q_{tn}$, the inherent loss quality factor $Q_{in}$, and the exterior loss quality factor $Q_{on}$ for the nth resonant mode adhere to the correlation $1/Q_{tn}=1/Q_{in}+1/Q_{on} \left(n=1, 2, 3, 4, 5\right)$, where $Q_{tn}=f/\Delta f$ represents the complete quality factor of the nth mode ($f$ stands for resonance frequency and $\Delta f$ signifies the full width at half maximum). The calculation of $Q_{in}$ involves the formulation $Q_{in}=Re\left(n_{eff}\right) / Im\left(n_{eff}\right)$ [40]. Abiding by the law of energy conservation, the correlation between the input wave and the output wave can be articulated as follows:

$$B_{+}^{\mathrm{in}}=A_{+}^{\mathrm{out}}\cdot \exp (i{\phi }_1),A_{-}^{\mathrm{in}}=B_{-}^{\mathrm{out}}\cdot \exp (i{\phi }_1),$$

(8)

$$C_{+}^{\mathrm{in}}=B_{+}^{\mathrm{out}}\cdot \exp (i{\phi }_2),B_{-}^{\mathrm{in}}=C_{-}^{\mathrm{out}}\cdot \exp (i{\phi }_2),$$

(9)

$$D_{+}^{\mathrm{in}}=C_{+}^{\mathrm{out}}\cdot \exp (i{\phi }_3),C_{-}^{\mathrm{in}}=D_{-}^{\mathrm{out}}\cdot \exp (i{\phi }_3),$$

(10)

$$E_{+}^{\mathrm{in}}=D_{+}^{\mathrm{out}}\cdot \exp (i{\phi }_4),D_{-}^{\mathrm{in}}=E_{-}^{\mathrm{out}}\cdot \exp (i{\phi }_4),$$

(11)

$$A_{+}^{out}=A_{+}^{in}-{\gamma }_{o1}^{1/2}\cdot a,A_{-}^{out}=A_{-}^{in}-{\gamma }_{o1}^{1/2}\cdot a,$$

(12)

$$B_{+}^{out}=B_{+}^{in}-{\gamma }_{o2}^{1/2}\cdot b,B_{-}^{out}=B_{-}^{in}-{\gamma }_{o2}^{1/2}\cdot b,$$

(13)

$$C_{+}^{out}=C_{+}^{in}-{\gamma }_{o3}^{1/2}\cdot c,C_{-}^{out}=C_{-}^{in}-{\gamma }_{o3}^{1/2}\cdot c,$$

(14)

$$D_{+}^{out}=D_{+}^{in}-{\gamma }_{o4}^{1/2}\cdot d,D_{-}^{out}=D_{-}^{in}-{\gamma }_{o4}^{1/2}\cdot d,$$

(15)

$$E_{+}^{out}=E_{+}^{in}-{\gamma }_{o5}^{1/2}\cdot e,E_{-}^{out}=E_{-}^{in}-{\gamma }_{o5}^{1/2}\cdot e,$$

(16)

where ${\varphi }_1$, ${\varphi }_2$, ${\varphi }_3$, and ${\varphi }_4$ denote the phase differences between A and B, B and C, C and D, D and E, respectively. Meanwhile, the values of ${\varphi }_n \left(n=1, 2, 3, 4\right)$ can be derived utilizing the association ${\varphi }_n=Re\left(\beta \right)d_n=\omega d_{n}Re\left(n_{eff}\right)/c$, where $d_n$ designates the distance between layers of graphene.

In addition, there is an absence of input wave propagating in the opposing direction in $E$, thus

$$E_{-}^{in}=0.$$

(17)

According to Equations (7)–(17), the complete system’s transmission coefficient can be represented as follows:

$$\begin{array}{c}t=\frac{E_{+}^{out}}{A_{+}^{in}}=1+\frac{{\gamma }_{\mathrm{o}5}}{{\gamma }_5}+\left(\sqrt{{\gamma }_{o1}}+\sqrt{{\gamma }_{o5}}\frac{x_{51}}{r_5}\right)t_1+\left(\sqrt{{\gamma }_{o2}}+\sqrt{{\gamma }_{o5}}\frac{x_{52}}{r_5}\right)t_2\\ +\left(\sqrt{{\gamma }_{o3}}+\sqrt{{\gamma }_{o5}}\frac{x_{53}}{r_5}\right)t_3+\left(\sqrt{{\gamma }_{o4}}+\sqrt{{\gamma }_{o5}}\frac{x_{54}}{r_5}\right)t_4\end{array}$$

(18)

where

$$t_1=-\frac{t_a}{t_c},$$

(19)

$$t_2=-\frac{t_b}{t_c},$$

(20)

$$t_3=\frac{lt+ps}{kp-lo}+\frac{ip-lm}{kp-lo}.\frac{t_a}{t_c}+\frac{jp-ln}{kp-lo}.\frac{t_b}{t_c},$$

(21)

$$t_4=\frac{o(lt+ps)-t(lo-kp)}{p(kp-lo)}+\frac{m(lo-kp)+o(ip-lm)}{p(kp-lo)}.\frac{t_a}{t_c}+\frac{n(lo-kp)+o(jp-ln)}{p(kp-lo)}.\frac{t_b}{t_c},$$

(22)

$$\begin{array}{ll}{t}_a& =\left[\left(hs+lr\right)\left(ch+dg\right)-\left(hq-dr\right)\left(hk-gl\right)\right]\times \left[\left(kp-lo\right)\left(fl+jh\right)-\left(jp-ln\right)\left(hk-gl\right)\right]\\ & +\left[\left(\mathrm{df}+bh\right)\left(hk-gl\right)+\left(ch+dg\right)\left(fl+jh\right)\right]\times \left[\left(kp-lo\right)\left(el+hi\right)-\left(ip-lm\right)\left(hk-gl\right)\right],\end{array}$$

(23)

$$\begin{array}{l}a={\gamma }_1{x}_{25}+x_{15}{x}_{21},b={\gamma }_2{x}_{15}+x_{12}{x}_{25},c=x_{15}{x}_{23}+x_{13}{x}_{25},d=x_{15}{x}_{24}+x_{14}{x}_{25},\\ e=x_{25}{x}_{31}+x_{21}{x}_{35},f={\gamma }_2{x}_{35}+x_{25}{x}_{32},g={\gamma }_3{x}_{25}+x_{23}{x}_{35},h=x_{25}{x}_{34}+x_{24}{x}_{35},\\ i=x_{35}{x}_{41}+x_{31}{x}_{45},j=x_{35}{x}_{42}+x_{32}{x}_{45},k={\gamma }_3{x}_{45}+x_{35}{x}_{43},l={\gamma }_4{x}_{35}+x_{34}{x}_{45},\\ m={\gamma }_5{x}_{41}+x_{45}{x}_{51},n={\gamma }_5{x}_{42}+x_{42}{x}_{52},o={\gamma }_5{x}_{43}+x_{45}{x}_{53},p={\gamma }_4{\gamma }_5+x_{45}{x}_{54},\\ q=\sqrt{{\gamma }_{o1}}{x}_{25}-\sqrt{{\gamma }_{o2}}{x}_{15},r=\sqrt{{\gamma }_{o2}}{x}_{35}-\sqrt{{\gamma }_{o3}}{x}_{25},s=\sqrt{{\gamma }_{o3}}{x}_{45}-\sqrt{{\gamma }_{o4}}{x}_{35},t=\sqrt{{\gamma }_{o4}}{\gamma }_5-\sqrt{{\gamma }_{o5}}{x}_{45},\end{array}$$

(24)

$$x_{mn}=\sqrt{{\gamma }_{om}{\gamma }_{on}}+i{\mu }_{mn}(m,n=1,2,3,4,5;m\ne n).$$

(25)

Theoretical transmission characteristics closely associated with the proposed graphene metamaterial structure can be mathematically expressed as $T=|t{|}^2$.

5. Analysis and Discussion

As is well known, graphene exhibits satisfactory dynamic tunability. The gate voltage dependency enables dynamic tuning of the transparency window without altering the geometric structure of graphene [45,46,47]. From Figure 10a,b, it becomes evident that the transmission spectra exhibit a nearly linear blue shift with the increase in the Fermi level from 0.6 eV to 0.72 eV. The correlation between the bias voltage and Fermi levels is described by the following equation [30,48]:

$$E_F=\hslash v_F{\left(\frac{\pi {\epsilon }_0{\epsilon }_d{V}_g}{d_{c}e}\right)}^{\frac{1}{2}}{}_{{}^{{}^{}}}$$

(26)

where $v_F$, $V_g$, ${\epsilon }_0$, ${\epsilon }_d$, and $d_c$ correspond to the Fermi carrier velocity, the gate voltage, the relative dielectric constant of air, the relative dielectric constant of silicon and the thickness of silicon, respectively.

The commonly used metrics for evaluating the performance of optical switch modulators are the MD [49] and the IL [50], which can be obtained as follows:

$$MD=\frac{|T_{\max }-T_{\min }|}{T_{\max }}\times 100\%$$

(27)

$$\mathrm{IL}=-{10\mathrm{lgT}}_{\max }$$

(28)

where $T_{\max }$ signifies the maximum transmission and $T_{\min }$ corresponds to the minimum transmission.

The transmission spectra curve, acquired via FDTD simulations, is depicted in Figure 10a. With an elevation of the graphene Fermi level from 0.6 eV to 0.72 eV, a noticeable blue shift becomes apparent along the red curve. The green dashed line represents the results fitted using the CMT theory, which closely matches the red curve. To observe the blue shift more clearly, the variation process with the increment of the $E_F$ is illustrated in Figure 10b. It is shown that the functional relationship between the frequency and transmittance of the structure changes as $E_F$ changes continuously within a certain range. The red curve represents the baseline transmission spectra of the quadruple PIT at 0.6 eV. Then, the overall graphene structure’s $E_F$ is adjusted to 0.72 eV. As demonstrated in Figure 10c, when the $E_F$ is 0.6 eV, the optical switch is considered “off” at the frequency points of 3.74 THz, 4.32 THz, 4.98 THz, and 5.94 THz, while it is “on” at the frequency points of 3.32 THz, 4.08 THz, 4.70 THz, and 5.41 THz. After modulation, the original “off” states at the frequency points of 3.74 THz, 4.32 THz, 4.98 THz, and 5.94 THz change to the “on” states. Conversely, at frequency points of 3.32 THz, 4.08 THz, 4.70 THz, and 5.41 THz, they change to the “off”. Before and after the modulation, the “on” and “off” states are completely reversed. Therefore, the influence of this optical switch is asynchronous.

The MD values and IL values in Figure 10c are 94.0%, 92.48%, 93.54%, 96.54%, 97.51%, 92.86%, 94.82%, and 88.20%, and 0.52 dB, 0.98 dB, 1.37 dB, 0.70 dB, 0.43 dB, 0.63 dB, 0.16 dB, and 0.17 dB, respectively. Here, to demonstrate the superior performance of our proposed modulator, we compare it with several graphene-based metamaterial structures used for optical switch modulators, as shown in

Table 1. Comparison of graphene-based optical switch performance.

Reference/Year	Material Structure	Modulation Mode	MD (%)	IL (dB and %)	Modulation Band
[49]/2020	Multilayer graphene metamaterials	Multiple frequency	87.8%	0.31 dB	Terahertz
			77.7%	0.20 dB
[44]/2020	Multilayer graphene metamaterials	Single frequency	83.3%	0.33 dB	Terahertz
[51]/2020	Multilayer graphene metamaterials	Single frequency	83.3%	7.2%	Terahertz
[52]/2020	Single-layer patterned graphene	Dual frequency	93.0%	0.32 dB	Terahertz
			85.0%	0.25 dB
[53]/2021	Bilayer graphene metamaterials	Triple frequency	86.1%	8.1%	Terahertz
This work	Double-layer patterned graphene	Octa-frequency	94.00%	0.52 dB	Terahertz
			92.48%	0.98 dB
			93.54%	1.37 dB
			96.54%	0.70 dB
			97.51%	0.43 dB
			92.86%	0.63 dB
			94.82%	0.16 dB
			88.20%	0.17 dB

. Clearly, the proposed octa-frequency switch modulator exhibits excellent MD and lower IL performances, providing important guidance for the switch modulation at terahertz frequencies. Overall, this structure exhibits dynamic and controllable spectral response characteristics. Through modulation of the entire structure’s Fermi level, an octa-frequency asynchronous switch can be achieved with an MD of up to 97.51% and excellent modulation performance.

6. Conclusions

In essence, our study introduces a dual-layer graphene terahertz metamaterial structure to realize an octa-frequency asynchronous switch, since quadruple PIT transparency windows have been generated. The mechanism of PIT formation is studied by investigating the electric field distribution. The results derived from the FDTD simulation and CMT fitting show high consistency. Through the manipulation of the Fermi level in graphene, dynamic and controllable multi-switch modulation functionality is achieved. For the proposed octa-frequency switch, the MD is generally greater than 88%, and the IL (0.16 dB < IL < 1.37 dB) is also remarkable. In the case of $E_F$ = 0.6 eV ($E_F$ = 0.72 eV), the maximum MD reaches 97.51%, corresponding to an IL of 0.40 dB. Moreover, the minimum IL of 0.16 dB corresponds to an MD of 94.82%. Compared to other structures based on PIT, our proposed dynamically controllable optical switch exhibits excellent modulation performance in terms of MD. Therefore, we believe that this structure offers a fresh avenue for exploring optical switches and modulators at terahertz frequency.