The Havriliak-Negami and Jurlewicz-Weron-Stanislavsky relaxation models revisited: memory functions based study

We provide a review of theoretical results concerning the Havriliak-Negami (HN) and the Jurlewicz-Weron-Stanislavsky (JWS) dielectric relaxation models. We derive explicit forms of functions characterizing relaxation phenomena in the time domain - the relaxation, response and probability distribution functions. We also explain how to construct and solve relevant evolution equations within these models. These equations are usually solved by using the Schwinger parametrization and the integral transforms. Instead, in this work we replace it by the powerful Efros theorem. That allows one to relate physically admissible solutions to the memory-dependent evolution equations with phenomenologically known spectral functions and, from the other side, with the subordination mechanism emerging from a stochastic analysis of processes underpinning considered relaxation phenomena. Our approach is based on a systematic analysis of the memory-dependent evolution equations. It exploits methods of integral transforms, operational calculus and special functions theory with the completely monotone and Bernstein functions. Merging analytic and stochastic methods enables us to give a complete classification of the standard functions used to describe the large class of the relaxation phenomena and to explain their properties.


I. INTRODUCTION
Researchers working in the fields of natural and engineering sciences are used to understand relaxation phenomena as processes which describe effects related to the delay between the application of an external stress to a system and its response. Phenomena which exhibit such behavior are observed practically everywhere in our surroundings: decay of induced dielectric polarization or magnetization usually does not follow switching off the external electromagnetic field instantaneously and various luminescence phenomena persist in the absence of illumination. Similarly, deformed viscoelastic materials return to their original shapes with shorter or longer delay, kinetics of chemical reactions in many cases is retarded with respect to changes of external conditions, and so on. The simplest mathematical model used to describe relaxation is the exponential, or the Debye [133], law which provides us with the time decay rule proportional to the exponential function exp (−λt), λ > 0. The Debye law works well for different physical systems but suffers more and more difficulties or even fails with the growing complexity of the system under investigation. The first more systematic observations of non-exponential relaxations and/or decays were performed more than one hundred seventy years ago, in the middle of the XIX-th century, by R. Kohlrausch [73] who carried out a series of relaxation experiments. At first he was focused on the mechanical creep but soon after he became interested in electrical phenomena and successfully measured the residual discharge current in the Leyden jar. Studying his own experimental results Kohlrausch found that the data were much better described by the stretched exponent exp (−λt α ), α ∈ (0, 1), which means that the relaxation is slower than predicted by the Debye law. During the next few decades, deviations from the Debye law were observed also in relaxation phenomena going far beyond mechanics of viscoelastic media and effects occurring in simple resistor-capacitor circuits. We mention the luminescence phenomena [10,11] or the response of dielectric material to the step input of a direct current voltage as well. Mathematical models tested to explain experimental results were, as a rule, based on the inverse power-like decay model. Examples include the Curievon Schweidler law of current decay or Nutting's equation for the stress-strain relation, both usually treated as purely phenomenological input proposed without any further justification. This phenomenology dominated period of relaxation studies ended in the context of physics of dielectrics around the 60s and 70s of the previous century. With the progress of experimental techniques, new materials were more and more common applied and investigated. Seminal results were obtained by A. K. Jonscher via analysis of a vast majority of available dielectric relaxation data. He found regularities hidden behind them [64][65][66]. Jonscher's discoveries, unified in the so-called Universal Relaxation Law (URL), boil down to the statement that the asymptotics of relaxation functions, both in the frequency and the time domains, are governed by fractional power-like functions. Also at that time, more than a hundred years after Kohlrausch's experiments, the stretched exponential function got back in the game and found wider interest among physicists thanks to the work of G. Williams and D. C. Watts who applied the stretched exponential to fit the dielectric relaxation data in the time domain [128]. It should be mentioned that despite its popularity among experimentalists, the use of the stretched exponent (now called also the Kohlrausch-Williams-Watts (KWW) function) has remained formal as it has never found a deeper theoretical explanation. Moreover, its use leads to an embarrassing situation: the KWW function, if transformed from the time to the frequency domain, is expressed by rather involved special functions, see Sect. IV. The latter are defined either by special contour integrals in the complex plane or through slowly convergent power series. This makes it difficult to read out correctly the dependence of dielectric permittivities on the frequency of applied harmonic electric field, even if these relations, nick-named spectral functions, are not only well-known phenomenologically from the broadband dielectric spectroscopy experiments, but do appear to be simple rational functions [18]. For that reason, the KWW model is often replaced by the Cole-Cole (CC) [18,30,61], Cole-Davidson (CD) [18,30,61], Havriliak-Negami (HN) [18,30,61], and Jurlewicz-Weron-Stanislavsky (JWS) [30,117,125] patterns. Their spectral functions are much simpler than the KWW model but the opposite situation appears for the time dependence functions.
Contemporary ongoing experiments aimed at collecting results essential for a better understanding of dielectric relaxation phenomena may be classified in a twofold way. The data which are measured in a suitable experimental setup concern either the time dependence of physically meaningful properties characterizing the system (e.g., the number of previously excited dipoles which survive some amount of time after sudden switching off the polarizing field or retreating mechanical deformation of the sample), or the frequencydependent response of the system perturbed by a steplike or alternating (usually harmonic) external electric field which plays the role of time-dependent external stress. In the first case, the results are traditionally encoded in the time-dependent relaxation n(t) (or the response φ(t) = − dn(t)/ dt) function, while in the second case, customarily used objects are complex-valued spectral functions φ(iω). In dielectric physics, the spectral functions φ(iω) mimic normalized complex dielectric permittivity [ ε ⋆ (ω)−ε ∞ ]/[ε 0 −ε ∞ ], whose real and imaginary parts are responsible for diffractive and absorptive properties of the medium, and the material constants ε 0 and ε ∞ denote the low (static) and high-frequency values of the dielectric permittivity. The time and frequency approaches are mutually connected by the Laplace transform (I.2) Equation (I.1) is crucial for our further consideration. It gives an operational rule which constitutes a bridge between the time and frequency descriptions of the relaxation phenomena and opens possibilities to compare them. This applies also if the data come from independent measurements performed in the time and frequency domains on the same, or twin-like prepared, sample.
We shall now give more details on the phenomenological models specified by CC, CD, HN, and JWS earlier.
In the frequency domain the phenomenological pattern relevant to the Debye model is equal to φ D (iω) = (1 + iωτ ) −1 (I. 3) and non-Debye relaxations are modelized through: (I. 4) Real numbers α, β ∈ (0, 1] are called the width and the symmetry parameters, respectively, and each time they are fitted to the experimental data together with the material-dependent characteristic time scale τ . It is easy to see that Eqs. (I.4) extend the Debye rule and that among them the HN and JWS patterns are the most general involving the CC and CD as their special cases [18,30]. The Debye and CC models emerge either from φ HN (iω) and φ JWS (iω) for α = β = 1 or for β = 1 and α ∈ (0, 1), respectively. The CD model comes out from the HN pattern for α = 1 and β ∈ (0, 1). On the other hand, measurements done in the time domain provide us with the data for the relaxation function which are usually fitted by the standard KWW function where parameters τ ′ and α ′ ∈ plane, observed in dielectric relaxation of polymers. Shortly after Havriliak and Negami's discovery, the HN pattern gained popularity because, although involving only three parameters, it appeared to be so flexible and universal that within a few years it became a standard model used to describe data obtained for a large range of relaxation and viscoelastic phenomena. Consequently, the calculation of its time domain counterpart, leading either to the response or to the relaxation function, became a serious challenger to routine parameterizations of their time evolution, based on the KWW functions or finite sums of exponential decays. This has motivated many authors to investigate the time dependence of the HN model and explore it using different mathematical tools. They include the analytical methods like the integral transforms theory and fractional calculus as well as an the approach worked out in the probability and stochastic processes theory and known as subordination methods. We shall show that these seemingly distant approaches merge entirely into the framework of memory-dependent kinetics.
The general expressions for the time-dependent counterparts of Eqs. (I.4) (as well as the frequency-dependent analogue of Eq. (I.5)) were found analytically 20 years ago by R. Hilfer [61,62] who, for arbitrary real values of α and β, expressed functions relevant for dielectric relaxation in terms of the Fox H-functions. Note that the applicability of this class of functions was used 10 years ealier for studies of other relaxation phenomena, especially in the research devoted to rheological models of viscoelastic materials [32,33,84,85,93,112,124]. In applications essential for the relaxation phenomena, the Fox H-functions are replaced by much better-known functions of the Mittag-Leffler family [30] and the Meijer G-function. In fact, its well-known representative, namely the standard Mittag-Leffler function, naturally enters the time relaxation for the CC model. Following and exploring this direction of theoretical investigations we shall show how to express the functions relevant to the HN and JWS relaxation models, like the response or the relaxation functions and the probability densities, in terms of the Meijer G-function and the generalized Mittag-Leffler functions. These families of functions are well-implemented in modern computer algebra systems (CAS) and may be effectively studied using methods of classical analysis much more accessible to a larger community. The crucial step to achieve this purpose is to replace the irrational values of the parameters sitting in the Fox H-functions with rational numbers. This simple, but physically well-justified [135] approximation, enables us, as the first step, to represent results in terms of the Meijer G-functions, and in the second step in terms of the generalized hypergeometric functions which specialize to the generalized, or multi-parameter, Mittag-Leffler functions.

