Non-Debye relaxations: The ups and downs of the stretched exponential vs Mittag-Leffler's matchings

Experimental data collected to provide us with information on the course of dielectric relaxation phenomena are got according to two distinct schemes: one can measure either the time decay of depolarization current or use methods of the broadband dielectric spectroscopy. Both sets of data are usually fitted by time or frequency dependent elementary functions which in turn may be analytically transformed among themselves using the Laplace transform and compared each other. This leads to the question on comparability of results got using just mentioned experimental procedures. If we would like to do that in the time domain we have to go beyond widely accepted Kohlrausch-Williams-Watts approximation and get acquainted with description using the Mittag-Leffler functions. To convince the reader that the latter is not difficult to understand we propose to look at the problem from the point of view of objects sitting in the heart of stochastic processes approach to relaxation. These are the characteristic exponents which are read out from the standard non-Debye frequency dependent patterns. Characteristic functions appear to be expressed in terms of elementary functions which asymptotic analysis is simple. This opens new possibility to compare behavior of functions used to describe non-Debye relaxations. Results of such done comparison are fully confirmed by calculations which use the powerful apparatus of the Mittag-Leffler functions.


INTRODUCTION
Description of physical phenomena which kinetics is influenced by complexity, disorder or randomness often requires a radical departure from theoretical methods established for analogous, but simpler, phenomena discussed in textbooks of general physics. Such a situation is met when we get interested in study of dielectric relaxations and encounter their time behavior different from the commonly expected exponential decay. Depending on the experimental setup empirical investigation of the relaxation phenomena and collecting the data is done by measuring either their time behaviour or frequency characteristics, i.e. the experiment provides us with data in the time or in the frequency domains. Typical example of dielectric relaxation phenomena is provided by a dipolar system which approaches the equilibrium being earlier driven out of it, i.e., polarized, by a step or alternating external electric field. Depolarization is usually described in terms of the relaxation or spectral functions [5]. The first just mentioned quantity, namely the time dependent relaxation function n(t), counts dipoles surviving depolarization during the time (0, t) ⊂ (0, ∞) and evolves form n(0+) = 1 to n(∞) = 0. The frequency dependent spectral functionφ(iω), describing diffractive and absorptive effects, results from the analysis of phenomenological data obtained as a response of the system when it is probed by the harmonic electric field. Defined as the normalized ratio of dielectric permittivities [ε(iω) − ε ∞ ]/[ε 0 − ε ∞ ], where ε ∞ = lim ω→∞ε (iω) and ε 0 = lim ω→0ε (iω), it is complex valued function which analytical properties stem from those obeyed byε(iω). (1) In typical non-Debye relaxation experiments data measured in the time domain are usually fitted using the stretched exponential or, in physicists' community language, the Kohlrausch-Williams-Watts (KWW) function with α > 0. In what follows we will consider only the case α ∈ (0, 1) which preserves the interpretation of the stretched exponential as the continuous sum of Debye exponential decays weighted by a probability distribution which appears to belong to the class of Lévy stable distributions [35,45]. Phenomenological functions usually used to fit the data in the frequency domain, called the standard non-Debye relaxation patterns, are the Cole-Cole (CC), Havriliak and simultaneously can be obtained from φ HN (iω) or φ JW S (iω) for β = 1. Also, if we set α = 1 and β ∈ (0, 1) in the HN model then we get the Cole-Davidson (CD) pattern. Important property of relaxation patterns (3) is that if transformed to the time domain they all lead to the relaxation functions n(t) expressed in terms of functions belonging to the family of the Mittag-Leffler functions (see Appendix ). Thus we arrive at a two-fold way how to analyse the non-Debye relaxation phenomenawe can take into account their modeling either in terms of the KWW function or to choose the Mittag-Leffler functions. None of these approaches is preferred by fundametal theoretical arguments and thus it is understood to treat them as challengers whose usefulness is to be determined by comparison with experimental data. In our opinion results of such comparison lead far beyond its instructive meaning as they may be used to clarify ambiguities coming from difficulties what experiments are facing in asymptotic (short/long time and high/low frequencies) regimes. Thus looking for arguments shedding light on choosing one of the just mentioned different approaches on experimentally observed data is worth attention and more systematic research. We shall begin comparison of the Mittag-Leffler family and KWW matchings with recalling relations between the KWW and the standard Mittag-Leffler function is responsible for the time behavior of the CC model. While the KWW function has been used in modeling physical processes mainly in the context of relaxations (e.g. the Curie-von Schweidler law) the CC pattern is by no means restricted to this class of phenomena. Taken for real argument, t ∈ R + , and α > 0, the CC pattern becomes an example of the generalized Cauchy-Lorentz (GCL) distributions used in numerous fields of basic and applied sciences. To attract attention on the utility of the GCL distributions recall that for α = 2 it found applications in optics long time ago [38,41] and much more recently in quantum mechanics [40] where it describes the so-called Maxwell's fish eye problem. Interesting application, coming from the interface of the basic and applied science, is using the CC pattern in electrochemistry [11,23], bioelectrochemistry [26][27][28][29] and photovoltaics [30]. Effects of distributed, i.e. non-Debye, relaxation processes, inhomogeneities of the system and possible deviations from the Gaussian diffusion spreading lead to non-ideal interfacial behaviour and cause that to model electrochemical response one has to go beyond simple models of electric circuits including capacity, like e.g. the Randles circuit [? ], and to necessity of modifying current-voltage relations introducing into them (sometimes ad hoc) additional time dependent factors given by the KWW or GCL functions [23]. Recent progress in investigations of electrochemical processes taking place in biological systems has shown that some results coming from the fractional calculus may be useful to push forward understanding of non-ideal interfacial capacitance. Working example is that if the so-called constant phase element (CPE) is mounted to replace the standard capacitor in the effective circuit then the differential relation which describes capacitor discharging becomes fractional. Thus we leave the realm of exponential decays (also generalized, like the KWW pattern is) and it becomes quite natural that functions characteristic for fractional calculus, like the Mittag-Leffler ones, come into play and replace exponential-like decay laws. Namely such an approach has been presented in investigations [26][27][28][29] where it has been also noticed that asymptotic properties of the (standard) Mittag-Leffler function for short times agree with those of the KWW function and that the Mittag-Leffler function interpolates between the KWW for short times and the power-like behavior for long times. Among consequences of this property it has been found that the (standard) Mittag-Leffler function appears useful not only in studies of short time effects characteristic for biochemical processes [27,29] but also in analysis of the long time phenomena occurring in the perovskite solar cells [28,30]. Here we want to turn the readers' attention and emphasize that attempts to fit the data by the (standard) Mittag-Leffler function instead of the KWW decay suggest to try the use of other function belonging to the Mittag-Leffler family, especially if one would be interested in search of matchings working beyond the leading order of small t asymptotics.
Let us suppose that we perform an experiment in which we are able to observe the relaxation process taking place in two samples, each prepared exactly in the same way, having exactly the same structure and put into the same experimental conditions. Measurements performed for the first sample are taken in the time domain and provide us with direct information on the time decay of polarization while for the second sample we collect spectroscopy data which are next transformed to the time domain. In such twin-like experimental setup it arises the question about agreement between the KWW function commonly used to fit the time data and the relaxation function(s) obtained using the Laplace transform of spectral functions Eq. (3) where the choice of suitable pattern emerges from the data analysis. Such comparison was first investigated numerically in Refs. [1, 2] for the KWW and HN models as the authors of analysis were unaware of the Laplace transform of the KWW function. Further study of the problem was announced a few years later in [25]. Currently, having in hands new mathematical tools, at that time unknown to the vast majority of physicists, we are going to show how to extend these results using contemporary knowledge, coming from sources far beyond the phenomenology, of the relaxation phenomena. Information expected to help us emerges from dynamics and evolution equations which govern the relaxation processes [10,42,55,58]. However, to get such equations from ab initio microscopic rules without implementing far going simplifications is extremely difficult, if possible at all. Thus some "effective" theoretical approach has to be used -in the case of relaxation processes suitable mathematical tools are provided by approaches rooted in the stochastic processes theory with the crucial role played by methods grown from the concepts of infinitely divisible distributions and subordination. To give very brief explanation -the most important property of nonnegative nondecreasing stochastic processes governed by infinitely divisible probability distributions is that they are uniquely characterized by functions called the characteristic (either Laplace or Lévy) exponents which carry on all information concerning distributions under consideration. This formalism adopted for studies of the relaxation phenomena leads to an unexpected result which merges basic, mathematical in fact, theory and pure phenomenology -characteristic exponents may be uniquely reconstructed from the knowledge of spectral function, i.e. experimentally obtained relaxation patterns. This provides us with a new tool to compare various schemes describing relaxation processes -as we mentioned a few lines earlier our goal is to study similarities and/or dissimilarities of the relaxation descriptions based on the KWW and Mittag-Leffler functions.
The content of our paper goes as follows. Sec. involves preliminaries concerning the characteristic exponents and stochastic approach to relaxations. The spectral and characteristic functions of the KWW model are computed in Sec.. Knowledge of the characteristic functions relevant for the standard non-Debye relaxation patterns and their asymptotics enables us to judge the challenge which of them is the best candidate to approximate the KWW model -we remark that results presented in Sec. are conclusive only for short times. The last section, Sec., summarizes properties of functions belonging to the Mittag-Leffler family, in particular their asymptotics and (fractional) equations which they obey. We collect in one place results of long and often cumbersome calculations in hope that experimentalists will find them useful in analyses of relaxation experiments. The paper is completed by Conclusions section and five appendices directly devoted to mathematical tools used throughout it.

