Numerical Accuracy of Finite-Difference Time-Domain Formulations for Magnetized Plasma

Article information

J. Electromagn. Eng. Sci. 2022;22(3):195-201
Publication date (electronic) : 2022 April 21
doi : https://doi.org/10.26866/jees.2022.3.r.77
1Department of Electronic Engineering, Hanyang University, Seoul, Korea
*Corresponding Author: Kyung-Young Jung (e-mail: kyjung3@hanyang.ac.kr)
Received 2021 June 30; Revised 2021 July 21; Accepted 2021 July 22.

Abstract

The finite-difference time-domain (FDTD) has been widely used to analyze electromagnetic (EM) wave propagation in complex dispersive media. Over the past three decades, a variety of FDTD approaches for the EM wave propagation in magnetized plasma has been presented. In this work, we perform a comprehensive study on the numerical accuracy of four FDTD formulations for magnetized plasma including the JE convolution (JEC) method, the exponential time differencing (ETD) method, the H–J collocated auxiliary differential equation (ADE) method, and the E-J collocated ADE method. Toward this purpose, the numerical permittivity tensor of magnetized plasma in the four FDTD formulations are derived and then we analyze them to determine which approach can provide the best accuracy. It is found that the E-J collocated ADE method can lead to the best accuracy. Numerical examples awere performed to validate our investigations.

I. Introduction

The finite-difference time-domain (FDTD) method [15] has been popularly used to study various electromagnetic (EM) wave problems due to its accuracy, robustness, and simplicity. Over the past three decades, the FDTD method has been extended to simulate anisotropic dispersive media, including magnetized plasma. There are various FDTD formulations for EM analysis of magnetized plasma, including the JE convolution (JEC) method [69], exponential time differencing (ETD) method [10, 11], and auxiliary differential equation (ADE) method [1214]. In the JEC method, recursive convolution is involved in the relation between the current density and electric field. The ETD method avoids the time-consuming recursive convolution based on an efficient first-order approximation. Simple arithmetic implementation is involved in the ADE method, which can also be straightforwardly extended to nonlinear dispersive media, unlike other methods [15, 16]. There are two particular implementations in the ADE method for EM analysis of magnetized plasma. First, in the H-J collocated ADE method, magnetic field (H), and current density (J) are collocated in the same time domain when discretizing J [12]. Second, in the E-J collocated ADE method, electric field (E), and J components are collocated simultaneously [13]. Unlike the H-J collocated ADE method, the stability condition of the E-J collocated ADE method is independent of the medium properties and remains the same as the Courant stability limit for free space [14].

We perform a comprehensive study on the numerical accuracy of four dispersive FDTD formulations for modeling magnetized plasma. For this purpose, we derive the numerical permittivity tensor of magnetized plasma in the different methods and compare them with the corresponding analytical counterpart. Numerical examples are employed to investigate the numerical permittivity tensor of the JEC, ETD, H-J collocated ADE, and E-J collocated ADE methods in detail.

II. Numerical Permittivity of Dispersive FDTD Formulations

In magnetized plasma, the governing equations are given by [14]

(1) ×H=ε0Et+J
(2) ×E=-μ0Ht
(3) Jt+vcJ=ε0ωp2E+ωb×J,

where νc is the collision frequency, ωp is the plasma frequency, ωb is the cyclotron frequency, and ɛ0 and μ0 are the permittivity and permeability of free space, respectively. Note that cyclotron frequency is a function of the static magnetic field. The cross-product terms in Eq. (3) can lead to anisotropy of plasma. Thus, EM wave behavior depends on the direction of the static magnetic field relative to the EM wave propagation direction. It is assumed that the external static magnetic field in Cartesian coordinates is parallel to the z-axis; then, the component equations of Eq. (3) can be written as

(4) Jxt×νcJx=ε0ωp2Ex-ωbJy
(5) Jyt×νcJy=ε0ωp2Ey+ωbJx
(6) Jzt+νcJz=ε0ωp2Ez,

