## Abstract

The one-step leapfrog alternating-direction-implicit finite-difference time-domain (ADI-FDTD) method is reformulated for simulating general electrically dispersive media. It models material dispersive properties with equivalent polarization currents. These currents are then solved with the auxiliary differential equation (ADE) and then incorporated into the one-step leapfrog ADI-FDTD method. The final equations are presented in the form similar to that of the conventional FDTD method but with second-order perturbation. The adapted method is then applied to characterize (a) electromagnetic wave propagation in a rectangular waveguide loaded with a magnetized plasma slab, (b) transmission coefficient of a plane wave normally incident on a monolayer graphene sheet biased by a magnetostatic field, and (c) surface plasmon polaritons (SPPs) propagation along a monolayer graphene sheet biased by an electrostatic field. The numerical results verify the stability, accuracy and computational efficiency of the proposed one-step leapfrog ADI-FDTD algorithm in comparison with analytical results and the results obtained with the other methods.

© 2013 Optical Society of America

## 1. Introduction

Due to its simplicity in directly discretizing Maxwell’s curl equations in time domain and its ability for obtaining a wideband response with a single run of simulation, the Finite-Difference Time-Domain method (FDTD) has been widely used for modeling various isotropic and anisotropic dispersive media and structures for RF to Optical applications [1–10]. In general, the FDTD schemes for modeling dispersive media can be divided into several categories: auxiliary differential equation (ADE) method [11], *Z*-transform technique [12], and discrete convolution of the dispersion relation (DC-DR) method [13]. Among them, the ADE method is simple and applicable for handling linear and nonlinear dispersive media, and its capability can be enhanced greatly by using the complex- conjugate pole-residue pair (CC-PR), critical points (CP), and rational-fraction dispersion (RFM) techniques developed recently [14–16]. Different from the ADE method, the DC-DR method proposed in [17,18] can model all dispersive media with a single general formulation. It is based on the discretization of the dispersion relation with the convolutions that include current density convolution (JEC), recursive convolution (RC), and piecewise linear recursive convolution (PLRC) [19,20]. However, their numerical implementations are more complex than that of the ADE with little accuracy improvement.

Many FDTD methods were developed in the past few years, such as EJ collocated FDTD [21], Runge–Kutta exponential time differencing formulation (RKETD) FDTD [22], matrix exponential (ME) FDTD [23], and exponential time differencing (ETD) FDTD [24], *etc*. Unfortunately, they cannot be easily used to model the general dispersive media.

On the other hand, in order to overcome the Courant–Friedrich–Levy (CFL) of the conventional FDTD method, one-step leapfrog ADI-FDTD method has been developed [25]. It is originated from the conventional ADI-FDTD but with no mid-time computation needed [26]. Therefore, it is similar to the conventional FDTD and has higher computational efficiency than the other unconditionally stable FDTD methods [27]. Recently, both convolutional perfectly matched layer (CPML) and Mur absorbing boundary condition for the leapfrog ADI-FDTD have been developed to simulate open structure problems [28,29]. While in [30], the leapfrog ADI-FDTD method for periodic structures was presented. More recently, the divergence properties of the leapfrog ADI-FDTD was studied in [31], and it has been further extended to simulating the lossy media [32]. In addition, the leapfrog ADI-FDTD has also been adapted to modeling magnetized plasma [33].

In this paper, we further extend the leapfrog ADI-FDTD method to simulating the general dispersive media based on the newly developed ADE method presented in [17,18]. Its successful applications are demonstrated by simulating the magnetized plasma, monolayer graphene sheet biased by a magnetostatic field, and surface plasmon polaritons (SPPs) propagating along the graphene sheet biased by an electrostatic field. Evidently, it can also be employed to simulate various graphene-based nonreciprocal RF, THz nanoelectonic and nanophotonic devices and components.

## 2. The proposed leapfrog ADI-FDTD method for general dispersive media

For a generalized electrically anisotropic media, the Ampere’s law in frequency domain can be written as

Where $E={\left[{E}_{x},{E}_{y}\text{,}{E}_{z}\right]}^{T}$, $H={\left[{H}_{x}\text{,}{H}_{y}\text{,}{H}_{z}\right]}^{T}$, ${\epsilon}_{0}$is the permittivity in free space, and ${\epsilon}_{re}$ is the relative permittivity tensor,*i.e.*