CHARACTERISTIC EXPONENTS AND STOCHASTIC DESCRIPTION -A BRIEF TUTORIAL
Characteristic exponents, Ψ(s)'s, s > 0, appear as basic objects reflecting properties of nonnegative infinitely divisible stochastic processes U (ξ) parametrized by a nonnegative nondecreasing random variable ξ. They are defined by the relation exp (−sU (ξ)) = exp (−ξ Ψ(s)) (4) and given by the Lévy-Khintchine formula [50,Eq. (1. 3)]. Among properties of the characteristic exponents the most essential is that they belong to the class of Bernstein functions (BFs) closely related to the class of completely monotone functions (CMFs), see Appendix . To make the notions of BFs and CMFs more intuitive one may understand BFs as "maximally regularly" increasing positive functions while CMFs as "maximally regularly" decreasing, but still non-negative, ones. Within the subordination approach to relaxation processes the characteristic exponents Ψ(s) are used to construct distribution functions which subordinate the Debye law assumed to depend on irregularly flowing stochastic operational time ξ. Subordination is realized by convoluting such Debye law with some infinitely divisible probability density function (PDF) g(t, ξ) which provides us with the probability density of finding the system at ξ if it is at the instant of time t measured by a laboratory clock. Having this in mind we can write down n(t) as the integral decomposition [4,9,55] where B(τ ) (B in short-hand-notation) denotes the material, time independent, transition rate characterizing the system. The pdf g(t, ξ) may be calculated from the cumulative distribution function of U (ξ) and its "inverse" process S(t) = inf{ξ : U (ξ) > t} and is uniquely determined by the pdf of U (ξ). In Ref. [51] it is shown that if the characteristic exponent Ψ(s) is the completely Bernstein function (CBF), see Appendix , then there exists its associated partner function Φ(s) = s/ Ψ(s) which also is CBF. The pair of Ψ(s) and Φ(s) satisfies the relation Ψ(s) Φ(s) = s which is called the Sonine property, mathematical condition which enables reformulation of integral equations in terms integro-differential equations and vice versa [24]. This unexpected duality has deeply meaningful consequences which have been noticed and discussed elsewhere [17,56]. Here they are only briefly mentioned in Sec. .

