# Newmark-FDTD Formulation for Modified Lorentz Dispersive Medium and Its Equivalence to Auxiliary Differential Equation-FDTD with Bilinear Transformation

Authors: Hongjin Choi; Jeahoon Cho; Yong Bae Park; Kyung-Young Jung
Journal: International Journal of Antennas and Propagation
 (2019)

---

## Abstract

The finite-difference time-domain (FDTD) method has been popularly utilized to analyze the electromagnetic (EM) wave propagation in dispersive media. Various dispersion models were introduced to consider the frequency-dependent permittivity, including Debye, Drude, Lorentz, quadratic complex rational function, complex-conjugate pole-residue, and critical point models. The Newmark-FDTD method was recently proposed for the EM analysis of dispersive media and it was shown that the proposed Newmark-FDTD method can give higher stability and better accuracy compared to the conventional auxiliary differential equation- (ADE-) FDTD method. In this work, we extend the Newmark-FDTD method to modified Lorentz medium, which can simply unify aforementioned dispersion models. Moreover, it is found that the ADE-FDTD formulation based on the bilinear transformation is exactly the same as the Newmark-FDTD formulation which can have higher stability and better accuracy compared to the conventional ADE-FDTD. Numerical stability, numerical permittivity, and numerical examples are employed to validate our work.


---

## Body

## 1. Introduction

The finite-difference time-domain (FDTD) method has been widely utilized to analyze various electromagnetic wave (EM) problems owing to its simplicity, robustness, and accuracy [1–3]. It is of great importance to choose well-fitted dispersion model to obtain accurate results in dispersive media. A variety of dispersion models were proposed such as Debye, Drude, Lorentz, quadratic complex rational function (QCRF), complex-conjugate pole-residue, and critical point models [4–14]. Note that the modified Lorentz dispersion model can unify the aforementioned dispersion models and has more degree of freedom than Debye, Drude, and Lorentz models [15, 16].Recently, the Newmark time-stepping algorithm was applied to the dispersive FDTD modeling of Debye, Drude, Lorentz, and QCRF dispersion models [17, 18]. The authors showed that the Newmark-FDTD method can lead to simple arithmetic implementation, higher stability, and better accuracy compared to the conventional auxiliary differential equation- (ADE-) FDTD method [19]. In this work, the Newmark-FDTD method is successfully extended to the modified Lorentz dispersion model. It is found that the formulation of the ADE-FDTD method based on the bilinear transformation (BT) [20–22] leads to the same formulation of the Newmark-FDTD method. Moreover, the ADE-FDTD formulations with BT for Debye, Drude, Lorentz, and QCRF dispersion models are also considered and it will be shown that the resulting FDTD formulations are equivalent to the Newmark-FDTD counterparts.The remainder of this paper is organized as follows. The Newmark time-stepping algorithm is reviewed and then the Newmark-FDTD formulation is derived for modified Lorentz medium. Numerical stability and numerical permittivity of the Newmark-FDTD formulation are discussed and the equivalence of the Newmark-FDTD formulation to the ADE-FDTD formulation based on the BT is also addressed. In the next section, numerical examples involving homogenous one-dimensional (1D) structure and three-dimensional (3D) plasmonic nanosphere are used to validate our work.


## 2. Formulations

