# An Extended Newmark-FDTD Method for Complex Dispersive Media

Authors: Yu-Qiang Zhang; Peng-Ju Yang
Journal: International Journal of Antennas and Propagation
 (2018)

---

## Abstract

Based on polarizability in the form of a complex quadratic rational function, a novel finite-difference time-domain (FDTD) approach combined with the Newmark algorithm is presented for dealing with a complex dispersive medium. In this paper, the time-stepping equation of the polarization vector is derived by applying simultaneously the Newmark algorithm to the two sides of a second-order time-domain differential equation obtained from the relation between the polarization vector and electric field intensity in the frequency domain by the inverse Fourier transform. Then, its accuracy and stability are discussed from the two aspects of theoretical analysis and numerical computation. It is observed that this method possesses the advantages of high accuracy, high stability, and a wide application scope and can thus be applied to the treatment of many complex dispersion models, including the complex conjugate pole residue model, critical point model, modified Lorentz model, and complex quadratic rational function.


---

## Body

## 1. Introduction

The finite-difference time-domain (FDTD) method has been widely used to study electromagnetic (EM) wave interaction with a wide variety of materials due to its robustness and its ability to calculate accurate broadband response via a single simulation [1, 2]. In the analysis of dispersive materials using the FDTD approach, it is required that the variation of dielectric parameters with frequency is modeled efficiently and accurately. In recent years, some novel dispersive models have been introduced, for example, the complex-conjugate pole-residue (CCPR) model [3], critical point (CP) model [4], modified Lorentz (m-Lo) model [5], and quadratic complex rational function (QCRF) model [6]. Some hybrid dispersion models, such as the Drude-Lorentz model [7, 8] and Drude-CP model [4, 9], were also usually adopted. However, the introduction of these new dispersion models and their hybrid models undoubtedly challenges the versatility of the existing FDTD methods. For the typical dispersive FDTD approaches, the recursive convolution (RC) method and Z transform (ZT) method are not convenient to deal with the high-order models such as m-Lo and QCRF. The auxiliary differential equation (ADE) method has more stringent stability conditions than the Courant-Friedrichs-Lewy (CFL) condition of Yee’s FDTD scheme [10–12].The Newmark algorithm was originally used for the numerical solution of structural dynamics and was introduced into the finite-element time-domain (FETD) method for electromagnetic simulations [13–15]. Then, it is introduced into the analysis of the dispersive FDTD computation in 2014 for the unified treatment of the Debye, Drude, and Lorentz medium [16–18]. In this paper, the Newmark-FDTD method is further extended to the high-order model so that it can be used as a general approach for dispersive media, not only to treat with Debye, Drude, and Lorentz media, but also to deal with complex dispersion models including CCPR, CP, m-Lo, and QCRF media. The scope of application of the algorithm is extended greatly. Then, the computational accuracy and stability are investigated from the two aspects of theory and numerical computation, and higher accuracy and stability are shown.


## 2. Extended Newmark-FDTD Method

