One-step leapfrog ADI-FDTD method for simulating electromagnetic wave propagation in general dispersive media

Xiang-Hua Wang; Wen-Yan Yin; Zhi Zhang (David) Chen

doi:10.1364/OE.21.020565

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

j ω ε_{0} ε_{r e} \cdot E (ω) = \nabla \times H (ω)

Where

E = {[E_{x}, E_{y}, E_{z}]}^{T}

,

H = {[H_{x}, H_{y}, H_{z}]}^{T}

,

ε_{0}

is the permittivity in free space, and

ε_{r e}

is the relative permittivity tensor, i.e.

ε_{r e} (ω) = [\begin{matrix} ε_{x x} (ω) & ε_{x y} (ω) & ε_{x z} (ω) \\ ε_{y x} (ω) & ε_{y y} (ω) & ε_{y z} (ω) \\ ε_{z x} (ω) & ε_{z y} (ω) & ε_{z z} (ω) \end{matrix}]

Therefore, the x-component of

ε_{0} ε_{r e} (ω) \cdot E (ω)

can be written as

\begin{array}{l} {[ε_{0} ε_{r e} (ω) \cdot E (ω)]}_{x} = ε_{0} E_{x} (ω) + ε_{0} (ε_{x x} - 1) E_{x} (ω) + ε_{0} ε_{x y} E_{y} (ω) + ε_{0} ε_{x z} E_{z} (ω) \\ = ε_{0} E_{x} (ω) + P_{x x} (ω) + P_{x y} (ω) + P_{x z} (ω) \end{array}

where

P_{x x} (ω)

,

P_{x y} (ω)

and

P_{x z} (ω)

are the polarization variables [18]. Then, by substituting Eq. (2) into Eq. (1a), we have

j ω ε_{0} E_{x} (ω) + j ω [P_{x x} (ω) + P_{x y} (ω) + P_{x z} (ω)] = j ω ε_{0} E_{x} (ω) + J_{x} (ω) = {[\nabla \times H (ω)]}_{x}

where

J_{x} (ω)

is the x-component of the polarization current density, and

J_{x} (ω) = j ω [P_{x x} (ω) + P_{x y} (ω) + P_{x z} (ω)]

Similarly, we also have

J_{y} (ω) = j ω [P_{y x} (ω) + P_{y y} (ω) + P_{y z} (ω)]

J_{z} (ω) = j ω [P_{z x} (ω) + P_{z y} (ω) + P_{z z} (ω)] .

Therefore, Eq. (1a) can be expressed as

j ω ε_{0} E (ω) = \nabla \times H (ω) - J (ω) .

Converting Eq. (5) into the time domain by using

j ω \to \partial / \partial t

, we have

d E / d t = [(A - B) H - J] / ε_{0}

Where

J = {[J_{x} J_{y} J_{z}]}^{T}

,

A = [\begin{matrix} 0 & 0 & \partial / \partial y \\ \partial / \partial z & 0 & 0 \\ 0 & \partial / \partial x & 0 \end{matrix}], B = [\begin{matrix} 0 & \partial / \partial z & 0 \\ 0 & 0 & \partial / \partial x \\ \partial / \partial y & 0 & 0 \end{matrix}] .

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

d H / d t = (B - A) E / μ_{0}

where μ₀ is the permeability in free space.

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,

E^{n + 1 / 2} = E^{n} + (0.5 Δ t / ε_{0}) (A H^{n + 1 / 2} - B H^{n} - J^{n + p_{1}})

H^{n + 1 / 2} = H^{n} + (0.5 Δ t / μ_{0}) (B E^{n + 1 / 2} - A E^{n})

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

E^{n + 1} = E^{n + 1 / 2} + (0.5 Δ t / ε_{0}) (A H^{n + 1 / 2} - B H^{n + 1} - J^{n + p_{2}})

H^{n + 1} = H^{n + 1 / 2} + (0.5 Δ t / μ_{0}) (B E^{n + 1 / 2} - A E^{n + 1})

where p₁ and p₂ 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

