A Review of Computational Electromagnetic Methods for Graphene Modeling

Graphene is a very promising optoelectronic material and has gained more and more attention. To analyze its electromagnetic properties, several numerical methods have been developed for graphene simulation. In this paper, a review of application of graphene in electronic and photonic device is provided, as well as some widely used computational electromagnetic algorithms for graphene modeling. The advantages and drawbacks of each method are discussed and numerical examples of these methods are given to illustrate their performance and application.


Introduction
Graphene consists of a monolayer of carbon atoms arranged on a two-dimensional honeycomb lattice which is made up of hexagons. Since its first discovery in 2004, graphene has attracted tremendous research interest in various fields due to its distinctive properties [1][2][3]. Graphene is a rising star not only in material science and condensed matter physics, but also in the electronic and photonic device communities [4].
Graphene has unique electronic band structure, and the electrons in it behave as massless Dirac fermion [5,6]. Graphene acquires a pronounced electric field effect which means carrier concentration can be tuned by electrostatic gating. It has been shown that, in bilayer graphene, the electronic gap between conduction bands and the valence can be tuned between zero and midinfrared energies, which makes bilayer graphene the only known semiconductor with a tunable energy gap [7,8]. Zhang et al. [9] presented a widely tunable electronic bandgap in electrically gated graphene, and they realized a gate-controlled, continuously tunable bandgap of up to 250 meV. Large-area graphene films of the order of centimeters on copper substrates were realized by chemical vapor deposition using methane [10], which facilitated the fabrication of graphene transistors. A good review of graphene transistors in both logic and radiofrequency applications is provided in [11].
Besides electron-device application, graphene has also been recognized as a novel optical material for photonic application. Compared to traditional metal, graphene exhibits several favorable properties. Particularly, surface plasmon polaritons (SPPs), which are electromagnetic waves coupled charge excitations on the surface of a conductor, can be excited in graphene. The plasmons in graphene are tightly confined and the volumes of SPPs in graphene can be 10 6 times smaller than those in free space [12]. This property leads to strong light-matter interaction in graphene [13]. Additionally, the dielectric properties of graphene can be electrically or chemically tuned by changing the charge carrier density and Fermi energy [13][14][15]. Brar et al. [12] created and probed plasmons in graphene with ≤ 0 /100 and resonant energies as high as 310 meV for 15 nm nanoresonators. Fei et al. [16] showed common graphenes/SiO 2 /Si back-gated structure support propagating surface plasmons and they altered both the amplitude and the wavelength of these plasmons by varying the gate voltage. The graphenebased plasmonics may enable the novel optic devices working in different frequency ranges-from terahertz to the visible range-with extremely high speed, low driving voltage, low power consumption, and compact sizes [17]. Vakil and Engheta in [18] showed that graphene can be made into a one-atom-thick platform for infrared metamaterials and transformation optical devices by manipulating spatially 2 International Journal of Antennas and Propagation inhomogeneous and nonuniform conductivity patterns. Ju et al. [19] used metamaterial made up of periodic graphene microribbon arrays for terahertz plasmon excitations and demonstrated that the plasmon resonance can be tuned over a broad terahertz frequency range by adjusting microribbon width and electrostatic doping.
Our group has been dedicated to the study of molecular sensors for a long time. Francescato et al. [20] presented a platform for broadband molecular spectroscopy based on the propagation of strongly confined antibonding plasmons supported by graphene sandwiches. This novel scheme measures the absorption spectrum of the molecule and shows a broadband capability and high sensitivity. Recently, Yang et al. [21] proposed a cylindrical graphene plasmon waveguide and investigated its application in molecular sensing.
The terahertz antenna based on graphene has promising application in wireless communications in nanosystems. The tunable property of graphene facilitates the designing of reconfigurable antenna. Huang et al. [22] presented a beam reconfigurable antenna based on a switchable high-impedance surface using single-layer graphene. The graphene-based antennas which have broad wavelength tuning range are proposed in [23,24].
In this paper, the electromagnetic simulation of graphene is introduced and some widely used computational electromagnetic methods for graphene modeling are reviewed. The advantages and drawbacks of each method are discussed and numerical examples of these methods are given to illustrate their performance and application.

