Ultra-Low-Loss Mid-Infrared Plasmonic Waveguides Based on Multilayer Graphene Metamaterials

Manipulating optical signals in the mid-infrared (mid-IR) range is a highly desired task for applications in chemical sensing, thermal imaging, and subwavelength optical waveguiding. To guide highly confined mid-IR light in photonic chips, graphene-based plasmonics capable of breaking the optical diffraction limit offer a promising solution. However, the propagation lengths of these materials are, to date, limited to approximately 10 µm at the working frequency f = 20 THz. In this study, we proposed a waveguide structure consisting of multilayer graphene metamaterials (MLGMTs). The MLGMTs support the fundamental volume plasmon polariton mode by coupling plasmon polaritons at individual graphene sheets over a silicon nano-rib structure. Benefiting from the high conductivity of the MLGMTs, the guided mode shows ultralow loss compared with that of conventional graphene-based plasmonic waveguides at comparable mode sizes. The proposed design demonstrated propagation lengths of approximately 20 µm (four times the current limitations) at an extremely tight mode area of 10−6 A0, where A0 is the diffraction-limited mode area. The dependence of modal characteristics on geometry and material parameters are investigated in detail to identify optimal device performance. Moreover, fabrication imperfections are also addressed to evaluate the robustness of the proposed structure. Moreover, the crosstalk between two adjacent present waveguides is also investigated to demonstrate the high mode confinement to realize high-density on-chip devices. The present design offers a potential waveguiding approach for building tunable and large-area photonic integrated circuits.


Introduction
Significant enhancements achieved in light-matter interactions, nonlinear optical effects, chemical and biological sensing sensitivities, and resolution in imaging and spectroscopy can further benefit from the localization of light waves in deep subwavelength sizes. To date, the most promising approach to achieve this goal is to excite surface plasmon polaritons (SPPs) [1] by coupling photons and free electrons at metal-dielectric interfaces. For visible and near-infrared (IR) wavelengths, noble metals have been used to build various nanoscale photonic circuits [2][3][4][5]. However, alternative metallic materials are necessary because of the weak confinement of SPP modes in conventional noble metals for mid-IR and terahertz (THz) bands [6,7]. Graphene [8][9][10][11][12] is an emerging 2D material with extraordinary electric, thermal, and optical properties that can be flexibly tuned by electrical gating or chemical doping. Graphene is considered a promising candidate for SPP waveguiding in the mid-IR [6,7,9,11,12] and THz [9,12] ranges because of its nearly pure imaginary surface conductivity and extreme light confinement. The real and imaginary parts of the surface conductivity of graphene characterize the ohmic loss and magnitude of the wavevector, respectively. At the spectral bands of mid-IR and THz, the real part of graphene's surface conductivity approaches zero, leading to a comparatively low ohmic loss. Because of their superior merits, many graphene-based optoelectronic and photonic devices, including polarizers [13][14][15][16], modulators [17][18][19][20][21], sensors [22][23][24], switches [25][26][27][28], and couplers [29], have been reported in recent years. To design functional optical devices well, detailed analyses of the mode properties of graphene plasmonic waveguides (GPWs) are of essential importance.
Many studies [30][31][32][33][34][35][36][37][38][39] have reported a variety of GPWs operating in the mid-IR range to improve SPP mode performance. For a plasmonic waveguide, several indices to evaluate the waveguiding performance of a device include the normalized mode area (A m ); the propagation length (L p = λ/[4πIm(n e )]), where λ is the working wavelength in free space and Im(n e ) is the imaginary part of the effective refractive index; and the figure of merit (FoM), which indicates the ratio between L p and mode size. For a dielectric loaded GPW, Xu et al. [31] theoretically investigated the dependence of the L p of fundamental graphene-based SPP (GSPP) modes on the geometry parameter, frequency (f ), Fermi energy (E F ), and carrier mobility (µ) of graphene. At f = 20 THz, they obtained L p < 3 µm at E F = 0.6 eV and µ = 1 m 2 /V·s (the same value is used for subsequent comparisons, unless stated otherwise). Note that the calculations of L p in Refs. [31][32][33][34][35][36][37][38][39] are unified by the same definition as mentioned above. Liu et al. [32] proposed a symmetrical long-range GSPP hybrid waveguide on a silica (SiO 2 ) buffer layer on a silicon (Si) substrate showing an Lp of approximately 5 µm and an extremely confined area, A m = 8.0 × 10 −7 , at f = 30 THz and E F = 0.8 eV. Considering a pattern-free suspended graphene sheet over a Si ridge substrate, Bahadori-Haghighi et al. [33] numerically demonstrated a low-loss GPW with L p = 9 µm at f = 30 THz and E F = 0.35 eV.
In addition to these GPWs [31][32][33], researchers have focused on another kind of GPW, called graphene-coated nanowire waveguide (GCNW)-related structures [34][35][36][37][38][39]. Gao et al. [34,35] showed that a single GCNW performs low propagation loss and is cutoff free for the fundamental mode. At f = 30 THz, a single GCNW with a radius of R = 100 nm and permittivity ε NW = 2.1 shows L p = 4 µm and A m = 9.0 × 10 −4 at E F = 0.6 eV for the fundamental mode. To improve the weak confinement of a radically polarized mode in a single GCNW [34,35], Teng et al. [36] adopted a GCNW dimer to couple two GCNW modes. The mode properties were moderately improved to L p = 4 µm and A m = 9.0 × 10 −5 at E F = 0.6 eV and f = 30 THz. To significantly enhance the mode localization, Liu et al. [37] reported a two-layer dielectric GCNW composed of a Si core surrounded by a SiO 2 layer and an outermost graphene sheet, forming a conventional hybrid plasmonic waveguide (HPW) structure [40]. They obtained an A m of approximately 1.5 × 10 −5 and an L p of approximately 1.5 µm at E F = 0.5 eV operating at a wavelength of λ = 7 µm (about f = 42.85 THz). The cost for improving energy confinement is reflected in a moderately shorter L p . Later, Liu et al. [38] extended their design by adding an extra graphene sheet between the Si and SiO 2 layers to form a symmetric long-range coupling mode with an L p of~10 µm and A m of~10 −5 at E F = 0.6 eV and f = 30 THz. Last year, Teng et al. [39] adopted a GCNW-loaded Si nano-rib (GCNWLSNR) structure based on the coupling between the GCNW plasmon mode and the Si nano-rib to significantly shrink the mode size to A m = 9.8 × 10 −7 while keeping L p at~9 µm at E F = 0.6 eV and f = 20 THz. The experimental fabrication approach to a GCNW, by rolling a graphene ribbon, was presented in refs. [41,42] to demonstrate practical feasibility.
In this study, we propose a high-performance GPW based on multilayer graphene metamaterials (MLGMTs) [43,44] on a Si nano-rib waveguide partially covered by a lowindex porous SiO 2 film. The fundamental volume plasmon polariton (VPP) modes, which are supported by MLGMT-coupling SPPs at individual graphene sheets [45], are coupled with a dielectric mode supported by a Si nano-rib waveguide. The resulting coupled mode performs not only deep subwavelength mode confinement but also ultra-low loss through introducing the high-conductivity MLGMTs compared with those of previously published results [32][33][34][35][36][37][38][39]. The proposed approach provides an additional degree of freedom, the number of graphene layers, to control the mode characteristics. Moreover, fabrication imperfections and spectral response are also addressed to evaluate the robustness and operating bandwidth of the proposed structure.

