Quadruple plasmon-induced transparency and tunable multi-frequency switch in monolayer graphene terahertz metamaterial

Graphene plasmonics, a vehicle for generating strong light and matter interactions, have been widely studied in optical fields such as ultrasensitive optical biosensing [1–3], photonic metamaterials [4], light harvesting [5, 6], optical nanoantennas [7] and quantum information processing [8]. They are an ideal substitute for noble-metal plasmons [9, 10]. In the mid-infrared and terahertz bands, the surface plasmon polaritons (SPPs) excited by graphene possess much tighter local field confinement effect [11]. Moreover, graphene SPPs not only can break through the traditional diffraction limit [12], but also achieve a relatively long propagation distance [13]. Importantly, by varying the Fermi level of graphene, the broadband of graphene SPPs exhibits the advantage of being highly tunable [14, 15]. Based on the above advantages, it has been favored by more and more researchers and widely explored in many applications, such as perfect absorbers [16–18], optical switches [19–24], optical storages [23–27] and optical biosensors [28].

Plasmon-induced transparency (PIT), a plasmonic analogue of classical electromagnetically induced transparency (EIT) [29], is generated by the destructive interference between the resonant modes excited in the SPPs metamaterials. It has been extensively studied in many reports, where the single-PIT can be realized in single-layer [18] or multi-layer structures [19, 22, 27] while the single-layer structure can also achieve dual-PIT [25, 30] and triple-PIT [21, 26]. In addition, the quadruple-PIT excited by five modes is rarely discussed. Further research indicates that the PIT based on graphene SPPs metamaterials can be dynamically regulated by a gate voltage [14] and can be used to design electro-optical switches and plasmonic logic gates [31–35]. Compared with the single-PIT, the multiple-PIT possesses more transparency windows and wider bandwidth, which can be used to realize the electro-optical switch with more modulated frequency points [23, 26]. To the best of our knowledge, there are few reports on multi-frequency switches.

In this paper, a monolayer graphene metamaterial composed of a graphene block and four graphene strips is proposed to generate an outstanding quadruple-PIT in terahertz frequency range. The chemical vapor deposition and the electron beam lithography can be employed to fabricate the patterned graphene of the proposed structure [36, 37]. For each PIT window, there are at least a bright mode and a dark mode present in the system, their strong destructive interference under the action of the incident light results in a PIT with four transparency windows. In addition, the distributions of electric field intensity and electric field vectors at the resonant frequencies further illustrate the forming physical mechanism of the quadruple-PIT. The spectral response characteristics of the proposed structure were investigated through finite-difference time-domain (FDTD) (Lumerical FDTD solutions) simulations and coupled mode theory (CMT) calculations, and a very satisfactory curve fitting effect was obtained. It is worth noting that six-frequency or five-frequency asynchronous switches can be achieved in the proposed structure. More interesting, a four-frequency synchronous switch has been achieved in the proposed structure, which is more prominent than the other graphene metamaterials with the same parameters [23, 26]. In particular, when E_f = 0.6 eV (E_f = 1.2 eV), the maximum modulation depth (MD) is up to 99.7%, which corresponds to the maximum extinction ratio (ER) of 25.04 dB and the minimum insertion loss (IL) of 0.19 dB. Finally, a significant performance of the slow-light effect in the proposed structure is obtained, the time delay is as high as 0.919 ps, and the corresponding group refractive index is up to 2755, which are superior to many other slow-light devices with the similar structures [23–27]. Thus, the proposed structure provides a guidance for the design of terahertz multi-frequency switches and slow light devices.

The designed three-dimensional schematic of the terahertz metamaterial structure is shown in figure 1(a). The monolayer graphene metamaterial is periodically distributed along the x- and y-directions at the x–y plane and deposited on a silicon substrate with the thickness of 100 nm. A linear polarized planar light is shined perpendicularly on the structure along the positive direction of the z axis to produce the excitation phenomenon. The thin metal wires are employed to connect the graphene strips between two adjacent graphene lattice units, the metal electrodes are added onto the graphene surface to modulate the Fermi levels of the graphene by changing the applied gate voltage V_g between the metal electrodes and the substrate, which is schematically shown in figures 1(a)–(c). The top view of the unit cell is schematically shown in figure 1(b). It can be seen that a graphene block and four graphene strips are deposited on the silicon substrate. The upper and lower graphene strips (ULGSs) are a horizontal double-rectangle graphene structure, the middle graphene block (MGB) is a vertical single-rectangle graphene structure, the left and right graphene strips (LRGSs) are a vertical double-rectangle graphene structure. The thickness of the monolayer graphene is 1 nm. The geometric parameters are as follows: l₁ = 5.4 μm, l₂ = 0.9 μm, l₃ = 0.1 μm, l₄ = 1.2 μm, l₅ = 3.6 μm, l₆ = 2 μm, l₇ = 0.6 μm, L_x = L_y = 6 μm.

Zoom In Zoom Out Reset image size

The Fermi level E_f can be altered by the applied gate voltage V_g according to the following equation [38, 39]:

$$\begin{equation}{E_{\text{f}}} = \hbar {V_{\text{F}}}{\left( {\frac{{\pi {\varepsilon _0}{\varepsilon _{{\text{si}}}}{V_{\text{g}}}}}{{{d_{{\text{sub}}}}e}}} \right)^{1/2}}\end{equation} \tag{ 1 }$$

