Subordination and memory dependent kinetics in diffusion and relaxation phenomena

The concept of subordination, originally introduced in the probability and stochastic processes theories, has also appeared in analysis of evolution equations. So it is not surprising that we meet it in physics of complex systems, in particular when study equations describing diffusion and dielectric relaxation phenomena. Grace to intuitively understood decomposition of complex processes into their simpler and better known components, called parent and leading processes, subordination formalism enables us to attribute physical interpretation to integral decompositions representing plethora of solutions to anomalous diffusion and relaxation problems. Moreover, it makes investigation of properties obeyed by these solutions far easier and more effective. Using the Laplace-Fourier transform method to solve memory-dependent evolution equations we show that subordination can be naturally implemented in their solutions. The key to achieve this goal is the use of operational calculus merged with the application of the Efros theorem [1]. Adopting exclusively methods of classical mathematical analysis we are able to derive the memory-stemmed origin of subordination and build a bridge connecting functional analysis/operator calculus based methods of solving the evolution equations with well established stochastic and probabilistic approaches. With such a developed general formalism in hands we apply it to several models of anomalous diffusion and relaxation phenomena.


Introduction
Simple models of kinetic phenomena, to mention those treating about standard diffusion, heat flow or dielectric relaxation, fail if one is going to apply them in order to describe physical systems which behavior shows significant deviations from idealized patterns. Those are exemplified by the Gaussian spreading of diffusing particles, the heat transfer generated by Fourier's law or the Debye exponential decay of mutually independent polarized dipoles. Such deviations, commonly nick-named anomalies, characterize systems which exhibit properties stemming from more and more complex structure and no longer can be considered as composed of non-interacting constituents placed in homogenous environments. Consequently, basic assumptions used to build up their theoretical modeling must be revisited. Among new concepts being introduced to solve this task the main role is played by modifications of evolution equations, sometimes proposed to be quite radical. Intuitive and physically justified approach leading beyond the standard evolution equations is to replace the latter ones by their time nonlocal generalizations expected to take into account physical consequences of the time delay with which the system responds to external stimuli. Called "memory" effects such phenomena are nothing strange in physics and are frequently met in continuum mechanics, macroscopic electrodynamics and various transport phenomena. Considering memory dependent kinetic equations and analysing their solutions we anticipate to expand our knowledge of complex systems and to understand much better how possible non-local time evolution effects influence their behavior. Illustration of this approach is provided by generalized diffusion equations in which the standard time derivatives are replaced by their fractional counterparts. A typical example, often encountered in studies of diffusion-like phenomena, is using the Caputo derivative involving, by definition, convolution of the standard time differentiation and the power-like time dependent kernel proportional to t −μ with μ > 0. So the evolution equations governed by the Caputo fractional time derivatives may be considered as resulting from the smearing the time derivatives with power-like memory functions. From the other side using the Caputo fractional time derivatives is important for self-consistent stochastic interpretation of diffusion-like phenomena considered as evolution patterns coming from the continuous time random walk model (CTRW). Recall that CTRW describes the motion of a diffusing object, traditionally named the walker, which jumps are neither in time nor in space regular but instead distributed stochastically in both domains. The randomness is attributed to the sojourn of the walker staying in the point (x, t) and contains information about the probability distributions of the jumps lenghts and the so-called waiting times, i.e. randomly long periods of time which the walker spends at rest at the point x [29]. Physically observed consequence of modeling anomalous diffusion using waiting time distributions (or the memory functions) given by a single power-like function t −μ is proportionality of the mean square value x 2 (t) of x(t) (MSD) to the power-like function of time t μ . The exponent μ, μ ∈ (0, 2], is used as a criterion for classification of the observed types of diffusion: for μ ∈ (0, 1) the process is called subdiffusion, for μ = 1 it is normal diffusion, for μ ∈ (1, 2) is the superdiffusion while for μ = 2 it bears the name of ballistic motion. Besides of the just listed physical implications using the power-like memory kernels has deep mathematical consequences which cause that stable stochastic processes and the concept of subordination do appear in diffusion/relaxation physics. Modern example is provided by advanced generalizations of the standard diffusion equation based on merging the CTRW stochastic concepts with usage the time and space-fractional derivatives coexisting in the same evolution equation [28,29,51]. This approach throws a bridge over the probability/stochastic methods, theory of stable stochastic processes and the Laplace-Fourier transform toolbox used to solve the evolution equations. The model explains how to incorporate the idea of subordination, grown out from the probability theory, into diffusion theory. It achieves this through studying processes initially governed separately by the time and space-fractional derivatives (random walks along the time and the space lines, respectively) and à posteriori joins them into one stochastic process, for comprehensive exposition see [29]. Here the power-like shape of the memory functions, sitting as kernels in the Caputo derivatives, turns out to be decisive to judge nonnegativity of obtained solutions and show the stable character of the subordination. Among another paths to generalize diffusion equations we mention giving up the constant character of the diffusion coefficient [67] or modeling diffusive processes using evolution equations with more than one fractional time derivative, including also the case of so-called distributed order fractional derivatives [20,21,37]. Models using the time dependent diffusion coefficient, inspired by experimental data, have been used to explain differences in behavior of MSDs for short and long times -relevant examples are provided by the scaled Brownian motion with diffusion coefficient D(t) ∝ t α−1 , where t ∈ R + and α > 0, for which the MSD grows linearly at short times, x 2 (t) ∝ t, whereas at long times behaves as x 2 (t) ∝ t α [17,42,68,69,81]. We note that similar duality, but of different origin, occurs for the Cattaneo-Vernotte equation (CVE) [18,19,82] and its generalizations [6,7,23,32]. For instance, the MSD calculated for the CVE is proportional to t for short times (like it takes place for the Brownian motion) while for long times it reads x 2 (t) ∝ t 2 which is typical for the ballistic motion. The effect results directly from the fact that CVE is governed by two time derivatives of different order which can be regarded as a trace of the time delayed Fourier's heat conduction law [52]. Note also that using CVE means changing the type of equation from parabolic to hyperbolic one and opens possibility to describe diffusion and heat transfer with finite propagation speed [84].
Huge amount of earlier unreachable data collected during recent few decades in diffusion/relaxation experiments forced theorists' community to push forward the development of mathematical tools providing the basis for effective models. The best among them have appeared those rooted in advanced methods of fractional calculus entangled with those coming from the probability/stochastic processes theory, in practical applications completed next with sophisticated computer science and numerical recipes. Merged together, these methods have formed a toolbox which usage has enabled us to solve equations of anomalous diffusion and non-Debye relaxations for a quite large number of problems [22,33,70,72,76]. However, regardless of the progress achieved so far, some questions, among them those concerning the properties of memory functions and their influence on behavior of physical systems, still wait for comprehensive physical and mathematical clarification. Our objective is to touch and clarify this problem using the Efros theorem [5,24,25,30,38,86], being itself almost 90 years old result of classical operational calculus and integral transforms theory. The usefulness, and at the same time simplicity, of the Efros theorem prompts us that it sheds new light on subordination and helps a lot to understand this probabilistic concept from the mathematical analysis point of view. We do believe that our exposition will make subordination concepts more affordable and easier to understand for physicists, also for those working with experiment. We also hope that our results will be useful to push forward physical interpretation of anomalous kinetics and explain frequently used, but a little bit naive, statement that subordination describes some standard evolution but governed by "internal" time in enigmatic way related to the physical one, measured by laboratory clocks.
We consider our paper as a review but emphasize that our presentation is entirely focused on the memory dependent evolution equations and their solutions which, under the form of integral decompositions, do hide subordination. Exploiting solutions written down in terms of integral decompositions we learn how to relate experimentally measured quantities with memory functions chosen to describe the system and how to determine requirements which they should obey. To illustrate the methods being used we shall concentrate ourselves on the time smeared 1 + 1-dim Fokker-Planck (FP) equation without drift for which relevant solutions may be expressed in analytic form. Thus the structure of the paper goes as follows: in the forthcoming Sect. 2 we recall the formalism of memory dependent evolution equations, discuss how they are related to the general class of Volterra equations and introduce concept of the Sonine pair. In Sect. 3 we shall demonstrate how the Efros theorem enables one to show that integral decompositions routinely exploited to represent solutions to the diffusion and relaxation equations are straightforward consequences of using methods coming from the operational calculus. Perspectives which the Efros theorem offers are quite wide -its formulation provides us with a general scheme which allows to construct and interpret integral decompositions for various evolution equations and shows how to relate them with probabilistic notions. In Sect. 4 we study conditions which, if required from the memory functions, imply the positivity of solutions to the equations met in studies of the anomalous diffusion. Our analysis is backgrounded in the properties of completely monotone (CMF), Bernstein (B) and complete Bernstein (CBF) functions and their relation to the infinitely divisible PDFs. In Sect. 5 we switch the subject of our investigations to the non-Debye relaxations. Applying methods presented in Sects. 3 and 4 we analyze from the subordination point of view popular relaxation models, namely the Cole-Cole, Havriliak-Negami and Jurlewicz-Weron-Stanislavsky ones. Section 6 is devoted to summing up discussion of general properties of integral decompositions as well as to the formulation of concluding remarks. The paper is completed by five appendices explaining details of calculations made throughout.