For the linear electric dispersive medium, the Maxwell’s curl equation can be written as(1)∇×H=∂D∂t+σE,(2)∇×E=−∂B∂t,where E, D, H, and B are the electric field strength, electric displacement, magnetic field strength, and magnetic flux density, respectively. σ is the electric conductivity. For simplicity, consider only the electric dispersion medium. The constitutive relation in a frequency domain can be characterized as
(3)B=μH,(4)Dω=εωEω=ε0ε∞+χωEω,where μ is the magnetic permeability. ε0 and ε∞ are the relative electric permittivity in free space and at the infinite frequency, respectively. εω and χω are the electric permittivity and electric susceptibility in frequency domain, respectively. In this paper, the electric susceptibility is adopted in the following rational function form:
(5)χω=∑q=1Qχqω=∑q=1Qb0,q+b1,qjω+b2,qjω2a0,q+a1,qjω+a2,qjω2,where ai,q and bi,q (i = 0, 1, 2) are both real constants. In this way, the three typical kinds of dispersive models, that is, Debye, Drude, and Lorentz, and some of the complex dispersive models, that is, CCPR, CP, m-Lo, and QCRF, can be viewed as special cases of (5). The higher degrees of freedom and wider range of application are possessed.The time-stepping equation ofH can be simply obtained by applying the Yee discrete scheme to (2). However, to obtain the time-stepping equation of E from (1), we must firstly solve (4) in a time-stepping manner. For this reason, the polarization vector P in the frequency domain is introduced
(6)Pω=χωEω,and let Pq=χqωEω; thus, (4) can further be represented as
(7)Dω=ε0ε∞Eω+ε0∑q=1QPqω.Substituting (5) into P, we can obtain
(8)a0,q+a1,qjω+a2,qjω2Pqω=b0,q+b1,qjω+b2,qjω2Eω.Converting (8) into a time domain by the relation jω→∂/∂t, we arrive at a two-order differential equation
(9)a2,qd2Pqtdt2+a1,qdPqtdt+a0,qPqt=b2,qd2Etdt2+b1,qdEtdt+b0,qEt.It can be observed that the right side of (9) not only contains Et but also contains the first-order derivative and second-order derivative of Et so that it cannot be solved directly using the Newmark algorithm employed in [16, 17]. For this reason, a temporary variable W is introduced and (9) is further represented as
(10)a2,qd2Pqtdt2+a1,qdPqtdt+a0,qPqt=W,b2,qd2Etdt2+b1,qdEtdt+b0,qEt=W.Then, the Newmark algorithm is applied to (10) with a unified time step Δt and to obtain higher stability, γ=0.5 and β=0.25 are adopted [13, 14]. Thus, we can obtain
(11)c0,q⋅Pqn+1−c1,q⋅Pqn+c2,q⋅Pqn−1=Δt24Wn+1+Δt22Wn+Δt24Wn−1,(12)d0,q⋅En+1−d1,q⋅En+d2,q⋅En−1=Δt24Wn+1+Δt22Wn+Δt24Wn−1,where
(13)c0,q=a2,q+12a1,qΔt+14a0,qΔt2,(14)c1,q=2a2,q−12a0,qΔt2,(15)c2,q=a2,q−12a1,qΔt+14a0,qΔt2.d0,q, d1,q, and d2,q are similar to c0,q, c1,q, and c2,q. It only needs replacing a0,q, a1,q, and a2,q in (13), (14), and (15) with b0,q, b1,q, and b2,q.Due to be a unifiedΔt, W in (11) and (12), the value of W at n+1, n, and n−1 must be the same, so we can obtain the Pq update formula
(16)Pqn+1=1c0,qc1,q⋅Pqn−c2,q⋅Pqn−1+d0,q⋅En+1−d1,q⋅En+d2,q⋅En−1.Finally, substituting (16) and (7) into the discretized form of (1), the FDTD-updating formula of E is given
(17)ε0ε∞+ε0∑q=1Qd0,qc0,q+σΔt2En+1=ε0ε∞+ε0∑q=1Qd1,qc0,q−σΔt2En+Δt⋅∇×H−ε0∑q=1Qd2,qc0,qEn−1+ε0⋅∑q=1Q1−c1,qc0,q⋅Pqn+∑q=1Qc2,qc0,q⋅Pqn−1.In this way, a complete FDTD time-domain iterative computation is constituted by (17) and (16) and the standard FDTD formula of H corresponding to (2). The FDTD computation can be realized.


## 3. Analysis of Accuracy and Stability