SPECTRAL FUNCTION FOR THE KWW PATTERN
According to our best knowledge the analytic expression for the spectral function of KWW model was found in Refs. [31,32] and is not quoted elsewhere. Numerically it was calculated in [2] employing Eq. (1). Calculations presented in Refs. [31,32] lead to representation of the KWW spectral in terms of the Fox H function: where the contour L omits the poles of Γ(1 + αξ) and Γ(−ξ). We remark that according to Eq. (1) in the Laplace space the Fox H function of Eq. (6) equals to s n KW W (s).
The Fox H function is complicated mathematical object difficult to apply in practice. Except of its general properties only a little information useful in calculations can be found in standard compendia dealing with special functions [8,20,48]. Also it is not implemented in the computer algebra systems like Mathematica and Maple. All this makes calculations involving the Fox H function difficult, time consuming and hard to be verified. To avoid this trouble we have found a way how to express φ KW W (iω) in terms of special functions which are analytically and numerically much more tractable. The solution goes as follows: begin with the observation that in any numerical calculation we are always restricted to using rational numbers such that the parameter α in Eq. (7) may be put equal to l/k. In such a case the Fox H in Eq. (7) can be expressed in terms of the Meijer G function G m,n p,q z| (ap) (bq) (see Appendix ). Setting ξ/k = −u in Eq. (7) we rewrite the latter as where we have used the Gauss multiplication formula, i.e. Γ(nz) = (2π) (1−n)/2 n nz−1/2 n−1 i=0 Γ(z + i/n), applied to Γ(1 − lu) and Γ(ku). The upper and lower list of parameters are denoted as ∆(l, 0) and ∆(k, 0), respectively, where ∆(n, a) is a sequence of numbers a/n, (a − 1)/n, . . . , (a + n − 1)/n. For l ≤ k, the Meijer G function in (8) can be expressed as the finite sum of the generalized hypergeometric functions (see Appendix ) by using [48,Eq. (8.2 For passing between Eqs. (9) and (10) we have used the series definition of generalized hypergeometric function p F q and the formula in which the sum of a r is split into k sums with the term a kr , a kr+1 , . . . a kr+k−1 , namely

COMPARISON OF CHARACTERISTIC EXPONENTS
According to [53,55] the characteristic (Laplace or Lévy) exponent may be retrieved from the knowledge the spectral function φ(s): In the remaining part of the paper we assume B(τ ) = 1 such that all dependence from τ is shifted to the spectral functions φ(s). The spectral functions for the CC, HN and JWS relaxation models are listed in Eq.
(3) while for the KWW model it is given by Eq. (8) and/or Eq. (9). Asymptotic behavior of all considered spectral functions for small and large frequencies confirms the experimentally established Jonscher's universal relaxation law (Jonscher's URL) [36]. From Refs. [31,32] where the parameters α and β belong to the range (0, 1]. We put reader's attention that contrary to the CC model the asymptotics of the HN and JWS spectral functions is governed by two different exponentials -for the CC model it is only α while for the HN and JWS models α and αβ = γ ≤ α ≤. This suggests to consider the CC and HN/JWS cases separately. To compare the characteristic exponents relevant for the above presented models we choose the KWW spectral function as the reference. In Refs. [31,32] it was shown that (13) which flows out also from the series form of φ KW W (s) given by Eq. (10).
(a) The characteristic exponents for CC relaxation model are given by the power-law functions which differ from the asymptotic behavior of Ψ KW W (s) and Φ KW W (s) For large s the leading asymptotic term of Ψ HN (s) is (sτ ) αβ . For small s te relevant asymptotics is got if we rewrite Ψ HN (s) as the series r≥0 Γ(1 + β)(sτ ) αr /[r!Γ(1 + β − r)] − 1 whose first two terms (i.e. the terms with r = 0 and r = 1) give the asymptotics of Ψ HN (s) proportional to β(sτ ) α . Gathered together the asymptotic behavior of Ψ HN (s) reads Ψ A;HN (s) ∼ (sτ ) αβ , sτ 1, and which determine the asymptotics of Φ HN (s) As in the previous example also here Ψ A;HN (s) and Φ A;HN (s) are CBFs for α, β ∈ (0, 1]. The power-law asymptotics given by Eq. (17) for sτ 1 shows that in order to match Ψ KW W (s) the relations of exponentials α HN β HN = α KW W has to be satisfied. It means that α HN may be chosen arbitrarily if simultaneously β HN = α KW W /α HN . Thus the small s asymptotics of Ψ A;HN (s) becomes incompatible with the asymptotics of Ψ A;KW W (s) and matches it only for β = 1 which is the condition reducing the HN pattern to the CC one. Fig.  2, with α KW W = 1/3, α HN = 5/6, and β HN = 2/5, shows that for large τ s Ψ HN (s) and Φ HN (s) fit well Ψ KW W (s) and Φ KW W (s), respectively, but the matching breaks down for small sτ .
Φ A;JW S (s) ∼ s 1−αβ , sτ 1; As in the previous case also here the asymptotics' presented by Eq. (20) are given by CBFs for α, β ∈ (0, 1]. The comparison between Ψ KW W (s) and Ψ JW S (s) as well as between Φ KW W (s) and Φ JW S (s) is shown in Fig. 3. It is seen that for large s Ψ JW S (s) and Φ JW S (s) match Ψ KW W (s) and Φ KW W (s) faster than for the CC and HN models. Nevertheless the matchings for small s remain disappointing although at the first glance they seem to be more acceptable than those resulting from the CC and HN models. This, however, may be treated as an artefact coming from the choice of parameters.
We complete this section with two remarks: 1. The leading order of large s, i.e. short t, asymptotics of all relaxation patterns being considered matches the KWW function.  According to it, for G(s) and q(s) being analytic functions, we have in which one immediately recognizes the structure of integral decomposition. In the probabilistic language, if h 1 (x, ξ) and h 2 (ξ, t) are independent probability distributions, Eq. (22) expresses the Bayes theorem and thus may be treated as a joint probability distributions. Namely this identification is made when stochastic methods are applied to the relaxation theory [9]. Within this approach the non-negative random variable ξ is interpreted as an "internal" or operational time which governs the evolution of the function h 1 (x, ξ) = L −1 [ h 1 (x, s); ξ] bearing the name of the parent process. The second component of the integral decomposition Eq. (22) h 2 (ξ, t) = L −1 [ G(s) e −ξ q(s) ; t] describes the mutual dependence of operational ξ and physical t times. Unlike regularly clocked physical time t the internal time ξ has the nature of a càdlàg (left continuous right limited) nonnegative and non-decreasing stochastic process. According to classification proposed in Ref. [52] and recently reconsidered in Refs. [6,19] both functions h i 's, i = 1, 2, can be either the "safe" or "dangerous" probability densities (PDFs). Sufficient condition to be the "safe" PDF is infinite divisibility, if it is not the case we may deal with an example of "dangerous" PDF. The working criterion to distinguish the "safe" and "dangerous" cases is their adherence to the class of Bernstein function. This guarantees that the features characterizing the "safe" PDFs, i.e. nonnegativity and infinite divisibility, are satisfied (Appendix ). Remember also that h 2 (ξ, t) should be normalized. For h 2 (ξ, t) being the "safe" PDF we can say that it subordinates h 1 (x, ξ). In the opposite case, i.e. for h 2 (ξ, t) being the "dangerous" PDF, we name h 2 (ξ, t) and h 1 (x, ξ) the constituents of integral decomposition only.

Integral decompositions as subordinations
In the case of relaxation theory we know from Eq. (5) that h 1 (x, ξ) is independent on x. It is equal to the Debye relaxation function, i.e., h 1 (ξ) ≡ n D (ξ) = exp(−Bξ). The latter is CBF as it is the nonnegative, normalized, and infinitely divisible with respect to ξ. Thus, n D (ξ) is the PDF of parent process. The function h 2 (ξ, t) ≡ h 2,Ψ (ξ, t) involves G(s) = Ψ(s)/s as well as q(s) = Ψ(s) and it is equal to L −1 {[ Ψ(s)/s] exp[−ξ Ψ(s)]; t}. Normalization of h 2,Ψ (ξ, t) in ξ is fulfilled automatically. Because Ψ(s) is the characteristic exponent, i.e. it belongs to the class of Bernstein functions, then using Appendix we can show that h 2,Ψ (ξ, t) is "safe" PDF and it subordinates the Debye relaxation process h 1 (ξ). With the help of Efross theorem Eq. (22) n(t) can be written as  5)), we know that h 2,Ψ (ξ, t) and h 2,Φ (ξ, t) are always the "safe" PDFs. Thus they may be used to subordinate the Debye relaxation as well as other "safe" distributions, e.g., the normal distribution, cf. Ref. [56]. We can also use another subordinators, e.g. L −1 {s αγ−β exp[−us α ]; t}, which subordinates the CD relaxation model. The example of using two kinds of subordinators is presented in Ref. [19].

