A Novel Piecewise Linear Recursive Convolution Approach for Dispersive Media Using the Finite-Difference Time-Domain Method

Two novel methods for implementing recursively the convolution between the electric field and a time dependent electric susceptibility function in the finite-difference time domain (FDTD) method are presented. Both resulting algorithms are straightforward to implement and employ an inclusive susceptibility function which holds as special cases the Lorentz, Debye, and Drude media relaxations. The accuracy of the new proposed algorithms is found to be systematically improved when compared to existing standard piecewise linear recursive convolution (PLRC) approaches, it is conjectured that the reason for this improvement is that the new proposed algorithms do not make any assumptions about the time variation of the polarization density in each time interval; no finite difference or semi-implicit schemes are used for the calculation of the polarization density. The only assumption that these two new methods make is that the first time derivative of the electric field is constant within each FDTD time interval.


I. INTRODUCTION
T HE finite-difference time domain method (FDTD) [1], [2] is a very popular numerical technique for solving Maxwell's equations for a wide range of different media, including materials with frequency dependent properties.In a number of practical applications of FDTD modelling, single and multi-pole Lorentz, Debye and Drude functions are widely used because they allow us to simulate the electric susceptibility of a range of materials amongst them: water [3], human tissues [4]- [7] cold plasma [8], [9], gold [10], soils [11]- [17].Electric susceptibility is defined as the function relating the polarization density to the electric field by .In essence, both the real and imaginary parts of a complex frequency depended electric permittivity variation are included in the electric susceptibility function.In this paper, for simplicity we define in the same way as used in [18] and [19] the electric susceptibility from the definition of electric flux density as .In such definition the effects of any frequency-independent parts of the medium's relative permittivity are included with the relative permittivity at infinite frequency term.
A number of different methodologies have been suggested for numerically implementing dispersive materials into the FDTD.The resulting methods can be roughly divided in three categories: 1) auxiliary differential equation (ADE) methods [20]- [25]; 2) Z-transform methods [26]- [28]; and 3) the recursive convolution methods [18], [19], [29]- [36].A systematic review of the methods related to ADE and recursive convolution can be found in [2] and [37].In addition, an interesting review about the modelling techniques used to simulate Lorentz media can be found in [32].
Inclusive algorithms, with a uniform implementation for a wide range of materials, are very attractive especially in situations where materials with different dispersion mechanisms need to be modelled.An inclusive ADE algorithm for modelling Lorentz, Debye, and Drude media is presented in [23].The algorithm needs two additional variables to be stored per pole for both Drude, Lorentz, Debye, and conductive term mechanisms.A more efficient ADE inclusive algorithm which is based on a complex-conjugate pole-residue method is presented in [22].This algorithm uses an inclusive electric susceptibility function which holds as special cases the Debye, Lorentz and Drude medium.
Piecewise linear recursive convolution (PLRC) [18] is an efficient method for dealing with dispersive materials.PLRC assumes that the electric field has a piecewise linear behavior and also uses a central difference scheme in order to calculate the derivative of the polarization density in time.In other words, it implicitly assumes that the polarization density can be accurately simulated by a second order polynomial in each time interval.PLRC is a widely used method, and one of the key reasons for its popularity, is that it is an accurate algorithm which can simulate materials with an inclusive susceptibility function which holds as special cases the Lorentz, Debye, and Drude media.This makes the implementation of dispersive materials in FDTD codes easy and practical and at the same time retains its computational efficiency [22], [36].Standard PLRC methodologies require two additional variables for each Lorentz medium and one additional variable for each Drude or Debye medium.Complex media like the ones described by the Havriliak-Negami equation can be approximated by multi-Debye functions [6], [17] and can be implemented very effectively using PLRC [38].Also, inclusive algorithms like PLRC and complex-conjugate pole-residuals method [22] that 0018-926X © 2014 IEEE.Personal use is permitted, but republication/redistribution requires IEEE permission.
can simulate susceptibility functions which hold as special cases the Debye, Lorentz, and Drude relaxations are more computationally efficient for modelling complex materials like Ag [22] and graphene [39].
In this work, we present two new novel recursive convolution based methodologies which in their development make only the assumption that the first time derivative of the electric field is constant in each time interval-in other words that the electric field is piecewise linear between time steps-and do not make any assumptions about the variation in time of the polarization density.Instead, an analytical calculation of the polarization density is used subject to the assumption that the electric field has a piecewise linear behavior in each time interval.As a result of reducing the number of the assumptions made in the derivation of the recursive algorithm, the proposed methods in this paper are found to be systematically improved compared with PLRC for Debye, Lorentz, and Drude media and at the same time they preserve the main advantages of recursive convolution techniques.ADE interpretations of the proposed methods are also given in the Appendixes.

