Consistent FDTD Modeling of Dispersive Dielectric Media via Multiple Debye Terms Derived Rigorously by Padé Approximants

An finite-difference time-domain (FDTD) scheme for simulating wave propagation inside Cole-Cole (CC), Havriliak-Negami (HN), or Raicu (R) dispersive media is proposed. The main difficulty in FDTD implementations for such media is the fact that the time-domain polarization relations are fractional integro-differential equations. To circumvent this difficulty, the dispersive susceptibility of the medium is modeled by Padé approximants. It is proven that under certain conditions the Padé approximant is equal to a weighted sum of Debye terms, while the physical restrictions for positive weights and relaxation times are fulfilled with certainty. Hence, a set of first-order auxiliary differential equations (ADEs), which approximate the original fractional polarization relation, is derived and directly embedded in the FDTD scheme. In contrast to previous works, the derivation of the Debye expansion of the original dispersive permittivity is carried out in a straightforward manner without involving nonlinear optimization techniques. Comparisons between analytical and FDTD results of the reflection and transmission coefficients as well as of the transfer functions inside the dispersive media illustrate the efficiency of the proposed scheme, over wideband frequency domains.


I. INTRODUCTION
T HE finite-difference time-domain (FDTD) method has been proven an accurate technique for simulating electromagnetic wave propagation inside media of various types [1], [2], [3]. In the case of dispersive media, the main objective of any FDTD scheme is the efficient and reliable discretization of the polarization relation between the polarization field and the electric field intensity. Given the relative complex permittivity of the dispersive medium, two basic approaches can be adopted to derive the updating relation for the polarization field. In the first approach, if the permittivity of the medium is expressed in the time domain, then the FDTD method deals with the appropriate discretization of the time-domain convolution integral. In the second approach, if the permittivity is given in the frequency domain, then the inverse Fourier transform of the polarization relation results in an auxiliary differential equation (ADE), which is then discretized by means of finite differences [3], [4]. Manuscript  When the frequency-domain permittivity involves integer powers of j ω (ω is the angular frequency), the derivation of the ADE is rather straightforward and results in a differential equation of integer order. Actually, this is the case with Drude, Debye, and Lorentz models of media [2], [3]. However, the permittivity of many dispersive media, such as biological tissues, polymers, liquid dielectrics, and diols [5], [6], [7], [8], [9], [10], [11], is accurately modeled by means of empirical formulae that involve non-integer powers of j ω. Among the most famous dispersive models are the Cole-Cole (CC) [5], the Davidson-Cole (DC) [6], the Havriliak-Negami (HN) [7], and the Raicu (R) model [8]. In these types of media, the time-domain polarization relation is actually a fractional-order integro-differential equation [12], [13] and, as a result, the application of FDTD schemes meets significant difficulties.
Several FDTD techniques have been proposed that are applicable to the simulation of wave propagation inside the aforementioned dispersive media. Some approaches are based on the approximation of the fractional differential operator, which appears in the time-domain polarization relation, by weighted sums of decaying exponentials whose coefficients are obtained either via nonlinear optimization [14] (for CC model) 1 or the Prony's method [15], [16] (both for CC), [17] (HN). In the FDTD methods proposed in [18] (CC) and [19] (HN and R models), the frequency-domain representation of the fractional differential operator is approximated by polynomial or fractional power expansions of j ω, respectively. It is noted that in [19], the enhanced weighted quantum particle swarm optimization technique is applied. Other FDTD schemes utilize empirical [20] (CC) or particle swarm optimization [21] (CC) techniques to approximate the distribution of relaxation time (DRT) functions associated with the dispersive models; these schemes end up in approximating the dispersive permittivity by weighted sums of Debye terms. The use of Debye expansions that approximate directly the frequency-domain susceptibility function of the dispersive medium has been proposed in [22] (HN), where the approximation is carried out using a hybrid particle swarm optimization approach. Alternatively, the frequency-domain susceptibility function can be modeled by Padé approximants [23], as proposed in [24]. Following this approach, the polarization fractional differential equations are substituted with ADEs of high integer order allowing the direct application of finite differences [25], [26], [27] (CC, DC, and HN). Finally, the FDTD scheme proposed in [28] involves the estimation of the time-domain susceptibility using the fast inverse Laplace transform and the Prony's method. As a result, the dispersive medium is approximated by weighted sums of Debye, and possibly, Lorentz terms. It is mentioned that the appearance of Lorentz terms is incompatible with the DRT representation of the medium susceptibility.
For the DC dispersive medium, it has been rigorously proven that under certain conditions the Padé approximant of the relative complex permittivity is equivalent to a weighted sum of Debye terms [26]. Also, these conditions ensure that the weights and the relaxation times of the corresponding Debye terms are positive, as required by the susceptibility representation via the DRT function. In the present study, the analysis presented in [26] is further extended to cover the cases of CC, HN, and R dispersive models. In particular, Padé approximants are utilized to express the frequency-domain susceptibility functions of CC, HN, and R models. The coefficients of the Padé approximants are easily computed by solving algebraic linear systems of equations. It is proven rigorously that if the degree of the numerator of the Padé approximant is one less than the degree of the denominator, then the Padé approximant is equal to a weighted sum of Debye terms, with strictly positive weights and relaxation times. Also, the parameters of the Debye terms are directly related to the poles of the Padé approximant, and can be computed in a straightforward manner. Thus, the whole approach does not involve the solution of any nonlinear multivariate optimization problem by means of deterministic or stochastic methods, which are not always reliable, and are time and memory consuming. To this end, the time-domain polarization relations of CC, HN, and R media, which are originally fractional integro-differential equations, are approximated by a set of first-order differential equations (one for each Debye term). The approximate polarization relation, through finite-differences discretization, is included in an FDTD scheme.
The rest of this article is organized as follows: In Section II, the models of the dispersive media investigated as well as the polarization relations are presented. The approximation of the susceptibility function by means of weighted sums of Debye terms, which are derived rigorously by Padé approximants, is presented in Section III. The proposed FDTD scheme, along with its memory-storage and computationaltime demands, is described in Section IV. Numerical results illustrating the accuracy of the Padé approximants and the validity of the proposed FDTD scheme are given in Section V. The conclusion is given in Section VI. Finally, the article is supported by three Appendixes, A, B, and C, where necessary definitions, theorems, propositions, and proofs are provided.

