From Random Times to Fractional Kinetics

In this paper we study the effect of the subordination by a general random time-change to the solution of a model on spatial ecology in terms of its evolution density. In particular on traveling waves for a non-local spatial logistic equation. We study the Cesaro limit of the subordinated dynamics in a number of particular cases related to the considered fractional derivative making use of the Karamata-Tauberian theorem.


Introduction
A study of a random time change in a Markov process X t was initiated by S. Bochner [Boch55] by considering an independent Markov random time ξ t . The resulting process Y t = X ξt is again Markov and is called a subordinated process. In a pioneering work [Kols], T. Kolsrud initiated the study of the general independent random time process. Later the concept of a random time change became an effective tool in the study of physical phenomena related to relaxation and diffusion problems in complex systems. We refer here to the section "Historical notes" in [MTM08]. An additional essential motivation for the random time change did appear in applications to biological models. The point is such that there exists a notion of biological time specific for each particular type of biological system and which is very different compared to the usual time scale employed in physics. One of the possibilities to incorporate this notion is related to a random time change. Moreover, this approach gives the chance to include in the model an effective influence of dynamical random environment in which our system in located.
An especially interesting situation appears for the case of an inverse subordinator ξ t . We will describe this framework roughly and leaving the details and more precise formulations in Subsection 2.2 below. The marginal distributions µ t of a Markov process X t describe an evolution of states in the considered systems and deliver an essential information for the study of the dynamics. We call that the statistical dynamics in contrast to the stochastic dynamics X t which contains more detailed information about the evolution of the system. The statistical dynamics may be formulated by means of the Fokker-Planck-Kolmogorov (FPK) evolution equation (weak sense) where L * is the (dual) generator of the Markov process on states. For an inverse subordinator ξ t and the process Y t = X ξt , denote ν t the corresponding marginal distributions. The key observation is such that the dynamics of ν t is described by the evolution equation where D ξ t denotes a generalized (convolutional) fractional derivative in time canonically associated with ξ t . This fractional Fokker-Planck-Kolmogorov equation (FFPK) gives the main technical instrument for the study of subordinated statistical dynamics. There is a well known particular case of the inverse to stable subordinators. In this case, such standard objects appear as Caputo-Djrbashian fractional derivatives, and all the well developed techniques of fractional calculus work perfectly. But in the case of general inverse subordinators we should think about proper subclasses, for which certain analytic properties of the related objects may be established, see e.g. [CKKW].
Note that there exist two possible points of view. We can start from an inverse subordinator and arrive in FFPK equation [MS06, MS08, Kol09, To2015]. Or, vice versa, we develop at first a notion of generalized fractional derivative and then search for a probabilistic interpretation of the solution. The latter was first realized in [Koc11]. For a detailed discussion of both possibilities see [CKKW].
The aim of this paper is to analyze the effects of random time changes on Markov dynamics for certain models of interacting particle systems in the continuum. For the concreteness, we will consider the important Bolker-Pacala model in the spatial ecology that is a particular case of general birthand-death processes in the continuum. The scheme of our study is the following. The FPK equation for the states of the model may be reformulated in terms of a hierarchical evolution of correlation functions. In a kinetic scaling limit this system of equations leads to a kinetic hierarchy for correlation functions and to a non-linear evolution equation for the density of the system. The latter is a Vlasov-type non-linear and non-local evolution equation. Considering a random time change by an inverse subordinator we arrive in the FFPK equation for correlation functions and to a fractional kinetic hierarchy. For a discussion of this approach and certain new properties of fractional kinetic hierarchy see [KK17].
A surprising feature appearing in the kinetic hierarchy is related to the dynamics of the density of the considered systems. The evolution of the density in the time changed kinetics is not the solution of a related Vlasovtype equation with a fractional time derivative as one may expect. In reality, this dynamics is a subordination of the solution to kinetic equation for the hierarchy in the initial physical time. A particular problem which does appear in this situation is related with an effect of the subordination for such special solutions as traveling waves which are known for the Bolker-Pacala model. In the special case of stable subordinators this question was studied in [DKT]. In the present paper we are dealing with certain classes of inverse subordinators for which the analysis of subordinated waves may be carried out.
We summarize our observation as follows: a heuristic consideration of a kinetic equation for the density with a fractional time derivative has no relation to the real dynamics in the kinetic limit of the time changed Markov evolution of the model. The correct behavior is given by a subordination of the solution to the kinetic equation in physical time.