In the FDTD method, Maxwell’s curl equations are replaced with a different equation and the dispersion relation in the discrete form is thus quite different from that in the continuous form. For the dispersive FDTD method, the accuracy of the FDTD method is related also to a discrete scheme applied to the dispersive model of the media in the time domain. So, the accuracy of FDTD method can clearly be understood by investigating the errors between the numerical dispersive model of media and analytical model of media. The smaller the errors, the higher the computational accuracy [19]. The accuracy of the extended Newmark-FDTD algorithm is discussed below. Moreover, for simplicity and without loss of generality, let Q=1 and ε∞=0; the permittivity of dispersive media becomes a rational fraction form similar to (5).To investigate the accuracy of the extended FDTD method, thez-domain electric susceptibility is given by applying the z transform to (16)
(18)χ~z=d0,qz−d1,q+d2,qz−1c0,qz−c1,q+c2,qz−1.Substituting c0,q, c1,q, c2,q, d0,q, d1,q, d2,q, and z=ejωΔt into (18) and then rearranging the resulting equation, we can arrive at
(19)χ~ejω=b0,qΔt2cos21/2ωΔt+jb1,qΔtsinωΔt+4b2,qsin21/2ωΔta0,qΔt2cos21/2ωΔt+ja1,qΔtsinωΔt+4a2,qsin21/2ωΔt.To derive a more compact form, the following double-angle formula of trigonometric identity is applied to (19):
(20)sinωΔt=2sin12ωΔtcos12ωΔt,tan12ωΔt=sin12ωΔtcos12ωΔt.The numerical electric susceptibility corresponding to (5) is derived as follows:
(21)χ~ω~=b0+b1jω~+b2jω~2a0+a1jω~+a2jω~2,where
(22)ω~=2Δttan12ωΔt.Comparing (21) with (5), it is shown that the forms of numerical electric susceptibility χ~ω~ and electric susceptibility χω are exactly the same, apart from replacing ω with ω~.By expanding (22) into the power series, we can obtain
(23)2tanωΔt2=ωΔt+ωΔt312+OωΔt5≈ωΔt.The error between ω~ and ω is not higher than ωΔt3. For practical FDTD computation, ωΔt restricted by CFL stable condition is always very small. Hence, χ~ω~ is very close to χω. In comparison with the ADE method given in [10, 19], the extended Newmark-FDTD method possesses higher accuracy.Then, the stability of the extended FDTD method is investigated, which is of great importance to fully employ this method for general dispersive media. The stability of the dispersive FDTD method depends not only on the dispersive model but also on the discretization scheme [11]. It can be carried out by the Von Neumann method combined with the Hurwitz-Routh (R-H) criterion [20]. Consider a time-dependent wave equation in a source-free homogeneous dispersive medium, we have
(24)z−12D~+4zε0v2E~=0,where
(25)v2=c0Δt2∑s=x,y,zsin2k~sΔs/2Δs2.In (24) and (25), D~ and E~ are complex amplitudes of field quantities, k~ss=x,y,z is a numerical wavenumber in s direction, Δs is the space step in s direction, and c∞=1/μ0ε∞. Substituting (18) into (4), we have
(26)ε~z=D~E~=ε0d0z2−d1z+d2c0z2−c1z+c2.The stability polynomial in z-plane can be given as follows from (24) and (26):
(27)SDEz=d0z4−d1+2d0−4c0v2z3+d2+2d1+d0−4c1v2z2−2d2+d1−4c2v2z+d2.It is known from the knowledge of signal processing that the root of SDEz must lie in the unit circle on the z-plane to make the FDTD algorithm stable. In order to avoid solving directly the root of SDEz, we can transform (27) into the s-plane by the relation z=s+1/s−1, then substituting c0, c1, c2, d0, d1, and d2, the following stability polynomial in the s-plane is obtained:
(28)SDEs=4v2a0Δt2s4+8v2a1Δts3+4b0−v2a0Δt2+16v2a2s2+8b1−v2a1Δts+16b2−v2a2.In this way, the problem of judgment of whether the root of SDEz is in the unit circle on the z-plane is converted into the problem of whether the root of SDEs lies in the left of the s-plane, which can be calculated according to the R-H criterion with the Routh table. Finally, the following effective stability conditions are obtained:
(29)a0≥0,(30)a1≥0,(31)a1b0−a0b1Δt2+4v2a1a2≥0,(32)b2−a2ε0v2≥0.Taking (32) into account, for a2>0, the last numerical stability condition of the extended Newmark-FDTD method can be written in terms of Δt:
(33)Δt≤b2a21c0∑s=x,y,z1Δs2−1/2.A comparison of the extended Newmark-FDTD algorithm with the double average scheme (DAS) of the ADE method as well as the Mobius transformation (MT) method indicates that a similar time-domain update formula is presented except for the coefficient difference with a constant in the three dispersive FDTD algorithms, although these three formulations are based on different approximation principles [11, 21–23]. The aforementioned numerical stability condition expressed by (33) is also the same as the stability condition of DAS and MT and is more relaxed than the stability condition of the direct scheme (DS) of the ADE algorithm [12]. Therefore, the extended Newmark-FDTD method also has the best stability [24, 25]. However, the derivation of the extended Newmark scheme is much simpler than that of the ADE and MT method, as the problem of the discrete scheme and time synchronism in the ADE algorithm as well as the tedious derivation in the MT algorithm can be avoided [22, 26].


## 4. Numerical Verification

