Advanced Numerical Methods for Graphene Simulation With Equivalent Boundary Conditions: A Review

1. Introduction

Graphene [1,2,3] is an allotropy of 3-D crystal graphite, a true two-dimensional material composed of a single layer of carbon atoms. Due to its unique electrical, electromagnetic, and optical characteristics, it has attracted widespread attention, leading to the design of many new systems and equipment [4,5,6]. In 2010y, Novoselov and Heim [1] won the Nobel Prize for their research and observation of graphene properties. Graphene exhibits semi metallic properties and has a strong bipolar electric field effect. And graphene has been shown to possess electrical properties similar to semiconductors (although with a zero band gap) [7,8,9]. Due to its excellent properties, many methods have been developed to fabricate graphene, such as micromechanical cleavage technique [10], chemical vapor deposition [11], solvent exfoliation [12,13], solvothermal method [14], etc. The development of preparation technology has made it easy to separate high-quality graphene, triggering a research boom in the 2-D material family. In fact, the optical properties generated by the unique electronic band structure of graphene are considered attractive features in the design of nanophotons and optoelectronic components [15,16]. Graphene has a strong interaction with electromagnetic fields, and its response to light is nonlinear, manifested by plasmonic characteristics. Graphene-based plasmonic can not only limit the electromagnetic field to a smaller transverse propagation range, but also modulate it over a wide frequency range through gating and chemical doping [17,18,19,20].

In the past few years, in addition to the physical research on graphene and graphene-based devices, another important topic has been the numerical simulation of graphene. Kubo formula puts forward the expression of graphene conductivity, which is a function formed by physical parameters such as wavelength, chemical potential and temperature [21,22,23,24]. A surface conductivity model is also proposed to describe graphene as an isotropic infinitesimal sheet. Although the model can provide problematic results for some tests, the electromagnetic simulation of specific graphene devices still remains a challenge when it comes to practical graphene problems. Fortunately, with the development of computational electromagnetics, many efficient and accurate numerical methods can have been proposed to analyze the electromagnetic response of graphene related devices. Among them, the most popular methods are the finite difference time domain method (FDTD), finite element method (FEM), and spectral element method (SEM) [25,26,27,28,29,30]. The numerical approaches of processing graphene include: (a) Taking graphene as a zero-thickness sheet [31,32,33,34]; (b) Graphene is regarded as a thin plate with finite thickness, and its surface conductivity is converted into a volumetric dielectric constant [35]. Due to the easy realization of (b), some published papers and commercial software will model graphene as a thin sheet with limited thickness [36,37,38]. However, direct discretization of graphene thin plates result in extremely fine grids and a large number of unknowns. Especially in time-domain simulation, extremely small time steps are required to ensure stability, which consumes an enormous amount of CPU time and memory costs. Therefore, equivalent boundary conditions are preferred methods used to eliminate thin plates in the computational domain, which include impedance transmission boundary condition (ITBC) [39,40], surface current boundary condition (SCBC) [41] impedance matrix boundary condition (IMBC) [42]. These numerical methods and equivalent boundary conditions have been proved that can greatly improve the design and manufacturing speed.

In this review, we will focus on the electromagnetic simulation of graphene, and select several numerical cases to show their serviceability in practical engineering applications.

2. Advanced Numerical Methods

With the development of computer performance and new efficient and accurate algorithms, electromagnetic numerical methods have become more applicable and efficient. In order to analyze the electromagnetic field interaction related to graphene and the shielding effect of graphene, several different numerical methods can be used, including mixed finite element method (Mixed FEM), the mixed spectral element method (Mixed SEM), auxiliary source method (MAS) and interior penalty discontinuous Galerkin time domain (IPDG), and equivalent boundary conditions such as ITBC, SCBC and IMBC are introduced. Regardless of the boundaries used by these methods, a mathematical model is required to describe the conductivity of graphene. The following subsection depicts the mathematical model employed in this review.

2.1. Mathematical Model

The numerical methods listed above are commonly used to simulate graphene as a thin sheet with a conductivity that is both complex and dependent on frequency. The graphene surface conductivity is usually expressed using the Kubo formula [43,44,45]:

$$\begin{array}{cc}\hfill {\sigma }_g=& \frac{j{e}^2(\omega -j2\Gamma )}{\pi {\hslash }^2}[-{\int }_0^{\infty }\epsilon \frac{f_d(-\epsilon )-f_d\left(\epsilon \right)}{{(\omega -j2\Gamma )}^2-4{\left(\frac{\epsilon }{\hslash }\right)}^2}d\epsilon \hfill \\ \hfill & +\frac{1}{{(\omega -j2\Gamma )}^2}{\int }_0^{\infty }\epsilon (\frac{\partial f_d\left(\epsilon \right)}{\partial \epsilon }-\frac{\partial f_d(-\epsilon )}{\partial \epsilon })d\epsilon ]\hfill \end{array}$$

(1)

where $\epsilon$ is the electron charge, $\omega$ is the radian frequent, $\Gamma$ is the phenomenological scattering rate, ℏ is the reduced Planck constant and e is the energy state. $f_d={(e^{(\epsilon -\left|{\mu }_c\right|)/k_{B}T}+1)}^{-1}$ is the Fermi-Dirac distribution, where ${\mu }_c$ is the chemical potential or Fermi level, T is the temperature and $k_B$ is the Boltzman constant. Equations (1) contains both intraband and interband contributions that can be expressed as ${\sigma }_g={\sigma }_{\mathrm{intra}}(\omega ,{\mu }_c,\Gamma ,T)+{\sigma }_{inter}(\omega ,{\mu }_c,\Gamma ,T)$, where:

$${\sigma }_{\mathrm{intra}}(\omega ,{\mu }_c,\Gamma ,T)=\frac{-2j{e}^{2}T}{{\hslash }^2\pi (\omega -j\Gamma )}\ln \left[2\cos \mathrm{h}\left(\frac{{\mu }_c}{2T}\right)\right]$$

(2)

$${\sigma }_{\mathrm{inter}}(\omega ,{\mu }_c,\Gamma ,T)=\frac{e^2}{4\hslash }\left[\frac{1}{2}+\frac{1}{\pi }\arctan \left(\frac{\hslash \omega -2{\mu }_c}{2T}\right)+\frac{j}{2\pi }\ln \left\{\frac{{(\hslash \omega +2{\mu }_c)}^2}{{(\hslash \omega -2{\mu }_c)}^2+{\left(2T\right)}^2}\right\}\right]$$

(3)

The interband conductivity, which is smaller than the intraband conductivity, is in the order of magnitude of $e^2/\hslash$. The dielectric constant of graphene is determined by $\widehat{\epsilon }=1-j{\sigma }_g/\left({\epsilon }_0\omega \tau \right)$. The electrons and holes close to the band edges "blocke" interband transitions at lower frequencies (i.e., in the THz range). Graphene exhibits the characteristics of a conductive film, and a straight forward Drude model is used to describe its conductivity (primarily caused by in-band contributions). The electrons and holes close to the band edges "blocke" interband transitions at lower frequencies

2.2. Mixed Finite Element Method

With the increasing complexity of microwave and optical waveguide structures, it is difficult to effectively conduct modal analysis of waveguide problem, which seriously delays the subsequent adjustment of geometric and material parameters in waveguide design. Scholars have developed different numerical methods to solve this problem, among which the FEM is a successful representative. The mathematical theory of the finite element method was completed by Feng Kang and others scholars. Compared with the FDTD, the FEM can utilize unstructured meshes to discretize the computational domain, which is convenient for solving problems in complex calculation regions. However, the traditional FEM based on scalar node basis functions inevitably generate spurious modes with non-zero eigenvalues. The FEM based on vector basis functions is a reliable solution, where there are no non-zero spurious solutions for solving Maxwell problems, but it will generate zero spurious solutions (also known as DC spurious modes) [46]. All modes can be suppressed by a mixed FEM based on a new weak form.

The mixed FEM is a powerful technique that not only maintains the flexibility of accurately modeling complex geometries shapes with finite elements, but also suppresses all spurious modes with new weak forms [47,48,49]. It applies systematic and rigorous mathematical techniques to the solution of equations and boundary conditions and uses meshes (for example, using triangle and quadrilateral elements) to discrete the computational domain involving geometrically complex structures that have potentially heterogeneous material properties. The medium in a waveguide structure is a two-dimensional variable with a transverse anisotropic tensor, and its longitudinal component is constant where the electromagnetic field is a three-dimensional function (assumed to be in the z-direction). Thus, the electric field in an arbitrary cross section is described as $\mathit{e}(x,y)e^{-j{k}_{z}z}$. This results in the expression of the vectorial Helmholtz equation and Gauss’ law can be written as:

$$\begin{array}{cc}\hfill {\nabla }_t\times {\mu }_{rz}^{-1}{\nabla }_t\times {\mathit{e}}_t+j{k}_z\widehat{z}\times & {\mathit{\mu }}_{rt}^{-1}\widehat{z}\times {\nabla }_t{e}_z-k_z^2\widehat{z}\times {\mathit{\mu }}_{rt}^{-1}\widehat{z}\times {\mathit{e}}_t-k_0^2{\epsilon }_{rt}{\mathit{e}}_t=0\hfill \\ \hfill {\nabla }_t\times {\mathit{\mu }}_{rt}^{-1}(\widehat{z}& \times {\nabla }_t{e}_z+j{k}_z\widehat{z}\times {\mathit{e}}_t)+k_0^2{\epsilon }_{rt}{e}_z\widehat{z}=0\hfill \\ \hfill & {\nabla }_t·{\epsilon }_{rt}{\mathit{e}}_t-j{k}_z{\epsilon }_{rz}{e}_z=0\hfill \end{array}$$

(4)

where ${\mu }_r=\mathrm{diag}({\mu }_{rt},{\mu }_{rz})$ and ${\epsilon }_r=\mathrm{diag}({\epsilon }_{rt},{\epsilon }_{rz})$ are the magnetic permeability and relative dielectric permittivity tensors with the transversely anisotropic medium ${\mu }_{rt}$ and ${\epsilon }_{rt}$, respectively. $k_0$ is the wavenumber in vacuum.

The graphene layer is assumed to have a finite thickness and a specific volumetric conductivity.

2.2.1. Impedance Transmission Boundary Condition

In the traditional FEM, due to its simplicity and feasibility, the graphene is assumed to be a layer of finite-thickness with a specific volumetric conductivity [50]. However, a graphene sheet having finite thicknesses results in exceedingly small meshes and a substantial number of unknowns when they are directly discretized. Therefore, to ensure numerical accuracy, a large amount of CPU time and memory costs are required. The impedance boundary condition (IBC) is a useful way for eliminating thin plates from the computational field, such as impedance network boundary condition (INBC), impedance transmission boundary condition (ITBC) and surface impedance boundary condition (SIBC) [51,52]. By using two-port network equations to describe the field inside the conductive sheet, INBC is effectively implemented in both FDTD and FEM. Meanwhile, Leontovich proposed a hypothetical SIBC for thin plates that the thickness is less than or equal to the skin depth. However, when the frequency of graphene reaches the terahertz band, the thickness of graphene is much lower than the skin depth. Therefore, ITBC is a better fit method to describe the interdependent tangential electromagnetic field on each side of graphene surface [53]. In order to clearly describe the implementation of ITBC, Figure 1 uses two triangular elements that share ITBC to illustrate.

The prerequisite for employing ITBC is that the wavelength and skin depth $\delta$ are larger than the conducting sheet’s thickness $\tau$. As shown in Figure 2, it can be observed that using line segments to represent thin conductive sheets with finite thickness in a 2-D waveguide cross-section model [50]. The condition for using ITBC is that the wavelength and the skin depth $\delta$ are greater than the conducting sheet’s thickness $\tau$, while avoiding sharp edges and angles. At this dot, the electromagnetic field is transmitted on either side of a sheet. Meanwhile, considering the interdependence of electromagnetic fields above and below the thin film, ITBC uses transmission line theory to represent the relationship between tangential electromagnetic fields above and below the line segment L:

$$\left[\begin{array}{c}{\mathit{n}}_1\times {\mathit{h}}_1\\ {\mathit{n}}_2\times {\mathit{h}}_2\end{array}\right]=\left[\begin{array}{cc}{Y}_{11}& Y_{12}\\ Y_{12}& Y_{11}\end{array}\right]\left[\begin{array}{c}{\mathit{n}}_1\times {\mathit{n}}_1\times {\mathit{e}}_1\\ {\mathit{n}}_2\times {\mathit{n}}_2\times {\mathit{e}}_2\end{array}\right]$$

(5)

