The HIE-FDTD Method for Simulating Dispersion Media Represented by Drude, Debye, and Lorentz Models

The hybrid implicit–explicit finite-difference time-domain (HIE-FDTD) method is a weakly conditionally stable finite-difference time-domain (FDTD) method that has attracted much attention in recent years. However due to the dispersion media such as water, soil, plasma, biological tissue, optical materials, etc., the application of the HIE-FDTD method is still relatively limited. Therefore, in this paper, the HIE-FDTD method was extended to solve typical dispersion media by combining the Drude, Debye, and Lorentz models with hybrid implicit–explicit difference techniques. The advantage of the presented method is that it only needs to solve a set of equations, and then different dispersion media including water, soil, plasma, biological tissue, and optical materials can be analyzed. The convolutional perfectly matched layer (CPML) boundary condition was introduced to truncate the computational domain. Numerical examples were used to validate the absorbing performance of the CPML boundary and prove the accuracy and computational efficiency of the dispersion HIE-FDTD method proposed in this paper. The simulated results showed that the dispersion HIE-FDTD method could not only obtain accurate calculation results, but also had a much higher computational efficiency than the finite-difference time-domain (FDTD) method.


Introduction
Nanomaterials are widely used in the field of electromagnetism because of their special structure. First, nanomaterials can be used as wave-absorbing materials. On one hand, nano-absorbing materials are widely used on stealth aircraft, tanks, and other military equipment to reduce their probability of detection and destruction and improve the survivability of these targets. On the other hand, nano-absorbing materials are also widely used in daily life such as in the windshields of vehicles and the doors and windows of some public places. Second, nanomaterials have been intensively researched in the field of electromagnetic radiation protection. With the wide application of electromagnetic radiation in communication, transportation, medical, military, and other fields, the human living environment is increasingly polluted by electromagnetic radiation. Research has shown that electromagnetic radiation can affect or even damage biological tissues through various mechanisms. Nanomaterials for electromagnetic radiation protection applications are mainly divided into two categories. One includes radiation protection nanomaterials such as nano-metal materials, nan-oxide materials, and so on. The other is nano protective fabric. In addition, graphene, as the thinnest and strongest new nanomaterial, has been widely used in microwave device design. Graphene is a two-dimensional material made of carbon atoms, so it is only one atom thick. Furthermore, graphene has unique electrical properties, excellent mechanical properties, and thermal stability. Nanomaterials are relatively difficult and costly to prepare, so the co-simulation of nanomaterials and electromagnetic devices is an important prior mean in engineering applications. boundary conditions. Furthermore, the existing HIE-FDTD method can only analyze the Drude, Debye, and Lorentz models separately, and cannot simulate the target composed of two or more dispersion models.
In this paper, we combined the Drude model, Debye model, and Lorentz model with Maxwell's equations, and used auxiliary differential equations to obtain the unified iterative formula of the HIE-FDTD method for these three dispersion models. Thus, the presented method can be used to calculate water, soil, plasma, biological tissue, optical materials, etc. In addition, we introduced the convolutional perfectly matched layer (CPML) boundary condition in the dispersion HIE-FDTD method, so the method can not only simulate closed-domain problems, but also calculate open-domain problems. The high absorbing performance of the CPML was verified, and numerical examples were used to prove the accuracy and computational efficiency of the dispersion HIE-FDTD method. This proves that the dispersion HIE-FDTD not only had accurate calculated results, but also had a much higher computational efficiency than the conventional FDTD method.

Methods and Formulations
The permittivity of many media varies with frequency such as plasma, water, biological tissue, radar absorbing materials, etc. Such media are called dispersive media. In this paper, it was assumed that the magnetic constitutive parameters of these media were independent of the frequency, and only the dielectric coefficients were related to the frequency. The frequency constitutive relation is Therefore, for different dispersive media, we only need to find the relationship between their polarization current and electric field, and then substitute this relationship into Equation (4). Finally the calculation formula of each electric field component can be obtained by using centered finite-difference expressions for the space derivatives in Equation (4). There are three typical frequency models for dispersive media: the Drude model, Debye model, and Lorentz model. Next, we will analyze these models, respectively.

