Theoretical Analysis of Terahertz Dielectric–Loaded Graphene Waveguide

The waveguiding of terahertz surface plasmons by a GaAs strip-loaded graphene waveguide is investigated based on the effective-index method and the finite element method. Modal properties of the effective mode index, modal loss, and cut-off characteristics of higher order modes are investigated. By modulating the Fermi level, the modal properties of the fundamental mode could be adjusted. The accuracy of the effective-index method is verified by a comparison between the analytical results and numerical simulations. Besides the modal properties, the crosstalk between the adjacent waveguides, which determines the device integration density, is studied. The findings show that the effective-index method is highly valid for analyzing dielectric-loaded graphene plasmon waveguides in the terahertz region and may have potential applications in subwavelength tunable integrated photonic devices.

In the THz band, Huang et al. [49] proposed a graphene-coated nanowire with a drop-shaped cross section to realize low loss waveguiding with an ultra-strong mode confinement. Additionally, a graphene-coated elliptical nanowire was suggested for ultra-deep We first employ EIM to give an approximate analytical model for the DLGSPWG and then verify the results by using the FEM simulation. In the EIM, the 2D cross-section of the DLGSPWG is divided into two 1D waveguide structures (Figure 1b,c). By individually solving the guided modes of these two 1D waveguides [58], one can finally obtain the  We first employ EIM to give an approximate analytical model for the DLGSPWG and then verify the results by using the FEM simulation. In the EIM, the 2D cross-section of the DLGSPWG is divided into two 1D waveguide structures (Figure 1b,c). By individually solving the guided modes of these two 1D waveguides [58], one can finally obtain the guided modes of the DLGSPWG. The proposed DLGSPWG (Figure 1a) could be obtained by narrowing the width of the top GaAs dielectric layer (Figure 1b). Then, we can calculate the effective mode index n co of the guided modes in such an equivalent four-layer structure. All the dimensions are infinite along the x axis, and the upper air and substrate dielectric layers are semi-infinite along the y direction. Here, n co is irrelevant to the width of the dielectric strip and is equivalent to that of the effective mode index of the SP mode supported by the air-dielectric-graphene-dielectric structure (Figure 1b), which is given as [31] n co = ε 0 (ε 1 + ε 2 )ci/σ g , where ε 0 is the permittivity in air, σ g is the optical conductivity of graphene, and c is the speed of light. Then, the GaAs strip serves as the core of a three-layer dielectric waveguide shown in Figure 1c. The effective mode index n cl is equivalent to that of the SP mode supported by the graphene-air interface and is given by [31] n cl = ε 0 (ε 1 + 1)ci/σ g . Finally, the eigen-equation of the equivalent dielectric planar waveguide for the N-th order guided transverse electric (TE) mode (i.e., TM plasmon mode) is given as [56]: where k 1 = k 0 n 2 co − n 2 neff , k 2 = k 0 n 2 neff − n 2 cl , k 0 = 2π/λ 0 , λ 0 is the wavelength of incident light in free space. w is the width of the dielectric strip, and n eff is the effective mode index of the propagating modes in DLGSPWG. By numerically solving Equation (1), the effective mode index n eff of the N-th order mode could be obtained. The real part of n eff is related to the dispersion, and the imaginary part Im(n eff ) is related to the modal loss. The power propagation length L P is defined as L P = 1/2α, where α is the loss factor related to Im(n eff ) and given by α = k 0 Im(n eff ). Therefore, the power propagation length could be calculated by L P = λ 0 /[4πIm(n eff )]. The cut-off condition of the guided modes is k 2 = 0. Then, the cut-off wavelength of the N-th order guided mode can be calculated by [58,60]: Here, the impact of the strip height h (assumed to be infinite) has not been taken into account in the EIM, since graphene plasmons are tightly concentrated at the interface of the graphene layer. In the FEM calculation, the height of the dielectric strip is h = 2 µm, which is enough to guarantee the accuracy. Furthermore, enlarging h will be much closer to the situation in EIM and will lead to a better match between the EIM and FEM results, especially for the fundamental mode. We will show this subsequently.
In the calculation, the graphene layer is modelled as an electric field-induced surface current J = σ g E on the surface of the substrate, thus neglecting the thickness of the graphene layer. Within the random-phase approximation, the complex optical conductivity of graphene [61][62][63] consists of the interband and intraband contributions in the THz range, that is σ g = σ intra + σ inter , with: Here, we set τ = 0.5 ps (unless otherwise mentioned) for the electron relaxation time and T = 300 K for the temperature. E F is the Fermi energy level of the graphene, ω is the angular frequency of incident light,h is the reduced plank constant, k B is the Boltzmann's constant, and e is the charge of the electron. The FEM results are obtained by use of the wave optics module of COMSOL Multiphysics. The eigenvalue solver is used to find modes of the waveguide. The calculation domain is 2λ 0 × 2λ 0 , and a perfectly matched layer (PML) is applied around the geometry to avoid the influence of reflection. A convergence analysis is also conducted to ensure that the numerical boundaries and meshing do not interfere with the solutions.
In the derivation of Equation (1), we assume that graphene can support and propagate TM plasmon modes. Actually, graphene can support TM plasmon modes when Im(σ g ) > 0 [32]. To illustrate this, we plot σ g with respect to the frequency and E F in Figure 2. Obviously, within the frequency region (2-10 THz, i.e., corresponding wavelength region 30-150 µm) and Fermi energy level region considered here, the imaginary parts of σ g (Im(σ g ) = σ i ) are always above zero, thus indicating that TM plasmon modes could be excited.
Here, we set τ = 0.5 ps (unless otherwise mentioned) for the electron relaxation time and T = 300 K for the temperature. EF is the Fermi energy level of the graphene, ω is the angular frequency of incident light, ћ is the reduced plank constant, kB is the Boltzmann's constant, and e is the charge of the electron.
The FEM results are obtained by use of the wave optics module of COMSOL Multiphysics. The eigenvalue solver is used to find modes of the waveguide. The calculation domain is 2λ0 × 2λ0, and a perfectly matched layer (PML) is applied around the geometry to avoid the influence of reflection. A convergence analysis is also conducted to ensure that the numerical boundaries and meshing do not interfere with the solutions.
In the derivation of Equation (1), we assume that graphene can support and propagate TM plasmon modes. Actually, graphene can support TM plasmon modes when Im(σg) > 0 [32]. To illustrate this, we plot σg with respect to the frequency and EF in Figure  2. Obviously, within the frequency region (2-10 THz, i.e., corresponding wavelength region 30-150 μm) and Fermi energy level region considered here, the imaginary parts of σg (Im(σg) = σi) are always above zero, thus indicating that TM plasmon modes could be excited.

