Frequency-Dependent FDTD Algorithm Using Newmark’s Method

According to the characteristics of the polarizability in frequency domain of three commonmodels of dispersivemedia, the relation between the polarization vector and electric field intensity is converted into a time domain differential equation of second order with the polarization vector by using the conversion from frequency to time domain. Newmark βγ difference method is employed to solve this equation.The electric field intensity to polarizability recursion is derived, and the electric flux to electric field intensity recursion is obtained by constitutive relation. Then FDTD iterative computation in time domain of electric and magnetic field components in dispersive medium is completed. By analyzing the solution stability of the above differential equation using central difference method, it is proved that this method has more advantages in the selection of time step. Theoretical analyses and numerical results demonstrate that this method is a general algorithm and it has advantages of higher accuracy and stability over the algorithms based on central difference method.

In the case of nondispersive medium, the spatial discretization of Maxwell curl equation can be determined by the requirement from numerical dispersion to stability, and the time discrete can be determined by Courant stability condition.So the calculation involves just the first order partial differential equation.The time domain expression of constitutive relation often can be written as a second order partial differential equation in the case of dispersive medium.To obtain an accurate solution of the partial differential equation, time step must be very small.In this paper, the general differential equation of dispersive medium is given and the Newmark method [17][18][19] is introduced to the treatment of the constitutive relation in dispersive medium.This new method has a higher accuracy and stability than the central difference method.The detailed process is as follows.
According to the polarizability characteristics in frequency domain of three common models (Debye, Drude, and Lorentz models) of dispersive medium, the relation in frequency domain between P and E is converted into a second order differential equation (has a general form to the three common models) by using the conversion  → / from frequency to time domain.Then the Newmark  method which is used in FETD is employed to solve the equation.
Hence the E → P recursion is derived, and the D → E recursion is obtained by constitutive relation.

2
International Journal of Antennas and Propagation

The General Newmark-FDTD Method
The electric constitutive relation is where () and () are permittivity coefficient and polarizability in frequency domain, respectively. ∞ is the relative permittivity when frequency is infinite.
The FDTD iterative computation of E → H → D can be accomplished by the discretization of (1) according to the standing Yee cell.And it is needed to deal with constitutive relation in the iterative computation of D → E.

The United Expression of Polarizability Function of Com-
mon Dispersive Medium Models.Three common dispersive medium models are Debye model, Drude model, and Lorentz model [2].Their polarizabilities can be written as [20] where  is numbers of poles,  represents the th pole, Δ  =  , −  ∞ ,  , is the relative permittivity at zero frequency,  0, is the relaxation time of pole, ] , is collision frequency,   is Drude frequency, and  0, is the natural frequency of oscillator.Let () of the three dispersive medium models in (3) can be written as a united form: where (6)

Constitutive Relations in Time Domain and Its Newmark's
Solution.Polarization vector P is introduced: The constitutive relation equation ( 2) can be rewritten as Let Then from ( 3) and ( 7) we know that From ( 4) and ( 9), we have After substituting ( 5) into (11) and using conversion  → / from frequency to time domain, we have E() in the above differential equation with second order time derivative is equivalent to a excitation source.In the year of 1959, a time-stepping algorithm for solving the differential equations like (12) is given by Newmark according to the time derivative difference approximate solution of Taylor series and the mean value theorem.This algorithm provides higher accuracy and stability [17].By the deriving method of Zienkiewicz, solution of ( 12) is [18] where Equation ( 13) is called two-step or Newmark  method [19].The calculation of P +1  involves the values of P   and P −1  .The coefficients in (13) include two variables  and .The stability and accuracy of (12) when selecting different  and  are also discussed in [19] in detail.Usually, 0 ≤  ≤ 1 and 0 ≤  ≤ 1/2.

International Journal of Antennas and Propagation 3
To assure the accuracy, it is needed to choose an appropriate  and  ≥ 0.25(0.5 + ) 2 .According to [19], unconditional convergence can be obtained by taking  = 0.25,  = 0.5.The time-stepping calculation from E +1 to P +1  can be finished by (13).
After converting the constitutive relation equation ( 8) into discrete time domain, we have Substituting ( 13) into ( 15), we know that Then the time-stepping calculation from D +1 , P   to E +1 can be accomplished.
In conclusion, the general Newmark-FDTD method in dispersive medium can be reduced as follows.
(1) Use the difference discrete expressions in curl (1) to get E → H → D.

