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

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.

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.
Nanomaterials 2023, 13, 2474 2 of 16 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.

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: 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.
ble 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]: ( ) 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]: where e, , k B , ω, 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.µ is carrier mobility, and v F = 10 6 ms −1 denotes the Fermi velocity.From τ = µE F /ev 2 F [39], the carrier relaxation time τ can be deduced.The ambient temperature is configured as T = 300 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, σ inter can be disregarded, so σ can be obtained [40,41]: In the proposed graphene metamaterial structure, with silicon as the substrate, the propagation constant β and effective refractive index n e f f are as follows [38,42]: ) where ε si , ε air , and ε 0 present the relative dielectric constants of silicon, air and vacuum dielectric, respectively.k 0 and η 0 correspond to the wave number in free space and the inherent impedance, respectively.

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 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.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

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 " " 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: 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.

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 in/out

±
, and E in/out ± " are respectively the input and output waves of the resonator.The subscripts "+" and "−" of "A in/out ± , B in/out ± , C in/out ± , D in/out ± , and E in/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: where µ mn denotes the mutual coupling coefficients that interconnect the five resonators (m, n = 1, 2, 3, 4, 5, m = n).The variable γ n (with n values of 1, 2, 3, 4, 5) is contingent upon the resonant angular frequency ω n and is explicitly characterized by the equation , where the inter-loss coefficient γ in is determined by γ in = ω n /(2Q in ) and the external loss coefficient γ on is obtained by γ on = ω 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 , where Q tn = f /∆ f represents the complete quality factor of the nth mode ( f stands for resonance frequency and ∆ f signifies the full width at half maximum).The calculation of Q in involves the formulation [40].Abiding by the law of energy conservation, the correlation between the input wave and the output wave can be articulated as follows: where φ 1 , φ 2 , φ 3 , and φ

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]: where F v , g V , 0 ε , d ε , and c d correspond to the Fermi carrier velocity, the gate voltage, the relative dielectric constant of air, the relative dielectric constant of silicon and the In addition, there is an absence of input wave propagating in the opposing direction in E, thus According to Equations ( 7)-( 17), the complete system's transmission coefficient can be represented as follows: x 53 where x mn = √ γ om γ on + iµ mn (m, n = 1, 2, 3, 4, 5; m = n).
Theoretical transmission characteristics closely associated with the proposed graphene metamaterial structure can be mathematically expressed as T =|t| 2 .

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]: where v F , V g , ε 0 , ε 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: IL = − 10lgT 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.graphene metamaterial structure.

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]: where F v , g V , 0 ε , d ε , and c d 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.
Nanomaterials 2023, 13, x FOR PEER REVIEW 13 of 17 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: where max T signifies the maximum transmission and min T corresponds to the mini- 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.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.

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.

Figure 1 .
Figure 1.(a) A 3D diagram presents the entirety of the proposed graphene structure.(b) An overhead perspective providing a top view of the graphene arrangement.(c) Diagram of the graphene structure on the upper layer in the z-direction of the structure.(d) Diagram of the graphene structure on the lower layer in the z-direction of the structure.

Figure 1 .
Figure 1.(a) A 3D diagram presents the entirety of the proposed graphene structure.(b) An overhead perspective providing a top view of the graphene arrangement.(c) Diagram of the graphene structure on the upper layer in the z-direction of the structure.(d) Diagram of the graphene structure on the lower layer in the z-direction of the structure.
curve.At the resonant frequencies of the bright modes ( 0 f , 1 f , 2 f , 3 f , and 4 f ), the originally bright modes transform into the states with transmission rates close to 1, as confirmed in Figure 2a,b.

Figure 2 .
Figure 2. (a) The transmission spectra of individual M1 (black line), M2 (green line), M3 (dark blue line), M4 (purple line), and M5 (light blue line).(b) The transmission spectra of the entire structure when the structure is under direct incident light excitation.Here, the EF = 0.6 eV.
represents the spatial distribution of the electric field within the entire structure at the first resonant peak a f .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 fre-

