Graphene-Modulated Terahertz Metasurfaces for Selective and Active Control of Dual-Band Electromagnetic Induced Reflection (EIR) Windows

Currently, metasurfaces (MSs) integrating with different active materials have been widely explored to actively manipulate the resonance intensity of multi-band electromagnetic induced transparency (EIT) windows. Unfortunately, these hybrid MSs can only realize the global control of multi-EIT windows rather than selective control. Here, a graphene-functionalized complementary terahertz MS, composed of a dipole slot and two graphene-integrated quadrupole slots with different sizes, is proposed to execute selective and active control of dual-band electromagnetic induced reflection (EIR) windows. In this structure, dual-band EIR windows arise from the destructive interference caused by the near field coupling between the bright dipole slot and dark quadrupole slot. By embedding graphene ribbons beneath two quadrupole slots, the resonance intensity of two windows can be selectively and actively modulated by adjusting Fermi energy of the corresponding graphene ribbons via electrostatic doping. The theoretical model and field distributions demonstrate that the active tuning behavior can be ascribed to the change in the damper factor of the corresponding dark mode. In addition, the active control of the group delay is further investigated to develop compact slow light devices. Therefore, the selective and active control scheme introduced here can offer new opportunities and platforms for designing multifunctional terahertz devices.


Introduction
Over the last few decades, the electromagnetic induced transparency (EIT) first observed in the atomic system has aroused a lot of interest because it allows the enhancement of transmission within a narrow frequency range along with strong dispersion properties and displays promising engineering prospects [1][2][3]. However, its implementation requires rigorous experimental conditions, severely hindering practical applications [4]. With the rapid development of metasurfaces (MSs), recently, the MS-based EIT-like effect has attracted tremendous attention owing to flexible design and easy fabrication [5]. Currently, the EIT-like effect in terahertz (THz) MSs has been widely investigated and experimentally demonstrated by bright-dark mode coupling [6,7] or bright-bright mode coupling [8,9]. Unfortunately, most of these structures can only realize a single transparency window, in which the resonance strength of the window is fixed after fabrication, restricting their application fields [10].
To overcome the above limitation, a variety of THz MSs with multi-band EIT windows were designed to actively modulate resonance properties by means of different approaches, including mechanical movement, phase change material, and graphene. For example, Devi et al. designed a concentrically coupled asymmetric THz MS to obtain a dual-band EIT effect; moreover, the dual-band EIT effect can be effectively modulated by rotating the element angle or by changing the asymmetric degree of the unit cell [11]. Chen et al. demonstrated an actively tunable dual-band EIT window in a VO 2 -integrated MSs [12]. In addition, the last decade saw the advent of graphene with flexible tunability, which has given rise to unprecedented progress in actively tunable microwave [13,14], terahertz [15,16], and optical [17] devices. According to similar principles, Gao et al. reported an actively tunable dual plasmon-induced transparency (PIT) MS based on graphenepatterned structures [18]. Meanwhile, three-band tunable EIT MSs (or those with more than three bands) were also investigated by stacking multilayer graphene/dielectric structures [19,20]. Among all of the previously mentioned MS structures, they primarily rely on uniform excitation to implement the global control of multi-band windows rather than selective control; that is, each one of the multi-band windows cannot be independently controlled or modulated. With the increasing demand for selective and tunable control systems, however, simple and multifunctional MSs are highly desirable for reconfigurable THz devices in wireless communications and data storage systems [21][22][23].
Here, we propose a scheme to obtain a dual-band EIR window in a graphenefunctionalized complementary THz MS consisting of two dolmen-like slot structures with different sizes and two graphene ribbons connecting to different top gates. The basic article structure is as follows: First, the structural design and parameter selection of the proposed MS are described in detail, as well as the simulation setting. Second, the formation process and the physical mechanism of two EIR windows are elucidated through the analysis of different element arrays in a unit cell and field distributions. Then, the selective and active control of two EIR windows are further investigated by tailoring the Fermi energy of the corresponding graphene ribbons via electrostatic doping. Moreover, the field distributions and classic three-oscillator model are employed to discover the selective and active control mechanism of two EIR windows. Finally, the active control of the group delay originating from strong phase dispersion is also investigated to develop compact slow light devices. Therefore, the results obtained for our devices could provide new opportunities and platforms for the design and development of multi-function terahertz devices. Figure 1 shows the schematic of the proposed graphene-functionalized complementary THz MS structure, in which the metal Al film, graphene layer, dielectric layer, and substrate are represented by yellow, blue, red, and green, respectively. The unit cell is composed of the metallic slot structures and two patterned monolayer graphene ribbons, and the metallic slot structures comprise a dipole slot serving as a bright resonator and two quadrupole slots of different sizes serving as two different dark resonators fabricated by etching the Al film, thus forming two dolmen-like slot structures with a communal dipole slot, as shown in Figure 1a,b. When the designed structure is excited by the x-polarized incident wave, two distinct EIR windows would be induced by tuning the near field coupling between the corresponding bright-dark modes. To realize selective control of the two EIR windows, two quadrupole slots are mutually separated by a small gap in the center of the dipole slot and are deposited on the graphene ribbons connected by different top gates, as shown in Figure 1c. Compared to traditional structures, the key feature of our structure is that it can realize the selective and active control of the two windows by adjusting the Fermi energy of the corresponding graphene ribbon via electric doping. To investigate the THz response of the proposed complementary MS, numerical calculations were conducted using the FDTD method, in which the boundary conditions in the x-and y-directions were set as the unit cell, and the z-direction was left open. The xpolarized incident light was normally illuminated on the surface of the designed complementary MS along the z-direction, as displayed in Figure 1a. Moreover, the size and number of the rational mesh were used to meet the accuracy requirements. In addition, the lossy metal Al with the DC conductivity of 3.56 × 10 7 S/m was employed as the slot resonators during the calculation [24], while the low-doped silicon with a permittivity of 11.7 and thickness of 300 μm and a SiO2 film with a permittivity of 3.9 and thickness of 30 μm were used as the substrate and dielectric layer, respectively. Since the graphene is composed of a one atom-thick layer of carbon atoms arranged in a hexagonal pattern, generally, the graphene is assumed to be a 1.0-nm-thick homogenous layer in order to facilitate simulation [25], while the effective permittivity is written as [26]:

Structures, Materials and Methods
Here, ω is the angular frequency, tg is the single-layer graphene thickness, and ε0 is the vacuum permittivity. In the terahertz range, the graphene conductivity σg is described by the following expression [27]: where e is the electron charge, is the reduced Planck's constant, kB is the Boltzmann's constant, T is the Kelvin temperature, and EF is the Fermi energy of graphene. τ is the carrier relaxation time described by = ( )/( 2 ) ; here, T = 300 K, = 3000 cm 2 /V • s, and = 1.1 × 10 6 m/s are used in our calculations [28]. According to Equation (2), the conductivity of graphene can be actively tuned by shifting the Fermi energy of graphene, while the Fermi energy of the graphene can be dynamically controlled by the voltages of the top gates (as shown in Figure 1c); that is, the conductivity of the single-layer graphene can be actively modulated by the top-gate voltages, which is different from the electro-intercalation way applied in multilayer graphenes [29], and the corresponding expression is as follows [30]:  Figure 1. The proposed graphene-functionalized EIR metasurface: (a) three-dimensional schematic of the EIR metasurface, (b) schematic of the unit cell; the geometric parameters are as follows: l 1 = 120 µm, l 2 = 109 µm, l 3 = 100 µm, w 1 = 36 µm, w 2 = 28 µm, w 3 = 35 µm, s 1 = 30 µm, s 2 = 20 µm, g 1 = 2 µm, g 2 = 3.4 µm, g 3 = 2 µm, P x = 30 µm and P y = 160 µm, and (c) cross-sectional view of unit cell.
To investigate the THz response of the proposed complementary MS, numerical calculations were conducted using the FDTD method, in which the boundary conditions in the x-and y-directions were set as the unit cell, and the z-direction was left open. The x-polarized incident light was normally illuminated on the surface of the designed complementary MS along the z-direction, as displayed in Figure 1a. Moreover, the size and number of the rational mesh were used to meet the accuracy requirements. In addition, the lossy metal Al with the DC conductivity of 3.56 × 10 7 S/m was employed as the slot resonators during the calculation [24], while the low-doped silicon with a permittivity of 11.7 and thickness of 300 µm and a SiO 2 film with a permittivity of 3.9 and thickness of 30 µm were used as the substrate and dielectric layer, respectively. Since the graphene is composed of a one atom-thick layer of carbon atoms arranged in a hexagonal pattern, generally, the graphene is assumed to be a 1.0-nm-thick homogenous layer in order to facilitate simulation [25], while the effective permittivity is written as [26]: Here, ω is the angular frequency, t g is the single-layer graphene thickness, and ε 0 is the vacuum permittivity. In the terahertz range, the graphene conductivity σ g is described by the following expression [27]: where e is the electron charge, = h/2π is the reduced Planck's constant, k B is the Boltzmann's constant, T is the Kelvin temperature, and E F is the Fermi energy of graphene. τ is the carrier relaxation time described by τ = (µE F )/ ev 2 F ; here, T = 300 K, µ = 3000 cm 2 /V·s, and v F = 1.1 × 10 6 m/s are used in our calculations [28]. According to Equation (2), the conductivity of graphene can be actively tuned by shifting the Fermi energy of graphene, while the Fermi energy of the graphene can be dynamically controlled by the voltages of the top gates (as shown in Figure 1c); that is, the conductivity of the single-layer graphene can be actively modulated by the top-gate voltages, which is different from the electro-intercalation way applied in multilayer graphenes [29], and the corresponding expression is as follows [30]: Here, ε r and d s are the relative permittivity and thickness of the SiO 2 dielectric layer. Thus, by tuning the conductivity of graphene through the gate voltage V g , the THz properties of the designed device could be effectively controlled.