Results and Discussion
By analytically solving Equation (1), we first investigate the effective mode indices of the SP modes supported by the DLGSPWG, and then verify the analytical results by using

Results and Discussion
By analytically solving Equation (1), we first investigate the effective mode indices of the SP modes supported by the DLGSPWG, and then verify the analytical results by using the FEM simulation. In the calculation, the width and height of the dielectric strip are w = 4 µm and h = 2 µm. As seen from Figure 3, the EIM results (symbols) are in good agreement with the FEM (solid lines) results when the wavelength is smaller than the cut-off wavelength. The fundamental plasmon mode (N = 0, TM 0 ) has a larger real part of the effective mode index (Re(n eff )) than the higher order modes (N = 1, 2, 3), which implies a shorter wavelength (λ SP = λ 0 /Re(n eff )), as shown in Figure 3a. All the mode indices are located between n co and n cl . When increasing the wavelength from 30 to 150 µm, Re(n eff ) decreases monotonically and the higher order modes are gradually cut off. The power propagation length increases monotonically, which suggests a smaller modal loss at longer wavelengths ( Figure 3b).
As the wavelength is close to the cut-off wavelength, the EIM results show slight deviations from the FEM results. This is due to the fact that the electromagnetic fields in the corner regions have not been taken into account in the EIM. At shorter wavelengths, the modal fields are perfectly localized in the GaAs strip. Therefore, the EIM results are highly consistent with the FEM simulation results. Figure 4a shows the modal fields (|E y |) of the first four order modes at the wavelength of 60 µm (5 THz), and one can see that the fields are perfectly localized below the GaAs strip. However, at longer wavelengths, the modal fields are not perfectly localized, which leads to slight difference between these two methods.
wavelength. The fundamental plasmon mode (N = 0, TM0) has a larger real part of the effective mode index (Re(neff)) than the higher order modes (N = 1, 2, 3), which implies a shorter wavelength (λSP = λ0/Re(neff)), as shown in Figure 3a. All the mode indices are located between nco and ncl. When increasing the wavelength from 30 to 150 μm, Re(neff) decreases monotonically and the higher order modes are gradually cut off. The power propagation length increases monotonically, which suggests a smaller modal loss at longer wavelengths (Figure 3b). As the wavelength is close to the cut-off wavelength, the EIM results show slight deviations from the FEM results. This is due to the fact that the electromagnetic fields in the corner regions have not been taken into account in the EIM. At shorter wavelengths, the modal fields are perfectly localized in the GaAs strip. Therefore, the EIM results are highly consistent with the FEM simulation results. Figure 4a shows the modal fields (|Ey|) of the first four order modes at the wavelength of 60 μm (5 THz), and one can see that the fields are perfectly localized below the GaAs strip. However, at longer wavelengths, the modal fields are not perfectly localized, which leads to slight difference between these two methods.  The fundamental plasmon mode is cut-off-free and has a lower modal loss compared with higher order modes and could thus be used for long-range propagation. Here, the waveguide could also support higher order modes, and we next study the single and multi-mode operation regions by solving Equation (2). As shown in Figure 4b, at a certain wavelength the numbers of guided modes (K) increase with an increasing width of the GaAs strip. The critical values of w for the single-mode operation are 0.3073 μm and 7.6937 μm for λ0 = 30 and 150 μm, respectively. As w increases, the single-mode operation region tends towards a longer wavelength. At a certain width, the number K decreases with an increase in the wavelength. Finally, single mode propagation could be achieved at a very small w and long wavelength.
As depicted in Figure 3, the EIM results (symbols) are in very good agreement with the FEM (solid lines) results when the wavelength is smaller than 120 μm. Therefore, the wavelength ranges from 30 to 120 μm (i.e., 2.5-10 THz) in what follows.
In Section 2, it is mentioned that the impact of the strip height h (assumed to be infinite) has not been taken into account in the EIM, while in the FEM calculation the height of the dielectric strip is h = 2 μm. To make the structure used in the numerical simulation The fundamental plasmon mode is cut-off-free and has a lower modal loss compared with higher order modes and could thus be used for long-range propagation. Here, the waveguide could also support higher order modes, and we next study the single and multi-mode operation regions by solving Equation (2). As shown in Figure 4b, at a certain wavelength the numbers of guided modes (K) increase with an increasing width of the GaAs strip. The critical values of w for the single-mode operation are 0.3073 µm and 7.6937 µm for λ 0 = 30 and 150 µm, respectively. As w increases, the single-mode operation region tends towards a longer wavelength. At a certain width, the number K decreases with an increase in the wavelength. Finally, single mode propagation could be achieved at a very small w and long wavelength. Figure 3, the EIM results (symbols) are in very good agreement with the FEM (solid lines) results when the wavelength is smaller than 120 µm. Therefore, the wavelength ranges from 30 to 120 µm (i.e., 2.5-10 THz) in what follows.