II. THE TIME DOMAIN DESCRIPTION: GENERAL PROPERTIES OF THE RELAXATION FUNCTIONS
Assume that the object of our investigation is the timedependent quantity J (t), like the relaxation or the response function is, which reflects the time behavior of some, suitably chosen, property P(t) of the system. For example, P(t) can be the number of its constituents which at an instant of time t are still excited (have survived the decay), or the depolarization current, or the intensity of luminescence light, etc. The normalized form of J (t/τ ), with the characteristic time scale τ put in, reads: where P 0 and P ∞ denote the values of P(t) for t = 0 and t = ∞, respectively. Parameterization given by Eq. (II.1) implies by construction that J (0) = 1 and J (∞) = 0. We also assume that J (t) is experimentally measurable (at least in the sense of Gedankenexperiment) and that in the system, there exist no internal factors which can revert the decay, i.e., J (t) never increases. Such requirements are too vague to specify J (t/τ ) satisfactorily and additional assumptions are needed to make any prospective mathematical modeling of J (t/τ ). Two commonly used suggestions are: which generalizes the standard exponential Debye decay law, dn(t)/ dt = −λn(t), by extending it to the case when the relaxation/decay rate λ becomes the time-dependent quantity, λ → r(t) instead of being a constant.
is interpreted as a continuous weighted sum of the Debye decays, provided the additional condition that p(x) is the probability distribution function (PDF), i.e., it is nonnegative and normalizable. If we focus attention on the relaxation phenomena and identify J (t/τ ) either with the relaxation n(t), or the response function φ(t), then both these functions may be represented as the Laplace integrals [136] If we have enough information on the time dependence of n(t) or φ(t), e.g., know that their time decay is governed by some KWW-like function, then Eqs. (II.5) and (II.6) may reduce to the well-defined problems of classical mathematical analysis whose solutions allow us to determine the function g(ξ), see e.g. [98]. The above propositions are strongly supported by simple physical intuition and, at first glance, look very attractive as a way to solve the problem. Indeed, imagine that the constituents of the relaxing system are under consideration, e.g., polarized dipoles, decay according to the Debye law but each constituent of the system does that with a different decay rate λ k ; we treat the latter quantity as a non-negative discrete or continuous random variable whose distribution g(λ) we somewhat know. We also assume that the randomization of elementary decays does not influence the Debye pattern dominating the process. Under such formulated conditions of independence the use of a joint probability formula for the fraction of constituents which survive the decay gives i.e., we have arrived at an expression which has an obtrusive interpretation to be treated as a weighted sum of exponential decays. Analogous reasoning is the basis of the majority of routine analyses whose goal is to represent decay curves describing experimental data as a weighted finite (or infinite) sum of exponentials [6,15,27,69,91]. However, despite their attractive simplicity and possible utility, the statements cited above are too naive. Equations (II.3) and (II.4) suffer from essential physical shortcomings and lack of mathematical precision. Equation (II.3), if treated as a source of experimental information concerning r(t), is difficult to be verified both for short and long times. Another objection is that r(t) calculated from exactly solvable fractional models of relaxation becomes singular for t → 0 and thus unphysical, e.g., [41]. The main criticism against Eq. (II.4) comes from the objection that the interpretation of factors exp (−ξt) as the Debye decays, each characterized by its own decay rate ξ, is formal and difficult to be explained. Does it mean that decay rates appearing in such a way are distributed according to some mysterious function p(τ ξ)? What is the origin of the time scale τ introduced ad hoc because of dimensional reasons but finally acquiring universal meaning? Nevertheless, so far proposed Eqs. (II.3) and (II.4), if understood in mathematically and physically correct ways, provide us with useful guideposts which show how to choose, arrange and push forward theoretical methods whose development will lead to a better understanding of relaxation phenomena. The mathematical condition which guarantees that Eq. (II.4) makes sense, comes from classical mathematical analysis -according to the Bernstein theorem (for a brief tutorial on it, on completely monotone (CM) and on Bernstein (B) functions see Sect. IV A) the sufficient and necessary condition to represent the function f (t) through the Laplace integral of the measure g(x) is that f (t) is CM. In Eq. (II.8) we deal with t ∈ R + . It means that we work in the realm of real-valued Laplace integral whose continuation to the complexvalued Laplace transform needs care. Merging Eqs. (II.3) and (II.4) and requiring complete monotonicity of J (t/τ ) represented by Eq. (II.2), one learns that the minus exponent, which is in its right-hand side (RHS), i.e., t 0 r(t/τ ) dt, must belong to the class of Bernstein functions, t 0 r(t/τ ) dt ∈ B. Thus, r(t) must be the CM and, consequently, must not be singular. As mentioned earlier, solvable models, e.g., the CC, show that it is not the case, which excludes them from further considerations. The distinguished role of the CM functions in the theory of relaxation phenomena has been noted a long time ago, starting with early observations concerning viscoelasticity [22] and rheology [9]. Additional links between the CM properties and relaxation phenomena have been provided by the studies of Green's functions [7] and of the extended family of the Mittag-Leffler-type functions [19,30,44,52,77,123], describing of the non-Debye relaxation.
Besides their well-established place in the classical mathematical analysis [14], the CM and B functions are objects which play an important role also in the probability and the theory of stochastic processes [113]. Thus, their appearance in the formalism signalizes existence of non-accidental relations between probabilistic tools and analytical methods of fractional and operational calculi, if applied to the relaxation theory. Requiring probabilistic interpretation of Eq. (II.4), besides of the indication that J (t) should be CM, leads to another important result. If we calculate the Laplace transform of J (t) with respect to t, J (z) = L [J (t), z], z ∈ C, then, after changing the order of integration, we get Under the condition that p(t) is integrable and nonnegative for all t ∈ R + , it means that J (z) belongs to the class of Stieltjes functions (S) which obeys well-known analytical properties, play important role in the Laplace transform theory [127] and, as we will see, turns out to be of a crucial importance in our further analysis. Notice that to get Eq. (II.9), we assume Re (z) ∈ R + . From the physical point of view, this assumption is too restrictive -to make our scheme compatible with physical data incorporated in Eq. (I.4), we do need to include also Re (z) = 0. The justification of this statement comes from the fact that the Laplace transform Eq. (II.9) for z = iω may be considered as the Fourier transform of a function vanishing on the negative semiaxis, i.e., the Fourier transform of a function Θ(t)f (t) being undefined for t = 0. The fundamental physical implication of the Fourier transform, namely the time-frequency correspondence, implies that J (z)| z=iω in such a way understood is closely related to the earlier introduced spectral functions being objects of direct physical meaning, as reconstructed from experimental data provided by methods of the broadband dielectric spectroscopy. But the spectral functions should be (and they are) defined also for ω = 0, in the case of a static field. Mathematically, it means that we have to point out the value J (0) to redefine the singularity of ∞ 0 Here the Plemelj-Sokhotski formula [137] enters the game with the Cauchy principal value being introduced. Consequently, we should extend Eq. (II.9) to the form which explains different shapes of the spectral functions relevant for the HN and JWS models. Summarizing the results obtained so far -using intuitive, and from time to time also hand-waving arguments, we can expect that functions admittable to describe non-exponential relaxation phenomena should be searched among functions belonging to the CM class or eventually its subclasses. Natural justification of such a conjecture is that the CM functions obey characteristic properties of relaxation functions, namely they are always non-negative and non-increasing. Infinite differentiability and sign-alternating derivatives, which in fact define complete monotonicity, are impossible to be reliably verified for functions known from phenomenological analyses only. Intuitively, to require the CM property shared with the exponential function used to describe the standard Debye decay, can shed more light on the problem. In Secs. V and VI we will see that it is really the case. Now, leaving aside the CM functions, we are going to convince the readers that the most promising factor which will enable us to go forward is to consider memorydependent kinetics giving dynamical laws governing the relaxation phenomena.

