Fractional Derivative Modification of a Drude Model

A modification of the Drude dispersive model based on fractional time derivative is presented. The dielectric susceptibility is calculated analytically and simulated numerically, showing a good agreement between theoretical description and numerical results. The absorption coefficient and wave vector -- key parameters describing the propagation of waves in such a medium -- are shown to follow a power law in the frequency domain, which is a common phenomenon in many real life applications. The introduction of two separate parameters provides a more flexible model than some other approaches found in literature and is well suited for numerical implementation.


I. INTRODUCTION
Factional calculus -branch of mathematics devoted to study non-integer order derivatives and integrals, becomes an increasingly popular tool for analysis of various problems in physics, ranging from quantum mechanics and cosmology [1] to electric circuits [2] and electromagnetic wave propagation [3]. An extensive review of recent developments is presented in [4].
The application of fractional derivative in theoretical [5] and numerical [6] description of mechanical waves has been extensively studied, including the propagation of acoustic waves [7]. One of the key motives behind the development of fractional models is the observation that wave attenuation in many systems follows a power law with a non-integer exponent [8,10], which cannot be described with standard time domain partial differential equations. Furthermore, many physical media exhibit hereditary features [11], where some physical property is dependent on the history of its previous values; fractional derivative provides a tool to analyse such systems that carry information about their present as well as past states [12][13][14].
Motivated by these developments, we propose an extension of the Drude model, which is one of the basic tools for describing the electric permittivity of metals. To better fit the experimental observations in complex media, several extensions of Drude model have been proposed, ranging from introduction of frequency-dependent parameters [15] to the use of fractional derivative [17]. The approach presented here introduces two fractional derivatives to the equations of motion describing the medium, resulting in a more general and flexible model. Moreover, the implementation of the medium description in Finite-Difference Time-Domain (FDTD) scheme [16] is presented. The accuracy of simulation results is discussed. * david.ziemkiewicz@utp.edu.pl

II. DRUDE MODEL
We assume that the medium (usually metal) contains some concentration n of free charges e, with effective mass m. The equation describing their motion under the influence of external field E is where r is the charge position and Γ describes dissipative processes. By introducing the polarization vector P = ne r, damping constant γ = Γ/m and plasma frequency ω 2 p = ne 2 /m, one obtains which, in the time domain, for harmonic wave E = E 0 e −iωt , P = P 0 e −iωt yields which is a standard expression for modelling electric susceptibility of metals [15]. In our modified model, we introduce the fractional derivative operator D to the Eq.
(2) in all places where time derivative is used. In particular, we transform the Eq. (2) into the form where D β denotes fractional derivative of the order β and α, β are real parameters. Standard relation is obtained for γ α = α = β = 1. In the frequency domain, the above relation leads to the susceptibility Note that the constants γ α , γ β ensure the proper dimensionality of the equation, e. g. γ α ω α+1 = ω 2 and γ β ω β = ω 2 , so that χ remains a dimensionless quantity.
In contrast to other extensions of Drude model, where dissipation constant and/or plasma frequency are frequency-dependent quantities [15], we introduce additional frequency dependence by directly modifying the exponent of the frequency.

III. NUMERICAL IMPLEMENTATION
The Finite-Difference Time-Domain method (FDTD) consists of dividing the simulation space into grid points and calculating the values of electric and magnetic field at those points with evolution equations derived directly from Maxwell's equations, with some finite time step ∆t [16]. Due to the fact that the algorithm is based on first principles, has a well-known sources of numerical errors and is easy to implement in parallel computing, it is one of the leading tools for analysis of complex optical and plasmonic systems. One of the methods to include complex media in FDTD simulations is ADE (Axillary Differential Equations) approach [9]; the medium polarization P is computed by numerically solving a partial differential equation describing its evolution in time. For regular Drude model, one can derive the evolution equation for P from Eq. (2). In the FDTD scheme, one has a set of discrete values of polarization P t , t = 1, 2, 3.... To obtain the first-order approximation of the first time derivative, one can use the relationṡ In a similar manner, one can define the second derivative asP In our modified approach, one has to first define the fractional derivative operator; from the standpoint of computation with discrete time step ∆t, the most convenient definition is the truncated Grunwald-Letnikov derivative [18] where in numerical implementation ∆t is set to some finite value and N is a suitably large number such that the components of the sum are negligibly small. Note that the left-sided derivative is used, where only previous values of f (t) are needed. In other words, the model is casual. In terms of discrete values of P , one obtains By applying the above definition and the averaging procedure in Eqs. (6) and (7) to the Eq. (4), one obtains the relations for derivatives (10) and the equation of motion By expanding and rearranging the terms, one obtains the evolution relation which relates the new value of P t+1 to the current value of P t , the electric field E t and the history P t−k . For the case of α = β = 1, α k = β k = 0 for k > 1 and the relation simplifies to the standard form presented in [9].

