On the Application of Fractional Derivative Operator Theory to the Electromagnetic Modeling of Frequency Dispersive Media

: Fractional derivative operators are finding applications in a wide variety of fields with their ability to better model certain phenomena exhibiting spatial and temporal nonlocality. One area in which these operators are applicable is in the field of electromagnetism, thereby modelling transient wave propagation in complex media. To apply fractional derivative operators to electromagnetic problems, the operator must adhere to certain principles, like the trigonometric functions invariance property. The Grünwald–Letnikov and Marchaud fractional derivative operators comply with these principles and therefore could be applied. The fractional derivative arises when modelling frequency-dispersive dielectric media. The time-domain convolution integral in the relation between the electric displacement and the polarisation density, containing an empirical extension of the Debye model, is approximated directly. A common approach is to recursively update the convolution integral by approximating the time series by a truncated sum of decaying exponentials, with the coefficients found through means of optimisation or fitting. The finite-difference time-domain schemes using this approach have shown to be more computationally efficient compared to other approaches using auxiliary differential equation methods.


Introduction
Fractional calculus is the theory of derivatives and integrals of an arbitrary order [1].In other words, it is a generalisation of classical integer order calculus.Therefore, it possesses all the capabilities of classical calculus with additional features.For example, fractional operators can better model certain phenomena with spatial and temporal nonlocality [2].
The theory of fractional calculus was postulated in 1695 shortly after classical calculus [3] and developed thereafter in 1823 [4].Its applications were limited, compared to classical calculus, due to the complexity of the equations, thereby making it extremely difficult or impossible to solve analytically.However, after the development of the computer, it became possible to compute fractional derivatives and integrals numerically.
Because of this, they are finding applications in a wide variety of fields: viscoelasticity, control theory, biology, physics, economics, electromagnetism, and many more [5,6].Most physical phenomena in these domains are governed by partial differential equations (PDEs).However, certain phenomena can be better modelled using fractional-order derivatives.
One such phenomenon is in the field of electromagnetism.Maxwell's equations are a set of coupled partial differential equations (PDEs), which, along with the Lorenz force, governs classical electrodynamics, although these equations only employ integer calculus.On the other hand, heuristic measurement-based models describing the relaxation of electromagnetic waves in frequency-dispersive materials fit very well with fractional calculus [7].
Frequency-dispersive media are often encountered in nature, so being able to model these materials well is essential to understanding and predicting their behaviour.Some prominent applications of this include the following: ground-penetrating radar [8,9], biological media [10][11][12][13], and nanophotonics [14,15].
Algorithms implementing frequency-dispersive media in the framework of the finitedifference time-domain (FDTD) method have been developed for this purpose [16].The most commonly used methods are recursive convolution (RC) methods [17,18], auxiliary differential equation (ADE) methods [19], and Z-transform methods [20], where the ADE and Z-transform methods are very similar in their approach.
The RC method approximates a convolution integral, in the relation between the electric displacement and the polarisation density, as a discrete summation solved using a recursive procedure.In the ADE method, the approach is to transform the complex permittivity from the frequency domain into the time domain using the inverse Fourier transform, which results in a differential equation in the time domain relating the electric displacement to the electric field.Then, the derivatives are approximated with finite differences based on Taylor series expansions [21].The Z-transform method converts the time-domain convolution integral into a multiplication in the z domain, thus leading to a direct FDTD implementation.Essentially, the Z-transform is a more general version of the discrete Fourier transform.Although the RC and ADE approaches have similar accuracy, the RC method can be more computationally efficient [22], while ADE methods are easier to formulate and directly implement in the framework of the FDTD.
To better model dielectric media in the FDTD framework, more accurate dispersion laws should be used.Such empirical dispersion models have been derived based on the Debye model.Namely, these include the Cole-Cole [23], , , and Raicu models [10].However, implementing these models in the FDTD framework is challenging.
For example, the RC method cannot be directly incorporated into the FDTD scheme when using these empirical models.This is because the RC methods need to transform the frequency-domain permittivity into a time-domain function in order to compute the convolution integral, which, when using these empirical models, leads to complicated timedomain expressions of the permittivity function containing fractional derivatives.There are analytical solutions for some fractional orders of the fractional differential equation [26].However, this is not the case for modelling transient wave propagation, which requires numerical treatment.
It is also more difficult to utilise empirical dispersion models in the framework of the FDTD using the ADE or Z-transform methods, as the fractional-order conversion schemes lead to nonrational Z formulae [27].
In this paper, an overview of fractional derivative operators is provided in Section 2, with an emphasis on the operators most suitable to the application of EM problems.A review of the different methods to implement these empirical dispersion models in the framework of the FDTD is given in Section 3, thereby focusing on methods using fractional derivatives.