Electromagnetic Simulation of Graphene
To investigate the electromagnetic property of graphene, Maxwell equations need to be solved to simulate the wave propagation in graphene. In most cases, Maxwell equations have no analytical solutions except for a few simple canonical problems. Therefore, numerical methods are crucial to understand the wave guiding and scattering by graphene. There are three popular computational electromagnetic methods to simulate graphene: finite-difference time-domain method (FDTD), finite element method (FEM), and method of moment (MoM). Each method has its own advantages and disadvantages depending on the specific problem. Many commercial software packages based on these methods are available, such as FDTD Solutions, CST, HFSS, COMSOL, IE3D, and FEKO. All the software packages can be used to model graphene; however, none of these are developed specifically for graphene which makes them less efficient in dealing with this monoatomic layer device.
Graphene is modeled as a thin surface with complex conductivity which can be tuned by external electric field bias or chemical doping. The surface conductivity of graphene is commonly described by Kubo formula [32]: where Γ is the phenomenological scattering rate, is the electron charge, is the energy state, and ℎ denotes the reduced Plank constant. The first term of (1) is corresponding to the intraband contribution while the second term is due to the interband contribution. That is, The intraband term in (2) can be evaluated by and the interband term in (2) is approximated by where is the chemical potential, is the temperature, and is the Boltzmann constant.
In electromagnetic simulation, the 2D graphene surface is usually approximated by a 3D dielectric slab whose 3D conductivity is evaluated by the following equation: where is the thickness of the graphene layer. The accuracy of 3D dielectric slab model degrades as increases. It can be shown that, using the 2D conductivity, the SPP wavelength is [25] The 3D dielectric slab can support even and odd modes. The odd mode has the wavelength of and (6) can be approximated by (7) only if the following three conditions are satisfied [25]: (8) Figure 1 shows the relative error of using the dielectric slab model for graphene with respect to normalized and ; more details can be found in [25].
International Journal of Antennas and Propagation

Finite-Difference Time-Domain Method
FDTD is time-domain method which has specific advantage in transient problem and wideband problem. The broadband results can be obtained in one simulation. Compared to frequency-domain method, it takes less computing memory and simulation time. Additionally, FDTD is an iterative scheme which eliminates solving large linear system; hence, it is simple, robust, and easy to implement.
There are three general approaches to model graphene in FDTD: (I) using standard FDTD method with very fine mesh within the graphene sheet; (II) using the subcell FDTD method; (III) using surface resistive boundary condition.
In the first approach, graphene is modeled as thin layer with finite thickness and the surface conductivity is transformed to volume permittivity. According to the Courant-Friedrichs-Lewy (CFL) stability condition of FDTD, fine mesh leads to small time step, so this approach is both memory-and time-consuming which reduces its efficiency. To alleviate the time step constraint, some unconditionally stable FDTD methods are developed for graphene simulation. In [33], the locally one-dimensional (LOD) FDTD is applied for the efficient simulation of graphene device. They demonstrated that the LOD-FDTD can be 60% faster than standard FDTD with reasonable accuracy. Wang et al. [26] proposed an unconditionally stable one-step leapfrog alternating-direction-implicit (ADI) FDTD to study the SPPs in graphene structure. In Figure 2, they showed the SPPs propagating along the spiral waveguide. The one-step ADI-FDTD method gives accurate result which is comparable to conventional FDTD. They also showed that the conventional FDTD took 1.7 × 10 5 s for this problem, while ADI-FDTD only took 0.8 × 10 5 s with CFLN = 10.
In the second approach, graphene is treated as thin layer occupying a fraction of FDTD cell [27]. The Yee cell for the subcell FDTD method is shown in Figure 3, where component is split into , and , in cells occupying graphene [27]. However, this scheme is complex in mathematics and needs special type of PML to model infinitely thin sheets [34].
In the third approach, graphene is modeled as a zerothickness conductive sheet over which the fields satisfy the surface boundary condition. Nayyeri et al. [28,35] presented this scheme in detail. As shown in Figure 4(a), a conductive sheet locates at + 1/2. The tangential component of field is discontinuous over the sheet, so is split into 1 and 2 , which satisfy the following discretized Faraday law: The surface boundary condition is from which ( + 1/2) can be written as Substituting (11) into (9a) and (9b), 1 and 2 can be derived as The expressions of 1 , 2 , 1 , and 2 are shown in [28]. Figure 4(b) is a 3D FDTD cell containing conductive sheet at + 1/2. In this case, tangential components of magnetic field ( , ) and normal component of electric field ( ) are discontinuous over the sheet, so they are split into two parts on the two sides of the surface. Their updating equations can be derived in a similar way to 1D case and are presented in [28]. The surface boundary condition approach has been validated [28]. Figure 5(a) shows a 2D problem in which an infinite line source radiates above an infinite graphene sheet. The pattern of at the wavelength of 100 m is shown in Figure 5(b). A strong agreement between the proposed FDTD result and that derived from semianalytic method is obtained. In addition, unlike the subcell FDTD method, the classical PML is applicable to the surface boundary condition scheme.
In time-domain simulation, the conductivities of graphene shown in (2)-(4) need to be converted from frequency domain into time domain. The intraband conductivity, which has a Drude-like expression, can be easily converted into time domain. However, the interband conductivity has a complex logarithmic form which needs a vector fitting technique. A summation of partial fractions in terms of complex conjugate pole-residues is applied to approximate the conductivity. The detailed description of the fitting technique can be found in [35][36][37].