A. Relative Complex Permittivities
The relative complex permittivities of the CC, DC, HN, and R models are given by where ε ∞ is the optical relative permittivity, ε is the relative dielectric increment, and σ is the ionic conductivity. The susceptibility function χ( j ω) is given by for the CC [5], DC [6], HN [7], and R [8] models, respectively, where τ 0 > 0 is the characteristic relaxation time. It should be noted that in the cases of the CC, DC, and HN models, ε = ε s −ε ∞ , where ε s stands for the relative static permittivity. It is evident that the R model compared to the other three models is more general. In (5), if γ = 0, then χ R ( j ω) = χ HN ( j ω), if γ = 0 and α = 1, then χ R ( j ω) = χ DC ( j ω), and finally, if γ = 0 and β = 1, then χ R ( j ω) = χ CC ( j ω). Also, if α = 1, then the CC model (2) corresponds to the Debye model with susceptibility function given by

B. Distribution of Relaxation Times
The susceptibility function χ( j ω) can be alternatively expressed as an integral of the form where function G(τ ) stands for the DRT [11]. The DRT function satisfies the conditions: G(τ ) ≥ 0 and ∞ 0 G(τ )dτ = 1. The physical interpretation of (7) is that the CC, DC, HN, and R media can be considered as superpositions of Debye media, which have continuously varying relaxation times, τ , and are weighed by an appropriate distribution function G(τ ). If χ( j ω) is known, then the DRT function is given by [29] By use of (8), the DRT functions for the CC, DC, HN, and R models are obtained, i.e., where ψ HN = arctan[(sin(απ))/((τ/τ 0 ) α + cos(απ))], and where respectively. It is mentioned that the DRT functions (9)- (12) are not used in the present study; they are given just for reference, because they have occasionally appeared in the literature with some misprints.

C. Polarization Relation
In a linear dispersive medium with relative complex permittivity ε r (ω), the phasor of the electric flux density D, is given by where E and P are the phasors of the electric field intensity and the polarization field, respectively. The frequency-domain polarization relation between E and P is written as If we consider any of the four models, namely CC, DC, HN, or R, then (14) is transformed in the time-domain into a fractional integro-differential equation [12], i.e., where D is a time-domain fractional integro-differential operator given by the inverse Fourier transform For example, in the case of the CC model, D = 1 + τ 0 (∂ α /∂t α ) and (15) takes the form of the fractional differential equation Because of the fractional order of the integro-differential polarization relation (15), the direct application of the FDTD method to the simulation of wave propagation in CC, DC, HN, or R media exhibits significant difficulties. However, in the case of the Debye model, the time-domain polarization relation, which is given by is a differential equation of integer order (first-order) and can be discretized easily via the FDTD method [2], [3].