Graphene nanoribbon sandwich waveguide [54] is considered which placed in vacuum background, with two identical graphene ribbons on both sides and a dielectric strip with a refractive index of $n_d$ = 1.4 in the middle as shown in Figure 3 [50]. In the model, the graphene is an isotropic medium which thickness is 0.5 nm and $\Gamma =0.1$ meV, ${\mu }_c=0.5$ eV, and $T=300$ K. The computing domain of 20 µm × 20 µm yields 11 212 elements, including 81 011 unknowns. Figure 4 shows the calculation results $\left|k_z/k_0\right|$ for the first twelve modes. We can clearly see that the mixed FEM-ITBC results and FEM solutions have good consistency. In addition to the accuracy of the mixed FEM-ITBC, Table 1 also shows that the mixed FEM-ITBC incurs lower memory costs and requires less CPU time than the traditional FEM.

2.2.2. Surface Current Boundary Condition

The previous ITBC is only applicable to model isotropic media. The modeling thin layers of anisotropic media requires the use of SCBC. In the right figure of Figure 2, one-dimensional lines are used to replace thin conductive sheets with finite thickness. When graphene thin layers are equivalent to zero thickness [33,55], fine mesh generation can be avoided. SCBC can be expressed as:

$$\mathit{n}\times ({\mathit{h}}_2-{\mathit{h}}_1)={\mathit{J}}_s={\mathit{\sigma }}_s\mathit{e}$$

(6)

where $\mathit{e}$ is the electric field, ${\sigma }_s$ is the surface conductivity of the nanomaterial sheet which can be anisotropic media. $\mathit{n}={\mathit{n}}_1=-{\mathit{n}}_2$ is the unit normal vector from ${\Omega }_1$ to region 2. ${\mathit{h}}_1,{\mathit{h}}_2$ are the magnetic fields of ${\Omega }_1$ and ${\Omega }_2$, respectively.

Figure 5 shows a graphene-coated dielectric nanowire waveguide (GNW), with an inner layer of dielectric nanowires and an outer layer covered with graphene sheet. The nanowire with a radius of $R=100$ nm has a relative permittivity of ${\epsilon }_{r1}=3$ and the external region has a relative permittivity of ${\epsilon }_{r2}=1$. The surface conductivity s of graphene ${\sigma }_s$ is calculated using Kubo’s formula with ${\mu }_c=0.5$ e V, $T=300$ K, and $\tau =200$ ns. By using PML which thickness is 0.5 µm to truncate the open GNW structure. The length and width of the calculation area are both 5 µm. Through comparing the results calculated by the mixed FEM-SCBC method with the analysis results [58], the field distributions of the first three propagation modes (m0, m1, m2), propagation lengths $L=1/\left[2\mathrm{Im}\left(k_z\right)\right]$, and the dispersion relations as a function of frequency is investigated. The proposed method has been proven to be accurate, as presented in Figure 6a,b, with maximum relative errors of 0.4 % and 0.15 %, respectively. The field distributions of the modes (m0, m1, m2) in Figure 6c are also in well-alignment with the field distribution of the analytical method. Furthermore, Table 2 demonstrates that the mixed FEM-SCBC is also a more valuable option to studying the propagation characteristics of plasmonic waveguides which is based on graphene, as it requires only 91 997 DoFs and avoids the demand for a dense mesh around the graphene. In contrast, the traditional FEM requires 384 737 DoFs to discretize thin graphene sheets with a thickness of 0.5 nm. In summary, this method offers both effectiveness and computational efficiency for the analysis of graphene-based plasmonic waveguides.

2.3. Mixed Spectral Element Method With SCBC

The SEM which combines spectral method with FEM extended by Patera. Due to its unique basis function construction method and node arrangement, it can achieve higher accuracy than the FEM, and consumes less memory costs and computational time. For solving Stokes equations in fluid mechanics, Patera took the lead in introducing spectral methods into calculations, resulting in better numerical accuracy and convergence compared to traditional methods [59]. Since then, the SEM has been widely used in numerical calculations simulation in the field of fluids. Later, pseudospectral method was proposed by Chebyshev, and Liu incorporated the pseudospectral method into the FEM [60]. The SEM has become a new and efficient tool for electromagnetic numerical computation.

The mixed SEM possesses both spectral accuracy and the capability to eliminate zero spurious eigenvalues [61]. The combination of SCBC and mixed SEM is utilized to analyze the modes of graphene based pasmonic waveguides due to their enormous potential in eigenvalue solver. This method employs a new variational formulation to suppress zero spurious modes by introducing Gauss’ law to the vectorial Helmholtz equation. Additionally, an equivalent SCBC is used in the computational domain to replace the nanoscale graphene sheets.