Evolution Mechanism of Two EIR Windows
To disclose the formation process of dual-band EIR effect, first, the reflection spectra of four different complementary MSs, consisting of a longer quadrupole slot array alone, a dipole slot array alone, and shorter quadrupole slot array alone as well as their combined structure were calculated for the x-polarized incident wave, as presented in Figure 2. For the individual dipole slot, a sharp reflection dip is observed at 0.93 THz due to direct coupling with the x-polarized incident wave, which serves as a bright mode, as shown by the shallow blue line in Figure 2b. For the individual longer or shorter quadrupole slots, by contrast, the reflection amplitude is almost close to 1.0, and no resonance appears in the interesting frequency range due to the structural symmetry, as displayed by the red line in Figure 2a or the orange line in Figure 2c. Thus, the quadrupole slot acts as a dark mode. When three slot elements are combined together to form the unit cell shown in Figure 1b, two distinct EIR windows with a resonance strength of more than 90% in the reflection spectrum appear at 0.86 THz and 0.94 THz due to the destructive interference caused by the near field coupling between the corresponding bright-dark mode, as represented by the deep blue line in Figure 2d.
Here, εr and ds are the relative permittivity and thickness of the SiO2 dielectric layer. Thus, by tuning the conductivity of graphene through the gate voltage Vg, the THz properties of the designed device could be effectively controlled.

