An Efficient Numerical Formulation for Wave Propagation in Magnetized Plasma Using PITD Method

A modified precise-integration time-domain (PITD) formulation is presented to model the wave propagation in magnetized plasma based on the auxiliary differential equation (ADE). The most prominent advantage of this algorithm is using a time-step size which is larger than the maximum value of the Courant–Friedrich–Levy (CFL) condition to achieve the simulation with a satisfying accuracy. In this formulation, Maxwell’s equations in magnetized plasma are obtained by using the auxiliary variables and equations. Then, the spatial derivative is approximated by the second-order finite-difference method only, and the precise integration (PI) scheme is used to solve the resulting ordinary differential equations (ODEs). The numerical stability and dispersion error of this modified method are discussed in detail in magnetized plasma. The stability analysis validates that the simulated time-step size of this method can be chosen much larger than that of the CFL condition in the finite-difference time-domain (FDTD) simulations. According to the numerical dispersion analysis, the range of the relative error in this method is 10−6 to 5×10−4 when the electromagnetic wave frequency is from 1 GHz to 100 GHz. More particularly, it should be emphasized that the numerical dispersion error is almost invariant under different time-step sizes which is similar to the conventional PITD method in the free space. This means that with the increase of the time-step size, the presented method still has a lower computational error in the simulations. Numerical experiments verify that the presented method is reliable and efficient for the magnetized plasma problems. Compared with the formulations based on the FDTD method, e.g., the ADE-FDTD method and the JE convolution FDTD (JEC-FDTD) method, the modified algorithm in this paper can employ a larger time step and has simpler iterative formulas so as to reduce the execution time. Moreover, it is found that the presented method is more accurate than the methods based on the FDTD scheme, especially in the high frequency range, according to the results of the magnetized plasma slab. In conclusion, the presented method is efficient and accurate for simulating the wave propagation in magnetized plasma.


Introduction
The simulations of the electromagnetic (EM) wave propagation in the magnetized plasma are attractive and have a wide range of applications, e.g., high frequency components, PCB design, microstrip antenna, and so on [1][2][3][4][5][6]. Recently, the finite-difference time-domain (FDTD) formulation is and the dispersion error are analyzed numerically. The stability analysis verifies that the numerical stability criterion of the presented method in magnetized plasma is much looser than the CFL stability condition of the FDTD methods so as to increase the maximum allowable time step, and the numerical dispersion errors are almost invariant when the time-step size is increased. Thus, this method has the potential to balance both the efficiency and the accuracy. The magnetized plasma slab and the magnetized plasma filled cavity are simulated to validate that the modified PITD method in the paper is reliable and efficient. The analyses of the numerical results indicate that the presented method can provide an evident reduction of the execution time by using a larger simulated time step, meanwhile, the computational error of the presented algorithm is also lower than those of the formulations based on the FDTD scheme, especially in the high frequency range.  Due to the advantages mentioned above, the PITD method is a promising approach to model the EM wave propagation in magnetized plasma efficiently. The resulting Maxwell's equations of magnetized plasma are firstly obtained by employing the auxiliary variables and equations related to the current density. Then, the second-order accurate finite-difference formulation is used to approximate the spatial derivative in the presented method, and several ordinary differential equations (ODEs) with respect to the time derivative are obtained directly. Finally, we use the PI scheme to solve the ODEs. After establishing the modified PITD method in magnetized plasma, the stability condition and the dispersion error are analyzed numerically. The stability analysis verifies that the numerical stability criterion of the presented method in magnetized plasma is much looser  Due to the advantages mentioned above, the PITD method is a promising approach to model the EM wave propagation in magnetized plasma efficiently. The resulting Maxwell's equations of magnetized plasma are firstly obtained by employing the auxiliary variables and equations related to the current density. Then, the second-order accurate finite-difference formulation is used to approximate the spatial derivative in the presented method, and several ordinary differential equations (ODEs) with respect to the time derivative are obtained directly. Finally, we use the PI scheme to solve the ODEs. After establishing the modified PITD method in magnetized plasma, the stability condition and the dispersion error are analyzed numerically. The stability analysis verifies that the numerical stability criterion of the presented method in magnetized plasma is much looser than the CFL stability condition of the FDTD methods so as to increase the maximum allowable time The rest of this paper is organized as follows. The formulation of the presented method is introduced in Section 2. The stability analysis and the dispersion analysis are discussed numerically in Section 3. Numerical results are simulated to verify the efficiency and the accuracy of the presented algorithm in Section 4. Finally, conclusions are drawn in Section 5.