Drude Model
The Drude model was proposed by Paul Drude in 1900 to explain the transport properties of electrons in matter, especially in metal. This model is based on four assumptions.
(1) There is no interaction or collision between free electrons, and the electrons can move independently of each other. (2) The ions in the crystal are stationary, while the electrons can move freely, and there is no electrostatic interaction between electrons and ions. (3) There are collisions between electrons and ions, and the collision is instantaneous. The motion velocity of electrons after the collision is only related to temperature, and has nothing to do with the velocity of the electrons before the collision. Through this collision, the electrons reach a thermal equilibrium with the surrounding environment. (4) The electron collisions follow the Poisson process and the average time interval between two collisions is called relaxation time. The permittivity of the Drude model is expressed as where ε ∞ is the relative permittivity at infinite frequency; ν c is the collision frequency; τ is the relaxation time and τ = 1/ν c ; ω p is the Drude frequency, and the polarizability is The Drude model is often used to describe the dispersion characteristics of plasmas, metals, etc.
By substituting Equation (6) into Equation (3), we obtain Equation (13) can be simplified, and it becomes ∆t where ∆t ∆t ε e f f .

Debye Model
Under the action of an applied electric field, the dielectric will be electrically polarized. Electric polarization is a microscopic process, but it can be comprehensively reflected by macroscopic physical quantity. This quantity is the relative permittivity. In addition, polarization relaxation is an important concept in the analysis of electric polarization phenomenon. Polarization relaxation means that when the external electric field acts on the dielectric, it takes a certain period of time for the polarization intensity P to reach a steady state value. The Debye model considers the polarization relaxation phenomenon in the polarization process and its permittivity is expressed as where ε s is the relative permittivity at static or zero frequency; ν c is the collision frequency; τ is the relaxation time; τ = 1/ν c . The polarizability χ(ω) in the Debye model is where ∆ε = ε s − ε ∞ . The Debye model is often used to describe the dispersion characteristics of media such as soil, water, and human tissue. By substituting Equation (16) into Equation (3), we obtain Equation (17) is the frequency-domain expression between the polarization current and the electric field of the Debye model. An efficient means to obtain the time-domain expression between the polarization current and the electric field is to first multiply both sides of Equation (17) by 1 + jωτ. This gives By applying the transition relation jω → ∂/∂t in Equation (18), we obtain Using the centered finite-difference expression for the time derivative in Equation (19), we obtain Equation (20) can also be expressed in the following form where k p1 = 2τ−∆t 2τ+∆t , β p1 = 2ε 0 ∆ε∆t 2τ+∆t . Equation (21) is the finite-difference equation between the polarization current and the electric field of the Debye model.

Substituting Equation (21) into Equation (4), we obtain
Equation (22) can be simplified, and it becomes ∆t where ∆t ∆t ∆t ε e f f .

Lorentz Model
In the Lorentz model, it is considered that the electron is no longer free and unrestrained. The electron and the ionic are bound to each other by the Coulomb force, and the motion of the electron is simple harmonic vibration. Considering that the linear oscillator has elastic restoring force and damping force, the expression for the permittivity of the Lorentz model is where ε s is the relative permittivity at static or zero frequency; ε ∞ is the relative permittivity at infinite frequency; ω 0 is the natural frequency of the oscillator; ν c is the collision frequency; τ is the relaxation time and τ = 1/ν c . The polarizability χ(ω) of Lorentz model is By substituting Equation (25) into Equation (3), we obtain where ∆ε = ε s − ε ∞ . Equation (26) is the frequency-domain expression between the polarization current and the electric field of the Lorentz model. An efficient means to obtain the time-domain expression between the polarization current and the electric field is to first multiply both sides of Equation (26) by ω 2 0 − ω 2 + 2jων c . This obtains By applying the transition relation jω → ∂/∂t in Equation (27), we obtain Using the centered finite-difference expression for the time derivative in Equation (28), we obtain Equation (29) can also be expressed in the following form (30) is the finite-difference equation between the polarization current and the electric field of the Lorentz model.
By substituting Equation (30) into Equation (4), we obtain Equation (31) can be simplified and becomes ∆t where ∆t ∆t

