#### 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.

#### 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 formwhere denotes the displacement [24]. The Taylor expansion of the velocity () about () iswhere denotes the value of at . For sufficiently smooth functions, the above expansion can be truncated aswhere Using the linear interpolation, (3) can be expressed aswhere Similarly, the Taylor expansion of the displacement about can be written aswhere Employing the linear interpolation to the above equation, we havewhere . When the governing equation (1) is expressed at time step , and , the following equations can be obtained:Similar procedure to (4) and (6) is applied to the value of :There are three nonderivative terms and six derivative terms (in seven equations (see (4) and (6)–(11)). After simple mathematical manipulation, we can finally obtain the update equation without the derivative terms:

Now, let us derive the Newmark-FDTD method for the modified Lorentz dispersion model. The modified Lorentz model [15, 16] can be expressed as , whereand 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]: Ampere’s law can be written as where , is the permittivity in the free space, and is the constitutive relation in the frequency domain. Inverse Fourier transform (IFT) is applied to the constitutive relation: where a temporary variable** W** is introduced to apply the Newmark time-stepping algorithm.

Therefore, we have two sets of differential equations: By applying the Newmark time-stepping algorithm to the above equations, the solution in terms of** P** and** E **can be obtained as follows:where Therefore, the** E** field can be updated by inserting of (19) into the CDS version of Ampere’s law (15):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):This equation can also be written in the same form as (19), where the coefficients are expressed as Now, let us consider the ADE-FDTD formulation based on the BT. The BT is basically an approximate of [20–22]. When the BT is applied to the constitutive relation, (16) can be expressed asThe above equation can also be expressed as the form of (19), where the coefficients arePlease note that the above coefficients are exactly the same as (20) with and .

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 two* Z*-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 By transforming the* Z*-domain into the* r*-domain, roots inside the unit circle in the* Z*-domain map onto the left half-plane in the* r*-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:where , , and 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

Next, the numerical permittivity can be derived inserting the numerical solution and into the constitutive relation [24, 28]. The numerical permittivity of the ADE-FDTD formulation can be obtained as follows:where , and The numerical permittivity of the ADE-FDTD formulation based on the BT is the same as that of the Newmark-FDTD formulation:

We also consider Debye, Drude, Lorentz, and QCRF dispersion models. The dispersion models are as follows:For Debye, Drude, and Lorentz dispersion models, can be written as , where* M*,* C*, and are constants for each model. By applying the Newmark-FDTD, one can obtainwhere 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 density** D** and electric field** E** 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 asAll things considered, the Newmark-FDTD formulation with and 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 , and The space step size and time step size , where is a classic Courant number [1] and 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 the* numerical* coefficients and are involved in the ADE-FDTD method in addition to the* numerical* angular frequency , differently from both the Newmark-FDTD and the ADE-FDTD with BT (see (28)-(29)).

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 the* Z*-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 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 the* Z*-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.

**(a)**

**(b)**

**(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, , and 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 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.

**(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.

#### Data Availability

The data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This research was supported in part by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (no. 2017R1D1A1B03034537) and in part by Ministry of Culture, Sports and Tourism (MCST) and Korea Creative Content Agency (KOCCA) in the Convergence Tourism Service Research and Business Development Program (no. SF0718106).