where, d_sub is the thickness of the silicon substrate, ℏ is the reduced Planck constant. ε₀, ε_si (=11.9) [23], V_F and e are the vacuum permittivity, the relative permittivity of dielectric silicon, the Fermi velocity and the electron charge, respectively. According to the Pauli exclusion principle [40], the interband contribution of the complex conductivity of the monolayer graphene can be left out, owing to ℏω ≪ 2E_f in the terahertz range. After that the Fermi level satisfies the condition of E_f ≫ k_B T (k_B is Boltzmann constant), as the temperature T is set to 300 K. The complex conductivity of the monolayer graphene can be further simplified into a Drude-like model [41]:

$$\begin{equation}{\sigma _{\text{g}}}\left( \omega \right) = \frac{{i{e^2}{E_{\text{f}}}}}{{\pi {\hbar ^2}\left( {\omega + i{\tau ^{ - 1}}} \right)}}\end{equation} \tag{ 2 }$$

here, the carrier relaxation time τ can be expressed as $\tau = \mu {E_{\textrm{f}}}/\left( {eV_{\textrm{F}}^2} \right)$. In this paper, the Fermi velocity V_F and the carrier mobility μ of the whole monolayer graphene are fixed to 10⁶ m s⁻¹ and 4 m² (V s⁻¹), respectively [42]. Then, the effective permittivity of graphene can be calculated by the following expression [43]:

$$\begin{equation}\varepsilon \left( \omega \right) = {\varepsilon _{\text{r}}} + \frac{{i{\sigma _{\text{g}}}\left( \omega \right)}}{{{\varepsilon _0}\omega {t_{\text{g}}}}}\end{equation} \tag{ 3 }$$

where _r = 1 is the relative permittivity of the surrounding air medium, ω is the angular frequency of the incident light, t_g = 1 nm is the thickness of the monolayer graphene. According to the above equations (2) and (3), the graphene permittivity as the Fermi level changes in the interested frequency range can be presented in figure 2. It can be seen from figure 2(a) that the real part of the permittivity of the monolayer graphene increases rapidly with the increasing of the terahertz frequency from 0.5 to 3 THz. When the terahertz frequency exceeds 3 THz, the real part of the permittivity of the monolayer graphene increases gently. In addition, the real part of the permittivity of the monolayer graphene decreases with the increasing of the Fermi level under the constant frequency. The imaginary part of the permittivity of the monolayer graphene decreases radically firstly and then keeps a constant value with the increasing of the terahertz frequency. However, it nearly keeps unchanging under different Fermi levels, as shown in figure 2(b).

Zoom In Zoom Out Reset image size

Since the patterned graphene structure locates between the silicon substrate and air, the dispersion relationship between them can be clearly explained by utilizing the Maxwell equations and electromagnetic field boundary conditions, whose function expression can be as below [26, 44]:

$$\begin{equation}\frac{{{\varepsilon _{{\text{Si}}}}}}{{\sqrt {{\beta ^2} - {\varepsilon _{{\text{Si}}}}k_0^2} }} + \frac{{{\varepsilon _{{\text{air}}}}}}{{\sqrt {{\beta ^2} - {\varepsilon _{{\text{air}}}}k_0^2} }} = - \frac{{i\sigma }}{{\omega {\varepsilon _{\text{0}}}}}.\end{equation} \tag{ 4 }$$

β is the propagation constant of graphene surface plasmon wave, k₀ is the wave vector of the x-polarized light, ε_air is the relative permittivity of air.

The designed monolayer graphene structure comprised of LRGSs, MGB and ULGSs can be simulated numerically by FDTD method and achieve a significant quadruple-PIT phenomenon. The periodic boundary conditions are set along the x- and y-directions, and the perfectly matched boundary condition is applied in the z-direction. The effective simulation area is divided into Yee cells with the uniformly mesh step of 0.01 μm. The transmission spectra of the individual ULGSs (olive line), MGB (violet line) and the combined structure of ULGSs and MGB (red line) are plotted in figure 3(a). The transmission spectra of the individual LRGSs (blue line), the combined structure of ULGSs and MGB (red line), the whole structure (black line) are plotted in figure 3(b). It can be seen that the transmission spectrum of the combined structure of ULGSs and MGB (red line) in figures 3(a) and (b) has two induced transparent windows, so we name it as double-PIT, the transmission spectrum of the whole structure (black line) in figure 3(b) presents the quadruple-PIT.

Zoom In Zoom Out Reset image size

As shown in figure 3(a), we divide the frequency region into two parts with a black dashed line. In the first frequency region, it can be observed from the olive and violet lines that there is a significant transmission valley in the transmission spectrum of the ULGSs structure, while the spectral transmittance of the MGB structure is always higher than 90%. Therefore, the ULGSs act as a bright mode, while the MGB acts as a dark mode. The interaction between the bright mode and the dark mode produces two remarkable transmission valleys at the frequency points of $f_1^{\prime}$ = 2.97 THz and $f_2^{\prime}$ = 4.51 THz. The distributions of electric field intensity and electric field vectors in ULGSs, MGB and the combined structure of ULGSs and MGB at the frequency of $f_1^{\prime}$ = 2.97 THz are respectively displayed in figures 4(a)–(c). It can be seen from figures 4(a) and (b) that the electric field intensity is mainly focused on the two ends of ULGSs and there is nearly not electric field intensity in and around MGB, demonstrating further the above ideas that the ULGSs acting as a bright mode and the MGB acting as a dark mode are correct. From the distribution of the electric field intensity of the combined structure of ULGSs and MGB in figure 4(c), one can see that the electric field energy is mainly gathered at the two ends of ULGSs and the gaps between ULGSs and MGB, indicating that the dark mode acted by the MGB is excited by the local field generated by the coupling between the bright mode acted by the ULGSs and the incident light. In addition, at the frequency of $f_2^{\prime}$ = 4.51 THz, the distributions of electric field intensity and electric field vectors in ULGSs, MGB and the combined structure of ULGSs and MGB are shown in figures 4(d)–(f), respectively. One can see from figures 4(d) and (e) that there is nearly not electric field intensity in and around ULGSs and MGB. The reason for this phenomenon is that both transmittances of ULGSs and MGB at this frequency are beyond 90%, which can be seen from figure 3(a). The distributions of the electric field intensity in the combined structure of ULGSs and MGB in figure 4(f) are mainly focused on the gaps between ULGSs and MGB, and weak electric field intensity is distributed at the two ends of ULGSs, which lead to the resonant dip at the frequency of $f_2^{\prime}$ = 4.51 THz.