III. MEMORY DEPENDENT KINETICS
As signalized in the Introduction, the characteristic feature of the relaxation phenomena is that the relaxing system exhibits delayed response to the changes of external conditions. A typical example of such a behavior is extended over time approach to equilibrium after switching off factors perturbing the system during its previous history. In other words -approaching equilibrium, even in the case of absent external influence, needs time and the process under consideration does not satisfy the instantaneous response principle. This means giving up the time-local evolution rules and strongly suggests looking for prospective challengers which would be used to replace the standard time-local evolution equations. The "first choice" solution of the problem is to take into account the time non-local analogues of basic equations, in particular to consider integro-differential equations with kernels suited to mimic memory effects occurring in the system. Thus, to study the non-Debye relaxation phenomena, as well as other kinetic phenomena taking place in complex systems, e.g., anomalous diffusion, we should to reformulate the standard time-local approach. The aim is to construct a new theoretical scheme which has the time non-locality built-in from the very beginning, and leads to new evolution laws modified by various memory effects. The cornerstone of such new formalism is classical works of R. Zwanzig [131,132] who, after paying attention to the fact that the formalism of instantaneous response did not work in many physically interesting situations, proposed the self-consistent construction of kinetic equations respecting causality and taking into account delayed response of the system. The essence of mathematical encoding such delays is to replace point-like operations, e.g., the usual multiplication of functions, by integral operators taken in the form of time convolutions of memory functions and solutions looked for. Thus, it is justified to say that such modified evolution equations are the "time-smeared" standard ones. Almost 40 years later, at the very beginning of the current century, it was I. M. Sokolov [116] who shed new light on the ideas underpinning Zwanzig's seminal works. Sokolov showed that the use of generalized non-Markovian Fokker-Planck (FP) equations with built-in memory kernels leads to considerable progress in understanding anomalous transport and non-Debye relaxations. Subsequent years brought new research techniques and huge amount of earlier unreachable experimental data. These, urgently needed to be theoretically explained, pushed forward theoreticians' efforts and soon led to the development of effective mathematical tools rooted either in mathematical analysis, e.g., provided by the fractional calculus, or coming from the probability/stochastic processes theory. Taken separately or merged together, these methods have formed a toolbox whose usage has enabled physicists to solve equations of anomalous diffusion and non-Debye relaxations for a quite large number of problems [21,40,107,108,118].