IV. RESULTS
A. Calculation of susceptibility Figure 1 shows the real part of susceptibility calculated from analytical solution (5) and obtained in FDTD simulation. The medium parameters are set to ω p = 0.3, γ β = 0.1, γ α = 1. On the Fig. 1 a) the value of β is set to 1 and 0 ≤ α ≤ 1.5. Like in a standard Drude model, the susceptibility is negative and for α ≤ 1, χ → −∞ when ω → 0. For higher values of α, a distinct minimum of susceptibility occurs. One can see that the accuracy is excellent in the medium frequency range, e. g. ω ≈ ω p ; at the low frequency limit, the large susceptibility corresponds to short wavelength, which approaches the finite spatial step of the simulation. Additionally, waves are highly absorbed in this region, further decreasing prediction accuracy. In the very high frequency regime, accuracy is limited by the fact that the wave period approaches the finite time step. In this region, χ → 0 and the value of α has a significant effect on how quickly this asymptote is approached. Fig. 1 b) shows the susceptibility calculated for 0 ≤ β ≤ 1 and α = 1. Here, one can observe a transition from pure Drude model β = 1 to a resonant model with resonance frequency of ω = ω p (β = 0). Again, the calculations become inaccurate in the regions of high absorption at ω → 0 and near the resonance. Overall, one can conclude that the model allows for a very significant alteration of the dispersion relation while retaining good stability and accuracy. An important factor for efficient implementation of the proposed model is the number of terms N in the sums in Eq. (12). The Fig. 2 shows the relative error of the numerical susceptibility as a function of frequency, calculated with various numbers of the memory terms P t−k . One can observe that in general, the error initially quickly decreases with frequency for ω < 0.5ω p , with a further, slower decrease at higher ω. The number of terms has a high impact on the low frequency accuracy. In the solutions with low number of terms, the error has a significant minimum; the calculated values below and above the frequency where the minimum occurs are overestimated and underestimated, accordingly. This can be attributed to the so-called numerical dispersion [16], which introduces small frequency-dependent term to the susceptibility regardless of the medium model. The effect is easily visible on the Fig. 1 a) for the case of α = 1.5; the numerical results (dots) overestimate the susceptibility for ω < ω p and underestimate it beyond that frequency. In majority of the spectral range, the error is greater than in the case of a standard Drude model (Fig. 2, green line) by a factor of 2-3. Apart from the large minimum, the error also exhibits a slight oscillatory behaviour, with amplitude of oscillations reducing with increasing number of terms. One may conclude that the minimum number of summation terms that provide satisfactory accuracy is N ∼ 15, with relatively little gain from increasing N further. Moreover, the smooth error function with no short-term changes indicates that the model can be easily tuned to maintain a very high accuracy within a chosen, limited spectral range; in the standard FDTD implementation, it is straightforward to add a frequency-independent term ∞ to the dielectric susceptibility (ω) = 1 + χ(ω) [16], regardless of the medium model introduced with ADE. By doing so, one can shift the spectral region where the numerical results match the theory exactly. Alternatively, additional Drude-like term can be introduced to ADE to counteract the numerical dispersion. Finally, due to the fact that parameters α and β can be changed continuously, one can optimize the accuracy by using slightly different values for theoretical and numerical calculations. The importance of the memory terms and hereditary properties of the medium decreases as the parameters α, β approach the limit of standard Drude model, e.g. α = β = 1; in such a case, the necessary number of memory terms for accurate computation is reduced. In many applications, only a small modification of the medium is considered [17].
In a standard, two-dimensional implementation of FDTD, one has to define three scalar field values at every grid point (for example, two components of magnetic field vector H x , H y and one component of electric field E z perpendicular to the simulation plane). [16] The inclusion of Drude or Drude-Lorentz dispersion model as described in [9] adds the medium polarization P . Specifically, 1 current (P t ) and 1 previous (P t−1 ) value of polarization is needed, resulting in 5 scalar values per grid point. The example implementation of fractional Drude model with 12-term memory adds another 10 scalar values, increasing the memory requirement by a factor of 3. However, it should be noted that the increase is needed only for the part of computational domain that contains the fractional model medium.
Finally, it should be noted that the coefficients α k , β k in the sums in Eq. (12) need to be computed once for any given values of α, β and the calculation of the weighted sum in Eq. (9) is essentially a discrete convolution operation, which can be subject to various numerical optimizations. An extensive discussion of numerical application of convolution to calculate a fractional derivative is presented in [19].
Another advantage of the presented model is its tunability. By allowing adiabatic changes of α and β in the time domain simulation, one achieve a dynamical tuning of the optical properties of the medium.