Standard Model for the Drude, Debye, and Lorentz Models
Based on Equation (14), Equation (23) and Equation (32), the magnetic field curl equation in the time-domain for the Drude, Debye, and Lorentz models can be described as follows: ∆t where Based on Equation (12), Equation (21), and Equation (30), the polarization current for the Drude, Debye, and Lorentz models, respectively, can be calculated as follows: 2∆t , B 4 = k P2 , and B 5 = k P3 . Equation (34) is the auxiliary differential equation that introduces the dispersion property of the dispersive media into the HIE-FDTD method.
The electric field curl equation in the time-domain is the same as that in the linear media, namely,

The HIE-FDTD Method in Dispersion Media
Without loss of generality, it is assumed that the fine structure is along the y direction, so the hybrid implicit-explicit difference technique was applied to the derivative term ∂y/∂t in the HIE-FDTD method. Thus, according to Equation (33), we obtained the formulas of the three electric field components in dispersive media as According to Equation (35), the calculations for the magnetic field are as follows It can be seen that E n+1 x and E n+1 z cannot be calculated directly by using Equations (37) and (38) because they all include the unknown magnetic field defined at the same time. Thus, modified equations were derived from the original equations. Substituting Equations (41) and (39) into Equations (37) and (38), respectively, the equations for E n+1 x and E n+1 z are given as 2ε e f f µ∆x∆y E n+1 where the iterative calculations of the polarization current Thus, the iterative process of the HIE-FDTD method in dispersion media are as follows:  (1). In [12], the weakly conditional stability of the HIE-FDTD method was well-validated. The time step size of the HIE-FDTD method in the y direction was only determined by two spatial increments ∆x and ∆z, namely, where c is the speed of light in the medium and c = 1/ √ εµ. This will significantly reduce the computation time when the simulated structure had fine-scale dimensions in one direction, which will be validated in a later section by using numerical examples.

The CPML Absorption Boundary in the Dispersion HIE-FDTD Method
To truncate the calculation area of the dispersive medium, the CPML absorption boundary can be applied. The CPML formula does not split the field components and the implementation is convenient [22][23][24].
Taking the calculations of E n+1 x and H n+1 z as examples, in the CPML absorption boundary, their formulas are expressed as The values of parameters κ m , a m , and b m (m = x, y, z) in the above calculations refer to the formulas in [22]. The updating equations of other field components in the CPML absorbing boundary can be derived in the same way.

Performance of Absorbing Boundary
To illustrate the CPML termination performance, the simulation of radiation from a small current source in a dispersive space was studied. The total computational domain was 45 × 45 × 45 cells, and a uniform mesh with a cell size ∆x = ∆y = ∆z = 0.015 m was applied. The volume of the dispersive media was 0.075 m × 0.075 m×0.075 m and located in the center of the computational domain. The parameters of the dispersive media are shown in Table 1. A 10-cell thick CPML boundary, which is ten cells far from the dispersive media, was used to terminate all six sides of the lattice.  Figure 1 shows the steady field distribution in the computational domain calculated by using the dispersion HIE-FDTD method. Here, a point source with a frequency of 2 GHz was located in the center of the dispersive media. The time step size was selected as ∆t = 1/ c (1/∆x) 2 + (1/∆z) 2 = 35 ps, which was the maximum time step to satisfy the time stability condition in the HIE-FDTD method. From Figure 1, we can see that in all three models, the radiation was annular distribution, and there was no obvious reflection from the absorbing boundary into the computational domain. This validates the effectiveness of the CPML absorbing boundary. The reflection error of the CPML boundary was quantitatively studied in this paper. Assuming that the excitation source is a plane wave, and the time dependence of the source is where t0 and t1 are constants, and both are equal to 1 × 10 -9 , thus the minimum wavelength of the source is about 0.15 m. The reflection error was computed at an observation point located at cell (30, 30, 30). A reference solution based on an extended lattice was computed in order to isolate the error of the CPML from the grid dispersion error. The relative error was computed as where E std (t) is the value of the electric field at the observation point computed by the HIE-FDTD method in the extended lattice. E(t) is the electric field calculated by the HIE-FDTD method truncated by the CPML boundary.
The reflection error with respect to the time step is shown in Figure 2. As can be seen from this figure, the reflection relative error of CPML was less than −80 dB in all three dispersion models. This validates the absorbing performance of the CPML quantitatively. The reflection error of the CPML boundary was quantitatively studied in this paper. Assuming that the excitation source is a plane wave, and the time dependence of the source is where t 0 and t 1 are constants, and both are equal to 1 × 10 −9 , thus the minimum wavelength of the source is about 0.15 m. The reflection error was computed at an observation point located at cell (30, 30, 30). A reference solution based on an extended lattice was computed in order to isolate the error of the CPML from the grid dispersion error. The relative error was computed as err(dB) = 20 log 10 where E std (t) is the value of the electric field at the observation point computed by the HIE-FDTD method in the extended lattice. E(t) is the electric field calculated by the HIE-FDTD method truncated by the CPML boundary. The reflection error with respect to the time step is shown in Figure 2. As can be seen from this figure, the reflection relative error of CPML was less than −80 dB in all three dispersion models. This validates the absorbing performance of the CPML quantitatively.