Overview of Operators
There are many different definitions of fractional derivative operators: Riemann-Liouville, Caputo, Hilfer, Grünwald-Letnikov, Caputo-Fabrizio [28], Atangana-Baleanu [29], Coimbra and many more [30].The reason so many fractional derivative operators exist is because there are no unique solutions.So, the choice depends on the context or application, more specifically, the functional space under consideration.Additionally, unlike the integer derivative operator, a fractional derivative is not a local operator, thereby leading to many different options.
Further increasing the number of operators, there are left-and right-sided fractional derivatives [30].Take the Riemann-Liouville derivatives as example.Suppose the function f (t) is defined on the interval [a, b], where a and b can be infinite; then, for t ≥ a, the Riemann-Liouville left-sided derivative [30] can be given as and for t ≤ b, the Riemann-Liouville right-sided derivative [30] reads Here, D α denotes the fractional derivative operator raised to the order α, γ(•) is the canonical gamma function [31], α ∈ C, and ℜ(α) ∈ (n − 1, n], where n ∈ Z ≥0 and ℜ(•) define the real part of a complex number.A physical interpretation can be given assuming t is time, then f (t) represents a process changing with time.Taking t as the present moment in time and letting τ < t, f (τ) describes the past states of f (left-sided derivative), whereas in letting τ > t, f (τ) describes the future states of f (right-sided derivative).Therefore, to reduce the number of definitions under consideration, it is reasonable to only consider left-sided derivatives when describing physical phenomena, as the dependence of present states on the future is generally not known (causality).
Since it is not convenient to have multiple definitions of fractional derivative operators, there have been studies working towards generalising them [32,33].They introduce a general framework whereby changing a kernel and normalisation constant gives a known fractional derivative.It is important to note that the generalisation reduces to two forms, with one being the generalised Riemann-Liouville type derivative: (GRL) [32]: and the other being the generalised Caputo-type derivative (GC) [32]: where D α represents the fractional derivative operator, 0 ≤ α ≤ 1, k(t, α) is the kernel function, N(α) is the normalisation function, and f (t) is some continuous, defined, and differentiable function.By filling in the kernel and normalisation functions given in Table 1 into the generalised framework, known fractional derivatives can be obtained [32].
Table 1.List of some kernel functions for fractional derivatives [32].