II. POLARIZATION DENSITY METHOD
Maxwell's equations for an isotropic, linear, non-magnetic, dispersive medium are given by (1) (2) (3) where and are the electric permittivity and the magnetic permeability of vacuum, respectively, is the relative electric permittivity for infinity frequency, is the magnetic field, is the electric field, is the polarization density for each dispersive pole, is the electrical conductivity, and is the total number of dispersive poles.The susceptibility function used here is an inclusive function which holds as special cases the Debye, Drude, and Lorentz relaxations [36] (4) We then define to be a complex function having as its real part the polarization density (5) where and can, in general, be both real and complex numbers; if they are both real then .For Debye and Drude media both and are real, while for a Lorentz medium is purely imaginary and is generally a complex number with both non-zero real and imaginary parts.In the proposed algorithms can also be complex with both non-zero real and imaginary parts and could be a purely imaginary or real number (similarly to the complex-conjugate pole-residue method [22]).
Equation ( 5) can be written as (6) which in discritized form is given by (7) Because , it can be easily shown that ( 7) can be rewritten as (8) To solve the integral in (8), subject to the constrain that the electric field has a piecewise linear behavior in each time interval, we apply the integration by parts rule instead of using the approach suggested in [18].In order to apply the integration by parts rule we rewrite the integral in ( 8) as (9) Applying the integration by parts rule in the integral of the right hand side of ( 9) and assuming that the first derivative of the electric field is constant over each time interval results in (10) The assumption of the linearity of the electric field which is directly implied by the assumption that its first derivative with respect to time is constant, is consistent with the development of an algorithm that is suited for a general second order in time FDTD scheme.
Substituting (10) into (9) and then using the result into (8) and using a central finite difference approximation to calculate the time derivative of the electric field in (10) results in (11) (12) (13) If we assume that the derivative of the polarization density can be accurately obtained using a central difference scheme (14) and by further substituting the real part of ( 14) into ( 1) and using a central difference scheme for the calculation of and a semi-implicit technique for the calculation of we arrive at Equations ( 16)-( 18) can be shown to be the same as the coefficients of the update equation of the electric field (15) using PLRC [18] (with a semi-implicit scheme to calculate in (1)) for solving recursively the convolution between the electric field and the electric susceptibility function.From (14) it is evident that a central difference scheme is used in order to calculate the time derivative of the polarization density.In (15) and (16) is used as a temporary holding value for the summation containing the terms and therefore there is no need for extra memory storage.
In order to improve upon the standard PLRC algorithm instead of using a central difference scheme as in (14), the time derivative of the polarization density at the required time instance is calculated analytically.To achieve this we have to rewrite (8) for ( ) as (19) The integral in ( 19) can be calculated using the integration by parts rule similarly to (8).In order to do this (19) is cast as (20) Therefore, the integral in (20) can be obtained by (21) Because of the assumption of the linearity of the electric field in each time interval the term in ( 21) is equal to (22) From ( 22), (21) and (20) it can be easily shown that equals (23) where (24) (25) For , and .The derivative of the polarization density at between two time intervals ( and ) can be written as (26) From ( 23), ( 24)-( 26) it can be shown that the analytic expression of the derivative subject to the constrain that the electric field is piecewise linear in each time interval equals to (27) The derivatives with respect to in (27) can be calculated analytically using (24) and (25).Substituting the real part of ( 27) into (1) and using the same discretization as in (15) yields the same equation ( 15) but with different coefficients (28) (29) (30) where (31) (32) The algorithm stores in a temporary memory variable the value of and then calculates from (15) using the coefficients given in ( 28)- (30), and .Subsequently is calculated from (23) using , , the coefficients given in ( 24) and ( 25), and the temporary memory variable which stores . is used to calculate from (28) which is used for the calculation of .When and are real numbers only one additional memory variable is needed to be stored, if or are complex numbers one complex variable is required (i.e., both real and imaginary parts need to be stored).In addition, an equivalent auxiliary differential equation (ADE) derivation of this method based on Laplace transformation is presented in Appendix A.

