Heat kernel asymptotics of the subordinator and subordinate Brownian motion

For a class of Laplace exponents, we consider the transition density of the subordinator and the heat kernel of the corresponding subordinate Brownian motion.We derive explicit approximate expressions for these objects in the form of asymptotic expansions: via the saddle point method for the subordinator’s transition density and via the Mellin transform for the subordinate heat kernel. The latter builds on ideas from index theory using zeta functions. In either case, we highlight the role played by the analyticity of the Laplace exponent for the qualitative properties of the asymptotics.

theoretical physics. The heat kernel appears typically as the fundamental solution of a partial differential equation, as the integral kernel of an operator semigroup or as the transition density of a stochastic process.
In its most basic form, the heat kernel refers to the Laplace operator on Euclidean space or the Laplace-Beltrami operator on a Riemannian manifold. Classical results link the short-time asymptotics of the heat kernel on a closed manifold to the geometry of the manifold [33], and it plays a significant role in index theory [2,5].
The heat kernel can be computed explicitly only in special cases, usually in the presence of high degrees of symmetry of the underlying manifold. We refer the reader to [10] for examples of techniques for explicitly constructing heat kernels. In all other cases, one must content oneself with bounds or (asymptotic) approximations of the heat kernel.
Recent attention has moved towards the heat kernels of nonlocal operators of the form f ( ) for suitable functions f . These operators also naturally appear in probability theory as the infinitesimal generators of subordinate Brownian motion with Laplace exponent f . Typical examples would be f (z) = z α or f (z) = (1 + z) α − 1. One is interested in the properties of the associated heat kernels on Euclidean space or Riemannian manifolds with various types of boundary conditions. Such operators and the heat kernels are also important from a practical point of view since they appear naturally in physics [4] or financial mathematics [16].
This paper continues the theme of [21] albeit on Euclidean space. To summarise our approach and results in a non-technical way, denote by B t a standard Brownian motion on R n and let X t be a subordinator with corresponding Laplace exponent f , i.e., X t is an almost surely increasing Lévy process that takes values in the nonnegative reals. It can be thought of as introducing an "operating time". Our key assumption is that the subordinators possess a Lévy density which has an asymptotic power series near the origin and is of rapid decay at infinity. We then derive asymptotics in t for both the transition density p t (τ ) of the subordinator X t and the heat kernel k t (x, y) of the subordinate Brownian motion B X t . Qualitatively, our results and related results can be summarised in Fig. 1 (note that we recast cited results in our notation, some authors write p t (x) instead of p t (τ ) or p t (x − y) instead of k t (x, y) for symmetric processes).
The asymptotics of p t (τ ) are obtained via an application of the saddle point method to an integral representation of the transition density. The difficulty here is that the saddle point depends on t and τ so that a detailed analysis is required. A key role is played by the analyticity of the Laplace exponent on a cut plane. The saddle point methods shows that the transition density decreases exponentially for t → ∞ and τ → 0.
On the other hand, the asymptotics of the subordinate heat kernel are obtained using ideas from index theory: we introduce a function ("zeta function") that is the Mellin transform of the heat kernel. Standard properties of the Mellin transform together with the fact that f is not entire then show that for small values of t, the heat kernel possesses an asymptotic power series in t. This is contrast to the rapid decay of k t which is observed for f entire. We give several characterisations of the zeta function and derive explicit lowest order terms for our class of Bernstein functions.
Closely related to our investigation of the subordinator's transition density is [31] (see also [32] for similar estimates in higher-dimensional Euclidean space). The authors employ the saddle point method to derive asymptotics of the transition densities of certain one-dimensional Lévy processes. Under the key assumption of the characteristic exponent (equivalently the symbol of the process) being an entire function, they show that p t (τ ) has exponential decay as t +τ → ∞. In a sense, this complements our analysis since the authors consider the "opposite" direction of the τ -variable. Also, the assumption of an everywhere analytic symbol leads to exponential decay. Note that in general the characteristic exponent of subordinate Brownian motion is not entire as it is given in terms of a Bernstein function. Except for special cases, these functions are analytic on a half-plane or on a cut plane. It is precisely the lack of analyticity that leads to an asymptotic power series in t for the heat kernel of subordinate Brownian motion as opposed to rapid decay.
Another closely related paper is the recent [34] in which the author uses scaling properties of the Laplace exponent f to derive upper and lower heat kernel bounds for subordinate Brownian motion. The estimates are expressed in terms of the inverse f −1 and an auxiliary function that also appears in our paper. The class of Laplace exponents considered covers our class.
An efficient and widely used method for deriving heat kernel estimates was developed in [11] that links Nash's inequality and Dirichlet forms. The connection with stochastic processes on Euclidean space was exploited for example in [3,12] or [13,15] and related papers by the same authors. Their approach relies on growth properties of the jumping intensity and its impact on the associated Dirichlet form. For the type of Laplace exponent considered by us, we recover parts of the cited results in the short-term asymptotics. Analytic methods to obtain heat kernel bounds are described in [18,23]. We do not even attempt to summarise the recent literature for the situation on bounded domains. (We do not claim completeness of this list and refer the reader to the references in these papers.) This paper is organised as follows. The subsequent section introduces some notation and collects some preliminary material. We state and comment on the key results in Sect. 3. These are proved in the subsequent two sections followed by a short appendix that calculates a particular line integral.