Fractional Derivatives k(t, α) N(α)
After the development of multiple definitions of fractional derivatives, there has been much discussion on the conditions classifying them [34][35][36], with many based on the general laws formulated by Volterra [37].There are properties that should be met by the operator, but generalising them is not trivial, because depending on the application, one type of kernel may be more useful than another [2].For example, a stretched exponential model was used to fit the data in a study on heterogeneous water diffusion in the human brain using magnetic resonance imaging (MRI) [38].So, a fractional derivative operator with a stretched exponential kernel would be more suited to this application than a power law kernel.
The power law kernels of the Riemann-Liouville and Caputo derivatives are widely used, as they can model many phenomena [39].But the power law kernel is an infinite system with a so-called heavy tail (slow decay rate); thus, it should not be used for systems which do not share this property [40].Another distinguishing property is that the kernel bears a singularity at one of the end points of the interval.However, singular kernels cannot accurately describe certain phenomena, like the those related to material heterogeneities [28] or dispersive phenomena [41].This is because the power law kernel is more affected by events further in the past.So, another type of kernel must be adopted for these situations.
One such nonsingular kernel is the Caputo-Fabrizio kernel, which uses an exponential function that does not contain any singularities [28].Exponential kernels are less affected by past events; in other words, they exhibit stronger fading memory than the power law kernel [42].In other words, the kernel is more local.Consequently, exponential kernels do not meet all the conditions required to be a fractional operator [43].However, it can still be used as a fractional operator in applications which do not require this nonlocality principle to be satisfied, like for the Gurtin-Pipkin theory of heat conduction with finite wave speeds [32].
Another nonsingular kernel is the Atangana-Baleanu kernel, which is based on the Mittag-Leffler function, E α and was developed to solve the nonlocality mathematical issue of the Caputo-Fabrizio kernel [29].
There has been some discussion on whether nonsingular kernels should be used [44][45][46].Disregarding the mathematical properties, these kernels have been shown to model certain phenomena well [47,48].Although certain kernels might not be appicable for frequency-dispersive media, these kernels might be very suitable for other problems.

Specific Operators for Electromagnetism
The Riemann-Liouville type can be seen as a derivative of the convolution of a given function with the power law kernel.Similarly, the Caputo type can be seen as the convolution of the derivative of a given function with the power law kernel.In some situations, one type might be more useful than the other.
For example, when solving initial value problems for fractional-order differential equations, the Caputo-type derivative could be preferred [49].It is possible to solve these equations using a Riemann-Liouville-type derivative, but the initial conditions required for a solution are themselves of noninteger order, which is physically unacceptable.However, if the initial conditions are zero, then a physical meaning can be attributed [50].And for this case, the Riemann-Liouville-type derivative and Caputo-type derivative are equivalent.
It should be noted that the generalisation in Ref. [32] is restricted to fractional orders of 0 < α < 1.A generalised Caputo-type fractional derivative that also holds for higher orders is given in Ref. [33]: where α ∈ [0, 1],n ≤ (n + α) ≤ n + 1, with n ∈ N. A generalised Riemann-Liouvilletype fractional derivative which holds for higher orders was excluded because it required physically unacceptable initial conditions [33].But, as discussed previously, depending on the application, physically acceptable initial conditions can be realised [50].So, a higherorder semigeneralised Riemann-Liouville-type fractional derivative should be possible within certain conditions [33].However, unless there is a way to derive physical meaning to initial conditions expressed as fractional derivatives, then a fully generalised framework for the Riemann-Liouville-type fractional derivatives does not make sense physically [33,50].
It has been suggested that the so-called Grünwald-Letnikov and Marchaud derivatives, which are equivalent for some cases [51], are more suited for EM problems than the Riemann-Liouville and Caputo derivatives with a finite base point [52].This is due to the fact that the Riemann-Liouville and Caputo derivatives do not satisfy all the properties deemed necessary when treating electromagnetic problems.These properties are linearity the index law and trigonometric functions invariance where j = √ −1, ω is the angular frequency, and α, β ∈ R. The so-called trigonometric function invariance property, regarding phasor representations of signals, is specific to electromagnetic theory [52].
Additionally, the Grünwald-Letnikov and Marchaud derivatives are a generalisation of the Riemann-Liouville derivative and can be applied to a wider class of functions [53].Both are based on the use of differences, with the Grünwald-Letnikov derivative being based on extending the binomial expansion of the Newtonian derivative to fractional orders.
To summarise, when dealing with problems in the field of electromagnetism, only fractional derivative operators that satisfy these aforementioned properties should be used.