Evolution Mechanism of two EIR Windows
To disclose the formation process of dual-band EIR effect, first, the reflection spectra of four different complementary MSs, consisting of a longer quadrupole slot array alone, a dipole slot array alone, and shorter quadrupole slot array alone as well as their combined structure were calculated for the x-polarized incident wave, as presented in Figure 2. For the individual dipole slot, a sharp reflection dip is observed at 0.93 THz due to direct coupling with the x-polarized incident wave, which serves as a bright mode, as shown by the shallow blue line in Figure 2b. For the individual longer or shorter quadrupole slots, by contrast, the reflection amplitude is almost close to 1.0, and no resonance appears in the interesting frequency range due to the structural symmetry, as displayed by the red line in Figure 2a or the orange line in Figure 2c. Thus, the quadrupole slot acts as a dark mode. When three slot elements are combined together to form the unit cell shown in Figure 1b, two distinct EIR windows with a resonance strength of more than 90% in the reflection spectrum appear at 0.86 THz and 0.94 THz due to the destructive interference caused by the near field coupling between the corresponding bright-dark mode, as represented by the deep blue line in Figure 2d. To further understand the physical mechanism of the two EIR windows, the field distributions of the designed structure were calculated when the incident wave was polarized along the x-direction, and the corresponding field distributions at three different resonance frequencies, marked as the symbols "A", "B", and "C" in Figure 2b,d, are displayed in Figure 3. For the dipole slot alone in Figure 3a,d, the magnetic fields are strongly focused on both ends of the dipole slot, while the electric fields are concentrated in the center, indicating the excitation of a dipole resonance, as presented in Figure 2b. For the combined structure, the magnetic and electric fields at two resonance peaks are confined to both ends and the center of the corresponding quadrupole slot in which the fields are To further understand the physical mechanism of the two EIR windows, the field distributions of the designed structure were calculated when the incident wave was polarized along the x-direction, and the corresponding field distributions at three different resonance frequencies, marked as the symbols "A", "B", and "C" in Figure 2b,d, are displayed in Figure 3. For the dipole slot alone in Figure 3a,d, the magnetic fields are strongly focused on both ends of the dipole slot, while the electric fields are concentrated in the center, indicating the excitation of a dipole resonance, as presented in Figure 2b. For the combined structure, the magnetic and electric fields at two resonance peaks are confined to both ends and the center of the corresponding quadrupole slot in which the fields are opposite to each other, while that of the dipole slot is fully restrained, as shown in Figure 3b-f. Such field distributions imply that there is a notable energy shift between the dipole slot and the quadrupole slot due to strong near field coupling; as a result, the two pronounced EIR windows are excited by the destructive interference caused by the near field coupling, as shown in Figure 2d.
opposite to each other, while that of the dipole slot is fully restrained, as shown in Figure  3b-f. Such field distributions imply that there is a notable energy shift between the dipole slot and the quadrupole slot due to strong near field coupling; as a result, the two pronounced EIR windows are excited by the destructive interference caused by the near field coupling, as shown in Figure 2d.

Selective and Active Modulation of Two EIR Windows
To examine the selective and active modulation of the two windows, we next investigated the tunable characteristics of the graphene-functionalized MS by tailoring the Fermi energies of two graphene ribbons, in which the Fermi energy of the graphene ribbons beneath the long and short quadrupole slots are denoted by Ef1 andEf2, respectively. Thus, by independently tuning any of Ef1 and Ef2 or by synchronously tuning Ef1 and Ef2, the graphene-functionalized MS can execute three different ways to control the resonance intensity of the two windows, as shown in Figure 4. As shown in the first row of Figure 4, two distinct EIR windows are observed at 0.86 THz and 0.94 THz without graphene. To facilitate expression, the low-and high-frequency windows are, described as Window I and Window II in the following. Once the graphene ribbons have been embedded in the complementary structure, however, it is noted that the profile features of two windows remain still similar to that of the one no graphene ribbon, with the exception of a slight decline in the strength (not shown here). The slight decrease in strength can be ascribed to the intrinsic conductivity of the graphene, which can shield the near field coupling between the two modes, and similar results have also been observed in previous reports [31]. As the Fermi energy increases in the graphene ribbon, the amplitude of the EIR window should further decline. Figure 4a shows the intensity modulation of Window I when Ef1 is individually tailored and when Ef2 is maintained at 0.1 eV. It can be clearly observed that when Ef1 increases from 0.1 eV to 0.5 eV, the intensity of Window I changes dramatically, while Window II remains unchanged. For example, as Ef1 increases from 0.1 eV to 0.4 eV, the amplitude of Window I decreases from 76% to 56%. With further increases to 0.5 eV, however, Window I is totally switched off and becomes a broad resonance valley. The extinction of Window I can be attributed to the increasing conductivity of the graphene. Thus, Window I can be independently controlled by merely tuning Ef1. To individually tune Ef2 from 0.1 eV to 0.5 eV, a similar change is also observed in Window II, with dependent control of Window II being obtained by Window I shifting Ef2 on its own, as shown in Figure 4b. The intensity change of the two windows that occurs when simultaneously tuning Ef1 and Ef2, however, is shown in Figure 4c. It can be seen that as Ef1 and Ef2 increase, the two windows weaken gradually and are annihilated, eventually becoming a considerable broad resonance dip, different from that of Figure 2b. As a result, the two windows can also realize an on-to-off switch by synchronously shifting Ef1 and Ef2. In addition, it is noteworthy that during the entire modulation process, the locations of the two