where Jx, Jy, and Jz are the current densities. Eqs. (4) and (5) indicate that two components of current density are coupled. Therefore, the update equations for polarization current density must be solved simultaneously. Moreover, for magnetized plasma, such as in the earth’s ionosphere, the electrons rotate about a steady magnetic field vector. Therefore, the plasma becomes nonreciprocal, and the scalar relationship between the electric flux density and electric field must be replaced by the tensor relation. The analytical permittivity tensor can be obtained using Eqs. (1), (4), and (5) as follows [17]:

(7) εxx(ω)=εyy(ω)=ε0[1-(ωp/ω)2{1-jνc/ω}{1-jνc/ω}2-(ωb/ω)2]
(8) εxy(ω)=εyx(ω)=ε0[j(ωp/ω)2(ωb/ω){1-jνc/ω}2-(ωb/ω)2].

Note that the z component of tensor permittivity is reciprocal because the external static magnetic field does not affect the wave behavior in that direction [15].

Due to the discrete nature of the FDTD technique, the numerical permittivity tensor of magnetized plasma in dispersive FDTD approaches is different from its analytical permittivity tensor. According to the standard FDTD method, the E field is defined at integer time steps and the H field is defined at half integer time steps. By applying the central difference scheme (CDS) to Eqs. (1) and (2), we have

(9) En+1=En+Δtε0(×H)n+1/2-Δtε0Jn+1/2
(10) Hn+1/2=Hn-1/2-Δtμ0(×E)n,

where Δt indicates the FDTD time step size and the superscript indicates the FDTD time step. In what follows, Eqs. (9) and (10) are employed, unless specified otherwise. In addition, the tilde characters indicate their numerical counterparts.

1. JEC method

For time-harmonic dependence, the following frequency-domain relation can be derived from Eqs. (4) and (5):

(11) Jx(ω)=σ(ω)Ex(ω)-ρ(ω)Jy(ω)
(12) Jy(ω)=σ(ω)Ey(ω)+ρ(ω)Jx(ω)

where

(13) σ(ω)=ε0ωp2jω+νc
(14) ρ(ω)=ωbjω+νc.

In the time domain, the above equations can be expressed by convolution [69]:

(15) Jx(t)=0t[σ(t-τ)Ex(τ)-ρ(t-τ)Jy(τ)]dτ
(16) Jy(t)=0t[σ(t-τ)Ey(τ)+ρ(t-τ)Jx(τ)]dτ
(17a) σ(t)=ε0ωp2e-νctu(t)
(17b) ρ(t)=ωbe-νctu(t)

where u(t) is the unit step function. Substitution Eq. (17) into Eqs. (15) and (16) yields the following equations:

(18) Jx(t)=ε0ωp2e-νct0teνcτEx(τ)dz-ωbe-νct0teνcτJy(τ)dz
(19) Jy(t)=ε0ωp2e-νct0teνcτEy(τ)dz+ωbe-νct0teνcτJx(τ)dz.

Let us consider the numerical permittivity tensor of magnetized plasma in the JEC method. To derive the numerical permittivity tensor, using Yee’s notation and t = (n+1/2)Δt in Eq. (18), we get

(20) Jxn+1/2=ε0ωp2e-νc(n+1/2)Δt0(n+1/2)ΔteνcτEx(τ)dz-ωbe-νc(n+1/2)Δt0(n+1/2)ΔteνcτJy(τ)dz.

At t = (n–1/2)Δt, we have

(21) Jxn-1/2=ε0ωp2e-νc(n-1/2)Δt0(n-1/2)ΔteνcτEx(τ)dz-ωbe-νc(n-1/2)Δt0(n-1/2)ΔteνcτJy(τ)dz.

Substituting Eq. (21) into Eq. (20), we obtain

(22) Jxn+1/2=e-νcΔtJxn-1/2+e-νc(n+1/2)Δt(n+1/2)Δt(n+1/2)Δtf(τ)dτ

where

(23) f(τ)=eνcτ[ε0ωp2Ex(τ)-ωbJy(τ)].

By using the Taylor series expansion of Eq. (22), we have

(24) (n-1/2)Δt(n+1/2)Δtf(τ)dτ=Δteνcτ[ε0ωp2Ex(τ)-ωbJy(τ)].

According to Eqs. (22) and (24), the following second-order approximation can be written as