Waveguide Structure and Methods
A 3D schematic diagram of the GPW and its front view are shown in Figure 1a,b, respectively. The structure consists of a porous SiO 2 [46,47] layer sandwiched by MLGMTs, which are formed by alternating graphene-dielectric layers, and a Si nano-rib waveguide deposited on a conventional SiO 2 substrate. Here, the same SiO 2 as the substrate is chosen as the dielectric layers in the MLGMTs. Following the HPW structure [40], the MLGMTs, porous SiO 2 , and Si nano-rib are considered the metal, low-index dielectric, and high-index dielectric layers, respectively. Therefore, the mode field of the present structure can be strongly squeezed in the nanoscale region between the Si nano-rib and the MLGMTs, thus achieving ultrasmall mode sizes. freedom, the number of graphene layers, to control the mode characteristics. Moreover, fabrication imperfections and spectral response are also addressed to evaluate the robustness and operating bandwidth of the proposed structure.

Waveguide Structure and Methods
A 3D schematic diagram of the GPW and its front view are shown in Figure 1a,b, respectively. The structure consists of a porous SiO2 [46,47] layer sandwiched by MLGMTs, which are formed by alternating graphene-dielectric layers, and a Si nano-rib waveguide deposited on a conventional SiO2 substrate. Here, the same SiO2 as the substrate is chosen as the dielectric layers in the MLGMTs. Following the HPW structure [40], the MLGMTs, porous SiO2, and Si nano-rib are considered the metal, low-index dielectric, and high-index dielectric layers, respectively. Therefore, the mode field of the present structure can be strongly squeezed in the nanoscale region between the Si nanorib and the MLGMTs, thus achieving ultrasmall mode sizes. The fabrication steps are schematically shown in Figure 2 and are described as follows: (1) Si and positive photoresist (PR) films with thickness hr (nm) are deposited on a conventional SiO2 substrate; (2) a mask is applied and followed by PR exposure and PR development to define a rectangular groove; (3) a Si layer is deposited and followed by the lift off of the PR to form a Si nano-rib waveguide; (4) the top of the Si nano-rib is rounded using e-beam lithography by carefully controlling the exposure time and  freedom, the number of graphene layers, to control the mode characteristics. Moreover, fabrication imperfections and spectral response are also addressed to evaluate the robustness and operating bandwidth of the proposed structure.