A. The basics
Restricting Zwanzig's approach to the problems depending on the time only, e.g., to the simplest analysis of dielectric relaxation [138], but holding on to the spirit of seminal paper [131] we can write down the evolution equation as where t ∈ R + . The linear operator O acts on the time variable but is not simple multiplication by a timedependent function, and is the object which represents a more general characterization of non-locality in time. If O is represented as an integro-differential operator with the kernel O(t−ξ)B, B = const, then Eq. (III.1) becomes i.e., has a typical form of the Volterra equation. It is convenient to rewrite Eq. (III.2) as coming directly from the integral master equation [120,121] If we take the Laplace transforms (with respect to t) of Eqs. (III.2) and (III.3) or (III.4) and symbolically write down these equations as is to be specified, cf. [21]. Looking for solutions of (III.2) -(III.4) we are interested in the functions n(t) which are normalizable and non-negative, i.e., fulfill a minimal set of requirements needed to interpret them as the relaxation functions. The crucial property to be shown is non-negativity -we will see that requiring this property allows to find conditions which have to be put on the memory functions. The theory of integral equations, more precisely the possibility of a dichotomic description of the system using either integral or differential equations, suggests that it would be interesting also for Eq. (III.4) to find its integro-differential analogue. To achieve this goal let us notice that if we have given Eq. (III.3) governed by the kernel function M (t) then we may ask for a function k(t) defined by the convolution which, if required to be satisfied for any t > 0, is known in the theory of integral equations as the Sonine relation [53,54,74,75,82]. Transformed to the Laplace domain Eq. (III.5) becomes k(z) M (z) = z −1 from which it follows: General properties required from M (t) in order to satisfy basic conditions of the Laplace transform theory guarantee that Eq. (III.5) defines the couple (M, k) uniquely. Moreover, the symmetry of Eq. (III.6) with respect to the interchange M (t) ↔ k(t) means that if M (t) is interpreted as a memory function then k(t) does share this property: if Eq. (III.5) holds then any suitable M (t) which governs Eq. (III.3) has its partner memory function k(t) entering the kernel of integro-differential equation Equation (III.7) is equivalent to Eq. (III.3) and in mathematics it is known as the Caputo-Djhrbashyan problem [71].  [39,43]. Duality in the description of the same process using either Eq. (III.3) or Eq. (III.7) was noticed and studied in many papers [39,43,107,108]. Those were first of all addressed to the anomalous diffusion (through studying various applications of the generalized Fokker-Planck equation [107,108]) but led to results by no means limited to this class of phenomena [39,43]. To name different points of view on the problem we list: i.) solvability of the Cauchy problem for Eq. (III.7) studied either using advanced partial differential equations methods [71] or operator calculus tools developed for the Volterra-type evolution equations [83] and leading to the appearance of subordination [12,13], ii.) investigations devoted to the role played by the Sonine relation and resulting duality of memories [39,41,43,45,53,122] and iii.) efforts oriented on understanding mutual relations between the so-called deterministic and probability/stochastic approaches to anomalous kinetic phenomena [40,107,108,120,121].
Recall that restricting the kinetic problem under consideration to be depending on the time only means that the action of the FP operator reduces to multiplication by a constant factor B. Thus we arrive at equations governing the relaxation phenomena which provide us with other examples of non-Markovian processes [118,121]. It has to be noted that Eq. (III.4) enables us to write down the spectral function, i.e., the Laplace transform φ(z) of the response function φ(t) = − dn(t)/ dt in terms of M (z): which gives the mutual relation between the spectral functions and evolution kernels responsible for the memory effects. This observation has far-going consequences. From the physicists point of view Eqs. (III.8) and (III.9), if taken for z = i ω with ω identified as the frequency of harmonic field used to polarize the sample, couple in a unique way evolution equations introduced on theoretical background with phenomenologically known spectral functions given by formulae which classify the standard non-Debye relaxations as the CC, CD, HN, and JWS patterns (cf. [39,41,43,45]).
B. The stochastic approach -a few general remarks As mentioned, Eqs. (III.3), (III.4), and (III.7) are the Volterra-type equations [49]. They may be considered as equations which, because of a possible introduction of "atypical" kernels M (t) and k(t), go beyond the fractional differential equations conventionally used to describe the anomalous diffusion and non-Debye relaxation phenomena. In turn, Eq. (III.4) may be identified as the master equation which governs some stochastic processes underlying mechanisms of anomalous kinetics [121]. Recall that the memory functions M (t) and k(t) play a dual role in our scheme -they determine the mathematical structure of equations under study and provide us with links to observational data fitted in the case of relaxation experiments by phenomenological spectral functions φ(iω). However, if we restrict ourselves to the experimental data only, then usually we are not able to acquire knowledge sufficient to determine basic processes underlying physics of considered phenomena. To proceed further we do need more detailed information, also coming from mathematics, expected to be obtained from analysis of properties which the memory functions should obey. Obvious requirements that the memory functions have to be non-negative and non-increasing are insufficient to judge the existence and physical applicability of solutions as well as to find their interpretation. Some extra conditions, like the Boltzmann fading memory concept, did not satisfy initial expectations due to the lack of mathematical precision [4,5]. Clarification and resolution of many doubts appearing in the statistical description of both anomalous transport and relaxation phenomena came with using advanced probabilistic methodology, first of all with identification of just mentioned phenomena with stable stochastic processes.
Stochastic point of view (cf. e.g., [80,81,118,120] and [121,Chs. 4.1,4.3]) led to methods which made it possible to analyse non-Debye relaxations (as well as other anomalous kinetics processes) using probabilistic concepts like subordination [16,17,28,114], infinitely divisible distributions [76] and applications the theory of B and S functions, see [14,113] and Sect. IV A for a brief tutorial. Simultaneously, investigations based on less popular methods of mathematical analysis (fractional and operational calculi) made the integro-differential equations (III.3) and (III.7) much better understood if they were governed by kernels whose the Laplace images belong to the S class [40,71,72]. Here we would like to turn once more the readers attention to the importance of the Sonine relation and the special role of the S functions hidden behind it. Indeed, if one requires that M (z) and k(z) of Eq. (III.6) belong to the same class of functions then the most natural and convincing possibility to satisfy this requirement is to put both of them in the S class. Additionally, belonging to the S class directly links the memory functions with the complete Bernstein (CB) functions and the Laplace exponents, objects which characterize PDFs relevant for infinitely divisible stochastic processes [41,43,45,114,121].