(25) Jxn+1/2=e-νcΔtJxn-1/2+Δte-νcΔt/2ε0ωp2Exn-ωbΔt2e-νcΔt/2Jyn-1/2-ωbΔt2e-νcΔt/2Jyn+1/2.

Using a similar procedure, we have

(26) Jyn+1/2=e-νcΔtJyn-1/2+Δte-νcΔt/2ε0ωp2Eyn+ωbΔt2e-νcΔt/2Jxn-1/2+ωbΔt2e-νcΔt/2Jxn+1/2.

By using plane-wave expansion [18, 19] and applying some mathematical manipulations, we have

(27) Jx0=ε0ωp2Δt(2sinh β)(2sinh β)2+{ωbΔtcos (ωΔt/2)}2Ex0-ε0ωp2Δt2ωbcos (ωΔt/2)(2sinh β)2+{ωbΔtcos (ωΔt/2)}2Ey0
(28) Jy0=ε0ωp2Δt(2sinh β)(2sinh β)2+{ωbΔtcos (ωΔt/2)}2Ey0+ε0ωp2Δt2ωbcos (ωΔt/2)(2sinh β)2+{ωbΔtcos (ωΔt/2)}2Ex0

with

β=νcΔt+jωΔt2.

Substituting Eqs. (27) and (28) into the plane-wave expansion version of Eq. (9) and then rearranging the resulting equation can lead to the following numerical permittivity tensor:

(29) ε˜r,xx=1-(ω˜p/ω˜)2(1-jν˜c/ω˜)(1-jν˜c/ω˜)2-(ω˜b/ω˜)2
(30) ε˜r,xy=j(ω˜p/ω˜)2(ω˜b/ω˜)(1-jν˜c/ω˜)2-(ω˜b/ω˜)2,

where

ω˜=2Δttan(ωΔt2)ν˜c=2Δttanh(νcΔt2)ω˜b=ωbcosh(νcΔt/2)ω˜p=ωpcos(ωΔt/2)1cosh(νcΔt/2).

The numerical permittivity converges to analytical permittivity as the time-step size approaches zero.

2. ETD method

In this subsection, we derive the numerical permittivity tensor of magnetized plasma using the ETD method. In the ETD method, the discrete form for Eq. (4) can be written as [10, 11]

(31) Jxn+1/2=e-νcΔtJxn-1/2+e-νcΔt/2-Δt/2Δt/2f(τ)dτ,

where

(32) f(τ)=eνcτ[ε0ωp2Ex(nΔt)-ωbJy(nΔt)]

and

(33) Jyn=12(Jyn+1/2+Jyn-1/2).

Substituting Eqs. (32) and (33) into (31) yields

(34) Jxn+1/2=e-νcΔtJxn-1/2+1νc(1-e-νcΔt)[ε0ωp2Exn-ωb2(Jyn+1/2+Jyn-1/2)].

Using a similar procedure

(35) Jyn+1/2=e-νcΔtJyn-1/2+1νc(1-e-νcΔt)[ε0ωp2Eyn+ωb2(Jxn+1/2+Jxn-1/2)].

By utilizing plane-wave expansion and applying some mathematical manipulations, we have

(36) Jx0=(νcsinh β)ε0ωp2sinh (νcΔt/2)(νcsinh β)2+[ωbsinh (νcΔt/2)cos (ωΔt/2)]2Ex0-ε0ωp2ωbcos(ωΔt/2)[sinh (νcΔt/2)]2(νcsinh β)2+[ωbsinh (νcΔt/2)cos (ωΔt/2)]2Ey0
(37) Jy0=(νcsinh β)ε0ωp2sinh (νcΔt/2)(νcsinh β)2+[ωbsinh (νcΔt/2)cos (ωΔt/2)]2Ey0+ε0ωp2ωbcos (ωΔt/2){sinh (νcΔt/2)}2(νcsinh β)2+[ωbsinh (νcΔt/2)cos (ωΔt/2)]2Ex0.

Substituting Eqs. (36) and (37) into the plane-wave expansion version of Eq. (9) and then rearranging the resulting equation, we have the same Eqs. (29) and (30) with