Resulting Maxwell's Equations of Magnetized Plasma
The curl Maxwell's equations for describing the magnetized plasma problem is firstly established by employing the auxiliary variables and equations related to the current density J(t). The resulting matrix form is shown as follows: where ω p and γ are the natural angular frequency and the collision frequency, respectively; ω b = eB 0 /m is the electron cyclotron angular frequency, wherein B 0 is the applied magnetic field, e is the electric quantity, and m is the mass of the electron. Assume that the applied magnetic field is z-direction in the following analysis, and Equation (1) can be expanded as follows: Here, we know that in the formulations based on the FDTD scheme, both the spatial and time derivatives are discretized by using the finite difference technique to obtain a set of algebraic equations from the Maxwell's equations. However, in the proposed PITD algorithm, only the spatial derivative is discretized, and several ODEs are obtained temporarily and written as matrix form: where Y = E x , E y , E z , H x , H y , H z , J x , J y , J z T , the coefficient matrix H is determined by both the EM parameters of the medium and the spatial step of the simulation, and g(t) is an inhomogeneous term related to the excitations.
The analytical and discrete solutions of Equation (3) are: and where Y n = Y(n∆t) is the discrete form of Y(t) and T = exp(H∆t) which can be calculated by the the PI technique.

PI Technique Review
The exponential matrix T is firstly reformulated as: where s = ∆t/l, l = 2 N , and N is a preselected arbitrary integer. If N is large enough, the interval of s will be extremely small. Then, the Taylor series expansion is employed to approximate exp(Hs) with a high precision as shown in the following: where: T a = (Hs) + (Hs) 2 2! + (Hs) 3 3! + (Hs) 4 4! , (8) and evidently, It should be noted that if T a is added to the identity matrix I directly, T a will be neglected because T a is extremely small, which leads to a precision reduction of the exponential matrix. Therefore, it is evident that T a should be operated in the process.
The matrix T is computed as follows: and Start with Equation (8) to compute T a and then run the following instruction, the exponential matrix T can be calculated: end do Equations (8), (12), and (13) constitute the whole process of the PI technique to calculate the exponential matrix. Relying on the previous work, the simulated time step of the PITD method is much larger than the maximum allowable value of the CFL stability condition of the FDTD formulation in the lossless or lossy problems, which is very significant in the full wave analysis. For the magnetized plasma material, we believe the application of the PI technique can achieve the same effect and the following stability analysis and numerical experimentations will prove this point.
Furthermore, for the integration on the right-hand side of Equation (5), a linear variation of the term g(t) is assumed within the interval (t n , t n+1 ), expressed as: Substitute the above expression in the integration, we have the following recursive form solution: In most cases, Equation (15) is unavailable directly since the matrix M is noninvertible. To mitigate this problem, the three-points Gauss integral formula is adopted to calculate the integration in Equation (5), and the recurrence formula is obtained as follows:

Stability Analysis
In this subsection, the amplitude of the eigenvalues of the exponential matrix T is used to discuss the numerical stability of the presented PITD algorithm in the magnetized plasma. Von Neumann stability criterion indicates that if the amplitudes of all the eigenvalues of the exponential matrix T are no larger than unity, the update equations of the presented PITD algorithm will be stable.
We consider the 2-D case in the following analysis. The preselected integer N is selected as 20, and l = 2 20 in the proposed method. The natural angular frequency of the magnetized plasma discussed is 2π × 28.7 × 10 9 rad/s, the collision frequency is 20 GHz, and the electron cyclotron frequency is set as 1.0 × 10 11 rad/s. Figure 3 graphs the comparison of the unit circle and the eigenvalues of the exponential matrix T in the complex plane when the time-step size is 10 6 ∆t CFL . Here, ∆t CFL is the upper limit time-step size of the conventional FDTD method. It is clear seen that all the eigenvalues are within or on the unit circle, which means that the presented PITD method of the magnetized plasma is stable under such a large time-step size. Therefore, the proposed formulation can use a time-step size much larger than the maximum value of the CFL stability condition to achieve the simulation of the magnetized plasma problems.
in Equation (5), and the recurrence formula is obtained as follows:

Stability Analysis
In this subsection, the amplitude of the eigenvalues of the exponential matrix T is used to discuss the numerical stability of the presented PITD algorithm in the magnetized plasma. Von Neumann stability criterion indicates that if the amplitudes of all the eigenvalues of the exponential matrix T are no larger than unity, the update equations of the presented PITD algorithm will be stable.
We consider the 2-D case in the following analysis. The preselected integer N is selected as 20, and 20 2 l = in the proposed method. The natural angular frequency of the magnetized plasma discussed is 28.7 10 rad/s π , the collision frequency is 20 GHz, and the electron cyclotron frequency is set as × 11 1.0 10 rad/s . Figure 3 graphs the comparison of the unit circle and the eigenvalues of the exponential matrix T in the complex plane when the time-step size is Here, CFL t Δ is the upper limit time-step size of the conventional FDTD method. It is clear seen that all the eigenvalues are within or on the unit circle, which means that the presented PITD method of the magnetized plasma is stable under such a large time-step size. Therefore, the proposed formulation can use a time-step size much larger than the maximum value of the CFL stability condition to achieve the simulation of the magnetized plasma problems.

Numerical Dispersion Analysis
In this subsection, the dispersion performance of the presented PITD method in magnetized plasma is discussed numerically by adopting the Fourier approach. The dispersion performance of the presented formulation is described by the differences between the numerical wave number k num and the analytical wave number k anal . Suppose c is the velocity of light in the vacuum, the analytic wave number of the left-hand circularly polarized (LCP) EM wave is: and the analytic wave number of the right-hand circularly polarized (RCP) EM waves is: Electronics 2020, 9, 1575 7 of 18 Assuming that the monochromatic plane wave propagates in the magnetized plasma, the field components are expressed as: where k is the wave number, X 0 and ω are the amplitude and the angular frequency of the EM wave, respectively. The discrete form of Equation (19) is obtained as follows: where m and n are the space index and the time index, respectively. Here, we consider the 1-D case, and the vector X includes the field components E x , E y , H x , H y , J x , and J y . Substituting the discrete form of the field components into Equation (2) for the 1-D case, a homogenous system of ODEs can be obtained and written in matrix form as: Here, the coefficient matrix H 1 is: where: The following discrete iterative formula is used to solve the ODEs Equation (21): where the exponential matrix T 1 is: Then, we have: Since it is true for any X 0 that is nonzero, the determinant of the coefficient matrix e jωt I − T should be zero: In the following analysis, the numerical dispersion error and the numerical dissipation error are defined to describe the precision of the presented PITD method in the magnetized plasma. The definition of the relative dispersion error is (Re(k num )−Re(k anal ))/Re(k anal ), and it is related to the phase error. The definition of the relative dissipation error is (Im(k num )−Im(k anal ))/Im(k anal ), and it is related to the amplitude error [27].
It is assumed that ∆z = 75 µm, ∆t = 0.125 ps, and ω b = 3.0 × 10 11 rad/s. The preselected integer N is chosen as 20, and l = 2 20 in the presented PITD method. The solutions of numerical wave number are computed by Equation (27).

Effect of Wave Frequency on Numerical Error
The natural angular frequency of the magnetized plasma is ω p = 2π × 50 × 10 9 rad/s. The collision frequency is 20 GHz. Figure 4a,b graphs the relative dispersion and relative dissipation errors with respect to the wave frequency ω for both the LCP and RCP waves, respectively. We can see that both the dispersion and dissipation errors increase monotonically with the wave frequency. Furthermore, the relative dispersion error is higher than the relative dissipation error in the LCP wave, and is lower than the relative dissipation error in the RCP wave. In addition, the relative error range of the proposed PITD method is 10 −6 to 5 × 10 −4 when the EM wave frequency is from 1 GHz to 100 GHz.