Waveguide Structure and Methods
A 3D schematic diagram of the GPW and its front view are shown in Figure 1a,b, respectively. The structure consists of a porous SiO2 [46,47] layer sandwiched by MLGMTs, which are formed by alternating graphene-dielectric layers, and a Si nano-rib waveguide deposited on a conventional SiO2 substrate. Here, the same SiO2 as the substrate is chosen as the dielectric layers in the MLGMTs. Following the HPW structure [40], the MLGMTs, porous SiO2, and Si nano-rib are considered the metal, low-index dielectric, and high-index dielectric layers, respectively. Therefore, the mode field of the present structure can be strongly squeezed in the nanoscale region between the Si nanorib and the MLGMTs, thus achieving ultrasmall mode sizes. The fabrication steps are schematically shown in Figure 2 and are described as follows: (1) Si and positive photoresist (PR) films with thickness hr (nm) are deposited on a conventional SiO2 substrate; (2) a mask is applied and followed by PR exposure and PR development to define a rectangular groove; (3) a Si layer is deposited and followed by the lift off of the PR to form a Si nano-rib waveguide; (4) the top of the Si nano-rib is rounded using e-beam lithography by carefully controlling the exposure time and  (1) Si and positive photoresist (PR) films with thickness h r (nm) are deposited on a conventional SiO 2 substrate; (2) a mask is applied and followed by PR exposure and PR development to define a rectangular groove; (3) a Si layer is deposited and followed by the lift off of the PR to form a Si nano-rib waveguide; (4) the top of the Si nano-rib is rounded using e-beam lithography by carefully controlling the exposure time and scanning speed; (5) a PR film is coated on the surface and is used to flatten the PR film using chemical mechanical polishing (CMP); (6) a mask is applied and followed by PR exposure and PR development to lift out the PR film and define the porous SiO 2 region; (7) a porous SiO 2 film is evaporated on the surface using the oblique deposition technique [46,47], the PR film is lifted off, and then CMP is used to flatten the porous SiO 2 layer; finally, (8) forming multilayer graphene-SiO 2 stacking by a chemical vapor deposition (CVD), a layer transfer method [48][49][50], or transfer-free, solution-phase deposition method [51,52]. Considering the high-technique step to form the MLGMTs, we describe the main fabrication processes as follows. The former approaches [48][49][50] include depositing graphene on a copper foil, transferring it to the substrate using the PMMA transfer technique, etching the copper foil using ammonium persulfate, doping the graphene by soaking the sample in a suitable solution, depositing the target dielectric layer by the atomic layer deposition being able to carefully control its thickness, and finally repeating the processes to achieve the desired number of layers. Note that the number of graphene layers in [50] can reach 11 layers. For the latter approach [51,52], the fabrication processes include the use of graphene-oxide (GO) in solution, deposition directly on a dielectric layer without requiring a transfer process, repetition to achieve the desired number of GO layers, and finally laser photoreduction from GO to graphene with removing the oxygen-functional groups, in which the bandgap of the graphene can be tuned by varying the lase power. Note that the number of graphene layers in [52] can reach 20 layers. Compared to the CVD transfer method [48][49][50], the advantages of the solution-phase method [51,52] are transfer-free, with quality independent of the number of layers, and controllable for the optical responses by tunning the bandgap of graphene. Consequently, the fabrication of the multilayer graphene metamaterials may be not that challenging and have complex techniques required and is moderately achievable using modern fabrication techniques, making the proposed waveguide structure practically realizable. Note that the conventional electrical gating on graphene layers uses a single voltage to a top contact [53,54], thus resulting in an inhomogeneous chemical potential of the graphene layers varying from layer to layer due to interlayer screening [55] in a multilayer graphene structure. The non-uniformity of chemical potential is more significant as the number of graphene layers increases. As a result, a potential scheme [49] can be adopted to achieve the required chemical potential in the proposed MLGHPW. This approach controls the Fermi energy levels of individual graphene layers by different gate voltages, making the carrier concentrations alter together in all layers.
In theory, the MLGMTs can be considered a coupled system with interacting multiple graphene sheets supporting multiple nondegenerate plasmon modes (called VPP modes). The VPP modes show hyperbolic isofrequency contours [45] and a large density of electromagnetic states, resulting in high-k guided modes. Among the VPP modes, the fundamental VPP mode (VPP 0 ) shows the lowest loss, although it possesses a comparably larger mode size. By coupling the VPP 0 with the dielectric mode of the Si nano-rib waveguide, a hybrid mode not only preserves the low-loss property of the VPP 0 , but also significantly benefits from the high mode confinement of the Si nano-rib. In the present structure, the geometry parameters are as follows: w g and t d are the width and thickness, respectively, of the dielectric (SiO 2 ) layer; N is the number of MLGMT layers between graphene sheets (N + 1 layers); and w r , h r , and t Si are the width, height, and bottom thickness of the Si nano-rib waveguide, respectively. Here, we set the radius of curvature to r = w r /2. The relative permittivities of Si, conventional SiO 2 , and porous SiO 2 at f = 20 THz are ε Si = 12.25 [56], ε SiO2 = 2.25 [56], and ε p-SiO 2 = 1.10 [46], respectively. Here, graphene is modeled as an infinitely thin sheet with a surface current density of J = σE in-plane, where E is the electric field vector and σ is the total surface conductivity of graphene with σ = σ intra + σ inter , including intraband (σ intra ) and interband (σ inter ) contributions, which can be calculated using the local random phase approximation [57]: where E F is the Fermi energy, τ = µE F /eV F 2 is the carrier relaxation lifetime, T is the temperature, k B is the Boltzmann constant,h is the reduced Planck constant, e is the electron charge, µ is the carrier mobility in graphene, and V F = 10 6 m/s is the Fermi velocity of electrons. Here, we consider a carrier mobility, µ = 1 m 2 /V·s, at T = 300 K, which is also used to study the performance of GPWs [34][35][36][37][38][39]. To evaluate the waveguiding performance of plasmonic waveguides, the propagation length (L p ), FoM (L p /2 √ A m /π [40]), and normalized mode area (A m = A e /A 0 ) are used, where A 0 = λ 2 /4 (λ is the working wavelength) is the diffraction-limited mode area and A e is the effective mode area given by Equation (3) includes the ratio of the total mode energy, W m , and the peak of the energy density, W(r), which is given by: where ω is the angular frequency, ε(r) is the profile of relative permittivity, µ 0 is the permeability in a vacuum, and |E| and |H| are the intensities of the electric and magnetic fields, respectively. The numerical results are calculated using the COMSOL Multiphysics software based on the rigorous finite element method. Figure 3 shows the mode properties of the present structure versus N for several values of t Si at the following parameters: w g = 200 nm, t d = 5 nm, w r = 10 nm, h r = 30 nm, E F = 0.6 eV, and f = 20 THz. We observe that the real part of n e , Re(n e ) sharply decreases from approximately 20 to 5 as N increases (see Figure 3a), but A m moderately increases from 5.34 × 10 −7 to 1.93 × 10 −6 for t Si = 5 nm (see Figure 3b). Note that the slopes of Re(n e ) and A m versus N become smaller when N is greater than five, meaning that Re(n e ) and A m are slightly influenced by larger N values. In contrast, L p linearly increases from 5.1 to 27.9 µm for t Si = 5 nm as N increases from 0 to 10 (see Figure 3c). Further increasing L p is only limited by the ability to fabricate an increasing number of N. The FoM increases from 6214 at N = 0 to 17,830 at N = 10 for t Si = 5 nm (see Figure 3d).