As depicted in
In Section 2, it is mentioned that the impact of the strip height h (assumed to be infinite) has not been taken into account in the EIM, while in the FEM calculation the height of the dielectric strip is h = 2 µm. To make the structure used in the numerical simulation much closer to that used in EIM, we have studied the effective mode indices with respect to the wavelength at different strip heights h for the TM 0 and TM 1 modes, as shown in Figure 5. We find that the electromagnetic field in the corner regions is only partly responsible for the difference between the EIM and FEM results. By enlarging h, these two results are in better agreement with each other, even at longer wavelengths. As seen from Figure 5, the FEM results with a larger strip height (blue (h = 4 µm) and black lines (h = 8 µm)) have smaller deviations compared with the EIM results (symbols) for the TM 0 mode. For instance, when h = 2 µm and the wavelength is 120 µm, the relative deviation of Re(n eff ) between FEM and EIM is about 9.9%, while the relative deviation is only 2.7% for h = 8 µm. As for the imaginary parts of n eff , the relative deviation is always less than 6%. Interestingly, the strip height seems to have little effect on the effective mode indices for the TM 1 mode, and the EIM results are always in good agreement with the FEM results, even at different h values. Next, we study the group velocity, defined as g , of the propagated THz plasmon modes. As shown in Figure 6, the group velocity of the TM0 mode gradually decreases with an increasing frequency, which is due to the fact that eff eff is much larger at a higher frequency. As for the TM1 mode, the group velocity first rapidly decreases and then approaches that of TM0. One can see that the EIM results are consistent with the FEM results.  Next, we study the group velocity, defined as v g = ∂ω/∂(k 0 n eff ), of the propagated THz plasmon modes. As shown in Figure 6, the group velocity of the TM 0 mode gradually decreases with an increasing frequency, which is due to the fact that Re(n eff ) +ω∂Re(n eff )/∂ω (that is, c/v g ) is much larger at a higher frequency. As for the TM 1 mode, the group velocity first rapidly decreases and then approaches that of TM 0 . One can see that the EIM results are consistent with the FEM results.
One of the outstanding properties of graphene is the tunability of the surface conductivity by employing a DC bias [64], and this therefore provides us with a feasible method for modulating the effective mode indices of the guided modes. Figure 7 shows the effective mode indices of the fundamental mode in DLGSPWG with respect to the wavelength at different E F values. As we have mentioned above, one could get a better fitting between the FEM and EIM results by enlarging the strip height at longer wavelengths. Hence, h is set to be 8 µm. The carrier density in graphene could reach a value as high as 10 14 cm −2 [65], and the corresponding E F is about 1.17 eV. In a relatively recent work [66], the Fermi energy level reached a value as high as 1.77 eV. Here, E F varies from 0.6 eV to 1.2 eV. It is worth mentioning that we set τ = 0.5 ps in the above investigations. Actually, τ is related to E F and is given by τ = nE F /eV F 2 [67,68], where n = 10,000 cm 2 /(V·s) is the carrier mobility of graphene and V F = 10 6 m·s −1 . As a result, τ is 0.9 (1.2) ps for E F = 0.9 (1.2) eV. From Figure 7a, one can see that at a fixed wavelength, Re(n eff ) decreases as E F increases. For E F = 0.9 eV, the analytical results of Re(n eff ) range from 38.49 to 7.24. Clearly, the EIM results are in good agreement with the FEM results. In terms of the imaginary parts, shown in Figure 7b, when enlarging E F , Im(n eff ) decreases dramatically, which is due to the fact that the interband contribution of σ g is reduced by increasing the Fermi energy level. Thus, the propagation length could be massively enhanced. However, the EIM results show a slight discrepancy with the FEM results at longer wavelengths. This is because the electromagnetic fields in the corner regions have not been taken into account in the EIM; thus, the modal loss is slightly larger than for the FEM results. Next, we study the group velocity, defined as g , of the propagated THz plasmon modes. As shown in Figure 6, the group velocity of the TM0 mode gradually decreases with an increasing frequency, which is due to the fact that eff eff is much larger at a higher frequency. As for the TM1 mode, the group velocity first rapidly decreases and then approaches that of TM0. One can see that the EIM results are consistent with the FEM results. One of the outstanding properties of graphene is the tunability of the surface conductivity by employing a DC bias [64], and this therefore provides us with a feasible method for modulating the effective mode indices of the guided modes. Figure 7 shows the effective mode indices of the fundamental mode in DLGSPWG with respect to the wavelength at different EF values. As we have mentioned above, one could get a better fitting between the FEM and EIM results by enlarging the strip height at longer wavelengths. Hence, h is set to be 8 μm. The carrier density in graphene could reach a value as high as 10 14 cm −2 [65], and the corresponding EF is about 1.17 eV. In a relatively recent work [66], the Fermi energy level reached a value as high as 1.77 eV. Here, EF varies from 0.6 eV to 1.2 eV. It is worth mentioning that we set τ = 0.5 ps in the above investigations. Actually, τ is related to EF and is given by τ = nEF/eVF 2 [67,68], where n = 10,000 cm 2 /(V·s) is the carrier mobility of graphene and VF = 10 6 m·s −1 . As a result, τ is 0.9 (1.2) ps for EF = 0.9 (1.2)  Figure 7b, when enlarging EF, Im(neff) decreases dramatically, which is due to the fact that the interband contribution of σg is reduced by increasing the Fermi energy level. Thus, the propagation length could be massively enhanced. However, the EIM results show a slight discrepancy with the FEM results at longer wavelengths. This is because the electromagnetic fields in the corner regions have not been taken into account in the EIM; thus, the modal loss is slightly larger than for the FEM results. Besides the modal properties, crosstalk between the adjacent waveguides, which determines the device integration density, is also an important factor in subwavelength photonic integration. In order to investigate the crosstalk, we set a coupling system with two DLGSPWG structures separated by a distance of S, as shown in Figure 8a. To evaluate the performance of the coupling system in the integrated THz circuit, crosstalk is characterized by using a coupling length LC [69], which is calculated by LC = π/(|βs − βas|), based on the coupled mode theory [69], where βs = k0neff,s and βas = k0neff,as are the complex propagation constants of the symmetric and antisymmetric modes at 5 THz, respectively. Figure  8b shows the normalized coupling length LC/LP with respect to S. The dependence of the coupling length LC on S indicates that the crosstalk would be weaker with the increase of S. To ensure a very low crosstalk, LC/LP needs to be large enough (usually exceeding 10). This means that the energy of the mode decays to 1/e of its original value before it is coupled to the adjacent structure; thus, the crosstalk between adjacent structures is negligible. Besides the modal properties, crosstalk between the adjacent waveguides, which determines the device integration density, is also an important factor in subwavelength photonic integration. In order to investigate the crosstalk, we set a coupling system with two DL-GSPWG structures separated by a distance of S, as shown in Figure 8a. To evaluate the performance of the coupling system in the integrated THz circuit, crosstalk is characterized by using a coupling length L C [69], which is calculated by L C = π/(|β s − β as |), based on the coupled mode theory [69], where β s = k 0 n eff,s and β as = k 0 n eff,as are the complex propagation constants of the symmetric and antisymmetric modes at 5 THz, respectively. Figure 8b shows the normalized coupling length L C /L P with respect to S. The dependence of the coupling length L C on S indicates that the crosstalk would be weaker with the increase of S. To ensure a very low crosstalk, L C /L P needs to be large enough (usually exceeding 10). This means that the energy of the mode decays to 1/e of its original value before it is coupled to the adjacent structure; thus, the crosstalk between adjacent structures is negligible. Here, even at a very small separation distance of S = 0.1 µm, L C /L P is as high as 26.89 for E F = 0.6 eV; thus, the crosstalk can almost be ignored. However, for E F = 1.2 eV, L C /L P is about 10 for S = 0.6 µm. The low-mode crosstalk between the adjacent waveguides shows its potential application in ultra-compact and high-density integration.