ω˜=2Δttan (ωΔt2)ν˜c=2Δttanh (νcΔt2)ω˜b=ωbtanh (νcΔt/2)νcΔt/2ω˜p=ωpcos (ωΔt/2)tanh (νcΔt/2)νcΔt/2.

Again, as the time step size approaches zero, the numerical permittivity converges to the analytical permittivity.

3. H-J collocated ADE method

In the H-J collocated ADE method, we can write [12]

(38) Jxn+1/2=2-vcΔt2+vcΔtJxn-1/2+2Δt2+vcΔt[ε0ωp2Exn-ωb2(Jyn+1/2+Jyn-1/2)]
(39) Jyn+1/2=2-vcΔt2+vcΔtJyn-1/2+2Δt2+vcΔt[ε0ωp2Eyn+ωb2(Jxn+1/2+Jxn-1/2)].

Now, we use the plane-wave expansion again and have

(40) Jx0=ε0ωp2cos (ωΔt/2)[j2Δttan(ωΔt2)+νc][j2Δttan (ωΔt2)+νc]2+ωb2Ex0-ε0ωp2ωbcos(ωΔt/2)[j2Δttan(ωΔt2)+νc]2+ωb2Ey0
(41) Jy0=ε0ωp2cos (ωΔt/2)[j2Δttan (ωΔt2)+νc][j2Δttan (ωΔt2)+νc]2+ωb2Ey0+ε0ωp2ωbcos(ωΔt/2)[j2Δttan (ωΔt2)+νc]2+ωb2Ex0.

Following the procedure described above, the numerical permittivity tensor is the same as Eqs. (29) and (30) with

ω˜=2Δttan (ωΔt2)ν˜c=νcω˜b=ωbω˜p=ωpcos (ωΔt/2).

Note that the collision and cyclotron frequencies are the same as the analytical ones.

4. E-J collocated ADE method

Unlike the other methods considered above, the current density vectors are collocated at the same time step and position of the electric field vectors in the E-J collocated ADE method. Therefore, the FDTD update equation for Ampere’s law can be written as

(42) En+1=En+Δtε0(×H)n+1/2-Δtε0(Jn+1+Jn2).

According to [13], the update equations of the current density in the E-J collocated ADE method can be written as

(43) Jxn+1=2-vcΔt2+vcΔtJxn+2Δt2+vcΔt[ε0ωp22(Exn+1+Exn)-ωb2(Jyn+1+Jyn)]
(44) Jyn+1=2-vcΔt2+vcΔtJyn+2Δt2+vcΔt[ε0ωp22(Eyn+1+Eyn)+ωb2(Jxn+1+Jxn)].

We use the plane-wave expansion again and have

(45) Jx0=ε0ωp2[j2Δttan(ωΔt2)+νc][j2Δttan(ωΔt2)+νc]2+ωb2Ex0-ε0ωp2ωb[j2Δttan(ωΔt2)+νc]2+ωb2Ey0
(46) Jy0=ε0ωp2[j2Δttan(ωΔt2)+νc][j2Δttan(ωΔt2)+νc]2+ωb2Ex0+ε0ωp2ωb[j2Δttan(ωΔt2)+νc]2+ωb2Ex0.

The numerical permittivity tensor of magnetized plasma in the E-J collocated ADE method can be obtained by a similar procedure, and we have the same as Eqs. (29) and (30) with

ω˜=2Δttan (ωΔt2)ν˜c=νcω˜b=ωbω˜p=ωp.

Note that only the angular frequency has a numerical value different from the analytical one, which implies that the E-J collocated ADE method can lead to better results than the other methods considered in this study.

III. Numerical Examples

In this section, we investigate the numerical accuracy of the JEC, ETD, H-J collocated ADE, and E-J collocated ADE methods. We assume that the plasma frequency ωp=2π×50×109 rad/s, cyclotron frequency ωb = 3×1011 rad/s [20], and vc = 5×1011 Hz. The simulation frequency ranges from 10 to 90 GHz, and the spatial cell is Δz = 75 μm. The FDTD time size is Δt = 1.111 ps, which is temporal points per period (PPP) equal to 10.