Numerical Approximations of Dispersive Media
Fractional derivative operators applied to the fractional partial differential Maxwell's equations cannot be solved analytically.Therefore, computational electromagnetics (CEMs) techniques must be employed.The most common techniques are the method of moments (MoM) [54], the finite element method (FEM) [55], and the FDTD method [56].These methods can be divided in two categories: the frequency domain (FEM and MoM) and the time domain (FDTD).
FEM divides the domain into finite elements and solves for polynomial approximations to the unknown functions on each element.It is especially useful when dealing with complicated geometries where higher-order approximations are desired.Another frequency domain method is MoM.This integral equation-based boundary element method reformulates the equations governing a boundary value problem into a matrix equation whose numerical solution can become computationally intensive and require significant memory and time as the size of the problem increases.This is particularly true for problems involving large structures or very fine details, where the number of unknowns can become very large.
The aforementioned frequency domain approaches provide a solution for a single frequency.Since one often has easy access to the properties of complex dispersive materials as a function of frequency, this allows for easily setting up the equations for a single frequency.However, the main drawback is that it is hard to find a solution to ultrawide bandwidth (UWB) excitations, since then the method would need to be iterated over many frequencies, which can be very computationally expensive.For UWB excitations, it is often beneficial to work in the time domain, in which an UWB excitation requires only one simulation.Therefore, the FDTD framework would be the most convenient choice, since it works in the time domain, and complex materials can also be conveniently modelled in the time domain [57].Therefore, the FDTD method is most commonly used.
For completeness, it should be noted that variants of FEM and MoM also exist for the time domain.In principle, the finite-element time-domain (FETD) method could also be employed.Even though the FETD method has shown potential for modelling dispersive media, it is more complicated to implement and is not as developed as the FDTD method when considering complex media [58].Additionally, work has also been done on the time-domain MoM [59,60], but this work is still far away from implementing dispersive media altogehter.Therefore, we focus on the FDTD method, since it is currently the most-used approach.

Finite-Difference Time-Domain
The FDTD method is based upon the application of the finite-difference approximation for both temporal and spatial discretization.Its explicit formulation, as introduced by Yee in Ref. [61], is second-order accurate, robust, simple to implement, and displays reduced computational cost.The numerical errors arising from the FDTD method are well known [62].More importantly, the fact that the FDTD method treats nonlinear and impulsive behaviour naturally is especially useful for the UWB EM characterisation of complex media.
When solving problems using the FDTD method, it was initially assumed that the medium is isotropic, nondispersive, linear, and time-invariant, as this was valid for many cases.However, as the FDTD method grew in popularity, algorithms were devised to solve Maxwell's equations in complex media [63,64], where these assumptions are not absolutely necessary [65].

Dispersive Dielectric Media
Dispersive dielectric media exhibit a permittivity function that is frequency-dependent when exposed to a time-harmonic electromagnetic field.This is because there is a slight delay between changes in the material polarisation and changes in the electric field.The intrinsic relationship between the electric displacement, D, and the electric field intensity, E, for a causal, linear, isotropic dispersive dielectric is given as where t represents time, ϵ 0 is the electric permittivity in a vacuum, and χ e is the electric susceptibility.And in the time-harmonic form, it is given as where ω is the angular frequency, and ϵ(ω) is the complex permittivity function.
There are many models which characterise this complex permittivity function, with the most simple being the Debye model, which can be expressed as where ϵ r (ω) is the relative complex permittivity function, ϵ r ∞ is the relative permittivity at infinite frequency, ϵ r s is the static (zero frequency) relative permittivity, j is the imaginary number, j 2 = −1, τ is the relaxation time, and σ is the ionic conductivity.Although this model is a suitable representation for some dielectrics, there are a number of dielectrics exhibiting exponential relaxations, which show significant deviations from the Debye model [66].To allow for more accurate fitting of the experimental data, empirical extensions of the Debye model have been developed.Namely, in increasing order of generalisation, the dielectric dispersion laws, for media exhibiting multiple relaxation processes, are defined in the Cole-Cole model [23] and the Raicu model [10] where 0 < α i , β i , and ν i < 1 are adjustable parameters that determine the shape of the complex permittivity curve, and i =1,2, . . ., N, where N is the number of relaxation processes transpiring in the dielectric material.