Selective and Active Modulation of Two EIR Windows
To examine the selective and active modulation of the two windows, we next investigated the tunable characteristics of the graphene-functionalized MS by tailoring the Fermi energies of two graphene ribbons, in which the Fermi energy of the graphene ribbons beneath the long and short quadrupole slots are denoted by E f1 andE f2 , respectively. Thus, by independently tuning any of E f1 and E f2 or by synchronously tuning E f1 and E f2 , the graphene-functionalized MS can execute three different ways to control the resonance intensity of the two windows, as shown in Figure 4. As shown in the first row of Figure 4, two distinct EIR windows are observed at 0.86 THz and 0.94 THz without graphene. To facilitate expression, the low-and high-frequency windows are, described as Window I and Window II in the following. Once the graphene ribbons have been embedded in the complementary structure, however, it is noted that the profile features of two windows remain still similar to that of the one no graphene ribbon, with the exception of a slight decline in the strength (not shown here). The slight decrease in strength can be ascribed to the intrinsic conductivity of the graphene, which can shield the near field coupling between the two modes, and similar results have also been observed in previous reports [31]. As the Fermi energy increases in the graphene ribbon, the amplitude of the EIR window should further decline. Figure 4a shows the intensity modulation of Window I when E f1 is individually tailored and when E f2 is maintained at 0.1 eV. It can be clearly observed that when E f1 increases from 0.1 eV to 0.5 eV, the intensity of Window I changes dramatically, while Window II remains unchanged. For example, as E f1 increases from 0.1 eV to 0.4 eV, the amplitude of Window I decreases from 76% to 56%. With further increases to 0.5 eV, however, Window I is totally switched off and becomes a broad resonance valley. The extinction of Window I can be attributed to the increasing conductivity of the graphene. Thus, Window I can be independently controlled by merely tuning E f1 . To individually tune E f2 from 0.1 eV to 0.5 eV, a similar change is also observed in Window II, with dependent control of Window II being obtained by Window I shifting E f2 on its own, as shown in Figure 4b. The intensity change of the two windows that occurs when simultaneously tuning E f1 and E f2 , however, is shown in Figure 4c. It can be seen that as E f1 and E f2 increase, the two windows weaken gradually and are annihilated, eventually becoming a considerable broad resonance dip, different from that of Figure 2b. As a result, the two windows can also realize an on-to-off switch by synchronously shifting E f1 and E f2 . In addition, it is noteworthy that during the entire modulation process, the locations of the two windows remain almost unchanged. To evaluate the modulation effect of the three control methods discussed above, the intensity modulation depth of the two windows is calculated by ∆R/R 0 = (R 0 − R g )/R 0 , where R 0 and R g are the reflectance amplitude of the window without and with graphene, respectively. Thus, the corresponding intensity modulation depths are 46% when shifting E f1 alone, 46% when shifting E f2 alone, and more than 40% when simultaneously shifting E f1 and E f2 . Therefore, the two windows can be selectively and actively modulated by tuning the Fermi energy of the corresponding graphene ribbons. Based on this scheme, a MSs with more windows can also be realized, in which each EIR window can be independently and actively controlled.
windows remain almost unchanged. To evaluate the modulation effect of the three control methods discussed above, the intensity modulation depth of the two windows is calculated by ΔR/R0 = (R0 − Rg) / R0, where R0 and Rg are the reflectance amplitude of the window without and with graphene, respectively. Thus, the corresponding intensity modulation depths are 46% when shifting Ef1 alone, 46% when shifting Ef2 alone, and more than 40% when simultaneously shifting Ef1 and Ef2. Therefore, the two windows can be selectively and actively modulated by tuning the Fermi energy of the corresponding graphene ribbons. Based on this scheme, a MSs with more windows can also be realized, in which each EIR window can be independently and actively controlled. To further demonstrate the claims that the above-mentioned two windows can be selectively and actively controlled by shifting the Fermi energy of the corresponding graphene ribbons, the field distributions at two window frequencies are calculated at different Fermi energies, as shown in Figure 5. Figure 5 presents the magnetic and electric field distributions for different Ef1 and Ef2 values obtained by simultaneously doping two graphene ribbons. At Ef1 = 0.10 eV and Ef2 = 0.08 eV, the graphene ribbon can be regarded as a semiconductor with low conductivity, and the confined magnetic and electric fields are strongly focused on the corresponding quadrupole slots (see Figure 5a,d), indicating the occurrence of two windows. However, when Ef1 and Ef2 are gradually increased up to 0.6 eV and 1.0 eV, respectively, the confined fields in the two quadrupole slots initially exhibit considerable attenuation and, ultimately, total annihilation; as a result, this leads to the two windows switching off, as shown in Figures 5b,c,e,f. The change in the magnetic and electric field strengths can be attributed to the increased conductivity of the graphene ribbon, which is able to shield the near field coupling between the bright and dark modes. Here, and in the following section, only the synchronous doping method for the two graphene ribbons is analyzed due to the similar mechanism, while other methods (selective doping) are not further discussed. In addition, it is noted that during the entire tuning process, the confined fields in the dipole slots are also suppressed due to the fact that it increasingly dampens itself, as demonstrated in the last row of Figure 4. To further demonstrate the claims that the above-mentioned two windows can be selectively and actively controlled by shifting the Fermi energy of the corresponding graphene ribbons, the field distributions at two window frequencies are calculated at different Fermi energies, as shown in Figure 5. Figure 5 presents the magnetic and electric field distributions for different E f1 and E f2 values obtained by simultaneously doping two graphene ribbons. At E f1 = 0.10 eV and E f2 = 0.08 eV, the graphene ribbon can be regarded as a semiconductor with low conductivity, and the confined magnetic and electric fields are strongly focused on the corresponding quadrupole slots (see Figure 5a,d), indicating the occurrence of two windows. However, when E f1 and E f2 are gradually increased up to 0.6 eV and 1.0 eV, respectively, the confined fields in the two quadrupole slots initially exhibit considerable attenuation and, ultimately, total annihilation; as a result, this leads to the two windows switching off, as shown in Figure 5b,c,e,f. The change in the magnetic and electric field strengths can be attributed to the increased conductivity of the graphene ribbon, which is able to shield the near field coupling between the bright and dark modes. Here, and in the following section, only the synchronous doping method for the two graphene ribbons is analyzed due to the similar mechanism, while other methods (selective doping) are not further discussed. In addition, it is noted that during the entire tuning process, the confined fields in the dipole slots are also suppressed due to the fact that it increasingly dampens itself, as demonstrated in the last row of Figure 4. To discover tunable mechanism, a classic three-oscillator coupled model was employed to describe the designed graphene-functionalized MS composed of a dipole slot and two graphene-integrated quadrupole slots. In this model, the dipole slot is denoted by the oscillator b being excited directly by the incident wave E, and the two quadrupole slots are denoted by the oscillators a and c, which are being excited indirectly by near field coupling with oscillator b. Thus, the coupled differential equations among the three oscillators are expressed as follows [12]: where (xa, γa), (xb, γb), and (xc, γc) are the resonance intensities and damping factors of three oscillators a, b, and c, respectively. ωa, ωb, and ωc are the resonance frequencies of the three resonators before coupling. κab denotes the coupling coefficient between resonator b and oscillator a, while κbc is the coupling coefficient between resonator b and oscillator c. q represents the coupling coefficient between resonator b and the incident wave. After solving Equations (3)-(5), the energy dissipation of the graphene-functionalized complementary MS as the functions of the frequency is given as follows: where Cj = ωj 2 + iωγj -ω 2 (j = a, b, c). According to Equation (6), the simulated reflection spectrum can be theoretically fitted, and the corresponding fitting results are shown in Figure 6. It can be seen from Figure 6a that for synchronous doping, the theoretical fitting results show reasonable agreement with the simulated results. Moreover, based on a similar fitting method, the fitting and simulated results can also exhibit very excellent agreement for selective doping (not shown here). To discover tunable mechanism, a classic three-oscillator coupled model was employed to describe the designed graphene-functionalized MS composed of a dipole slot and two graphene-integrated quadrupole slots. In this model, the dipole slot is denoted by the oscillator b being excited directly by the incident wave E, and the two quadrupole slots are denoted by the oscillators a and c, which are being excited indirectly by near field coupling with oscillator b. Thus, the coupled differential equations among the three oscillators are expressed as follows [12]: ..
where (x a , γ a ), (x b , γ b ), and (x c , γ c ) are the resonance intensities and damping factors of three oscillators a, b, and c, respectively. ω a , ω b , and ω c are the resonance frequencies of the three resonators before coupling. κ ab denotes the coupling coefficient between resonator b and oscillator a, while κ bc is the coupling coefficient between resonator b and oscillator c. q represents the coupling coefficient between resonator b and the incident wave. After solving Equations (4)-(6), the energy dissipation of the graphene-functionalized complementary MS as the functions of the frequency is given as follows: where C j = ω j 2 + iωγ j − ω 2 (j = a, b, c). According to Equation (7), the simulated reflection spectrum can be theoretically fitted, and the corresponding fitting results are shown in Figure 6. It can be seen from Figure 6a that for synchronous doping, the theoretical fitting results show reasonable agreement with the simulated results. Moreover, based on a similar fitting method, the fitting and simulated results can also exhibit very excellent agreement for selective doping (not shown here). To quantitatively explain the active modulation of the two windows in the synchronous doping situation shown in Figure 4c, the fitting parameters κab, κbc, γa, γb, and γc as the functions of Ef1 and Ef2 are extracted during the fitting calculations and are represented in Figure 6b. During synchronous doping, the coupling coefficients κab and κbc stay nearly constant with Ef1 and Ef2, while the amplitudes of the damping factors γa and γc exhibit a linear increase. When Ef1 and Ef2 are increased up to 0.6 eV and 1.0 eV, respectively, γa and γc become large enough to suppress the excitation of two dark elements, as a result, whichleads to the complete disappearance of the two windows. In addition, it is observed that the damping factor γa arises gradually with the increase of Ef1 and Ef2, which is consistent with the field distributions shown in Figure 5. Therefore, the active modulation of the reflection window can be attributed to the increased damping factor of the dark element arising from the increasing conductivity of the doped graphene ribbon.