General Facts and Notation
Let B(R d ) be the family of all Borel sets in R d , d ≥ 1 and let B b (R d ) denote the system of all bounded sets in B(R d ).
The space of n-point configurations in an arbitrary Y ∈ B(R d ) is defined by where | · | the cardinality of a finite set. We also set Γ (0) (Y ) := {∅}. As a set, Γ (n) (Y ) may be identified with the symmetrization of The configuration space over the space R d consists of all locally finite subsets (configurations) of R d , namely, (2.1) The space Γ is equipped with the vague topology, i.e., the minimal topology for which all mappings Γ ∋ γ → x∈γ f (x) ∈ R are continuous for any continuous function f on R d with compact support. Note that the summation in x∈γ f (x) is taken over only finitely many points of γ belonging to the support of f . It was shown in [KK06] that with the vague topology Γ may be metrizable and it becomes a Polish space (i.e., a complete separable metric space). Corresponding to this topology, the Borel σ-algebra B(Γ) is the smallest σ-algebra for which all mappings It follows that one can introduce the corresponding Borel σ-algebra on Γ (n) (Y ), which we denote by B(Γ (n) (Y )). The space of finite configurations in an arbitrary Y ∈ B(R d ) is defined by This space is equipped with the topology of disjoint unions. Therefore one can introduce the corresponding Borel σ-algebra B(Γ 0 (Y )). In the case of Y = R d we will omit Y in the notation, thus Γ 0 : The restriction of the Lebesgue product measure (dx) n to Γ (n) , B(Γ (n) ) will be denoted by m (n) , and we set m (0) := δ {∅} . The Lebesgue-Poisson measure λ on Γ 0 is defined by For any Λ ∈ B b (R d ), the restriction of λ to Γ(Λ) := Γ 0 (Λ) will be also denoted by λ. The space Γ, B(Γ) is the projective limit of the family of spaces (Γ(Λ), B(Γ(Λ))) Λ∈B b (R d ) . The Poisson measure π on Γ, B(Γ) is given as the projective limit of the family of measures {π Λ } Λ∈B b (R d ) , where π Λ := e −m(Λ) λ is the probability measure on Γ(Λ), B(Γ(Λ)) . Here m(Λ) is the Lebesgue measure of Λ ∈ B b (R d ).
For any measurable function f : Then, by (2.2), for f ∈ L 1 (R d , dx) we obtain e λ (f ) ∈ L 1 (Γ 0 , dλ) and We will make use of the following classes of functions on Γ 0 : (i) L 0 ls (Γ 0 ) is the set of all measurable functions on Γ 0 which have local support, i.e., H ∈ L 0 is the set of bounded measurable functions with bounded support, i.e., H ↾ Γ 0 \B = 0 for some bounded B ∈ B(Γ 0 ).
In fact, any B(Γ 0 )-measurable function H on Γ 0 is a sequence of functions )-measurable function on Γ (n) . On Γ we consider the set of cylinder functions F cyl (Γ). These functions are characterized by the relation F (γ) = F ↾ Γ Λ (γ Λ ).
The following mapping from L 0 ls (Γ 0 ) into F cyl (Γ) plays a key role in our further considerations: where H ∈ L 0 ls (Γ 0 ). See, for example, [KK02] and references therein for more details. The summation in (2.5) is taken over all finite sub-configurations η ∈ Γ 0 of the (infinite) configuration γ ∈ Γ; this relationship is represented symbolically by η ⋐ γ. The mapping K is linear, positivity preserving, and invertible, with (2.6) Here and in the sequel, inclusions like ξ ⊂ η hold for ξ = ∅ as well as for ξ = η. We denote the restriction of K onto functions on Γ 0 by K 0 . A probability measure µ ∈ M 1 f m (Γ) is called locally absolutely continuous with respect to (w.r.t.) a Poisson measure π if for any Λ ∈ B b (R d ) the projection of µ onto Γ(Λ) is absolutely continuous w.r.t. the projection of π onto Γ(Λ). By [KK02], there exists in this case a correlation functional κ µ : Γ 0 → R + such that the following equality holds for any H ∈ B bs (Γ 0 ): (2.7) Restrictions κ