III. APPROXIMATION OF THE SUSCEPTIBILITY FUNCTION
Exploiting the physical interpretation of (7), it is reasonable to approximate the polarization relation (14) using sums of Debye terms. In particular, our intention is to approximate the susceptibility function χ( j ω) by means ofχ( j ω) given bŷ where the parameters τ m and a m need to be estimated. For (19) to be consistent, the parameters τ m (relaxation times) and a m should be positive. In this context, several approaches have been proposed where the parameters of the Debye terms are derived either empirically [20] or by solving nonlinear optimization problems using deterministic [21] or evolutionary algorithms [22].
In the present study, χ( j ω) is approximated by means of Padé approximants. In the following, we focus on CC, HN, and R media, because the case of DC media has already been examined in [26]. In contrast to the DC model, the susceptibility functions χ(s) of the CC, HN, and R models are not analytic at s = 0. Thus, the Padé approximant of χ(s) is based on derivatives of χ(s) at s = ξ > 0 where the derivatives of χ(s) are bounded. In particular, the Padé approximant of χ( j ω) is written aŝ where p n (n = 0, 1, . . . , N) and q m (m = 0, 1, . . . , M) are the Padé coefficients. In (20), [N/M](·) denotes the Padé approximant, which is a rational function whose numerator and denominator are polynomials of Nth and Mth degree, respectively. By substituting (20) into (1), we derive the approximation of the relative complex permittivitŷ The derivation of (20), ]. The particular choice, N = M − 1, will be clarified later. Typically and without loss of generality, we can set ξ = 1. This choice is related to the normalization of the actual frequencies ω with respect to the characteristic frequency of the medium ω 0 = 1/τ 0 [25].
In (21), the selection of the degrees, N and M, of the numerator and denominator polynomials, respectively, has to be consistent with the physical restrictions set by the original relative complex permittivity ε r (ω). For instance, if N > M, then the real and/or the imaginary part ofε r (ω) will be unbounded when the frequency tends to infinity. Also, (20) is actually equal to a weighted sum of Debye terms (19), i.e., where the physical restrictions τ m > 0 and a m > 0 are fulfilled with certainty. It is mentioned that the relaxation times τ m (m = 1, 2, . . . , M) are obtained analytically from the poles of the Padé approximant. This fact is proved in Appendix B. By substituting the original susceptibility, χ( j ω), with its approximation,χ( j ω) (19), into the polarization relation (14), we have where In the time-domain, (24) is transformed into which means that the original fractional integro-differential equation (15) is approximated by M first-order differential equations that can easily be discretized via an FDTD scheme.
IV. FDTD SCHEME The partial time-derivatives that appear in Maxwell's equations as well as in the polarization relations (25) are discretized by use of central differences. For the time instance t = (i + 1/2)t (t is the time step and i is the time-instance index), the finite-difference formulation of (25) is written as where r m = (2τ m − t)/(2τ m + t) and c m = b m t/ (2τ m +t). Then, by substituting (26) into the finite-difference formulation of the Ampere's law we derive the updating formula for the electric field, which is given by If we define auxiliary fields as Q m = (r m − 1)P m and set d m = c m (r m − 1), then (26) is written as while the electric field updating formula (28) is written as To summarize, assuming that the electric field, E, and the auxiliary fields, Q m (m = 1, 2, . . . , M), are known at the time instance i t as well as the magnetic field, H , at the instance (i + 1/2)t, then a single iteration of the FDTD scheme consists of the following steps.
3) Finally, the updated magnetic field is given by which is derived by the Faraday's law. In the proposed FDTD scheme, the memory storage locations required per field spatial component and per computational cell are equal to M + 3; two for E field [current and previous values are needed in (29)], M for the auxiliary fields Q m , and one for the H field. The number of memory locations for storing the parameters, η, ν, κ, r m , and d m is not considered here, and it actually depends on whether the dispersive medium properties vary in space or not. Concerning the computational time associated with (29)-(31), the number of multiplications required for each field spatial component in a cell is equal to 2M + 4; three in (30), 2M in (29), and one in (31).
Finally, the FDTD computational boundary of the dispersive medium is truncated by convolutional perfectly matched layer (CPML), as in [28]. The formulation of the CPML is based on the electric flux density, whose updating formula, at the interface between the medium and the CPML, is given by