B. Wave propagation in fractional medium
Dispersive properties of the medium described by the function χ(ω) influence the propagation of electromagnetic waves through the material. It directly affects the permittivity = 1 + χ, refraction index n = √ and the value of the wave vector k = ωn/c. Assuming a harmonic wave with frequency ω and wavevector k, from the relation (5) one obtains The above relation is nontrivial for real parameters α, β due to the fact that both terms in denominator introduce separate, frequency-dependent contributions to both real and imaginary part of k. One can simplify the problem by assuming that the medium is a small modification of Drude model, with α = 1 and β ∼ 1; in such a case, in the high frequency limit one obtains Thus, the imaginary part of k is In a similar manner, one can derive the limit for α = 0 which is Im k ∼ ω 0 . As mentioned in [7] and [10], there is a demand on dispersion models where the attenuation (which is proportional to Im k) follows a frequency power law. Fig. 3 shows numerically calculated imaginary part of wave vector in the whole spectral range. The discussed limiting cases are shown with color lines, while the transitional spectra with fractional values of α, β are shown in gray. There are two distinct regions ω < ω p and ω > ω p ; below plasma frequency, the susceptibility is negative and correspondingly ≈ 0, Re k ≈ 0; the wavevector is thus purely imaginary, which corresponds to highly absorbed, evanescent waves. The parameter β has a negligible impact on the absorption in this range. Above the ω p , the absorption follows a power law (straight line), with the results consistent with Eq. (15). The dependence on the parameter α is more complicated; the steep reduction of absorption at ω ≈ ω p for α = 1 becomes more gradual for smaller values of α. In the limit of α → 0, the absorption becomes almost constant. In contrast to β, the parameter α has a significant impact on the imaginary part of k in the region ω < ω p . While waves with such frequencies are highly attenuated when propagating through the medium, the negative value of permittivity (see Fig. 1) allows for formation of surface plasmons. These collective electron oscillations are highly sensitive to changes of medium properties [20], which makes them a particularly promising field of study where the proposed model could be applied. One of the prospective applications of fractional derivative model and resulting power-law dissipation are metallic nanoparticle chains [10]. The above mentioned power-law dependence extends to other material functions such as conductivity In the paper [17] the author proposed a fractionalderivative model based on Eq. (2), which is transformed to a dimensionless, first-order differential equation describing particle velocity. Our model is consistent with results in [17] when one sets β = 1 and α ≈ 1, resulting in a small modification of the standard Drude model. Calculation results for such parameters are shown on the Fig. 4. One can notice that while a small modification of α has a little impact on the σ for ω ∼ ω p , it dramatically changes the high-frequency behaviour of the medium. As the order of fractional derivative decreases, the slope of the high frequency asymptote increases, in agreement with [17]. The impact of β is much smaller in the high frequency range, but more significant for ω < ω p .

V. CONCLUSIONS
A novel, two-parameter modification of the Drude model based on Grunwald-Letnikov fractional derivative has been presented. The analytical formulas for basic optical functions such as susceptibility and wave vector have been derived and discussed in the context of wave propagation in the medium. A numerical implementation in a FDTD method has been realized, expanding a wellknown ADE method and taking advantage of efficient discrete convolution computation. The numerical complexity and accuracy of the approach was discussed. The results indicate that the proposed model is highly flexible and applicable to a wide variety of optical and plasmonic systems, allowing for modelling of other modified Drude models as well as dynamic modulation of medium properties.