Effect of the Natural Angular Frequency on Numerical Error
The EM wave frequency is set as 50 GHz. The collision frequency of the magnetized plasma is 20 GHz. Figure 5a,b graphs the relative dispersion and relative dissipation errors with respect to the natural angular frequency ωp for both the LCP and RCP waves, respectively. For the LCP wave, the relative dispersion error curve has lower-peak error when ωp/2π is 18 GHz, and the relative dissipation error curve has lower-peak errors when ωp/2π are 4 GHz and 21 GHz. For the RCP wave, the relative dispersion error curve has no lower-error peak, and the relative dissipation error curve has lower-peak error when ωp/2π is 6 GHz. Moreover, both the relative dispersion and relative dissipation errors increase monotonically with the natural frequency when ωp/2π is larger than the

Effect of the Natural Angular Frequency on Numerical Error
The EM wave frequency is set as 50 GHz. The collision frequency of the magnetized plasma is 20 GHz. Figure 5a,b graphs the relative dispersion and relative dissipation errors with respect to the natural angular frequency ω p for both the LCP and RCP waves, respectively. For the LCP wave, the relative dispersion error curve has lower-peak error when ω p /2π is 18 GHz, and the relative dissipation error curve has lower-peak errors when ω p /2π are 4 GHz and 21 GHz. For the RCP wave, the relative dispersion error curve has no lower-error peak, and the relative dissipation error curve has lower-peak error when ω p /2π is 6 GHz. Moreover, both the relative dispersion and relative dissipation errors increase monotonically with the natural frequency when ω p /2π is larger than the frequency of the lower-peak error.

Effect of the Plasma Collision Frequency on Numerical Error
The EM wave frequency is set as 50 GHz, and the natural angular frequency is = × × 9 2 50 10 rad/s p ω π . Figure 6a,b graphs the relative dispersion and relative dissipation errors with respect to the collision frequency γ for both the LCP and RCP waves, respectively. It is found that the relative dispersion and relative dissipation errors are slightly decreased when the collision frequency is increased.

Effect of the Plasma Collision Frequency on Numerical Error
The EM wave frequency is set as 50 GHz, and the natural angular frequency is ω p = 2π × 50 × 10 9 rad/s. Figure 6a,b graphs the relative dispersion and relative dissipation errors with respect to the collision frequency γ for both the LCP and RCP waves, respectively. It is found that the relative dispersion and relative dissipation errors are slightly decreased when the collision frequency is increased.

Effect of the Time-Step Size on Numerical Error
It is assumed that the natural angular frequency is × × 9 2 50 10 rad/s π , and the collision frequency is 20 GHz. Figure 7 graphs the relative dispersion and relative dissipation errors with respect to the wave frequency ω for both the LCP and RCP waves under different Courant number S, respectively. It is clear that all the curves are in agreement. Figure 8 shows the relative dispersion and relative dissipation errors with respect to the Courant number when the EM wave frequency is 50 GHz. Figure 8 indicates that the relative dispersion and relative dissipation errors are almost invariant when the Courant numbers is increased. These mean that the proposed method can maintain a lower numerical dispersion error in the simulations when the time step is increased.

Effect of the Time-Step Size on Numerical Error
It is assumed that the natural angular frequency is 2π × 50 × 10 9 rad/s, and the collision frequency is 20 GHz. Figure 7 graphs the relative dispersion and relative dissipation errors with respect to the wave frequency ω for both the LCP and RCP waves under different Courant number S, respectively. It is clear that all the curves are in agreement. Figure 8 shows the relative dispersion and relative dissipation errors with respect to the Courant number when the EM wave frequency is 50 GHz. Figure 8 indicates that the relative dispersion and relative dissipation errors are almost invariant when the Courant numbers is increased. These mean that the proposed method can maintain a lower numerical dispersion error in the simulations when the time step is increased. Electronics 2020, 9,