A. Accuracy of Padé Approximants
The accuracy of the Padé approximants, when applied to approximate the relative complex permittivities of CC, HN, and R dispersive media is investigated. For all three cases of dispersive media, their common parameters are chosen to be the same, namely ε ∞ = 2, ε = 48, τ 0 = 153 ps, σ = 0, and α = 0.8, while β = 0.6 for the HN and R media,  1, 2, . . . , 10). The relative mean square (rms) error, when the original permittivity ε r = ε r − j ε r is approximated byε r =ε r − jε r (21), is given by where the integrations are carried out along the frequency domain of interest [F L , F H ]. The permittivity approximation error for a wide frequency range (from F L = 0.1 GHz up to F H = 10 GHz), versus the degree of the Padé approximants is presented in Fig. 1. It is evident that the permittivity approximation error is reduced when the degree of the Padé approximants increases, while for M ≥ 6 the rms errors, in all media cases, is less than 1%. The real and the imaginary part of the exact relative permittivities and the Padé approximated ones, when M = 7, are depicted in Fig. 2.

B. Validation of the FDTD Scheme
To validate the efficiency of the proposed FDTD scheme, we consider 1-D problems where the dispersive medium (CC, HN, or R) occupies the region z ≥ 0, whereas the remaining space is free. In each case, the medium is excited by a plane wave, which propagates along the positive z-direction having the electric field polarized in the x-axis. We consider a wideband excitation, where the incident electric field is a modulated Gaussian pulse given by where f c is the central frequency, parameter A controls the bandwidth of the pulse, and u(t) is the step function. The half-power bandwidth of the pulse is HPBW = A √ 2 ln(2)/π.
By post-processing the FDTD simulation results of the electric field, we compare the estimates of the magnitude of the reflection, |R|, and the transmission, |T |, coefficients with their analytical values given by respectively. Furthermore, during the FDTD simulation, the electric field is recorded at two positions inside the dispersive medium, which are separated by distance D along the z-axis. By utilizing the Fourier transform of the two time series, it is possible to estimate the transfer function, i.e., and compare it to the analytically known one given by

1) CC Medium Example:
In the first application, we consider the case of a CC medium with parameters: ε ∞ = 2.  Fig. 3(a). In particular, the rms error (33) between the Padé approximated and the exact relative complex permittivity is e(ε r , ε r ) = 0.12%.
The time step and the grid spacing applied in the proposed FDTD scheme are set equal to t = 0.085 ps and z = 0.05 mm, respectively. The magnitudes of the reflection and transmission coefficients, at the boundary between free space and the CC medium, have been computed by post-processing the electric field time series obtained by the FDTD method. The numerical results and the ones given analytically by (35) and (36) are presented in Fig. 3(b). It should be noted that the rms errors between the analytical and the FDTD-simulated values of the reflection and the transmission coefficient magnitudes are e(|R|, |R|) = 0.22% and e(|T |, |T |) = 0.36%, respectively.
Finally, the simulated result of the transfer function, T(D, ω) (for D = 70z), inside the CC medium, is presented in Fig. 3(c) and compared to its analytical value (38). The rms error between the analytical and the FDTD-simulated transfer function is e(T, T) = 0.52%.
2) HN Medium Example: As a second example, we consider an HN medium with parameters: ε ∞ = 2.61, ε = 10.27, σ = 0.01 Sm −1 , τ 0 = 39.98 ps, α = 0.73, and β = 0.66. The characteristic frequency of the medium is f 0 ≈ 4 GHz. The excitation Gaussian pulse (34) has f c = 10 GHz and HPBW = 50 GHz. A [4/5] Padé approximant (M = 5) is utilized to approximate the relative complex permittivity in the frequency range 1-100 GHz. The exact and the Padé approximated relative complex permittivity of the medium are presented in Fig. 4(a), while the corresponding rms error of the approximation is e(ε r , ε r ) = 0.41%. In the proposed FDTD scheme, t = 0.1 ps and z = 0.06 mm are set. The analytical and the FDTD-simulated values of the reflection and the transmission coefficient magnitudes are presented in Fig. 4(b), while the corresponding rms errors are e(|R|, |R|) = 0.38% and e(|T |, |T |) = 0.42%. Furthermore, the analytical and the FDTD-simulated values of the transfer function T(D, ω), for D = 200z, are presented in Fig. 4(c), while e(T, T) = 0.28%.
In the FDTD scheme, t = 0.05 ps and z = 0.03 mm are set. The analytical and the FDTD-simulated values of the reflection and the transmission coefficient magnitudes are presented in Fig. 5(b), while e(|R|, |R|) = 0.07% and e(|T |, |T |) = 0.16%. Also, the analytical and the FDTD-simulated values of the transfer function T(D, ω), for D = 100z, are illustrated in Fig. 5(c), while e(T, T) = 0.22%.
The numerical results obtained by applying the FDTD method in the three types of dispersive media (CC, HN, and R) are almost identical to the analytical results. Hence, the proposed FDTD scheme is characterized by remarkable accuracy in simulating wave propagation in these media. Furthermore, the agreement between the numerical and the analytical results of the transfer function inside the dispersive media, over a wide frequency range, clearly indicates the absence of any significant numerical dispersion introduced by the FDTD scheme.