A tunable multilayer nanoring waveguide consisting of six layers of graphene and seven layers of dielectric is firstly considered. The graphene and dielectric layers are alternately distributed [62]. Figure 7 depicts the cross-section and three-dimensional structure of the waveguide. The innermost ring radius is 140 nm, and the total ring radius is 210 nm. The innermost and outermost layers have a thickness of 10 nm, while the rest layers and graphene sheets are 9.5 nm and 0.5 nm thick, respectively. The dielectric layer has a relative permittivity of 2.1. In this situation, the environment temperature is $T=300$, $\tau =0.5$ ps and ${\mu }_c=0.5$ eV. The relative dielectric constant of graphene is expressed as ${\epsilon }_r=2.5-j{\sigma }_s/\left({\epsilon }_0\omega d\right)$ [63], where d represents the thickness of graphene.

The multilayer nanoring open waveguide was truncated using a PML with 100 nm thickness, and the entire simulation domain size was 1.2 µm × 1.2 µm. To obtain accurate simulation results for the graphene sheets, conventional FEM requires a dense finite element mesh with 127 620 triangular elements, yielding 893 421 unknowns. Instead, the mixed SEM-SCBC method substitutes SCBC for graphene sheets in the computational domain. As a result, only 1449 quadrilateral elements are produced across the entire domain, and only 17 528 of them are unknown. Figure 8a displays the extremely fine meshes of graphene sheet, where at least one element is discretized in the graphene sheet, while Figure 8b displays the quadrilateral elements that result from SCBC, which are relatively coarse. It is evident that using SCBC instead of graphene sheets can significantly decrease on the number of elements and degrees of freedom (DOF) required for numerical calculations.