Memory dependent equations: their origin and properties
The need for describing kinetic phenomena taking place in complex physical systems by using equations containing memory functions materialized in early 60's of the previous century when Robert Zwanzig in seminal papers [87,88] first noted that the formalism of instantaneous response does not work in many physically interesting situations. Next he proposed and developed the selfconsistent construction of kinetic equations respecting causality and taking into account delayed response of the system. The standard mechanism of encoding such delays is to replace point-like operations, e.g. 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 a way modified evolution equations are the time smeared standard ones. Almost 40 years later it was Igor Sokolov [75] who shed new light on ideas underpinning Zwanzig's seminal works. Sokolov showed that the use of generalized non-Markovian Fokker-Planck equations with built in memory kernels leads to considerable progress in understanding anomalous transport and non-Debye relaxations.
Following the spirit of [87] we write down the generalized FP equation as which is typical form of the Volterra equation. For our further considerations it is convenient to rewrite Eq. (2.2) in the form coming directly from the integral master equation is to be specified, cf. [22]. Looking for solutions of (2.2)-(2.4) we are interested in the functions p(x, t) which fulfill minimal set of requirements needed to interpret them as probability density functions (PDFs) demanded to be normalizable and nonnegative. The widely known example is the case of the memory kernel M(t) present in Eq. (2.3) assumed to be proportional to the power-law function. Such obtained equation, for L = D∂ 2 x , D > 0, was in details discussed in [70,72] in the context of the CTRW approach to anomalous diffusion. From the second of just quoted papers, cf. [72,Eq. (3.8)], we learn that the solution to this model may be assigned the MSD x 2 (t) given by The Laplace image of the memory functionM(z) determines also the remaining moments: for even n = 2m the moments x n (t) are proportional to L −1 {[z −1M (z)] m ; t} while for odd n's they vanish. At this point we want to pay the readers' attention that the model under consideration provides us with an example which directly links moments with the memory function. As will be seen from Eqs. (3.2) and (4.1) in Sects. 3 and 4 such relations appear naturally and allow to suppose that using methods of the classical moment problem may be helpful to recover information concerning conditions which the memory function M(t) should satisfy to guarantee the existence of physically meaningful p(x, t) [2,61,74]. The interpretation of the memory function in terms of the CTRW can be found in Ref. [89].
Given Eq. (2.3) governed by the kernel function M(t) one 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 [40,41,47,48,57]. Transformed to the Laplace domain Eq. (2.6) becomesk(z)M(z) = z −1 from whicĥ  [34,35]. Thus it is seen how important for completing and consistent of our scheme is the Sonine relation. From one side it unifies in mathematical sense memory dependent effects coming from the integral and integro-differential equations, namely Eq. (2.8) and Eq. (2.3), respectively. From the other, physical side, it emphasizes the equivalence of memory dependent, time smeared, response of the system (Eq. (2.8)) with the time smeared evolution law (Eq. (2.3)).
Duality in description of the same process using either Eq. (2.3) or Eq. (2.8) was noticed and studied in many papers treating the problem from different points of view -e.g., i) solvability of the Cauchy problem for Eq. (2.8) studied using either standard methods of differential equations theory [45] or abstract operator/semigroup formalism based search for solutions of the Volterra equations, developed in [66] and applied to diffusion problems by Emilia Bazhlekova [9][10][11] whose approach has allowed to introduce and explain subordination on the basis of mathematical analysis considerations [12][13][14], ii) investigations devoted to the role played by the Sonine relation and resulting duality of memories [32, 34-36, 40, 79] and iii) efforts oriented on understanding mutual relations between the so-called deterministic and probability/stochastic approaches to anomalous kinetic phenomena [33,70,72,77,78].
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 another examples of non-Markovian processes [76,78] which differential form is The relaxation function n(t) counts dipoles which survive depolarization during the time (0, t) ∈ (0, ∞) and if normalized evolves from n(0+) = 1 to n(∞) = 0. Using the Sonine relation we can transform Eq.
It has to be marked that Eq. (2.10) enables us to write down the Laplace transform φ(z) of the spectral function φ(t) = − dn(t)/ dt in terms ofM(z) Looking at Eq. (2.11) from the physicists point of view we would like to note that this equation, if taken for z = i ω with ω identified as the frequency of harmonic field used to polarize the sample, couples theoretical formulae (cf. [32,[34][35][36]) with phenomenological models of the spectral functionsφ(iω), namely the standard non-Debye relaxations: the Cole-Cole, Havriliak-Negami and Jurlewicz-Weron-Stanislavsky patterns, which are obtained as fits to experimental data. Basic equations (2.3), (2.8) and (2.10) are the Volterra type equations [39,66]. They may be considered as equations which, because of permissible introduction of "atypical" kernels M(t) and k(t), go beyond the fractional differential equations conventionally used to describe the anomalous diffusion and relaxation phenomena within the framework of fractional calculus approach. 2 In turn, Eqs. (2.4) and (2.9) may be identified as master equations which govern some stochastic processes underlying mechanisms of anomalous kinetics [78]. Understanding properties of memory functions M(t) or k(t) plays important role in mathematical analysis of equations under study and, as it was mentioned earlier, provides us with links to observational data, either to the MSD (and eventually higher moments) in the diffusion case or to the spectral functionφ(iω) for the case of relaxation experiments. 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. We do need additional information practically available only from detailed analysis of properties which the memory functions M(t) and k(t) should obey. Obvious requirements that memory functions have to be nonnegative and nonincreasing are insufficient to judge the existence and physical applicability of solutions as well as to find their interpretation. Some extra conditions, like the fading memory concept, were tried to be put on the memory functions but due to the lack of mathematical precision they did not satisfy initial expectations [3,4]. The progress has come with using seemingly distant methods connected by the concept of subordination -from the one side linked with the probabilistic methodology, first of all identification of anomalous kinetic phenomena with stable stochastic processes (cf. e.g. [55,56,76,77,85] and [78,Chs. 4.1,4.3]) and from the other side with the operator theory equipped with ensuing functional analysis tools like resolvents [9][10][11]. Both approaches lead to methods which make it possible to analyse and solve anomalous diffusion and non-Debye relaxation problems. For probability-oriented researchers the crucial step is using concepts of stable and infinitely divisible distributions [49] merged with subordination technique of [16] and applications of the Bernstein (BF) and Stieltjes (SF) functions theory, see [15,73] and Appendix 1 for a brief tutorial. Investigations based on operator methods and functional analysis tools made the integro-differential equations (2.3) and (2.8) better and better understood especially if they were governed by kernels which Laplace images are power-like or, more generally, belong to the class of SFs, the latter being close relatives of BFs [9-11, 33, 45, 46]. Here we would like to pay readers' attention to two facts: -The Operator methods clarify the meaning of subordination as an universal method which enables one to find solutions to the evolution equations, or, more precisely, to solve complicated Cauchy problems, using solutions known for their much simplified analogues. 3 -In analysis of the Volterra-like evolution equations special care has to be paid for the Sonine relation and the role of the SFs hidden behind it. Indeed, if one requires thatM(z) andk(z) of Eq. (2.7) belong to the same class of functions then the overriding possibility to satisfy this requirement is to put both of them in the Stieltjes class. Additionally, belonging to the Stieltjes class directly links the memory functions with the BFs and Laplace exponents, objects which characterize PDFs relevant for infinitely divisible stochastic processes [34][35][36]78].

Integral decompositions: how the Efros theorem hides subordination
The standard method to solve either Eqs. inverting the algebraic expressions forp(κ, z) orn(z). We shall show that completing this procedure leads to the so-called integral decompositions used for almost 30 years by physicists working in the field of kinetic processes [8,26,55,56,75] and constructing their mathematical description which at that time arose from merging the Langevin-type evolution equations with the CTRW approach. Leaving aside, but not for a long time, the approach based on the probabilistic/stochastic processes methods we shall focus our attention on the exposition promoted in the Introduction, i.e., based on the study of evolution equations. To begin with we propose the leading role in the forthcoming investigations to be played by generalized FP equations and check if solutions p(x, t) of FP equations (2.3), (2.4) or (2.8) may be represented as integral decompositions given by the formula in which functional forms of h(x, ξ) and f (ξ, t) are to be determined from the primary equation. The keynote for physical interpretation of Eq. (3.1) is that it separates out the x and t dependence sitting in p(x, t) and replaces it, through a convolution type relation, by entangling the functions h(x, ξ) and f (ξ, t) temporarily playing roles of auxiliary objects not required to satisfy any special conditions. Under the additional assumption that h(x, ξ) and f (ξ, t) are PDFs the convolution Eq. (3.1) reflects the composition of independent probabilities, the so-called Bayes formula. Possibility to assign p(x, t) the interpretation of a probability distribution related to the process which emerges as composition of stochastic processes has become the crown argument supporting the importance of the probability/stochastic approach to kinetic phenomena and has opened possibility for the subordination formalism to get into the game [26,55]. We remark that although known and applied in practice for many years the origin, interpretation and possible constraints, either put on constituents of integral decompositions separately or demanded from Eq. (3.1) to be a PDF as a whole, are still the subject of ongoing research [22,33].
To derive Eq. (3.1) and to study its mathematical and physical content we begin either with Eq. (2.3) or Eq. (2.8) treated as equivalent under the provision of the Sonine relation. Assume also that the FP operator equals to D∂ 2 x which will enable us to present results in closed form. Taking the FL transform of Eqs.
Next, we can proceed two-fold way -either to invert first the Fourier transform κ ÷ x or to begin with inverting the Laplace transform z ÷ t. The inverse Fourier transform of Eq.
with h(x, ξ) being the inverse Fourier transform ofh(κ, ξ ) which in the (κ, ξ )-space The above shows that the integral decompositions appear as straightforward results of the Efros theorem. However, because of ambiguous choice of functionsĥ(x, z), q(z) andĜ(z) the representation of p(x, t) may be given, for the samek(z) orM(z), in different functional forms. It should be also noted that the Efros theorem scheme, Eqs. In what follows we will explain how to relate the set ofà priori chosen functions {ĥ(x, z),q(z),Ĝ(z) } with physical models.

Example 1: the Gaussian-like propagation
Our goal is to solve Eq. (3.2) using the Efros theorem. For that purpose we rename from which p N (x, t) is to be computed taking the inverse Fourier transform of and in a natural way links it to Eq.
completed with the initial condition N (x, 0) = δ(x). Equation (3.10) is the standard diffusion equation and N (x, ξ) is nothing else than its fundamental solution  N (x, ξ) represents a Markovian process. But there is no reason to claim that the evolution parameter ξ governing this process is the physical time t. Our only knowledge concerning this parameter is that it is nonnegative and appears as a partner of the Laplace variable z ∈ C. Moreover, in expressions expected to be physically meaningful ξ is integrated over and as such becomes a dummy variable which does not enter description of physical effects. From the other point of view the role played by ξ , or more precisely by the map f : t → ξ , is crucial for existence and properties of the integral decomposition Eq. (3.7). Recall that f (ξ, t) is the only object depending, through Eq. (3.6), on the memory functions. Having this in mind we are free to suppose that f (ξ, t) relates the physical time t to a variable ξ which enters the formalism as some "internal" time-like variable whose properties areà priori unknown and may be quite strange, e.g., its timing and flow are no longer uniformly distributed but represent random variables. If such an idea is accepted then the functions f (ξ, t) may serve as patterns how the "internal" time-like variables, whose origin and properties come from the memory effects, are functionally related to the physical time t measured by laboratory clocks. Obviously, not all patterns encoded by functions f (ξ, t) are eligible as mathematically and physically correct. Necessary mathematical condition is the existence of the inverse Laplace transform which may be deducted from analytical and asymptotic properties of the Laplace imagesM(z) andk(z). The minimal set of physical requirements demanded from ad hoc introduced quantity ξ to be called "new" time is that we have to deal with nondecreasing nonnegative valued function t → ξ which composition rule respects causality. However, if the "new" time is allowed to be stochastic then its properties may significantly differ from those of physical time, e.g., continuity of its flow may be relaxed and replaced by the so-called càdlàg condition defining functions (in maths classified as belonging to the Skorokhod space) which are right continuous with left limits -this allows e.g., to introduce the timing with jumps [22]. Note that for the memoryless case, i.e., M(t) = 1 ⇔M(z) = z −1 , the "new" time should reduce to that measured by laboratory clocks. Such "consistency test" we can formulate as follows: if M(t) = 1 (equivalent toM(z) = z −1 and k(t) = δ(t) andk(z) = 1, respectively) then f (ξ, t) = δ(ξ − t). It is easy to see that Eq. (3.6) satisfies this test and Eq. (3.7) leads to p N (x, t) = N (x, t) as it should be.
To proceed further in our search for the meaning of integral decompositions let us reanalyse previous steps. To get Eq. (3.7) we made two assumptions. First we adopted the condition (3.8) which led us to Eq.  4). However, the question concerning nonnegativity of this solution, i.e. its interpretation as a PDF, is pending. Here the help of probabilistic/stochastic tools appears useful. We have mentioned previously that integral decompositions which we are dealing with come out as results of convolution-like theorems. Simultaneously, if constituents of the integral decomposition admit interpretation as independent well-defined PDFs, then the integral decomposition Eq. (3.7) is well-defined PDF as well and may be read out as composition of two processes, in particular stochastic ones. Moreover, if f (ξ, t) is an infinitely divisible PDF then we can make the previous statement much stronger and say that f (ξ, t), called the leading process, subordinates the other process, called the parent one. For just presented example the parent process is determined by the shape ofh(κ, z) and for the case under study appears Gaussian. However, the choice ofh(κ, z) is by no means unique. To look closer on this we postpone to Sect. 5 further discussion of conditions under which subordination procedures do work and now will focus our attention on the problem how different choices of functionsh(κ, z),q(z) andĜ(z) influence factorization of the process into the parent and leading ones.