VI. CONCLUSION
An alternative time-domain technique, for simulating wave propagation inside CC, HN, and R dispersive media, has been investigated. The objective was to circumvent the difficulties that arise because of the fractional differential polarization relations. It has been shown that it is possible to approximate the susceptibility functions of the media by means of weighted sums of Debye terms, without invoking previously proposed empirical techniques, deterministic/stochastic optimization methods, or sophisticated evolutionary algorithms. The latter are time consuming and do not impose the restriction of positive relaxation times, in a physical way. The present method is based on Padé approximants of the susceptibility function, where the degree of the Padé numerator is one less than the degree of the denominator. It has been proven rigorously that such Padé approximants are actually equivalent to weighted sums of Debye terms with positive relaxation times. Hence, the approximation of the susceptibility functions is derived by a straightforward and sound analysis, without the solution of any nonlinear optimization problem. Using the proposed approximation, the FDTD scheme ends up simulating multiple Debye-term media; its application verifies the accuracy of the method even when broadband field excitations are considered. Finally, we have to underline a significant feature of the presented rigorous derivation of multiple Debye terms. The method is valid not only for the cases of CC, HN, and R media, but it is applicable in any dispersive dielectric medium whose frequency-domain susceptibility is a Stieltjes function with DRT integral representation (7). This is essentially the only restriction imposed on the basis of the method (see Appendix B, and Propositions 1 and 2).

APPENDIX A PADÉ APPROXIMANTS OF χ( j ω)
The objective is to derive Padé approximants of the susceptibility function χ( j ω). To simplify the analysis, we introduce the normalized frequency = ωτ 0 , with respect to the characteristic relaxation time τ 0 , and we define the normalized function χ ( j ) such as For example, χ ( j ) for the R model is given by It is evident that, for CC, HN, and R models, χ (s) is not analytic at s = 0, because its derivatives at s = 0 are not bounded. As a result, it is not possible to derive Padé approximants based on Maclaurin expansion of χ (s). However, all derivatives of χ (s) at s = ξ ∈ R + are bounded. Thus, it is possible to introduce a Padé approximant of the form where q 0 = 1 and all other coefficients p n and q m are derived algebraically [23]. In particular, the evaluation of coefficients p n and q m is based on settingχ (s) equal to the truncated Taylor expansion of χ (s) (at s = ξ ), up to the term of degree 2M − 1. Let be the derivatives of χ (s), which are evaluated at s = ξ . Also, for simplicity, we define

APPENDIX B FROM PADÉ APPROXIMANTS TO MULTIPLE DEBYE TERMS
Given the integral representation (7) of the susceptibility function χ( j ω), it can be shown that the normalized susceptibility function χ ( j ) (A.1) has an integral representation given by In the following, we will prove that the Padé approximant (A.3) of χ (s) is equal to a weighted sum of Debye terms. Since χ (s) is analytic at s = ξ ∈ R + , but not at s = 0, we define h(s) = χ (s + ξ). (B.3)

From (B.1), it is concluded that h(s) is represented by the integral
Obviously, the derivatives of h(s) at s = 0 are equal to the derivatives of χ (s) at s = ξ , which means that they are bounded, i.e., where Proposition 1: The function h(s), which is given by (B.7), is a Stieltjes function and it has a Stieltjes series (see Definitions 1 and 2 in Appendix C). Proof: Based on the fact that H (u) is a positive continuous function, it is concluded that φ(u) is a bounded nondecreasing function taking infinitely many different values on 0 ≤ u ≤ 1/ξ . Thus, φ(u) fulfills the first condition of Definition 1 (see Appendix C). Concerning the moments related to φ(u), we have It is easily shown that, at v = i /ξ , g i (v) takes its maximum value while, from (B.10), we have is defined in the cut s plane, cut along −∞ < s < −λ −1 .
In this case, from Theorem 1, the following corollary is derived.