Fig. 1 shows the numerical relative permittivity tensor of magnetized plasma. As shown in the figure, large differences between the analytical and numerical permittivity tensors in the JEC and ETD methods are observed because all four numerical parameters (ω̃,c,ω̃b,ω̃p) are not the same as the analytical ones. The H-J collocated ADE method has two numerical values (ω̃,ω̃p) different from the analytical ones, but the E-J collocated ADE method has only one numerical value (ω̃) different from the analytical one. Therefore, the E-J collocated ADE method can yield the best numerical accuracy among the four FDTD methods considered in this study. Next, let us investigate the root-mean-square (RMS) error of the numerical relative permittivity in the four FDTD methods versus the FDTD time step size. The RMS error is defined as

Fig. 1

Real and imaginary parts of ε̃rxx and ε̃rxy.

(47) RMS error=fafb|ε˜r-εr|2dffafb|εr|2df.

Here, ε̃r,ɛr,fa, and fb indicate the numerical relative permittivity, analytical relative permittivity, minimum frequency, and maximum frequency in the frequency range of interest.

Figs. 2 and 3 show the RMS error of the numerical relative permittivity tensor. As the time step size decreases (PPP increases), the RMS error decreases. Again, the E-J collocated ADE method yields the best accuracy.

Fig. 2

RMS error of ε̃rx versus PPP.

Fig. 3

RMS error of ε̃rxy versus PPP.

IV. Conclusion

We investigated the numerical accuracy of various FDTD formulations for magnetized plasma. The exact expressions of the numerical permittivity tensor were derived, and the E-J collocated ADE method was found to yield the best accuracy. In addition, numerical examples were used to validate our investigation.

Acknowledgments

This research was supported by an Institute of Information & Communications Technology Planning & Evaluation (IITP) grant funded by the Korea government (MSIT) (No. 2019-0-00098, Advanced and integrated software development for electromagnetic analysis) and by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2020R1F1A1055444).

References

1. Taflove A, Hagness SC. Computational Electrodynamics: The Finite-Difference Time-Domain Method 3rd edth ed. Norwood, MA: Artech House; 2005.
2. Chilton R, Jung KY, Lee R, Teixeira FL. Frozen modes in parallel-plate waveguides loaded with magnetic photonic crystals. IEEE Transactions on Microwave Theory and Techniques 55(12):2631–2641. Dec. 2007;
3. Chung H, Jung KY, Tee XT, Bermel P. Time domain simulation of tandem silicon” solar cells with optimal textured light trapping enabled by the quadratic complex rational function. Optics Express 22(S3):A818–A832. May. 2014;
4. Jung K-Y, Teixeira FL. Numerical study of photonic crystals with a split band edge: Polarization dependence and sensitivity analysis. Physics Review A 78(4):043826. Oct. 2008;
5. Park J, Jung K-Y. Numerical stability of modified Lorentz FDTD unified from various dispersion models. Optics Express 29(14):21639–21654. July. 2021;
6. Chen Q, Katsurai M, Aoyagi PH. An FDTD formulation for dispersive media using a current density. IEEE Transactions on Antennas and Propagation 46(10):1739–1746. Oct. 1998;
7. Liu S, Mo J, Yuan N. FDTD analysis of electromagnetic reflection by conductive plane covered with magnetized inhomogeneous plasmas. International Journal of Infrared and Millimeter Waves 23(12):1803–1815. Dec. 2002;
8. Xu L, Yuan N. FDTD formulations for scattering from 3-D anisotropic magnetized plasma objects. IEEE Antennas and Wireless Propagation Letters 5:335–338. 2006;
9. Chen W, Guo L, Li J, Liu S. Research on the FDTD method of electromagnetic wave scattering characteristics in time-varying and spatially nonuniform plasma sheath. IEEE Transactions on Plasma Science 44(12):3235–3242. Dec. 2016;
10. Huang SJ, Li F. FDTD implementation for magnetoplasma medium using exponential time differencing. IEEE Microwave and Wireless Components Letters 15(3):183–185. Mar. 2005;
11. Huang SJ. Exponential time differencing FDTD formulation for plasma. Microwave and Optical Technology Letters 49(6):1363–1364. Jun. 2007;
12. Young JL. A full finite difference time domain implementation for radio wave propagation in a plasma. Radio Science 29(6):1513–1522. Nov.-Dec. 1994;
13. Yu Y, Simpson J. An E-J collocated 3-D FDTD model of electromagnetic wave propagation in magnetized cold plasma. IEEE Transactions on Antennas and Propagation 58(2):469–478. Feb. 2010;
14. Cho J, Park M, Jung K-Y. Perfectly matched layer for accurate FDTD for anisotropic magnetized plasma. Journal of Electromagnetic Engineering and Science 20(4):277–284. Oct. 2020;
15. Choi H, Kim Y-H, Baek J-W, Jung K-Y. Accurate and efficient finite-difference time-domain simulation compared with CCPR model for complex dispersive media. IEEE Access 7:160498–160505. Nov. 2019;
16. Choi H, Baek J-W, Jung K-Y. Comprehensive study on numerical aspects of modified Lorentz model based dispersive FDTD formulations. IEEE Transactions on Antennas and Propagation 67(12):7643–7648. Dec. 2019;
17. Hunsberger F, Luebbers R, Kunz K. Finite-difference time-domain analysis of gyrotropic media-I: Magnetized plasma. IEEE Transactions on Antennas and Propagation 40(2):1489–1495. Dec. 1992;
18. Choi H, Baek J-W, Jung K-Y. Numerical stability and accuracy of CCPR-FDTD for dispersive medi. IEEE Transactions on Antennas and Propagation 68(11):7717–7720. Nov. 2020;
19. Kim Y-J, Jung K-Y. Accurate and efficient finite-difference time domain formulation of dusty plasma. IEEE Transactions on Antennas and Propagation 10.1109/TAP.2021.3069542. Apr. 2021;
20. Liu S, Mo J, Yuan N. Piecewise linear current density recursive convolution FDTD implementation for anisotropic magnetized plasmas. IEEE Microwave and Wireless Components Letter 14(5):222–224. May. 2004;