Tunable Slow-Light Applications
As it is widely known, one remarkable feature of the EIT effect is that it is accompanied by strong phase dispersion slowing down the light speed, thus exploiting slow light devices. Generally, slow light effect can be quantitatively characterized by a group delay (τg = -dφ /dω, in which φ and ω = 2πf are the phase shift and the angular frequency) [32], as shown in Figure 7. Figure 7a shows the reflection phase spectra for different Ef1 and Ef2. As expected, there is a notable phase dispersion around two EIR windows; moreover, the dispersion in phase shift progressively weakens with Ef1 and Ef2. Figure 7b shows the group delays of the corresponding phases. When there is not graphene ribbon, the group delays near the two EIR windows can be up to 2.6 ps and 3.4 ps, respectively, indicating that light waves propagate through the designed MS sample with a slower group velocity than they do through air. Such a large group delay can be ascribed to the steepest phase shift change within the reflection window. Once the graphene ribbons have been integrated to form a graphene-functionalized MS, however, it can be easily discerned that when Ef1 and Ef2 are 0.30 eV and 0.27 eV, respectively, the two group delays are suppressed to 0.36 ps and 0.32 ps, owing to the phase dispersion being weakened. As Ef1 and Ef2 are further increased up to 0.60 eV and 1.0 eV, the group delay totally disappears and becomes a negative value. Moreover, the selective control of the two delay groups can be also observed by doping the corresponding graphene ribbon (not shown here). By doping the To quantitatively explain the active modulation of the two windows in the synchronous doping situation shown in Figure 4c, the fitting parameters κ ab , κ bc , γ a , γ b, and γ c as the functions of E f1 and E f2 are extracted during the fitting calculations and are represented in Figure 6b. During synchronous doping, the coupling coefficients κ ab and κ bc stay nearly constant with E f1 and E f2 , while the amplitudes of the damping factors γ a and γ c exhibit a linear increase. When E f1 and E f2 are increased up to 0.6 eV and 1.0 eV, respectively, γ a and γ c become large enough to suppress the excitation of two dark elements, as a result, whichleads to the complete disappearance of the two windows. In addition, it is observed that the damping factor γ a arises gradually with the increase of E f1 and E f2 , which is consistent with the field distributions shown in Figure 5. Therefore, the active modulation of the reflection window can be attributed to the increased damping factor of the dark element arising from the increasing conductivity of the doped graphene ribbon.