Conclusions
We extend the concept of dielectric-loaded plasmon waveguide from the near-and mid-infrared band to the THz band. The SP modes of the DLGSPWG are studied by using the EIM and FEM. The tunability of the fundamental mode, effective mode index, propagation length, and the cut-off wavelength of higher order modes are investigated. The EIM results are in good agreement with the FEM results in the THz band. Our findings show that the electromagnetic field in the corner regions is only partly responsible for the difference between the EIM and FEM results. Enlarging the strip height will lead to a better fitting between the EIM and FEM results, especially for the fundamental mode. Besides this, we show that the DLGSPWG shows very small crosstalk between the adjacent structures, even at a very small separation distance of 0.1 μm. The DLGSPWG may have potential applications in tunable subwavelength terahertz photonic devices.

Conclusions
We extend the concept of dielectric-loaded plasmon waveguide from the near-and mid-infrared band to the THz band. The SP modes of the DLGSPWG are studied by using the EIM and FEM. The tunability of the fundamental mode, effective mode index, propagation length, and the cut-off wavelength of higher order modes are investigated. The EIM results are in good agreement with the FEM results in the THz band. Our findings show that the electromagnetic field in the corner regions is only partly responsible for the difference between the EIM and FEM results. Enlarging the strip height will lead to a better fitting between the EIM and FEM results, especially for the fundamental mode. Besides this, we show that the DLGSPWG shows very small crosstalk between the adjacent structures, even at a very small separation distance of 0.1 µm. The DLGSPWG may have potential applications in tunable subwavelength terahertz photonic devices.