*x*-component of ${\epsilon}_{0}{\epsilon}_{re}(\omega )\cdot E(\omega )$can be written as

*x*-component of the polarization current density, andSimilarly, we also have Therefore, Eq. (1a) can be expressed asConverting Eq. (5) into the time domain by using$j\omega \to \partial /\partial t$, we haveWhere $J={\left[{J}_{x}{J}_{y}{J}_{z}\right]}^{T}$,

Similarly, the Faraday’s induction law can also be rewritten as

where*μ*is the permeability in free space.

_{0}By applying the ADI process to Eq. (6a) and Eq. (7), we can obtain the ADI-FDTD formulations for the dispersive media as:

For the sub-step of *n to n* + 1/2,

For the sub-step of *n* + 1/2 to *n* + 1,

*p*

_{1}and

*p*

_{2}are the time index within single time step, and they are within the ranges of [0, 1/2] and [1/2, 1], respectively. With the algebraic manipulations used in [25], the one-step leapfrog ADI-FDTD method with the polarization current density can be derived as follows.

First, by substituting Eq. (8b) into Eq. (8a), we have

*n*with

*n*-1, we obtain

Similarly, an iterative equation for the magnetic field can be derived as follows. From Eq. (8a), it has

By substituting Eq. (12) into Eq. (8b), we obtainIn practice, it is not easy to compute the cross term $\Delta {t}^{2}/4{\mu}_{0}{\epsilon}_{0}A({J}^{n+{p}_{2}}-{J}^{n+{p}_{1}})$ in Eq. (15). To simplify it, we set *p*_{1} = *p*_{2} = 1/2 as

The above derivations of Eqs. (16a) and (16b) appear quite tedious. In fact, both of them can be reformulated in the form similar to that of the explicit FDTD method by two auxiliary variables **e** and **h** such that

The above leapfrog ADI-FDTD formulas are applicable for various dispersive media as demonstrated below.

#### 2.1 Lossy media

A lossy medium can be considered as a simple dispersive medium as studied in [32]. For simplicity, we choose its conductivity and equivalent permittivity as

It should be noted that Eq. (20) and Eq. (16b) have the same form as those in [32], and its unconditionally stable property has been studied in [34].

#### 2.2 Magnetized plasma

Assume that the applied biasing magnetic field is in the *z*-axis direction, we then have [18]

*ω*is the plasma frequency,

_{p}*ω*is the cyclotron frequency, and

_{b}*ν*is the electron collision rate. By substituting Eqs. (21a)-(21d) into Eq. (4), we obtain

_{c}*i.e. ${J}_{x}^{n}={J}_{xx}^{n}+{J}_{xy}^{n}$*. Both ${J}_{xx}(\omega )$ and ${J}_{xy}(\omega )$ have the general form as

*n*-1, we obtain

*J*is derived as

#### 2.3 Graphene biased by a magnetostatic field

The monolayer graphene is biased by a magnetostatic field **B _{o}** in the

*z*-axis direction, and it is placed in the

*x*-

*y*plane. Then, its volume conductivity$\sigma $can be calculated from its surface conductivity ${\sigma}_{s}$by using$\sigma ={\sigma}_{s}/\Delta z$. In the frequency up to 10 THz, its ${\sigma}_{s}$can be described by a 2-D Drude model [35], with${\epsilon}_{re}(\omega )=I+\sigma /(j\omega {\epsilon}_{0})$, and

*v*is the Fermi velocity,

_{F}*k*is the Boltzmann constant, ${\mu}_{c}$ is the chemical potential,

_{B}*T*is the operating temperature, and $\hslash $ is the reduced Planck’s constant.

According to Eqs. (27a)-(27d), it is found that the${\epsilon}_{re}(\omega )$of graphene sheet can be obtained by replacing the parameters of magnetized plasma in Eq. (21) with${\omega}_{p}^{2}$$\to $${\sigma}_{0}$, ${\nu}_{c}^{}$$\to $$\nu $, ${\omega}_{b}$$\to $${\omega}_{c}$, and${\epsilon}_{z}(\omega )=0$. Therefore, the magnetized graphene can also be simulated using Eq. (16a), Eq. (16b), and Eq. (26).