Waveguiding Performance Dependence on the Number of MLGMT Graphene Layers
The obtained results show that the proposed design can achieve both an ultralong L p (up to~28 µm) and an extremely small A m (on the order of 10 −6 ). Comparing these values with A m = 9.8 × 10 −7 and L p = 3.5 µm of the GCNWLSNR [39], the present structure improves L p by more than three times while maintaining a comparable order of magnitude of A m . For observing the mode profiles, Figure 4a-c depict the |E| values near the nanoscale region between the MLGMTs and Si nano-rib for N = 1, 5, and 10 at t Si = 5 nm. Clearly, increasing N results in a stronger peak of |E| but looser mode confinement, leading to the moderate increase in A m (see Figure 3b). The obtained results show that the proposed design can achieve both an ultralong Lp (up to ~28 µm) and an extremely small Am (on the order of 10 −6 ). Comparing these values with Am = 9.8 × 10 −7 and Lp = 3.5 µm of the GCNWLSNR [39], the present structure improves Lp by more than three times while maintaining a comparable order of magnitude of Am. For observing the mode profiles, Figure 4a-c depict the |E| values near the nanoscale region between the MLGMTs and Si nano-rib for N = 1, 5, and 10 at tSi = 5 nm. Clearly, increasing N results in a stronger peak of |E| but looser mode confinement, leading to the moderate increase in Am (see Figure 3b).   Figure 5a) and vertical dashed line V (inset of Figure 5b), respectively, to more clearly observe the relative mode profiles for different N values. The ultrasmall Am = 10 −6 of the present design shows full widths at half maximums of approximately 1 and 0.1 nm along the x and y directions, respectively. The Lp can be attributed to the mode profile being significantly shifted from the region of MLGMTs to the Si nano-rib as N increases (see Figure 5b). Note that the conditions EF = 0.6 eV, N = 10, and f = 20 THz are used in subsequent analyses unless stated otherwise.  The obtained results show that the proposed design can achieve both an ultralong Lp (up to ~28 µm) and an extremely small Am (on the order of 10 −6 ). Comparing these values with Am = 9.8 × 10 −7 and Lp = 3.5 µm of the GCNWLSNR [39], the present structure improves Lp by more than three times while maintaining a comparable order of magnitude of Am. For observing the mode profiles, Figure 4a-c depict the |E| values near the nanoscale region between the MLGMTs and Si nano-rib for N = 1, 5, and 10 at tSi = 5 nm. Clearly, increasing N results in a stronger peak of |E| but looser mode confinement, leading to the moderate increase in Am (see Figure 3b).      Figure 5b), respectively, to more clearly observe the relative mode profiles for different N values. The ultrasmall A m = 10 −6 of the present design shows full widths at half maximums of approximately 1 and 0.1 nm along the x and y directions, respectively. The L p can be attributed to the mode profile being significantly shifted from the region of MLGMTs to the Si nano-rib as N increases (see Figure 5b). Note that the conditions E F = 0.6 eV, N = 10, and f = 20 THz are used in subsequent analyses unless stated otherwise.
To clearly elucidate how the dielectric Si nano-rib changes the VPP 0 mode, we analyze the properties of the VPP 0 mode supported by the proposed structure without the nano-rib. Figure 6 shows the mode properties of the present structure with and without the nano-rib versus N at the same parameters as used in Figure 3 with t Si = 5 nm. To clearly elucidate how the dielectric Si nano-rib changes the VPP0 mode, we analyze the properties of the VPP0 mode supported by the proposed structure without the nano-rib. Figure 6 shows the mode properties of the present structure with and without the nano-rib versus N at the same parameters as used in Figure 3 with tSi = 5 nm. We observe that Re(ne) and Lp of the two structures are really close, but Am of the structure without the nano-rib is one to two orders of magnitude larger than that of the present structure as shown in Figure 6b, making the FoM of the structure without the nano-rib about one order of magnitude smaller. The above result reveals that the nano-rib can make the field distribution of the VPP0 mode significantly concentrated around the nanoscale region between the nano-rib and the MLGMTs. For observing the effect, Figure  7a,b show the |E| distributions of the present structure with and without the nano-rib, respectively, along with the zoomed-in view (see inset of Figure 7b) of the |E| distribution around the nano-rib. Without the nano-rib, the field profile spreads stronger out of the MLGMTs than that of the present design. In addition, the field of the present structure is significantly enhanced and focused mainly around the nano-rib, effectively shrinking its Am. To clearly elucidate how the dielectric Si nano-rib changes the VPP0 mode, we analyze the properties of the VPP0 mode supported by the proposed structure without the nano-rib. Figure 6 shows the mode properties of the present structure with and without the nano-rib versus N at the same parameters as used in Figure 3 with tSi = 5 nm. We observe that Re(ne) and Lp of the two structures are really close, but Am of the structure without the nano-rib is one to two orders of magnitude larger than that of the present structure as shown in Figure 6b, making the FoM of the structure without the nano-rib about one order of magnitude smaller. The above result reveals that the nano-rib can make the field distribution of the VPP0 mode significantly concentrated around the nanoscale region between the nano-rib and the MLGMTs. For observing the effect, Figure  7a,b show the |E| distributions of the present structure with and without the nano-rib, respectively, along with the zoomed-in view (see inset of Figure 7b) of the |E| distribution around the nano-rib. Without the nano-rib, the field profile spreads stronger out of the MLGMTs than that of the present design. In addition, the field of the present structure is significantly enhanced and focused mainly around the nano-rib, effectively shrinking its Am. We observe that Re(n e ) and L p of the two structures are really close, but A m of the structure without the nano-rib is one to two orders of magnitude larger than that of the present structure as shown in Figure 6b, making the FoM of the structure without the nano-rib about one order of magnitude smaller. The above result reveals that the nanorib can make the field distribution of the VPP 0 mode significantly concentrated around the nanoscale region between the nano-rib and the MLGMTs. For observing the effect, Figure 7a,b show the |E| distributions of the present structure with and without the nano-rib, respectively, along with the zoomed-in view (see inset of Figure 7b) of the |E| distribution around the nano-rib. Without the nano-rib, the field profile spreads stronger out of the MLGMTs than that of the present design. In addition, the field of the present structure is significantly enhanced and focused mainly around the nano-rib, effectively shrinking its A m .