Accuracy and Efficiency
To validate the computational accuracy and efficiency of the dispersion HIE-FDTD method, the scattered field of a multilayer complex dielectric plate was analyzed. The specific structure is shown in Figure 3. A four-layer dielectric plate, whose length, width, and height were 0.3 m, 0.3 mm and 0.045 m, respectively, was separated by three dispersive dielectric layers. The thickness of each dispersive medium layer was 3 mm, and was composed of three kinds of dispersive media: Drude, Debye, and Lorentz, respectively. The parameters of the dispersion media were same as those in Table 1. A plane wave propagating in the -z direction was vertically incident on the dielectric plate. The polarization of the electromagnetic wave was along y direction, and the time dependence of the wave is shown in Equation (52). Two observation points were located before and behind the plate at 0.03 m and 0.06 m, respectively. A 10-cell thick CPML layer, which was 10 cells far from the plate, was used to terminate all six sides of the lattice.
The electric field Ey(t) at the observation points were calculated by using the FDTD method and the HIE-FDTD method. The minimum spatial cell Δymin was 0.0005 m in the dispersion media. In the other computational domain, spatial cells Δx = Δy = Δz = 0.015 m were used. According to the CFL stability condition, the maximum time step size of the FDTD method was 1.66 ps. However, in the HIE-FDTD method, the time step size was not related to the minimum cell size, so its time step size could reach 35.35 ps. The calculated results of these two methods are shown in Figure 4. To validate the accuracy of the methods, the results simulated by using the typical software CST Studio Suite 2018 were also plotted in Figure 4. It can be seen from Figure 4 that the calculated results agreed well with each other.

Accuracy and Efficiency
To validate the computational accuracy and efficiency of the dispersion HIE-FDTD method, the scattered field of a multilayer complex dielectric plate was analyzed. The specific structure is shown in Figure 3. A four-layer dielectric plate, whose length, width, and height were 0.3 m, 0.3 mm and 0.045 m, respectively, was separated by three dispersive dielectric layers. The thickness of each dispersive medium layer was 3 mm, and was composed of three kinds of dispersive media: Drude, Debye, and Lorentz, respectively. The parameters of the dispersion media were same as those in Table 1. A plane wave propagating in the -z direction was vertically incident on the dielectric plate. The polarization of the electromagnetic wave was along y direction, and the time dependence of the wave is shown in Equation (52). Two observation points were located before and behind the plate at 0.03 m and 0.06 m, respectively. A 10-cell thick CPML layer, which was 10 cells far from the plate, was used to terminate all six sides of the lattice.