Microscopic Spatial Ecological Model
Let us consider the spatial ecological model a.k.a. the Bolker-Pacala model, for the introduction and detailed study of this model see [BP97,FM04,FKK10,FKK12,FKKK15]. Below we formulate certain results from these papers concerning the Markov dynamics and mesoscopic scaling in the Bolker-Pacala model.
The heuristic generator L in this model is defined on a space of functions over the configuration space (2.8) Here m > 0 is the mortality rate, a − and a + are competition and dispersion kernels, respectively. See Section 5 for the conditions on these kernels in the present paper. We assume that the initial distribution in our model is a probability measure µ 0 ∈ M 1 (Γ) and the corresponding sequence of correlation functions κ 0 = (κ (n) 0 ) ∞ n=0 , see e.g. [KK02]. Then the evolution of the model at time t > 0 is the measure µ t ∈ M 1 (Γ), and κ t = (κ (n) t ) ∞ n=0 its correlation functions. If the evolution of states (µ t ) t≥0 is determined by the heuristic Markov generator L, then µ t is the solution of the forward Kolmogorov equation (or Fokker-Plank equation FPE), where L * is the adjoint operator of L. In terms of the time-dependent correlation functions (κ t ) t≥0 corresponding to (µ t ) t≥0 , the FPE may be rewritten as an infinite system of evolution equations The expression for the operator L △ is obtained from the operator L via combinatoric calculations (cf. [KK02]).
The evolution equation (2.10) is nothing but a hierarchical system of equations corresponding to the Markov generator L. This system is the analogue of the BBGKY-hierarchy of the Hamiltonian dynamics, see [Bog46].
We are interested in the Vlasov-type scaling of stochastic dynamics which leads to the so-called kinetic description of the considered model. In the language of theoretical physics we are dealing with a mean-field type scaling which is adopted to preserve the spatial structure. In addition, this scaling will lead to the limiting hierarchy, which possesses a chaos preservation property. In other words, if the initial distribution is Poisson (non-homogeneous) then the time evolution of states will maintain this property. We refer to [FKK10] for a general approach, other examples, and additional references.
There exists a standard procedure for deriving the Vlasov scaling from the generator in (2.10). The specific type of scaling is dictated by the model in question. The process leading from L △ to the rescaled Vlasov operator L △ V produces a non-Markovian generator L V since the positivity-preserving property fails. Therefore instead of (2.9) we consider the following kinetic FPE, and observe that if the initial distribution satisfies µ 0 = π ρ 0 , then the solution is of the same type, i.e., µ t = π ρt , t > 0. In terms of correlation functions, the kinetic FPE (2.11) gives rise to the following Vlasov-type hierarchical chain (Vlasov hierarchy) This evolution of correlations functions exists in a scale of Banach spaces, cf. [FKK12]. Let us consider the Lebesgue-Poisson exponent, defined in (2.3) as the initial condition. Such correlation functions correspond to the Poisson measures π ρ 0 on Γ with the density ρ 0 . The scaling L △ V should be such that the dynamics κ 0 → κ t preserves this structure, or more precisely, κ t should be of the same type (2.13) The relation (2.13) is known as the chaos propagation property of the Vlasov hierarchy. Under certain assumptions on the mortality m and the kernels a ± , the density ρ t corresponding to the spatial ecologic model, the equation (2.13) implies, in general a non-linear differential equation where the initial condition ρ 0 is a bounded function. Equation (2.14) is called Vlasov-type kinetic equation for ρ t , see [FKK10], [FKKK15] and references therein for more details and [Sor06] for important applications of this model in various areas of science.
In general, if one does not start with a Poisson measure, the solution will leave the space M 1 (Γ). To have a bigger class of initial measures, we may consider the cone inside M 1 (Γ) generated by convex combinations of Poisson measures, denoted by P(Γ).
Remark 1. Below we discuss the concept of a fractional Fokker-Plank equation and the related fractional statistical dynamics, which is still an evolution in the space of probability measures M 1 (Γ) on the configuration space Γ. The mesoscopic scaling of this evolutions leads to a fractional kinetic FPE. The subordination principle provides the representation of the solution to this equation as a flow of measures that is a transformation of a Poisson flow for the initial kinetic FPE, see Sections 3 and 4 below.
The inverse subordinator E = {E t , t ≥ 0} (also called the first hitting time process for the subordinator S) is defined by and its density we denote by ̺ t (τ ), that is Then the subordination process Y t := X Et , t ≥ 0 is such that the onedimensional distribution ν t is given by Let k be the kernel defined by K as its Laplace transform With the help of the kernel k we define the general fractional derivative (GFD) developed in [Koc11] which plays a basic role in this paper (2.15) In Subsection 2.3 we study in more details the derivative D (k) t . The natural question about the subordination process Y is: What type of "differential" equation does the distribution ν t of Y t satisfies? The answer is as follows: The distribution ν t satisfies the following GFD equation As a result we shall consider the fractional Fokker-Plank equation with the GFD (2.15) D (k) The corresponding evolutions of the correlation functions for the Vlasov scaling is D which is a non-Markov evolution. We would like to study some properties of the evolution κ t,k . The general subordination principle gives which is a relation to all orders of the correlation functions. The kernel ̺ t and its properties are studied in Section 3 below. In particular, the density of "particles" is given by The general subordination principle (2.16) gives From this representation we should be able to derive an effect of the fractional derivative onto the evolution of the density, see Sections 4 and 5.
Remark 2. Certain heuristic motivations in physics are leading to the following non-linear equation with fractional time derivative: Note that the subordinated density dynamics has no relation with the solution to this equation. Both evolutions will coincide only in the case of a linear operator in the right hand side.
It is reasonable to study the properties of subordinated flows in (2.17) from a more general point of view, when the evolution of densities ρ t (x) is not necessarily related to a particular Vlasov-type kinetic equation, this is realized in Sections 4 and 5 below.