Mode Characteristic Dependence on Geometric Parameters
To fully assess the waveguiding performance of the proposed structure, we analyze the geometrical dependence of the mode properties. Figure 8 shows the mode properties versus wr for several hr values at wg = 200 nm, tSi = 5 nm, and td = 5 nm. We observe that Re(ne) slightly depends on wr and hr (see Figure 8a) but Am shows a substantial dependence on wr and hr (see Figure 8b). For example, Am varies from 1.6 × 10 −6 to 6.0 × 10 −6 , while wr changes from 5 to 35 nm at hr = 30 nm. This effect can be attributed to a larger wr, obviously leading to a looser mode confinement. On the other hand, larger hr attains tighter Am due to the increase in the area of a low-index region, making the mode field more concentrated towards the region between the MLGMTs and the Si nano-rib. However, Lp significantly increases as hr increases, although Am shrinks. This is because the mode field shifts toward the Si nano-rib side as hr increases. Interestingly, Lp slightly increases for the smaller hr = 10 nm but moderately decreases for the larger hr > 30 nm as wr increases (see Figure 8c). For example, Lp varies from 29.2 (25.8) at wr = 5 nm to 26.8 (26.2) µm at wr = 35 nm for hr = 50 (10) nm. Note that smaller wr or larger hr (>20 nm here) can effectively improve the FoM of the proposed waveguide structure (see Figure 8d).  Width of Si nano-rib, w r (nm)

Mode Characteristic Dependence on Geometric Parameters
To fully assess the waveguiding performance of the proposed structure, we analyze the geometrical dependence of the mode properties. Figure 8 shows the mode properties versus w r for several h r values at w g = 200 nm, t Si = 5 nm, and t d = 5 nm. We observe that Re(n e ) slightly depends on w r and h r (see Figure 8a) but A m shows a substantial dependence on w r and h r (see Figure 8b). For example, A m varies from 1.6 × 10 −6 to 6.0 × 10 −6 , while w r changes from 5 to 35 nm at h r = 30 nm. This effect can be attributed to a larger w r , obviously leading to a looser mode confinement. On the other hand, larger h r attains tighter A m due to the increase in the area of a low-index region, making the mode field more concentrated towards the region between the MLGMTs and the Si nano-rib. However, L p significantly increases as h r increases, although A m shrinks. This is because the mode field shifts toward the Si nano-rib side as h r increases. Interestingly, L p slightly increases for the smaller h r = 10 nm but moderately decreases for the larger h r > 30 nm as w r increases (see Figure 8c). For example, L p varies from 29.2 (25.8) at w r = 5 nm to 26.8 (26.2) µm at w r = 35 nm for h r = 50 (10) nm. Note that smaller w r or larger h r (>20 nm here) can effectively improve the FoM of the proposed waveguide structure (see Figure 8d).

Mode Characteristic Dependence on Geometric Parameters
To fully assess the waveguiding performance of the proposed structure, we analyze the geometrical dependence of the mode properties. Figure 8 shows the mode properties versus wr for several hr values at wg = 200 nm, tSi = 5 nm, and td = 5 nm. We observe that Re(ne) slightly depends on wr and hr (see Figure 8a) but Am shows a substantial dependence on wr and hr (see Figure 8b). For example, Am varies from 1.6 × 10 −6 to 6.0 × 10 −6 , while wr changes from 5 to 35 nm at hr = 30 nm. This effect can be attributed to a larger wr, obviously leading to a looser mode confinement. On the other hand, larger hr attains tighter Am due to the increase in the area of a low-index region, making the mode field more concentrated towards the region between the MLGMTs and the Si nano-rib. However, Lp significantly increases as hr increases, although Am shrinks. This is because the mode field shifts toward the Si nano-rib side as hr increases. Interestingly, Lp slightly increases for the smaller hr = 10 nm but moderately decreases for the larger hr > 30 nm as wr increases (see Figure 8c). For example, Lp varies from 29.2 (25.8) at wr = 5 nm to 26.8 (26.2) µm at wr = 35 nm for hr = 50 (10) nm. Note that smaller wr or larger hr (>20 nm here) can effectively improve the FoM of the proposed waveguide structure (see Figure 8d).  Width of Si nano-rib, w r (nm) Next, we consider the effects of w g and t d of the MLGMTs on mode properties. At w r = 10 nm, h r = 30 nm, and t Si = 5 nm, Figure 9 shows the mode properties versus w g for several t d values. We observe that Re(n e ), A m , and L p moderately depend on w g and t d . We also observe that A m and L p decrease as w g decreases. Differing from increasing h r , which leads to a smaller A m (see Figure 8b), increasing w g increases A m , although they all increase the low-index region. This is because increasing w g also increases the width of the MLGMTs, leading to a looser mode field. The compensation between A m and L p makes the FoM almost constant for different values of w g . As for the t d , the smaller value leads to stronger coupling between SPP modes at individual graphene sheets. Therefore, smaller t d concurrently achieves lower loss and smaller mode size, breaking the trade-off between A m and L p to effectively improve the waveguiding performance of the proposed structure. Next, we consider the effects of wg and td of the MLGMTs on mode properties. At wr = 10 nm, hr = 30 nm, and tSi = 5 nm, Figure 9 shows the mode properties versus wg for several td values. We observe that Re(ne), Am, and Lp moderately depend on wg and td. We also observe that Am and Lp decrease as wg decreases. Differing from increasing hr, which leads to a smaller Am (see Figure 8b), increasing wg increases Am, although they all increase the low-index region. This is because increasing wg also increases the width of the MLGMTs, leading to a looser mode field. The compensation between Am and Lp makes the FoM almost constant for different values of wg. As for the td, the smaller value leads to stronger coupling between SPP modes at individual graphene sheets. Therefore, smaller td concurrently achieves lower loss and smaller mode size, breaking the trade-off between Am and Lp to effectively improve the waveguiding performance of the proposed structure.