The mixed SEM-SCBC is employed to simulate the numerical dispersion relationship of graphene waveguides. Figure 9 shows the relationship between the real part of the effective refractive index (real ($N_{\mathrm{eff}}$, $N_{\mathrm{eff}}=k_z/k_0$) and the operating frequency. When m = 0 is a fixed, the real part of the effective refractive index varies with n, with mode $(0,1)$ being the lowest among modes $(0,n)$. The propagation lengths, represented by $L_p={\lambda }_0/4\pi \left[\mathrm{Im}\left(N_{\mathrm{eff}}\right)\right]$ in Figure 10, are calculated based on the operating frequency. It can be observed that higher order modes experience more losses compared to lower order modes, which is in line with the characteristics of propagation modes. In this figure, the numerical results of traditional FEM are presented by –, – –, and · lines. The proposed mixed SEM-SCBC method has great consistency with the traditional FEM with graphene sheet, but with significantly lower CPU time and memory costs (only 3% and 2%, respectively), as shown in Table 3. These results demonstrate that the mixed SEM-SCBC can enhance computing speed and save computer resources.

2.4. The Method of Auxiliary Sources with IMBC

The electromagnetic waves are incident on a conductive layer that separates the layered medium in a normal shielding configuration. The conductive layer can block the incident field and act as an electromagnetic shield. The planar, cylindrical, and spherical interfaces are the most common layered media models used to demonstrate electromagnetic shielding. However, for the high conductivity of interior materials, the conventional MAS [67,68] approach will yield divergent results. In Figure 11, if the standard MAS described in [69,70] is used to calculate the shielding effect, a convergent field cannot be obtained within the highly conductive layer. Therefor, a modified version of MAS was proposed in [71] to improve by associating the field at each point of boundary a with the electric field at the corresponding point of boundary b. This effect can be achieved by replacing the conductive layer with the IMBC. The shielding effectiveness of a single layer of graphene covering a cylindrical structure can be calculated using this modified method.

By superposing the fields emitted by $N_1$ and $N_2$ filaments set on the auxiliary surfaces $C_{aux}^1$ and $C_{aux}^2$, it is possible to replicate the internal field in region 2 and the scattering field in region 1. The $N_1$ and $N_2$ filaments, which radiate in the dielectric-filled unbounded space of regions 1 and 2, respectively, and carry unknown currents $I_{1l}\left(K_{1l}\right),l=0,...,N_1-1$ and $I_{2l}\left(K_{2l}\right),l=0,...,N_2-1$, respectively. The field of the l-th source produced in region $m(m=1,2)$ is expressed as:

$$E_{ml}(x,y)=-\widehat{\mathit{z}}\frac{k_m{Z}_m}{4}{I}_{ml}{H}_0^{\left(2\right)}\left(k_m{\rho }_{ml}\right)$$

(7)

$$H_{ml}(x,y)=\frac{k_m}{4i}{I}_{ml}\frac{\widehat{\mathit{x}}(y_{ml}-y)+\widehat{\mathit{y}}(x-x_{ml})}{{\rho }_{ml}}\times H_1^{\left(2\right)}\left(k_m{\rho }_{ml}\right)$$

(8)

for the ${\mathrm{TM}}_z$ case, and

$$H_{ml}(x,y)=-\widehat{\mathit{z}}\frac{k_m}{4Z_m}{K}_{ml}{H}_0^{\left(2\right)}\left(k_m{\rho }_{ml}\right)$$

(9)

$$E_{ml}(x,y)=\frac{k_m}{4i}{K}_{ml}\frac{\widehat{\mathit{x}}(y_{ml}-y)+\widehat{\mathit{y}}(x-x_{ml})}{{\rho }_{ml}}\times H_1^{\left(2\right)}\left(k_m{\rho }_{ml}\right)$$

(10)

for ${\mathrm{TE}}_z$ polarization, where $(x_{ml},y_{ml})$ denotes filament coordinates in region m, ${\rho }_{ml}$ denotes the separation between the observation point $(x,y)$ and the l-th filament in region m, $I_{ml}$ and $K_{ml}$ denote the current and magnetic current amplitudes of the l-th filament in region m, $k_m=\omega \sqrt{{\mu }_m/{\epsilon }_m}$ and $Z_m=\sqrt{{\mu }_m/{\epsilon }_m}$ denote the wave number and impedance of region m, respectively.

Then, the 2-D-IMBC approach is applied to M matching pairs. This is done to associate the fields at specific points on boundary a with the fields at corresponding points on boundary b.

$$\left[\begin{array}{c}{\widehat{\mathit{n}}}_a\times ({\widehat{\mathit{n}}}_a\times {\mathit{E}}_a)\\ {\widehat{\mathit{n}}}_b\times ({\widehat{\mathit{n}}}_b\times {\mathit{E}}_b)\end{array}\right]=\left[\begin{array}{cc}-Z_S& \sqrt{\frac{r_b}{r_a}}{Z}_{tr}\\ -\sqrt{\frac{r_a}{r_b}}{Z}_{tr}& Z_S\end{array}\right]\left[\begin{array}{c}{\widehat{\mathit{n}}}_a\times {\mathit{H}}_a\\ {\widehat{\mathit{n}}}_b\times {\mathit{H}}_b\end{array}\right]$$

(11)

where $Z_{tr}=-i{Z}_3\csc (k_3▵)$ and $Z_s=-i{Z}_3\cot (k_3▵)$ are the transfer and surface impedance, $k_3$ and $Z_3$ denote the wavenumber and characteristic impedances of the conductive layer, the curvature radii on boundaries a and bare $r_a$ and $r_b$, respectively. The appendix of [73] shows that ${\widehat{\mathit{n}}}_a=\widehat{\mathit{n}}={\widehat{\mathit{n}}}_b$ and $r_a=▵+r_b$ for uniformly thick conductive layers.

Figure 12 studies a multi-layer cylindrical medium with a circular square cross-section. Wherein region 3 is a graphene monolayer, region 2 is a silicon dioxide (SiO2) layer, air is filled in both regions 1 and 4; thus, ${\epsilon }_1={\epsilon }_4={\epsilon }_0$ and ${\epsilon }_2=3.9{\epsilon }_0$. The surface conductivity of graphene with $E_F=0.5$ eV, $\tau =6.582$ ps, $T=300$ K and $▵=0.5$ nm is utilized. Boundary b has collocation points with the following coordinates:

$$\begin{array}{cc}\hfill & x_b=15\mathrm{sgn}\left(\cos \theta \right){\left|\cos \theta \right|}^{\frac{3}{2}}\hfill \\ \hfill & y_b=15\mathrm{sgn}\left(\sin \theta \right){\left|\sin \theta \right|}^{\frac{3}{2}}\hfill \end{array}$$

(12)

A TMz or TEz plane wave incident the computational domain. There are four sets of auxiliary sources $N_1=N_2=N_3=N_4=110$, each placed on coordinate sources $(x_{a\_1},y_{a\_1})=0.8(x_b,y_b),(x_{a\_2},y_{a\_2})=2(x_b,y_b),(x_{a\_3},y_{a\_3})=0.8(x_c,y_c)$ and $(x_{a\_4},y_{a\_4})=1.2(x_c,y_c)$ of $C_{aux}^1,C_{aux}^2,C_{aux}^3$, and $C_{aux}^4$. Applying the boundary condition Equation (10) to a and b on the $M=110$ matching pair. The norm of the fields in regions 1 and 4 with $(x_{o\_1},y_{o\_1})=2(x_b,y_b)$ and $(x_{o\_2},y_{o\_2})=0.6(x_c,y_c)$ are described in Figure 13. It is also shown that the modified MAS in good agreement with the finite element calculation results of COMSOL. The relationship between the shielding effect with $E_F=0.5$ and f, and the relationship between the shielding effect with $f=3$ GHz and $E_F$ are presented in Figure 14 and Figure 15, respectively. The modified MAS not only matches the finite element calculation results, but also has advantages in memory usage. The shielding effect is shown in Figure 16.

2.5. Discontinuous Galerkin Time-Domain Method With SIBC

The discontinuous Galerkin time-domain method [34,74] is a versatile approach for solving differential equations in numerous fields such as computational science, engineering, and physics. It incorporates the advantages of the Finite Volume Method (FVM) [75,76] and the FEM, enabling mesh discretization of the computational domain. The spatial DGTD operations, like FVM, are localized, and the global mass matrix can be transformed and divided into a block diagonal mass matrix. The mass matrix block’s inversion and storage [77] are performed before initiating time marching. This makes the solver of DGTD very compact, especially when using explicit integration methods. These characteristics make DGTD an excellent method for simulating multi-scale electromagnetic fields in two-dimensional materials.

Since the jumping depth of graphene is much greater than the thickness of the graphene layer, the SIBC can be utilized to replace the graphene layer, which can be expressed as:

$$\begin{array}{cc}\hfill & \widehat{\mathit{n}}\times ({\mathit{E}}^2-{\mathit{E}}^1)=0\hfill \\ \hfill & {\sigma }_g^{-1}·[\widehat{\mathit{n}}\times ({\mathit{H}}^2-{\mathit{H}}^1)]={\sigma }_g^{-1}·\mathit{J}=\mathit{E}\hfill \end{array}$$

(13)

where ${\mathit{E}}^1$ and ${\mathit{E}}^2$ indicate the total electric field intensity on each side of the graphene layer, respectively. ${\mathit{H}}^1$ and ${\mathit{H}}^2$ represent the total magnetic field intensity on each side of the graphene layer, respectively. $\widehat{\mathit{n}}$ indicates the unit normal vector above the graphene layer, pointing outwards from the plane. $\mathit{J}$ is the induced polarization current density in the graphene layer and $\mathit{E}$ represents the intensity of the electric field.

A graphene nano-ribbon transmission line is simulated and the electromagnetic field interaction on graphene nano-ribbon (GNR) located on silicon substrate is analyzed. As is shown in Figure 17, the GNR dimension is $0.3$ µm × 3 µm × 0.05 µm. Port 1 is the current source, and both non-local and local conductivity models are simulated. Time-voltage is recorded at ports 1 and 2 of the nanobelt transmission line in two simulation models (see Figure 18). Figure 19 and Figure 20 show the relationship between $Z_{11}$ and $Z_{21}$ with frequency. Results in the time domain differ significantly from those in the frequency domain, especially in the high-end terahertz band.

2.6. Interior Penalty Discontinuous Galerkin-Time Domain Method With ITBC

The IPDG, a DG method, is also a good method for solving electromagnetic problems, including modeling graphene as an infinitely thin impedance surface using ITBC [79], which minimizes memory usage and computation time.

Equation (1) expresses the surface conductivity as a frequency-domain expression, which is difficult to solve in time domain. To overcome this, the vector-fitting technique [80] is employed to approximate ${\sigma }_g$ as the sum of partial fractions of the extreme residual pairs, either in the form of real or complex conjugate, on the frequency band. Equation (1) is modified as:

$${\sigma }_g=\sum _{l=1}^N{\sigma }_l=\sum _{l=1}^N\frac{r_l}{j\omega -p_l}$$

(14)

where N is the total number of poles and $p_l$ and $r_l$ are residuals and poles, respectively. Assuming a computational domain Grammar which has boundary ${\Gamma }_b$. The entire computational domain is discretized into multiple non-overlapping tetrahedrons, with ${\Gamma }_{interior}$ representing all inner surfaces and ${\Gamma }_{ITBC}$ representing the inner surfaces attached to the graphene surface. The enhanced IPDG method using ITBC can be reformulated as:

$$\begin{array}{cc}\hfill & {\int }_{\Omega }\epsilon \frac{{\partial }^2\tilde{\mathit{E}}}{\partial t^2}·{\mathit{N}}_{i}d\mathit{r}+{\int }_{\Omega }{\mu }^{-1}\nabla \times \tilde{\mathit{E}}·\nabla \times {\mathit{N}}_{i}d\mathit{r}-{\int }_{{\Gamma }_{ITBC}}{\mathit{J}}_S·\left\{{\mathit{N}}_i\right\}dS\hfill \\ \hfill & -{\int }_{{\Gamma }_{interior}}\left\{{\tilde{\mathit{E}}}^{*}-\tilde{\mathit{E}}\right\}·{\left[\left[{\mu }^{-1}\nabla \times {\mathit{N}}_i\right]\right]}_{T}dS-{\int }_{{\Gamma }_{interior}}\left\{{({\mu }^{-1}\nabla \times \tilde{\mathit{E}})}^{*}\right\}·{\left[\left[{\mathit{N}}_i\right]\right]}_{T}dS\hfill \\ \hfill & +{\int }_{{\Gamma }_{interior}}{\left[\left[{\tilde{\mathit{E}}}^{*}-\tilde{\mathit{E}}\right]\right]}_T·\left\{{\mu }^{-1}\nabla \times {\mathit{N}}_i\right\}dS+{\int }_{{\Gamma }_{interior}}{\left[\left[{({\mu }^{-1}\nabla \times \tilde{\mathit{E}})}^{*}\right]\right]}_T·\left\{{\mathit{N}}_i\right\}dS\hfill \\ \hfill & +{\int }_{{\Gamma }_b}\mathit{n}\times {({\mu }^{-1}\nabla \times \tilde{\mathit{E}})}^{*}·{\mathit{N}}_{i}dS+{\int }_{{\Gamma }_b}\mathit{n}\times ({\tilde{\mathit{E}}}^{*}-\tilde{\mathit{E}})·{\mu }^{-1}\nabla \times {\mathit{N}}_{i}dS=0\hfill \end{array}$$

(15)

where $\mathit{n}$ indicates the unit vector, which points towards the graphene sheet. $\tilde{\mathit{E}}={\partial }_t^{-1}\mathit{E}$ is expanded with vector basis function ${\mathit{N}}_i$. ${\left[\left[·\right]\right]}_T$ represents tangential jump, $\left\{·\right\}$ represents the average value on the face f of an element, and the superscript "*" represents numerical flux.

Assuming an infinite graphene sheet, the x direction boundary is taken as PEC while PMC is considered as the boundary in the y direction. The truncated boundary is a first-order ABC in the z direction . The entire domain size is 10.5 µm $\times 2.5$ µm $\times 2.5$ µm, and a plane wave source is created using total/scattered field boundary conditions. Unfolding unknown fields using full two vector basis functions, with a time step size of $▵t$ is $5e^{-17}$ s. The graphene parameter ${\mu }_c$ is estimated to be 0.15 eV.

3. Conclusions

In the past decade, due to its excellent optical properties, thermal, mechanical and electrical, widespread attention has been attracted on graphene, papers on the physical mechanisms and applications of graphene in various fields have been constantly updated. In this work, we reviewed some of the most commonly used algorithms for simulating graphene in the CEM field. The mixed FEM-ITBC and mixed FEM-SCBC and mixed SEM-SCBC are employed to solve the graphene-based plasmonic waveguide problem. The improved auxiliary source method with IMBC is used to calculate the shielding effect, which can obtain a convergence field in a high conductive layer; The DGTD method with SIBC is also an effective method for exploring the spatial dispersion characteristics of graphene. The IPDG-ITBC to model graphene, saving memory and usage time. All these equivalent boundary conditions have been proposed and significantly saved computational time and memory costs. With all of these performances and improvements, we anticipate that equivalent boundary conditions are promising techniques to simulate electromagnetic response of graphene.