Zoom In Zoom Out Reset image size

In the second frequency region, both the transmission spectra of ULGSs and MGB show two apparent transmission valleys with different transmission depths. The violet line has a transmission valley closed to zero transmittance at the frequency range, but the olive line has a transmission valley with a transmittance of about 75% in this frequency region, so the MGB acts as a bright mode and the ULGSs act as a dark mode. The distributions of electric field intensity and electric field vectors of ULGSs, MGB and the combined structure of ULGSs and MGB at the frequency of $f_3^{\prime}$ = 6.48 THz are shown in figures 4(g)–(i), respectively. It can be seen from figure 4(g) that there is nearly not electric field energy in and around ULGSs, the reason is that the transmittance of ULGSs at the frequency of $f_3^{\prime}$ = 6.48 THz is close to 100% in figure 3(a), which indicates that the ULGSs acting as a dark mode is correct. What is more, one can see from figure 4(h) that the distribution of electric field intensity is mainly gathered at the four angles of MGB, demonstrating that the MGB acting as a bright mode is correct. The reason why the electric field intensity is weak is that the transmittance of MGB at the frequency of $f_3^{\prime}$ = 6.48 THz is up to about 70%, which can be seen from figure 3(a). Finally, it is easy to see from figure 4(i) that the electric field intensity of the combined structure of ULGSs and MGB is mainly focused on the gaps between ULGSs and MGB, indicating that the dark mode acted by the ULGSs is excited by the coupled light field between the incident light and the bright mode acted by the MGB. The interaction between the bright mode and the dark mode under the irradiation of the incident light produces the resonant dip at the frequency of $f_3^{\prime}$ = 6.48 THz.

As shown in figure 3(b), there are five resonant dips in the transmission spectrum of the whole structure (black line). We name them as dip1, dip2, dip3, dip4 and dip5 (from left to right), the corresponding frequencies are f₁ = 2.60 THz, f₂ = 4.37 THz, f₃ = 5.11 THz, f₄ = 5.62 THz and f₅ = 6.38 THz, respectively. The distributions of electric field intensity and electric field vectors at five resonant dips are plotted separately in figures 5(a)–(e). In order to better analyze the forming mechanism of the quadruple-PIT, we divide the frequency range into three frequency regions with two black dashed lines in figure 3(b). In the first frequency region, the blue spectral line is always closed to 1, so the LRGSs act as a dark mode, while the red spectral line produces a significant transmission valley closed to zero. Combined with figure 3(a), the ULGSs acts as a bright mode, the MGB and the LRGSs act as two dark modes, and the three structures together produce a resonant dip at the frequency point of f₁ = 2.60 THz under the action of the incident light. Furthermore, it can be observed from figure 5(a) that the electric field intensity is mainly distributed at the left and right edges of the ULGSs and the interface between the ULGSs and LRGSs, indicating that the resonant dip is mainly ascribed to the bright mode of ULGSs and the dark mode of LRGSs.

Zoom In Zoom Out Reset image size

Similar to the first frequency region, in the second frequency region, the red spectral line also produces a significant transmission valley, while the transmittance of the blue spectral line decreases with the frequency increasing and is still higher than that of the red spectral line, so the ULGSs still act as the bright mode, the MGB and LRGSs act as two dark modes. The interaction between the bright mode and two dark modes produces two resonant dips at the frequency points of f₂ = 4.37 THz and f₃ = 5.11 THz. The corresponding distributions of electric field intensity and electric field vectors at these two resonant frequencies are shown in figures 5(b) and (c), respectively. It can be observed from figure 5(b) that the electric field intensity is mainly gathered at the interface between the bright mode of ULGSs and the dark mode of MGB, illustrating that the formation of the resonant dip at the frequency of f₂ = 4.37 THz is mainly attributed to ULGSs and MGB. Similarly, it is easy to find from figure 5(c) that the formation of the resonant dip at the frequency of f₃ = 5.11 THz is mainly ascribed to the bright mode of ULGSs and the dark mode of LRGSs.

