Open Access
Add to BibSonomy
Abstract
Guided modes of two-dimensional (2D) material-based plasmonic waveguides are applied in photonic devices because of their strong light–matter interaction within atomically thin layers and unique optical characteristics. Numerical simulations and experiments both play crucial roles in exploring unexpected phenomena at the sub-nanoscale of these materials. To efficiently and precisely compute mode characteristics, a multi-domain pseudospectral method (MPM) exhibiting high accuracy and fast convergence is proposed to study 2D material-based plasmonic waveguides in this study to alleviate the highly computational load of the widely used finite difference time domain or finite element method, as they demand extremely fine grid points or meshes around 2D materials. Models of graphene- and black phosphorus-based waveguides demonstrate that the MPM preserves exponential accuracy at relatively low computational cost, compared with the analytical characteristic equation and FEM, respectively. We believe that the proposed MPM offers a highly efficient and accurate approach to the study of 2D material-based photonics devices.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Two-dimensional (2-D) materials exhibit drastically different phenomena from their three- dimensional (3-D) counterparts. Recently, they have been emerging as a promising platform to build nanoscale devices [1]. Besides their extraordinary physical properties, 2-D materials including graphene, transition metal dichalcogenides (TMDs), hexagonal boron nitride, and black phosphorus (BP) are tuned by simple gating or doping to adapt to applications that are more extensive. Two-dimensional materials support various highly confined surface polariton modes in sub-nanoscale films. Among the members of the 2-D-material family, graphene [2], a monolayer of carbon atoms packed into a honeycomb crystal lattice, exhibits potential for application in the design of a variety of optoelectronic and photonic devices [3,4], which fuels extensive research on other 2-D materials. Particularly, rather than using noble metals to produce surface plasmon polaritons (SPPs) in the near-infrared (IR) wavelength regime [5,6], graphene with intrinsic tunable plasmons inspired the research on graphene-based plasmonics [7–9], enabling the manufacture of novel optical devices operating in the frequency regime ranging from terahertz to middle-IR (mid-IR). In graphene plasmonics, light is localized in the vicinity of an atomic thin sheet, leading to a significant enhancement of optical fields and thus offering a novel platform to manipulate light in a 2-D manner. Consequently, besides experimental studies, an accurate and efficient numerical approach is practically important to precisely investigate the localized mode characteristics based on 2-D materials.
Presently, the commonly applied numerical methods to study graphene plasmonics are the finite difference time domain (FDTD) using Lumerical FDTD solutions [10–12] and the finite element method (FEM) using COMSOL Multiphysics [13–17]. The FDTD is a versatile numerical technique used for broadband modeling by directly discretizing Maxwell’s curl equations. The desired response at any frequency is obtained by executing a single FDTD simulation with a wideband temporal source. Thus, FDTD has superior performance in the search for resonant frequencies that are not known in advance or in the prediction of broadband responses. Because FDTD discretizes space derivatives using second-order finite difference approximations, the grid size for the smallest geometrical feature (i.e., graphene sheet) must be extremely fine to guarantee numerical accuracy. Furthermore, a very small time-step is required to satisfy the Courant–Friedrich–Levy stability condition due to the use of the fine grid, thus significantly increasing computational costs [18]. The FEM is a prominent and versatile numerical scheme used to solve boundary value problems with complex geometry due to the use of arbitrary meshes to discretize the space domain. Inevitably, the atomically thin graphene sheets must be finely meshed in the FEM, as well as in the FDTD, because of the use of low-order linear or quadratic elements to capture the drastically varying fields. To alleviate the massive computational loads of the FDTD and FEM, an equivalent surface boundary condition [19,20], converting thin layers into conducting lines with complex conductivity, is assumed to model the graphene sheet. However, particular attention must be paid to discontinuous tangential magnetic components across the conducting lines in the FDTD and FEM, both of which benefit from avoiding volumetric discretization.
To overcome the extremely fine grid requirement due to the use of low-order approximations in FDTD and FEM, a multi-domain pseudospectral method (MPM) with high efficiency and fast convergence is applied to analyze the 2-D material-based plasmonic waveguides. The MPM was originally developed to address fluid dynamics problems; however, it was limited to its single-domain version [21,22]. The MPM was previously proposed by our group to successfully deal with waveguide problems [23,24], and the results demonstrate that the MPM outperforms FDTD and FEM regarding both accuracy and convergence rate, because of the use of well-matched high-order basis functions to represent unknown field profiles. In recent years, several groups likewise devoted to developing the variations of the MPM to solve various waveguide problems, including dielectric waveguides [25–28], diffraction gratings [29,30], photonic crystal waveguides [31], and leaky waveguides [32]. Although the MPM has not gained the extensive popularity of the FDTD and FEM, it complements the deficiency of the FDTD and FEM, particularly in addressing problems of high accuracy and lower computational loads. In this study, we aim to deal with atomically thin 2-D materials with the smallest geometrical feature at the sub-nanometer scale that is smaller than by a factor of 10−3 of the largest one. To our best knowledge, no reports to date have applied the MPM to analyze problems with significantly high-contrast geometry features. Examples including graphene- and BP-based plasmonic waveguides are analyzed thoroughly to demonstrate the suitability of the present scheme. We believe that the MPM offers a good approach to study a variety of 2-D material-based optoelectronic and photonic devices, provided that its’ superior computational performances are preserved.
2. Formulations of the wave equations and interfacial conditions
For a monochromatic light wave with a time dependence of exp(jωt) propagating along the z-direction in a medium with permittivity tensor [ε], the wave equation based on the magnetic field vector H is given by
To preserve the exponential convergence of the MPM, Eq. (2) is only used in each subdomain with homogeneous refractive index excluding the interface positions between different materials. The MPM then assembles all subdomains by considering the physical interfacial conditions including divergence-free of magnetic field vector $\nabla \cdot {\textbf{H}} = 0$ and Ampere’s law $\nabla \times $H = jω[ε]E, where E is the electric field vector. Note that explicitly imposing the constraint of $\nabla \cdot {\textbf{H}} = 0$ in the present formulations prevents the spurious modes from appearing in the solutions [33]. Using the continuous conditions of Hz and Ez at the interfaces between materials 1 and 2, the physical interfacial conditions are given by
Here, only TM modes are considered due to their better confinement for the design of compact 2D material-based optoelectronic and photonic devices comparing to that of TE counterpart [34].
3. Implementation of the MPM
First, the computational domain is divided into several subdomains with a uniform refractive index. This avoids discontinuous unknown fields by applying continuous basis functions in a single domain, giving rise to the Gibbs phenomenon [21,22], which preserves the high accuracy and exponential convergence of the MPM. In each subdomain, Hx and Hy are expanded by a product of suitable Lagrange-type interpolating basis functions Θ(x) with mx + 1 terms in the x-direction and Φ(y) with my + 1 terms in the y-direction multiplying the corresponding field components denoted as Hx,rs and Hy,rs at the (mx + 1)${\times} $(my + 1) collocation points. These collocation points are located at different positions for the chosen basis functions as follows [21–23]:
The pattern of the global matrix is a block diagonal matrix composed of the contributions from all subdomains. Finally, the interfacial points shared by adjacent subdomains must be replaced by the physical interface conditions of Eqs. (4a) and (4b) to connect the communications of unknown fields between subdomains. Subsequently, the degree of sparsity of the final matrix, which affects the computational time, to be solved here is below 40%, and the effective Arnoldi iteration method [35] is applied to solve the presented eigenvalue problems. The remainder issue of the present scheme is the choice of suitable basis functions. Note that the subdomains are classified into the interior with finite interval and exterior subdomains with semi-infinite extension. For an interior subdomain, the unknown fields are expanded by the Chebyshev polynomials (CPs), which are robustness for nonperiodic problems. By contrast, the unknown guided fields in exterior subdomains with exponential decay extending to infinity are expanded by the Laguerre-Gaussian functions (LGFs) showing exponential decay behavior with adjustable decay rate. Therefore, the LGFs can effectively model the profiles of evanescent waves of highly confined SPP modes for various 2-D materials. The best beneficial advantage of using the LGFs is that no absorbing boundary conditions (for example, perfectly matched layers) to enclose the computational boundaries are needed. The explicit form of the Lagrange-type CPs [22] is given by
4. Numerical results
In this section, 2D material-based plasmonic waveguides are investigated to demonstrate the high accuracy and exponential convergence of the MPM. First, we solve the symmetric and asymmetric SPP modes of a BP waveguide coupler for the light propagating along the zigzag (ZZ) and armchair (AC) crystal directions of BP. The results obtained in the example are compared with exact analytical solutions to evaluate the convergence of the present MPM. Furthermore, we verify the computational cost of the MPM by solving a 3D dielectric-loaded graphene plasmonic waveguide (DLGPW), compared with that of the FEM using the commercial COMSOL Multiphysics. Furthermore, the modal spectrum responses of the DLGPW are analyzed by controlling the Fermi energy levels, carrier mobilities, and the geometrical sizes of the dielectric ridge.
4.1 Black Phosphorus Slab Waveguide Coupler
A monolayer BP, also called phosphorene [36], is an elemental crystalline semiconductor comprising only phosphorus atoms in which each phosphorus atom covalently bonds with three adjacent phosphorus atoms forming a puckered hexagonal crystal structure as shown in Fig. 1(a). Layered BP exhibits high carrier mobility comparable to graphene and high on-off ratios that graphene with zero bandgap is unable to achieve, and it has great potential for mid-IR nanophotonic devices due to its highly in-plane anisotropic crystal structure with zigzag (ZZ) and armchair (AC) principal axes. In the example, we consider a phosphorene-based waveguide coupler consisted of two monolayer BP sheets with refractive index of nm separated by a gap g, and covered by substrate and cladding regions with refractive indices of ns and nc, respectively as shown in Fig. 1(b).
Fig. 1. Schematic diagrams of (a) the crystal structure of a monolayer BP and (b) a five-layer phosphorene-based waveguide coupler for the light propagating along the ZZ (z-) directions.
Here, a monolayer BP can be modeled as a uniaxial anisotropic permittivity by assuming that it has a finite thickness of tb = 0.53 nm, in which the equivalent permittivities of out-of-plane and in-plane permittivities can be expressed as ε┴ = 5.76 and ε‖ = 5.76 – jσin/(ωε0tb) [37], respectively, where σin is the in-plane conductivity. The σin shows highly in-plane anisotropy along the ZZ and AC crystal directions given by the formulas [37,38]:
In this example, we choose the relative permittivities nm = ns = nc = 1, the gap between BP sheets g = 50 nm, and the operating mid-IR wavelength λ = 30 µm. The permittivity tensor of the layered BP is [ε] = ε0 diag (εZZ [εAC], 5.76, εZZ [εAC]) for the ZZ (AC) direction aligning to the propagating (z-) direction of light. By the MPM, the computational domain is divided into five subdomains, in which three interior subdomains (the middle gap and the two BP sheets) are expanded by CPs and the other two exterior subdomains (substrate and cladding regions) are expanded by LGFs. The computational domains chosen are 1.2 and 2.5 µm for computing the symmetric and asymmetric SPP modes, respectively. Comparing with the analytical results [39] ($n_e^{\textrm{AN}}$ = 119.727923456506509 – j19.644659384139020 for the symmetric mode and $n_e^{\textrm{AN}}$ = 53.409161786092627 – j16.984227605081103 for the asymmetric mode) for the ZZ direction aligning to the z-direction, the relative errors of ne for the symmetric and asymmetric modes are shown in Figs. 2(a) and 2(b), respectively.
Fig. 2. Relative errors of the Re(ne) and Im(ne) of the MPM comparing with the analytical solutions [39] for the (a) symmetric and (b) asymmetric SPP modes of a BP-based plasmonic waveguide, while the ZZ direction aligning to the direction of light propagating.
We observe that the MPM shows the exponential convergence for both the symmetric and asymmetric SPP modes of the BP-based waveguide coupler. As N increases from 5 to 9 (only 41 [9 × 5–4] unknowns in total, where 5 and –4 denote five subdomains and 3 interfaces, respectively), the relative error of the Re(ne) (Im[ne]), exponentially decreases from the order of 10−2 (10−2) to 10−13(10−12) for the symmetric and asymmetric modes. Differing significantly from the FDTD and FEM, no extremely fine grid is required to discretize the region around phosphorene sheet for the MPM due to using high-order basis functions. Note that N = 6 (i.e., 26 unknowns) is sufficient if we require the accuracy of ne of 10−3. The field profiles (Hx) obtained by this work and the analytical solutions are shown in Figs. 3(a) and 3(b) for the symmetric and asymmetric SPP modes, respectively. In the zoomed-in view, the field profiles within the phosphorene sheet indicated by dashed lines clearly exhibit excellent agreement between those obtained by the MPM (symbols, ‘o’) and the analytical characteristic equation (solid line), demonstrating the high accuracy of the MPM. Note that the symmetric mode clearly exhibits better confinement than the asymmetric one in Fig. 3.
Fig. 3. (a) Symmetric and (b) asymmetric SPP mode profiles along the y-direction with the corresponding zoomed-in views of a BP-based slab waveguide coupler with g = 50 nm for the ZZ direction aligning to the propagating direction obtained by the MPM and the analytical solutions [39].
For the AC crystal direction aligning to the z-direction, the similar convergent behaviors of the ne as well as that in the ZZ condition are also obtained by the MPM as shown in Figs. 4(a) and 4(b) for the symmetric and asymmetric modes, respectively.
Fig. 4. Relative errors of the Re(ne) and Im(ne) of the MPM comparing with the analytical solutions [39] for the (a) symmetric and (b) asymmetric SPP modes of a BP-based plasmonic waveguide, while the AC direction aligning to the direction of light propagating.
The field profiles (Hx) obtained by the MPM and the analytical solutions are shown in Figs. 5(a) and 5(b) for the symmetric and asymmetric SPP modes, respectively, also verifying the high accuracy of the MPM as those of the ZZ condition. The weaker confinements of the AC condition than that of the ZZ one can be observed by comparing the field profiles in Figs. 3 and 5.
Fig. 5. (a) Symmetric and (b) asymmetric SPP mode profiles along the y-direction with the corresponding zoomed-in views of a BP-based slab waveguide coupler with g = 50 nm for the AC direction aligning to the propagating direction obtained by the MPM and the analytical solutions [39].
To analyze the mode couplings of the BP-based waveguide coupler, we show the Re(ne) of the symmetric and asymmetric modes for the ZZ and AC crystal directions as a function of g obtained by the MPM and the analytical solutions in Fig. 6(a). The calculated results obtained by the MPM show excellent agreement with the analytical solutions. We observe that decoupling gap between the symmetric and asymmetric SPP modes is longer for the AC condition than that for the ZZ one, and the phenomenon is consistent with the weaker mode confinements of the AC condition as shown in Fig. 5. In addition, the Im(ne) of the AC condition exhibits lower energy loss than the ZZ one as shown in Fig. 6(b). The results demonstrate once again that the present MPM can accurately model the plasmonic waveguides involving various 2D materials.
Fig. 6. (a) Re(ne) and (b) Im(ne) as a function of gap, g, between the two BP sheets obtained by the MPM and the analytical solution [39] for a layered BP-based waveguide coupler.
4.2 Dielectric loaded graphene plasmonic waveguide (DLGPW)
A 3D DLGPW structure is shown in Fig. 7(a), in which a graphene layer with thickness of tg is sandwiched by the substrate and dielectric strip with width of w and height of h. The refractive indices of the substrate and the dielectric strip are ns and nd, respectively. Here, we model a graphene layer as a uniaxial anisotropic permittivity by assuming that the graphene layer has a finite thickness tg [40,41]. Therefore, an equivalent permittivity of graphene can be expressed as ε┴ = 2.5 and ε‖ = 2.5 – jσ/(ωε0tg) for the out-of-plane and in-plane permittivities, respectively. The optical conductivity σ of graphene includes intraband and interband contributions, σ = σintra + σinter, which can be calculated by the local random-phase approximation [42]:
Fig. 7. Schematic diagrams of (a) a dielectric loaded graphene plasmonic waveguide and (b) its front view divided into 12 subdomains, where the graphene sheet is labelled by subdomains 7 to 9 (zoomed in to the top).
In this example, the subdomain extensions are classified into four types: Type I (subdomains 5 and 8) is consisted of finite extensions in both x- and y-directions; Type II (subdomains 2 and 11) is consisted of a semi-infinite extension in y-direction and a finite extension in x-direction; Type III (subdomains 4, 6, 7 and 9) is consisted of a semi-infinite extension in x-direction and a finite extension in y-direction; and Type IV (subdomains 1, 3, 10, and 12) is consisted of semi-infinite extensions in both x- and y-directions. The Hx and Hy are expanded by CPs in both x- and y-directions for Type I; CPs in the x-direction and LGFs in the y-directions for Type II; LGFs in the x-direction and CPs in the y-directions for Type III; and LGFs in both x- and y-directions for Type IV. Figures 8(a) and 8(b) show the convergences of the Re(ne) and Im(ne) versus N (we set the condition of N = mx+1 = my+1 for simplicity, although mx ≠ my is allowed) using the MPM along with the FEM solutions using the COMSOL Multiphysics. For the first three modes, the convergent effective indices (dashed lines in Fig. 11) of the mth mode obtained by the FEM with degree of freedoms (DoFs) 710610 and memory storage (RAM) 5.36 GB are $n_e^0$ = 64.878 – j1.0191 (m = 0), $n_e^1$ = 55.651 – j1.0836 (m = 1), and $n_e^2$ = 43.508 – j0.8623 (m = 2) using a workstation with the Intel Xeon CPU E5-2687W 3.40 GHz spending the calculating time 86 seconds.
Fig. 8. (a) Re(ne) and (b) Im(ne) versus N obtained by the MPM along with the reference values (dashed lines) obtained by the FEM using the COMSOL Multiphysics.
We observe that N = Nx = Ny = 8 is sufficiently to achieve excellent agreement with the results obtained by the FEM. Meanwhile, the DoFs and RAM of the present MPM are only 1276 (2 × [8 × 4–3] × [8 × 3–2]) and 72 MB spending the calculating time 2.4 seconds, respectively. In contrast, the FEM requires much more DoFs = 710610 and RAM = 5.36 GB. The significant savings of the computational efforts demonstrate the highly efficient MPM even for solving the 2D material-based plasmonic devices. Interesting, differing from the SPPs built in noble metals that tighter mode confinement accompanies by the higher energy loss [5,6], the fundamental mode (m = 0) with the strongest confinement has a lower energy loss than that of the first-order mode (m = 0). To clearly illustrate the mode distributions, the Hx field contours of the three modes, containing 50 isolines from normalized amplitude of −1 to 1, obtained by the MPM are shown in Fig. 9. Mode 0 with larger Re(ne) indeed shows the better mode confinement than other high-order modes with smaller Re(ne).
Fig. 9. The magnetic field (Hx) profiles of the (a) fundamental (m = 0), (b) first-order (m = 1), and (c) second-order (m = 2) of the DLGPW obtained by the MPM with N = Nx = Ny = 8 at the condition of w = 200 nm, h = 100 nm, EF = 0.5 eV, µ = 1 m2/V·s, and λ = 10 µm.
After verifying the efficiency and accuracy of the MPM, we further investigate the wavelength dependence on the graphene-based SPP (GSPP) mode properties. The Re(ne) and the propagation length Lp (= λ/[4πIm(ne)]) for the first four modes using the terms of Nx = Ny = 8 as a function of λ are shown in Figs. 10(a) and 10(b), respectively. Similar to the general dielectric waveguides [45], the Re(ne) of the GSPP modes increases as the wavelength decreases, and the high-order modes possess their corresponding cutoff wavelengths except the fundamental mode [13]. We observe that Lp of the fundamental GSPP mode with no cutoff wavelength linearly increases as λ increases. For other higher order modes, Lp’s also linearly increase as λ increases but dramatically increase as λ approaches the cutoff wavelengths because of the extremely loose mode distributions. For example, Lp of the mode 1 (2) increases from 1.07 (0.49) to 1.95 (1.12) µm from λ = 12.5 (8.2) to 14.5 (10.2) µm, respectively, in contrast, the corresponding Re(ne) decreases moderately from 36.49 (66.98) to 28.80 (41.88).
Fig. 10. (a) Re(ne) and (b) Lp of the first four modes of the DLGPW structure as a function of wavelength λ.
To analyze the EF dependence on the mode properties, we focus on the fundamental mode because of the preferred single-mode operation in most photonic devices. For different values of EF, the Re(ne) and Lp versus λ are shown in Figs. 11(a) and 11(b), respectively. We observe that Re(ne) increases as EF decreases because of the tighter mode confinement, however, the price reflects in the shorter Lp.
Fig. 11. (a) Re(ne) and (b) Lp of the fundamental mode of the DLGPW structure as a function of wavelength λ for different values EF.
In addition to tuning EF, Lp of the GSPP modes can be improved by controlling the carrier mobility μ in graphene. Experimentally, μ can be increased by reducing defects in lattices with advanced fabrication processes or building suspended graphene environment [44,46]. Considering the presently achievable range of μ [44], the Re(ne) and Lp as a function of λ for several values of μ are shown in Figs. 12(a) and 12(b), respectively, at the condition of EF = 0.5 eV. We observe that the Re(ne) has nothing occurred for different values of μ because μ devotes merely to the carrier collision frequency, but Lp is significantly increased by the larger μ. At the wavelength of λ = 16 µm, Lp is increased from ∼2 µm to ∼6 (8) µm by improving μ from 1 to 6 (12) m2/V·s. Therefore, increasing μ can effectively increase Lp of the GSPP modes without paying the penalty of weaker mode confinement, especially for the longer wavelength.
Fig. 12. (a) Re(ne) and (b) Lp of the fundamental modes of the DLGPW structure versus wavelength λ for different carrier mobilities μ.
Differing from tuning the material properties of EF and μ in graphene sheets, we tailor the mode properties by varying the device geometry. Figures 13(a) and 16(b) show the Re(ne) and Lp versus the width, w, of the dielectric ridge for different values of h at the condition of EF = 0.5 eV, μ = 1 m2/V·s, and λ = 10 µm. It can be seen that Re(ne) gradually increases as w increases and reaches to saturation for w >300 nm or h >50 nm because most optical power is confined in the dielectric ridge. By contrast, Lp decreases as w increases and also reduces to saturation for w >150 nm. At a fixed w while w > 150 nm, Lp moderately increases as h increases because of pulling a small amount of the mode profile away from the interface of graphene and dielectric ridge, leading to the reduction of energy loss.
Fig. 13. (a) Re(ne) and (b) Lp of the fundamental mode of the DLGPW structure versus width w of the dielectric ridge for several heights of h at the condition of EF = 0.5 eV, μ = 1 m2/V·s, and λ = 10 µm.
5. Summary
To deeply[ study the mode properties of photonic devices involving 2D materials, numerical analyses are vital and indispensable before performing experimental investigates. In this work, a multidomain pseudospectral method (MPM), which was previously reported by our group to efficiently solve general dielectric waveguides with comparable geometry dimensions, is extended to investigate the mode characteristics of 2D material-based plasmonic waveguides with the dimension contrast of the order of 10−3. To demonstrate the applicability of the MPM, we solve a BP waveguide coupler with zigzag and armchair crystal directions by the MPM, and the numerical results demonstrate the exponential accuracy comparing with the exact analytical solutions. Furthermore, we further analyzed a 3D dielectric-loaded graphene plasmonic waveguide with no analytical solutions. Comparing with the prominent FEM requiring 710610 unknowns and 5.3 GB memory storage to achieve convergent results, the present MPM requires only 1276 unknowns and 72MB memory storage. The significant reduction of the computational costs verifies that the MPM can still preserve fast convergence and high accuracy, as well as those in solving the general dielectric waveguides previously, for solving the graphene plasmonic waveguides. In addition, the mode properties of the graphene are also addressed in detailed by varying the Fermi energy levels, carrier mobilities, and geometry dimensions of the dielectric strip, exhibiting the flexible tunability of the graphene. The present MPM can be directly applied to another 2D materials without any difficulty, making the research topic of 2D material-based devices move forward swiftly.
Funding
Ministry of Science and Technology, Taiwan (109-2112-M-005-005).
Acknowledgments
The authors would like to thank Enago (www.enago.tw) for the English language review.
Disclosures
The author declares no conflicts of interest.
References
1. K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. C. Neto, “2D materials and van der Waals heterostructures,” Science 353(6298), aac9439 (2016). [CrossRef]
2. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, M. I. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, “Two-dimensional gas of massless Dirac fermions in graphene,” Nature 438(7065), 197–200 (2005). [CrossRef]
3. F. Bonaccorso, Z. Sun, T. Hasan, and A. C. Ferrari, “Graphene photonics and optoelectronics,” Nat. Photonics 4(9), 611–622 (2010). [CrossRef]
4. B. Deng, R. Frisenda, C. Li, X. Chen, A. C. Gomez, and F. Xia, “Progress on black phosphorus photonics,” Adv. Opt. Mater. 6(19), 1800365 (2018). [CrossRef]
5. W. L. Barnes, A. Dereux, and T. W. Ebbesen, “Surface plasmon subwavelength optics,” Nature 424(6950), 824–830 (2003). [CrossRef]
6. D. K. Gramotnev and S. I. Bozhevolnyi, “Plasmonics beyond the diffraction limit,” Nat. Photonics 4(2), 83–91 (2010). [CrossRef]
7. A. N. Grigorenko, M. Polini, and K. S. Novoselov, “Graphene plasmonics,” Nat. Photonics 6(11), 749–758 (2012). [CrossRef]
8. P. A. D. Gonçalves and N. M. R. Peres, An Introduction to Graphene Plasmonics, (World Scientific, 2016).
9. Y. Fan, N. H. Shen, F. Zhang, Q. Zhao, H. Wu, Q. Fu, Z. Wei, H. Li, and C. M. Soukoulis, “Graphene plasmonics: a platform for 2D optics,” Adv. Opt. Mater. 7(3), 1800537 (2019). [CrossRef]
10. J. P. Liu, X. Zhai, L. L. Wang, H. J. Li, F. Xie, S. X. Xia, X. J. Shang, and X. Luo, “Graphene-based long-range SPP hybrid waveguide with ultra-long propagation length in mid-infrared range,” Opt. Express 24(5), 5376–5386 (2016). [CrossRef]
11. J. P. Liu, W. L. Wang, F. Xie, X. Luo, X. Zhou, M. Lei, Y. J. Yuan, M. Q. Long, and L. L. Wang, “Efficient directional coupling from multilayer-graphene-based long-range SPP waveguide to metal-based hybrid SPP waveguide in mid-infrared range,” Opt. Express 26(22), 29509–29520 (2018). [CrossRef]
12. S. B. Haghighi, R. Ghayour, and M. H. Sheikhi, “Design and analysia of low loss plasmonic waveguide and directional coupler based on pattern-free suspended graphene sheet,” Carbon 129, 653–660 (2018). [CrossRef]
13. W. Xu, Z. H. Zhu, K. Liu, J. F. Zhang, X. D. Yuan, Q. S. Lu, and S. Q. Qin, “Dielectric loaded graphene plasmon waveguide,” Opt. Express 23(4), 5147–5153 (2015). [CrossRef]
14. X. Q. He, T. G. Ning, L. Pei, J. J. Zheng, J. Li, and X. D. Wen, “Tunable hybridization of graphene plasmons and dielectric modes for highly confined light transmit at terahertz wavelength,” Opt. Express 27(5), 5961–5972 (2019). [CrossRef]
15. D. Teng, K. Wang, Z. Li, and Y. Zhao, “Graphene-coated nanowire dimers for deep subwavelength waveguiding in mid-infrared range,” Opt. Express 27(9), 12458–12469 (2019). [CrossRef]
16. Z. Wu, T. Ning, J. Li, M. Zhang, H. Su, I. L. Li, and H. Liang, “Tunable photonic-like modes in graphene-coated nanowires,” Opt. Express 27(24), 35238–35244 (2019). [CrossRef]
17. D. Teng, K. Wang, Q. Huan, W. Chen, and Z. Li, “High-performance light transmission based on graphene plasmonic waveguides,” J. Mater. Chem. C 8(20), 6832–6838 (2020). [CrossRef]
18. K. Niu, P. Li, Z. Huang, L. J. Jiang, and H. Bagci, “Numerical methods for electromagnetic modeling of graphene: a review,” IEEE J. Multiscale Multiphys. Comput. Tech. 5(4), 44–58 (2020). [CrossRef]
19. V. Nayyeri, M. Soleimani, and O. M. Ramahi, “Modeling graphene in the finite-difference time-domain method using a surface boundary condition,” IEEE Trans. Antennas Propag. 61(8), 4176–4182 (2013). [CrossRef]
20. N. Liu, G. Cai, L. Ye, and Q. H. Liu, “The efficient mixed FEM with the impedance transmission boundary condition for graphene plasmonic waveguides,” J. Lightwave Technol. 34(23), 5363–5370 (2016). [CrossRef]
21. C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang, Spectral Methods in Fluid Dynamics, Springer Series in Computational Physics, (Springer-Verlag, 1988).
22. J. P. Boyd, Lecture Notes in Engineering, Chebyshev and Fourier Spectral Methods, (Springer-Verlag, 2000).
23. C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Select. Topics Quantum Electron. 11(2), 457–465 (2005). [CrossRef]
24. C. C. Huang, “Numerical investigation of mode characteristics of nanoscale surface plasmon-polaritons using a pseudospectral scheme,” Opt. Express 18(23), 23711–23726 (2010). [CrossRef]
25. R. Song, J. Zhu, and X. Zhang, “Full-vectorial modal analysis for circular optical waveguides based on the multidomain Chebyshev pseudospectral method,” J. Opt. Soc. Am. B 27(9), 1722–1730 (2010). [CrossRef]
26. D. Song and Y. Y. Lu, “Pseudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32(8), 1624–1630 (2014). [CrossRef]
27. A. Abdrabou, A. M. Heikal, and S. S. A. Obayya, “Efficient rational Chebyshev pseudo-spectral method with domain decomposition for optical waveguides modal analysis,” Opt. Express 24(10), 10495–10511 (2016). [CrossRef]
28. S. Wu and J. Xiao, “A pseudospectral reflective beam propagation method for optical waveguides,” IEEE Photonics Technol. Lett. 29(5), 435–438 (2017). [CrossRef]
29. Y. P. Chiou, W. L. Yeh, and N. Y. Shih, “Analysis of highly conducting lamellar gratings with multidomain pseudospectral method,” J. Lightwave Technol. 27(22), 5151–5159 (2009). [CrossRef]
30. D. Song, L. Yuan, and Y. Y. Lu, “Fourier-matching pseudospectral modal method for diffraction gratings,” J. Opt. Soc. Am. A 28(4), 613–620 (2011). [CrossRef]
31. P. J. Chiang, C. P. Yu, and H. C. Chang, “Analysis of two-dimensional photonic crystals using a multidomain pseudospectral method,” Phys. Rev. E 75(2), 026703 (2007). [CrossRef]
32. D. Song and Y. Y. Lu, “Analyzing leaky waveguide modes by pseudospectral modal method,” IEEE Photonics Technol. Lett. 27(9), 955–958 (2015). [CrossRef]
33. A. J. Kobelansky and J. P. Webb, “Eliminating spurious modes in finite-element waveguide problems by using divergence-free fields,” Electron. Lett. 22(11), 569–570 (1986). [CrossRef]
34. X. Y. He and R. Li, “Comparison of graphene-based transverse magnetic and electric surface plasmon modes,” IEEE J. Select. Topics Quantum Electron. 20(1), 62–67 (2014). [CrossRef]
35. R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly re-started Arnoldi iteration,” SIAM J. Matrix Anal. & Appl. 17(4), 789–821 (1996). [CrossRef]
36. H. Liu, A. T. Neal, Z. Zhu, Z. Luo, X. Xu, D. Tomanek, and P. D. Ye, “Phosphorene: an unexplored 2D semiconductor with a high hole mobility,” ACS Nano 8(4), 4033–4041 (2014). [CrossRef]
37. T. Low, R. Roldan, H. Wang, F. Xia, P. Avouris, L. M. Moreno, and F. Guinea, “Plasmons and screening in monolayer and multilayer black phosphorus,” Phys. Rev. Lett. 113(10), 106802 (2014). [CrossRef]
38. Z. Liu and K. Aydin, “Localized surface plasmons in nanostructure monolayer black phosphorus,” Nano Lett. 16(6), 3457–3462 (2016). [CrossRef]
39. X. Lin, G. Cai, H. Chen, N. Liu, and Q. H. Liu, “Modal analysis of 2-D material-based plasmonic waveguides by mixed spectral element method with equivalent boundary condition,” IEEE Photonics Technol. Lett. 38(14), 3677–3686 (2020). [CrossRef]
40. E. H. Hwang and S. D. Sarma, “Dielectric function, screening, and plasmons in two-dimensional graphene,” Phys. Rev. B 75(20), 205418 (2007). [CrossRef]
41. L. A. Falkovsky, “Optical properties of graphene,” J. Phys.: Conf. Ser. 129, 012004 (2008). [CrossRef]
42. L. A. Falkovsky and S. S. Pershoguba, “Optical far-infrared properties of a graphene monolayer and multilayer,” Phys. Rev. B 76(15), 153410 (2007). [CrossRef]
43. L. Ju, B. S. Geng, J. Horng, C. Girit, M. Martin, Z. Hao, H. A. Bechtel, X. G. Liang, A. Zettl, Y. R. Shen, and F. Wang, “Graphene plasmonics for tunable terahertz metamaterials,” Nat. Nanotechnol. 6(10), 630–634 (2011). [CrossRef]
44. K. I. Bolotin, K. J. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, and H. L. Stormer, “Ultrahigh electron mobility in suspended graphene,” Solid State Commun. 146(9-10), 351–355 (2008). [CrossRef]
45. S. J. Orfanidis, Electromagnetic waves and antennas, (online book, 2016).
46. A. S. Postigo, J. G. W. Pérez, J. S. Penadés, A. O. Moñux, M. Nedeljkovic, R. Halir, F. E. M. Mimun, Y. X. Cheng, Z. Qu, A. Z. Khokhar, A. Osman, W. Cao, C. G. Littlejohns, P. Cheben, G. Z. Mashanovich, and Í. M. Fernández, “Mid-infrared suspended waveguide platform and building blocks,” IET Optoelectron. 13(2), 55–61 (2019). [CrossRef]