Biography

Jeahoon Cho received his B.S. degree in Communication Engineering from Daejin University, Pocheon, Rep. of Korea, and his M.S. and Ph.D. degrees in electronics and computer engineering from Hanyang University, Seoul, Rep. of Korea, in 2004, 2006, and 2015, respectively. From 2015 to August 2016, he was a postdoctoral researcher at Hanyang University. Since September 2016, he has worked at Hanyang University, where he is currently a research professor. His current research interests include computational electromagnetics and EMP/EMI/EMC analysis.

Min-Seok Park received his B.S. degree in Department of Electrical Engineering from Myongi University, Yongin, Rep. of Korea, in 2015, and his M.S. degree in electrical engineering from Hanyang University, Seoul, Rep. of Korea, in 2017. From 2018 to 2020, he participated in EM-Tech’s circuit and product development. He is currently pursuing a Ph.D. degree in electrical and computer engineering. His current research interests include computational electromagnetics, wave propagation, and multi-physics.

Kyung-Young Jung received his B.S. and M.S. degrees in electrical engineering from Hanyang University, Seoul, Rep. of Korea, in 1996 and 1998, respectively, and his Ph.D. degree in electrical and computer engineering from The Ohio State University, Columbus, USA, in 2008. From 2008 to 2009, he was a postdoctoral researcher at the Ohio State University, and from 2009 to 2010, he was an assistant professor at the Department of Electrical and Computer Engineering, Ajou University, Suwon, Rep. of Korea. Since 2011, he has worked at Hanyang University, where he is now a professor in the Department of Electronic Engineering. His current research interests include computational electromagnetics, bioelectromagnetics, and nanoelectromagnetics. Dr. Jung was a recipient of the Graduate Study Abroad Scholarship from the National Research Foundation of Korea, the Presidential Fellowship from The Ohio State University, the HYU Distinguished Teaching Professor Award from Hanyang University, and the Outstanding Research Award from the Korean Institute of Electromagnetic Engineering Society.

Article information Continued

Fig. 1

Real and imaginary parts of ε̃rxx and ε̃rxy.

Fig. 2

RMS error of ε̃rx versus PPP.

Fig. 3

RMS error of ε̃rxy versus PPP.