Example 2: the Gaussian vs Cattaneo-Vernotte propagation
To discuss this example we takek(z) →K(z) = τ z[γ (z)] 2 +γ (z) whereK(z) is a "composite" memory function which forms the Sonine pair with a suitably chosen memory functionM(z). Next assume that: i) Eq. (3.8), i.e.,h(κ, z) = (z + Dκ 2 ) −1 , is kept and the auxiliary functions are: ii)q 1CV (z) = zK(z) = τ [zγ (z)] 2 + zγ (z) and iii)Ĝ 1CV (z) =q 1CV (z)/z. A little bit overloaded notation is introduced to distinguish this example from another ones discussed later on. Withh(κ, z) = (z + Dκ 2 ) −1 and the substitutions ii) and iii) the inverse Laplace transform z ÷ t of Eqs. (3.2) and (8.4) reads with initial conditions p(τ ; x, 0) = δ(x) and ∂ t p(τ ; x, t)| t=0 = 0 and subject to the conditionη(z) = γ (z) 2 . The latter, as been shown in [32,33], guarantees that the smearings of the first and second order time derivatives are self-consistent. Comparing the FL transform of the fundamental solution to Eq. (3.14) with the integrand of the second inverse Laplace transform in Eq. (3.12) we see that they are the same. So the integral decomposition presented in Eq. (3.12) is the solution of the GCVE, Eq. (3.14). At this moment a little bit disturbing observation appears. According to Eqs. (3.12) such obtained solution of GCVE is expressed in terms of the convolution which involves the Gaussian (diffusion-like) propagator describing the infinite speed propagation. This fact remains undisturbed by any shape of f 1CV (τ ; ξ, t) and thus seems to be shared by the integral decomposition Eq. (3.12) [33]. So what about the case This argument stops to work if we deal with functions whose asymptotics z → ∞ reduces to products of exp [−ξ z α j ] with all α j ∈ (0, 1). For such a case the relevant inverse Laplace transforms exist and are expressed in terms of the Lévy stable distributions [30,31,62,63]. This is realized if the memory functions are modeled by the power-like monomialsγ (α; z) = z α−1 , α ∈ (0, 1). From the example discussed in Sect. 4, we will learn that the construction (and the positivity condition) works if α ∈ (0, 1/2] but is broken for α ∈ (1/2, 1]. Such conclusion is understood because for α ∈ (0, 1/2] the GCVE which we have begun with is the generalized diffusion equation with two Caputo fractional time derivatives which maximal order is 1. Alternative result of using the Efros theorem to find the integral decomposition relevant for Eq. (3.12) comes from writing down the latter in terms of h 2CV (τ ; κ, z) = τ z + 1 τ z 2 + z + Dκ 2 (3.15) and substitutions z →q 2CV (z) = zγ (z) ≡q N (z) andĜ 2CV (z) =q 2CV (z)/z ≡ G N (z). This time, if compared to the previously discussed case, we change the form ofh(κ, z) but leave simple linear relations linking auxiliary and memory functions.