Tunable Slow-Light Applications
As it is widely known, one remarkable feature of the EIT effect is that it is accompanied by strong phase dispersion slowing down the light speed, thus exploiting slow light devices. Generally, slow light effect can be quantitatively characterized by a group delay (τ g = −dϕ /dω, in which ϕ and ω = 2πf are the phase shift and the angular frequency) [32], as shown in Figure 7. Figure 7a shows the reflection phase spectra for different E f1 and E f2 . As expected, there is a notable phase dispersion around two EIR windows; moreover, the dispersion in phase shift progressively weakens with E f1 and E f2 . Figure 7b shows the group delays of the corresponding phases. When there is not graphene ribbon, the group delays near the two EIR windows can be up to 2.6 ps and 3.4 ps, respectively, indicating that light waves propagate through the designed MS sample with a slower group velocity than they do through air. Such a large group delay can be ascribed to the steepest phase shift change within the reflection window. Once the graphene ribbons have been integrated to form a graphene-functionalized MS, however, it can be easily discerned that when E f1 and E f2 are 0.30 eV and 0.27 eV, respectively, the two group delays are suppressed to 0.36 ps and 0.32 ps, owing to the phase dispersion being weakened. As E f1 and E f2 are further increased up to 0.60 eV and 1.0 eV, the group delay totally disappears and becomes a negative value. Moreover, the selective control of the two delay groups can be also observed by doping the corresponding graphene ribbon (not shown here). By doping the corresponding graphene ribbon, the graphene-functionalized MS can also achieve the ability to selectively control the group delay, which has potential applications for the development of compact slow light devices.
corresponding graphene ribbon, the graphene-functionalized MS can also achieve the ability to selectively control the group delay, which has potential applications for the development of compact slow light devices.