Numerical Experiments
The performance of the presented PITD method are verified by two typical magnetized plasma examples which are also solved by the analytical formulas and the formulations based on the FDTD method, respectively, for comparison.

Magnetized Plasma Slab
As the first example, a magnetized plasma slab is simulated to validate the efficiency and the accuracy of the modified PITD algorithm in this paper. The diagram of the infinite magnetized plasma slab in infinite free space is shown in Figure 9. The numerical reflection coefficient and transmission coefficient are computed by the presented method, the JEC-FDTD method and the ADE-FDTD method, respectively. The results are also compared with the analytical solution.

Numerical Experiments
The performance of the presented PITD method are verified by two typical magnetized plasma examples which are also solved by the analytical formulas and the formulations based on the FDTD method, respectively, for comparison.

Magnetized Plasma Slab
As the first example, a magnetized plasma slab is simulated to validate the efficiency and the accuracy of the modified PITD algorithm in this paper. The diagram of the infinite magnetized plasma slab in infinite free space is shown in Figure 9. The numerical reflection coefficient and transmission coefficient are computed by the presented method, the JEC-FDTD method and the ADE-FDTD method, respectively. The results are also compared with the analytical solution. (30) Figures 10-13 show the magnitude and the phase of the complex reflection coefficient and transmission coefficient of the magnetized plasma slab calculated by the JEC-FDTD method, the ADE-FDTD method, the proposed PITD method, and the analytical solution, respectively. We can clearly see that the computational results of the proposed PITD method are coincident with the analytical solutions on both the magnitude and the phase. The magnetized plasma slab is 1.5 cm thick, and divided into 200 cells, i.e., the space step is 75 µm. The main computing region is composed of 600 free space cells (the space indexes are from 1 to 300, and 501 to 800) and 200 magnetized plasma cells (the space indexes are from 301 to 500). Perfectly matched layer (PML) is employed as the absorption boundary condition to eliminate the reflection error. The time-step size of the methods based on the FDTD formulation is set to ∆t FDTD = 0.125 ps, and the time-step size of the proposed PITD method is 5 times of ∆t FDTD (i.e., ∆t ProposedPITD = 0.625 ps). The parameters of the magnetized plasma are shown as follows:  According to Figure 11, Figure 12a and Figure 13, it should be noted that the solutions of the proposed PITD method is more accurate than those of the JEC-FDTD method and the ADE-FDTD method, especially in the higher frequency range. Meanwhile, it is also found that larger errors occur in the stopband of the transmission coefficient for the RCP wave. Figure 13 shows that larger errors occur from 13 GHz to 38 GHz for the FDTD methods, and from 13 GHz to 27 GHz for the proposed method. This means that the low computational error range of the presented PITD method is larger than the FDTD methods.
The CPU time of the three methods are also recorded. The execution time of the JEC-FDTD method, the ADE-FDTD method, and the presented method are 8.50 s, 8.09 s, and 4.52 s, respectively. It can be concluded that a larger simulated time step leads to an appreciably reduction of the CPU time.

2-D Magnetized Plasma Filled Cavity
The second example is a 2-D cavity filled with the magnetized plasma as shown in Figure 14. According to Figures 11, 12a and 13, it should be noted that the solutions of the proposed PITD method is more accurate than those of the JEC-FDTD method and the ADE-FDTD method, especially in the higher frequency range. Meanwhile, it is also found that larger errors occur in the stopband of the transmission coefficient for the RCP wave. Figure 13 shows that larger errors occur from 13 GHz to 38 GHz for the FDTD methods, and from 13 GHz to 27 GHz for the proposed method. This means that the low computational error range of the presented PITD method is larger than the FDTD methods.
The CPU time of the three methods are also recorded. The execution time of the JEC-FDTD method, the ADE-FDTD method, and the presented method are 8.50 s, 8.09 s, and 4.52 s, respectively. It can be concluded that a larger simulated time step leads to an appreciably reduction of the CPU time.

