Relaxation Equations: Fractional Models

The relaxation functions introduced empirically by Debye, Cole-Cole, Cole-Davidson and Havriliak-Negami are, each of them, solutions to their respective kinetic equations. In this work, we propose a generalization of such equations by introducing a fractional differential operator written in terms of the Riemann-Liouville fractional derivative of order $\gamma$, $0<\gamma \leq 1$. In order to solve the generalized equations, the Laplace transform methodology is introduced and the corresponding solutions are then presented, in terms of Mittag-Leffler functions. In the case in which the derivative's order is $\gamma=1$, the traditional relaxation functions are recovered. Finally, we presente some 2D graphs of these function.


Introduction
Fractional calculus has greatly developed during the last years and has established itself as a generalization of classical differential and integral calculus. Over the last three decades the investigation on this subject has intensified and applications were found in many fields of science such as physics, mathematics, chemistry, financial systems, economy and engineering [1]. As examples we may mention the following applications: the fractional control of engineering systems; variational calculus and optimum control for fractional dynamic systems; numerical and analytical tools and techniques; fundamental explorations of mechanical, electric and thermal constitutive relations; investigation of properties of several engineering materials such as viscoelastic polymers, foams, gels and animal tissues; diffusion phenomena; bioengineering and biomedical applications; thermal modeling for engineering systems such as breaks and tool-machines and, finally, imaging and signal processing [2,3,4,5,6,7,8,9,10].
The mathematical investigation of relaxation processes in dielectrics has been conducted with the use of fractional calculus tools [11,12,13,14]. The properties of dielectric materials are usually represented by the empirical susceptibility functions as they have been initially modeled by Debye (D) in 1929, then by the Cole brothers (C-C) and, later, in the works of Cole-Davidson (C-D) and Havriliak-Negami (H-N). With the help of the concept of memory developed in the formalism of Mori-Zwanzig, it is possible to deduce the kinetic equations for the relaxation processes described by those empirical functions [15]. The equations obtained are expressed in terms of the fractional derivative of Riemann-Liouville and their solutions are given in terms of Mittag-Leffler functions [16].
The Mittag-Leffler functions of a real variable t, with one, two and three parameters, are convenient to study anomalous processes such as those occurring in dielectrics, including the aforementioned classical cases C-C, C-D and H-N. It is crucial to emphasize that these functions are completely monotonic for t > 0 [11,17].
In this work we study fractional kinetic equations with a derivative of order γ, 0 < γ ≤ 1, proposed as generalizations of the classic kinetic equations associated with relaxation processes in dielectrics, in such a way that the traditional kinetic equations become particular cases of the fractional equations when γ = 1.
This work is organized as follows: in the first section a brief theoretical introduction is presented of the four empirical expressions of complex susceptibility (in the frequency domain) related to Debye's model and to the anomalous relaxation processes C-C, C-D e H-N. On the basis of these expressions of complex dielectric permittivity, we present in the second section the construction of the kinetic equations for the Debye and the anomalous relaxation functions. The construction of these kinetic equations is carried with the introduction of an integral memory function. At the end of the second section the equations are solved and the relaxation functions (in the time domain) are written in terms of Mittag-Leffler functions; we also show their graphic representations. In the third section, the previously mentioned fractional kinetic equations are solved and the graphs of such solutions are shown. Finally, the fourth section, dedicated to the conclusions, presents a comparative analysis of the results obtained and to the discussion of outline a proposal for future study.