Conclusions
In conclusion, we have numerically demonstrated the selective and active control of two EIR windows in a graphene-functionalized complementary MS consisting of a dipole slot and two graphene-integrated quadrupole slots of different sizes. In the structure, two EIR windows are induced by destructive interference originating from near field coupling between the dipole and the quadrupole slots and can be selectively or synchronously modulated by doping the corresponding graphene ribbons. The theoretical model and field distributions discovered that such active control can be ascribed to the increased damping factor of the dark element. In addition, the group delay is also actively controlled by doping the corresponding graphene ribbons, which is well-suited for the development and application of compact slow light devices. Based on this scheme, the selective and active control of multiple windows can further extend to other functionalized MSs by implanting different dark elements integrated with active materials, which offer new opportunities and platforms for the design of multifunctional terahertz devices.

Conclusions
In conclusion, we have numerically demonstrated the selective and active control of two EIR windows in a graphene-functionalized complementary MS consisting of a dipole slot and two graphene-integrated quadrupole slots of different sizes. In the structure, two EIR windows are induced by destructive interference originating from near field coupling between the dipole and the quadrupole slots and can be selectively or synchronously modulated by doping the corresponding graphene ribbons. The theoretical model and field distributions discovered that such active control can be ascribed to the increased damping factor of the dark element. In addition, the group delay is also actively controlled by doping the corresponding graphene ribbons, which is well-suited for the development and application of compact slow light devices. Based on this scheme, the selective and active control of multiple windows can further extend to other functionalized MSs by implanting different dark elements integrated with active materials, which offer new opportunities and platforms for the design of multifunctional terahertz devices.

Conflicts of Interest:
The authors declare no conflict of interest.