Fabrication Tolerance, Material Parameters of Graphene, and Spectral Response
In experiments, fabrication imperfections lead to a reduction in waveguiding performance. Among the geometric parameters, the strictest part of the structure that should be precisely fabricated is the dimension of the Si nano-rib. Figure 10a Figure 10c,d, respectively. At r = 10 nm, Am moderately varies from 4.2 × 10 −6 to 8.5 × 10 −6 while Δr/r deviates from 0 to 0.5, but Lp remains constant, although Δr/r = 0.5. The analyzed results verify that the proposed structure possesses high fabrication tolerance on mode properties.

Fabrication Tolerance, Material Parameters of Graphene, and Spectral Response
In experiments, fabrication imperfections lead to a reduction in waveguiding performance. Among the geometric parameters, the strictest part of the structure that should be precisely fabricated is the dimension of the Si nano-rib. Figure 10a,b show the dependence of A m and L p on the relative fabrication error, ∆x/w r , at h r = 30 nm, t Si = 5 nm, t d = 5 nm, and w g = 200 nm, where ∆x is the fabrication error. Evidently, A m and L p are almost invariant in the range between ∆x/w r = 0 and 0.5. Similarly, the dependence of A m and L p on the relative fabrication error ∆r/r, where ∆r is the fabrication error, are shown in Figure 10c,d, respectively. At r = 10 nm, A m moderately varies from 4.2 × 10 −6 to 8.5 × 10 −6 while ∆r/r deviates from 0 to 0.5, but L p remains constant, although ∆r/r = 0.5. The analyzed results verify that the proposed structure possesses high fabrication tolerance on mode properties. Considering a high doping level of graphene that leads to reducing the carrier mobility (µ), we investigated the mode properties versus µ for several EF's, as shown in Figure 11. We observe that Re(ne) and Am are nearly independent on µ; however, Lp linearly reduces as µ decreases due to increasing the ohmic loss significantly for a higher doping level of EF > 0.5 eV. At EF = 0.4 eV, Lp varies from 10.7 to 7.4 µm at µ = 1 and 0.6 (m 2 /V·s), respectively. On the other hand, a lower EF attains a moderately higher Re(ne) and slightly smaller Am showing a tighter mode confinement, but significantly leads to a shorter Lp (see Figure 11c). For example, Lp varies from 19.5 to 10.8 µm at EF = 0.5 and 0.4 eV, respectively, for the condition µ = 1 (m 2 /V·s). The results reveal that a lower EF or µ leads to a reduction in Lp for general GPWs, mainly due to the higher ohmic losses. Considering a high doping level of graphene that leads to reducing the carrier mobility (µ), we investigated the mode properties versus µ for several E F 's, as shown in Figure 11. We observe that Re(n e ) and A m are nearly independent on µ; however, L p linearly reduces as µ decreases due to increasing the ohmic loss significantly for a higher doping level of E F > 0.5 eV. At E F = 0.4 eV, L p varies from 10.7 to 7.4 µm at µ = 1 and 0.6 (m 2 /V·s), respectively. On the other hand, a lower E F attains a moderately higher Re(n e ) and slightly smaller A m showing a tighter mode confinement, but significantly leads to a shorter L p (see Figure 11c). For example, L p varies from 19.5 to 10.8 µm at E F = 0.5 and 0.4 eV, respectively, for the condition µ = 1 (m 2 /V·s). The results reveal that a lower E F or µ leads to a reduction in L p for general GPWs, mainly due to the higher ohmic losses.
To fully study the tunability of graphene within a bandwidth range, we address the spectral response of mode properties. At w g = 200 nm, w r = 10 nm, h r = 30 nm, t Si = 5 nm, and t d = 5 nm, the results for several values of E F are shown in Figure 12. As the working frequency, f, increases from 10 to 30 THz, Re(n e ) and A m moderately increase, but L p significantly decreases. Considering the effect of E F , Re(n e ) increases, but A m and L p decrease because of the enhanced mode localization as E F decreases. Note that A m is slightly influenced while L p is significantly influenced by varying E F . At f = 10 (30) THz, our design achieves unprecedented long L p values of 12.1 (4.8) and 38.9 (19.4) µm while maintaining ultrasmall A m values of 6.8 × 10 −7 (1.8 × 10 −6 ) and 7.8 × 10 −7 (2.6 × 10 −6 ) for E F = 0.4 and 0.6 eV, respectively. Exceptionally, the obtained FoM values are higher than 10 4 (3.0 × 10 3 ) within f = 10-30 THz, even while operating at E F = 0.6 (0.3) eV.
reduces as µ decreases due to increasing the ohmic loss significantly for a higher doping level of EF > 0.5 eV. At EF = 0.4 eV, Lp varies from 10.7 to 7.4 µm at µ = 1 and 0.6 (m 2 /V·s), respectively. On the other hand, a lower EF attains a moderately higher Re(ne) and slightly smaller Am showing a tighter mode confinement, but significantly leads to a shorter Lp (see Figure 11c). For example, Lp varies from 19.5 to 10.8 µm at EF = 0.5 and 0.4 eV, respectively, for the condition µ = 1 (m 2 /V·s). The results reveal that a lower EF or µ leads to a reduction in Lp for general GPWs, mainly due to the higher ohmic losses. To fully study the tunability of graphene within a bandwidth range, we address the spectral response of mode properties. At wg = 200 nm, wr = 10 nm, hr = 30 nm, tSi = 5 nm, and td = 5 nm, the results for several values of EF are shown in Figure 12. As the working frequency, f, increases from 10 to 30 THz, Re(ne) and Am moderately increase, but Lp significantly decreases. Considering the effect of EF, Re(ne) increases, but Am and Lp decrease because of the enhanced mode localization as EF decreases. Note that Am is slightly influenced while Lp is significantly influenced by varying EF. The conventional electrical gating on graphene layers uses a single voltage to a top contact [50,51], thus resulting in an inhomogeneous chemical potential of the graphene layers varying from layer to layer due to interlayer screening [52] in a multilayer graphene structure. The non-uniformity of chemical potential is more significant as the number of graphene layers increases. As a result, a potential scheme [49] can be adopted to achieve the required chemical potential in the proposed MLGHPW. This approach controls the Fermi energy levels of individual graphene layers by different gate voltages, making the carrier concentrations alter together in all layers.