#### 2.4 Graphene biased by an electrostatic field

If the monolayer graphene sheet is only biased by an electrostatic field${E}_{0}$, its relative permittivity can be obtained by setting ${\omega}_{c}=0$ in Eq. (27), *i.e.*

*J*and

_{x}*J*can be obtained easily.

_{y}In our simulation below, the chemical potential ${\mu}_{c}$of graphene is calculated by

## 3. Numerical results and discussion

To demonstrate the capability of our developed leapfrog ADI-FDTD algorithm for treating general dispersive media, three numerical experiments were performed: electromagnetic wave propagation in a PEC waveguide loaded with a magnetized plasma slab [33], electromagnetic wave transmission through a magnetized two-dimensional monolayer graphene sheet, and the surface plasmon polaritons along a two-dimensional graphene sheet.

#### 3.1 The PEC Waveguide Loaded with Magnetized Plasma

Figure 1 shows the rectangular waveguide of 40mm × 40mm, inserted with a magnetized plasma slab of 5 mm in width. The waveguide is terminated by the 10-layer CPMLs at both ends [33]. In such a structure, both propagating and evanescent modes exist.

In our simulation, we set ∆x = ∆y = ∆z = ∆ = 1mm. A modulated Gaussian pulse of$\mathrm{cos}(2\pi {f}_{0}t)\mathrm{exp}[-4\pi {(t-{t}_{0})}^{2}/{\tau}^{2}]$was used to excite electromagnetic field, where *f*_{0} = 10 GHz,$\Delta t=\Delta /\sqrt{3}{c}_{0}\cdot \text{CFLN}$,$\tau =160\Delta t/\text{CFLN}$, and ${t}_{0}=\tau $. CFLN is the CFL number which measures the ratio of the time step to the CFL limit. The TE_{10}-mode was excited and its cutoff frequency is 3.75GHz. The plasma frequency *ω _{p}* = 2π × 10 Grad/s,

*ν*= 10 GHz, and

*ω*= 10 Grad/s. An observation point is set at the waveguide center, and it is ten cells away from the magnetized plasma slab. Figure 2 shows the recorded

_{b}*E*-component at the observation point verses time with different CFLNs. It is seen that for CFLN = 1, the recorded

_{z}*E*–component agrees well with that of the conventional FDTD. As the CFLN increases from 1 to 3, the errors are increased slightly, but they are in an acceptable range. To check the stability of our developed algorithm, we run the simulation to one million time steps with the time step CFLN = 2 and 3, respectively; the recorded field was found to be always stable.

_{z}#### 3.2 Magnetized graphene sheet

When a graphene sheet is biased by a magnetostatic field **B**_{0}, it produces gyrotropy [35]. By applying the proposed method, we can predict the transmission properties of an infinite graphene sheet at both THz and optical bands.

In our simulation, a two-dimensional graphene sheet was placed at the center of *z*-direction in the *x*-*y* plane as in Fig. 3(a). We chose the maximum frequency *f _{max}* =

*c*

_{0}/

*λ*= 2 THz, the FDTD cells ∆x = ∆y = ∆z = ∆ =

_{min}*λ*/20,

_{min}*ν*= 20THz,

*T*= 300K,

*μ*= 0.1 eV, and $\Delta t=\Delta /(\sqrt{3}{c}_{0})$. The overall computational domain consists of 40 × 40 × 180 cells, including 10-layer CPML [28]. The normally incident plane wave is polarized in the

_{c}*x*-direction and described by

*x*-

*y*plane and 4 cells away from the graphene sheet. The analytical transmission coefficient of the infinite graphene sheet can be calculated as [35]:

*FFT*represents the fast Fourier transformation.

Figure 3(b) shows the computed transmission coefficient using the leapfrog ADI-FDTD with B_{0} = 1T and CFLN = 1, 2, and 3, respectively. The analytical result is also plotted for comparison. It is seen that, when CFLN = 1, the transmission coefficient |T| and the analytical solution agrees very well up to 2 THz. When CFLN = 2, the value of |T| still agrees well with the analytical solution at low frequencies. When CFLN = 3, the errors become much larger at high frequencies. This means that the leapfrog ADI-FDTD is quite suitable for computing the transmission coefficient at low frequencies.

