### I. Introduction

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.

### II. Numerical Permittivity of Dispersive FDTD Formulations

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

where

*ν**is the collision frequency,*_{c}*ω**is the plasma frequency, ω*_{p}*is the cyclotron frequency, and*_{b}*ɛ*_{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 aswhere

*J**,*_{x}*J**, and*_{y}*J**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]:*_{z}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 havewhere Δ

*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):

where

where

*u*(*t*) is the unit step function. Substitution Eq. (17) into Eqs. (15) and (16) yields the following equations: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 getAt

*t*= (*n*–1/2)Δ*t*, we havewhere

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

Using a similar procedure, we have

##### (27)

##### (28)

with

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:

where

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]

where

and

Using a similar procedure

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

##### (36)

##### (37)

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

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]

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

##### (40)

##### (41)

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

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

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

We use the plane-wave expansion again and have

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

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

*ω**=2π×50×10*_{p}^{9}rad/s, cyclotron frequency*ω**= 3×10*_{b}^{11}rad/s [20], and*v**= 5×10*_{c}^{11}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}*ω̃**) 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 (*_{p}*ω̃*) 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 asHere,

*ε̃**,*_{r}*ɛ**,*_{r}*f**, and*_{a}*f**indicate the numerical relative permittivity, analytical relative permittivity, minimum frequency, and maximum frequency in the frequency range of interest.*_{b}### 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.