Waveguide Crosstalk
In addition to the size of A m , the crosstalk of the modes in adjacent waveguides complements to describe the degree of mode confinement and examines the feasibility for high integration of photonic integrated circuits. Figure 13 shows a coupled waveguide consisting of two parallel waveguides with a center-to-center separation s. significantly decreases. Considering the effect of EF, Re(ne) increases, but Am and Lp decrease because of the enhanced mode localization as EF decreases. Note that Am is slightly influenced while Lp is significantly influenced by varying EF. At f = 10 (30) THz, our design achieves unprecedented long Lp values of 12.1 (4.8) and 38.9 (19.4) µm while maintaining ultrasmall Am values of 6.8 × 10 −7 (1.8 × 10 −6 ) and 7.8 × 10 −7 (2.6 × 10 −6 ) for EF = 0.4 and 0.6 eV, respectively. Exceptionally, the obtained FoM values are higher than 10 4 (3.0 × 10 3 ) within f = 10-30 THz, even while operating at EF = 0.6 (0.3) eV.  Operating frequency, f (THz) The conventional electrical gating on graphene layers uses a single voltage to a top contact [50,51], thus resulting in an inhomogeneous chemical potential of the graphene layers varying from layer to layer due to interlayer screening [52] in a multilayer graphene structure. The non-uniformity of chemical potential is more significant as the number of graphene layers increases. As a result, a potential scheme [49] can be adopted to achieve the required chemical potential in the proposed MLGHPW. This approach controls the Fermi energy levels of individual graphene layers by different gate voltages, making the carrier concentrations alter together in all layers.