III. CURRENT DENSITY METHOD
In this alternative method the dispersive properties are expressed by using apparent current density sources.Therefore (1) is written as (33) where (34) We define as a complex function with its real part being the current density (35) The identity [40] illustrates the difference between the polarization density and the current density method.In the case of the polarization density method the convolution between the electric susceptibility function and the electric field is recursively calculated and subsequently the derivative-left-hand side of the identity-is analytically obtained.In the second case of current density method the convolution between the electric susceptibility and the first time derivative of the electric field-right-hand side of the identity-is calculated.Both methods assume that the electric field is linear between time steps.
In the same way that (20) was derived, assuming that the derivative of the electric field is constant in each time interval it can be shown that (36) Using a central difference scheme to calculate the derivative of the electric field in time, (36) becomes (37) Substituting (37) to (34) and subsequently into (33) for and discretizing as was done in ( 15) yields ( 15) with (38) (39) (40) where (41) The algorithm works in the same way as the polarization density method described previously.To find for the next iteration we calculate (37) for .It is evident that no semi-implicit approximations are used to calculate , the current density term is calculated analytically subject to the constrain that the electric field has a piece-wise linear behavior.From ( 37), ( 38)- (40) it is evident that one less pre-calculated variable is needed to be stored compared with PLRC and the polarization density method.An alternative auxiliary differential equation (ADE) interpretation of this algorithm based on a power series method is given in Appendix B.

A. Debye Relaxation
For a Debye medium the susceptibility function is equal to [19] (42) where , , and are the relative permittivity for zero and infinite frequency, respectively, and is the relaxation time.From ( 4) and ( 42) it is evident that (43) (44) In order to test the two proposed algorithms we compare numerically obtained reflection coefficients to analytical ones from a simple 1-D problem of a plane wave normally incident from vacuum onto a two-pole Debye medium.For the numerical calculations an 1-D FDTD algorithm was used with 0.3 mm and 1 ps.The multi-Debye medium characteristics were , , , 271 ps, 10.8 ps, and .Fig. 1 presents the comparison between the analytical and the numerical reflection coefficients and relative permittivity for the polarization density method.The analytical and the numerical results are in excellent agreement showing the accuracy of the proposed algorithm.Results obtained by the current density method are almost identical to the ones obtained from the polarization density method although the two formulations are not.The polarization density method is slightly better by an almost negligible amount.
The two proposed algorithms in this paper are compared with the PLRC [18], the complex-conjugate pole-residue method [22], and the Trapezoidal Recursive Convolution [31] method.These are all algorithms that can implement an inclusive susceptibility function which holds as special cases the Debye, Lorentz, and Drude media relaxations which make them very attractive for modeling complex materials like graphene [39].A 1-D FDTD model, as used previously, is employed to test all the other methods mentioned above and the errors between the analytical and the numerical reflection coefficients are shown in Fig. 2. The overall errors for all the methods have the same order of magnitude because all of them are computationally efficient algorithms with approximately second order accuracy.Regarding the amplitude error the differences between the methods proposed in this paper are negligible, for the phase error as shown in Fig. 2 the new methods proposed here performed better when compared to the PLRC, TRC and complex-conjugate method.

B. Lorentz Relaxation
For a Lorentz medium the susceptibility function is equal to [29] (45) where is the resonant frequency, is the damping factor, and .From ( 4) and ( 45) yields To assess the accuracy of the formulations we compare numerically obtained reflection coefficients to analytical ones from a plane wave incident from vacuum onto a two-pole Lorentz medium.The 1-D FDTD parameters were 30 pm and 0.1 as.The properties of the Lorentz medium were set to , , PHz, , , and .Fig. 3 shows the comparison between the analytical and the numerical reflection coefficients and relative permittivity obtained using the current density method proposed in this paper.The same example is repeated using the polarization density, TRC, PLRC  and the complex-conjugate method and the errors are shown in Fig. 4.Both the amplitude and the phase related errors are reduced for the new algorithms with the polarization density to be slightly better when compared to the current density method.
Comparisons between reflection coefficients of plane waves incident normally from vacuum onto a dispersive medium have been widely used in the literature to evaluate the accuracy of different algorithms.This raises the question of how compre-   fractional bandwidth equals to 5. The received field is sampled at a location away from the source's position along -axis.Fig. 5 shows the geometry of the Lorentzian waveguide model and presents a number of snapshots that illustrate how the pulse is propagating inside it.
Fig. 6 compares the numerical-using the current density method-and the analytical calculated at the receiver position.In order to further compare the errors between the methods presented here and other approaches the same test was repeated using polarization density, PLRC, TRC, and the complex-conjugate method.By summing the absolute values of the differences between the analytical and numerical results and subsequently normalizing these errors to the maximum error (in this case to the error from complex-conjugate method), a comparison could be made of the relative performance of each approach.Fig. 7 presents these normalized errors obtained from all the methods above.The methods proposed in this paper show reduced overall errors when compared with PLRC, TRC, and the complex-conjugate (C-C) methods.Again it has been found that the polarization density method is marginally better when compared with the current density method.