Figures 4(a) and 4(b) show the snapshots of *E*_{y} – component at time *t* = 2.5*t*_{0} for CFLN = 1 and 3, respectively. The field is recorded at 10 cells away from the CPML in the *y*-*z* plane. It is seen that the field distribution computed with CFLN = 3 is slightly different from that of CFLN = 1. On the other hand, it is noted that, as the incident EM wave is polarized in the *x*-axis direction, the *E _{y}* – component remains zero before it interacts with the graphene sheet. Therefore, Fig. 4 shows the gyrotropic property of the magnetized graphene sheet.

#### 3.3 SPPs propagating along the graphene sheet

It is known that when$\mathrm{Re}[{\epsilon}_{d}(\omega )]<0$, the graphene sheet acts as a thin metal film, and a TM-mode SPPs can be excited. The properties of SPPs have been studied using the finite-difference frequency-domain (FDFD) and the conventional FDTD [36,37], respectively. However, as it is an atomically thin structure, both methods may take very long simulation time because of the small cell sizes and small time steps used. Here, we fast solve it using the proposed method with a large CFLN.

In our simulation, we chose ∆y = ∆z = ∆ = 50 nm, *ν* = 1 THz, *T* = 300 K, *μ _{c}* = 0.15 eV,and $\Delta t=\Delta /(\sqrt{3}{c}_{0})$. The overall computational domain, as shown in Fig. 5, consists of 820 × 180 cells, including 10-layer CPML [31]. A single point source, placed one cell over the graphene sheet and one cell away from the CPML, is used to excite the SPPs, and the exciting magnetic current source is${J}_{m,x}^{}(t)=\mathrm{sin}(2\pi {f}_{0}t)$with

*f*

_{0}= 6 THz. Therefore, the grid resolution is

*λ*

_{0}/∆ =

*c*

_{0}/(

*f*

_{0}∆) = 1000, and for such a fine structure the large CFLN is needed so as to reduce the simulation time.

Figures 6(a)-6(d) show snapshots of the *E _{z}* – and

*E*–components at time

_{y}*t*= 2.887 ps. In Figs. 6(a) and 6(b), the numerical results were obtained with the leapfrog ADI-FDTD with CFLN = 10, and in Figs. 6(c) and 6(d), they were computed using the conventional ADE-FDTD [17] with CFLN = 1. It is seen that the results obtained with the leapfrog ADI-FDTD agree well with that of the conventional FDTD. Its computational time is about 46 seconds for 3000 time steps with CFLN = 10, while the conventional FDTD uses 142 seconds for 30000 time steps. Our simulation platform is the Lenovo PC with Intel Dual Core P8700 of 2.53 GHz and RAM of 8GB.

## 4. Conclusion

In this paper, the one-step leapfrog ADI-FDTD method is adapted for simulating general dispersive media. In the method, dispersive properties are modelled with equivalent polarization currents which are then solved by the newly developed ADE formulations. Numerical experiments have shown that the proposed method can simulate the conventional magnetized plasma as well as magnetized graphene sheet accurately. In particular, it can be used to simulate the SPPs on a monolayer graphene sheet and the simulation time was reduced greatly without sacrificing the accuracy.

## Acknowledgments

This work was supported in part by the NSFC under Grant 61171037, the Specialized Research Fund for the Doctoral Program of Higher Education under Grant 20100101110065, the State Key Lab of MOI of Zhejiang University, the National Basic Research Program under Grant 2009CB320204 of China, and the Natural Science and Engineering Research Council of Canada through its Discovery Grant program.

## References and links

**1. **A. Taflove and S. C. Hagness, *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, 3rd ed. (Artech House, 2005).

**2. **C. Lundgren, R. Lopez, J. Redwing, and K. Melde, “FDTD modeling of solar energy absorption in silicon branched nanowires,” Opt. Express **21**(9), 329–400 (2013). [PubMed]

**3. **E. H. Khoo, I. Ahmed, R. S. M. Goh, K. H. Lee, T. G. G. Hung, and E. P. Li, “Efficient analysis of mode profiles in elliptical microcavity using dynamic-thermal electron-quantum medium FDTD method,” Opt. Express **21**(5), 5910–5923 (2013). [CrossRef] [PubMed]