Classic Empirical Models
The most important phenomenon associated with a dielectric material is its polarization, which consists in the change of the distribution of its molecular and atomic charges when it is subjected to the action of an electric field. Thus, when an electric field is applied to a dialectric it produces a very small electric current called dielectric loss, and its constituent particles, ions or molecules suffer small dislocations or rearrangements, thus altering their equilibrium positions [18]. These molecular parts do not leave nor reach their state of equilibrium instantaneously: a variable ammount of time is necessary for this change of positions to take place. This time lapse necessary for the material to respond to the electric field applied is called relaxation time. Thus, when submitted to an electric excitation, the dielectric will respond to this action in an attempt to reestablish its equilibrium during and after the electric stimulus. The polarization of each dielectric material depends on the nature of its molecular and atomic chemical bonds, and there is presently no universal model which can explain the polarization phenomenon in all materials [19]. Debye's response function was the first theoretical model for the dielectric behavior of some substances [20]. However, due to its limitations, this model is incapable of describing in detail the dielectric response of a large number of solids and liquids.
Therefore, in the years following the appearance of Debye's theory, several other response functions were proposed in the literature to serve as models for describing the dielectric relaxation of many materials. Among them, we here focus our attention on the model proposed by the Cole brothers [21,22], and the one by Cole-Davidson [23,24], both of which emerged in attempts to adjust the response function to the experimental behavior of some dielectric materials. There is also the model of H-N [25,26] which can be considered a generalization of the latter models. These last three models, called anomalous models, are the most relevant in the literature; however, there are other models which approach the phenomena from a different perspective, as the models by Weron and collaborators [27,28,29,30], Hilfer [13,31], Hanyga and Seredynska [32].
During the last years, some of these authors and others not mentioned above have addressed the anomalous relaxation processes of types C-C, C-D and H-N with the help of fractional calculus, for example, by associating the response functions to the Mittag-Leffler functions (which turn out to be related with the complex susceptibility s = iω through the Laplace transform) [11]; by modeling relaxation processes through equations with fractional derivatives, whose solutions are also given in terms of Mittag-Leffler functions; and establishing other theoretical connections discovered through the use of fractional calculus tools.
In 1929, Debye conceived a simple model for the relaxation process in which he supposed a unique relaxation time for all molecules, obtaining the following expression [18]: whereε(s) D is the complex susceptibility, ǫ * is the complex permittivity of the dielectric, the real value ǫ 0 is the low frequency dielectric constant, the real value ǫ ∞ is the high frequency dielectric constant and σ is the constant associated to the dipole's characteristic relaxation time.
Cole and collaborators [21,22] formulated a more complete model based on experimental data on dielectrics. This new formulation contains a parameter α which can assume values between 0 and 1 and is given byε In their studies of the dielectric relaxation process, C-D [24] proposed a modified equation for the dielectric liquids glycerol and propylene-glycol. They also studied n-propanol and confirmed that, in its case, the process is described by Debye's equation. The difference from the former two liquids lies in the broadening of the dispersion range under higher frequencies [21].
The works of C-D [23,24] allowed for an almost complete determination of the dielectric properties and a more accurate quantitative description. The empirical expression for the complex susceptibilityε(s) is given by:ε where 0 < β ≤ 1.
In their works, H-N [25,26] studied the complex dielectric behavior of twenty-one polymers and noticed that they had approximately the same form. They then arrived at an empirical expression which generalizes the dispersion models of Debye, C-C and C-D. This is the expression representing the relaxation process:ε It is important to observe that for β = 1 we have the C-C model, while α = 1 takes to C-D model and, finally, with α = β = 1 we recover the first model proposed by Debye.

Kinetic Equations
The properties of dielectric materials are usually described by two constants ε ′ and ε ′′ which are called dielectric constants (or loss factors). They can be combined in a complex dielectric constant given byε This constant, known as complex dielectric permittivity, is given by the following superposition relation [33]:ε where L [f (t)](iω) =f (iω) is the Laplace transform of f (t) in variable iω. We have that ϕ(t) is the normalized polarization decay function when a macroscopic electric field is removed from its medium. Function ϕ(t) contains only the contributions from the relaxation process and we have chosen ϕ(0) = 1 [34].
In the case of linear approximation response, the polarization changes caused by thermal motion are the same as for the macroscopic function dipole relaxation induced by the electric field [35]. Therefore, the laws governing the dipole correlation function φ(t) are directly related to the kinetic properties and macroscopic structures of the dielectric system, represented by function ϕ(t). Thus, it is possible to equate the relaxation function ϕ(t) to the macroscopic dipole correlation function φ(t) as follows: where M (t) is the fluctuating macroscopic dipole moment.
In the projection operator formalism developed by Mori [36] and Zwanzig [37], function φ(t) is called temporal correlation function, as the dipole correlation function defined above is a specific case of the temporal correlation function. Thus, function φ(t), now called temporal correlation function, has the following form in the approximation context [38]: This is an integro-differential equation which takes into account the effects of memory. Therefore, introducing the concept of an integral memory function given by M (t) = t 0 K(x)dx and using the fact that the relaxation function ϕ(t) also satisfies Eq. (8), it is possible to obtain the following relation involving function ϕ(t): where * denotes a convolution product. The relations given by Eqs. (9) and (6) can now be used to calculate the integral memory function M (t). Applying the Laplace transform to Eq.(9) we obtaiñ whereM (iω)) is the Laplace transform of M (t). Substituting Eq.(10) into Eq. (6), we obtain the following relation: Comparing the relation expressed by Eq.(11) with the classical empirical laws (1)-(4) (where s = iω), it is possible to obtain the corresponding memory functions: Applying the inverse Laplace transform to functions (12)- (15) we find that the memory functions in the time variable are respectively given by: with 0 < α ≤ 1, 0 < β ≤ 1 and where E c a,b (·) is the Mittag-Leffler function with three parameters. This function contains as particular cases the Mittag-Leffler with two parameters (c = 1) and one parameter (b = c = 1).
The Mittag-Leffler functions with two and three parameters are defined, respectively, by and where Re(a)> 0, Re(b)> 0, Re(c)> 0 and (c) k = Γ(c+k) Γ(c) is the Pochhammer symbol. Substituting the memory functions given by Eqs.(16)- (19) into Eq.(9), it is posssible to obtain the kinetic equations associated with the models of Debye, C-C, C-D and H-N, respectively: Expression D γ t f (t) denotes the Riemann-Liouville fractional derivative, defined by The solutions of kinetic equations (22)-(25) are given respectively by: The graphic representation of solutions (27)- (30) can be seen in Figures 1, 2, 3 and 4, respectively 1 .
1 In all graphs we have σ = 1.

