Numerical Stability of Modified Lorentz FDTD Unified From Various Dispersion Models

Jaesun Park; Kyung-Young Jung

doi:10.1364/OE.428656

1. Introduction

The finite-difference time-domain (FDTD) method [1] which can yield a broadband response in one simulation has been popularly used to study the various electromagnetic wave problems due to its accuracy, robustness and simplicity [2–5]. The FDTD method has been widely employed for the electromagnetic analysis of complex media [6–8]. When investigating the interaction of electromagnetic waves and dispersive media, a variety of dispersion models are used, including Debye [9,10], Drude [11,12], Lorentz [9,13] models. Besides these classical dispersion models, new dispersion models has been also widely used: complex-conjugate pole-residue (CCPR) [14], quadratic complex rational function (QCRF) [15], and modified Lorentz (mLor) [16,17]. The CCPR model was used for the material dispersion modeling of silver [14], fused silica [18], human tissues [19], fluoride glass [19], cold plasma [19], gold [19], graphene [20,21], and P3HT:PC61BM [22]. Note that the QCRF and mLor dispersion models has more degree of freedoms than the conventional Debye, Drude, Lorentz dispersion models and thus they can be widely employed for a variety of electromagnetic wave propagation and scattering problems in complex media. For example, the QCRF dispersion model was utilized for the electromagnetic analysis of human skin and organs [23–25], concrete materials with various water/cement ratios [24,26,27], amorphous and crystalline silicon [28], silver [28,29], cadmium telluride based alloys [30], and an aquifer [31]. The mLor dispersion model was employed to study electromagnetic wave propagation in silver [32–34], gold [32–34], silicon [33], cold plasma [35], human tissues [36], and graphene [37]. Among them, the mLor dispersion model is promising for FDTD simulations of complex dispersive media, because it can unify the other dispersion models mentioned above [32,38,39]. Therefore, based on the mLor dispersion model, a one unified FDTD code for complex dispersive media can be simply developed without including various dispersive FDTD formulations for other dispersion models. In this unified dispersive FDTD framework, the electromagnetic wave interaction with the complex dispersive media can be straightforwardly analyzed by utilizing various existing dispersion model parameters. Note that the mLor FDTD formulation can also produce good performance in computational efficiency [36,37].

In the FDTD method, the time step size should be limited to ensure the numerical stability. In other words, the FDTD method is conditionally stable due to the explicit solution of the time-dependent Maxwell curl (partial differential) equations. In the nondispersive FDTD methods, only one condition is considered, i.e., the Courant number is equal or less than one. In the dispersive FDTD methods, however, additional conditions should be determined to guarantee the numerical stability. Numerical stability studies were successfully conducted for the dispersive FDTD formulations based on Debye [9], Lorentz [9], CCPR [22], QCRF [24], and mLor [17,39]. Various dispersive FDTD formulations include recursive convolution, piecewise linear recursive convolution, $Z$ transform and auxiliary differential equation (ADE) methods [40–42]. Note the ADE method involves simpler arithmetic operations than other methods and it can be also easily extended to other complex media, as pointed out in [32,39]. In the ADE FDTD method, the double averaging scheme to 0th-order derivatives (equivalent to the bilinear transformation) can yield better performance in the aspects of numerical stability and numerical accuracy, compared to the direct scheme [9,23,24,39]. The ADE method with the double averaging scheme is preferable and thus in this work we discuss dispersive FDTD formulations based the ADE technique with the double averaging scheme.

To fully utilize the mLor FDTD formulation unified from the other dispersion models, it is of necessity to study its numerical stability in terms of the original (Debye, Drude, Lorentz, CCPR, or QCRF) model parameters. In this work, we investigate comprehensively the numerical stability conditions of the unified mLor FDTD formulation in aspects of the other dispersion model parameters. The remainder of this paper is organized as follows. The dispersive FDTD formulation (based on mLor, Debye, Drude, Lorentz, CCPR, or QCRF) and its numerical stability conditions are revisited. And then, the numerical stability study of the unified mLor FDTD formulation is performed by using the model parameter equivalence. It is elucidated that the numerical stability conditions of the unified mLor FDTD formulation are equivalent to those of the original model-based dispersive FDTD formulation. Albeit with this, for the QCRF model, a proper Courant number that satisfies the numerical stability of its unified mLor FDTD formulation should be chosen, different from the Debye, Drude, Lorentz, and CCPR models. Otherwise, its unified mLor FDTD simulation may be numerically unstable. Various numerical examples are employed to validate our investigations. Finally, concluding remarks are provided.

2. Numerical stability conditions of dispersive FDTD formulations

An $e^{j\omega t}$ time dependence is assumed. First of all, we briefly review the mLor FDTD formulation and its numerical stability conditions.

The mLor model [16] can be written as $\varepsilon _r(\omega )=\varepsilon _{r,\infty }+\chi (\omega )$, where

(1)$$\chi(\omega)=\frac{a_0+a_1(j\omega)}{b_0+b_1(j\omega)+b_2(j\omega)^2},$$

$\varepsilon _{r,\infty }$ is the relative permittivity at the infinite frequency. It should be noted that $b_{2}=1$ in above leads to the generalized dispersion model (GDM) which has been successfully implemented in both FDTD and discontinuous Galerkin time domain (DGTD) solvers [43–45]. To obtain the mLor FDTD update equation, the constitutive relationship between the electric field $(\textbf {E})$ and the current density $(\textbf {J})$, i.e., $\textbf {J}(\omega )=j\omega \varepsilon _{0}\chi (\omega )\textbf {E}(\omega )$, is considered, where $\varepsilon _{0}$ denotes the permittivity in the free space. In what follows, we employ the central finite difference operator $\delta _{t}f^{n}\equiv f^{n+1/2}-f^{n-1/2}$ and the central average operator $\mu _{t}f^{n}\equiv 0.5(f^{n+1/2}+f^{n-1/2})$. By applying the inverse Fourier transform (IFT) and the central difference scheme (CDS), we can derive

