Marching-on-in-Degree Time-Domain Integral Equation Solver for Transient Electromagnetic Analysis of Graphene

: The marching-on-in-degree (MOD) time-domain integral equation (TDIE) solver for the transient electromagnetic scattering of the graphene is presented in this paper. Graphene’s dispersive surface impedance is approximated using rational function expressions of complex conjugate pole-residue pairs with the vector ﬁtting (VF) method. Enforcing the surface impedance boundary condition, TDIE is established and solved in the MOD scheme, where the temporal surface impedance is carefully convoluted with the current. Unconditionally stable transient solution in time domain can be ensured. Wide frequency band information is obtained after the Fourier transform of the time domain solution. Numerical results validate the proposed method.

The TDIE method has been more and more popular in solving transient electromagnetic problems. The marching-on-in-degree (MOD) scheme [19][20][21][22][23][24][25] uses weighted Laguerre polynomials (WLPs) as the temporal basis functions, and involves no late-time instability, which may be encountered in the marching-on-in-time (MOT) scheme [26][27][28]. Generally, when handling resonant structures with a long tail current, many more degrees of WLPs are required. In this case, the frequency domain methods are recommended, but they may suffer from the ill-conditioned impedance matrix as well. When handling the potential internal resonance problems, the MOD scheme [19] and augmented electric field integral equation [22] might be the remedy.
In this paper, the MOD TDIE solver for the transient electromagnetic scattering of the graphene sheet is developed. The TDIE is based on the surface impedance boundary condition. The surface impedance is approximated by a vector fitting (VF) method, which has since its first introduction in 1999 widely applied for its robustness and efficiency [11][12][13][14][15][16][17][18][29][30][31][32][33][34]. Convolution of the temporal surface impedance and current is derived from properties of Laguerre polynomials.
In the next section, the formulation of VF and MOD TDIE is deduced. Section 3 presents numerical results to demonstrate the proposed method. Section 4 delineates conclusions.

Modeling of Graphene and VF Method
The Kubo formula-based graphene sheet dispersive frequency domain intraband and interband conductivity are as follows: where ω is the angular frequency, µ c is the chemical potential, Γ is the phenomenological scattering rate (or relaxation time), T is the temperature, e is the electron charge, k B is the Boltzmann constant, and h is the Planck constant. VF computes a rational approximation from the frequency domain data, using expressions in terms of complex conjugate pole-residue pairs. Graphene's frequency domain surface impedance is as follows: where a l and c l are poles and residues, respectively. The real parts of a l should be negative for causality and stability. Calculating ρ(ω) at several frequencies with a set of starting poles presents the over-determined linear problem, which can be solved with the least squares method. It is repeated using new poles as starting poles in an iterative procedure until the convergence is achieved [29][30][31].
Graphene's time domain surface impedance is [17] F −1 ( ) denotes inverse Fourier transform, and u(t) is the unit step function.

TDIE and Temporal Convolution
Enforcing the surface impedance boundary condition, the TDIE is | tan stands for tangential electric fields, S' is the surface of graphene sheet illuminated by the transient electric field E inc (r,t), r and r' are field and source point position vectors, R = |r − r |, J is the first derivative of current J with respect to time t, * is the temporal convolution, and c, ε 0 , and µ 0 are light velocity, permittivity and permeability in free space, respectively. The j-th degree Laguerre polynomial is The weighted Laguerre polynomial (WLP) is s is the temporal scaling factor, and ϕ j (t) converges to zero. The highest degree of WLP is t 0 is the time delay of the signal, and B is the bandwidth of the incident pulse. Given real (a i ) < 0, we have [23] N L is the highest degree of WLP, b l = a l /s − 1/2, and a l and c l are poles and residues in the VF method, respectively. In the computer model, the infinitely thin graphene sheet is built as a plane (zero-thickness) with a triangle mesh. The induced current of graphene is expressed by both spatial and temporal basis functions as N S is the number of inner edges of the triangles, J n,j is the unknown coefficient, and f n (r) is the RWG basis function defined by the triangle mesh. Using WLPs as temporal basis functions eliminates the late-time oscillation. With Equations (9) and (10), where Together with the property we get Coatings 2017, 7, 170