2-D Magnetized Plasma Filled Cavity
The second example is a 2-D cavity filled with the magnetized plasma as shown in Figure 14.
The main computing region is divided into 20 × 20 cells with a space step 75 µm. The time-step size of the ADE-FDTD method is ∆t = 0.1 ps. For the proposed PITD method, the time-step size is 6∆t in this example. Therefore, the simulations are executed by 3000 iterative steps in the ADE-FDTD method and 500 iterative steps in the presented PITD method. The parameters of the magnetized plasma filled in the cavity are ω p = 2π × 28.7 × 10 9 rad/s, γ = 10.0 GHz, and ω b = 3.0 × 10 11 rad/s.   Figure 15 graphs the time-domain waveforms of the electric field Ex simulated by the presented PITD method and the ADE-FDTD method, respectively. It is shown that good agreements are achieved between the two methods on the simulated waveform. Table 1 lists the numerical resonant frequencies and the execution time of the presented PITD method and the ADE-FDTD method, respectively. It can be found that the calculated resonant frequencies are also coincident, moreover, the CPU time of the presented method is at least 1/3 less than that of the ADE-FDTD method. The simulations of both the FDTD methods and the PITD method in above analysis are achieved by MATLAB and performed on Intel(R) Core(TM) i3 CPU M370 2.40 GHz PC (Intel, Santa Clara, CA, USA).    Figure 15 graphs the time-domain waveforms of the electric field E x simulated by the presented PITD method and the ADE-FDTD method, respectively. It is shown that good agreements are achieved between the two methods on the simulated waveform. Table 1 lists the numerical resonant frequencies and the execution time of the presented PITD method and the ADE-FDTD method, respectively. It can be found that the calculated resonant frequencies are also coincident, moreover, the CPU time of the presented method is at least 1/3 less than that of the ADE-FDTD method. The simulations of both the FDTD methods and the PITD method in above analysis are achieved by MATLAB and performed on Intel(R) Core(TM) i3 CPU M370 2.40 GHz PC (Intel, Santa Clara, CA, USA).   Figure 15 graphs the time-domain waveforms of the electric field Ex simulated by the presented PITD method and the ADE-FDTD method, respectively. It is shown that good agreements are achieved between the two methods on the simulated waveform. Table 1 lists the numerical resonant frequencies and the execution time of the presented PITD method and the ADE-FDTD method, respectively. It can be found that the calculated resonant frequencies are also coincident, moreover, the CPU time of the presented method is at least 1/3 less than that of the ADE-FDTD method. The simulations of both the FDTD methods and the PITD method in above analysis are achieved by MATLAB and performed on Intel(R) Core(TM) i3 CPU M370 2.40 GHz PC (Intel, Santa Clara, CA, USA).    In conclusion, according to the numerical experiments above, the efficiency of the modified PITD method in this paper is higher than the algorithms based on the FDTD scheme for modeling the magnetized plasma. Meanwhile, the solutions of the magnetized plasma slab validate that the presented method is more accurate than the JEC-FDTD method and the ADE-FDTD method, especially in the high frequency range.

Conclusions
Based on both the auxiliary differential equation and the PI technique, a modified PITD method has been proposed for modeling the EM wave propagation through magnetized plasma in this paper. The analyses of the numerical stability and dispersion are discussed respectively, and the superior performance of the proposed method has been confirmed numerically. It is found that the numerical stability criterion of the proposed method is much looser than the CFL stability condition of the FDTD methods so as to increase the allowable simulated time step, and the numerical dispersion error and the dissipation error are almost invariant when the time step is increased. The numerical results validate the efficiency and accuracy of the presented algorithm. In the numerical experiments above, the proposed method can use a time-step size much larger than the value allowed by the CFL limit of the FDTD method which leads to a reduction of the CPU time in the simulation. Meanwhile, the accuracy performance of the presented PITD method is better than the JEC-FDTD method and the ADE-FDTD method, especially in higher frequency range. In conclusion, the proposed algorithm is a strong numerical tool to solve the EM wave problems in the magnetized plasma.