**4. **D. T. Nguyen and R. A. Norwood, “Label-free, single-object sensing with a microring resonator: FDTD simulation,” Opt. Express **21**(1), 49–59 (2013). [CrossRef] [PubMed]

**5. **I. R. Çapoğlu, A. Taflove, and V. Backman, “Computation of tightly-focused laser beams in the FDTD method,” Opt. Express **21**(1), 87–101 (2013). [CrossRef] [PubMed]

**6. **S. Buil, J. Laverdant, B. Berini, P. Maso, J. P. Hermier, and X. Quélin, “FDTD simulations of localization and enhancements on fractal plasmonics nanostructures,” Opt. Express **20**(11), 11968–11975 (2012). [CrossRef] [PubMed]

**7. **C. M. Dissanayake, M. Premaratne, I. D. Rukhlenko, and G. P. Agrawal, “FDTD modeling of anisotropic nonlinear optical phenomena in silicon waveguides,” Opt. Express **18**(20), 21427–21448 (2010). [CrossRef] [PubMed]

**8. **G. Singh, K. Ravi, Q. Wang, and S. T. Ho, “Complex-envelope alternating-direction-implicit FDTD method for simulating active photonic devices with semiconductor/solid-state media,” Opt. Lett. **37**(12), 2361–2363 (2012). [CrossRef] [PubMed]

**9. **G. Singh, E. L. Tan, and Z. N. Chen, “Split-step finite-difference time-domain method with perfectly matched layers for efficient analysis of two-dimensional photonic crystals with anisotropic media,” Opt. Lett. **37**(3), 326–328 (2012). [CrossRef] [PubMed]

**10. **S. Zhao, “High-order FDTD methods for transverse electromagnetic systems in dispersive inhomogeneous media,” Opt. Lett. **36**(16), 3245–3247 (2011). [CrossRef] [PubMed]

**11. **W. H. Weedon and C. M. Rappaport, “A general method for FDTD modeling of wave propagation in arbitrary frequency-dispersive media,” IEEE Trans. Antenn. Propag. **45**(3), 401–410 (1997). [CrossRef]

**12. **D. Sullivan, “Nonlinear FDTD formulations using Z transforms,” IEEE Trans. Microw. Theory Tech. **43**(3), 676–682 (1995). [CrossRef]

**13. **R. J. Luebbers and F. Hunsberger, “FDTD *N*th-order dispersive media,” IEEE Trans. Antenn. Propag. **40**(11), 1297–1301 (1992). [CrossRef]

**14. **M. Han, R. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” IEEE Microw. Wirel. Compon. Lett. **16**(3), 119–121 (2006). [CrossRef]

**15. **P. G. Etchegoin, E. C. Le Ru, and M. Meyer, “An analytic model for the optical properties of gold,” J. Chem. Phys. **125**(16), 164705 (2006). [CrossRef] [PubMed]

**16. **L. Han, D. Zhou, K. Li, X. Li, and W. P. Huang, “A rational-fraction dispersion model for efficient simulation of dispersive material in FDTD method,” IEEE J. Light. Tech. **30**(13), 2216–2225 (2012). [CrossRef]

**17. **M. Alsunaidi and A. Al-Jabr, “A general ADE-FDTD algorithm for the simulation of dispersive structures,” IEEE Photon. Technol. Lett. **21**(12), 817–819 (2009). [CrossRef]

**18. **A. Al-Jabr, M. Alsunaidi, T. Khee, and B. S. Ooi, “A simple FDTD algorithm for simulating EM-wave propagation in general dispersive anisotropic material,” IEEE Trans. Antenn. Propag. **61**(3), 1321–1326 (2013). [CrossRef]

**19. **L. J. Xu and N. C. Yuan, “FDTD formulations for scattering from 3-D anisotropic magnetized plasma objects,” IEEE Antennas Wirel. Propag. Lett. **5**(1), 335–338 (2006). [CrossRef]

**20. **D. F. Kelley and R. J. Luebbers, “Piecewise linear recursive convolution for dispersive media using FDTD,” IEEE Trans. Antenn. Propag. **44**(6), 792–797 (1996). [CrossRef]