Definitions
Motivated by the preceding considerations, we recall the general concept of fractional derivative developed in [Koc11] which plays a basic role in this paper. The basic ingredient of the theory of evolution equations, [KST06,EIK04] is to consider, instead of the first time derivative, the Caputo-Djrbashian fractional derivative of order α ∈ (0, 1) , t > 0.
Further details on fractional calculus may be found in [KST06] and references therein.
More generally, it is natural to consider differential-convolution operators where k ∈ L 1 loc (R + ) is a nonnegative kernel. As an example of such operator, we consider the distributed order derivative D Remark 3. 1. The Caputo-Djrbashian fractional derivative (2.19) are widely used in physics, see [MK00,MK04,Mai10], for modeling slow relaxation and diffusion processes. In this case the power-like decay of the mean square displacement of a diffusive particle appears instead of the classical exponential decay.
Considering the general operator (2.20), it is natural to investigate the conditions on the kernel k ∈ L 1 loc (R + ) such that the operator D (k) t possess a right inverse (a kind of a fractional integral) and produce, a kind of a fractional derivative, equations of evolution type. In particular, it means that where λ > 0, has a unique solution u λ , infinitely differentiable for t > 0 and completely monotone, u λ ∈ CM, see Appendix A for this and other classes of functions in what follows.
(2.23) where w 0 is a bounded globally Hlder continuous function, that is |w 0 (x) − w 0 (y)| ≤ C|x − y| θ , 0 < θ ≤ 1, for any x, y ∈ R d , has a unique bounded solution. In addition, the equation (2.23) possesses a fundamental solution of the Cauchy problem, a kernel which is a probability density.