C. Drude Relaxation
For Drude medium the susceptibility function is equal to [41] (51) From ( 1 [41].Fig. 8 shows the numerical and the analytical reflection coefficients and the numerical and analytical relative permittivity calculated using the polarization density method.The current density method performed equally well.The same numerical test is repeated using TRC, PLRC, and the Complex-Conjugate methods and the errors are shown in Fig. 9. Similarly to the results obtained from the Lorentz relaxation medium numerical test, both the amplitude and the phase related errors for the Drude medium have been found to be less for the new proposed algorithms with the polarization density to perform again slightly better when compared to the current density method.V. CONCLUSION Two novel and elegant ways to implement dispersive media into the FDTD are presented.In both of the proposed methods the polarization density is calculated analytically subject to the constraint that the time derivative of the electric field is constant in each time interval, no central difference scheme like PLRC or TRC neither a semi-implicit technique is used to calculate the polarization density.Therefore, the proposed methods are more accurate than PLRC and TRC with the same computational efficiency while the major advantage of recursive convolution techniques-being the implementation of an inclusive electric susceptibility function which holds as special cases Debye, Drude, and Lorentz media-is been preserved.The superiority in accuracy of the proposed algorithms compared with TRC [31], PLRC [18], and complex-conjugate pole-residuals method [ From Laplace transform theory [40] we know that (61) (63) Transforming (60) back to time domain using (61)-( 63) we obtain (23).Therefore, one can easily arrive at exactly the same update equations as presented in Section II for the polarization density method.
In order to solve this equation subject to the assumption that the derivative of the electric field is constant in each time interval and without making any assumptions about the nature of , we use a power series method [42].First we define in the interval as an infinity order polynomial (i.e., a power series) (67) where are vectors and .Substituting (67) to (66) yields (68) Equation ( 68) can be written in a more practical form as (69) This method is applied in each time interval separately so, the right term can be considered as constant.The initial condition of this differential equation in each time interval is for .From that initial condition and (67) can be easily shown that .Repeating the process in a similar way for all the coefficients we derive the inclusive formula for (75) Substituting ( 75) to (67) and using the initial condition (i.e., ) yields (76) The Maclaurin series for is given by ( 77) Using ( 77), (76) and a central difference scheme for the derivative of the electrical field in time in (76) we derive (37).As a result, one can then easily arrive at exactly the same update equations as presented in Section III for the current density method.From (67), it is clear that no simplifications are made about the order of in the interval .

Fig. 1 .
Fig. 1.Multi-Debye medium.A) Comparison between the numerical (polarization density method) and analytical reflection coefficients.B) Comparison between the numerical (polarization density method) and the analytical relative permittivity.

Fig. 3 .
Fig. 3. Multi-Lorentz medium.A) Comparison between numerical (current density method) and analytical reflection coefficients.B) Comparison between numerical (current density method) and analytical relative permittivity.
hensive 1-D examples are.Therefore, in this paper numerical and analytical results are compared from a 2-D Lorentzian waveguide.The characteristics of the Lorentz medium in the waveguide are , , Ps and .The FDTD characteristics of the model are 30 pm and 70.711 zs (Courant limit).The waveguide's width is .An electric line source is placed in the center of the waveguide and the excitation pulse is a Gaussian-modulated sinusoidal function with central frequency 100 PHz and

Fig. 5 .
Fig. 5. Pulse propagating into a Lorentzian waveguide.The contours illustrate the electric field .

Fig. 6 .
Fig.6.Comparison between numerical and analytical results using the current density method.

Fig. 7 .
Fig. 7. Normalized summation of the absolute error between numerical and analytical results using PLRC, TRC, Complex-Conjugate (C-C), and current density method.
previous examples regarding the Debye and Lorentz media we compare the numerical and the analytical reflection coefficients from a plane wave incident onto a Drude medium excitation pulse did not have low frequency content and overcomes the obstacle explained in

Fig. 8 .
Fig. 8. Drude medium.A) Comparison between the numerical (polarization density method)l and analytical reflection coefficients.B) Comparison between the numerical (polarization density method)l and analytical relative permittivity.
transforming (65) to time domain a first order differential equation is obtained From (69) by equating the coefficients of the same power of on the left and right side yields 22] (i.e., with the methods that use an inclusive susceptibility function) is demonstrated through numerical examples.