Figure 2 . 17 Figure 3 .
Figure 2. (a) The transmission spectra of individual M1 (black line), M2 (green line), M3 (dark blue line), M4 (purple line), and M5 (light blue line).(b) The transmission spectra of the entire structure when the structure is under direct incident light excitation.Here, the E F = 0.6 eV.Nanomaterials 2023, 13, x FOR PEER REVIEW 7 of 17

Figure 3 .
Figure 3. (a) Transmission spectra for the scenarios involving only M2, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M2 and DVRG at the resonant frequency f 1 .(d) Electric field distribution of the entire structure at f a .

Figure 3 .
Figure 3. (a) Transmission spectra for the scenarios involving only M2, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M2 and DVRG at the resonant frequency f1.(d) Electric field distribution of the entire structure at fa.

Figure 4 .
Figure 4. (a) Transmission spectra for the scenarios involving only M3, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M3 and DVRG at the resonant frequency f2.(d) Electric field distribution of the entire structure at fb.

Figure 4 . 17 Figure 5 .
Figure 4. (a) Transmission spectra for the scenarios involving only M3, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M3 and DVRG at the resonant frequency f 2 .(d) Electric field distribution of the entire structure at f b .Nanomaterials 2023, 13, x FOR PEER REVIEW 8 of 17

Figure 5 .
Figure 5. (a) Transmission spectra for the scenarios involving only M4, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M4 and DVRG at the resonant frequency f 3 .(d) Electric field distribution of the entire structure at f c .

Figure 5 .
Figure 5. (a) Transmission spectra for the scenarios involving only M4, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M4 and DVRG at the resonant frequency f3.(d) Electric field distribution of the entire structure at fc.

Figure 6 .
Figure 6.(a) Transmission spectra for the scenarios involving only M5, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M5 and DVRG at the resonant frequency f4.(d) Electric field distribution of the entire structure at fd.

Figure 6 . 17 Figure 7 .
Figure 6.(a) spectra for the scenarios involving only M5, DVRG, and the entire structure.(b,c) The spatial distribution of the electric field within structures M5 and DVRG at the resonant frequency f 4 .(d) Electric field distribution of the entire structure at f d .Nanomaterials 2023, 13, x FOR PEER REVIEW 9 of 17

Figure 7 .
Figure 7. (a) Transmission spectra for the scenarios involving only M1, DVRG, and the entire structure.(b-d) The spatial distribution of the electric field within structures M4, DVRG, and the entire structure at the resonant frequency f 0 .

Figure 7 .
Figure 7. (a) Transmission spectra for the scenarios involving only M1, DVRG, and the entire structure.(b-d) The spatial distribution of the electric field within structures M4, DVRG, and the entire structure at the resonant frequency f0.

Figure 8 .
Figure 8.The spatial distribution of the electric field within structure M1 at 5.64 THz.

Figure 8 .
Figure 8.The spatial distribution of the electric field within structure M1 at 5.64 THz.
4 denote the phase differences between A and B, B and C, C and D, D and E, respectively.Meanwhile, the values of φ n (n = 1, 2, 3, 4) can be derived utilizing the association φ n = Re(β)d n = ωd n Re n e f f /c, where d n designates the distance between layers of graphene.Nanomaterials 2023, 13, x FOR PEER REVIEW 12 of 17

Figure 9 .
Figure 9. Theoretical schematic depicting the coupling among the resonant modes in the proposed graphene metamaterial structure.

Figure 9 .
Figure 9. Theoretical schematic depicting the coupling among the resonant modes in the proposed graphene metamaterial structure.

Figure 10 .
Figure 10.(a) The simulated transmission spectrum obtained from FDTD at different EF.(b) The theoretical transmission spectrum as EF undergoes continuous variation within the range of 0.6 eV to 0.72 eV in FDTD simulation.(c) The octa-frequency switch of the transmission spectrum under EF values of 0.60 eV and 0.72 eV, respectively.

Figure 10 .
Figure 10.(a) The simulated transmission spectrum obtained from FDTD at different E F .(b) The theoretical transmission spectrum as E F undergoes continuous variation within the range of 0.6 eV to 0.72 eV in FDTD simulation.(c) The octa-frequency switch of the transmission spectrum under E F values of 0.60 eV and 0.72 eV, respectively.

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