To verify the advantages of extended Newmark-FDTD method, the numerical analysis of accuracy and stability of this method is further conducted.Firstly, to verify the computational accuracy of the extended Newmark-FDTD algorithm, consider a 1-D human dry skin tissue modeled in the frequency range of 400 MHz–3 GHz by a QCRF with the coefficients ofQ=1, b0=188.85, b1=7.49×10−8, b2=1.04×10−18, a0=1.0, a1=1.80×10−9, and a2=4.38×10−20 [19]. The error is measured using the relative error between ε~ω~ and εω defined as
(34)error=εω−ε~ω~εω.The results atdt=dtmax and dt=0.5⋅dtmax are shown in Figure 1, where dtmax=Δz/c0 and Δz are one tenths of the minimum wavelength in incident wave. As a comparison, the relative error of the direct scheme (DS) of the ADE method is also given. It is seen from Figure 1 that both the relative errors of the extended Newmark algorithm and ADE algorithm are very small and decrease with the decrease of the discrete interval Δt. A comparison of the two methods shows that the extended Newmark algorithm has a smaller magnitude of relative error than the ADE method, so it has higher computational accuracy.Figure 1
The relative error ofε~ω~ and εω for the extended Newmark algorithm.The stability of the extended Newmark-FDTD method is also investigated by the root locus of the stability polynomial inz-plane as (27). Consider a human muscle tissue modeled in the frequency range of 400 MHz–3 GHz by QCRF as (5) with the following parameters [21]: Q=1, ε∞=0, b0=455.72, b1=2.50×10−7, b2=1.98×10−18, a0=1, a1=4.47×10−9, and a2=7.56×10−20. Let λmin=0.0138 m and Δz=λmin/10, then the maximum time step satisfying the CFL stability limit of 1-D FDTD in free space is Δt0=Δz/c0=4.6021 ps and the time step of 1-D FDTD in human muscle tissue is Δt=v⋅Δt0, where 0<v≤vmax=b2/a2=5.1177, vmax is computed from (33). Then, the root locus of the stability polynomial of an extended Newmark-FDTD implementation is computed and shown in Figure 2(a) for 0<v≤vmax. At the same time, the root locus of the stability polynomial of the direct scheme of the ADE algorithm is given in Figure 2(b).Figure 2
Root locus of the stability polynomial inz-domain.
(a)
Extended Newmark-FDTD method
(b)
Direct scheme of ADE methodIt is seen from Figure2(a) that for 0 < v < 5.1177, all the roots of the stability polynomial of the extended Newmark-FDTD implementation are inside the unit circle that shows the extended Newmark-FDTD method is stable for Δt satisfied (33). This is consistent with the theoretical analysis. On the other hand, it is seen from Figure 2(b) that one root of the stability polynomial of the DS-ADE method escapes outside the unit circle for 5.0422<v<5.1177. This indicates that the direct scheme of the ADE method becomes unstable. This result is also consistent with the results from [21]. The comparison of Figures 2(a) and 2(b) shows that the stability condition of the extended Newmark scheme is more relaxed than that of the ADE direct implementation scheme.To validate further the theoretical stability results, the FDTD simulation of electromagnetic wave propagation in a 1-D domain is considered. The simulation domain is composed of a human muscle tissue slab with the same parameters as described above and sandwiched between two semi-infinite dielectrics withεr=A2/B2+1. The size of the simulation domain is limited to 8000 space cells and terminated by the convolutional perfectly matched layer (CPML) with a thickness of 10 cells. The muscle slab is located between the nodes 4000 and 5000 of the domain. A modulated Gaussian pulse given by Ex=exp−t/τ−32sin2πfct is introduced at the S point that is 10 Δz away from the left boundary of the simulation domain (where fc=1.5 GHz and τ=0.5198 ns). The electric field Ex can be observed at point O that is 3000 Δz away from the left boundary of the simulation domain [21]. The simulation domain geometry is shown in Figure 3.Figure 3
Structure of simulation domain (where S and O are, resp., the source point and observation point).If the spatial step isΔz=0.00138 m, the maximum time step of the 1-D FDTD computation in free space that satisfied the Courant stability condition is Δt0=Δz/c0=4.6021 ps and the upper bound of the time step of 1-D FDTD in the human muscle tissue for the discretization of extended Newmark FDTD is Δtmax=A2/B2Δt0=5.1177Δt0=23.5521 ps obtained from (33). The maximum time step of the 1-D FDTD computation in the surrounding semi-infinite dielectrics is Δt=Δz/c0/εr=23.997 ps. In the first test, the FDTD simulation is performed with a time step limit Δtmax. This time step selection meets the stability limit in the human muscle tissue. Figure 4(a) shows the time-domain evolution of Ex as recorded at the observation point O. It can be seen that the field Ex remains stable during the whole period of time. The simulation is repeated by using Δt=1.01Δtmax=23.7877 ps that is bigger than the stability limit in the human muscle tissue. Figure 4(b) shows the time-domain evolution of Ex as recorded at observation point O. It can be seen that the field begins to increase without being bound in the late-time response, and therefore, it is unstable.Figure 4
Waveform ofEx at the observation point.
(a)
Δt=Δtmax
(b)
Δt=1.01Δtmax


## 5. Conclusion

In this paper, an extended Newmark-FDTD approach for the complex dispersive medium is presented in terms of polarizability in the form of a quadratic complex rational function, in which the time-stepping equation of the polarization vector is obtained by the Newmark algorithm used widely in finite-element time-domain (FETD) computation. The high accuracy and stability of this method are demonstrated from the two aspects of theoretical analyses and numerical results. It is also indicated that the application scope of this approach is extended greatly. In addition, this approach can be further used for the analysis of some complex dispersive models and some hybrid models, for example, CCPR, CP, m-Lo, QCRF, Drude-Lorentz, and Drude-CP.



---

*Source: 1573512-2018-04-10.xml*