When
3. The asymptotic properties of E α for real arguments are given by, see for example [GKMS14].
which resembles the classical case α = 1, E 1 (z) = e z . On the other hand Here and below C denotes a positive constant which changes from line to line. This slow decay property is at the origin of a large variety of applications of fractional differential equations.
A more general choice of the weight function µ leads to other type of decay patterns.
The conditions upon k guaranteeing a solution to (A) and (B) were given in [Koc11]. The sufficient conditions are as follows.

Asymptotic Properties
As the kernel k and its Laplace transform K are among the objects which play a major role in what follows, here we collect some of its asymptotic properties which depends on the kind of fractional derivative considered. Two cases are studied, the distributed order derivative with k given by (2.21) and the general fractional derivative (2.20) for which K is a Stieltjes function.
Distributed order derivatives. The following two propositions refers to the special case of distributed order derivative, we refer to [Koc08b] for the details and proofs . (2.30) Notice that (2.29) implies that k ∈ L 1 ([0, T ]), however k / ∈ L q ([0, T ]) for any q > 1.

Solutions of the Evolution Equations
Let L be a heuristic Markov generator defined on functions u(t, x), t > 0, x ∈ R d . We have in mind the Bolker-Pacala model and the related non-linear equation, see Subsection 2.2 for details. Consider the evolution equations of the following type with the same operator L acting in the spatial variables x with the same initial conditions The solutions of equations (3.1) and (3.2) typically satisfy the subordination principle [Baz00], that is there exists a nonnegative density kernel function ̺ t (s), s, t > 0, such that ∞ 0 ̺ t (s) ds = 1 and The appropriate notions of the solutions of (3.1) and (3.2) depend on the specific setting, they were explained • in [Koc11] for the case where L is the Laplace operator on R n , • in [Baz00, Baz01, Baz15] with abstract semigroup generators for special classes of kernels k, • in [Pr93] for abstract Volterra equations.
There is also a probabilistic interpretation of the subordination identities (see, for example, [Kol11,Sat99]). In the models of statistical dynamics we deal with a subordination of measure flows that will give a weak solution to the corresponding fractional equation.
In the above relation (3.3), the subordination kernel ̺ t (s) does not depend on L and can be found as follows [Koc11]. Example 8 (α-stable subordinator). Let S be a α-stable subordinator, α ∈ (0, 1) with Laplace exponent L(p) = p α and the corresponding Lévy measure In this case K(p) = p α−1 and the kernel k is given by The associated general fractional derivative D t , see (2.19). As for the density ̺ t (τ ), it follows from Corollary 3.1 in [MSK04] that where g α is the density function of S 1 , that is its Laplace transform is given byg α (p) = e −p α . In addition, it was shown in Proposition 1(a) in [Bin71], see also Theorem 4.3 in [BKS96], that E t has a Mittag-Leffler distribution, that is It follows from the asymptotic of the Mittag-Leffler function E α in (2.24) that̺