Waveguide Crosstalk
In addition to the size of Am, the crosstalk of the modes in adjacent waveguides complements to describe the degree of mode confinement and examines the feasibility for high integration of photonic integrated circuits. Figure 13 shows a coupled waveguide consisting of two parallel waveguides with a center-to-center separation s. According to the coupled mode theory (CMT) [58], energy exchange between the adjacent waveguides is due to the field coupling of evanescent tails of two normal modes. To evaluate the coupling strength, the coupling length of a lossless coupled waveguide system Lc = λ/[2(ns − na)] determines the length required to completely transfer power from one waveguide to another, where ns and na are the ne of the symmetric and asymmetric modes, respectively. For a lossy waveguide system such as plasmonic waveguides, we adopted a more suitable criterion, the normalized coupling length Lc/Lave [59], which considers both the power attenuation and maximum power transfer to measure the crosstalk, where Lave is the average Lp of the symmetric and asymmetric modes. The maximum transfer power ρmax between waveguides is only a function of Lc/Lave, and the adjacent waveguides can be considered as nearly isolated (ρmax = 0.33%) if Lc/Lave > 10 is reached. This is because the transfer power from one channel to the other is relatively weak at the distance of Lave. At the frequency f = 30 THz and the geometry parameters of wg = 180 nm, wr = 20 nm, hr = 30 nm, tSi = 25 nm, and td = 25 nm, the results of Lc/Lave along with ρmax versus s at µ = 1 m 2 /V.s for several EFs are shown in Figure 14a and those at EF = 0.4 eV are shown in Figure 14b for several µs. The results show that decreasing EF or µ leads to weaker coupling strength, and the dependence of coupling strength on EF is stronger than that on µ. We observe that the separations for negligible couplings between waveguides are s = 0.52, 0.64, and 0.72 µm for EF = 0.35, 0.40, and 0.45 eV, respectively (see s Figure 13. Schematic of a coupled waveguide consisting of two parallel MLGHPWs with a center-tocenter separation s. According to the coupled mode theory (CMT) [58], energy exchange between the adjacent waveguides is due to the field coupling of evanescent tails of two normal modes. To evaluate the coupling strength, the coupling length of a lossless coupled waveguide system L c = λ/[2(n s − n a )] determines the length required to completely transfer power from one waveguide to another, where n s and n a are the n e of the symmetric and asymmetric modes, respectively. For a lossy waveguide system such as plasmonic waveguides, we adopted a more suitable criterion, the normalized coupling length L c /L ave [59], which considers both the power attenuation and maximum power transfer to measure the crosstalk, where L ave is the average L p of the symmetric and asymmetric modes. The maximum transfer power ρ max between waveguides is only a function of L c /L ave , and the adjacent waveguides can be considered as nearly isolated (ρ max = 0.33%) if L c /L ave > 10 is reached. This is because the transfer power from one channel to the other is relatively weak at the distance of L ave .
At the frequency f = 30 THz and the geometry parameters of w g = 180 nm, w r = 20 nm, h r = 30 nm, t Si = 25 nm, and t d = 25 nm, the results of L c /L ave along with ρ max versus s at µ = 1 m 2 /V.s for several E F s are shown in Figure 14a and those at E F = 0.4 eV are shown in Figure 14b for several µs. The results show that decreasing E F or µ leads to weaker coupling strength, and the dependence of coupling strength on E F is stronger than that on µ. We observe that the separations for negligible couplings between waveguides are s = 0.52, 0.64, and 0.72 µm for E F = 0.35, 0.40, and 0.45 eV, respectively (see Figure 14a), and are s = 0.57, 0.61, and 0.64 µm for µ = 0.5, 0.75, and 1 m 2 /V.s, respectively (see Figure 14b). The results demonstrated that the coupling strength can be tuned by varying E F or µ, respectively. For instance, ρ max achieves about 20% at E F = 0.45 eV and 5% at E F = 0.35 eV for the condition µ = 1 m 2 /V.s and s = 0.4 µm. The small values of s (about 0.5 to 0.7 µm) make the proposed waveguide capable of building high-density graphene-based photonic-integrated circuits operating in the mid-infrared range.

Comparison of Waveguiding Performance
To demonstrate the superior waveguiding performance of the proposed design, we compared the reported results [32,34,36,38,39] in Table 1 with our mode properties at the following parameters: wg = 200 nm, wr = 10 nm, hr = 30 nm, tSi = 5 nm, td = 5 nm, EF = 0.6 eV, and N = 10. All the results in Table 1 were calculated at f = 30 THz, µ = 1 m 2 /V·s, and EF = 0.6 eV, except the results in ref. [32] with EF = 0.8 eV. If EF = 0.8 eV is decreased to 0.6 eV, Am and Lp will be further reduced. Liu et al. [32] achieved an extremely small area, Am = 8.0 × 10 −7 , by adopting a symmetrical HPW. Teng et al. [36] proposed a GCNW dimer to cause the coupling of the fundamental modes between two GCNWs to improve the Am of the single GCNW [34] by one order of magnitude while keeping the same Lp. Liu et al. [38] extended their previous report [32] to add an extra graphene sheet between Si and SiO2 layers and roll it into a cylindrical waveguide. It can be inferred that the two-layer graphene structure effectively improves the Lp but at the cost of a larger Am. Last year, Teng et al. [39] proposed a GCNWLSNR structure composed of a single GCNW deposited on a Si nano-rib, which showed a performance of Lp = 3.5 µm and Am = 2.0 × 10 −6 . In comparison with these published results [32,34,36,38,39], the proposed structure achieves an unprecedented waveguiding performance of Lp = 19.4 µm and Am = 2.6 × 10 −6 , thus obtaining an extremely high value of FoM = 10,612

Comparison of Waveguiding Performance
To demonstrate the superior waveguiding performance of the proposed design, we compared the reported results [32,34,36,38,39] in Table 1 with our mode properties at the following parameters: w g = 200 nm, w r = 10 nm, h r = 30 nm, t Si = 5 nm, t d = 5 nm, E F = 0.6 eV, and N = 10. All the results in Table 1 were calculated at f = 30 THz, µ = 1 m 2 /V·s, and E F = 0.6 eV, except the results in ref. [32] with E F = 0.8 eV. If E F = 0.8 eV is decreased to 0.6 eV, A m and L p will be further reduced. Liu et al. [32] achieved an extremely small area, A m = 8.0 × 10 −7 , by adopting a symmetrical HPW. Teng et al. [36] proposed a GCNW dimer to cause the coupling of the fundamental modes between two GCNWs to improve the A m of the single GCNW [34] by one order of magnitude while keeping the same L p . Liu et al. [38] extended their previous report [32] to add an extra graphene sheet between Si and SiO 2 layers and roll it into a cylindrical waveguide. It can be inferred that the two-layer graphene structure effectively improves the L p but at the cost of a larger A m . Last year, Teng et al. [39] proposed a GCNWLSNR structure composed of a single GCNW deposited on a Si nano-rib, which showed a performance of L p = 3.5 µm and A m = 2.0 × 10 −6 . In comparison with these published results [32,34,36,38,39], the proposed structure achieves an unprecedented waveguiding performance of L p = 19.4 µm and A m = 2.6 × 10 −6 , thus obtaining an extremely high value of FoM = 10,612. Table 1. Comparisons of the modal properties of A m , L p , and FoM.

Reference
A m L p (µm) FoM E F (eV) [32] 8.0 × 10 −7 5 4951 0.8 [34] 9.0 × 10 −4 4 118 0.6 [36] 9.0 × 10 −5 4 374 0.6 [38] 1.0 × 10 −5 10 2803 0.6 [39] 2.0 × 10 −6 3.5 2193 0.6 This work 2.6 × 10 −6 19.4 10,612 0.6 Currently, the obstacles for experimentally fabricating integrated graphene waveguides include three major points, despite the mature developments of both silicon photonics and graphene industries. (1) Excitation of extremely high-k SPP modes supported by monolayer graphene; (2) high-efficiency coupling between the nanoscale mode sizes of high-k SPP modes and submicron-scale dielectric waveguide modes or micron-scale optical fiber modes; (3) sufficient long propagation length beyond hundreds of micrometers or even millimeters. To relieve the first obstacle, employing multilayer graphene structures can efficiently control the field localizations and effective refractive index by varying the number of the graphene layers, thus overcoming the difficulty of exciting the high-k SPP modes. Next, the coupling efficiency can be improved by designing perfectly adiabatic metallic gratings, reducing field scattering during the coupling process. The final one is a common limit, as the SPP modes supported by general noble metals operating in the near-IR and visible light bands. In addition to discovering new low-loss materials, designing a novel waveguide structure by combining multiple guiding mechanisms together has been the most effective solution to decrease the propagation losses.

Conclusions
This work reported a mid-IR waveguiding structure based on MLGMTs on a Si nano-rib waveguide structure covered by a porous SiO 2 layer. By coupling the low-loss fundamental VPP mode of the MLGMTs, which is formed by coupling the SPP modes at individual graphene sheets and the dielectric mode of a Si nano-rib, the hybrid mode of the present design achieves an ultralong propagation length L p = 19.4 µm with A m = 2.6 × 10 −6 at E F = 0.6 eV operating at f = 30 THz. Compared with the reported results, the L p of our structure is five times greater than those reported at a comparable A m . Even for the looser A m , previously reported L p values at E F = 0.6 eV and f = 30 THz were still limited to below 10 µm. The MLGMTs provide a high-conductivity graphene structure that significantly increases L p with an increasing number of graphene layers. Therefore, the increased degree of L p is mainly restricted by the modern fabrication technique. In addition, the crosstalk between two adjacent waveguides demonstrates that the proposed structure is beneficial in realizing high-integration photonic devices operating in the mid-IR band. Our design is expected to pave the way for potential applications in building ultralow loss and compact and tunable mid-IR photonic devices and can be extended to other extraordinary 2D materials.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.