Fractional Kinetic Equations
As mentioned above, the generalization of the laws describing anomalous relaxation phenomena in dielectrics has been the object of studies by several important researchers. Here we turn our attention to a recent study [39] in which the following initial value problem has been considered: with u 0 and λ positive constants. Parameters α and β are subject to the conditions Fractional equation (31) appears as a generalization of the well-known model by Kohlrausch (Kohlrausch-Willians-Watts) based on the "stretched" exponential function [40,41,42] and its solution is represented in terms of a function introduced by Kilbas and Saigo [43,44], which is given by the series with α, m, l ∈ R, α > 0, m > 0 and α(im + l) = −1, −2, −3, . . . Without loss of generality, we can take λ = 1 in problem (31); the solution is then The conditions given in Eq.(32) guarantee the existence and complete monotonicity of solution u(t) for t ≥ 0. We now mention some particular cases. For β = 0 we have the particular case u(t) = E α,1,0 (−t α ) = E α (−t α ), where E α (·) is the classical Mittag-Leffler function [16]. The second particular case is associated with the wellknown "stretched" exponential function and is given by where α = 1 and −1 < β ≤ 0.
Still supposing β = 0 one finds the solution u(t) = E 1,1,0 (−t) = exp(−t), which is the exponential solution to the relaxation equation. With an analogous generalization purpose, Garra et al. [14], using the general theory of the hyper-Bessel operator [3], solved a generalized relaxation equation which was called fractional differential equation with time-variant coefficient. In that equation, operator t θ d dt α replaces the derivatives in the relaxation equation.
Here, the aim of the paper, we will consider the following fractional differential equation: In this equation, operator D γ t is the Riemann-Liouville fractional derivative, defined in Eq. (26). Imposing the normalization condition D γ−1 t ϕ(0) = 1, it is possible to obtain the fractional kinetic equations with the help of the memory functions explained in the previous section.
First, memory function Eq.(16) is substituted into Eq. (35) in order to obtain: Applying the Laplace transform to it, we obtain Then, isolatingφ(s) we getφ Calculating its inverse Laplace transform, Debye's fractional relaxation function turns out to be For γ = 1 we recover Debye's exponential solution ϕ D (t) = e −t/σ . In the same way, substituting Eq.(17) into Eq.(35) we have Again, applying Laplace transform to Eq.(40) we obtain the following expression: whence emerges the C-C fractional relaxation function: For α = 1 we recover Debye's fractional function Eq. (39). For γ = 1 one has the C-C model given by Eq. (28). In order to generalize the C-D model, Eq.(35) can be substituted into Eq. (18), in order to obtain Applying the Laplace transform to the last expression, we obtain: Expanding the exponential function and the Mittag-Leffler function with two parameters in power series, we get Multiplying and dividing by Γ(β + j + βk) and rearranging factors, this gives: which results in or, isolatingφ(s),φ Thus, using the inverse Laplace transform we finally get Now, taking β = 1 in Eq.(49) and using the following relation involving the Mittag-Leffler function with two parameters: we recover the function of Debye's fractional model, given by Eq. (39). For γ = 1 we recover the C-D function, given by Eq. (29).
Finally, we substitute H-N's memory function, given by Eq. (19), in Eq.(35), we have: (51) Applying the Laplace transform and using the following equality given in [11], we obtain or, equivalently, From this expression it follows that Applying the inverse Laplace transform, we obtain

Conclusion
This study has shown that certain fractional kinetic equations with a parameter 0 < γ ≤ 1, whose solutions are given in terms of Mittag-Leffler functions, can be used as generalized models of the classical relaxation kinetic equations. The classical solutions turned out to correspond to the particular case γ = 1.
In all the cases discussed, the results of the corresponding classical models (response functions) were recovered, as well as the empirical dielectric complex functions (in variable s).
The graphs of the fractional relaxation functions and of the solutions of the fractional kinetic equations were constructed. The analysis of fractional relaxation function graphs and their algebraic expressions given by Eqs. (39), (42), (49) and (56) has shown that the functions's signals and their complete monotonicity depend on the parameters involved. These aspects deserve a more careful analysis and will be the subject of a future work [45].