Subordination of Moving Step Function
Let u 0 (t, x) be a solution to the kinetic evolution equation, say for the Bolker-Pacala model from Subsection 2.2. The subordinated dynamics is given by As a simple example we take a traveling step function (this is a toy example of a traveling wave in fact) u 0 (t, x) = 1 (−∞,0] (x − tv), where x ∈ R, t, v ∈ R + and consider the subordination of the moving step function. In this case We are interested in studying this dynamics, one possibility is to study the Cesaro limit The Cesaro limit (4.2) may be realized in a number of particular cases related to the fractional derivative considered. Below we investigate three cases corresponding to the α-stable subordinator, the distributed order derivative and general fractional derivative with K fulfilling (H).
Remark 9. The Cesaro limit of the initial step function u 0 (·, x), for fixed x, is given by In addition, for the moving x(t) = ct β with β < 1 we obtain the same asymptotic. Note that the assumption β < 1 is needed to ensure that x(t)/t −→ 0 as t → ∞.
In order to study the Cesaro limit M t (u) in (4.2), at first we compute the Laplace transform of u(t, x) in t using Fubini's theorem and the equality (3.7) yields The following three cases are distinguished.
1. α-stable subordinator. It follows from Example 8 that K(p) = p α−1 which, for fixed x ∈ R, implies where L(x), x > 0 is a slowly varying function, that is see also Definition 26-2. It follows from Karamata Tauberian theorem (cf. Theorem 28-(i) with ρ = 1) that t 0 u(τ, x) dτ ∼ tL(t), as t → ∞ and this implies that the Cesaro limit The same results holds if ,instead of a fixed x ∈ R we take the moving x(t) = ct β with 0 < β < α < 1.
2. The distributed order derivative. It follows from the asymptotic at the origin in (2.32) that , p ≥ 0 is a slowly varying function. Then an application of the Karamata Tauberian theorem (cf. Theorem 28-(i) with ρ = 1) yields which implies the Cesaro limit M t (u) −→ 1 as t → ∞, for any fixed x ≥ 0. For the moving x(t) = ct β , for any β > 0, we obtain Note that the motion of the point x in time with any positive power turns the Cesaro limit to vanish.
3. General fractional derivative. If the measure σ in (2.33) is absolutely continuous with respect to the Lebesgue measure and the density ϕ satisfies (2.36), then Lemma 7 implies that K(p) ∼ Cp θ−1 , as p → 0 andũ and L(x), x > 0 is a slowly varying function. Once again by Karamata's Tauberian theorem (cf. Theorem 28-(i) with ρ = 1) we obtain the asymptotic for M t (u), namely and again we have M t (u) −→ 1 as t → ∞. For the moving x(t) = ct β for any β > θ, we obtain The motion of the point x in time with a positive power β such that β > θ turns the Cesaro limit to vanish.

Traveling Waves
Now we would like to consider a realistic dynamics u 0 (t, x) which is presented by a traveling wave for the non-local spatial logistic equation. This evolution equation appeared as the kinetic equation in the Bolker-Pacala ecological model, see Subsection 2.2 and [BP97, FM04, FKK09, FKK10, FKK12, FKKK15] for more details.
To avoid certain technical details, we will assume the following concrete relations between the mortality m, competition a − and dispersion a + kernels on the generator L (2.8), see [FKT18] for more details.
A traveling wave u(t, x), t ≥ 0, x ∈ R with velocity v > 0 is defined by a profile function ψ : R −→ [0, 1], that is a continuous monotonically decreasing function such that For a fixed x ∈ R the traveling wave u 0 (t, x) = ψ(x − vt) as a function of t is monotonically increasing and has the following properties: As a characteristic of this dynamics, we again consider the Cesaro limit It follows from (5.1)-5.3 the following upper bound for M t (u 0 ) On the other hand, a bound from below of M t (u 0 ) is obtained by Putting together, we have due to arbitrary δ > 0, the Cesaro limit lim t→∞ M t (u) = 1.
Remark 11. Note that for the moving x(t) = ct β , 0 < β < 1 the asymptotic for M t (u 0 ) will be the same. Now we will consider the subordinated dynamics and study the Cesaro limit To this end, at first we rewrite u(t, x) as the sum of three terms. Denoting (5.6) We study each term separately.
I 1 (t, x): It follows from (5.1) that The asymptotic for the integral on the rhs follows as in Section 4 (with ζ − instead of x v ) for all three cases of α-stable subordinator (corresponding to Caputo-Djrbashian fractional derivatives) distributed order derivative, general fractional derivative, we have I 2 (t, x) : In order to study the behavior in t of The Laplace transform of f is equal tõ where g 1 (p) = 1 p e −ζ + pK(p) , g 2 (p) = 1 p e −ζ − pK(p) .
In all the three cases of fractional derivative we have considered (e.g., for the general fractional derivative we use our hypothesis (2.27)) we have pK(p) → 0 as p → 0, therefore e −ζ + pK(p) → 1 and e −ζ − pK(p) → 1 as p → 0. An application of Karamata's Tauberian theorem, see Theorem 28, yields Therefore, for the term I 2 (t, x), we have I 3 (t, x) : Finally we investigate the Cesaro limit of I 3 (t, x), that is It follows from (5.2) the estimates was studied in Section 4 with ζ + = x /v and its Cesaro limit, for all the three types of fractional derivatives considered in Subsection 2.3, was shown to be Hence, we have From the arbitrary of δ > 0 we obtain Putting all together, the Cesaro limit for the subordinated dynamics by the density ̺ t (τ ) gives This result is true for the three type of fractional derivatives considered in Subsection 2.3.