In the third frequency region, the blue spectral line generates a significant transmission valley with a transmission rate of 3.44%. The transmittance of the red spectral line is always higher than that of the blue spectral line throughout the frequency region. Combined the transmission spectra of ULGSs and MGB in figure 3(a), it can be seen that LRGSs and MGB act as two bright modes in this frequency region, while ULGSs act as a dark mode. The interaction between two bright modes and the dark mode produces two significant resonant dips at the frequency points of f₄ = 5.62 THz and f₅ = 6.38 THz, the corresponding distributions of electric field intensity and electric field vector are respectively displayed in figures 5(d) and (e). One can see from figure 5(d) that the electric field intensity is mainly distributed at the outer edges of the LRGSs, indicating that the resonant dip at the frequency of f₄ = 5.62 THz is mainly affected by the bright mode of the LRGSs. A small part of the electric field intensity is distributed at the four ends of the ULGSs, the reason is that the dark mode is excited by the coupled light field between the incident light and LRGSs. In addition, it can be seen from figure 5(e) that the electric field intensity is mainly gathered at the two edges of LRGSs and MGB, a fraction of electric field intensity is distributed at the gaps between ULGSs and MGB, illustrating that the dark mode of ULGSs is excited by the coupled light field between the incident light and MGB. The formation of the resonant dip at the frequency of f₅ = 6.38 THz is mainly attributed to the interaction between the two bright modes and a dark mode under the irradiation of the incident light.

Through the above analysis, it can be seen that ULGSs, MGB and LRGSs act as different modes at the three frequency regions, their interaction under the irradiation of the incident light forms the quadruple-PIT.

In order to further explore the physical mechanism of the quadruple-PIT phenomenon in the proposed graphene metamaterial structure, CMT is used to fit the transmission spectra produced by the FDTD method. The theoretical coupling diagram based on CMT is shown in figure 6, five excitation states of the hypothetical resonator modes denoted as A¹, A², A³, A⁴ and A⁵ are closely related to each other. $A_{ \pm ,{\text{in}}}^{\left( n \right)}$ and $A_{ \pm ,{\text{out}}}^{\left( n \right)}$ (n = 1, 2, 3, 4, 5) are represented as the nth mode with amplitude a_n , where the subscripts 'in' and 'out' respectively represent the incoming and outgoing streams of oscillatory resonator modes, the subscripts '+' and '−' respectively denote the positive and negative directions of the propagation streams. According to the basic CMT [45, 46], the time-varying equations among the five resonant modes can be obtained mathematically as below:

$$\begin{align}\frac{{{\text{d}}{a_1}}}{{{\text{d}}t}} & = \left( { - j{\omega _1} - \frac{1}{{{\tau _{{\text{o}}1}}}} - \frac{1}{{{\tau _{{\text{i}}1}}}}} \right){a_1} + A_{ + ,{\text{in}}}^{(1)}\sqrt {\frac{1}{{{\tau _{{\text{o}}1}}}}} \nonumber\\ & \quad + A_{ - ,{\text{in}}}^{(1)}\sqrt {\frac{1}{{{\tau _{{\text{o}}1}}}}} - j{\mu _{12}}{a_2} - j{\mu _{13}}{a_3} - j{\mu _{14}}{a_4} - j{\mu _{15}}{a_5},\end{align} \tag{ 5 }$$

$$\begin{align}\frac{{{\text{d}}{a_2}}}{{{\text{d}}t}} & = \left( { - j{\omega _2} - \frac{1}{{{\tau _{{\text{o}}2}}}} - \frac{1}{{{\tau _{{\text{i}}2}}}}} \right){a_2} + A_{ + ,{\text{in}}}^{(2)}\sqrt {\frac{1}{{{\tau _{{\text{o}}2}}}}} \nonumber\\ & \quad + A_{ - ,{\text{in}}}^{(2)}\sqrt {\frac{1}{{{\tau _{{\text{o}}2}}}}} - j{\mu _{21}}{a_1} - j{\mu _{23}}{a_3} - j{\mu _{24}}{a_4} - j{\mu _{25}}{a_5},\end{align} \tag{ 6 }$$

$$\begin{align}\frac{{{\text{d}}{a_3}}}{{{\text{d}}t}} & = \left( { - j{\omega _3} - \frac{1}{{{\tau _{{\text{o}}3}}}} - \frac{1}{{{\tau _{{\text{i}}3}}}}} \right){a_3} + A_{ + ,{\text{in}}}^{(3)}\sqrt {\frac{1}{{{\tau _{{\text{o}}3}}}}} \nonumber\\ & \quad + A_{ - ,{\text{in}}}^{(3)}\sqrt {\frac{1}{{{\tau _{{\text{o}}3}}}}} - j{\mu _{31}}{a_1} - j{\mu _{32}}{a_2} - j{\mu _{34}}{a_4} - j{\mu _{35}}{a_5},\end{align} \tag{ 7 }$$

$$\begin{align}\frac{{{\text{d}}{a_4}}}{{{\text{d}}t}} & = \left( { - j{\omega _4} - \frac{1}{{{\tau _{{\text{o}}4}}}} - \frac{1}{{{\tau _{{\text{i}}4}}}}} \right){a_4} + A_{ + ,{\text{in}}}^{(4)}\sqrt {\frac{1}{{{\tau _{{\text{o}}4}}}}} \nonumber\\ & \quad + A_{ - ,{\text{in}}}^{(4)}\sqrt {\frac{1}{{{\tau _{{\text{o}}4}}}}} - j{\mu _{41}}{a_1} - j{\mu _{42}}{a_2} - j{\mu _{43}}{a_3} - j{\mu _{45}}{a_5},\end{align} \tag{ 8 }$$

$$\begin{align}\frac{{{\text{d}}{a_5}}}{{{\text{d}}t}} & = \left( { - j{\omega _5} - \frac{1}{{{\tau _{{\text{o}}5}}}} - \frac{1}{{{\tau _{{\text{i}}5}}}}} \right){a_5} + A_{ + ,{\text{in}}}^{(5)}\sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} \nonumber\\ & \quad + A_{ - ,{\text{in}}}^{(5)}\sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} - j{\mu _{51}}{a_1} - j{\mu _{52}}{a_2} - j{\mu _{53}}{a_3} - j{\mu _{54}}{a_4}.\end{align} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