C. Integral decompositions and subordination
To keep the forthcoming construction as general as possible let us observe that for the relaxation phenomena n(z) from Eq. (III.8) may be rewritten in the form (III.10) To get Eq. (III.10) the Schwinger parametrization was used and Ψ(z) = [ M (z)] −1 introduced. Ψ(z) is called the characteristic or the Laplace-Lévy exponent [120,121] which can be connected to the spectral function φ(z) via Eq. (III.9) : Taking the inverse Laplace transform of Eq. (III.10) we express it in the form of integral decomposition: the integrand being a product of the Debye relax- presenting the relation between the time t and integral variable ξ. The functions n D (ξ) and f Ψ (ξ, t) are given by normalized and non-negative functions which in probability theory means that they are the PDFs of "parent" [139] and "leading" processes, respectively. If these PDFs are independent Eq. (III.12) turns out to be the subordination. Within the latter scheme the integral variable ξ is interpreted as the internal time with randomize the physical time t.
The same results can be obtained by using the Efros theorem (Subsect. IV B) giving more flexibility than an application of the Schwinger parametrization. The use of the Efros theorem allows one to preset n(t) as where h(u) denotes the relaxation function of basic model. It is subordinated by the PDF f (u, t) which expresses the conversion between physical t and internal u times. The simplest realization of Eq. (III. 13) is for Other examples of Eq. (III.13) are presented in Subsecs. V F and VI F where for the HN and JWS models setting separately h(u) and f (u, t) we obtain dual representations of the same n(t).

IV. INTERLUDE -ELEMENTS OF THE MATHEMATICAL TOOLBOX
The wide applicability of the relaxation pattern attracts numerous authors to investigate it from various points of view. The goal of mathematically oriented research was to find the analytic expressions for the family of relaxations in the time domain. We mentioned in the Introduction that the problem was solved by R. Hilfer [61,62] who showed that solutions looked for may be expressed in terms of the Fox H-functions. This result, important for theoretical considerations, regrettably did not find deeper practical interest among experimentalists. It is not difficult to understand that as the Fox H-functions are special functions which for general values of parameters are defined symbolically either via the contour integrals of the Mellin-Barnes type [101] or by slowly convergent series. All together it makes relevant calculations difficult to understand, time-consuming and prone to mistakes emerging during hand-made manipulations. Thus, it is justified to emphasize that the use of the Fox H-functions formalism is not mandatory to invert analytically the Laplace transform in Eq. (I.4) if its RHS specializes in specifically chosen sets of parameters. As is shown in a series of research papers and books, e.g., [35,100] in such a case the suitable Laplace transform essentially simplifies and results in functions of the Mittag-Leffler family, much better known and easier to be handled both analytically and numerically. Moreover, relevant Mittag-Leffler functions obey properties placing them in special classes of analytic functions which for the positive value of the argument, besides of being nonnegative, also boils down to the completely monotone functions [14,113].