**21. **Y. Yu and J. Simpson, “An EJ collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma,” IEEE Trans. Antenn. Propag. **58**(2), 469–478 (2010). [CrossRef]

**22. **S. Liu and S. Liu, “Runge-Kutta exponential time differencing FDTD method for anisotropic magnetized plasma,” IEEE Antennas Wirel. Propag. Lett. **7**, 306–309 (2008). [CrossRef]

**23. **D. Y. Heh and E. L. Tan, “FDTD modeling for dispersive media using matrix exponential method,” IEEE Microw. Wirel. Compon. Lett. **19**(2), 53–55 (2009). [CrossRef]

**24. **S. Huang and F. Li, “FDTD implementation for magnetoplasma medium using exponential time differencing,” IEEE Microw. Wirel. Compon. Lett. **15**(3), 183–185 (2005). [CrossRef]

**25. **S. J. Cooke, M. Botton, T. M. Antonsen, and B. Levush, “A leapfrog formulation of the 3D ADI-FDTD algorithm,” Int. J. Numer. Model. **22**(2), 187–200 (2009). [CrossRef]

**26. **F. Zheng, Z. Chen, and J. Zhang, “Toward the development of a three dimensional unconditionally stable finite-different time-domain method,” IEEE Trans. Microw. Theory Tech. **48**(9), 1550–1558 (2000). [CrossRef]

**27. **E. L. Tan, “Fundamental schemes for efficient unconditionally stable implicit finite-difference time-domain methods,” IEEE Trans. Antenn. Propag. **56**(1), 170–177 (2008). [CrossRef]

**28. **X. H. Wang, W. Y. Yin, Y. Q. Yu, Z. Chen, J. Wang, and Y. Guo, “A convolutional perfect matched layer (CPML) for one-step leapfrog ADI-FDTD method and its applications to EMC problems,” IEEE Trans. Electromagn. Compat. **54**(5), 1066–1076 (2012). [CrossRef]

**29. **T. H. Gan and E. T. Tan, “Mur absorbing boundary condition for 2-D leapfrog ADI-FDTD method,” in *proceedings of IEEE Conference on Asia-Pacific Antennas and Propagation*, 3–4 (2012). [CrossRef]

**30. **Y. F. Mao, B. Chen, J. L. Xia, J. Chen, and J. Z. Tang, “Application of the leapfrog ADI-FDTD method to periodic structures,” IEEE Antennas Wirel. Propag. Lett. **12**, 599–602 (2013). [CrossRef]

**31. **T. H. Gan and E. L. Tan, “Analysis of the divergence properties for the three-dimensional leapfrog ADI-FDTD method,” IEEE Trans. Antenn. Propag. **60**(12), 5801–5808 (2012). [CrossRef]

**32. **T. H. Gan and E. L. Tan, “Unconditionally stable leapfrog ADI-FDTD method for lossy media,” Progress In Electromagnetics Research M **26**, 173–186 (2012).

**33. **X. H. Wang, W. Y. Yin, and Z. Chen, “One-step leapfrog ADI-FDTD method for anisotropic magnetized plasma,” in *proceedings of IEEE International Microwave Symposium*, 1–4 (2013).

**34. **J. Y. Gao and H. X. Zheng, “One-Step leapfrog ADI-FDTD method for lossy media and its stability analysis,” Progress In Electromagnetics Research Letters **40**, 49–60 (2013).

**35. **D. L. Sounas and C. Caloz, “Gyrotropy and nonreciprocity of graphene for microwave applications,” IEEE Trans. Microw. Theory Tech. **60**(4), 901–914 (2012). [CrossRef]

**36. **B. Wang, X. Zhang, X. Yuan, and J. Teng, “Optical coupling of surface plasmons between graphene sheets,” Appl. Phys. Lett. **100**(13), 131111 (2012). [CrossRef]

**37. **H. Lin, M. F. Pantoja, L. D. Angulo, J. Alvarez, R. G. Martin, and S. G. Garcia, “FDTD modeling of graphene devices using complex conjugate dispersion material model,” IEEE Microw. Wirel. Compon. Lett. **22**(12), 612–614 (2012). [CrossRef]