where, ω_n (n = 1, 2, 3, 4, 5) is the resonant angular frequency of the nth resonant mode, 1/τ_in = ω_n /(2Q_in) and 1/τ_on = ω_n /(2Q_on) are respectively the internal loss coefficient and the external loss coefficient of the nth resonant mode, while μ_mn (m, n = 1, 2, 3, 4, 5, m ≠ n) is the coupling coefficient between different resonant modes. Furthermore, the total quality factor Q_tn, the internal loss quality factor Q_in, and the external loss quality factor Q_on of the nth resonant mode can be satisfied with the relationship: 1/Q_tn = 1/Q_in + 1/Q_on (n = 1, 2, 3, 4, 5), where the total quality factor is expressed as Q_tn = f_n /Δf_n (f_n and Δf_n are the resonant frequency and the full width of half maximum (FWHM) of the nth resonant mode, respectively) and Q_in can be obtained by the expression: Q_in = Re(n_eff)/Im(n_eff) [47]. n_eff, given by n_eff = β/k₀, is the effective refractive index of graphene. Based on the energy conservation law, the relationship between the input and the output waves can be satisfied by:

$$\begin{align}A_{ - ,{\text{in}}}^{(n - 1)} & = A_{ - ,{\text{out}}}^{(n)} \cdot {{\text{e}}^{j{\phi _{n - 1}}}}(n = 2,{ }3,{ }4,{ }5);{ }A_{ + ,{\text{in}}}^{(n + 1)}\nonumber\\ & = A_{ + ,{\text{out}}}^{(n)} \cdot {{\text{e}}^{j{\phi _n}}}(n = 2,{ }3,{ }4),\end{align} \tag{ 10 }$$

$$\begin{equation}A_{ \pm ,{\text{out}}}^{(n)} = A_{ \pm ,{\text{in}}}^{(n)} - \sqrt {\frac{1}{{{\tau _{{\text{on}}}}}}} \cdot {a_n}(n = 1,{ }2,{ }3,{ }4,{ }5),{ }A_{ - ,{\text{in}}}^{\left( n \right)} = 0.\end{equation} \tag{ 11 }$$

ϕ_n is the phase difference between the (n+ 1)th and nth resonant modes. Since the five hypothetical resonators are located at the same plane, ϕ_n is taken to be zero. Thus, the transmission coefficient of the quadruple-PIT system can be expressed as:

$$\begin{align} t & = \frac{{A_{ + ,{\text{out}}}^{(5)}}}{{A_{ + ,{\text{in}}}^{(1)}}} = 1 + \frac{1}{{{\gamma _5}{\tau _{{\text{o}}5}}}} \nonumber\\ & \quad + \left( {\sqrt {\frac{1}{{{\tau _{{\text{o}}1}}}}} + \sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} \frac{{{x_{51}}}}{{{r_5}}}} \right){t_1} + \left( {\sqrt {\frac{1}{{{\tau _{{\text{o}}2}}}}} + \sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} \frac{{{x_{52}}}}{{{r_5}}}} \right){t_2} \nonumber\\ & \quad + \left( {\sqrt {\frac{1}{{{\tau _{{\text{o}}3}}}}} + \sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} \frac{{{x_{53}}}}{{{r_5}}}} \right){t_3} + \left( {\sqrt {\frac{1}{{{\tau _{{\text{o}}4}}}}} + \sqrt {\frac{1}{{{\tau _{{\text{o}}5}}}}} \frac{{{x_{54}}}}{{{r_5}}}} \right){t_4} \end{align} \tag{ 12 }$$

where,

$$\begin{equation}{t_1} = - {t_{\textrm{a}}}/{t_{\textrm{c}}},\end{equation} \tag{ 13 }$$

$$\begin{equation}{t_2} = - {t_{\textrm{b}}}/{t_{\textrm{c}}},\end{equation} \tag{ 14 }$$

$$\begin{equation}{t_3} = \frac{{LT + PS}}{{KP - LO}} + \frac{{IP - LM}}{{KP - LO}} \cdot \frac{{{t_{\text{a}}}}}{{{t_{\text{c}}}}} + \frac{{JP - LN}}{{KP - LO}} \cdot \frac{{{t_{\text{b}}}}}{{{t_{\text{c}}}}},\end{equation} \tag{ 15 }$$

$$\begin{align}{t_4} & = \frac{{O(LT + PS) - T(LO - KP)}}{{P(KP - LO)}} \nonumber\\ & \quad + \frac{{M(LO - KP) + O(IP - LM)}}{{P(KP - LO)}} \cdot \frac{{{t_{\text{a}}}}}{{{t_{\text{c}}}}} \nonumber\\ & \quad + \frac{{N(LO - KP) + O(JP - LN)}}{{P(KP - LO)}} \cdot \frac{{{t_{\text{b}}}}}{{{t_{\text{c}}}}},\end{align} \tag{ 16 }$$