Accuracy and Efficiency
To validate the computational accuracy and efficiency of the dispersion HIE-FDTD method, the scattered field of a multilayer complex dielectric plate was analyzed. The specific structure is shown in Figure 3. A four-layer dielectric plate, whose length, width, and height were 0.3 m, 0.3 mm and 0.045 m, respectively, was separated by three dispersive dielectric layers. The thickness of each dispersive medium layer was 3 mm, and was composed of three kinds of dispersive media: Drude, Debye, and Lorentz, respectively. The parameters of the dispersion media were same as those in Table 1. A plane wave propagating in the -z direction was vertically incident on the dielectric plate. The polarization of the electromagnetic wave was along y direction, and the time dependence of the wave is shown in Equation (52). Two observation points were located before and behind the plate at 0.03 m and 0.06 m, respectively. A 10-cell thick CPML layer, which was 10 cells far from the plate, was used to terminate all six sides of the lattice.
The electric field Ey(t) at the observation points were calculated by using the FDTD method and the HIE-FDTD method. The minimum spatial cell Δymin was 0.0005 m in the dispersion media. In the other computational domain, spatial cells Δx = Δy = Δz = 0.015 m were used. According to the CFL stability condition, the maximum time step size of the FDTD method was 1.66 ps. However, in the HIE-FDTD method, the time step size was not related to the minimum cell size, so its time step size could reach 35.35 ps. The calculated results of these two methods are shown in Figure 4. To validate the accuracy of the methods, the results simulated by using the typical software CST Studio Suite 2018 were also plotted in Figure 4. It can be seen from Figure 4 that the calculated results agreed well with each other.  The electric field E y (t) at the observation points were calculated by using the FDTD method and the HIE-FDTD method. The minimum spatial cell ∆y min was 0.0005 m in the dispersion media. In the other computational domain, spatial cells ∆x = ∆y = ∆z = 0.015 m were used. According to the CFL stability condition, the maximum time step size of the FDTD method was 1.66 ps. However, in the HIE-FDTD method, the time step size was not related to the minimum cell size, so its time step size could reach 35.35 ps. The calculated results of these two methods are shown in Figure 4. To validate the accuracy of the methods, the results simulated by using the typical software CST Studio Suite 2018 were also plotted in Figure 4. It can be seen from Figure 4 that the calculated results agreed well with each other.  To complete the simulation above, the time step size of the FDTD method was 1.66 ps, and its total simulation time was 81.7 minutes. In the HIE-FDTD method, the time step size was 35.35 ps, so the total simulation time of the HIE-FDTD method was only 5.4 minutes, which was only 1/15 of that of the FDTD method. This validates that the HIE-FDTD method had a much obvious higher computational efficiency than that of the FDTD method.
For the multilayer complex dielectric plate, as shown in Figure 3, when the thickness of each dispersive medium layer was 3 nm, the electric field Ey(t) at the observation points calculated by using the dispersion HIE-FDTD method proposed in this paper is shown in Figure 5. In this simulation, the minimum spatial cell Δymin was 0.5 nm in the dispersion media. According to the CFL condition, the time step size of the FDTD method was only 1.66 × 10 -6 ps. It can be deduced that the simulation time of the FDTD method was about 8.17 × 10 7 minutes, which was about 56736 days. The time step size of the HIE-FDTD method was still 35.35 ps because it was not limited by the grid size in the direction of the media layer thickness, and the simulation time was still 5.4 minutes. It can be seen that the simulation time required by the HIE-FDTD method was much shorter than that of the FDTD To complete the simulation above, the time step size of the FDTD method was 1.66 ps, and its total simulation time was 81.7 minutes. In the HIE-FDTD method, the time step size was 35.35 ps, so the total simulation time of the HIE-FDTD method was only 5.4 minutes, which was only 1/15 of that of the FDTD method. This validates that the HIE-FDTD method had a much obvious higher computational efficiency than that of the FDTD method.
For the multilayer complex dielectric plate, as shown in Figure 3, when the thickness of each dispersive medium layer was 3 nm, the electric field E y (t) at the observation points calculated by using the dispersion HIE-FDTD method proposed in this paper is shown in Figure 5.  In this simulation, the minimum spatial cell ∆y min was 0.5 nm in the dispersion media. According to the CFL condition, the time step size of the FDTD method was only 1.66 × 10 −6 ps. It can be deduced that the simulation time of the FDTD method was about 8.17 × 10 7 minutes, which was about 56736 days. The time step size of the HIE-FDTD method was still 35.35 ps because it was not limited by the grid size in the direction of the media layer thickness, and the simulation time was still 5.4 minutes. It can be seen that the simulation time required by the HIE-FDTD method was much shorter than that of the FDTD method, and the computational efficiency of the HIE-FDTD method was greatly improved.

Conclusions
In this paper, to solve the dispersive media effectively, the Maxwell equations and three dispersion models were combined, and the hybrid explicit-implicit difference technique was used to obtain the HIE-FDTD method. By introducing the CPML absorption boundary, the method was extended to solve open-domain dispersion problems. Numerical examples were used to verify the absorption performance of the absorption boundary, the calculation accuracy, and calculation efficiency of the algorithm. Since the presented method had a higher computational efficiency than the traditional FDTD method, it could be well applied to solve dispersive media with fine structures such as water, soil, plasma, biological tissue, optical materials, etc.