The finite-difference time-domain (FDTD) method [1–5] 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 [6–9], exponential time differencing (ETD) method [10, 11], and auxiliary differential equation (ADE) method [12–14]. 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.

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

where ν_c is the collision frequency, ω_p is the plasma frequency, ω_b is the cyclotron frequency, and ɛ₀ and μ₀ 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

where J_x, J_y, and J_z 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]:

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

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.

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×10⁹ rad/s, cyclotron frequency ω_b = 3×10¹¹ rad/s [20], and v_c = 5×10¹¹ 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

Here, ε̃_r,ɛ_r,f_a, and f_b 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.

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.