Subordinations as signposts leading to evolution equations
The subordination approach allows one to find the evolution equations which govern the behavior of subordinated parent process [52]; this we have named the relaxation function n(t). To obtain suitable equations we need two building blocks. The first of them is the standard evolution equation for n D (t), considered in the Laplace space with s replaced by Ψ(s). The second one is the relation n D ( Ψ(s)) = s n(s)/ Ψ(s) derived from Eq. (23).
The evolution equation (with respect to the time t) of n D (t) is well-known. It readsṅ D (t) = −B n D (t) and in the Laplace space equals to s n D (s) − n D (0) = −B n D (s) (25) with the initial condition n D (0) = 1. After replacing s by Ψ(s) we have Setting n D ( Ψ(s)) = s n(s)/ Ψ(s) and substituting it into Eq. (26) we get which can be rewritten as Taking the inverse Laplace transform of Eqs. (27) and (28) in which we use the Laplace convolution we get, respectively, (29) where and The functions M (t) and k(t) are nonnegative and interpreted as memory functions or kernels. Because Ψ(s) is CBF then M (s) is the Stieltjes function (SF) and s/ Ψ(s) is also CBF, see Appendix . Denoting the latter as Φ(s) we learn that k(s) is also SF [15,17,18]. Using the relation Ψ(s) Φ(s) = s shows that M (s) k(s) = s −1 which means that the memory functions M (t) and k(t) fulfill the Sonine equation Eqs. (29) express the time smeared evolution, either of the relaxation function or its derivative. Detailed discussion of mutual relation between these equations has been presented in [17,18]. From the physical point of view the first equation present in the pair (29) is known as the master equation and may be considered as general modeling of the memory dependent linear evolution scheme. Mathematically Eqs. (29) are both the Volterra type equations [22] which shape and utility goes beyond much more popular fractional differential equations introduced in the framework of fractional calculus approach to the relaxation phenomena [10,15].