E^{n + 1 / 2} = E^{n} + \frac{Δ t}{2 ε_{0}} (A H^{n} + \frac{Δ t}{2 μ_{0}} A B E^{n + 1 / 2} - \frac{Δ t}{2 μ_{0}} A A E^{n} - B H^{n} - J^{n + p_{1}}) .

By substituting Eq. (8d) into Eq. (8c), and replacing n with n-1, we obtain

E^{n} = E^{n - 1 / 2} + \frac{Δ t}{2 ε_{0}} (A H^{n} - \frac{Δ t}{2 μ_{0}} A B E^{n + 1 / 2} + \frac{Δ t}{2 μ_{0}} A A E^{n} - B H^{n} - J^{n + p_{2} - 1}) .

Adding Eq. (9) and Eq. (10) on the both corresponding sides, we obtain the iterative equation for the electric field in the leapfrog ADI-FDTD method as

(I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) E^{n + 1 / 2} = (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) E^{n - 1 / 2} + \frac{Δ t}{ε_{0}} (A H^{n} - B H^{n}) - \frac{Δ t}{2 ε_{0}} (J^{n + p_{1}} + J^{n + p_{2} - 1}) .

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

E^{n} = E^{n + 1 / 2} - (0.5 Δ t / ε_{0}) (A H^{n + 1 / 2} - B H^{n} - J^{n + p_{1}}) .

By substituting Eq. (12) into Eq. (8b), we obtain

H^{n + 1 / 2} = H^{n} + \frac{Δ t}{2 μ_{0}} (B E^{n + 1 / 2} + \frac{Δ t}{2 ε_{0}} A A H^{n + 1 / 2} - \frac{Δ t}{2 ε_{0}} A B H^{n} - A E^{n + 1 / 2} - \frac{Δ t}{2 ε_{0}} A J^{n + p_{1}}) .

By substituting Eq. (8c) into Eq. (8d), we also obtain

H^{n + 1} = H^{n + 1 / 2} + \frac{Δ t}{2 μ_{0}} (B E^{n + 1 / 2} - \frac{Δ t}{2 ε_{0}} A A H^{n + 1 / 2} + \frac{Δ t}{2 ε_{0}} A B H^{n + 1} - A E^{n + 1 / 2} + \frac{Δ t}{2 ε_{0}} A J^{n + p_{2}}) .

Then, add Eq. (13) and Eq. (14) on their both sides, the update equation for the magnetic field in the leapfrog ADI-FDTD method can be obtained as

\begin{array}{l} (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) H^{n + 1} = (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) H^{n} + \frac{Δ t}{μ_{0}} (B E^{n + 1 / 2} - A E^{n + 1 / 2}) \\ + \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A (J^{n + p_{2}} - J^{n + p_{1}}) . \end{array}

In practice, it is not easy to compute the cross term $Δ t^{2} / 4 μ_{0} ε_{0} A (J^{n + p_{2}} - J^{n + p_{1}})$ in Eq. (15). To simplify it, we set p₁ = p₂ = 1/2 as

(I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) E^{n + 1 / 2} = (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) E^{n - 1 / 2} + \frac{Δ t}{ε_{0}} (A H^{n} - B H^{n}) - \frac{Δ t}{ε_{0}} J^{n}

(I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) H^{n + 1} = (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B) H^{n} + \frac{Δ t}{μ_{0}} (B E^{n + 1 / 2} - A E^{n + 1 / 2}) .

Here the Taylor series approximation

J^{n - 1 / 2} + J^{n + 1 / 2} = 2 J^{n} + O (Δ t^{2})

is used. Equations (16a) and (16b) are the equations of the leapfrog ADI-FDTD method applicable for modeling general dispersive media. The method is efficient since no mid-time step computation is needed with less memory consumed.

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

E^{n + 1 / 2} = E^{n - 1 / 2} + \frac{Δ t}{ε_{0}} (A H^{n} - B H^{n}) - \frac{Δ t}{ε_{0}} J^{n} + e^{n}

H^{n + 1} = H^{n} + \frac{Δ t}{μ_{0}} (B E^{n + 1 / 2} - A E^{n + 1 / 2}) + h^{n + 1 / 2}

where