This choice leads to the integral decomposition
with a = (D/τ ) 1/2 and the standard notation been applied: δ(·) symbolizes the δ-Dirac distribution, (·) means the Heaviside step function while the modified Bessel function of the first kind are marked as I 0 (·) or I 1 (·). Remark that the solution Eq. (3.18) is localized on the bounded support aξ > |x| and has been used for many years in studies of diffusion-like phenomena obeying the finite propagation velocity -see, e.g., [33,50,59,84].

Counterexample: the "wave-like" propagation
We remarked previously that for τ = 0 Eqs ( Thus the integral decomposition Eq. (3.1), with f W (ξ, t) given by which for the memoryless caseM(z) = z −1 gives i.e. the Gaussian spreading. It explicitly contradicts the result expected from Eq. (3.18) [33]. Inspection of the calculation indicates the mistake: the formal substitution q W (z) = [M(z)] −1/2 = √ z is not allowed because the inverse Laplace transform L −1 [ √ z; t] does not exist and the Efros construction breaks down. We end up this section with a global remark. The Efros theorem, much easier to follow and be applied than subordination schemes coming out from operator calculus methods, is a mathematical tool which enables us to represent solutions to the wide class of evolution equation in the form of integral decompositions. However, similarly to many other methods rooted in operational calculus, it leaves unsettled problems concerning both its mathematical rigor as well as physical consequences that result from it. Rigorous mathematics learns us that the use of operational formulae which appear the cornerstone of our approach must be treated with care every time -the last quoted example warns us about that. Another caution comes from physical point of view -the quantities that determine the correctness and physical utility of the results are to be interpreted as the memory functions. Nevertheless, in order to avoid premature conclusions one has to remember that even subtle changes of parameters characterizing the memory functions (for examples illustrating such situations see [22,33]) may lead to incorrect statements concerning the existence of solutions obeying the probabilistic interpretation.

Properties of memory functions
In this section we shall focus our attention on nonnegativity of functions given by integral decompositions of the form Eqs (3.5) and (3.6). 5 To make deeper insight to this problem we are going to look for possibly well-defined conditions which the memory functions M(t) and k(t) should satisfy in order to obtain normalizable and nonnegative solutions either to Eqs (2.3), (2.8) or GCVE. These solutions are given by p N (x, t), p 1CV (x, t) and p 2CV (x, t) and in what follows denoted by a shorthand as p j (x, t), j = N , 1CV, and 2CV. The normalization of p j (x, t)s is fulfilled automatically because of R p j (x, t) dx = L −1 [z −1 ; t] = 1 for all j's and t ≥ 0. More effort is needed to show that p j (x, t)'s are nonnegative. We will achieve that with the help of the Bernstein theorem which uniquely joints, by the Laplace integral, the nonnegative function p j (x, t) with the completely monotone function (CMF)p j (x, s), s > 0. We would like to pay the readers' attention that, due to [39,Th. 2.6], we can change the Laplace variable z ∈ C onto s > 0. So from now all memory functions are the Laplace integrals (being by themselves nonnegative valued functions of the real variable) of their time dependent archetypes. Special classes of such functions, namely the CMFs and BFs (see Appendix 1 for a brief tutorial to the subject or consult [73] for comprehensive review of the properties and mutual relations of CMF, SF and CBF) and their theory based on the Bernstein theorem has been shown to be very useful to prove positivity of solutions to numerous diffusion and relaxation problems [33,70,72]. where j = N , 1CV, 2CV. Nonnegativity ofp j (x, s) for all x is easily seen so we will focus our attention on its nonnegativity with respect to t which must be unraveled from properties in the Laplace domain. Because the product of CMFs is CMF (see the property (a1) in Appendix 1) and exponentiation of minus completely Bernstein function (CBF) is CMF (the property (a3)) then the simplest (but pay attention that sufficient only) assumption which promises to be effective reads From the property (a4) we learn that in such a case q j (s)/s belongs to the SF class, being itself the subclass of CMF. Hence, for (4.2) satisfied, we can call p j (x, t) the PDF. The requirement (4.2) implies that [q j (s)/D] 1/2 = λ −1 j (α; s) is nonnegative. Simple example is to takeM(α; s) = s −α , α ∈ (0, 1] which gives λ 1 (α; s) > 0 for all s > 0. Then, for the power law memory functions, we obtain that sp 1 (α; x, s) is given by the Laplace distribution.
As mentioned, the assumption (4.2) is a sufficient condition only. This means that if (4.2) is true thenp j (x, s) of Eq. (4.1) is CMF and thus p j (x, t) is nonnegative. Unluckily we must not be extrapolate (4.2) to the analogous criterion forq j (s). The property (a2), which says that the composition of CBFs is CBF, is one-sided theorem only. So (4.2) does not imply the completely Bernstein character ofq j (s). As an example let us takek(α; s) = s α−1 for which we get q N (α; s) = s α/2 and q 1CV (α; s) = (τ s 2α + s α ) 1/2 . From Table 3 it is seen thatq N (α; s) is CBF. Simultaneously, in [32, Sect. IV B] it has been shown that q 1CV (α; s) is CBF as well. Now, for α ∈ (0, 1) we know thatq N (α; s) = s α is CBF but forq 1CV (α; s) we must not claim that. This is becauseq 1CV (α; s) = τ s 2α + s α is CBF only for α ∈ (0, 1/2] and it is not CBF for α ∈ (1/2, 1) [32,33]. This example warns us that the assumption (4.2), although illustrative and helpful, does not give information enough to be critical if one looks for (at least) sufficient conditions to be put on the memory functions M(t) and k(t). Useful properties of the memory functions are not encoded in the square roots ofq j (s), j = N ,1CV, 2CV but in the functionsq j (s) so we do need strongest requirement than Eq. (4.2) is.

The second hint: infinite divisibility
To find more information aboutq j (z) we will go back to the integral decomposition (3.1) and look more carefully on the examples of functions h(x, ξ) and f (ξ, t) found in Sect. 4. The first of them is either the Gaussian N (x, ξ) or the fundamental solution to the CVE h 2CV (τ ; x, ξ). Both these functions are normalized, nonnegative, and, what is much less known, infinitely divisible. For the Gaussian (normal) distribution relevant proofs are obvious [49]. For h 2CV (τ ; x, ξ) its normalization and nonnegativity are demonstrated in [32] while infinite divisibility can be deduced from the fact that h 2CV (τ ; x, ξ) can be derived from the central limit theorem if we consider its behavior of |x| larger than O( √ t) [44]. The reduction of h 2CV (τ ; x, ξ) to N (x, ξ) was presented in [33, Sect. IIIA]. In the language of subordination we can call each of these functions as the PDF of some parent process: once it is the Brownian motion and another time it is the process described by the continuous persistent random walk model [32,43,53,54,83]. These PDFs represent the Markovian processes and thus do not contain any information on the memory functions. The only objects which hide this information are, according to the Efros theorem, functions f j (ξ, t), j = N, 1CV, 2CV. In the Laplace space they readf As presented in [22] the functions f j (ξ, t) and g j (ξ, t), are connected to direct and parametric integral decompositions, respectively. Checking their probabilistic interpretation we see that the normalization of g j (ξ, t) as well as f j (ξ, t) with respect to ξ is fulfilled automatically while their nonnegativity comes out from the Bernstein theorem and can be ensured by the only one requirement, namelŷ Indeed, this condition guarantees the completely monotone character ofĝ j (ξ, s) and f j (ξ, s) and thanks to the Bernstein theorem g j (ξ, t) and f j (ξ, t) are nonnegative. The requirement (4.6) implies also that which is seen from [73, Lemma 5.8 (or Lemma 9.2)]. So we can callĝ j (ξ, s) the infinitely divisible PDF. Basic information concerning infinitely divisible functions, their connection to CMF and CBF as well as just quoted Lemmas 5.8/9.2 are given in Appendix 4. Here we recall that according to [22] the above constitutes the sufficient condition for identification the integral decomposition with the subordination approach. In the time t domain the functionsĝ j (ξ, s) satisfy for μ = μ 1 + μ 2 the relation which comes from the infinite divisibility and means the semigroup composition property of g j (ξ, t).
Concluding this section we can state that the condition Eq. (4.6) which locate the memory functionsq j (s) in the CBF class joins in a synergetic frame theoretical schemes stemmed from stochastic, probabilistic and mathematical analysis methods. Functions belonging to the CBF class entangle seemingly distant objects: infinitely divisible stochastic processes characterized by distributions coming from the Lévy-Khintchine formula, the Laplace exponents, the semigroup composition law, the functions of the Pick-Nevalinna class and methods of the classical moment problem. In our opinion understanding of their mutual relations should be treated as signposts which guide future development of various branches of kinetic theory -to support this hypothesis we will show how the scheme developed so far is applied to the relaxation phenomena.

Integral decompositions: Applications to relaxation phenomena
In what follows we consider evolution equations describing the relaxation phenomena.  [49,62,63]. Note that the stability property (5.5) may be used to obtain the fractional differential equation [31].  [49, p. 108] this function is infinitely divisible PDF. Calculations leading to Eq. (5.6) are based on the Efros theorem; they are rather laborious and because of that we shift them to Appendix 5.

Discussion and conclusions
The Efros theorem, originally derived as a little bit tricky mathematical tool inspired by operational calculus and used to solve integral equations, appears a versatile method applicable to study kinetic phenomena described by the time-nonlocal evolution equations. Such type of equations we identify as modeling processes usually classified as the memory dependent phenomena. The latter class of physical effects breaks down the instantaneous response paradigm and for diffusion and dielectric relaxation leads to the time behavior rules different from the Gaussian spreading or exponential decay. We have shown that using the Efros theorem is simple and easily applicable method which enables us to write down solutions to the time smeared evolution equations as integral decompositions. The latter are convolutions admitting, under well understood Table 2 Examples of h(ξ ) and f (ξ, t) in the integral decompositions (3.1) discussed in Sect. 5 physical requirements, consistent probabilistic interpretation in terms of processes governed by some probabilistic/stochastic rules. Presented approach explains how starting from deterministic evolution equations, like the integral (master) or integrodifferential ones, we arrive at the area of statistical physics ruled by the probability and stochastic mechanisms. These links are by no means restricted to formal similarities and observations. Significance and deep meaning of the Efros theorem is determined by the fact that expressions it provides, if merged with subordination schemes and stochastic approach, open possibility to identify physical quantities, the memory functions in the case, with objects of purely probabilistic origin like the infinitely divisible processes, subordinators and Laplace exponents are. This way revealed relations shed new light on mutual connections between memory functions and subordinators as we can expect their applicability for memory functions much general than commonly used as kernels of fractional derivatives. In this context our approach goes beyond the subordination scheme introduced in the framework of functional analysis/operator calculus approach focused on rigorous solving the Cauchy problem for fractional diffusionlike equations and showing nonnegativity of relevant solutions [12][13][14]. Considering physical interpretation of our approach we see that links to the memory or spectral functions equip subordination with a meaning more general than ad hoc introduced non-measurable, and thus physically vague, "internal" time. Novelty worthy for further studies is paying attention to the fact that choosing factorization of the integral decompositions we select the splitting of the process into the parent and leading ones. Such procedure is by no means unique. Solvable examples show that the parent and leading processes may be different for the same PDF assumed to represent the experimental data. Mathematically it is quite clear, also physically is not unexpected -physical processes looking the same at first glance in fact may be rooted in different primary laws.
In Tables 1 and 2 we recall dualities of integral decompositions found and discussed in Secs. 3 and 5. Note also that separating out the parent and leading processes we restricted ourselves to simple models dependent on a single memory function. However, we are convinced that the methods presented so far will work also for much complicated models involving larger number of memory functions or governed by "nested" processes. Such situations seem to be typical situation for various transport and relaxation phenomena taking place in complex systems. Our last remark is that the subordination interpretation of the integral decomposition (3.1) flows out from the Efros theorem under the strong requirement that the components of integral decomposition (3.1), i.e. f (ξ, t) and h(x, ξ) (in diffusion) or h(ξ ) (in relaxation) are normalizable and nonnegative. Both these properties are ensured if, according to (4.6), the memory functionsq(s)'s are CBF. Although it means only a sufficient condition such statement has far going implications -grace to it we can interpret any suitable memory function as the characteristic (Laplace-Lévy) exponent related to some nonnegative stochastic process U (ξ ) such that its PDFĝ(ξ, s) g(ξ, s) = e −sU (ξ ) = e −ξq(s) , ξ > 0 is infinitely divisible and satisfies the semigroup property (4.8).
Among the important properties of CMFs we have (a1) the product of two CMFs is also CMF.
Note that Eq. (7.1) is the real-valued Laplace integral of a nonnegative function g(t) and in order to deal with the Laplace transform its argument has to be complex. Hence, G(s) and the Laplace transform of g(t), denoted asĝ(z) are different objects: the first of them is real function of s > 0 while the second is complex valued and depends on z ∈ C\R − . Knowledge of analytic continuation of G(s) →ĝ(z) is important because these are special analytic properties ofĝ(z) (known as Herglotz conditions, [2,15]) which determine, e.g., according to Theorem 2.6 of Ref. [39] quoted as Theorem of Ref. [60], conditions under whichĝ(z) is representable as the Laplace transform of a nonnegative measure defined on positive semiaxis. In majority of probabilistic applications we may restrict considerations to the real variable s and treatĝ(s) as the real function G(s) but to find the inverse Laplace transform ofĝ(z) we must have the variable z complex. The next class of functions needed in our considerations is that of complete Bernstein functions (CBF) [71,73]: c(s) is CBF, s > 0, if c(s)/s is the Laplace transform of CMF restricted to the positive semiaxis, or, equivalently, in the same way restricted Stieltjes transform of a positive function named also the Stieltjes function (SF). Note that all SFs are completely monotone i.e., SFs are subclass of CMF.
The following properties of CBF will be used in what follows: (a2) the composition of CBFs is CBF; (a3) the composition of CMF and CBF is another CMF; (a4) if c(s) is CBF then c(s)/s is SF.
CBFs form a subclass of the Bernstein functions (BF). These are defined as nonnegative functions whose derivative is CMF [71,73]: h(s) > 0 is BF if The example which illustrates similarities and differences between CMF, SF, BF, and CBF is the power function. This is presented in Table 3.

Theorem 1 IfĜ(z) andq(z) are analytic functions, and
as well as From the Efros theorem it appears that

The relation between the infinitely divisible distributions and the Bernstein-class functions: BF as well as CBF
The relation between the CBF and infinitely divisible function is expressed by [

Duality of integral decompositions
To show the passage between Eqs. (46) and (49)  To calculate the inverse Laplace transform in the integrand of Eq. (11.9) we apply the Efros theorem once more, this time withĜ 2 (s) andq 2 (s) as given below Eq. (11.2). Thus the RHS of (11.9) becomes which leads to Eq. (50).