Research Article | Open Access

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

**Academic Editor:**Rodolfo Araneo

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

#### References

- A. Taflove and S. C. Hagness,
*Computational Electrodynamics: The Finite-Difference Time-Domian Method*, Artech House, Norwood, Mass, USA, 3rd edition, 2000. View at: MathSciNet - S.-G. Ha, J. Cho, J. Lee, B.-W. Min, J. Choi, and K.-Y. Jung, “Numerical study of estimating the arrival time of UHF signals for partial discharge localization in a power transformer,”
*Journal of Electromagnetic Engineering and Science*, vol. 18, no. 2, pp. 94–100, 2018. View at: Publisher Site | Google Scholar - J. Baek, D. Kim, and K. Jung, “Finite-difference time-domain modeling for electromagnetic wave analysis of human voxel model at millimeter-wave frequencies,”
*IEEE Access*, vol. 7, pp. 3635–3643, 2019. View at: Publisher Site | Google Scholar - J. A. Pereda, L. A. Vielva, A. Vegas, and A. Prieto, “Analyzing the stability of the FDTD technique by combining the von Neumann method with the Routh-Hurwitz criterion,”
*IEEE Transactions on Microwave Theory & Techniques*, vol. 49, no. 2, pp. 377–381, 2001. View at: Publisher Site | Google Scholar - K.-Y. Jung, F. L. Teixeira, and R. M. Reano, “Au/SiO2 nanoring plasmon waveguides at optical communication band,”
*Journal of Lightwave Technology*, vol. 25, no. 9, pp. 2757–2765, 2007. View at: Publisher Site | Google Scholar - J.-H. Kweon, M.-S. Park, J. Cho, and K.-Y. Jung, “FDTD analysis of electromagnetic wave propagation in an inhomogeneous ionosphere under arbitrary-direction geomagnetic field,”
*Journal of Electromagnetic Engineering and Science*, vol. 18, no. 3, pp. 212–214, 2018. View at: Publisher Site | Google Scholar - J. L. Young, “Propagation in linear dispersive media: finite difference time-domain methodologies,”
*IEEE Transactions on Antennas and Propagation*, vol. 43, no. 4, pp. 422–426, 1995. View at: Publisher Site | Google Scholar - S.-G. Ha, J. Cho, J. Choi, H. Kim, and K.-Y. Jung, “FDTD dispersive modeling of human tissues based on quadratic complex rational function,”
*IEEE Transactions on Antennas and Propagation*, vol. 61, no. 2, pp. 996–999, 2013. View at: Publisher Site | Google Scholar | MathSciNet - J. Cho, S. Ha, Y. B. Park, H. Kim, and K. Jung, “On the numerical stability of finite-difference time-domain for wave propagation in dispersive media using quadratic complex rational function,”
*Electromagnetics*, vol. 34, no. 8, pp. 625–632, 2014. View at: Publisher Site | Google Scholar - E.-K. Kim, S.-G. Ha, J. Lee, Y. B. Park, and K.-Y. Jung, “Three-dimensional efficient dispersive alternating-direction-implicit finite-difference time-domain algorithm using a quadratic complex rational function,”
*Optics Express*, vol. 23, no. 2, pp. 873–881, 2015. View at: Publisher Site | Google Scholar - S.-M. Park, E.-K. Kim, Y. B. Park, S. Ju, and K.-Y. Jung, “Parallel dispersive FDTD method based on the quadratic complex rational function,”
*IEEE Antennas and Wireless Propagation Letters*, vol. 15, pp. 425–428, 2016. View at: Publisher Site | Google Scholar - M. Han, R. W. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,”
*IEEE Microwave and Wireless Components Letters*, vol. 16, no. 3, pp. 119–121, 2006. View at: Publisher Site | Google Scholar - K. A. Michalski, “On the low-order partial-fraction fitting of dielectric functions at optical wavelengths,”
*IEEE Transactions on Antennas and Propagation*, vol. 61, no. 12, pp. 6128–6135, 2013. View at: Publisher Site | Google Scholar - A. Vial, T. Laroche, M. Dridi, and L. Le Cunff, “A new model of dispersion for metals leading to a more accurate modeling of plasmonic structures using the FDTD method,”
*Applied Physics A: Materials Science & Processing*, vol. 103, no. 3, pp. 849–853, 2011. View at: Publisher Site | Google Scholar - A. Deinega and S. John, “Effective optical response of silicon to sunlight in the finite-difference time-domain method,”
*Optics Expresss*, vol. 37, no. 1, pp. 112–114, 2012. View at: Publisher Site | Google Scholar - K. P. Prokopidis and D. C. Zografopoulos, “A unified FDTD/PML scheme based on critical points for accurate studies of plasmonic structures,”
*Journal of Lightwave Technology*, vol. 31, no. 15, pp. 2467–2476, 2013. View at: Publisher Site | Google Scholar - B. Wei, L. Cao, F. Wang, and Q. Yang, “Frequency-dependent FDTD algorithm using Newmark’s method,”
*International Journal of Antennas and Propagation*, vol. 2014, Article ID 216763, 6 pages, 2014. View at: Publisher Site | Google Scholar - Y.-Q. Zhang and P.-J. Yang, “An extended Newmark-FDTD method for complex dispersive media,”
*International Journal of Antennas and Propagation*, vol. 2018, Article ID 1573512, 7 pages, 2018. View at: Publisher Site | Google Scholar - M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple treatment of multi-term dispersion in FDTD,”
*IEEE Microwave and Guided Wave Letters*, vol. 7, no. 5, pp. 121–123, 1997. View at: Publisher Site | Google Scholar - J. A. Pereda, Á. Vegas, and A. Prieto, “FDTD modeling of wave propagation in dispersive media by using the Mobius transformation technique,”
*IEEE Transactions on Microwave Theory and Techniques*, vol. 50, no. 7, pp. 1689–1695, 2002. View at: Publisher Site | Google Scholar - Z. Lin, Y. Fang, J. Hu, and C. Zhang, “On the FDTD formulation for modeling wideband Lorentzian media,”
*IEEE Transactions on Antennas and Propagation*, vol. 59, no. 4, pp. 1338–1346, 2011. View at: Publisher Site | Google Scholar | MathSciNet - K. P. Prokopidis and D. C. Zografopoulos, “Efficient FDTD algorithms for dispersive Drude-critical points media based on bilinear z-transform,”
*IEEE Electronics Letters*, vol. 49, no. 8, pp. 534–536, 2013. View at: Publisher Site | Google Scholar - O. C. Zienkiewicz, “A new look at the Newmark, Houbolt and other time stepping formulas: a weighted residual approach,”
*Earthquake Engineering & Structural Dynamics*, vol. 5, no. 4, pp. 413–418, 1977. View at: Publisher Site | Google Scholar - W. L. Wood, “A further look at Newmark, Houbolt, etc., time-stepping formulae,”
*International Journal for Numerical Methods in Engineering*, vol. 20, no. 6, pp. 1009–1017, 1984. View at: Publisher Site | Google Scholar | MathSciNet - A. H. Panaretos and R. E. Díaz, “On the stability of the finite-difference time-domain modeling of lorentz media,”
*IEEE Microwave and Wireless Components Letters*, vol. 21, no. 6, pp. 283–285, 2011. View at: Publisher Site | Google Scholar - K. P. Prokopidis and D. C. Zografopoulos, “Investigation of the stability of ADE-FDTD methods for modified lorentz media,”
*IEEE Microwave and Wireless Components Letters*, vol. 24, no. 10, pp. 659–661, 2014. View at: Publisher Site | Google Scholar - O. Ramadan, “Stability considerations for the direct FDTD implementation of the dispersive quadratic complex rational function models,”
*IEEE Transactions on Antennas and Propagation*, vol. 64, no. 11, pp. 4929–4932, 2016. View at: Publisher Site | Google Scholar - S.-G. Ha, J. Cho, E.-K. Kim, Y. B. Park, and K.-Y. Jung, “FDTD dispersive modeling with high-order rational constitutive parameters,”
*IEEE Transactions on Antennas and Propagation*, vol. 63, no. 9, pp. 4233–4238, 2015. View at: Publisher Site | Google Scholar - J. Robinson and Y. Rahmat-Samii, “Particle swarm optimization in electromagnetics,”
*IEEE Transactions on Antennas and Propagation*, vol. 52, no. 2, pp. 397–407, 2004. View at: Publisher Site | Google Scholar | MathSciNet - H. Chung, S. Ha, J. Choi, and K. Jung, “Accurate FDTD modelling for dispersive media using rational function and particle swarm optimisation,”
*International Journal of Electronics*, vol. 102, no. 7, pp. 1218–1228, 2015. View at: Publisher Site | Google Scholar - D. Andreuccetti, R. Fossi, and C. Petrucci,
*Dielectric Properties of Body Tissues*, http://niremf.ifac.cnr.it/tisspro. - Z. Han, E. Forsberg, and S. He, “Surface plasmon bragg gratings formed in metal-insulator-metal waveguides,”
*IEEE Photonics Technology Letters*, vol. 19, no. 2, pp. 91–93, 2007. View at: Publisher Site | Google Scholar - J. A. Roden and S. D. Gedney, “Convolution PML (CPML): an efficient FDTD implementation of the CFS-PML for arbitrary media,”
*Microwave and Optical Technology Letters*, vol. 27, no. 5, pp. 334–339, 2000. View at: Publisher Site | Google Scholar - K.-Y. Jung, B. Donderici, and F. L. Teixeira, “Transient analysis of spectrally asymmetric magnetic photonic crystals with ferromagnetic losses,”
*Physical Review B*, vol. 74, Article ID 165207, pp. 1–11, 2006. View at: Google Scholar - M. Born and E. Wolf,
*Principles of Optics*, Macmillan Co., NewYork, NY, USA, 2nd edition, 1964.

#### Copyright

Copyright © 2019 Hongjin Choi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.