Empirical Dispersion Laws Applied to FDTD
The need for fractional derivatives arises when the procedure is based on approximating the time-domain convolution integral in (9).Most methods try to recursively update the convolution integral by approximating the time series by a truncated sum of decaying exponentials [67], with the coefficients found by some means of optimisation or fitting, like enhanced weighted quantum particle swarm optimisation [22] or the damped least-squares method [68], for example.But, other methods also exist, thereby approximating the complex permittivity function using a general series expansion [69], a binomial series [70], or by using a fractional polynomial series [7].
Although fractional derivatives are more computationally expensive to solve due to their nonlocality and difficulties in discretisation, there have been advancements to overcome these drawbacks [22,69].For example, a formulation was developed using the Riemann-Liouville derivative to find the optimal truncation of the binomial series relevant to the fractional derivative [70], meaning the method now requires less memory, which makes the scheme more computationally efficient, as the relation between the accuracy of the method and the length of the weighted summation of exponentials is known.Additionally, in Ref. [71], a technique was developed which reduced the computational storage for the Havriliak-Negami model from O(N) to O(log N).
These approaches are based on the time-domain empirical dispersion laws in the FDTD method by applying fractional derivatives.The other methods deal with the frequencydomain permittivity function resulting in an ADE, which is appropriately discretised.This escapes directly dealing with the fractional derivative.It does this by transforming the time-domain response into the frequency or Z domain, like in Ref. [72].Similar to the previous approach, a means to find the coefficients of the FDTD update equations are used, for example, using a least-squares fitting method [73].
The Z-transform method was implemented for fractional derivatives in Ref. [74] by expanding the nonrational Z formulae using a truncated Taylor series, which is converted into a difference equation.This truncation, which is a form of the "short memory" principle [5], causes a clipping of the "history" of the function so that only the "recent memory" is considered.However, the stability and accuracy of the method is not well defined, so it is unclear at which point to truncate the power series to achieve a certain level of accuracy.
Many FDTD algorithms constructed for Debye media can efficiently accommodate multiple Debye functions.In view of this, a method was developed to approximate the frequency-domain complex permittivity using a weighted sum of Debye functions by utilising a hybrid particle swarm optimisation approach [75].Even though the stability and accuracy properties are the same as the schemes incorporating the Debye dielectric model [76], due to the band-limited nature of the approximation and not taking the singularity in the fractional relation into account, the method is unable to compute the impulse or step response of spatially complex dispersive media.
The Pade approximant is another technique to approximate a function using a rational function [77].It often gives a better approximation of a function compared to truncating its Taylor series, and it could be used where the Taylor series does not converge.Accordingly, instead of using a truncated Taylor series, like in Ref. [74], a method which uses the Pade approximation of the complex frequency-dependent permittivity in the Z domain was developed [78,79].However, the relation between the order of the polynomials and the accuracy of the results is not stated.This dependency is important because there is a tradeoff between accuracy and memory requirements.The higher order polynomials give more accurate results but at the cost of memory storage.Additionally, the method in Ref. [79] derives a dispersion relation that is dependent on the model parameters, meaning the stability analysis must be performed every time the model parameters change.
Another technique which used the Pade approximant found that under certain conditions, it is equal to a weighted sum of Debye terms [80], thereby building on the work in Ref. [81].This means that there is no need to use any nonlinear multivariate optimisation, thus saving time and memory.They started by filling in the susceptibility function into the polarisation relation between the electric field and the polarisation density, P, and taking the inverse Fourier transform (IFT), which, for the Cole-Cole media, results in the equation from Ref. [80]: To get around dealing with this fractional derivative, Ref. [80] approximated χ e (ω) as a sum of weighted Debye terms: where M is the number of relaxation times, and a m and τ m are parameters that can be estimated empirically [82] by using nonlinear optimisation techniques [83] or by using Pade approximants [80].This approximates ( 16) by M discretized first-order differential equations: thus making it easier to discretize and directly implement in the framework of the FDTD.Finally, the update equations are derived by applying finite difference to Ampere's law and Faraday's law.The limits of this method are not restricted to only the models presented in Section 3.2, but instead can be used for any dispersive dielectric media where χ e (ω) is a Stieltjes function with the integral representation of the distribution of relaxation times (DRTs) of the form defined in Ref. [80]: where G(τ) is the DRTs function [84].
Using a slightly different approach to the other ADE methods, Ref. [72] incorporated ideas from digital filter theory.By transforming (10) into the Z domain, they derived the following: χ e (z) becomes the transfer function from E(z) and D(z).Therefore, χ e (z) must be known.First, Ref. [72] used the fast inverse Laplace transform (FILT) to transform χ e (ω) into the time domain at a specific time step.Then, Prony's method [85] is used to find the signal parameters to derive χ e (z): where A l , p l , B k , C k , r k , and q k are real values calculated using a procedure formulated in Ref. [72], N l is the number of real poles, and N k is the number of complex pole pairs.From this, the update equations can be found by using the Z-transform method [86], thereby applying finite differences to Ampere's law and Faraday's law and substituting and rearranging.Plus, the methods can be used for any dispersive dielectric media where the inverse Laplace transform can be solved numerically.However, it is interesting that the generalised pencil-of-function (GPOF) method is not used over Prony's method, as it is more computationally efficient and robust [87], which should result in a more efficient and robust FDTD scheme.The downside is that the proposed method is more computationally expensive compared to other ADE methods [88] and may be prone to numerical instability, as it relies heavily on the FILT and Prony's method, which are approximate methods.
To summarise, the ADE method eliminates the need to formulate the time-domain representation of the fractional derivatives, meaning the application of finite differences within the FDTD framework is easy and straightforward.Despite this advantage, the method may require storing the field quantities for quite a few previous time steps for each Yee cell edge.This may lead to large memory requirements and computation costs compared to methods which approximate the time-domain fractional derivative, as they store the field quantities relevant only to the previous time step [83].Although dealing with fractional derivatives leads to more complex formulations, it has been shown to be worth the effort, thereby displaying more accurate results than those using the ADE method [22].Furthermore, there was limited knowledge regarding the optimum truncation point for the ADE methods, whereas the time-domain methods have explored this, so the computational costs could be estimated [89].
Therefore, the fractional derivative methods have the potential to be the better choice over the ADE method.A more detailed description of how they are implemented in the FDTD follows.