A. Non-negatively definite functions
The completely monotone (CM) functions are introduced as non-negative and infinitely differentiable (i.e., belonging to the C ∞ class) functions g(s) on R + whose all derivatives alternate in sign for any s ∈ R + : Notice that any CM function is a non-increasing and convex function. According to the Bernstein theorem [113], we can connect in a unique way CM and non-negative functions: The important property of the CMs is that the product of two CM functions is CM as well, i.e., CM · CM ⊂ CM .
The Stieltjes (S) functions are denoted as k(s) and form a subclass of the CM functions. On the real semiaxis they admit the integral representation can be interpreted as the Laplace transform of the CM function, i.e., the Laplace transform of the Laplace transform of the non-negative function [14,113]. This fact was pointed out and commented in the Introduction.
The Bernstein (B) functions are these non-negative infinitely differentiable functions b(s) on R + whose first derivative is CM. It means that (−1) n−1 b (n) (s) ≥ 0, n = 1, 2, . . . . denoted as c(s), s ∈ R + , are CB functions such that c(s)/s is the Laplace transform of the CM function restricted to the positive semi-axis, or, equivalently, the same way restricted Stieltjes transform of a positive function. In the paper we will use the following properties of CB functions: (a1) the composition of two CB functions is CB function, i.e., CB(CB) ⊂ CB; (a2) S functions can be obtained by making the composition of another S function with a CB function or the composition of a CB function with another S function, i.e., S(CB) ⊂ S and CB(S) ⊂ S; (a3) reciprocal (algebraic inverse) of an S function is a CB function. The inverse property is true as well: 1/ k(s) is CB function and 1/ c(s) is a S function for k ∈ S and c ∈ CB.
The illustrative example of differences between CM, S, B, and CB functions is provided by the power law function. This is presented in Tab. I. Two theorems of classical mathematical analysis are important for the mathematical toolbox used in the relaxation theory. First of them, which we abbreviate as Theorem 1, allows one to make analytic continuation of completely monotone functions to the complex domain. Thus, if its conditions are satisfied, we can freely jump between the complex variable z, e.g., z = iω and real s > 0. The second one is the Efros theorem [8,24,26,36,47,79,130] which generalizes the Borel convolution theorem for the Laplace transform. We will see that if applied to our purposes, the Efros theorem justifies integral decompositions and can be interpreted as a guidepost leading to the subordination approach [40,46]. Just mentioned theorems state as follows. Characterization of the Laplace transform of a CM, locally integrable function ( [19,49]) is given by Theorem 1: The Laplace transform f (z) of a function f (s), s > 0, that is locally integrable on R + and completely monotone, has the following properties: (a) f (z) has an analytical extension to the cut complex plane C \ R − . Conversely, every function f (z) that satisfies (a)-(c) together with (d) or (e), is the Laplace transform of a function f (s), which is locally integrable on R + and completely monotone on (0, ∞).
Proof. The proof of this theorem is presented in Ref. [49].
Remark 1. The Laplace transforms of CM functions satisfying the assumptions of Theorem 1 are defined for z ∈ C \ R − and because of conditions (d) or (e) they can be called the Nevanlinna functions [3,14,113]. For positive argument s they are, according to [113], the Stieltjes functions.
The Meijer G-functions.
The Meijer G-functions [94,101] are defined through the Mellin-Barnes integrals introduced as inverses in the sense of the Mellin transform of expressions being ratios of certain Γ functions. The definition of the Meijer Gfunction reads Parameters in Eq. (IV.7) are subject to conditions (i) They form a set closed under the Laplace transform. It means that its (direct and/or inverse) Laplace transform leads to another Meijer Gfunction: and . (ii) Below we quote two identities which will be used in our calculations: see [94,101].
(iii) Important procedure often applied in calculations is to convert the Meijer G-function G m,n p,q into the finite sum of the generalized hypergeometric functions p F q . This formula can be found in, e.g., [  The three parameter generalization of the Mittag-Leffler function [35] in the literature devoted to the fractional calculus and its applications is known as the Prabhakar function. It was introduced by T. R. Prabhakar [100] and reads: (IV.14) where the Pochhammer symbol (ν) r is given below Eq. (IV.6). For µ = ν = 1 it reduces to the (standard, called also one-parameter) Mittag-Leffler function E α (z) being the relaxation function of the CC model, which was first considered in [87][88][89]. For ν = 1 it becomes the two-parameters Mittag-Leffler function E α,µ (z) studied in 1905 by Wiman [129].
Remark 3. From the definition (IV.14) it follows that where introducing the δ-Dirac function allows us to avoid the pole at x = 0.
The comprehensive review of the family of the Mittag-Leffler functions can be found in [29,35,55]. Here, we shall quote a few formulae employed throughout the paper. Lemma 1. The asymptotics of E ν α,µ (−x α ) for small and large x is for x ≫ 1.
For x ∈ R and α = l/k such that 0 < l < k, the threeparameters Mittag-Leffler function E ν l/k,µ (x) can be represented as a finite sum of k generalized hypergeometric functions [39]: where ∆(n, a) is a special list of parameters defined below Eq. (IV.10). The Meijer G-representation of E ν α,µ (x) for α = l/k can be obtained by employing the Laplace integral of .
Example 1. The first of them is obtained from Eq. (IV.19) by taking µ = ν = 1. Due to the definition of the Meijer G-function in the numerator of (IV.7) the only non-vanishing terms are equal to m i=1 Γ(b i + s) whereas in the denominator survive only the terms p i=n+1 Γ(a i + s) which read ∆(l, 1) = 1 l , 2 l , . . . , l−1 l , 1 and ∆(k, 1) = 1 k , 2 k , . . . , k−1 k , 1, respectively. Delating "1" in both these sequences and adding "0" allows us to shift the upper and lower lists to ∆(l, 0) = 0, 1 l , . . . , l−1 l and ∆(k, 0) = 0, 1 k , . . . , k−1 k . In consequence, we get the integral representation of the Mittag-Leffler function first obtained by H. Pollard in Ref. [99] and independently derived in Refs. [37,126]: Example 2. As the second example, we take µ = αβ and ν = β. Then Φ β α,αβ (y) = αy −αβ Φ α (y)/Γ(β) for which we have [38,Eq. (11)], this is The substitution of Eq. (IV.19) into Eq. (IV.18) enables us to calculate the integral in Eq. (IV.10) and express it in the form (IV.21) is crucially important for solving problems which appear in the relaxation theory, e.g., it allows us to calculate the response function φ(t) from the experimental data encoded in phenomenological spectral functions φ(s) whose forms, see Eq. (I.4), match the RHS of Eq. (IV.21). Both heuristic, and rigorous arguments quoted in Sect. II suggest that looking for the relaxation functions we should focus our attention on the functions belonging to the CM class. Here it has to be said that physicists' interest in Mittag-Leffler functions as elements of the CM class is by no means limited to dielectric relaxation phenomena [19,30,52,77,123] but is also well-grounded in studies of other phenomena exhibiting the memory-dependent time evolution, to mention hereditary mechanics and viscoelasticity [51,78,[103][104][105]. Mathematicians' interest to study CM property of the Mittag-Leffler functions stems from the classical works of H. Pollard and W. R. Schneider who proved it for E α (−x) [99] and E α,β (−x) [115]. Currently the research in the field is inspired by problems and applications of the special functions theory [86,110,111], fractional calculus and fractional differential equations [35,78]. Recent results include proofs of the CM character obeyed by the Prabhakar function [77], see the lemma below, and by E ν α,µ (−at α ) itself if 0 < α ≤ 1, µ ≥ αν, and γ > 0, as shown in [44]. CM property of E ν α,µ (−at α ) explains restriction β < 1 put on the parameter β in the Prabhakar function to preserve its CM as a result of the multiplication of two CM functions.
Proof. Different proofs of this lemma were presented in [44,77,123].
Proof. Lemma 3 comes from direct calculations:

V. THE HAVRILIAK-NEGAMI MODEL AND ITS PHYSICAL CONTENT
In this section, we shall illustrate the general theoretical approach presented in Secs. I-III on the example of the HN relaxation model and its special cases, like the Debye, CC and CD models. The HN relaxation pattern was introduced in [56] to parametrize experimental data describing the frequency dependence of the complex dielectric permittivity ε ⋆ (ω) measured in polymers. Despite its purely phenomenological origin and apparent simplicity, the applicability of the HN model went far beyond its initial implementation. The model has appeared well-working for a much larger plethora of dielectric phenomena and has become the "first choice" method to analyse experimental relaxation data for various relaxation phenomena taking place in different complex systems, by no mere limited to those of condensed matter physics and materials science. Unexpected examples of its utility include the use of the HN function for monitoring the contamination in sandstone [106], or investigations of complex systems representing plant tissues of fresh fruits and vegetables, for which the HN relaxation in the frequency range 10 7 − 1.8 × 10 9 Hz was shown to be an useful tool of analysis [90].
Applying Lemma 1 to Eq. (V.15) we find the asymptotics of n HN (α, β; t) which is proportional to as is given in [30,45,62]. According to the general formalism developed in Subsect. III A, the relaxation function n HN (α, β; t) satisfies the integro-differential equations (III.4) and (III.7). They contain the memory functions M HN (α, β; t) and k HN (α, β; t) connected by the Sonine equation (III.5). The memory functions M HN and k HN in the Laplace domain can be determined algebraically using Eqs. (III.9) and (III.6): The The pseudo-operator on the LHS above belongs to the class of Prabhakar-like integral operators, and for the case under consideration is defined as [30, Eq. (B.23)] where we keep the notation of [30, Appendix B] and denoteṅ HN (α, β; t) = dn HN (α, β; t)/dt. It has to be pointed out that the pseudo-operator (V.21) must be distinguished from the pseudo-operator ( c D α + τ −α ) β considered in [30] as an alternative to the LHS of Eq. (V.21), with the Caputo derivative c D α (see [97]) sitting in. The difference is evidently seen if we take β = 1, which simplifies HN relaxation to the CC model. Then, it can be concluded that C (D α + τ −α )n CC (α; t) = ( c D α + τ −α )n CC (α; t) − τ −α , where n CC (α; t) = n HN (α, 1; t).
F. Two variants of the subordination approach to nHN (α, β; t) Let us start with the first type of subordination, namely Eq. (III.12), according to which we have represents the probability density distribution of (random) internal time ξ with respect to the laboratory time t. It is found by using Eq. (III.11). Inserting of Eq.
where the term exp(Bξ) in Eq. (V.23) canceled n D (ξ) = exp(−Bξ). The inverse Laplace transform L −1 [−, ·] in the integrand of (V.24) can be calculated by employing the Efros theorem give the Eq. (IV.5) of it with G 2 (z) = z −1 and q 2 (z) = τ −α + z α . That gives where z ÷ t and z ÷ u constitute the Laplace pairs for the first and second inverse Laplace transforms, respectively. The only term which depends on ξ is equal to exp(−ξBτ αβ z β ). Thus, inserting Eq. (V.25) into Eq. (V.24) and changing the order of integration we get The real function φ CD (β; u) ≡ φ HN (1, β; u) is the response function for the CD model which determines φ CD (β; u) = L −1 [(1 + T s) −β ; u] with T = τ α . Having this in mind we present Eq. (V.26) as where we used taken with minus sign the standard relation φ CD (β; u) = −ṅ CD (β; u) between the response function and the time derivative of relaxation function. Then, we employ the Leibniz formula which allowed us to shift the derivative over u from the relaxation function where f (α; u, t) can be expressed in terms of the onesided Lévy stable distribution Φ α (σ) with α ∈ (0, 1) and That means that we arrived at the second type of subordination, in which the CD relaxation appears, which plays the role of the parent process and is subordinated by the Lévy-like process described by the PDF f (α; u, t). This result explicitly confirms our previous claim [40] that the subordination description of the same process may be realized in different ways which eventually may be understood as an effect of "nested" processes.

VI. THE JURLEWICZ-WERON-STANISLAVSKI MODEL AND ITS PHYSICAL CONTENT
The next illustration which shows of applicability of the so far developed universal scheme presented in Secs. I-III is the Jurlewicz-Weron-Stanislavsky (JWS) relaxation model. The JWS relaxation model (for its brief description see [30,120]) complements and modifies the HN model leaving unchanged its general structure. The model was introduced in [67,117,125] to explain discrepancies emerging in some experiments between the results coming out from the Jonscher URL and those described by the HN pattern, see [120, Fig. 1]. (According to statistical enumeration measuring applicability of different models, the JWS pattern fits and reproduces the data for approximately 20% of relaxation experiments [65,119].) Therefore it is legitimate to analyse this model applying theoretical framework and methods developed in the current work.
To rephrase it in terms of functions implemented in the CAS we use Eq. (VI.2) and the Mittag-Leffler representation of the response function φ JWS (α, β; t). This enables one to express φ JWS (α, β; t) for rational α = l/k in terms of Meijer G-function and the finite sums of generalized hypergeometric functions. With the help of Eqs. (IV. 20) and (V.3) we have The asymptotic behavior of the Mittag-Leffler functions presented by Lemma 1 result in [30,Eq. (3.45)], that is Because φ JWS (α, β; s), s > 0, has the CM character for β ≤ 1/α, then thanks to the Bernstein theorem (Sect. IV A) we know that φ JWS (α, β; t) is a non-negative function in the same range of parameters α and β. That is confirmed in Fig. 4. C. The probability density gJWS(α, β; ξ) By analogy with to Subsect. V C we can represent the spectral function φ JWS (α, β; z), z = iω, as follows which for z = s ∈ R + is a S function [45]. The Schwinger parametrization enables us to write φ JWS (α, β; z) = L [φ JWS (α, β; t); z], where from Eq. (II.6) it follows that φ JWS (α, β; t) = L [p JWS (α, β; ξ); t] with p JWS (α, β; ξ) = ξ g JWS (α, β; ξ)/τ . Hence, Derivation of the exact form of g JWS (α, β; ξ) for α = l/k relies on calculating the inverse Laplace transform sitting in Eq. (VI.6). To do this we employ the Meijer G-representation of the JWS response function. That gives The first inverse Laplace transform in Eq. (VI.7), L −1 [δ(t); ξ], can be obtained by using the so-called limits representation of δ-distribution, namely δ(t) = lim ǫ→+0 ǫ t ǫ−1 /2 for t ∈ R + [95], and next changing the order of the inverse Laplace transform and the limit. That leads to To calculate the second inverse Laplace transform we use Eq. (IV.10): The last formula becomes more readable if the Meijer G-function is expressed in terms of the finite sum of the hypergeometric functions p F q according to Eq. (IV.13).
(completed with a suitable initial condition). The pseudo-operator in the above equation reads (VI. 18) It was introduced in [118,119] and discussed in Appendix B of [30]. For β = 1 and α ∈ (0, 1) it gives the evolution equation written in terms of fractional derivative in the sense of Riemann-Liouville derivative D α , α ∈ (0, 1) [97], used instead of the fractional derivative in the Caputo sense. Using the link between these derivatives, i.e., , we get the same formula.
F. Two kinds of subordination approaches of nJWS (α, β; t) As the first type of subordination we take the Debye subordination, i.e., Eq. (III.12) with fΨ(ξ, t) given now by After the cancelation of the Debye relaxation function by exp(Bξ) coming from Eq. (VI.19) we get To calculate the inverse Laplace transform in the integrand of Eq. (VI.20) we apply once more the Efros theorem but this time with G 2 (z) and q 2 (z) given below Eq. (V.24). Thus, the RHS of Eq. (VI.20) becomes Inserting it into Eq. (VI.20) and changing the order of integration we can simplify the calculations. It results in The next observation is where we applied Eq. (VI.1) for α = 1. Then, it turns out from that φ MCD (β; z) = 1 − (zτ α ) β φ CD (β; z). Since L −1 [1; u] = δ(u) and the inverse Laplace transform of the spectral function is equal to the response function, then Eq. (VI.22) is expressed as δ(u)− φ MCD (β; u). Consequently, n JWS comes out as (VI.23) The last steps to complete the calculation are: the use φ MCD (β; u) = −ṅ MCD (β; u), the Leibniz formula, and, finally, rewrite the RHS of Eq. (VI.21) as n JWS (α, β; x, t) = with f (α; u, t) given by Eq. (V.29). Equation (VI. 24) can be interpreted as the second type of subordination in which n MCD (β; u) is subordinated by f (α; u, t). We see that the HN and the JWS relaxation models lead to at least two types of subordinations: one described by Eqs. (V.22) and (V.28) for the HN model, as well as second describe by Eqs. (VI.20) and (VI.24) for the JWS model. Looking for physical interpretation of this fact enables us to suspect that with the growing complexity of the relaxing system, the simple partition of the process into two components, namely the parent and leading processes is not enough to reflect and understand all properties of the system, in particular it may be necessary to take into account the fact that the parent and leading processes can have a complex structure on their own.

VII. SUMMARY
The broadband dielectric spectroscopy allows us to obtain experimental data which, if fitted to the spectral functions φ (·) (iω), enable us to classify the latter as corresponding to one among of the standard relaxation models. Simpler of them, the CC, CD, and MCD models, emerge as reductions of three parameter ones, namely the HN and JWS patterns. HN and JWS models reduce to the CC relaxation for α ∈ (0, 1) and β = 1. For α = 1 and β ∈ (0, 1) the HN model goes to CD pattern, whereas the JWS pattern leads to the MCD model. All these spectral functions were tabulated in Tab. II from which one sees that φ JWS (α, β; iω) + φ HN (−α, β; iω) = 1. (VII.1) The knowledge of phenomenologically found spectral functions is crucial for our investigations. From one side, the methods of dielectric relaxation theory and the Laplace transform enable us to recover the response φ(t) and relaxation n(t) functions, see Eqs. (I.1) and (I.2), from the knowledge of spectral function. Obviously, the reverse procedure is also justified -we can transform φ(t) and/or n(t) into φ(iω). For the readers convenience we itemized φ(t) and n(t) in Tabs. III and IV which give φ JWS (α, β; t) + φ HN (−α, β; t) = δ(t) (VII.2) and n JWS (α, β; t) + n HN (−α, β; t) = 1.

(VII.3)
From another side the importance of spectral function φ(iω) comes from the fact that its knowledge enables us to find the memory functions k(t) and M (t). The latter are basic objects used to determine the time evolution of the system under study, see Eq. (III.9). To develop this idea we took into account that the memory functions M (t) and k(t) are mutually related by the Sonine equation. This implied that the relevant evolution equations, see Eqs. (III.3)/(III.4) and (III.7), led to the same results. Thus we concluded that the search for principles which govern the evolution may be done in two free chosen ways. The simplest is to consider integro-differential equations interpreted as memory dependent. For standard relaxation function these equations are presented in Tab. V. Here, we would like to pay the readers attention that this approach prefers the evolution equations for the HN and JWS models in the form which involves the pseudo-operators C (D α + τ −α ) β and (D α + τ −α ) β . The second observation concerning the spectral functions is that they are involved in the definition of characteristic functions Ψ(z), named also the Laplace-Lévy exponents. This correspondence goes through using the memory function M (t) and appears to be essential if one links the memory dependent evolution schemes and subordination approach, see Eq. (III.11). Our approach introduces an essential, previously unnoticed, noveltyin the dielectric relaxation theory we can distinguish, at least two, subordination patterns. Besides of the subordination which comes from the Schwinger parametrization (called by us the basic one), see Eq. (III.12), we found its alternative coming from the Efros theorem. Both these subordinations are connected to the various choice of "internal" timing, see Tab. VI. For instant, the HN relaxation function n HN (t) can be build for from the Debye relaxation in which we introduce internal time ξ. The relation between the physical time t and internal time ξ is given by the PDF of leading process, here f Ψ(z) (ξ, t).
Another possibility to get n HN (t) is from the CD relaxation function in which the timing characterized by the PDF f (α; ξ, t).
φHN φJW S φCC φCD φMCD  nD(ξ) f (α; ξ, t) nCC (ξ) δ(t − ξ) CD nD(ξ) fCD(β; ξ, t) = fHN (1, β; ξ, t) nCD(β; ξ) δ(t − ξ) MDC nD(ξ) fMCD(β; ξ, t) = fJWS (1, β; ξ, t) nMCD(β; ξ) δ(t − ξ) imental data. These models satisfactorily describe relaxation phenomena characterized by one-peak (unimodal) behavior in the frequency domain. However there exist materials which for frequencies 10 5 Hz exhibit either slower decay or multi-peak behavior of polarizability. This means that so far studied simple relaxation models do not cover all possibilities thoroughly enough. As prospective challenger models used to describe experimentally more complex phenomena we mention the excess wing (EW) model [42,63,92] or models involving sums of standards relaxation patterns taken with different parameters [68]. We remark that the EW model preserves complete monotonicity in the time domain while the second approach may lead to the lack of this property and consequently to abandon mathematical methods which result in the complete monotonicity concepts and reflect seemingly imposing, but far from complete, physical interpretation of relaxation in terms of summing up the Debye decays.
We also emphasize that the two-parameters KWW model mentioned in the Introduction (see Eq. (I.5)) must not be discarded as historical and old-fashioned. If we reduce it to the short-time, i.e., power-law, its asymptotics can be successfully used to describe discharge of atypical capacitors in modeling special electric circuits important for electrochemistry [31,50] and biochemical processes [58][59][60]. It also serves to be a starting point comparison between exponential-and power-like time behaviors of relaxation functions [1,2,45].