e^{n} = \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B (E^{n + 1 / 2} - E^{n - 1 / 2}), h^{n + 1 / 2} = \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B (H^{n + 1} - H^{n}) .

Equations (17a) and (17b) show that the leapfrog ADI-FDTD can be regarded as the conventional FDTD method, with the second-order perturbation terms

e^{n}

and

h^{n + 1 / 2}

added.

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

σ = [\begin{matrix} σ_{x} & 0 & 0 \\ 0 & σ_{y} & 0 \\ 0 & 0 & σ_{z} \end{matrix}], ε_{r e} (ω) = [\begin{matrix} 1 + σ_{x} / j ω ε_{0} & 0 & 0 \\ 0 & 1 + σ_{y} / j ω ε_{0} & 0 \\ 0 & 0 & 1 + σ_{z} / j ω ε_{0} \end{matrix}] .

By substituting Eq. (18) into Eq. (2) and Eq. (4), we have

J^{n} = σ E^{n} = σ (E^{n + 1 / 2} + E^{n - 1 / 2}) / 2.

Furthermore, by substituting Eq. (19) into Eq. (16a), we obtain

(I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B + \frac{Δ t}{2 ε_{0}} σ) E^{n + 1 / 2} = (I - \frac{Δ t^{2}}{4 μ_{0} ε_{0}} A B - \frac{Δ t}{2 ε_{0}} σ) E^{n - 1 / 2} + \frac{Δ t}{ε_{0}} (A H^{n} - B H^{n}) .

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]

ε_{r e} (ω) = [\begin{matrix} ε_{d} (ω) & j ε_{g} (ω) & 0 \\ - j ε_{g} (ω) & ε_{d} (ω) & 0 \\ 0 & 0 & ε_{z} (ω) \end{matrix}]

ε_{d} (ω) = 1 + \frac{ω_{p}^{2} + ν_{c} ω_{p}^{2} / (j ω)}{(ν_{c}^{2} + ω_{b}^{2}) + 2 j ν_{c}^{} ω - ω^{2}}, ε_{g} (ω) = \frac{- ω_{b} ω_{p}^{2} / (j ω)}{(ν_{c}^{2} + ω_{b}^{2}) + 2 j ν_{c}^{} ω - ω^{2}}

ε_{z} (ω) = 1 - \frac{ω_{p}^{2}}{j ν_{c}^{} ω + ω^{2}}

where ω_p is the plasma frequency, ω_b is the cyclotron frequency, and ν_c is the electron collision rate. By substituting Eqs. (21a)-(21d) into Eq. (4), we obtain

J_{x} (ω) = \frac{ε_{0} ω_{p}^{2} j ω + ε_{0} ν_{c} ω_{p}^{2}}{(ν_{c}^{2} + ω_{b}^{2}) + 2 j ν_{c}^{} ω - ω^{2}} E_{x} + \frac{- ε_{0} ω_{b}^{} ω_{p}^{2}}{(ν_{c}^{2} + ω_{b}^{2}) + 2 j ν_{c}^{} ω - ω^{2}} E_{y} .

As discussed in [18], the first and second terms on the right-hand side of Eq. (22) can be replaced by two auxiliary variables

J_{x x}

and

J_{x y}

, i.e. $J_{x}^{n} = J_{x x}^{n} + J_{x y}^{n}$ . Both

J_{x x} (ω)

and

J_{x y} (ω)

have the general form as

J (ω) = \frac{a_{0} + a_{1} j ω}{b_{0} + b_{1} j ω + b_{2} {(j ω)}^{2}} E (ω)

where

a_{0}

,

a_{1}

,

b_{0}

,

b_{1}

and

b_{2}

are the real coefficients. Converting Eq. (23) into the time domain, we have

b_{0} J + b_{1} \frac{d J}{d t} + b_{2} \frac{d^{2} J}{d t^{2}} = a_{0} E + a_{1} \frac{d E}{d t} .

By discretizing and approximating the derivatives around the time step n-1, we obtain

b_{0} J^{n - 1} + b_{1} \frac{J^{n} - J^{n - 2}}{2 Δ t} + b_{2} \frac{J^{n} - 2 J^{n - 1} + J^{n - 2}}{Δ t^{2}} = a_{0} \frac{E^{n - 1 / 2} + E^{n - 3 / 2}}{2} + a_{1} \frac{E^{n - 1 / 2} - E^{n - 3 / 2}}{Δ t} .