$$\begin{align} {t_{\text{a}}} & = \left[ {\left( {HS + LR} \right)\left( {CH + DG} \right) - \left( {HQ - DR} \right)\left( {HK - GL} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {KP - LO} \right)\left( {FL + HJ} \right) - \left( {JP - LN} \right)\left( {HK - GL} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad + \left[ {\left( {LT + PS} \right)\left( {HK - GL} \right) - \left( {KP - LO} \right)\left( {HS + LR} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {DF + BH} \right)\left( {HK - GL} \right) + \left( {CH + DG} \right)\left( {FL + JH} \right)} \right], \nonumber\\[4pt] {t_{\text{b}}} & = \left[ {\left( {HQ - DR} \right)\left( {HK - GL} \right) - \left( {CH + DG} \right)\left( {HS + LR} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {KP - LO} \right)\left( {EL + HI} \right) - \left( {IP - LM} \right)\left( {HK - GL} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad - \left[ {\left( {HS + LR} \right)\left( {KP - LO} \right) - \left( {LT + PS} \right)\left( {HK - GL} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {AH - DE} \right)\left( {HK - GL} \right) - \left( {CH + DG} \right)\left( {EL + HI} \right)} \right], \nonumber\\[4pt] {t_{\text{c}}} & = \left[ {\left( {AH - DE} \right)\left( {HK - GL} \right) - \left( {CH + DG} \right)\left( {EL + HI} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {KP - LO} \right)\left( {FL + JH} \right) - \left( {JP - LN} \right)\left( {HK - GL} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad + \left[ {\left( {DF + BH} \right)\left( {HK - GL} \right) + \left( {CH + DG} \right)\left( {FL + JH} \right)} \right]{\text{ }} \nonumber\\[4pt] &\quad \times \left[ {\left( {KP - LO} \right)\left( {EL + HI} \right) - \left( {IP - LM} \right)\left( {HK - GL} \right)} \right], \end{align} \tag{ 17 }$$

$$\begin{align} A & = {\gamma _1}{x_{25}} + {x_{15}}{x_{21}},{\text{ }}B = {\gamma _2}{x_{15}} + {x_{12}}{x_{25}},{\text{ }}C = {x_{15}}{x_{23}} - {x_{13}}{x_{25}}, \nonumber\\ D & = {x_{15}}{x_{24}} - {x_{14}}{x_{25}},{\text{ }}E = {x_{25}}{x_{31}} - {x_{21}}{x_{35}},{\text{ }}F \! = \! {\gamma _2}{x_{35}} + {x_{25}}{x_{32}}, \nonumber\\ G & = {\gamma _3}{x_{25}} + {x_{23}}{x_{35}},{\text{ }}H = {x_{25}}{x_{34}} - {x_{24}}{x_{35}},{\text{ }}I = {x_{35}}{x_{41}} - {x_{31}}{x_{45}}, \nonumber\\ J & = {x_{35}}{x_{42}} - {x_{32}}{x_{45}},{\text{ }}K = {\gamma _3}{x_{45}} + {x_{35}}{x_{43}},{\text{ }}L = {\gamma _4}{x_{35}} + {x_{34}}{x_{45}}, \nonumber\\ M & = {\gamma _5}{x_{41}} + {x_{45}}{x_{51}},{\text{ }}N = {\gamma _5}{x_{42}} + {x_{45}}{x_{52}},{\text{ }}O = {\gamma _5}{x_{43}} + {x_{45}}{x_{53}}, \nonumber\\ P & = {\gamma _4}{\gamma _5} - {x_{45}}{x_{54}}, \end{align} \tag{ 18 }$$

$$\begin{align} Q & = \sqrt {1/{\tau _{{\text{o}}1}}} {x_{25}} - \sqrt {1/{\tau _{{\text{o}}2}}} {x_{15}},{\text{ }}R = \sqrt {1/{\tau _{{\text{o}}2}}} {x_{35}} - \sqrt {1/{\tau _{{\text{o}}3}}} {x_{25}}, \nonumber\\ S & = \sqrt {1/{\tau _{{\text{o}}3}}} {x_{45}} - \sqrt {1/{\tau _{{\text{o}}4}}} {x_{35}},{\text{ }}T = \sqrt {1/{\tau _{{\text{o}}4}}} {\gamma _5} + \sqrt {1/{\tau _{{\text{o}}5}}} {x_{45}}, \end{align} \tag{ 19 }$$

$$\begin{equation}{x_{mn}} = \sqrt {\frac{1}{{{\tau _{{\text{o}}m}}{\tau _{{\text{o}}n}}}}} + j{\mu _{mn}}(m,n = 1,{ }2,{ }3,{ }4,{ }5;{ }m \ne n).\end{equation} \tag{ 20 }$$

Thus, the theoretical transmission of the quadruple-PIT system can be obtained by T = |t|². Then, we can further have a more in-depth study of the spectral response of the whole system.

The transmission spectra of the proposed graphene metamaterial structure under different Fermi levels of the monolayer graphene are given in figure 7(a). The blue solid and red dashed lines are the results obtained by the FDTD and CMT methods, respectively. It can be seen that the numerical simulated results by the FDTD method are well consistent with the theoretical calculated results by the CMT method, indicating that our results obtained by the theoretical calculation is correct. In addition, the transmission spectra take a blue shift with the increasing of the Fermi levels. This behavior can be explained by the resonance condition. The wave vector of surface plasmons along the graphene nanostructure can be expressed as k_spp ∝ 1/L and k_spp = ℏω²/(2α₀ E_f c), where α₀ = e²/(ℏc) represents the fine nanostructure constant, L is the length of the graphene nanostructure, and f is the resonant frequency. Thus, the resonant frequency can be expressed as f ∝ (α₀ E_f/(2π²ℏcL))^1/2 [48], which is the mathematical analysis of the blue shift that occurs in the transmission spectra of the graphene metamaterial structure when the Fermi level of graphene increases. In addition, the fitting quality factors of Q_tn and Q_in used in the CMT calculation are depicted in figure 7(b). As the graphene Fermi level grows in increments of 0.1 eV, Q_tn basically remains the same level and Q_in steadily rises. Figure 7(c) illustrates the three-dimensional theoretical evolution of the transmission spectra when the Fermi level of graphene is continuously modulated from 0.6 to 1.2 eV. It is easy to see that the transmission spectra of the proposed graphene metamaterial structure have an obvious blue shift as the Fermi level increases continuously, which is in accord with the results in figure 7(a).

Zoom In Zoom Out Reset image size

Given that the graphene Fermi level can be dynamically adjusted by changing the gate voltage V_g, the design exhibits dynamically tunable spectral response characteristics. As shown in figure 8, we have selected the transmission spectra corresponding to the different Fermi levels and divided them into three pairs. As shown in figure 8(a), when E_f is set to be 0.6 and 1.2 eV, it can be used to develop an excellent hexa-frequency asynchronous switch at the frequency points of [f_a1, f_a2, f_a3, f_a4, f_a5, f_a6] = [2.82, 3.60, 4.60, 5.32, 5.63, 6.49] THz and a quatary-frequency synchronous switch at the frequency points of [f_a1, f_a3, f_a4, f_a6] = [2.82, 4.60, 5.32, 6.49] THz. When E_f is set to be 0.8 and 1.0 eV, as shown in figure 8(b), it can be used to achieve the first penta-frequency asynchronous switch at the frequency points of [f_b1, f_b2, f_b3, f_b4, f_b5] = [4.38, 4.77, 5.11, 5.31, 5.62] THz and the first triple-frequency synchronous switch at the frequency points of [f_b1, f_b3, f_b5] = [4.38, 5.11, 5.62] THz. And when E_f is set to be 0.9 and 1.1 eV, it can also be used to achieve the second penta-frequency asynchronous switch at the frequency points of [f_c1, f_c2, f_c3, f_c4, f_c5] = [4.52, 4.95, 5.23, 5.48, 5.63] THz and the second triple-frequency synchronous switch at the frequency points of [f_c1, f_c3, f_c5] = [4.52, 5.23, 5.63] THz, as shown in figure 8(c). In order to obtain excellent switching performance, the transmission peaks and transmission valleys are used as the 'on' and 'off' states of the switch, respectively, which correspond to '1' and '0' in binary in the digital circuit, respectively. Tr, on/off and binaries corresponding to the frequency point of each switch-on/off at different Fermi energy levels are listed in tables 1–3, respectively. It can be clearly observed from table 1 that we can obtain the binary array (1 0 1 1 0 1) at the frequency points of [2.82, 3.60, 4.60, 5.32, 5.63, 6.49] THz when E_f = 0.6 eV, the binary array (0 1 0 0 1 0) can be achieved at the corresponding frequency points when E_f = 1.2 eV. Table 2 shows that the binary array (1 0 1 0 1) can be obtained at the frequency points of [4.38, 4.77, 5.11, 5.31, 5.62] THz when E_f = 0.8 eV, the binary array (0 1 0 1 0) can be achieved at the corresponding frequency points when E_f = 1.0 eV. From table 3, we can observe that the binary array (1 0 1 0 1) can be obtained at the frequency points of [4.52, 4.95, 5.23, 5.48, 5.63] THz when E_f = 0.9 eV, the binary array (0 1 0 1 0) can be achieved at the corresponding frequency points when E_f = 1.1 eV.

Zoom In Zoom Out Reset image size

The performance of the switch can be assessed with several inter-played parameters, such as modulation speed, MD, IL, energy consumption, ER, and dephasing time [49, 50]. In this paper, we mainly analyze the MD, the IL, and the ER. The MD is used to describe the degree of amplitude variation, which can be defined (in %) as MD = |T_on − T_off |/T_on × 100% [51]. The ER and the IL, which are two critical parameters to characterize the performance of the switch, can be defined (in dB) as ER = 10 lg(T_on/T_off) and IL = −10 lg(T_on), respectively [50, 52]. Here, T_on and T_off refer to the transmission peaks and valleys, respectively, during a full modulation cycle. When E_f is set to be 0.6 and 1.2 eV, the maximum MD is high up to 99.7%, which corresponds to the maximum ER of 25.04 dB and the minimum IL of 0.19 dB. When E_f = 0.8 eV (E_f = 1.0 eV) and E_f = 0.9 eV (E_f = 1.1 eV), the MD respectively gets the maximum values of 98.3% and 97.8%, respectively corresponding to the maximum ERs of 17.64 and 16.58 dB, which are excellent for most applications.

The sensitivity characteristics of the switch are mainly influenced by the dephasing time, which can be defined (in ps) as T = 2ℏ/FWHM [53–55]. For the hexa-frequency asynchronous switch, the dephasing time of the switch at the six frequencies are 3.22, 7.01, 6.24, 11.75, 4.49 and 5.71 ps, respectively. For the first penta-frequency asynchronous switch, the dephasing time of the switch at the five frequencies are 6.24, 11.10, 11.10, 11.10 and 11.75 ps, respectively. For the second penta-frequency asynchronous switch, the dephasing time of the switch at the five frequencies are 6.24, 11.10, 11.75, 11.10 and 12.11 ps, respectively. The detailed calculation results of MD, ER, IL and DT are listed in the following three tables.

In addition, the PIT effect, which shows an excellent dispersion property, has another potential application in optical storage for the proposed system. The slow light, an inherent property of EIT and PIT, has ability to slow down the light velocity and enhance interaction time between light and matter [56], making more light information to be stored. The time delay τ_g and the group refractive index n_g, which are two crucial parameters to characterize the performance of the slow light effect, are related to the transparency windows of the proposed system. They can be expressed as follows [57, 58]:

$$\begin{equation}{\tau _{\text{g}}} = \frac{{{\text{d}}\phi }}{{{\text{d}}\omega }},\end{equation} \tag{ 21 }$$

$$\begin{equation}{n_{\text{g}}} = \frac{c}{{{\nu _{\text{g}}}}} \cdot = \frac{c}{l} \cdot {\tau _{\text{g}}},\end{equation} \tag{ 22 }$$

where v_g is the group velocity, l is the thickness of the silicon substrate, and ϕ is the phase shift of the amplitude transmittance, which is calculated from equation (12) with ϕ = arg(t).

Figures 9(a) and (b) show the evolution of the phase shift and time delay when the Fermi level is modulated from 0.8 to 1.2 eV, respectively. The corresponding transmission spectra are shown by the blue solid line. As shown in figure 9(a), the phase shift corresponding to the red solid line can be obtained by the mathematical calculation. On the one hand, for the same frequency region, the transmission spectra gradually take a blue shift with the Fermi level increasing, and the phase shift also produces a blue shift. On the other hand, for a fixed Fermi level, the phase shift oscillates up and down with the frequency increasing, and shows most dramatically at the frequency points of the transmission valleys. Furthermore, it is interesting to note that the wider the transmission valley, the larger the FWMH, and the larger the amplitude of the phase shift. When the Fermi energy level is fixed to 1.2 eV, the phase shift achieves a maximum value of 1.24 and a minimum value of −1.22 near the first transmission valley. As shown in figure 9(b), it clearly demonstrates the connection between the induced transparency window and the group time delay, a stronger dispersion phenomenon is exhibited at the transmission valley. In this paper, we only discuss the slow light effect. When the Fermi level is modulated from 0.8 to 1.2 eV, the maximum time delay near each transmission valley increases as the Fermi level increases. In other words, the slow light effect near each transmission valley also becomes more significant as the Fermi level increases. When the Fermi level is fixed at 1.2 eV, the maximum time delays near the five transmission valleys (from left to right) are 0.373, 0.605, 0.816, 0.919 and 0.667 ps, respectively. Meanwhile, the corresponding maximum group refractive indexes near the five transmission valleys are 1117.71, 1813.27, 2446.63, 2755.24, and 2000.35, respectively. Therefore, our proposed graphene structure is able to obtain a maximum time delay of 0.919 ps and a maximum group refractive index of 2755, which are superior to the slow light effect in the other graphene metamaterial structures with the similar parameters [23, 26].

Zoom In Zoom Out Reset image size

A comprehensive comparison on the properties of our proposed graphene metamaterial structure with that of other similar graphene-based structures [19–27] is given in table 4. Table 4 clearly illustrates the excellent properties of our proposed structure. On the one hand, the MD, ER and group time delay (or group refractive index) of the proposed structure are superior to some previous works. On the other hand, compared with other structures the proposed structure is simpler and has realized the four induced transparent windows in a single-layer graphene structure, making the realization of multi-frequency synchronous and asynchronous switch possible. Therefore, the proposed structure has a significant advantage over other structures in terms of overall performance.

In summary, we have proposed a monolayer graphene metamaterial composed of a graphene block and four graphene strips and achieved an outstanding quadruple-PIT. For each PIT window, there is at least a bright mode and a dark mode present in the system, so a strong destructive interference between bright mode and dark mode under the action of incident waves results in a PIT with four transparency windows. In addition, the analysis on the distributions of electric field intensity and electric field vectors makes this physical mechanism more fully plausible. The spectral response characteristics of the proposed structure were investigated through numerical simulations and theoretical calculations, and a very satisfactory curve fitting effect was obtained. It is worth noting that a hexa-frequency asynchronous switch (including a quatary-frequency synchronous switch) and two penta-frequency asynchronous switches (including two triple-frequency synchronous switches) have been achieved in the proposed structure. For all switches, the MD is generally greater than 91%, the ER (10.67 dB < ER < 25.04 dB) and the IL (0.19 dB < IL < 0.76 dB) are also prominent. In particular, when E_f = 0.6 eV (E_f = 1.2 eV), the maximum MD is high up to 99.7%, which corresponds to the maximum ER of 25.04 dB and the minimum IL of 0.19 dB. Finally, a significant performance of the slow light effect in the proposed structure is obtained. The time delay is as high as 0.919 ps, the corresponding group refractive index is up to 2755. Thus, the proposed structure provides a significant guidance for the design of terahertz multi-frequency switches and slow light devices.

This work was supported in part by National Natural Science Foundation of China (NSFC) (61605018, 61535015, 11904032), Xi'an Scientific and Technological Projects (2020KJRC0013), Hubei Natural Science Foundation of China (2016CFC767), the Doctoral Scientific Research Startup Foundation of Yangtze University (801080010128). The authors would like to thank Professor Yang Shen (School of Materials Science and Engineering, Sun Yat-sen University, China) for the help and discussion on numerical simulation by using of Lumerical FDTD solutions.

The data that support the findings of this study are available upon reasonable request from the authors.

The authors declare that there are no conflicts of interest related to this article.