Before proceeding with the Newmark-FDTD method, it is worth reviewing the Newmark time-stepping algorithm briefly [23, 24]. Let us consider the second-order differential equation of the form(1)Mx′′+Cx′+Kx=f,where x denotes the displacement [24]. The Taylor expansion of the velocity (x′) about (n+1)Δt is(2)xn+1′=xn′+Δtxn′′+Δt22!xn′′′+…,where xn denotes the value of x at nΔt. For sufficiently smooth functions, the above expansion can be truncated as(3)xn+1′=xn′+Δtxq′′,where n≤q≤n+1. Using the linear interpolation, (3) can be expressed as(4)xn+1′=xn′+1-γ∆txn′′+γ∆txn+1′′,where 0≤γ≤1. Similarly, the Taylor expansion of the displacement about (n+1)Δt can be written as(5)xn+1=xn+∆txn′+∆t22!xr′′,where n≤r≤n+1. Employing the linear interpolation to the above equation, we have(6)xn+1=xn+∆txn′+1-2β∆t22!xn′′+2β∆t22!xn+1′′,where 0≤2β≤1. When the governing equation (1) is expressed at time step n+1,n, and n-1, the following equations can be obtained:(7)Mxn+1′′+Cxn+1′+Kxn+1=fn+1,(8)Mxn′′+Cxn′+Kxn=fn,(9)Mxn-1′′+Cxn-1′+Kxn-1=fn-1.Similar procedure to (4) and (6) is applied to the value of nΔt:(10)xn′=xn-1′+1-γ∆txn-1′′+γ∆txn′′,(11)xn=xn-1+∆txn-1′+1-2β∆t22!xn-1′′+2β∆t22!xn′′.There are three nonderivative terms (xn-1,xn,xn+1) and six derivative terms (xn-1′,xn′,xn+1′,xn-1′′,xn′′,xn+1′′)in seven equations (see (4) and (6)–(11)). After simple mathematical manipulation, we can finally obtain the update equation without the derivative terms:(12)Kβ∆t2+Cγ∆t+Mxn+1+K0.5-2β+γ∆t2-C2γ-1∆t-2Mxn+K0.5+β-γ∆t2+Cγ-1∆t+Mxn-1=β∆t2fn+1+0.5-2β+γΔt2fn+0.5+β-γ∆t2fn-1.Now, let us derive the Newmark-FDTD method for the modified Lorentz dispersion model. The modified Lorentz model [15, 16] can be expressed as εr(ω)=εr,∞+χ(ω), where(13)χω=a0+a1jωb0+b1jω+b2jω2and εr,∞ is a relative permittivity at the infinite frequency. Note that the update equation of magnetic field can be implemented by using standard central difference scheme (CDS) to Faraday’s law [1]: (14)∇×E=-μ0∂H∂t.Ampere’s law can be written as (15)∇×H=ε∞∂E∂t+∂P∂t,where ε∞=ε0εr,∞, ε0 is the permittivity in the free space, and P(ω)=ε0χ(ω)E(ω) is the constitutive relation in the frequency domain. Inverse Fourier transform (IFT) is applied to the constitutive relation: (16)b0Pt+b1∂Pt∂t+b2∂2Pt∂t2=a0ε0Et+a1ε0∂Et∂t=Wt,where a temporary variableW is introduced to apply the Newmark time-stepping algorithm.Therefore, we have two sets of differential equations:(17)b0Pt+b1∂Pt∂t+b2∂2Pt∂t2=Wt,(18)a0ε0Et+a1ε0∂Et∂t=Wt.By applying the Newmark time-stepping algorithm to the above equations, the solution in terms ofP andEcan be obtained as follows:(19)CaPn+1+CbPn+CcPn-1=CdEn+1+CeEn+CfEn-1,where (20)Ca=b0β∆t2+b1γ∆t+b2,Cb=b00.5-2β+γ∆t2-b12γ-1∆t-2b2,Cc=b00.5+β-γ∆t2+b1γ-1∆t+b2,Cd=ε0a0β∆t2+a1γ∆t,Ce=ε0a00.5-2β+γ∆t2-a12γ-1∆t,Cf=ε0a00.5+β-γ∆t2+a1γ-1∆t.Therefore, theE field can be updated by inserting Pn+1 of (19) into the CDS version of Ampere’s law (15):(21)Caε∞+CdCa∆tEn+1=Caε∞-CeCa∆tEn-CfCa∆tEn-1+∇×Hn+1/2+Cb+CaCa∆tPn+CcCa∆tPn-1.In what follows, we choose β=0.25 and γ=0.5 the same as in [17, 18] because this Newmark-FDTD can have higher stability and better accuracy compared to the conventional ADE-FDTD method [18].Next, let us derive the conventional ADE-FDTD method. Toward this purpose, CDS is applied to (16):(22)b0Pn+b1Pn+1-Pn-12∆t+b2Pn+1-2Pn+Pn-1∆t2=a0ε0En+a1ε0En+1-En-12∆t.This equation can also be written in the same form as (19), where the coefficients are expressed as (23)Ca=0.5b1∆t+b2,Cb=b0∆t2-2b2,Cc=-0.5b1∆t+b2,Cd=ε00.5a1∆t,Ce=ε0a0∆t2,Cf=ε0-0.5a1∆t.Now, let us consider the ADE-FDTD formulation based on the BT. The BT is basically an approximate of jω≈2/∆t1-Z-1/1+Z-1 [20–22]. When the BT is applied to the constitutive relation, (16) can be expressed as(24)b0Pn+1+2Pn+Pn-14+b1Pn+1-Pn-12∆t+b2Pn+1-2Pn+Pn-1∆t2=a0ε0En+1+2En+En-14+a1ε0En+1-En-12∆t.The above equation can also be expressed as the form of (19), where the coefficients are(25)Ca=0.25b0∆t2+0.5b1∆t+b2,Cb=0.5b0∆t2-2b2,Cc=0.25b0∆t2-0.5b1∆t+b2,Cd=ε00.25a0∆t2+0.5a1∆t,Ce=ε00.5a0∆t2,Cf=ε00.25a0∆t2-0.5a1∆t.Please note that the above coefficients are exactly the same as (20) with β=0.25 and γ=0.5.It is worth comparing the numerical stability of the Newmark-FDTD formulation, the conventional ADE-FDTD formulation, and the ADE-FDTD formulation based on the BT. The numerical stability conditions for FDTD formulations can be obtained by using von Neumann method [4, 25–27].Z transform is applied to the wave equation and the constitutive relation. For the twoZ-transformed equations, the determinant must be zero to get nonzero solutions. To be stable, all the roots of the polynomial must be inside or on the unit circle. Instead of finding the roots of the polynomial, let us transform Z=r+1/r-1. By transforming theZ-domain into ther-domain, roots inside the unit circle in theZ-domain map onto the left half-plane in ther-domain. Next, Routh-Hurwitz (R-H) criterion can be used to derive the numerical stability conditions [4]. To be stable, all components of the first column on the Routh table must be equal to or larger than zero. Following the above procedure, the numerical stability conditions of the conventional ADE-FDTD formulation are as follows:(26)b0≥0,b1≥0,Q∆t2+b1εr,∞ν24b2-b0∆t2≥0,Qa1+b1εr,∞∆t2+4a1b1b2εr,∞ν2≥0,4b2-b0∆t21-ν2εr,∞-a0∆t2≥0,where Q=a0b1-a1b0,ν2=(c∞∆t)2∑α=x,y,zsin2(kα~∆α/2)/∆α2, c∞=1/μ0ε∞, and kα~ are the numerical wavenumbers in the α direction. The numerical stability conditions of ADE-FDTD formulation based on the BT are the same as those of the Newmark-FDTD formulation, since both methods have the same formulation. Their numerical stability conditions are(27)b0≥0,b1≥0,Q∆t2+4b1b2εr,∞ν2≥0,Qa1+b1εr,∞1-ν2∆t2+4a1b1b2εr,∞ν2≥0,b21-ν2≥0.Next, the numerical permittivity can be derived inserting the numerical solutionP0ejωn∆t and E0ejωn∆t into the constitutive relation [24, 28]. The numerical permittivity of the ADE-FDTD formulation can be obtained as follows:(28)εr~ω=εr,∞+a0~+a1jω~b0~+b1jω~+b2jω~2,where ω~=2tan(ω∆t/2)/∆t,a0~=a0/cos⁡2(ω∆t/2), and b0~=b0/cos⁡2(ω∆t/2). The numerical permittivity of the ADE-FDTD formulation based on the BT is the same as that of the Newmark-FDTD formulation: (29)εr~ω=εr,∞+a0+a1jω~b0+b1jω~+b2jω~2.We also consider Debye, Drude, Lorentz, and QCRF dispersion models. The dispersion models are as follows:(30)εrDebyeω=εr,∞+∆ε1+jωτ,εrDrudeω=εr,∞+ω02jωγ+jω2,εrLorentzω=εr,∞+∆ε·ω02ω02+2jωδ+jω2,εrQCRFω=A0+A1jω+A2jω2B0+B1jω+B2jω2.For Debye, Drude, and Lorentz dispersion models, 1/χ(ω) can be written as M(jω)2+C(jω)+K, whereM,C, and K are constants for each model. By applying the Newmark-FDTD, one can obtain(31)SPn+1=w1Pn+w2Pn-1+u0En+1+u1En+u2En-1,where (32)S=M+0.5C∆t+0.25K∆t2,w1=2M-0.5K∆t2,w2=-M+0.5C∆t+0.25K∆t2,u0=0.25ε0∆t2,u1=0.5ε0∆t2,u2=0.25ε0∆t2.Note that this update equation is the same as the ADE-FDTD formulation based on the BT. For QCRF dispersion model, the constitutive relation between electric flux densityD and electric fieldE is simpler. The update equation for the Newmark-FDTD method is the same as the ADE-FDTD formulation based on the BT and it can be expressed as(33)0.25B0∆t2+0.5B1∆t+B2Dn+1+0.5B0∆t2-2B2Dn+0.25B0∆t2-0.5B1∆t+B2Dn-1=0.25A0ε0∆t2+0.5A1ε0∆t+A2ε0En+1+0.5A0ε0∆t2-2A2ε0En+0.25A0ε0∆t2-0.5A1ε0∆t+A2ε0En-1.All things considered, the Newmark-FDTD formulation with β=0.25 and γ=0.5 is the same as the ADE-FDTD formulation based on the BT.