and then, one iterative equation for J is derived as

J^{n} = \frac{4 b_{2} Δ t - 2 b_{0} Δ t^{2}}{b_{1} Δ t + 2 b_{2}} J^{n - 1} + \frac{b_{1} Δ t - 2 b_{2}}{b_{1} Δ t + 2 b_{2}} J^{n - 2} + \frac{a_{0} Δ t^{2} + 2 a_{1} Δ t}{b_{1} Δ t + 2 b_{2}} E^{n - 1 / 2} + \frac{a_{0} Δ t^{2} - 2 a_{1} Δ t}{b_{1} Δ t + 2 b_{2}} E^{n - 3 / 2} .

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 $σ$ can be calculated from its surface conductivity $σ_{s}$ by using $σ = σ_{s} / Δ z$ . In the frequency up to 10 THz, its $σ_{s}$ can be described by a 2-D Drude model [35], with $ε_{r e} (ω) = I + σ / (j ω ε_{0})$ , and

ε_{r e} (ω) = [\begin{matrix} ε_{d} (ω) & j ε_{g} (ω) & 0 \\ - j ε_{g} (ω) & ε_{d} (ω) & 0 \\ 0 & 0 & 0 \end{matrix}]

ε_{d} (ω) = 1 + \frac{σ_{0} + ν σ_{0} / (j ω)}{(ν_{}^{2} + ω_{c}^{2}) + 2 j ν ω - ω^{2}}

ε_{g} (ω) = \frac{- ω_{b} σ_{0} / (j ω)}{(ν_{}^{2} + ω_{c}^{2}) + 2 j ν ω - ω^{2}}

σ_{0} = \frac{2 e^{2} τ k_{B} T ν}{π ℏ^{2} ε_{0} Δ z} \ln (2 \cos h \frac{μ_{c}}{2 k_{B} T})

where

ν

is the phenomenological scattering rate,

ω_{c}

is the cyclotron frequency and

ω_{c} = e B_{0} v_{F}^{2} / μ_{c}

,

e

is the electron charge, v_F is the Fermi velocity, k_B is the Boltzmann constant,

μ_{c}

is the chemical potential, T is the operating temperature, and

ℏ

is the reduced Planck’s constant.

According to Eqs. (27a)-(27d), it is found that the $ε_{r e} (ω)$ of graphene sheet can be obtained by replacing the parameters of magnetized plasma in Eq. (21) with $ω_{p}^{2}$ $\to$ $σ_{0}$ , $ν_{c}^{}$ $\to$ $ν$ , $ω_{b}$ $\to$ $ω_{c}$ , and $ε_{z} (ω) = 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 $ω_{c} = 0$ in Eq. (27), i.e.

ε_{d} (ω) = 1 + [σ_{0} / (j ω)] / (ν + j ω)

and

ε_{g} (ω) = 0

. With the iterative equation Eq. (26) used, both J_x and J_y can be obtained easily.

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

E_{0} = [e / (ε_{0} π ℏ^{2} v_{F}^{2})] \int_{0}^{\infty} ε [f_{d} (ε) - f_{d} (ε + 2 μ_{c})] d ε

where

f_{d} (ε)

is the Femi-Dirac distribution function.

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.

Fig. 1 A rectangular waveguide of 40mm × 40mm loaded with a magnetized plasma slab of 5 mm in width.

Download Full Size | PDF

In our simulation, we set ∆x = ∆y = ∆z = ∆ = 1mm. A modulated Gaussian pulse of $\cos (2 π f_{0} t) \exp [- 4 π {(t - t_{0})}^{2} / τ^{2}]$ was used to excite electromagnetic field, where f₀ = 10 GHz, $Δ t = Δ / \sqrt{3} c_{0} \cdot CFLN$ , $τ = 160 Δ t / CFLN$ , and $t_{0} = τ$ . CFLN is the CFL number which measures the ratio of the time step to the CFL limit. The TE₁₀-mode was excited and its cutoff frequency is 3.75GHz. The plasma frequency ω_p = 2π × 10 Grad/s, ν = 10 GHz, and ω_b = 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 E_z -component at the observation point verses time with different CFLNs. It is seen that for CFLN = 1, the recorded E_z –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.