(2)$$\left[b_{0}\mu_{t}^{2}+b_{1}\frac{\delta_{t}\mu_{t}}{\Delta t}+b_{2}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{J}^{n}=\varepsilon_{0}\left[a_{0}\frac{\delta_{t}\mu_{t}}{\Delta t}+a_{1}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{E}^{n},$$

where $\Delta t$ indicates the FDTD time step size. In addition to the constitutive relation, the Ampere’s law should be also considered:

(3)$$\nabla\times\textbf{H}(t)=\varepsilon_{\infty}\frac{\partial \textbf{E}(t)}{\partial t}+\textbf{J}(t).$$

By utilizing the CDS, one can have

(4)$$\varepsilon_{\infty}\frac{\delta_{t}}{\Delta t}\textbf{E}^{n+1/2}=\nabla\times\textbf{H}^{n+1/2}-\mu_{t}\textbf{J}^{n+1/2},$$

where $\varepsilon _{\infty }=\varepsilon _{0}\varepsilon _{r,\infty }$. Note that the magnetic field (H) update equation can be obtained using the Faraday’s law

(5)$$\nabla\times\textbf{E}(t)={-}\mu_{0}\frac{\partial \textbf{H}(t)}{\partial t},$$

where $\mu _{0}$ means the permeability in the free space. $\textbf {H}$ can be updated by applying the CDS to the Faraday’s law:

(6)$$-\mu_{0}\frac{\delta_{t}}{\Delta t}\textbf{H}^{n}=\nabla\times\textbf{E}^{n}.$$

To obtain the numerical stability conditions, the von Neumann method with the Routh-Hurwitz (R-H) criterion [9] is utilized. A time-dependent wave equation in a source-free homogeneous dispersive medium can be derived by equating Eqs. (3) and (5),

(7)$$c_{\infty}^{2}\varepsilon_{\infty}\nabla^{2}\textbf{E}-\frac{\partial\textbf{J}}{\partial t}-\varepsilon_{\infty}\frac{\partial^{2}\textbf{E}}{\partial t^{2}}=0,$$

where $c_{\infty }=1/\sqrt {\mu _{0}\varepsilon _{\infty }}$. By applying the CDS and the Z-transform to the wave equation, we can obtain

(8)$$\varepsilon_{\infty}(Z-1)^{2}\textbf{E}_{0}+4\varepsilon_{\infty}\nu^{2}Z\textbf{E}_{0}+\frac{Z^{2}-1}{2}\Delta t\textbf{J}_{0}=0,$$

where

(9)$$\nu^{2}=(c_{\infty}\Delta t)^{2}\sum_{\alpha=x,y,z}{\frac{\textrm{sin}^{2}(\widetilde{k_{\alpha}}\Delta_\alpha/2)}{(\Delta_{\alpha})^{2}}},$$

$\widetilde {k_{\alpha }}$ is the numerical wavenumber in the $\alpha$ direction, $\textbf {J}_{0}$ and $\textbf {E}_{0}$ are the complex amplitudes of field quantities, Z is an amplification factor, and $\Delta _{\alpha }$ indicates the FDTD space step size. By solving the simultaneous equation composed of Eq. (8) and the Z-transformed equation of Eq. (2), we can derive the following stability polynomial in the Z-domain:

(10)$$\begin{aligned}S(Z) = & Z^4\big{[}(b_0\varepsilon_{\infty}+a_0\varepsilon_{0})\Delta t^2 + 2(b_1\varepsilon_{\infty}+a_1\varepsilon_{0})\Delta t + 4b_2\varepsilon_{\infty}\big{]}\\ &+ Z^3\big{[}4b_0\varepsilon_{\infty}\nu^2\Delta t^2 -4(b_1\varepsilon_{\infty}+a_1\varepsilon_0)\Delta t + 8b_1\varepsilon_{\infty}\nu^2\Delta t -16b_2\varepsilon_{\infty}(1-\nu^2)\big{]}\\ &+ Z^2\big{[}-2(b_0\varepsilon_{\infty}+a_0\varepsilon_{0})\Delta t^2 + 8b_0\varepsilon_{\infty}\nu^2\Delta t^2+24b_2\varepsilon_{\infty}-32b_2\varepsilon_{\infty}\nu^2\big{]}\\ &+ Z\big{[}4b_0\varepsilon_{\infty}\nu^2\Delta t^2 + 4(b_1\varepsilon_{\infty}+a_1\varepsilon_{0})\Delta t-8b_1\varepsilon_{\infty}\nu^2\Delta t-16b_2\varepsilon_{\infty}(1-\nu^2)\big{]}\\ &+ \bigg{[}(b_0\varepsilon_{\infty}+a_0\varepsilon_{0})\Delta t^2-2(b_1\varepsilon_{\infty}+a_1\varepsilon_{0})\Delta t+4b_2\varepsilon_{\infty}\bigg{]}=0.\end{aligned}$$

Based on $Z=(\mathit {r}+1)/(\mathit {r}-1)$, the stability polynomial in the $r$-plane can be written as

(11)$$\begin{aligned} S(r) =&r^4\big{[}b_0\varepsilon_{\infty}\nu^2\Delta t^2\bigg{]}+r^3\bigg{[}2b_1\varepsilon_{\infty}\nu^2\Delta t\bigg{]}\\ = & r^2\bigg{[}(b_0\varepsilon_{\infty}+a_0\varepsilon_{0})\Delta t^2-b_0\varepsilon_{\infty}\nu^2\Delta t^2+4b_2\varepsilon_{\infty}\nu^2\bigg{]}\\ &+r\bigg{[}2(b_1\varepsilon_{\infty}+a_1\varepsilon_{0})\Delta t-2b_1\varepsilon_{\infty}\nu^2\Delta t\bigg{]}+\bigg{[}4b_2\varepsilon_{\infty}(1-\nu^2)\bigg{]}. \end{aligned}$$

In order to be stable, all roots should be located on the left half or imaginary axis on the $r$-plane. Using the R-H criterion, the numerical stability conditions of the mLor FDTD formulation can be finally obtained [39], as shown in Table 1. Note that $\varepsilon _{r,\infty }$ is positive, although it is not explicitly shown in Table 1. In above, the numerical stability conditions for the mLor FDTD formulation are determined based on the constitutive relation between $\textbf {E}$ and $\textbf {J}$. However, it is also available to develop the mLor FDTD formulation considering the constitutive relation between $\textbf {E}$ and the electric flux density ($\textbf {D}$). Note that these two different FDTD formulations lead to the same numerical stability conditions [39]. In what follows, the Courant number is defined as $C_{n}=c_{\infty } \Delta t / \Delta z$ [17,39] for a one-dimensional (1fD) wave propagation, unless specified otherwise.

Table 1. Numerical Stability Conditions of mLor FDTD

View Table | View all tables in this article

Next, we summarize the FDTD update equations for various dispersion models (which are unified into the mLor dispersion model) and the corresponding numerical stability conditions. And then, based on those numerical stability conditions, we examine the numerical stability of the unified mLor FDTD formulation in a systematic fashion.

2.1 Debye model

The relative permittivity of the Debye model [9] can be expressed as

(12)$$\varepsilon_{r}{(\omega)}{=}\varepsilon_{r,\infty}{+}\frac{\varepsilon_{s}-\varepsilon_{r,\infty}}{1+j\omega\tau},$$

where $\varepsilon _{s}$ is the static relative permittivity, and ${\tau }$ indicates the relaxation time constant. As mentioned previously, the mLor dispersion model can cover various dispersion models such as Debye, Drude, Lorentz, CCPR, and QCRF. The coefficients of each dispersion model can be equivalent to those of the mLor dispersion model. The equivalent coefficients of the mLor dispersion model for the Debye dispersion model are $a_{0}=\varepsilon _{s}-\varepsilon _{r,\infty }, a_{1}=0, b_{0}=1, b_{1}=\tau , b_{2}=0$. According to [9], the constitutive relationship between $\textbf {E}$ and $\textbf {D}$, $\textbf {D}(\omega )=\varepsilon _{0}\varepsilon _{r}(\omega )\textbf {E}(\omega )$, is considered. By applying the IFT and the CDS, we can obtain

(13)$$\left[\mu_{t}+\tau\frac{\delta_{t}}{\Delta t}\right]\textbf{D}^{n}=\varepsilon_{0}\left[\varepsilon_{s}\mu_{t}+\tau\varepsilon_{r,\infty}\frac{\delta_{t}}{\Delta t}\right]\textbf{E}^{n}.$$

The Ampere’s law is written as

(14)$$\nabla\times\textbf{H}(t)=\frac{\partial\textbf{D}(t)}{\partial t}.$$

D can be updated by applying the CDS to the above Ampere’s law:

(15)$$\frac{\delta t}{\Delta t}\textbf{D}^{n+1/2}=\nabla\times\textbf{H}^{n+1/2}.$$

In this case, a wave equation in a source-free homogeneous dispersive medium can be derived by using Eqs. (5) and (14)

(16)$$\frac{\partial^{2}\textbf{D}}{\partial t^{2}}-\varepsilon_{\infty}c_{\infty}^{2}\nabla^{2}\textbf{E}=0.$$

By applying the CDS and the Z-transform to the wave equation, we can derive

(17)$$(Z-1)^{2}\textbf{D}_{0}+4\varepsilon_{\infty}\nu^{2}Z\textbf{E}_{0}=0,$$

where $\textbf {E}_{0}$ and $\textbf {D}_{0}$ are the complex amplitudes of field quantities. By considering Eq. (17) and the Z-transformed equation of Eq. (13) simultaneously, we can derive the stability polynomial in the Z-domain:

(18)$$\begin{aligned} S(Z)= & Z^3\left[\frac{\varepsilon_{s}}{\varepsilon_{r,\infty}}+2\frac{\tau}{\Delta t}\right]+Z^2\left[4\nu^2\left(1+2\frac{\tau}{\Delta t}\right)-\frac{\varepsilon_{s}}{\varepsilon_{r,\infty}}-6\frac{\tau}{\Delta t}\right]\\ & + Z\left[4\nu^2\left(1-2\frac{\tau}{\Delta t}\right)-\frac{\varepsilon_{s}}{\varepsilon_{r,\infty}}+6\frac{\tau}{\Delta t}\right]+\left[\frac{\varepsilon_{s}}{\varepsilon_{r,\infty}}-2\frac{\tau}{\Delta t}\right]=0. \end{aligned}$$

Finally, using the von Neumann method with the R-H criterion, the numerical stability conditions of the Debye FDTD formulation are obtained [9] and they are summarized on the left column of Table 2.

Table 2. Numerical Stability Conditions of Debye FDTD and its unified mLor FDTD

View Table | View all tables in this article

Now, let us investigate the numerical stability of the mLor FDTD unified from the Debye model. Toward this purpose, the equivalent mLor parameters for the Debye model are plugged into Table 1 and the resulting numerical stability conditions are shown in the right column of Table 2. A close look at the right column of Table 2 leads to the following simple numerical stability conditions: $\tau \geq 0$, $\varepsilon _{s}\geq \varepsilon _{r,\infty }$, $\nu ^{2}\leq 1$ (as mentioned previously, $\varepsilon _{r,\infty } > 0$). Therefore, the numerical stability conditions of the mLor FDTD formulation unified from the Debye model are equivalent to those of the Debye FDTD formulation. Therefore, one can simply implement the unified mLor FDTD simulation by using the parameter equivalence as long as the numerically stable Debye parameters are adopted.

2.2 Drude model

The relative permittivity of the Drude dispersion model [11] can be written as

(19)$$\varepsilon_{r}{(\omega)}{=}\varepsilon_{r,\infty}{+}\frac{\omega_{0}^{2}}{\gamma(j\omega)+(j\omega)^{2}},$$

where $\omega _{0}$ is the resonant frequency and ${\gamma }$ indicates the inverse of the relaxation time. Equivalent mLor parameters for the Drude model give $a_{0}=\omega _{0}^{2}, a_{1}=0, b_{0}=0, b_{1}=\gamma , b_{2}=1$. By applying the IFT and the CDS to the constitutive relationship $\textbf {D}(\omega )=\varepsilon _{0}\varepsilon _{r}(\omega )\textbf {E}(\omega )$ with the Drude dispersion, we can have

(20)$$\left[\gamma\frac{\delta_{t}\mu_{t}}{\Delta t}+\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{D}^{n}=\varepsilon_{0}\left[\omega_{0}^{2}\mu_{t}^{2}+\varepsilon_{r,\infty}\gamma\frac{\delta_{t}\mu_{t}}{\Delta t}+\varepsilon_{r,\infty}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{E}^{n}.$$

Following the similar procedure in the previous subsection, the stability polynomial for the Drude FDTD formulation is as follows

(21)$$\begin{aligned} S(Z) = & Z^4\left[\varepsilon_{0}\omega_{0}^2\Delta t^2+2\varepsilon_{\infty}\gamma\Delta t+4\varepsilon_{\infty}\right]+Z^3\left[-4\varepsilon_{\infty}\gamma(1-2\nu^2)\Delta t-16\varepsilon_{\infty}(1-\nu^2)\right]\\ & + Z^2\left[-2\varepsilon_{0}\omega_{0}^2\Delta t^2+24\varepsilon_{\infty}-32\varepsilon_{\infty}\nu^2\right]+Z\left[4\varepsilon_{\infty}\gamma(1-2\nu^2)\Delta t-16\varepsilon_{\infty}(1-\nu^2)\right]\\ & + \left[\varepsilon_{0}\omega_{0}^2\Delta t^2-2\varepsilon_{\infty}\gamma\Delta t+4\varepsilon_{\infty}\right]=0. \end{aligned}$$

In addition, the numerical stability conditions of the Drude FDTD formulation are derived and they are shown in the first column of Table 3.

Table 3. Numerical Stability Conditions of Drude FDTD and its Unified mLor FDTD

View Table | View all tables in this article

In order to examine the numerical stability of the mLor FDTD formulation unified from the Drude model, its mLor parameter equivalence is applied to Table 1. The second column of Table 3 lists the numerical stability conditions of the unified mLor FDTD formulation for the Drude model. The numerical stability conditions in the second column can be concisely expressed as $\gamma \geq 0$, $\nu ^{2}\leq 1$. Again, equivalence is observed in terms of the numerical stability conditions between the Drude FDTD formulation and its unified mLor FDTD formulation.

2.3 Lorentz model

The relative permittivity of the Lorentz dispersion model [13] is expressed as

(22)$$\varepsilon_{r}{(\omega)}{=}\varepsilon_{r,\infty}{+}\frac{(\varepsilon_{s}-\varepsilon_{r,\infty})\omega_{0}^{2}}{\omega_{0}^{2}+2\delta(j\omega)+(j\omega)^{2}},$$

where $\delta$ is the damping coefficient. Equivalent mLor parameters for the Lorentz model are $a_{0}=(\varepsilon _{s}-\varepsilon _{r,\infty })\omega _{0}^{2}, a_{1}=0, b_{0}=\omega _{0}^{2}, b_{1}=2\delta , b_{2}=1$. The constitutive equation for the Lorentz model in the discrete FDTD time domain is written as

(23)$$\left[\omega_{0}^{2}\mu_{t}^{2}+2\delta\frac{\delta_{t}\mu_{t}}{\Delta t}+\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{D}^{n}=\varepsilon_{0}\left[\varepsilon_{s}\omega_{0}^{2}\mu_{t}^{2}+2\delta\varepsilon_{r,\infty}\frac{\delta_{t}\mu_{t}}{\Delta t}+\varepsilon_{r,\infty}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{E}^{n}.$$

The stability polynomial for the Lorentz FDTD formulation is

(24)$$ \begin{aligned} S(Z) &= Z^4\left[\varepsilon_{s}\omega_{0}^2\Delta t^2+4\varepsilon_{r,\infty}\delta\Delta t+4\varepsilon_{r,\infty}\right]\\ & + Z^3\left[4\omega_{0}^2\nu^2\Delta t^2-8\delta(\varepsilon_{r,\infty}-2\nu^2)\Delta t-16(\varepsilon_{r,\infty}-\nu^2)\right]\\ &+ Z^2\left[-2\omega_0^2(\varepsilon_{s}-4\nu^2)\Delta t^2+24\varepsilon_{r,\infty}-32\nu^2\right]\\ & + Z\left[4\omega_0^2\nu^2\Delta t^2+8\delta(\varepsilon_{r,\infty}-2\nu^2)\Delta t-16(\varepsilon_{r,\infty}-\nu^2)\right]\\ & + \left[\varepsilon_{s}\omega_{0}^2\Delta t^2-4\varepsilon_{r,\infty}\delta\Delta t+4\varepsilon_{r,\infty}\right]=0. \end{aligned}$$

The numerical stability conditions of the Lorentz FDTD formulation [9] are listed in the first column of Table 4. The second column of Table 4 shows the numerical stability conditions of the unified mLor FDTD formulation for the Lorentz model and they can be equivalently expressed as $\delta \geq 0$, $\nu ^{2}\leq 1$, $\varepsilon _{s}\geq \varepsilon _{r,\infty }$. It is found out that the numerical stability conditions for both FDTD formulations in Table 4 are equivalent to each other.

Table 4. Numerical Stability Conditions of Lorentz FDTD and its Unified mLor FDTD

View Table | View all tables in this article

2.4 CCPR model

The relative permittivity of the CCPR dispersion model [14] can be expressed as $\varepsilon _{r}(\omega )=\varepsilon _{r,\infty }+\chi (\omega )$, where

(25)$$\chi(\omega)=\frac{r}{j\omega-p}{+}\frac{r^{*}}{j\omega-p^{*}},$$

$p$ and $r$ denote the pole and residue of the CCPR model, and $*$ means the complex conjugate. Equivalent mLor coefficients for the CCPR model are $a_{0}=-2{\rm Re}[rp^{*}], a_{1}=2{\rm Re}[r], b_{0}=\left | p \right |^{2}, b_{1}=-2{\rm Re}[p], b_{2}=1$. Using the IFT and the CDS, we can have the FDTD update equation for the CCPR model:

(26)$$\left[\left|p\right|^{2}\mu_{t}^{2}-2{\rm Re}(p)\frac{\delta_{t}\mu_{t}}{\Delta t}+\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{J}^{n}=\varepsilon_{0}\left[{-}2{\rm Re}(rp^{*})\frac{\delta_{t}\mu_{t}}{\Delta t}+2{\rm Re}(r)\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{E}^{n}.$$

According to [22], the stability polynomial can be written as

(27)$$\begin{aligned} S(Z) = & \big{[}Z^4+1\big{]}\big{[}(|p|^2\varepsilon_\infty-2\text{Re}(rp^*)\varepsilon_0)\Delta_t^2+4\varepsilon_\infty\big{]}-\big{[}Z^4-1\big{]}\\ &\times \big{[}4(\text{Re}(p)\varepsilon_\infty+\text{Re}(r)\varepsilon_0)\Delta_t\big{]}+4\big{[}Z^3+Z\big{]}\\ &\times \big{[}|p|^2\varepsilon_\infty\nu^2\Delta_t^2- 4\varepsilon_\infty(1-\nu^2)\big{]}+8\big{[}Z^3-Z\big{]}\big{[}(\text{Re}(p)\varepsilon_\infty-\text{Re}(r)\varepsilon_0)\Delta_t - 2\text{Re}(p)\varepsilon_\infty\nu^2\Delta_t\big{]}\\ & + 2Z^2\big{[}(2\text{Re}(rp^*)\varepsilon_0-|p|^2\varepsilon_\infty)\Delta_t^2 + 4|p|^2\varepsilon_\infty\nu^2\Delta_t^2+12\varepsilon_\infty-16\varepsilon_\infty\nu^2\big{]}=0. \end{aligned}$$

The numerical stability conditions of the CCPR FDTD formulation can be derived similarly and they are listed in the left column of Table 5. In addition, the right column shows the numerical stability conditions of the mLor FDTD formulation unified from the CCPR model. Note that $Q$ in the second column is twice the first-column $Q$. It is demonstrated that the numerical stability conditions of the CCPR FDTD formulation are exactly same as the unified mLor FDTD counterparts.

Table 5. Numerical Stability Conditions of CCPR FDTD and its Unified mLor FDTD

View Table | View all tables in this article

2.5 QCRF model

The relative permittivity of the QCRF model [15] can be expressed as

(28)$$\varepsilon_{r}{(\omega)}{=}\frac{A_{0}+A_{1}(j\omega)+A_{2}(j\omega)^{2}}{B_{0}+B_{1}(j\omega)+B_{2}(j\omega)^{2}}.$$

In above, positive QCRF parameters are assumed. Different from the Debye, Drude, Lorentz, and CCPR dispersion models, $\varepsilon _{r,\infty }$ is not explicitly expressed. Here, it can be represented by $\varepsilon _{r,\infty }={A_2}/{B_2}$. Equivalent mLor parameters for the QCRF model can lead to $a_{0}=A_{0}-A_{2}B_{0}/B_{2}, a_{1}=A_{1}-A_{2}B_{1}/B_{2}, b_{0}={B_{0}}, b_{1}={B_{1}}, b_{2}={B_{2}}$. The FDTD update equation for the QCRF model can be expressed as

(29)$$\left[B_{0}\mu_{t}^{2}+B_{1}\frac{\delta_{t}\mu_{t}}{\Delta t}+B_{2}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{D}^{n}=\varepsilon_{0}\left[A_{0}\mu_{t}^{2}+A_{1}\frac{\delta_{t}\mu_{t}}{\Delta t}+A_{2}\frac{\delta_{t}^{2}}{\Delta t^{2}}\right]\textbf{E}^{n}.$$

According to [24], a wave equation in a source-free homogeneous QCRF medium is expressed as

(30)$$\frac{\partial^{2}\textbf{D}}{\partial t^{2}}-\varepsilon_{0}c_{0}^{2}\nabla^{2}\textbf{E}=0,$$

where $c_{0}$ is the speed of light in free space, $c_{0}=1/\sqrt {\mu _{0}\varepsilon _{0}}$. By applying the CDS and the Z-transform to the wave equation, we have [24]

(31)$$(Z-1)^{2}\textbf{D}_{0}+4\varepsilon_{0}\nu_{q}^{2}Z\textbf{E}_{0}=0$$

where

(32)$$\nu_{q}^{2}=(c_{0}\Delta t)^{2}\sum_{\alpha=x,y,z}{\frac{\textrm{sin}^{2}(\widetilde{k_{\alpha}}\Delta_\alpha/2)}{(\Delta _{\alpha})^{2}}}.$$

By considering Eq. (31) and the Z-transformed equation of Eq. (29) simultaneously, we have the stability polynomial:

(33)$$ \begin{aligned} S(Z) =& Z^4\left[A_0\Delta t^2+2A_1\Delta t+4A_2\right]\\ &+ Z^3\left[4B_0\nu_q^2\Delta t^2-4(A_1-2B_1\nu_q^2)\Delta t-16(A_2-B_2\nu_q^2)\right]\\ &+ Z^2\left[-2(A_0-4B_0\nu_q^2)\Delta t^2+24A_2-32B_2\nu_q^2\right]\\ &+ Z\left[4B_0\nu_q^2\Delta t^2+4(A_1-2B_1\nu_q^2)\Delta t-16(A_2-B_2\nu_q^2)\right]\\ &+\left[A_0\Delta t^2-2A_1\Delta t+4A_2\right]=0. \end{aligned}$$

Applying the von Neumann method with the R-H criterion to the above equation leads to the numerical stability conditions of the QCRF FDTD formulation [24], which are listed in the first column of Table 6. Note that the parameter $\nu _{q}$ in the QCRF FDTD formulation is different from the previous dispersive FDTD cases. For the 1D QCRF FDTD formulation, the Courant number is defined as $C_{n}=c_{0} \Delta t / \Delta z$ [24].

Table 6. Numerical Stability Conditions of QCRF FDTD and its Unified mLor FDTD

View Table | View all tables in this article

Let us investigate the numerical stability of the mLor FDTD unified from the QCRF model. Following the similar procedure, the equivalent mLor parameters for the QCRF model are plugged into Table 1 and the resulting numerical stability conditions are shown in the second column of Table 6. At the first glance, both numerical stability conditions seem to be different to each other. However, the numerical stability conditions of the QCRF FDTD are exactly same as those of its unified mLor FDTD. To show this, let us derive the relation between $\nu$ of the mLor FDTD and $\nu _q$ of the QCRF FDTD. By comparing Eq. (9) by Eq. (32), it can be shown that

(34)$$\nu^2=\frac{(c_0 \Delta t)^2}{\varepsilon_{r,\infty}} \sum_{\alpha=x,y,z}{\frac{\textrm{sin}^{2}(\widetilde{k_{\alpha}}\Delta_\alpha/2)}{(\Delta_{\alpha})^{2}}}=\frac{1}{\varepsilon_{r,\infty}}\nu_q^2.$$

Note that $\varepsilon _{r,\infty }$ of mLor model can be represented as ${A_2}/{B_2}$ of the QCRF parameters and thus we rewrite

(35)$$ \nu^2=\frac{B_2}{A_2}\nu_q^2, $$

(36)$$ \nu_q^2=\frac{A_2}{B_2}\nu^2. $$

Plugging Eq. (36) into the first column of Table 6, it is observed that both the QCRF FDTD and its unifed mLor FDTD have the same numerical stability conditions. In the mLor FDTD formulation unified from the Debye, Drude, Lorentz, or CCPR models, the same Courant number can be used as in the original stable dispersive FDTD formulation. However, in the mLor FDTD formulation unified from the QCRF model, a proper Courant number should be chosen to ensure its numerical stability, which may be different from the Courant number in the stable QCRF FDTD formulation. This will be illustrated in detail in the next section.

3. Numerical examples

In this section, numerical examples are provided to investigate the numerical stability of the original dispersive FDTD formulations and their unified mLor FDTD formulation.

First, let us consider the Debye dispersion model. For liquid water [10], the Debye parameters are $\varepsilon _{r,\infty }=5.285$, $\varepsilon _{s}{-}\varepsilon _{r,\infty }=74.789$, and $\tau =9.352$ ps/rad. The FDTD time step size is indicated by $\Delta t=C_{n}\Delta z/c_{\infty }$, where $\Delta z=\lambda _{\rm min}/15$ in the frequency range of DC to 100 GHz. The root locus of the stability polynomial in the Z-domain is plotted in Fig. 1(a). In what follows, $C_{n}=0 \sim 1$ is considered for the root locus plot. It is observed that numerical stability conditions are satisfied for this Debye FDTD simulation (as referred in the left column in Table 2). As can be seen in the figure, it is numerically stable since all roots are located inside the unit circle. The root locus in the mLor FDTD formulation unified from this Debye model is also plotted in Fig. 1(b), indicating that the mLor FDTD formulation is also numerically stable. Therefore, the mLor FDTD algorithm unified from the Debye model can be simply implemented by using the mLor parameter equivalence.

Fig. 1. Stable root locus of the Debye model. (a) Debye. (b) Corresponding unified mLor.

Download Full Size | PDF

For the second example, we consider the Drude FDTD simulation with $\varepsilon _{r,\infty }=1$, $\omega _{0}=11.96\times {10}^{15}$ rad/s, $\gamma =80.52\times {10}^{12}$ rad/s, and $\Delta z=1.7$ nm. As can been from Table 3, the numerical stability conditions are satisfied for this Drude FDTD simulation and its unified mLor FDTD simulation. The root loci of the stability polynomial in the Z-domain are plotted in Fig. 2. It is observed that both simulations are numerically stable. Again, similar to the Debye case, the unified mLor FDTD algorithm for the Drude model is simply implemented by using the model parameter equivalence.

Fig. 2. Stable root locus of the Drude model. (a) Drude. (b) Corresponding unified mLor.

Download Full Size | PDF

Third, let us take a numerical example of the Lorentz dispersion model. The considered Lorentz parameters are $\varepsilon _{r,\infty }=1$, $\varepsilon _{s}=2.25$, $\omega _{0}=4\times {10}^{16}$ rad/s, $\delta =0.07\omega _{0}$ with $\Delta z=0.2$ nm [13]. Plugging these parameters into Table 4 shows that the numerical stability conditions are satisfied for the Lorentz and its unified mLor FDTD simulations. The root loci in Fig. 3 indicates that both FDTD simulations are numerically stable, as expected.

Fig. 3. Stable root locus of the Lorentz model. (a) Lorentz. (b) Corresponding unified mLor.

Download Full Size | PDF

Next, let us investigate the numerical stability of the CCPR FDTD formulation. For the electromagnetic wave analysis of silver [22], the CCPR-FDTD parameters are $\varepsilon _{r,\infty }=5.2830$, $p=-4.8113\cdot {10}^{12}+j6.8512\cdot {10}^{13}$, and $r=4.9581\cdot {10}^{14}-j1.5110\cdot {10}^{18}$ with $\Delta z=40$ nm, which are satisfied for the numerical stability conditions of the CCPR and its unified mLor FDTD formulation (Refer Table 5.). As shown in Fig. 4, all roots exist on or inside the unit circle. Note that the mLor FDTD algorithm unified from the CCPR model can be straightforwardly implemented by utilizing the mLor coefficient equivalence.

Fig. 4. Stable root locus of the CCPR model. (a) CCPR. (b) Corresponding unified mLor.

Download Full Size | PDF

Let us discuss the numerical stability conditions of the QCRF model in detail. Let us consider the QCRF FDTD simulation with $A_{0}=455.72$, $A_{1}=2.5\times {10}^{-7}$, $A_{2}=4.31\times {10}^{-18}$, $B_{0}=1$, $B_{1}=4.47\times {10}^{-9}$, $B_{2}=7.56\times {10}^{-20}$, and $\Delta z=1.38$ mm. The root locus for its FDTD simulation is plotted in Fig. 5(a) and it can be seen that it is numerically stable. Note that the numerical stability is also confirmed by the left column of Table 6. By using the mLor parameter equivalence to the QCRF model, we have the unified mLor parameters: $a_{0}=398.7094$, $a_{1}=-4.8373\times {10}^{-9}$, $b_{0}=1$, $b_{1}=4.47\times {10}^{-9}$, $b_{2}=7.56\times {10}^{-20}$, and $\varepsilon _{r,\infty }=57.0106$. For the unified mLor FDTD formulation, some roots are located outside the unit circle as shown in Fig. 5(b) and thus it may be numerically unstable. This can be also inferred from the right column of Table 6. Next, let us consider the QCRF parameters with $A_{2}=7.56\times {10}^{-20}$ and $B_{2}=1.98\times {10}^{-18}$ (while other parameters are unchanged). Its unified mLor parameters are $a_{0}=455.6818$, $a_{1}=2.4983\times {10}^{-7}$, $b_{0}=1$, $b_{1}=4.47\times {10}^{-9}$, $b_{2}=1.98\times {10}^{-18}$, and $\varepsilon _{r,\infty }=0.0382$. Different from the previous case, as can be seen from Table 6 and Fig. 6, the numerical stability of the QCRF FDTD simulation is not satisfied but its unified mLor FDTD simulation is numerically stable.

Fig. 5. Root locus of the QCRF model. (a) QCRF (stable). (b) Corresponding unified mLor (unstable).

Download Full Size | PDF

Fig. 6. Root locus of the QCRF model. (a) QCRF (unstable). (b) Corresponding unified mLor (stable).

Download Full Size | PDF

To confirm the root locus results, we also perform actual 1D FDTD simulations with $C_{n}=1$. The whole computational domain of 140 cells is modelled as the dispersive medium and it is truncated by the ten layers of the perfectly matched layer (PML) [1]. The Gaussian-modulated sinewave with the center frequency of 1.7 GHz and the half-power fractional bandwidth of 65 % is excited and the observation point is 50 cells away from the source point. Figure 7 shows the FDTD simulations corresponding to Fig. 5. As shown in the figure, the QCRF FDTD simulation is stable but its unified mLor FDTD simulation diverges. On the contrary, the QCRF FDTD simulation is unstable but its unified mLor FDTD simulation converges, as shown in Fig. 8 corresponding to Fig. 6. For the same $C_n$, the unified mLor FDTD simulation can be numerically unstable even though the QCRF FDTD simulation is numerically stable, vice versa.

Fig. 7. 1D FDTD simulation corresponding to Fig. 5. (a) QCRF. (b) Unified mLor.

Download Full Size | PDF

Fig. 8. 1D FDTD simulation corresponding to Fig. 6. (a) QCRF. (b) Unified mLor.

Download Full Size | PDF

As alluded previously, the numerical stability conditions are same to each other for the QCRF FDTD and its unified mLor FDTD. This indicates that (un)stable QCRF FDTD simulations lead to (un)stable mLor FDTD simulations. To validate this observation, we use the relationship Eq. (35) between $\nu$ and $\nu _q$ to find an appropriate Courant number that satisfies the numerical stability conditions of the mLor FDTD. We consider the same problem in Fig. 7(b) again for the mLor FDTD. Using Eq. (35), the proper Courant number $C_n=0.1324$ is found out. Figure 9 shows that the mLor FDTD result agrees perfectly with the QCRF FDTD result.

Fig. 9. 1D FDTD simulation corresponding to Fig. 7. In this case, $C_n$=1 for the QCRF FDTD simulation and $C_n=0.1324$ for the mLor FDTD simulation.

Download Full Size | PDF

As a final example, we consider a human fat in 3D. The uniform FDTD space step size is used as $\Delta s=\Delta x=\Delta y=\Delta z=4.37$ mm. The fat cube consisting of 20 cells is surrounded by a host medium with $\varepsilon _{r} = 4$. The QCRF parameters are $A_0=23.40, A_1=2.15\times {10}^{-8}, A_2=3.40\times {10}^{-19}, B_0=1, B_1=3.89\times {10}^{-9}$ and $B_2=8.66\times {10}^{-20}$ [24]. By using the parameter equivalence, we have the unified mLor parameters: $a_0=19.4739, a_1=6.2275\times {10}^{-9}, b_0=1$, $b_1=3.89\times {10}^{-9}, b_2=8.66\times {10}^{-20}$, and $\varepsilon _{r,\infty }=3.9261$. The time step size of the QCRF FDTD simulation is set to as $\Delta t=C_n \Delta s/\left (\sqrt {3}c_0 \right )$ with $C_n=1$. Using the relationship Eq. (35) again, the proper Courant number of the mLor FDTD simulation is found out as $C_n=0.5047$. The Gaussian-modulated sinewave with the center frequency of 1.7 GHz and the half-power fractional bandwidth of 65 % is excited and the ten-layer PML is used to terminate the computational domain. Figure 10 shows the results of QCRF FDTD and mLor FDTD at the center of the fat cube. It is found that the both two FDTD solutions are numerically stable and they are same to each other. When integrating the QCRF model into the mLor FDTD formulation, one should choose the appropriate Courant number that satisfies its numerical stability conditions, using the relationship Eq. (35).

Fig. 10. 3D FDTD simulation for fat.

Download Full Size | PDF

4. Conclusion

The mLor FDTD formulation can unify the Debye, Drude, Lorentz, CCPR, and QCRF models by harnessing mLor model parameter equivalence. In this work, we have studied the numerical stability of the unified mLor FDTD formulation for other dispersion models. Its numerical stability conditions has been investigated in aspects of the Debye, Drude, Lorentz, CCPR, and QCRF model parameters. It has been demonstrated that the numerical stability conditions of the unified mLor FDTD formulation are equivalent to those of the original-model-based dispersive FDTD formulations. Numerical examples has shown that the unified mLor FDTD formulation is numerically stable for the stable FDTD formulation based on the Debye, Drude, Lorentz, or CCPR model but the unified mLor FDTD formulation may suffer from numerical instability for the stable QCRF FDTD formulation. We have stressed that, for the QCRF model, an appropriate Courant number should be used for the unified mLor FDTD simulation, using the relationship between $\nu$ and $\nu _q$ presented in this work. Finally, it should be noted that most of dispersive media in practical applications are usually described by multiple dispersion terms. Our discussion can be extended to the multiple dispersion terms by considering each dispersion term and then adding it.

Funding

National Research Foundation of Korea (2020R1F1A1055444).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. A. Taflove and S. C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain method, (Artech House, Norwood, MA, 2005), 3rd ed.

2. K.-Y. Jung, S. Ju, and F. L. Teixeira, “Application of the modal CFS-PML-FDTD to the analysis of magnetic photonic crystal waveguides,” IEEE Microw. Wireless Compat. Lett. 21, 179–181 (2011). [CrossRef]

3. G. Alsharahi, A. Faize, C. Maftei, M. Bayjja, M. Louzazni, A. Driouach, and A. Khamlichi, “Analysis and modeling of GPR signals to detect cavities: Case studies in Morocco,” J. Electromagn. Eng. Sci. 19, 177–187 (2019). [CrossRef]

4. J.-W. Baek, D.-K. Kim, and K.-Y. Jung, “Finite-difference time-domain modeling for electromagnetic wave analysis of human voxel model at milimeter-wave frequencies,” IEEE Access 7, 3635–3643 (2019). [CrossRef]

5. J. Cho, M.-S. Park, and K.-Y. Jung, “Perfectly matched layer for accurate FDTD for anisotropic magnetized plasma,” J. Electromagn. Eng. Sci. 20, 277–284 (2020). [CrossRef]

6. R. J. Luebbers and F. Hunsberger, “FDTD for Nth-order dispersive media,” IEEE Trans. Antennas Propag. 40, 1297–1301 (1992). [CrossRef]

7. F. L. Teixeira, “Time-domain finite-difference and finite-element methods for Maxwell equations in complex media,” IEEE Trans. Antennas Propag. 56, 2150–2166 (2008). [CrossRef]

8. J. Shibayama, H. Kawai, J. Yamauchi, and H. Nakano, “Analysis of a 3D MIM waveguide-based plasmonic demultiplexer using the TRC-FDTD method,” Opt. Commun. 452, 360–365 (2019). [CrossRef]

9. 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 Trans. Microw. Theory Tech. 49, 377–381 (2001). [CrossRef]

10. Z. Lin, W. Qiu, J. Pu, Y. Fang, and J. Hu, “Accuracy and von Neumann stability of several highly accurate FDTD approaches for modelling Debye-type dielectric dispersion,” IET Microw. Antennas Propagat. 12, 211–216 (2018). [CrossRef]

11. M. Okoniewski and E. Okoniewska, “Drude dispersion in ADE FDTD revisited,” Electron. Lett. 42, 503–504 (2006). [CrossRef]

12. K.-Y. Jung, F. L. Teixeira, and R. M. Reano, “Au/SiO₂ nanoring plasmon waveguides at optical communication band,” J. Lightwave Technol. 25, 2757–2765 (2007). [CrossRef]

13. Z. Lin, P. Ou, Y. Jia, and C. Zhang, “A Highly Accurate FDTD Model for Simulating Loretnz Dielectric Dispersion,” IEEE Photonics Technol. Lett. 21, 1692–1694 (2009). [CrossRef]

14. M. Han, R. W. Dutton, and S. Fan, “Model dispersive media in finite-difference time-domain method with complex-conjugate pole-residue pairs,” IEEE Microw. Wireless Compat. Lett. 16, 119–121 (2006). [CrossRef]

15. 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 Trans. Antennas Propag. 61, 996–999 (2013). [CrossRef]

16. A. Deinega and S. John, “Effective optical response of silicon to sunlight in the finite-difference time-domain method,” Opt. Lett. 37, 112–114 (2012). [CrossRef]

17. K. P. Prokopidis and D. C. Zografopoulos, “Investigation of the stability of ADE-FDTD methods for modified Lorentz media,” IEEE Microw. Wireless Compat. Lett. 24, 659–661 (2014). [CrossRef]

18. J. J. Hu, P. Shum, C. Lu, and G. Ren, “Generalized finite-difference time-domain method utilizing auxiliary differential equations for the full-vectorial analyis of photonic cyrstal fibers,” IEEE Photonics Technol. Lett. 19, 1970–1972 (2007). [CrossRef]

19. I. Udagedara, M. Premaratne, I. D. Rukhlenko, H. T. Hattori, and G. P. Agrawal, “Unified perfectly matched layer for finite-difference time-domain modleing of dipsersive optical materials,” Opt. Express 17, 21179–21190 (2009). [CrossRef]

20. H. Lin, M. F. Pantoja, L. D. Angulo, J. Alvarez, R. G. Martin, and S. G. Garcia, “FDTD modeling of graphene devices using complex conjugate dispersion material model,” IEEE Microw. Wireless Compat. Lett. 22, 612–614 (2012). [CrossRef]

21. H. Lin, D. Xu, M. F. Pantoja, S. G. Garcia, and H.-L. Yang, “Characterization of graphene-based photonic crystal in THz spectrum with finite-difference time domain method,” Chinese Phys. B 23, 094203 (2014). [CrossRef]

22. H. Choi, J.-W. Baek, and K.-Y. Jung, “Numerical stability and accuracy of CCPR-FDTD for dispersive medi,” IEEE Trans. Antennas Propag. 68, 7717–7720 (2020). [CrossRef]

23. K.-Y. Jung, “On the numerical accuracy of finite-difference time-domain dispersive modeling based on a complex quadratic rational function,” Electromagn. 34, 380–391 (2014). [CrossRef]

24. J. Cho, S.-G. Ha, Y. Park, H. Kim, and K.-Y. Jung, “On the numerical stability of finite-difference time-domain for wave propagation in dispersive media using quadratic complex rational function,” Electromagn. 34, 625–632 (2014). [CrossRef]

25. H. Chung, S.-G. Ha, J. Choi, and K.-Y. Jung, “Accurate FDTD modelling for dispersive media using rational function and particle swarm optimization,” Int. J. Electron. 102, 1218–1228 (2015). [CrossRef]

26. H. Chung, J. Cho, S.-G. Ha, S. Ju, and K.-Y. Jung, “Accurate FDTD dispersive modeling for concrete materials,” ETRI J. 35, 915–918 (2013). [CrossRef]

27. 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 Wirel. Propagat. Lett. 15, 425–428 (2016). [CrossRef]

28. H. Chung, K.-Y. Jung, X. T. Tee, and P. Bermel, “Time domain simulation of tandem silicon solar cells with optimal textured light trapping enabled by the quadratic complex rational function,” Opt. Express 22, A818–A832 (2014). [CrossRef]

29. 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,” Opt. Express 23, 873–881 (2015). [CrossRef]

30. H. Chung, C. Zhou, X. T. Lee, K.-Y. Jung, and P. Bermel, “Hybrid dielectric light trapping designs for thin-film CdZnTe/Si tandem cells,” Opt. Express 24, A1008–A1020 (2016). [CrossRef]

31. Z.-H. Lai and J.-F. Kiang, “Dispersive FDTD scheme and surface impledance bounary condition for modeling pulse propagation in very lossy medium,” IEEE Trans. Antennas Propag. 68, 3060–3067 (2020). [CrossRef]

32. K. P. Prokopidis and D. C. Zografopoulos, “A unified FDTD/PML scheme based on critical points for accurate studies of plasmonic structures,” J. Lightwave Technol. 31, 2467–2476 (2013). [CrossRef]

33. K. P. Prokopidis and C. Kalialakis, “Physical interpretation of a modified Lorentz dielectric function for metals based on the Lorentz-Dirac force,” Appl. Phys. B 117, 25–32 (2014). [CrossRef]

34. K. P. Prokopidis and D. C. Zografopoulos, “An ADI-FDTD formulation with modified Lorentz dispersion for the study of plasmonic structures,” IEEE Photonics Technol. Lett. 26, 2267–2270 (2014). [CrossRef]

35. K. P. Prokopidis and D. C. Zografopoulos, “Generalised 3D ADI-FDTD algorithm for the modeling of wave propagation in dispersive media,” Electron. Lett. 53, 1242–1244 (2017). [CrossRef]

36. H. Choi, Y.-H. Kim, J.-W. Baek, and K.-Y. Jung, “Accurate and efficient finite-difference time-domain simulation compared with CCPR model for complex dispersive media,” IEEE Access 7, 160498–160505 (2019). [CrossRef]

37. Y.-H. Kim, H. Choi, J. Cho, and K.-Y. Jung, “FDTD modeling of the accurate electromagnetic wave analysis of graphene,” J. Electrical Eng. Technol. 15, 1281–1286 (2020). [CrossRef]

38. K. P. Prokopidis and D. C. Zografopoulos, “Modeling plasmonic structures using LOD-FDTD methods with accurate dispersion models of metals at optical wavelengths,” J. Lightwave Technol. 35, 193–200 (2017). [CrossRef]

39. H. Choi, J.-W. Baek, and K.-Y. Jung, “Comprehensive study on numerical aspects of modified Lorentz model-based dispersive FDTD formulations,” IEEE Trans. Antennas Propag. 67, 7643–7648 (2019). [CrossRef]

40. D. F. Kelley and R. J. Luebbers, “Piecewise linear recursive convolution for dispersive media using FDTD,” IEEE Trans. Antennas Propag. 44, 792–797 (1996). [CrossRef]

41. D. M. Sullivan, “Z-transform theory and the FDTD method,” IEEE Trans. Antennas Propag. 44, 28–34 (1996). [CrossRef]

42. M. Okoniewski, M. Mrozowski, and M. A. Stuchly, “Simple treatment of multi-term dispersion in FDTD,” IEEE Microw. Guid. Wave Lett. 7, 121–123 (1997). [CrossRef]

43. L. J. Prokopeva, J. D. Borneman, and A. V. Kildishev, “Optical dispersion models for time-domain modeling of metal-dielectric nanostructures,” IEEE Trans. Magn. 47, 1150–1153 (2011). [CrossRef]

44. M. D. Thoreson, J. Fang, A. V. Kildishev, L. J. Prokopeva, P. Nyga, U. K. Chettiar, V. M. Shalaev, and V. P. Drachev, “Fabrication and realistic modeling of three-dimensional metal-dielectric composites,” J. Nanophotonics 5, 051513 (2011). [CrossRef]

45. Q. Ren, H. Bao, S. D. Campbell, L. J. Prokopeva, A. V. Kildishev, and D. H. Werner, “Continuous-discontinuous Galerkin time domain (CDGTD) method with generalized dispersive material (GDM) model for computational photonics,” Opt. Express 26, 29005–29016 (2018). [CrossRef]

Debye	Unified mLor for Debye
$τ \geq 0$	$1 \geq 0$
	$τ \geq 0$
	$0 \geq 0$
$ε_{s} \geq ε_{r, \infty}$	$(ε_{s} - ε_{r, \infty}) τ Δ t^{2} \geq 0$
$ν^{2} \leq 1$	$(ε_{s} - ε_{r, \infty}) τ^{2} ε_{r, \infty} (1 - ν^{2}) Δ t^{2} \geq 0$

Drude	Unified mLor for Drude
$γ \geq 0$	$0 \geq 0$
$γ \geq 0$	$γ \geq 0$
$ν^{2} \leq 1$	$1 - ν^{2} \geq 0$
	$ω_{0}^{2} γ Δ t^{2} + 4 γ ε_{r, \infty} ν^{2} \geq 0$
	$ω_{0}^{2} γ^{2} ε_{r, \infty} (1 - ν^{2}) Δ t^{2} \geq 0$

Lorentz	Unified mLor for Lorentz
$δ \geq 0$	$ω_{0}^{2} \geq 0$
$δ \geq 0$	$2 δ \geq 0$
$ν^{2} \leq 1$	$1 - ν^{2} \geq 0$
$ε_{s} \geq ε_{r, \infty}$	$(ε_{s} - ε_{r, \infty}) ω_{0}^{2} 2 δ Δ t^{2} + 8 δ ε_{r, \infty} ν^{2} \geq 0$
$ε_{s} \geq ε_{r, \infty}$	$4 (ε_{s} - ε_{r, \infty}) ω_{0}^{2} δ^{2} ε_{r, \infty} (1 - ν^{2}) Δ t^{2} \geq 0$

CCPR	Unified mLor for CCPR
$R e (p) \leq 0$	${\| p \|}^{2} \geq 0$
$R e (p) \leq 0$	$- 2 R e (p) \geq 0$
$ν^{2} \leq 1$	$1 - ν^{2} \geq 0$
$Q Δ t^{2} - 4 R e (p) ε_{r, \infty} ν^{2} \geq 0$	$Q Δ t^{2} - 8 R e (p) ε_{r, \infty} ν^{2} \geq 0$
$Q [R e (r) - R e (p) ε_{r, \infty} (1 - ν^{2})] Δ t^{2}$ $- 4 R e (r) R e (p) ε_{r, \infty} ν^{2} \geq 0$	$Q [2 R e (r) - 2 R e (p) ε_{r, \infty} (1 - ν^{2})] Δ t^{2}$ $- 16 R e (r) R e (p) ε_{r, \infty} ν^{2} \geq 0$
$Q = 2 R e (p) R e (r p^{*}) - {\| p \|}^{2} R e (r)$	$Q = 4 R e (p) R e (r p^{*}) - 2 {\| p \|}^{2} R e (r)$

QCRF	Unified mLor for QCRF
$B_{0} \geq 0$	$B_{0} \geq 0$
$B_{1} \geq 0$	$B_{1} \geq 0$
$A_{2} - B_{2} ν_{q}^{2} \geq 0$	$B_{2} (1 - ν^{2}) \geq 0$
$(A_{0} B_{1} - A_{1} B_{0}) Δ t^{2} + 4 ν_{q}^{2} B_{1} B_{2} \geq 0$	$(A_{0} B_{1} - A_{1} B_{0}) Δ t^{2} + 4 ν^{2} B_{1} A_{2} \geq 0$
$\begin{aligned} (A_{0} B_{1} - A_{1} B_{0}) (A_{1} - B_{1} ν_{q}^{2}) Δ t^{2} \\ + 4 (A_{1} B_{2} - A_{2} B_{1}) B_{1} ν_{q}^{2} \geq 0 \end{aligned}$	$\begin{aligned} (A_{0} B_{1} - A_{1} B_{0}) (A_{1} B_{2} - A_{2} B_{1} ν^{2}) Δ t^{2} \\ + 4 (A_{1} B_{2} - A_{2} B_{1}) B_{1} A_{2} ν^{2} \geq 0 \end{aligned}$

Numerical Stability of Modified Lorentz FDTD Unified From Various Dispersion Models

Abstract

1. Introduction

2. Numerical stability conditions of dispersive FDTD formulations

2.1 Debye model

2.2 Drude model

2.3 Lorentz model

2.4 CCPR model

2.5 QCRF model

3. Numerical examples

4. Conclusion

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (10)

Tables (6)

Equations (36)

Optics Express