Examples
As noticed, Eqs. (23) and (24) yield the same results as governed by equivalent Eqs. (29). Thus, the relaxation function n(t) can be derived with the help either of Eqs. (23) or (24). Without loss of generality we take Eq. (23), the characteristic function Ψ(s) presented in Sec. ?? and more convenient for our purposes this one of Eqs. (29) which describes the evolution of n(t).
(iii) Eq. (23) for the JWS model gives (46) which after using once again the Efross theorem where G(s) = s αβ−1 and q(s) = s α + τ −α can be represented as Calculating the integral over ξ (as it was done in the example (ii)) we can simplify the first inverse Laplace transform. It enables us to write down the above equation in the form To show that L −1 [s αβ−1 e −us α ; t] = Γ(β)t/(αu β )g β α,1 (u, t) we employ Eq. (33). Next, using Eq. (43) we express n JW S (t) as ; t] du, (48) which means that to obtain the relaxation function of the JWS model we can use two kinds of subordination approaches. Using the different subordinators we can subordinate the Debye or CD process. Thus, the processes which lead to the JWS relaxation model can be obtained using the various approaches. The asymptotic behaviour of n JW S (t) can be given by , t τ, and Its evolution equation reads Appendix contains the definition of the pseudo-operator (D α + τ −α ) β . Eq. (50) has been discussed in [54,Eq. (4.3)] and justified by under-and overshooting subordination technique applied for the anomalous diffusion.
Analyzing the above examples we see that the integral decomposition Eq. (22) can be interpreted as an alternative form of Eq. (1) obtained from Eq. (5) by employing Eq. (11). It provides us also missing long time asymptotics of the KWW relaxation model Eq.
(2). The relevant asymptotic behavior for short and long times is obtained in Refs. [31,32] .
Looking at the asymptotics of CC and JWS relaxation functions we see that for short times they have similar power-law behavior as n KW W (t) albeit they differ by a constant. The HN relaxation function differs more significantly -the exponential involves both parameters and to agree the asymptotics we have to put restrictive condition α HN β HN = α KW W . Analogous, but reverse, situation we met when want to agree the asymptotics for large t. This leads to the conclusion going beyond the short time asymptotics one has to be careful with choosing one of the Mittag-Leffler functions as an object suitable to replace the KWW function -to make the proper choice it is necessary to have information concerning the middle and long time behavior of the relaxation function.