Fig. 2 The Recorded E_z –component versus time with CFLN = 1, 2, and 3; the result obtained with the conventional FDTD method is also plotted for comparison.

Download Full Size | PDF

3.2 Magnetized graphene sheet

When a graphene sheet is biased by a magnetostatic field B₀, 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₀/λ_min = 2 THz, the FDTD cells ∆x = ∆y = ∆z = ∆ = λ_min/20, ν = 20THz, T = 300K, μ_c = 0.1 eV, and $Δ t = Δ / (\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 x-direction and described by

E_{x}^{i n c} (t) = \exp [- 4 π {(t - t_{0})}^{2} / τ^{2}]

where

τ = 80 Δ t

and

t_{0} = 2 τ

. The observation point was set at the center of 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]:

T_{a n a l y t i c} (f) = \frac{2 \sqrt{{| 2 + j ω [ε_{d} ω - 1] η_{0} |}^{2} + {| ε_{0} ε_{g} ω_{}^{2} η_{0} |}^{2}}}{| {(2 + j ω [ε_{d} ω - 1] η_{0})}^{2} + {(ε_{0} ε_{g} ω_{}^{2} η_{0})}^{2} |}

where

η_{0}

is the wave impedance in the free space. The transmission coefficient in the FDTD simulation was computed with

T_{l e a p f r o g} (f) = \frac{\sqrt{{| F F T (E_{x, o b} (t)) |}^{2} + {| F F T (E_{y, o b} (t)) |}^{2}}}{| F F T (E_{i n} (t)) |}

where FFT represents the fast Fourier transformation.

Fig. 3 (a) The computational domain where a plane wave normally incident on an infinite graphene sheet. (b) The computed transmission coefficient |T| of the plane wave normally incident on an infinite graphene sheet biased by a magnetostatic field B₀ = 1 T. The analytical result is also plotted for comparison.

Download Full Size | PDF

Figure 3(b) shows the computed transmission coefficient using the leapfrog ADI-FDTD with B₀ = 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.5t₀ 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.

Fig. 4 E_y-component snapshot at t = 2.5t₀ for (a) CFLN = 1 and (b) CFLN = 3, respectively. The field is recorded at 10 cells away from the CPML in the y-z plane.

Download Full Size | PDF

3.3 SPPs propagating along the graphene sheet

It is known that when $Re [ε_{d} (ω)] < 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 $Δ t = Δ / (\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) = \sin (2 π f_{0} t)$ with f₀ = 6 THz. Therefore, the grid resolution is λ₀/∆ = c₀/(f₀∆) = 1000, and for such a fine structure the large CFLN is needed so as to reduce the simulation time.

Fig. 5 Computational domain for the simulation of SPPs on the graphene sheet.

Download Full Size | PDF

Figures 6(a)-6(d) show snapshots of the E_z – and E_y –components at time 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.

Fig. 6 Snapshots of the normalized electric field at time t = 2.887 ps. (a) E_z: leapfrog ADI-FDTD with CFLN = 10; (b) E_y: leapfrog ADI-FDTD with CFLN = 10; (c) E_z: conventional ADE-FDTD with CFLN = 1 [17]; and (d) E_y: conventional ADE-FDTD with CFLN = 1 [17].

Download Full Size | PDF

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 Nth-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]

One-step leapfrog ADI-FDTD method for simulating electromagnetic wave propagation in general dispersive media

Abstract

1. Introduction

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

2.1 Lossy media

2.2 Magnetized plasma

2.3 Graphene biased by a magnetostatic field

2.4 Graphene biased by an electrostatic field

3. Numerical results and discussion

3.1 The PEC Waveguide Loaded with Magnetized Plasma

3.2 Magnetized graphene sheet

3.3 SPPs propagating along the graphene sheet

4. Conclusion

Acknowledgments

References and links

Cited By

Figures (6)

Equations (47)

Optics Express