The Advantages of Algorithm in This Paper
For the second order differential equations like After using central difference method, we have The time domain solution stability of the above differential equation can be analyzed by using the  transform.To a time series of () ( = 0, 1, 2...), its  transform is defined as () is the  transform of the sampling sequence (). transform is proposed to (18); we know that According to the final value theorem, time domain function () ( = 0, 1, 2...) has time stable solutions when the poles of () located in a unit circle on complex plane in  domain.An analysis of the solution stability of second order homogeneous differential equation without first order derivative is made in [21].In this paper, the case of including first order derivative is analyzed.Equation ( 20) is considered as a quadratic equation with one unknown using  as dependent variable.The equation can be solved and let −1 ≤  ≤ 1; we know that The corresponding cubic equation with one unknown to the above inequality has three solutions: Let , ,  be all positive quantities.The derivative of the cubic equation with one unknown is 8 + 8 √  at Δ 3 , which is an increasing function satisfying the requirements of the inequality.Hence, when central difference method is used to solve (17), the time-stepping is limited to According to (23), the requirements of time-stepping when central difference method is used to solve the second order time domain differential equations corresponding to the three models of dispersive medium in (6) are That is to say, for Debye medium, corrected results cannot be obtained whatever the value of the time-stepping takes when central difference method is used to solve (17).For Drude medium, corrected results can be obtained whatever the value of the time-stepping takes.But for Lorentz medium, the solution stability is relevant to its resonance frequency at poles.However, the solution of ( 12) is unconditionally stable when using Newmark difference algorithm.Hence, the selection of time-stepping is more flexible.Meanwhile, a higher accuracy is obtained because the Newmark  method has a combination of Taylor expansions and means value theorem is employed.
To analyze the difference between using Newmark  method and central difference method, the results obtained from these two methods and Runge-Kutta solutions [20] are given.Figure 1 is the absolute error comparing the classical Runge-Kutta solutions of using these two methods, in which the solid line and the circle are the result from Newmark  and central difference method, respectively.As it is shown in in (17).The theoretical analysis and numerical results show that the algorithm presented in this paper has higher steady state and calculation accuracy than central difference method.

Numerical Results
To demonstrate the correctness of our algorithm, the backscatter of three different dispersive media sphere is given and compared with Mie series solution.Let  = 0.25,  = 0.5, Δ = /2, where  is spatial discrete grid-scale of Yee cell.Consider where  = 60Δ and  0 is the pulse peak time.
Example 1 (the backscatter of Debye sphere).Let the radius of the sphere be 0.25 m and filled with single pole Debye medium.The complex relative permittivity is where  ,1 = 1.16,  ∞ = 1.01,  0,1 = 6.497 × 10 −10 s.The backscatter RCS is shown in Figure 2. The solid line and circle are the results from this method and Mie series solution, respectively.It is shown that these two results are in good  Example 2 (the backscatter of plasma sphere).The radius of the sphere is 3.75 mm and filled with single pole Drude medium.The complex relative permittivity is where  1 = 1.8 × 10 11 rad/s, ] ,1 = 2.0 × 10 10 Hz.The backscatter RCS is shown in Figure 3.The solid line and circle are the result from this method and Mie series solution, respectively.It is shown that these two results are in good agreement with each other.Consider  = 5.0 × 10 −2 mm in FDTD calculation.

Conclusions
In this paper, a general Newmark-FDTD algorithm is given to deal with the electromagnetic problems in dispersive medium.This new method combines the Newmark  difference method which is widely used in FETD calculation in dispersive medium.Theoretical analyses and numerical results demonstrate that this algorithm has advantages of higher accuracy and stability over the algorithms based on central difference method.

Figure 1 :
Figure 1: Comparison of stability and calculation error between central difference and Newmark method.