CONCLUSIONS
It is known that the KWW function does not describe properly many relaxation phenomena, much better (and still friendly in use) is to use the CC model. But also this model breaks down in many physically interesting casethis was the reason of introducing the HN and JWS models which as phenomenological schemes fitting the data in the frequency domain. If transformed to the time domain both these models lead to the time decay laws given in terms of multiparameter Mittag-Leffler functions, long time unfamiliar to the physicists' community. Simultaneously to get information on the time decay of dielectric polarization is often much more needed for practical applications than data obtained the spectroscopy experiments although the latter are more precise and cover much larger range of the frequency involving 10 or even more orders of magnitude. Unfortunately, a fear of using unpopular and scary looking special functions of the Mittag-Leffler family effectively discourages a vast majority of physicists (first of all the experimentalists) and causes that they consider fitting the data by the stretched exponential function not as a routine coming from the long-time habit but as a method which is the only doable procedure. We consider this situation perplexing and propose to give up this long-time habit. Having in mind that functions describing the relaxation phenomena are widely unknown we propose to begin with analysis of the characteristic exponents expressed in terms of easily calculable functions. Properties of these functions illustrate the problems which we face when compare various theoretical schemes used in the relaxation theory, in particular the problem of choosing the most suitable relaxation pattern. As a benchmark for comparing different relaxation patterns we took the stretched exponential and collated it, one by one, with well established models: Cole-Cole, Havriliak-Negami and Jurlewicz-Weron-Stanislavsky. These models are significantly different and overlap only in the asymptotic regime of large s, i.e., short times. Nevertheless we do not consider this result valueless -it gives a warning that to choose properly we have to look for additional information, e.g. the higher order short time asymptotics of relaxation functions or their long-time behavior. Tools to be used in order to push forward such investigations are collected and listed in the Sec. . Thus our main conclusion is that both theoretical studies as well as the time domain measurements of polarization decays should go beyond the stretched exponential fit and should use more extensively the results obtained by spectroscopy methods translated to the time domain. Here we would like to emphasize that Mittag-Leffler functions are well manageable with standard computer mathematics packages, like Mathematica, Matlab and Maple, and it is not a great problem to familiarize with them. The Lévy stable distributions take distinguished place in our considerations [? ] because they enter the integral representations of basic functions used to describe non-Debye relaxations, namely the stretched exponential exp(−ap α ) [46] and the one-parameter Mittag-Leffler function E α (−ap α ) [47,58], both with a > 0 and α ∈ (0, 1]. Introducing for x, y > 0 modified functions g α (x, y) = x −1/α g α (y x −1/α ) we rewrite the standard Pollard definitions [46,47] as Equations (52) have the form of the Laplace integrals but the variable u enters them in different ways: either through g α (a, u) = a −1/α g α (ua −1/α ) or through g α (u, p) = u −1/α g α (pu −1/α ). Worthy to note is also the different shape of differential relations held for the stretched exponential and the one-parameter Mittag-Leffler function d dp e −ap α = −aαp α−1 e −ap α and The first equality above is a differential relation which introduces the so-called Weibull distribution αp α−1 exp(−ap α ) [57] while in the second equality we deal with an integro-differential relation of the eigenequation form in which the operator c D α p denotes the Caputo fractional derivative, see Appendix . The one-parameter Mittag-Leffler function constitutes a generalization of exponential function as its series expansion differs from the usual exponential function series only by a parameter α present in the argument of the Γ function settled in denominator. Moreover, the α parameter indicates that the eigenequation in Eq. (53) is a generalization of the differential equation for the exponential function obtained when we put α = 1 in any of the equations Eq. (53). The one-parameter Mittag-Leffler function provides us with the representative of a class of functions which generalize the exponential function and are widely met in studies of anomalous kinetic phenomena. Another frequently used function belonging to this class is the threeparameter Mittag-Leffler function E γ α,β (x) whose series form is , α, β, γ > 0 while its integral representation reads E γ α,β (−ap α ) = ∞ 0 e −ua p αu g γ α,β (u, p) du, g γ α,β (u, p) = u −1/α g γ α,β (pu −1/α ), (55) with a generalization of the one-sided Lévy stable distribution g α (u, p), denoted by g γ α,β (u, p), being used [17]. For rational α = l/k ∈ (0, 1) the function g γ α,β (u, p) can be expressed through the Meijer G function (see [17, Corollary 3] and Appendix ) .
For A j = B j = 1 the Fox H functions reduces to the Meijer G function: where conditions listed in Eq. (59) become conditions for [a p , 1] = (a p ) = a 1 , · · · , a p and [b q , 1] = (b q ) = b 1 , · · · , b q . For a full description of the integration contour L H and L G and its properties as well as special cases for the H and G functions, see [48,Secs. 8.2 and 8.3].
The generalized hypergeometric function p F q is defined as follows where (a) r is the Pochhammer symbol (rising factor) equals to Γ(a + r)/Γ(a) = a(a + 1) · · · (a + r − 1).

The completely monotone and completely Bernstein functions
The Bernstein functions (BFs) are non-negative functions on R + , differentiable there infinitely many times and satisfying for s ∈ R + and n ∈ N 0 the conditions (−1) n h (n+1) (s) ≥ 0 everywhere in their domain.
The function h(s), s ∈ R + , is a completely Bernstein function (CBF) if it is BF and h(s) and s/ h(s) have the representation given by the Stieltjes transform [51]. Alternative criterion says that h(s) is CBF if s/ h(s) is CBF [51].
The completely monotone functions (CMF) H(s) are non-negative function of a non-negative argument whose all derivatives exist and alternate on R + , i.e.
The Bernstein theorem uniquely and mutually connects CMF with the non-negative function defined on R + by the Laplace transform: where H(s) is CMF [37,45,59].
We can say that f the Stieltjes function (SF) if, and only if, 1/f is a CBF [51].

Relation between the infinitely divisible distribution and the Bernstein-class functions
The relation between the CBF and infinitely divisible function is expressed by [51,Lemma 9.2]. It say that the measure g on [0, ∞) is infinitely divisible iff L[g; λ] = exp[−f (λ)] where f is CBF.