Fractional Derivatives in FDTD
Methods involving fractional derivatives in the framework of the FDTD all follow a similar procedure, only differing in the way in which the fractional derivative is approximated.First the definition of the complex permittivity function or electric susceptibility function is given.For example, Ref. [11] suggested a formulation for the Cole-Cole media: This is plugged into the relation between the auxiliary displacement current density, J, and the time derivative of the electric field or the relation between the electric field and the induced dielectric polarisation density, P, thereby giving the equation from Ref. [11]: The equation is rearranged and transformed into the time domain using the IFT.Since the auxiliary displacement current density is needed for the update equations, methods using the relation between P and E need to differentiate both sides of the equation, as At this point, the methods differ from each other, thereby dealing with the fractional derivative caused by the IFT in different ways.
For example, in Ref. [11], the resulting equation is Discretizing the formulation using the finite difference method and expanding the vector terms using a semiimplicit approximation [62] leads to the update equation for the auxiliary displacement current density: where and where Due to the fractional derivative, the update of J n+1 requires the storage of all previous auxiliary displacement current density terms in time.Therefore, a recursive update scheme was used with the parameters found through an optimisation technique, like least squares for example.This approximation significantly reduces the computational memory, as only the previous time instant, J n , is required.The electric field update equation was found by substitutions and rearranging Ampere's law, with the curl of the magnetic field approximated with the standard Yee equations.In Ref. [69], the electric susceptibility function was defined as where A p is a real coefficient based on the material parameters.Following the same steps as previously stated, this leads to where ⌈•⌉ is the ceiling function, meaning −1 < α p − ⌈α p ⌉ < 0. This was discretized using a second-order accurate finite-difference scheme evaluated at t = m∆t.To solve the fractional derivative, the principle of the FDTD was used, which assumes constant electromagnetic quantities over each time interval, ∆t, to be equal to the value of the quantity at the centre of ∆t: Using a binomial expansion of the Newtonian derivative gives the following: The update equations can be found by discretising, substituting, and rearranging Ampere's law, along with Faraday's law.
In Ref. [7], the complex permittivity is approximated using a general fractional polynomial series approximation: where a k , α k , b l , and β l are real-valued parameters, which leads to the resulting equation [7]: where D β l t and D 1+α k t are the fractional derivative operators of order β l and 1 + α k , respectively.The equation is discretised using a second-order finite-difference scheme and evaluated using a semiimplicit approximation, with the discretized fractional derivatives approximated as in Ref. [7]: where Ψ l,q is the auxiliary current density vector, S l is a summation of coefficients, and s l,q is a coefficient following an exponential expansion.A similar expression was obtained for D 1+α k t E m [7].In short, this expression was obtained by applying the central finitedifference approximation to the fractional derivative using the principle of FDTD in (31) and solving the integral part of the discretized fractional derivative, I l : where an an exponential expansion of order Q l is applied from [90]: r l,q e −s l,q p , ( with the coefficients r l,q and s l,q found by minimising a mean square error function.Now, it can be rewritten into the form shown in (35).And from this point, the update equations can be found.
A very similar approach is used in Refs.[7,22,70,91].However, they differ in certain aspects, with one being the way in which the optimum truncation point of a series approximation of the denominator of the complex permittivity is found.For example, Ref. [70] uses a truncated binomial series to approximate the complex permittivity denominator, while Ref. [22] uses a fractional power expansion.Another difference is the method used to calculate the coefficients in (35).While Refs.[7,70] minimise the mean square error function, Ref. [22] uses enhanced weighted quantum particle swarm optimisation (EWQPSO) on a nonlinear minimisation function.
It should be noted that Ref. [22] reduced the computational costs of the scheme in Ref. [70] to O(N 2 ), as the formulation could be written in matrix form, with the inverse of a coefficient matrix only needing to be computed once.
Thus far, the fractional derivative methods have used the Riemann-Liouville derivative due to its power law kernel.However, the Grünwald-Letnikov derivative has been used to model Cole-Cole media [89] and Havriliak-Negami media [92].Following Ref. [92], they define the permittivity function as in (14), plug it into (10), and perform the IFT, thus resulting in the equation from Ref. [92]: where ∆ϵ is the numerator of the fraction in (14), and where F −1 represents the IFT.In a similar way, since there was no closed form solution available, D t is approximated by a sum of time fractional derivatives, like in (34): where q i and r i are optimisation variables found using the Nelder-Mead nonlinear optimisation algorithm [93], thus leading to the following polarisation relation [92]: From this point, the Grünwald-Letnikov derivative is employed, with the r i th fractional derivative of the discreitized P(t) given as in Ref. [92]: where w (r i ) k are coefficients found using a recursive relation.However, (42) requires all previous values of P n to be stored, meaning the memory requirements increase with time.
To reduce this burden on the memory, the authors of Ref. [92] applied Prony's method to approximate w (r i ) k using a weighted sum of truncated decaying exponentials [94]: where k ≥ 1, a i,j < 0, and b i,j > 0. Substituting ( 43) into ( 42) and again into (41) gives where v n i,j are the auxiliary vectors.So, only N • M auxiliary vectors need to be stored in the memory, and the fact that they can be computed recursively means that the memory requirements are heavily reduced.From this, the update equation for the polarisation density can be found.And discretizing Ampere's law and Faraday's law gives the electric field and magnetic field update equations, respectively.It is important to note that there is no mention of the numerical stability of the method.Regardless, the method showed good agreement between the estimated and analytical solutions over a wide frequency range of the reflection coefficient and transfer function.
To evaluate the effectiveness of the fractional representation of complex permittivity functions as outlined in (33) and the corresponding fractional FDTD algorithms, a comparison with the Padé approximation model ( 21) that is used in traditional FDTD approaches, like the ADE and Z-transform methods, was conducted in Ref. [22].This benchmarking revealed that the formulations based on fractional derivatives enable a significantly enhanced numerical accuracy.Moreover, unlike the Padé approximation, they lead to more compact representations of permittivity functions.This results in notable advantages, including accelerated convergence rates, reduced memory usage, and simplified integration into the framework of existing FDTD codes.

Conclusions
This study has presented an in-depth analysis of the electromagnetic modelling of dispersive materials using the FDTD technique, thereby focusing on three prominent methodologies: the ADE, the Z-transform, and fractional derivative-based methods.Each of these techniques provides a unique framework for accurately modelling the complex frequency-dependent behaviour of dispersive materials, which is crucial for the design and analysis of electromagnetic devices and systems.
The ADE method facilitates the incorporation of material dispersion by coupling Maxwell's equations with additional differential equations that describe the material's polarisation response.However, the main drawbacks of the ADE method include its complexity in formulating the auxiliary equations for different types of dispersive materials and the potential for numerical instability under certain simulation conditions.These limitations can lead to increased computational resources and longer simulation times, particularly for materials with arbitrary dispersive properties.
Similarly, the Z-transform approach offers a powerful technique to model dispersive materials by converting the time-domain FDTD equations into the Z domain, thereby allowing for the direct incorporation of material dispersion into the FDTD grid.Despite its effectiveness, the Z-transform method is criticised for its inherent approximation errors and the need for inverse Z-transform to return to the time domain, which can introduce inaccuracies in the simulation results.
In contrast, fractional derivative-based methods emerge as a superior alternative for modelling dispersive materials in FDTD simulations.These methods leverage the concept of fractional calculus to more accurately represent the power law behaviour of material dispersion over a wide frequency range.The benefits of fractional derivative-based methods are manifold.Firstly, they offer a more generalised framework that can be applied to a broad spectrum of dispersive materials without the need for material-specific auxiliary equations.This universality significantly simplifies the modelling process and reduces the setup time for simulations.Secondly, fractional derivative-based methods have been shown to produce more accurate results when compared to traditional integer-order differential equations, especially in capturing the long tail effects of material dispersion.This accuracy is paramount in applications where precise electromagnetic behaviour prediction is critical.Thirdly, these methods exhibit better stability and efficiency in numerical simulations, thereby enabling faster convergence and reduced computational overhead.
In summary, while the ADE and Z-transform methods have played a pivotal role in advancing the electromagnetic modelling of dispersive materials, the fractional derivativebased approach offers a promising pathway forward.Its ability to deliver higher accuracy, efficiency, and stability in FDTD simulations makes it an appealing choice for tackling the complex challenges associated with the electromagnetic response of dispersive materials.Future research should focus on further refining fractional derivative models, exploring their integration with more advanced higher-order FDTD schemes, and extending their application to a wider range of electromagnetic phenomena.
The fractional derivative-based techniques presented herein offer a versatile tool for tackling complex electromagnetic challenges in bioengineering and remote sensing.They enable the study of pulse wave propagation and biological processes, including electroporation within the human body.Additionally, these methods find application in the domain of ground-penetrating radar, thereby facilitating the investigation of electromagnetic scattering by objects, such as utilities and landmines buried in realistic multilayered dispersive soils.