## 3. Results and Discussion

In this section, numerical examples are used to validate our study. First, we consider human blood from 300 MHz to 3 GHz. The modified Lorentz parameters are extracted by using the particle swarm optimization [29, 30] for the complex relative permittivity data [31]. They are a0=6.9379·1021,a1=1.5057·1012,b0=6.1637·1018,b1=4.5425·1010,b2=1, and εr,∞=31.1662. The space step size ∆z=1mm and time step size ∆t=Cn∆z/c∞, where Cn is a classic Courant number [1] and Cn= 0.99 is used in this work. The relative errors of the numerical permittivity of the three FDTD formulations are shown in Figure 1. The accuracy of the ADE-FDTD with BT is the same as the Newmark-FDTD and they have higher accuracy than the ADE-FDTD. This result can be explained from the fact that thenumerical coefficients a0~ and b0~ are involved in the ADE-FDTD method in addition to thenumerical angular frequency ω~, differently from both the Newmark-FDTD and the ADE-FDTD with BT (see (28)-(29)).Figure 1
Relative error of the numerical permittivity.Next, actual FDTD simulations are performed. A sinewave with the frequency of 300 MHz is excited in 1D homogenous modified Lorentz medium. The 10-layer perfectly matched layer (PML) [1] is used to truncate the computational domain of 10,000 FDTD cells. First, we consider the aforementioned modified Lorentz parameters that satisfy the numerical stability conditions of (26) and (27) for all possible values of ν. As expected, all FDTD simulations are stable for 10,000 time steps in Figure 2(a). As shown in Figure 3(a), all roots of the stability polynomial in theZ-domain are inside or on the unit circle for all cases. The memory requirement and the central processing unit (CPU) time of all three FDTD formulations are the same as 6.6 MB and 39.6 s, respectively, because the number of auxiliary variables and the number of arithmetic operations are completely the same for all three formulations. Now, only b2 is changed to 0.8, where the stability conditions of (26) are not satisfied. For the ADE-FDTD method, instability can be found in Figure 2(b). It can be observed that some roots of stability polynomial in the Z-domain are outside the unit circle for the ADE-FDTD method but all roots of stability polynomial in theZ-domain are inside or on the unit circle for the others, as shown in Figure 3(b). As can be seen in these FDTD simulations, both the Newmark-FDTD and the ADE-FDTD with BT have the same accuracy and they are better than the ADE-FDTD in terms of the numerical stability.Figure 2
1D FDTD simulations for different sets of modified Lorentz parameters. (a) All FDTD formulations are stable. (b) Only ADE-FDTD is unstable. Legends are the same as Figure1.
(a) (b)Figure 3
Root locus of the stability polynomial in the Z-domain for Figure1. Legends: Black circle: Newmark-FDTD/ADE-FDTD with BT; blue-gray square: ADE-FDTD.
(a) (b)As a final example, a 3D Ag nanosphere with the radius of 40 nm is considered. In this case, the parameters in [32] are used, a0=1.9136·1032,a1=0,b0=0,b1=2.7362·1013,b2=1, and εr,∞=3.7. The x-polarized Gaussian-modulated sinewave plane wave is excited and the computational domain is terminated by a 10-layer complex frequency shifted- (CFS-) PML [33, 34]. The FDTD computational domain comprised 260x260x260 cells and the total number of time marching steps is 40,000. To compare FDTD results with Mie theory [35], we calculate the magnitude of Ex at the centre of the nanosphere normalized by the magnitude of the incident electric field. As shown in Figure 4, the Newmark-FDTD is the same as the ADE-FDTD with BT and they have higher accuracy than the ADE-FDTD. For all 3D FDTD simulations, the memory requirement is 4.5 GB and the CPU time is 6 hours.Figure 4
3D FDTD simulations for an Ag nanosphere with the radius of 40 nm. (a) Spectral response. (b) Relative error.
(a) (b)


## 4. Conclusions

In this work, the Newmark-FDTD method is applied to modified Lorentz dispersion model which can systematically unify various existing dispersion models. It is found that the Newmark-FDTD is equivalent to the ADE-FDTD with the BT in terms of update formulation, numerical stability, and numerical accuracy. In addition, it is figured out that both methods can yield better accuracy and higher stability against the conventional ADE-FDTD method. Moreover, Debye, Drude, Lorentz, and QCRF dispersion models are extended to the Newmark-FDTD method and its equivalence to the BT-based ADE-FDTD counterpart is also observed. Numerical permittivity and the 1D blood example are used to illustrate that both the Newmark-FDTD method and the ADE-FDTD method based on the BT are better than the conventional ADE-FDTD method in terms of numerical accuracy and numerical stability. A further numerical example involving 3D plasmonic nanosphere is presented to demonstrate that both the Newmark-FDTD simulation and the BT-based ADE-FDTD simulation are in good agreement with the Mie solution and they are superior to the conventional ADE-FDTD method.



---

*Source: 4173017-2019-06-20.xml*