Non-Debye relaxations: The characteristic exponent in the excess wings model

The characteristic (Laplace or L\'evy) exponents uniquely characterize infinitely divisible probability distributions. Although of purely mathematical origin they appear to be uniquely associated with the memory functions present in evolution equations which govern the course of such physical phenomena like non-Debye relaxations or anomalous diffusion. Commonly accepted procedure to mimic memory effects is to make basic equations time smeared, i.e., nonlocal in time. This is modeled either through the convolution of memory functions with those describing relaxation/diffusion or, alternatively, through the time smearing of time derivatives. Intuitive expectations say that such introduced time smearings should be physically equivalent. This leads to the conclusion that both kinds of so far introduced memory functions form a"twin"structure familiar to mathematicians for a long time and known as the Sonine pair. As an illustration of the proposed scheme we consider the excess wings model of non-Debye relaxations, determine its evolution equations and discuss properties of the solutions.


I. INTRODUCTION
Typical example of dielectric relaxation is provided by a dipolar system which approaches the equilibrium being earlier driven out of it by a step or alternating external electric field. The phenomenon is usually described in terms of the relaxation function n(t) which counts dipoles surviving depolarization during the time (0, t) ⊂ (0, ∞) and, if normalized, evolves form n(0+) = 1 to n(∞) = 0. The function n(t) comes out as the solution of macroscopic differential equatioṅ n(t) = −r(t, τ )n(t). (1) The non-negative quantity r(t, τ ) is the transition rate of the system and besides of the time depends on properties characterizing the medium among which a material constant called the relaxation, or characteristic, time τ is the most important. Solution to Eq. (1) is easily got as n(t) = exp[− ∞ 0 r(ξ, τ ) dξ] but it remains of very limited physical utility because the knowledge of r(t, τ ), especially for short and long times t, is far insufficient except of the Debye case for which r(t, τ ) = τ −1 = const. Data provided by the broadband dielectric spectroscopy (encoded in the so-called spectral functions) extrapolated to the full frequency range and next transformed to the time domain, do not help very much -using them to calculate the ratioṅ(t)/n(t) usually leads to cumbersome formulae [15,17,28,29], in addition singular at the origin, which * katarzyna.gorska@ifj.edu.pl † andrzej.horzela@ifj.edu.pl ‡ poganj@pfri.hr makes their experimental verification rather impossible for the time laps close to the origin. These difficulties have prompted efforts to look for mesoscopic description of relaxation phenomena based on dynamical rules being non-local in time and leading to the evolution equations which from the very beginning take into account the memory effects. The simplest way to mimic the memory is to introduce the time smearing which may proceed two-fold: either one smears the left hand side of Eq. (1), i.e., the time derivative inṅ(t) = −r(t, τ )n(t) or rewrites Eq. (1) in the integral form n(t) = 1 − t 0 r(ξ, τ )n(ξ) dξ and uses the smearing r(ξ, τ ) → r(t − ξ, τ ). Keeping the relaxation time τ explicitly separated out from the other material depending parameter's vector p, say, this leads to t 0 k(t − ξ)ṅ(ξ) dξ = −B(τ ; p) n(t) , where B(τ, p) denotes the universal, time independent, transition rate and k(t) stands for the memory kernel responsible for smearing the time derivative. In turn, for r(t − ξ, τ ) = B(τ ; p)M (t − ξ), we get with M (t − ξ) being another memory kernel, a priori not connected to k(t − ξ). Mathematically Eqs. (2) and (3) are both the Volterra type equations [23] which utility goes beyond more popular fractional differential equations introduced in the framework of fractional calculus approach to the relaxation phenomena [19]. If we require physically justified equivalence of Eqs. (2) and (3) then the memory kernels become mutually related and form a coupled (Sonine) pair which appearance and properties we shall discuss a bit later. The integro-differential equation (2) is mathematically very well understood [32]. In fact it is the equation which for special choices of k(t − ξ) reduces to equations with fractional derivatives (more precisely, various types of them) for a long time proposed to investigate the relaxation phenomena. Simultaneously, both Eqs. (2) and (3) have the form of kinetic equations which are the starting point to describe relaxation in the subordination framework [44], [46], [47,Chs. 4.1,4.3] developed as a general scheme within the stochastic processes approach to the relaxation phenomena and anomalous diffusion.
The cornerstone of the stochastic processes based approach to relaxation (as well as to anomalous diffusion if one adopts a suitably reinterpreted language) is the assumption that the transition rate r(t, τ ) introduced in Eq. (1) takes on the meaning of a non-negative stochastic quantity parametrized by the randomized characteristic time τ . The latter does not denote any longer one among material properties of the relaxing medium and becomes physically meaningful variable which shape of postulated randomization strongly influences, or even determines, modeling the relaxation. The choice of stochastic processes proposed to investigate relaxation phenomena is dominated by choosing those which are non-negative, nondecreasing and have distributions which are infinitely divisible. The last means that relavant distributions functions f ( β) are representable as N → ∞ limit of distributions obeyed by random variables β (N ) = 1 N N i=1 λ i where all λ i are independent identically distributed random variables [50]. Randomization of the characteristic time τ (in the Debye systems assumed to be the same and fixed for all dipoles forming the system) means that we are going to change description of the system -instead of looking for deterministic evolution in the time t measured by a laboratory clock we search for stochastic evolution in terms of the "internal" time τ (t) whose dependence on t is hidden in some probability distribution f (τ, t).
Any non-negative stochastic process whose distribution is infinitely divisible, herewith denoted as U (τ ), satisfies the relation where Ψ(s) bears the name of characteristic (either Laplace or Lévy) exponent and is uniquely given by the Lévy-Khintchine formula [41, Eq. (1. 3)] where µ(dx) (subject to some additional conditions) is called the Lévy measure while λ is named the drift parameter. For s > 0 the relation (5) places all functions Ψ(s) in the class of Bernstein functions (BFs), i.e. nonnegative functions on R + , differentiable infinitely many times and satisfying for s > 0 and n ∈ N 0 the conditions (−1) n f (n+1) (s) ≥ 0 everywhere in their domain [51]. We remark that the BFs are close relatives to the completely monotone functions (CMFs), also being nonnegative on R + , differentiable there infinitely many times and satisfying (−1) n f (n) (s) ≥ 0, s > 0, n ∈ N 0 . To make these notions more intuitive one may understand BFs as "maximally regularly" increasing positive functions while CMFs as "maximally regularly" decreasing ones [52]. The deep mutual relation between infinitely divisible distributions, BFs and CMFs is encoded as follows: for h : (0, ∞) → (0, ∞) the following statements are equivalent: (i) h is CMF and it is infinitely divisible with h(0+) ≤ 1; and (ii) h = exp(− Ψ) with Ψ being BF [42,p. 52,Lemma 5.8]. Coming back to the relaxation phenomena we remind that the relaxation function n(t) (which provides us with the information on the number of relaxation centers which did not decay during the time (0, t)) if calculated from the spectroscopic data appears to be CMF for a vast majority of commonly used phenomenological models [6,12,20,21,24,35,49]. This fact merged with the just mentioned theorem strongly suggests that characteristic exponents are inextricably linked with investigation of the relaxation processes. Research which sheds light on this problem is the leitmotif of our paper.
We present and discuss a number of arguments which clarify the role played by characteristic exponents in description of the relaxation phenomena, in particular provide the reader their interpretation as memory functions. The methods which we advocate are general and, as recently demonstrated in [22,48], applicable to various phenomenological models of relaxation. In what follows we focus our attention on the excess wings model [26,27] which goes beyond the Jonscher universal relaxation law (URL) [31] and is less popular among experimentalists if compared with models of the Havriliak-Negami family. General considerations of Sec. II show how the characteristic exponent enters the spectral, relaxation, and memory functions. Also we explain the physical interpretation of Ψ(s) and demonstrate that required equivalence of Eqs. (2) and (3) inevitably leads to the concept of the Sonine pair which non-negligible role in theoretical studies of relaxation, viscoelasticity and anomalous diffusion was recently noticed, analysed and developed [13,14,25]. Starting from Sec. III we investigate the excess wings model. Using its spectral function we recover suitable characteristic exponent whose knowledge enables us to find the appropriate relaxation function. In Sec. IV we use the characteristic exponent to introduce two coupled memory kernel functions which form the Sonine pair. Thus we arrive at a pair of evolution equations which involve either the smearing of the relaxation function or its time derivative and should give the same excess wings relaxation function. Both equations are solved in Sec. V where also requirements demanded from their solutions are checked. The paper is concluded in Sec. VI.

II. CHARACTERISTIC EXPONENTS AS CONSTITUTIVE ELEMENTS OF THE RELAXATION THEORY
As sygnalized in the Introduction the first step in the construction of stochastic approach to relaxation phenomena, see [44] and [46,47] for recent exhaustive reviews, is to assume that they are underpinned by randomization of the characteristic time τ and that the stochastic processes U (τ ) emerging from such a randomization have non-negative infinitely divisible distributions. In the majority of physically meaningful applications these distributions are realized as heavy tailed α-stable Lévy ones related to various variants of random walks. The second step, essential for making the method effective in modeling physical applications, is to use the subordination formalism [4] within which the parent process, usually the Debye law dependent on operational time τ , is subordinated by a directing process τ (t) which links τ and t in a random relation encoded in a probability density (pdf) f (τ, t). Intuitively, employing the subordination scheme means to replace a process described in terms of the laboratory clock measured time t by a composed random process governed by an irregular nondecreasing flow of randomized time τ given by a stochastic process t → τ = τ (t). Physically it is expected that properties of τ (t) may shed light on the internal structure of the system or provide us with some hints how its macroscopic behaviour is influenced by many-body effects.
According to Eq. (4) the probability theory introduces the characteristic exponent Ψ(s) in terms of the mean value of the exponentiated non-negative stochastic process exp (−sU (τ )). Suppose that f (τ, t) (which says how to find the system in the operational time τ if it is in the laboratory time t) is also the infinitely divisible pdf of U (τ ). Then, Assumption that the process U (τ ) results from t → τ : U (τ ) ≤ t with the pdf f (τ, t) opens the possibility to ask for the "inverse" process S(t) = inf{τ : U (τ ) > t} and its pdf g(t, τ ). The latter may be calculated from the cumulant distribution functions of U (τ ) and S(t) (see e.g. [8,44]) Taking the Laplace transform of Eq. (7) and using Eq. (6) leads to [53] g(s, τ ) = Ψ(s) s e −τ Ψ(s) .
Alternatively, any non-Debye relaxation process may be seen as summing up effects of multichannel exponential decays with each channel characterized by some randomly distributed relaxation time θ. Under this assumption the relaxation function n(t) counting the fraction of objects which have survived the decay in the laboratory time interval (0, t) boils down to the weighted average of exponential decays where µ(dθ) denotes the probability with which the random relaxation time θ occurs. In the framework of the subordination approach the same quantity n(t) comes from weighted average of the Debye law expressed in the operational time τ and the pdf g(t, τ ). Thus Eq. (9) may be rewritten as the integral decomposition [11] Using Eq. (10) enables us to calculate the response (called also spectral) function defined in the frequency domain as φ(iω) = L −1 [−ṅ(t); iω]. Because of Eqns. (8) and (10) it is uniquely expressed in terms of the characteristic exponent Here we point out that Eq. (11) explicitly determines the relation between purely phenomenological object which is the spectral function φ obtained as a fit to experimental data and Ψ, a mathematical quantity one to one related to the stochastic process being assumed to underlie physical phenomenon under consideration but of origin rather loosely supported by specific physical properties of the system. To look for physical justification of so far presented construction notice that the relation and Eq. (11) implies As recalled in the Introduction the time evolution equations involving memory effects may be obtained by modeling memory effects through the time smearing, either ofṅ(t) like it has taken place in Eq. (2) or of r(t, τ )n(t) like has been done in Eq.
while for Eq. (3) we get where k(s) = L[k(t); s] and M (s) = L[M (t); s]. Physical equivalence of the above approaches requires that Eqs. (14) and (15) describe the same situation, i.e., the memory effects influencing the behaviour ofṅ(t) and r(t, τ )n(t) should yield the same results for observed properties of n(t). The equality of Eqs. (14) and (15), i.e., n k (s) = n M = n(s), if compared with Eq. (13), gives

III. THE EXCESS WINGS MODEL
The spectral function which corresponds to the simplest version of the excess wings model is It depends on two characteristic times τ 1 > 0 and τ 2 > 0 and so does not fit to the Jonscher's URL [12,31,45] involving only a single characteristic time τ . Despite this reservation the excess wings model appears useful in analysis of experimental data as it successfully describes the relaxation phenomena in the high frequency regime when the frequency of applied electric field is of the order 10 5 − 10 10 Hz [5,9,10,26,27]. For α = 1 the spectral function φ α (iω) is proportional to the Debye spectral function φ D (iω) = [1 + (τ 1 + τ 2 ) iω] −1 with the characteristic time τ 1 + τ 2 , i.e., φ 1 (iω) = φ D (iω) + iωτ 2 φ D (iω). Comparing the spectral function (17) with (11) we find that the characteristic exponent Ψ formally reads if we set B(τ, p) = B(τ 1 , τ 2 , α) = τ α 2 /τ 1 = τ . But some doubt arises: is the construction described in Sec. II legitimate if we have two characteristic times -which of f (s) CMF BF (s + b) µ , b ≥ 0 µ ≤ 0 µ ∈ (0, 1) (s ν + b) µ , b > 0 µ ≤ 0 and ν ∈ (0, 1] µ ∈ (0, 1) and ν ∈ (0, 1)  [6,Theorem]. It leads to the crucially important result -namely enables us to represent in an unique way Ψ(z) as the Laplace transform of a nonnegative function. Furthermore, restricting the argument z of Ψ(z) to the positive semiaxis, i.e., z = s > 0, we can identify Ψ(s)| s∈R+ as a BF and make use of a plethora of results concerning CMFs and BFs. In the first step notice that for s > 0 the function Ψ(s) is non-negative while its first derivative is CMF. Indeed, from Tab. I we see that for non-negative τ 2 and α ∈ (0, 1) this expression is a convex sum of CMFs and hence it is CMF as well. Thus, the characteristic exponent Ψ(s) itself is BF. This is the result which we do need and which for the case under consideration is by no means obvious from the stochastic point of view since we lack the information concerning the infinite divisibility of underlying stochastic process. Needed result, which obviously confirms infinite divisibility, is obtained from the completely different sources, namely from the phenomenology merged with mathematical analysis. We would also like to remark that within the stochastic approach we deal with functions of real variables exemplified by those being CMFs and BFs. Starting from the spectral function treated as a complex function of the complex variable we avoid the path marked out by principles of the stochastic approach. Equipped with tools of the complex analysis we can leave aside the probability rooted description of relaxation phenomena and may understand much better results not once or twice hidden behind paradigms of the real functions approach [21].
The relaxation function n(t) Eq. (10) with substituted Eq. (8) reads where τ = τ α 2 /τ 1 was introduced a few lines above. To get Eq. (19) we changed the order of integration over ξ ∈ [0, ∞) in the inverse Laplace transform which reduced the integral over ξ to the elementary one. The inverse Laplace transform in the lower line of Eq. (19) can be calculated by virtue of the formula [40,p. 10,Eq. (1.38)] from which we get As shown in [34] this function, known as the binomial (multivariable) Mittag-Leffler function, is defined by the double power series ,x, y ∈ R, and is non-negative for λ 1 , λ 2 ≥ 0, β ∈ (0, 1), and α 1 , α 2 ≥ β − 1. Thus, n(t) given by Eq. (20) is non-negative for α ∈ (0, 1). Notice that in Eq. (21) the infinite sum over k is followed by sums over l 1 and l 2 constrained by l 1 + l 2 = k. As a consequence the double sum in l 1 and l 2 can be represented two-fold: (a) l 1 = 0, 1, . . . , k and l 2 = k − l 1 or (b) l 2 = 0, 1, . . . k and l 1 = k − l 2 . Without loss of generality we consider the case (a). In such a case Eq. (20) becomes .
Using the definition of Mittag-Leffler polynomials (A6) we get The same expression as in Eq. (22) will be obtained if we change k≥0 k l1=0 into l1≥0 k≥l1 . Making this change and setting r = k − l 1 we transform the series and the sum sitting inside Eq. (22) into two independent series Treating once the series over l 1 and another time the series over r as the definition (A1) of the three parameter Mittag-Leffler (or Prabhakar) function [54] we express Eq. (24) in two equivalent forms, namely Calculations made for (a) can be repeated for (b) with l 2 written instead of l 1 ; thus, r = k − l 2 . The formulae (25) and (26) , we set µ = α, ν = (1 − α)k + 1, and γ = −k. That allows us to rewrite Eq. (23) in the form (27) where we change the summation index in the first series by setting k + 1 = r. From Eqs. (A3) and (A4) it comes out that Eq. (27) can be expressed as is the fractional derivative in the Riemann-Liouville sense for α ∈ (0, 1) whereas (I ν 0 f )(x) (given by Eq. (A5)) is the Riemann-Liouville fractional integral for ν ∈ (0, 1). We point out that the time operator in square bracket of (28) is equivalent to [30,Eq. (5) for β = 0]. Acting with I 1 0 on both sides of Eq. (28) we get represented also in the form That leads to Eq. (3) with B(τ 1 , τ 2 , α) = τ = τ α 2 /τ 1 and which is interpreted as power-like smearing of r(t, τ )n(t) in Eq. (1). The related Laplace transform becomes From the above and Eq. (16) we restore the characteristic exponent described by Eq. (18).

B. Coupled memories
The explicit form of the characteristic exponent Ψ(s) enables us to find memories M (t) and k(t) responsible for the time smearing of Eq. (1). Recall that the memory M (t) reflects the smearing of n(t) whereas k(t) is related to the smearing of the time derivativeṅ(t) and that the memory M (t) and its Laplace form are given by Eqs. (29) and (30). Using the coupled pair M (s) k(s) = 1/s we find that k(s) and its Laplace form k(t) yield (31) The singularity of M (t) and k(t) at t = 0 is controlled by the parameter α. In the example quoted just below Corollary 4.1 in [25] it is pointed out that M (t) and k(t) are the so-called Sonine functions and the coupled pair (k, M ) is the Sonine pair [25, pp. 213-4]; at a moment we conclude that they are only Sonine functions k and M . Such functions are locally integrable non-decreasing functions which satisfy Thus, [25,Theorem 3.1] is revealed. According to the philosophy of the coupled memories M (t) is linked to Eq. (28) and Hence, the smearing of the relaxation function n(t) can be changed into the smearing of its first time derivativė n(t) like it is done in Eq. (32).

V. THE SERIES FORM OF SOLUTIONS TO (32)
General conditions of solvability Eq.(32) are precised in [32,Theorem 2]. It guarantees the uniqueness of the solution, its continuity, differentiability, and completely monotone character on (0, ∞). From Eq. (31) the asymptotics of k(s) turns out to be so we reconstruct the conditions listed in [32,Theorem 2] except of the first of them: k(s) does not tend to infinity with s → 0 but to the constant τ α 2 instead. This clearly suggests the existence of a solution to (32) which differs from (20).

VI. CONCLUSIONS
We have shown that the kinetic equations (2) and (3) assumed to govern the relaxation phenomena and stemmed from the time smearing of either LHS or RHS in non-Debye evolution equationṅ(t) = −r(t, τ )n(t) determine their stochastic interpretation. The crucial role in the presented approach is played by the characteristic exponent Ψ which provides us with a bridge connecting kinetic equations and stochastic methods. Moreover, for a large set of relaxing systems Ψ obeys well-defined properties which put it in the class of Bernstein functions and open new ways to push forward mathematical and physical understanding of the relaxation phenomena.
To illustrate our methods we went beyond the family of the Havriliak-Negami models and considered the excess wings model of relaxation. We identified the characteristic exponent related to it and derived and solved kinetic equations which reflect two ways of introducing the memory effects -the time smearing ofṅ(t) or r(t, τ )n(t) reflected in Eqns. (2) and (3), respectively. Natural assumption that both approaches lead to the same physical results allowed us to claim that the memory functions, M (t) and k(t), responsible for both variants of smearing, form the Sonine pair, i.e., their transforms to the Laplace domain satisfy M (s) k(s) = 1/s. Results of the paper complete and, in a sense, unify so-called deterministic and stochastic processes based investigations of the non-Debye relaxation phenomena. We show that both these approaches are not only mutually related but realize a correspondence principle which joins different, but in fact equivalent views on the same physical problem.