Finite Element Method
FEM is a full-wave numerical technique for electromagnetic boundary-value problems. The basic principle of this method is to discretize the whole computing domain with a finite number of subdomains in which the unknown function is expanded by simple interpolation functions with unknown coefficients. Then, a system of algebraic equations of these unknown coefficients is obtained by using Ritz variational or Galerkin's method. Finally, by solving this linear system, the approximated solution of the entire domain can be obtained. In FDTD, the computing space is discretized by orthogonal grid; this staircase approximation will reduce the modeling fidelity when it comes to complex geometry, while in FEM, where the triangles or tetrahedral elements are applied, arbitrary geometries can be modeled accurately. As a result, FEM has advantage in complex and inhomogeneous problems. Furthermore, FEM is a frequency domain method which makes it efficient in dealing with narrow-band problems.
Brar et al. [12] solved Maxwell's equation by FEM and modeled graphene as a thin sheet with 0.1 nm thickness; the results suggested that graphene can increase light-matter interactions at infrared energies. Software packages such as HFSS and COMSOL are based on FEM, and they are used widely in graphene simulation. Recently, we used COMSOL to investigate the transmission properties of a cylindrical graphene plasmon waveguide [21]. Tamagnone et al. [29] simulated a reconfigurable graphene antenna with HFSS; the structure of the antenna is shown in Figure 6(a); the input impedance of the antenna can be tuned by changing the chemical potential as shown in Figure 6(b).

Method of Moments
In contrast to FDTD and FEM which solve differential equations, MoM is a technique used to solve electromagnetic surface or volume integral equations in the frequency domain. In MoM, the quantities of interest are not the fields but the electromagnetic sources (surface or volume current), so only the surface of the geometry needs to be discretized. The surface current is discretized into wire segments and/or surface patches. A linear system can be constructed by the method of weighted residuals and the results of the linear equations give the surface current. The far-field result can be derived from the surface current by Green's function. Because it only requires calculating the boundary values instead of the values throughout the space, MoM is highly efficient for electrically large objects and is widely used in solving radiation and scattering problems. However, when applied to complex inhomogeneous cases, it will be very computationally expensive and less efficient.
The analytical expressions of dyadic Green's function for graphene are derived in [38]. Shapoval et al. [39] proposed integral equations based on surface-impedance boundary  condition to analyze plane wave scattering and absorption by graphene-strip gratings. The method of moments for graphene nanoribbons was developed in [40][41][42][43][44], in which the issue of nonlocality of graphene conductivity was taken into account. Nonlocal effect arises from spatial dispersion of graphene which is nonnegligible when dealing with slow modes supported by graphene nanoribbons. The spatially dispersive intraband conductivity tensor was derived in [41].
Software packages such as IE3D and FEKO are based on MoM. IE3D is applied to simulate the microwave propagation in a coplanar waveguide over graphene from 40 MHz to 110 GHz [45]. Cabellos-Aparicio et al. [30] used FEKO to study the radiated power of a graphene plasmonic antenna fed by photoconductive source. The antenna structure is shown in Figure 7(a) and the radiated power with respect to frequency is shown in Figure 7(b). The detailed parameters of the photoconductive antenna can be found in [30].

Discontinuous Galerkin Time-Domain Method
The graphene involved problems are often multiscale. Graphene is monoatomic and its thickness is much smaller than wavelength, so it is an electrically fine structure. In contrast, the substrate belongs to electrically coarse structure because its dimension is much greater than wavelength. As we have mentioned before, time-domain methods have the advantage that the broadband characterization can be obtained with only a single simulation. However, FDTD method has a serious efficiency problem in multiscale simulation because it requires high discretization density to model electrically fine structure due to its cartesian grid. The finite element time-domain (FETD) method is capable of modeling complex and fine structures and achieving highorder accuracy with high-order basis functions. The major drawback is a global linear system of equations that needs to be solved at each time step. The multiscale problems usually contain a large number of unknowns; FETD will be very computationally expensive in this case. Discontinuous Galerkin time-domain (DGTD) method is promising in multiscale problems [46]. DGTD allows for domain decomposition. A multiscale structure is divided into several subdomains and each subdomain can be discretized separately. All operations in DGTD are local, so large global matrix is split into several smaller matrices. Unlike FETD, the matrices of DGTD are inverted and stored before time marching, and different time integration scheme can be used in different subdomain. Additionally, DGTD is naturally adapted to parallel computing. Recently, DGTD has been used in nanophotonics field and considered as a viable alternative to the well-established FDTD and FETD methods [47]. Li et al. [31] proposed DGTD method with resistive boundary condition for the electromagnetic analyzing of graphene. They also applied this method to the magnetized graphene from microwave to THz range, where the anisotropic and disperse surface conductivity is involved [48]. the results of DGTD agree very well with the theoretical data. Figure 8(b) is the normalized extinction cross section of a freestanding graphene patch; good consistency is achieved between DGTD and integral equation method. More numerical examples can be found in [31].

Conclusion
Due to its intriguing properties, graphene is used more and more widely in electronic and photonic community. This paper gives a brief review of application of graphene, as well as the computational methods which can be used to investigate its electromagnetic properties. Although a number of famous commercial software packages are available for graphene simulation, none of them are developed specifically for graphene; therefore, they suffer from inefficiency. Each computational method has its own advantages and drawbacks; one should choose the appropriate method according to the specific problem. FDTD is simple and easy to implement; however, it is less efficient in modeling fine structure. FEM is flexible in geometry modeling but solving large global system makes it computationally expensive. DGTD is very suitable for multiscale problem and we believe it will be more and more widely used in modeling and designing graphene-assisted device.