Preliminaries
We introduce some notation and collect various prerequisites. Landau O-notation Let f and g be two functions (0, ∞) → C. Then we characterise the growth of this function as follows.
Asymptotic series This definition closely follows Chapter 2 in [17]. Let g : (0, ∞) → R be a function whose asymptotic behaviour at 0 we wish to characterise. (Analogous definitions apply for the asymptotics at ∞.) A sequence of functions {ϕ k } is called an asymptotic sequence at 0 if ϕ k+1 = o(ϕ k ) at 0. We then say that for every N ≥ 0 as x → 0 + . Typical choices of ϕ k we will encounter are ϕ Slowly varying functions We say that a function f : (0, ∞) → C is slowly varying at ∞ if it is nonzero for large enough arguments and f (λz)/ f (z) → 1 as z → ∞ for all λ > 0, cf. [43] and Appendix 1 of [7]. Likewise for z → 0.
Bernstein functions Recall (cf. [37,Definition 3.1]) that a function f : Any Bernstein function can be represented in Lévy-Khintchin form as for constants a, b ≥ 0 and μ a measure on (0, ∞) such that ∞ 0 t ∧ 1 μ(dt) < ∞. This means in particular that it is smooth on (0, ∞), can be extended to an analytic function on the half-plane {z ∈ C|Re z > σ 0 } for some σ 0 ≤ 0 and is continuous on the axis σ 0 + iR. If the measure μ has a density m with respect to Lebesgue measure, we call m the Lévy density.
The Mellin transform Given a locally integrable function f : dt for any z ∈ C for which the integral converges. The transform exchanges growth of f at 0 and ∞ with complex differentiability in the following sense, cf. [8,Chapter 4] for details. For our purposes, it suffices to consider the case where f has the following asymptotic expansions: Then the Mellin transform of f will be analytic on the strip α < Re z < β where the boundaries are determined explicitly: In the cases of finite α, β one can meromorphically extend the Mellin transform with simple poles that are given in terms of the asymptotic expansions of f , cf. [ Correspondingly, if q = 0, the extension holds in the half-plane Re z < α with at most simple poles at z = − a m and residue p m . The converse of these claims also holds, cf. [24]. Of particular importance will be the Plancherel formula which reads with c ∈ R in the intersection of the strips of analyticity of the transforms, cf. [41].

Key results
We state the assumptions, key results and comment on them.
3.1. The transition density of the subordinator Let X t be a subordinator, i.e., a càdlàg stochastic process taking only nonnegative values, X 0 = 0 a.s. and having stationary and independent increments. The reader is referred to [1,6,37] for further details of these processes. The characteristic function of X t can be expressed as for a negative-definite function ψ, the characteristic exponent or alternatively by its moment generating function with Laplace exponent f . The relationship between ψ and f is ψ(λ) = f (−iλ). The function f is a Bernstein function.
We assume that the subordinator X t has a transition density p t (·). The characteristic function is the Fourier transform of p, i.e., The density can be recovered by Fourier inversion as or by inverting the moment generating function using the Bromwich integral Here c ∈ R is chosen so that all singularities of the integrand are to the left of the vertical axis c + iR. cf. [41]. The integral can only be evaluated in closed form in special cases and hence must be approximated in general. We expect the growth of f at ∞ to determine the behaviour of the transition density near τ = 0. To illustrate some possible situations, we consider two examples. Although the first example is discussed in the literature, for example [39,Proposition 5.29], we sketch an alternative proof that sheds light on the relationship between the growth of f at ∞ and the growth of the transition density at τ = 0. EXAMPLE 3.1. Let f (z) = log(1 + z α ) for α ∈ (0, 1). Then for t > 0 fixed and Proof. ( with a complex logarithm defined on the cut plane C\(− ∞, 0] and arguments in (− π, π]. 2. Now compute the Mellin transforms. The Mellin transform of e ±iτ λ is explicitly given in terms of the Gamma function as M[e ±iτ λ ; z] = τ −z e ±πiz/2 (z). This is originally defined for 0 < Re z < 1 but has an obvious meromorphic extension to the whole of C. The Mellin transform of (1 + (±iλ) α ) −t could also be computed explicitly in terms of the Beta function, cf. [19, equation 6.2.(30)]. However, to illustrate the general case we consider the growth/decay of (1 + (±iλ) α ) −t as λ → 0 and λ → ∞ since that yields the desired properties of the Mellin transform. Indeed, we have where the dots indicate terms of higher (lower) order in λ. By standard properties of the Mellin transform, M[(1+(±iλ) α ) −t ; z] is analytic in the strip 0 < Re z < αt and can be meromorphically extended to the whole complex plane. It has a simple pole at z = αt with residue (±i) −αt = e −αt (±πi/2) . 3. The Plancherel formula for the Mellin transform then yields where 0 < c < αt. 4. To obtain the asymptotic behaviour of p t (τ ) for τ → 0 we now apply Cauchy's theorem and move the contour to the right across the poles of the integrands. Since αt < 1, the first pole we encounter is at αt and it is due to M[(1 + (±iλ) α ) −t ; z]. We find where the last line follows from the reflection formula of the Gamma function. The remainder term is a line integral with powers in τ strictly greater than αt −1, cf.
Step 3 of the proof of Proposition 4.1. One can repeat this and move the contour across further poles. Each pole yields a term in the asymptotic expansion of p t (τ ) in powers of τ .
We contrast this with the case when the Laplace exponent grows like z α instead of log z. It is clear that the method just outlined will not yield a meaningful asymptotic expansion since e −t f (ξ ) is of exponential decay so its Mellin transform will have no poles in the right half-plane. The formally obtained asymptotic power series in τ is zero and does not yield much detailed information about the decay of the transition density.
τ 3/2 e −t 2 /4τ +t−τ , as can be seen by evaluating the Fourier integral (3). In particular, this is of rapid decay for τ → 0 with t fixed (and t → ∞ with τ fixed).
We are thus led to approximate the Bromwich integral (4) by the saddle point method (cf. for example [17,Section 8]) since we expect this method to extract an exponential term in the asymptotic expansion of p t (τ ). This requires us to find points where the derivative of the exponent vanishes, which clearly happens at the point Standard arguments using the Cauchy-Riemann equations show that this is a saddle point of the integrand. The saddle point is a function of τ and t so that it moves when the parameters change, which necessitates a careful analysis.
3. Since f is a Bernstein function, Eq. (5) has a unique solution for τ/t real. Indeed, the inverse function is smooth also. In general, there are further complexvalued solutions to (5) with each leading to a saddle point that should be considered in the approximation. We cannot expect to say in the general case that the real-valued solution is the "dominant" saddle point in the sense that the approximations due to the other solutions can be neglected. It turns out, however, that the real-valued solution gives a meaningful approximation of the density with explicit error bounds.
We shift the integration contour so that it goes through the saddle point. Ideally, we want to find a contour of steepest descent, i.e., one on which the imaginary part of the exponent in the Bromwich integral is constant. However, to paraphrase [26, Chapter 3.3], it is not necessary in practice to find the path of steepest descent since any path descending from the saddle point will yield the correct answer.
The easiest way is to integrate along the vertical line ξ s + iR. Changing coordinates and defining ξ = ξ s + iη, we find The saddle point method suggests that we approximate p t by the Gaussian integral We will make this precise and also address the difficulty that the saddle point is not fixed but depends on τ/t. We are now ready for the key assumption on the Laplace exponent. HYPOTHESIS 3.4. We assume that f is a complete Bernstein function with Lévy density m that has the following properties.
(i) Behaviour at 0. The density m is defined on a neighbourhood of the origin in C.
Moreover, it has an asymptotic expansion of the form for λ → 0 that is valid along all rays arg λ = const. Here, a j ∈ R, and exponents satisfy α 0 = −α − 1 for some α ∈ (0, 1) and α j ↑ ∞. (ii) Behaviour at ∞. The density decays exponentially fast in the sense that the exponent is strictly negative.
REMARK 3.5. Some remarks on this assumption and its consequences.
(i) The requirement that the asymptotic expansion of m be valid in a neighbourhood of the origin in the complex plane leads to uniform growth estimates of the Bernstein function in all directions in the complex plane, cf. Lemma 4.2. (ii) Lemma 4.2 also shows that the Bernstein function grows like z α so that the integrals in (3) and (4) converge absolutely. Necessarily, the transition density of the subordinator must be of this form so that we see its existence. The reader is also referred to [25] for sufficient conditions ensuring the existence of a density. (iii) Note that any complete Bernstein function can be analytically continued to the cut complex plane C\(− ∞, 0] by [37, Theorem 6.2 (v)]. The decay of m at infinity ensures that the function f is analytic on the cut plane C\(− ∞, σ 0 ]. (iv) In a sense, Hypothesis 3.4 is a generalisation of Hypothesis 1 in [20], where we required a j = − α + ( j − 1) so that f ( ) is a classical pseudodifferential operator. This required the constant "spacing" of the exponents. Here, we relax this assumption so that we allow any asymptotic expansion with increasing exponents. In particular, all of the examples exhibited in [20] satisfy the above hypothesis. (v) The parameter α has a probabilistic interpretation as the Blumenthal-Getoor index of the subordinate Brownian motion whose subordinator is defined in terms of the Bernstein function f satisfying Hypothesis 3.4. This is explained in Theorem 1 of [20]. We use the method of steepest descents to derive an approximation of the transition density that is asymptotically valid for small τ and large t.
. REMARK 3.7. Two brief remarks on the theorem.
(i) We note that the exponents in the above theorem are related to [28] and [34]. For certain stochastic processes, the authors derive bounds on the transition density in terms of the auxiliary function so that H is crucial for the estimates in Theorem 3.6, too. (ii) Under the assumption of the characteristic exponent of certain processes being entire, the authors of [31] show that the corresponding transition density p t (τ ) is of exponential decay as t + τ → ∞. They also use the saddle point method to identify the exponential term that causes the decay. Their result complements our analysis in the senses that they consider the case of large τ .
An important corollary concerns the exponential decay of the transition density.
In particular, the transition density decays exponentially as t → ∞ or τ → 0. EXAMPLE 3.9. For the relativistic α-stable subordinator we have f (z) = (1 + z) α − 1 with σ 0 = 1. Neglecting higher-order terms, this leads to the approximation of the transition density given by which equals the exact transition density if α = 1/2.

The heat kernel of subordinate Brownian motion
The key idea is to use the Mellin transform to translate the short-time asymptotic expansion of the heat kernel into the pole structure of a suitable meromorphic function that is the product of the Gamma function and a "zeta function" to be defined. We mimic the construction of the operator zeta function from index theory (see for example [22,Chapter 1.12]) but dispense with some of the restrictions inherent in the calculus of pseudodifferential operators.
We denote the heat kernel of the subordinate Brownian motion with Laplace exponent f by k t (x, y) if there is no confusion about f , and if we want to emphasize the dependence on the Bernstein function we write k[ f ](t; x, y). This may seem cumbersome but the reader will see immediately why this is sensible. The heat kernel is given by Here, h τ is the heat kernel of a standard Brownian motion and p is the transition density of the subordinator from the Bromwich integral (4).
so that ζ is the Mellin transform of k[ f ] with respect to t up to the factor 1/ (s).
REMARK 3.11. We comment on the motivation and suitability of this definition.
(i) The definition is motivated by the classical operator zeta function based on the relation If we are given an operator A acting on a Hilbert space with an orthonormal basis of eigenfunctions of A and eigenvalues λ n , then formally Recalling that λ −s n = Trace (A −s ) = ζ(s), the operator zeta function of A, and e −λ n t = Trace (e −At ) the heat trace, we find Moving from the trace of the complex powers A −s to their integral kernel yields the above definition of the zeta function ζ(s; x, y). We refer the reader also to [21] where this is done in the context of pseudodifferential operators and ζ(s; x, y) is indeed the integral kernel of the complex powers A −s . (ii) It is not obvious that Definition 3.10 makes sense as the integral may not converge. Indeed, for general f we cannot expect (7) to exist other than in the sense of an oscillatory integral, cf. Chapter 7.8 of [27] or Chapter I.1.2 of [38]. However, in our class of Bernstein functions if we choose a > 0 suitably large, then the integral for ζ [ f + a](s; x, y) converges absolutely as shown in Lemma 5.1.
Here, f + a stands for the Bernstein function z → f (z) + a. The reason why we are free to choose such an a is as follows: From the integrals (4) and (6), we see that which also justifies that we mention the Bernstein function in the notation. As we are interested in the short-time asymptotics of the heat kernel k[ f ], we may as well derive the asymptotics of k[ f + a] and scale by e at .
We characterise the zeta function from three perspectives. (i) Choose c a real number such that σ 0 < c < 0 and fix an a > 0 such that The zeta function can then be written as where ν = n/2 − 1 and . Here, K ν is a modified Bessel function of the second kind. (ii) We can also express the zeta function as where g s (r ) = f (σ 0 + re −iπ ) + a −s − f (σ 0 + re iπ ) + a −s and a is as in (i).

(iii) Denote by u s the solution of the Poisson equation on R n (Euclidean Klein-Gordon equation)
Here, 1 A is the indicator function for the set A ⊆ R n , = n i=1 ∂ 2 x i the Laplacian on R n , and vol(S n−1 ) denotes the volume of the unit sphere in R n . Then ζ(s; x, y) = u s (0). REMARK 3.13. A few remarks on Theorem 3.12 and Corollary 3.14.
(i) For the representation (8), we have the following comments.
(a) The integral is defined based on two cuts in the complex plane that allow the construction of a complex logarithm: To define [ f (ξ ) + a] −s we cut along (− ∞, σ 0 ], with respect to this cut, we allow the argument of a complex number to be in the interval (− π, π]. To define the square root (−ξ) 1/2 and its powers, we cut along [0, ∞) and with respect to this cut, we allow the argument of a complex number to be in the interval [0, 2π). The assumption σ 0 < 0 is crucial for our argument as it allows the contour c + iR to pass between the two cuts. If σ 0 were zero, this would not be possible. (b) If f were entire, then the line integral would vanish by Cauchy's theorem so that the zeta function would be identically zero. The consequence would be that all heat kernel coefficients vanish so that the heat kernel decays faster than any polynomial as t → 0. Hence, the polynomial nature of k t for small t is due to f not being entire. This is the case considered in Corollary 8 of [32] where an entire characteristic exponent leads to exponential decay of the heat kernel as t → 0. We refer to Remark 3.7 for a comparison with [31] where the assumption of entire characteristic exponents is important and leads to exponential decay of the transition density (albeit for t → ∞ and large distances d(x, y), but our analysis again complements this.) We also refer the reader to a related discussion of the decay of the transition density and the analyticity of the Laplace exponent in [9,Chapter 3] in the context of modelling real-world phenomena. (ii) We can view the representation (9) as a Bessel or Hankel transform. (iii) The representation in (3.12) is useful for numerical computations of the heat kernel coefficients. Moreover, the generator of the subordinate Brownian motion is the operator f ( ) in the sense of a Dunford integral (cf. the discussion in Section 4 of [36]). The Dunford integral requires the full resolvent family, i.e., we must know (− + λ) −1 for any λ > 0. However, the heat kernel coefficients for given x, y as expressed by ζ(s; x, y) need only a single resolvent, namely d(x, y)) −1 . When expanding the heat kernel of the nonlocal operator f ( ) in powers of t, we obtain coefficients that are determined in terms of the local (namely differential) operator .
The relation with the short-time asymptotics of the subordinate heat kernel is an immediate consequence.
for t → 0. This is independent of the choice of a.
We will use asymptotic methods analysis to evaluate this integral, of course exploiting the asymptotic expansion of m near the origin. The integral should, however, also serve as the basis in more general situations when the restrictions placed on f (and its density m) are less severe; this is the subject of current research.
The aim is now to find an approximation to the first heat kernel coefficient. We will do this in the form of an asymptotic series that is valid for small values of t and d.
as t → 0 in lowest orders in t and the Euclidean distance d.
which is obtained by setting V (r ) = r n and φ(r ) = r 2α in the notation of [15].
An immediate consequence concerns the heat kernel expansion of the relativistic α-stable process.
as t → 0 in lowest orders in t and the Euclidean distance d. It is instructive to compare this result with the literature: it agrees with [21] translated to the case of a flat manifold. Moreover, the bounds for the relativistic stable process (α = 1/2) obtained in [ for constants c 1 , c 2 and t ∈ (0, t 0 ] for t 0 fixed. This obviously agrees with our result for sufficiently small t.

Proof of the approximate transition density of the subordinator
In this section, we provide the proof of the approximation result in Theorem 3.6 and of some auxiliary results. We first collect growth properties of our class of Bernstein functions in a separate subsection and then prove the main result.

A class of Bernstein functions
We collect growth properties of the Bernstein functions satisfying Hypothesis 3.4.
The key property of our class Bernstein functions concerns the growth along parallels to the imaginary axis. This will be applied in the proofs of Proposition 4.11 and Lemma 5.1. (i) If x ≥ 0 and y ∈ R, then we can write Here, R is a function such that |y| −α R(x, y) → 0 as |y| → ∞ uniformly in x.
We note that this result is related to the familiar property of complete Bernstein functions preserving sectors in the complex plane, cf. [37, Corollary 6.6].
Proof. The asymptotics in y could be shown by a simple application of Watson's Lemma. However, we need a bound on the remainder term and its dependence on x so that a more detailed analysis is required.
We first prove assertion (i).
1. We use the representation (2) to write by the Plancherel formula for the Mellin transform for c ∈ R to be determined. with c is such that max(−α 1 , α) < c < 1 + α. The next pole is at −α 1 so that c > −α 1 would be sufficient here but we will require c > α in Step 5 below. 4. It remains to estimate the remainder term independently of x and show that |y| −α R(x, y) → 0 as |y| → ∞. Indeed, If we can show that the Mellin transform is bounded independently of x, then we are done: the integral of the Gamma function along parallels to the imaginary axis converges by the decay estimates in [ We note that the integral from 1 to ∞ is an entire function of z due to the rapid decay of m at ∞. Moreover, we can simply bound this independently of x: since x ≥ 0.
The integral from 0 to 1 is rewritten as By the asymptotic expansion of m, the first integral converges absolutely for z = c +iu with −α 1 < c < 1 + α. We bound e −λx above by 1 to obtain a bound independent of x.
In the second integral, we integrate by parts to obtain the desired analytic continuation: The integral term converges absolutely for c > α. The terms can clearly be bounded independently of x so that the proposition is proved.
To deduce assertion (ii) on the pointwise limit of R, we note that the integral on the left hand side of (13) can be bounded for x > σ 0 by the rapid decay of the density m. Recall that the Bernstein function has abscissa of convergence σ 0 so that multiplication by e −λx with x ∈ (− σ 0 , 0) does not affect the convergence of (13). Moreover, the terms in (14) can be bounded in terms of x, albeit not uniformly.
The second description of the growth of the Bernstein functions follows from Watson's Lemma where the key point is that the asymptotic expansion holds in all directions in which z can approach ∞ in the complex plane.
as z → ∞ along any ray in the complex plane.
Proof. This is simply the generalised Watson's Lemma of [45,Chapter I.5] applied to (2). The asymptotic property irrespective of the direction follows since the Lévy density m has the same asymptotic expansion when approaching 0 in any direction.
In the analysis of the transition density of the subordinator, we also need an alternative characterisation of the growth of the Bernstein function. We phrase this in terms of functions of slow variation as the most natural setting, although this is somewhat more general than we need here (all function will have limits so are naturally slowly varying).
for Re z > 0 where l 0 is a slowly varying function at infinity with lim z→∞ l 0 (z) = 1.
Proof. Using the asymptotics of m at 0 and ∞ we decompose the Lévy density additively as m(λ) = λ −1−α e σ 0 λ + n(λ) for some integrable function n(λ). For Re z > 0 we have by (2) that where g(z) = ∞ 0 (1 − e −λz )n(λ)dλ. Directly evaluating the first integral leads to We can thus write Based on the asymptotics of n at 0, Watson's Lemma shows that l 0 (w) → 1 as w → ∞ in the right half-plane. Slow variation is immediate.
This assumption also implies the growth of the derivatives of f which is key for the saddle point method.
Proof. This follows from the arguments leading to equation A1.3 of [7] exploiting the analyticity of the Bernstein function f in the half-plane Re z > σ 0 .
Another auxiliary result concerns the inverse of f on the real line which defines the saddle points, cf. Remark 3.3.
Another lemma concerns the composition of f and its derivatives with f −1 .
Proof. The claims follow by direct calculations using Lemmas 4.4 and 4.5.

Proof of the approximation of the transition density
The proof of Theorem 3.6 depends on a series of propositions that are given subsequently. Recall that we defined ϕ(η) = τ (ξ s + iη) − t f (ξ s + iη).
Proof of Theorem 3.6. The transition density is given as We split the integral into four terms and treat each summand separately. Define η 0 = τ β/(1−α) t γ /(1−α) for β ≤ 0 and γ ≥ 0 to be determined later. Then where we employed a Taylor expansion of ϕ around η = 0 with remainder term R, see Proposition 4.10. This can be further decomposed as .
Each proposition entails a growth estimate in τ and t so that we can write Here, c α is a positive constant depending only on α. By Lemma 4.6 and the definition of η 0 we have In assertion (i) we fix t. We want the first term in the brackets to dominate all other terms in the limit τ → 0.
In assertion (ii), we fix τ and set β = 0. We want the first term in the brackets to dominate all other terms in the limit t → ∞.
This completes the proof.
The first proposition merely recalls a standard result.
Proof. This follows from a direct calculation using ϕ (0) = t f (ξ s ) and the Gaussian integral The second proposition estimates a truncated Gaussian integral in terms the integration limits. PROPOSITION 4.9. We have the bound for any η 0 ≥ 0.
after a change of variables. Now write this as and consider the integral.

In the integral we make a change of variables
Note that f (ξ s ) < 0 since f is a Bernstein function so that e t f (ξ s )η 0 ν ≤ 1. This means that the integral is bounded above by giving the desired result after dividing by 2π .
The third integral entails a Taylor remainder term of ϕ on a finite range.
for any t > t 0 . Here, l 7 is the function of slow variation at 0 from Lemma 4.6.
The analogous result holds if we are given t: there is a τ 0 (t) such that the inequality holds for all τ < τ 0 . Proof.

5.
Overall we obtain the bound 1 2π The final proposition covers the most intricate integral.
1. We consider the integral for on [η 0 , ∞) in detail, the integral on (− ∞, η 0 ] being treated analogously. Using the definition of ϕ we rewrite the integral as Since ξ s → ∞ as τ/t → 0, we have ξ s > 0 for τ/t suitably small. From Proposition 4.1, we thus have where R is a function such that |η| −α R(ξ s , η) → 0 as |η| → ∞ uniformly in ξ s . Thus, we have the bound with c α = − a 0 (−α) cos( απ 2 ) which is positive. Recall that by its definition, η 0 can be made as large as one likes by choosing τ/t sufficiently small. Since η −α R(ξ s , η) tends to zero as η → ∞ uniformly in ξ s , there is an η 0 (and hence a suitable combination of τ and t) such that for η > η 0 , we have |1 + c −1 α η −α R(ξ s , η)| < 1/2. This means that for τ suitably small (or t suitably large), we obtain the bound which one can make explicit. 3. In the last integral, we make the change in variables u = 1 2 c α tη α which leads to The integral on the right hand side is the incomplete Gamma function typically denoted by (s, x) = ∞ x u s−1 e −u du. For the limit x → ∞ one has the asymptotics . In our situation, s = 1/α and x = 1 2 c α tη α 0 we have proving the claim upon collecting powers in t.
Proof of Corollary 3.8. We must compute the function e ϕ(0) / − 2π t f (ξ s ). For the exponent, we find by Corollary 4.7 that for l 8 a slowly varying function at 0 with lim x→0 l 8 (x) = 1. By Lemma 4.6 (ii) we compute the denominator as Combining the exponent and denominator, we arrive at the claim.

Proof of the heat kernel asymptotics
This section contains the proofs of the different characterisations of the zeta function and the general heat kernel asymptotics.

Definition of the zeta function and relation with the heat kernel
We first show that the definition of the zeta function make sense using the growth estimates from Proposition 4.1.
As a consequence, we may change the order of integration in (15) so that in particular the definition of the zeta function in (7) makes sense. REMARK 5.2. Note that given c such an a exists as Proposition 4.1 (iii) guarantees that the function η → Re f (c + iη) is eventually positive when |η| becomes large so that it can only be negative on a bounded set. By continuity this function it must assume a minimum value on that set, and this determines a lower bound for a.
Proof. To show that the absolute value of the integrand is integrable note that The right hand is integrable, and we address this in some detail. By the Fubini-Tonelli Theorem, if the triple integral converges for some ordering of the integrals, then it converges for any ordering. We thus consider the integral Now observe the following: • The τ -integration can be done without problems as we assumed c < 0.
• The t-integration is a standard Mellin integral that gives cf. formula 6.3(1) of [19] whenever Re s > 0 as Re f (c + ξ) + a > 0 by assumption.
• Thus we end up considering Now estimate Re f (c + iη) + a in terms of η. By the choice of a, the integrand is everywhere strictly positive. Moreover, by Proposition 4.1 (ii) we can write for a function R such that |η| −α R(c, η) → 0 as |η| → ∞. Here, c α is a positive constant that depends on α. So the η-integral converges if αRe s > 1.
This proves the claim.

Different characterisations of the zeta function
In Theorem 3.12, we gave three characterisations of the zeta function.
(i) Suppose first that Re s > 1/α. Plugging (4) and (6)  3. The extension to all s ∈ C follows from asymptotic properties of the modified Bessel functions of the second kind: recall that K ν (z) ∼ z −ν at z = 0 and K ν (z) ∼ z −1/2 e −z as z → ∞ with | arg(z)| < 3π/2, cf. [44, Chapter 7.23, equation (1)]. Note here that the square root in the argument of K ν ensures that when integrating along a parallel to the imaginary axis, the real part of (−ξ) 1/2 tends to infinity when |ξ | → ∞, so that the exponential decay is effective.
(ii) The expression for the zeta function in (8) is a line integral along a parallel to the imaginary axis. We deform the contour to a keyhole contour that encloses the cut (− ∞, σ 0 ]. This is allowed by the exponential decay of the integrand and the fact that the integrand is holomorphic on the left half-plane minus the cut.
Let > 0. The integral along the keyhole contour can be decomposed into the sum of three integrals along curves A, , B, and C, which are given as follows: • A, covers the "upper" part of the cut from σ 0 to −∞. It is parametrised as σ 0 + re iπ for r ∈ [ , ∞). • B, covers the circular arc parametrised as σ 0 + e iθ with θ ∈ (− π, π]. • C, covers the "lower" part of the cut from −∞ to σ 0 . It is parametrised as σ 0 + re −iπ for r ∈ [ , ∞).
We simplify or estimate the integrals separately.

For
Since f is bounded in a neighbourhood of σ 0 and the other factors in the integrand are analytic there, the integral converges to 0 as → 0. 3. For C, we find that The claim follows once we let → 0.
5. It remains to show that the series for F α,l and G α,l converge absolutely. We present the argument for F α,l , the claim for G α,l is proved similarly. It suffices to estimate the growth of ( n 2 + αl − m). By the reflection formula for the Gamma function (z) (1 − z) = π sin π z we find with z = 1 − n 2 − αl + m that .
The summands are positive and the growth of the Gamma function ensures convergence.
We now address the general case. For our class of subordinators, the asymptotic expansion of the Lévy density at 0 ensures that the relativistic α stable process is the basic building block and gives the lowest order heat kernel coefficients. 1. We decompose the Laplace exponent into a part corresponding to the relativistic α-stable subordinator and remainder terms that can be controlled. Using (2) we write where m(λ) = m(λ) − a 0 λ −1−α e σ 0 λ . Since by construction m is integrable, we can write with b = − a 0 (− σ 0 ) α + ∞ 0 m(λ)dλ andf (z) = ∞ 0 e −λz m(λ)dλ. 2. We now derive an upper bound forf in terms of z. To this end we apply the generalised Watson's Lemma from [45,Chapter I.5]. Since we assumed that the asymptotic expansion of m around 0 is valid in all directions, we have an