A Bernstein, Complete Bernstein and Stieltjes Functions
In this appendix we collect certain notions of functions theory needed throughout the paper. Namely, the classes of completely monotone, Stieltjes, Bernstein functions and complete Bernstein functions. They are used in connection with the properties of the Laplace transform (LT). More details on these classes may be found in [SSV12]. The family of all completely monotone functions will be denoted by CM.
The function [0, ∞) ∋ τ → e −τ t , 0 ≤ t < ∞ is a prime example of a completely monotone function. In fact, any element ϕ ∈ CM can be written as an integral mixture of this family. This is precisely the contents of the next theorem, due to Bernstein, on the characterization of the class CM in 2. Conversely, wheneverφ(p) < ∞, ∀p > 0, the function [0, ∞) ∋ p → ϕ(p) is completely monotone, that is ϕ belongs to the class CM.
Remark 14. The class CM of completely monotone functions is easily seen to be closed under pointwise addition, multiplication and convergence. However, the composition of elements of the class CM is, in general, not completely monotone.
Stieltjes functions. A subclass of the completely monotone functions is the, so called Stieltjes functions, and they play a central role in the study of complete Bernstein functions, defined below.
where a, b ≥ 0 and σ is a Borel measure on (0, ∞) such that The family of all Stieltjes functions we denote by S.
is a completely monotone function whose LTh(p) exists for any p > 0.
In particular, we see that S ⊂ CM and S consists of all ϕ ∈ CM such that its representation measure (from Theorem 13) has a completely monotone density on (0, ∞), for ϕ ∈ S is of the form Example 17. The following are examples of Stieltjes functions for any τ, t > 0 ϕ 1 (τ ) = 1, ϕ 2 (τ ) = 1 τ , ϕ 3 (τ ) = (τ + t) −1 , ϕ 4 (τ ) = 1 + t τ + t , Bernstein functions. Now we introduce the class of Bernstein functions which are closely related to completely monotone function. Bernstein functions are also known in probabilistic terms as Laplace exponents. The Lévy triplet (a, b, µ) determines ϕ uniquely and vice versa. In particular, 3. The class of Bernstein function will be denoted by BF .
The following structural characterization theorem of Bernstein functions is due to Bochner, see [SSV12,Thm 3.7] for the proof.
The class of complete Bernstein functions we denote by CBF .
2. A function ϕ ≡ 0 is a complete Bernstein function if, and only if, 1/ϕ is a non-trivial Stieltjes function.

S • S ⊂ CBF .
We conclude this subsection with some examples of elements in the class CBF .

B The Karamata Tauberian Theorem
Tauberian theorems deals with the deduction of the asymptotic behavior of functions from a certain class (regular varying in the original of Karamata [Kar33]) from the asymptotic behavior of their transforms (e.g. their Laplace-Stieltjes transforms). We refer to [Sen76, Sec. 2.2] and [BGT87] for more details and proofs. Let A > 0 be given and denote by F + (A) the class of positive measurable functions defined on [A, ∞). Then, if ρ ≥ 0 and L is a slowly varying function x /Γ(ρ + 1) as x → 0 + .