MOD Scheme
If the temporal derivative and integration terms [20] of the coefficients J n,j and Equation (16) is substituted into Equation (5), we get where Following the Galerkin's method with spatial testing functions f m (r), Equation (17) can be rewritten as where If a temporal testing procedure with ϕ i (t) is performed, we obtain the matrix equation of the i-th degree: where Coatings 2017, 7, 170 and R mn is the distance between centroids of the two triangles related to f m (r) and f n (r'). The matrix equation is solved degree by degree (0 ≤ i ≤ N L ), and the surface current can be obtained by Equations (10) and (28). After that, electromagnetic parameters, e.g., far field scattering, can be computed from the current.

Results and Discussion
Consider a 1 mm × 1 mm graphene sheet of µ c = 0.01 eV and Γ = 5 × 10 12 s −1 at T = 300 K. Following Equations (1)-(4), poles a l and residues c l of graphene's surface impedance computed with VF when p = 2 are listed in Table 1. Negative real parts of a l are required to ensure casualty and stability. Graphene's surface impedances are computed and compared in Figure 1, where the results of VF are identical to those computed by the Kubo formula.  (30) and Rmn is the distance between centroids of the two triangles related to fm(r) and fn(r').
The matrix equation is solved degree by degree (0 ≤ i ≤ NL), and the surface current can be obtained by Equations (10) and (28). After that, electromagnetic parameters, e.g., far field scattering, can be computed from the current.

Results and Discussion
Consider a 1 mm × 1 mm graphene sheet of μc = 0.01 eV and Γ = 5 × 10 12 s −1 at T = 300 K. Following Equations (1)-(4), poles al and residues cl of graphene's surface impedance computed with VF when p = 2 are listed in Table 1. Negative real parts of al are required to ensure casualty and stability. Graphene's surface impedances are computed and compared in Figure 1, where the results of VF are identical to those computed by the Kubo formula.  The graphene sheet is illuminated by the modulated Gaussian pulse The central frequency f0 and frequency bandwidth fbw are 0.204 THz and 0.408 THz, respectively, and ˆc t     rk / , k is the unit vector of the incidence direction and along the -z direction in this example, tp = 3.5η, η ==6/(2πfbw).
The transient electromagnetic problem is solved using the MOD TDIE solver presented in Section 2. The graphene surface current is expressed by 313 RWGs and 80 WLPs. The scaling factor s is 9.0 × 10 11 .
Two end points of a randomly chosen inner edge are (−0.413391 × 10 −3 , −0.254972 × 10 −3 , 0) and (−0.326874 × 10 −3 , −0.205698 × 10 −3 , 0). Current across this inner edge is shown in Figure 2. The dispersion of graphene leads to differences in comparison with the perfect electric conductor (PEC) plate of the same size. The unit of the lateral axis is light meter (lm). One light meter is the time that the EM wave propagates 1 m in free space. The current in late time is stable and converges to zero.  The graphene sheet is illuminated by the modulated Gaussian pulse The central frequency f 0 and frequency bandwidth f bw are 0.204 THz and 0.408 THz, respectively, and τ = t − r ×k/c,k is the unit vector of the incidence direction and along the -z direction in this example, t p = 3.5η, η = 6/(2πf bw ).
The transient electromagnetic problem is solved using the MOD TDIE solver presented in Section 2. The graphene surface current is expressed by 313 RWGs and 80 WLPs. The scaling factor s is 9.0 × 10 11 .
Two end points of a randomly chosen inner edge are (−0.413391 × 10 −3 , −0.254972 × 10 −3 , 0) and (−0.326874 × 10 −3 , −0.205698 × 10 −3 , 0). Current across this inner edge is shown in Figure 2. The dispersion of graphene leads to differences in comparison with the perfect electric conductor (PEC) plate of the same size. The unit of the lateral axis is light meter (lm). One light meter is the time that the EM wave propagates 1 m in free space. The current in late time is stable and converges to zero. The graphene sheet is illuminated by the modulated Gaussian pulse The central frequency f0 and frequency bandwidth fbw are 0.204 THz and 0.408 THz, respectively, and ˆc t τ = − × r k / , k is the unit vector of the incidence direction and along the -z direction in this example, tp = 3.5η, η ==6/(2πfbw). The transient electromagnetic problem is solved using the MOD TDIE solver presented in Section 2. The graphene surface current is expressed by 313 RWGs and 80 WLPs. The scaling factor s is 9.0 × 10 11 .
Two end points of a randomly chosen inner edge are (−0.413391 × 10 −3 , −0.254972 × 10 −3 , 0) and (−0.326874 × 10 −3 , −0.205698 × 10 −3 , 0). Current across this inner edge is shown in Figure 2. The dispersion of graphene leads to differences in comparison with the perfect electric conductor (PEC) plate of the same size. The unit of the lateral axis is light meter (lm). One light meter is the time that the EM wave propagates 1 m in free space. The current in late time is stable and converges to zero.   Wide frequency band information, e.g., bistatic radar cross section (RCS), is obtained after the Fourier transform. Examples are chosen at 0.102, 0.204, and 0.374 THz. The results show good agreement with those obtained via MOM in Figure 3. Another 12.5 µm × 25 µm (xoy plane) graphene sheet (µ c = 0.01 eV, Γ = 5 × 10 12 s −1 , T = 300 K) and frequency up to 12 THz are under discussion. The poles a l and residues c l of graphene's surface impedance computed with VF are listed in Table 2. Graphene's surface impedances are computed and compared in Figure 4, where the results of VF agree well with those computed by the Kubo formula. Another 12.5 μm × 25 μm (xoy plane) graphene sheet (μc = 0.01 eV, Γ = 5 × 10 12 s −1 , T = 300 K) and frequency up to 12 THz are under discussion. The poles al and residues cl of graphene's surface impedance computed with VF are listed in Table 2. Graphene's surface impedances are computed and compared in Figure 4, where the results of VF agree well with those computed by the Kubo formula.   For graphene, the total-scattering cross section (TSCS), the absorption cross section (ACS), and the extinction cross section (ECS) are other interesting parameters. The ECS at frequency f is computed according to the electromagnetic optical theorem [35] as:  For graphene, the total-scattering cross section (TSCS), the absorption cross section (ACS), and the extinction cross section (ECS) are other interesting parameters. The ECS at frequency f is computed according to the electromagnetic optical theorem [35] as: where E inc and E sca are the incident and scattering amplitude in the forward direction, and Im( ) is the imaginary part. The graphene surface current is expressed by 123 RWGs and 80 WLPs. The scaling factor s is 3.6 × 10 13 . The ECS above 1 THz computed with MOD TDIE show good agreement with those obtained via MOM in Figure 5. where E inc and E sca are the incident and scattering amplitude in the forward direction, and Im( ) is the imaginary part. The graphene surface current is expressed by 123 RWGs and 80 WLPs. The scaling factor s is 3.6 × 10 13 . The ECS above 1 THz computed with MOD TDIE show good agreement with those obtained via MOM in Figure 5. Based on the MOD TDIE solver developed in this paper, transient electromagnetic modeling of multilayer graphene with metallic and/or dielectric substrate will be studied with PMCHWT formulation in the future by the authors. The memory and time complexity of the standard MOD TDIE method are O (NS 2 NL) and O (NS 2 NL 2 ), respectively, where NS is the number of inner edges, NL is the highest degree of WLPs. The memory requirement and computation time can decrease by combining fast algorithms or parallelization codes.

Conclusions
TDIE for the transient electromagnetic scattering of the graphene is solved following the MOD scheme. TDIE is established enforcing the surface impedance boundary condition. The convolution of the temporal surface impedance and current is deduced. The temporal surface impedance is obtained after the inverse Fourier transform of the frequency domain counterparts, which are approximated using the VF method. Stable and accurate solution can be ensured. Based on the MOD TDIE solver developed in this paper, transient electromagnetic modeling of multilayer graphene with metallic and/or dielectric substrate will be studied with PMCHWT formulation in the future by the authors. The memory and time complexity of the standard MOD TDIE method are O (N S 2 N L ) and O (N S 2 N L 2 ), respectively, where N S is the number of inner edges, N L is the highest degree of WLPs. The memory requirement and computation time can decrease by combining fast algorithms or parallelization codes.

Conclusions
TDIE for the transient electromagnetic scattering of the graphene is solved following the MOD scheme. TDIE is established enforcing the surface impedance boundary condition. The convolution of the temporal surface impedance and current is deduced. The temporal surface impedance is obtained after the inverse Fourier transform of the frequency domain counterparts, which are approximated using the VF method. Stable and accurate solution can be ensured.