Heat kernel for higher-order differential operators and generalized exponential functions

We consider the heat kernel for higher-derivative and nonlocal operators in $d$-dimensional Euclidean space-time and its asymptotic behavior. As a building block for operators of such type, we consider the heat kernel of the minimal operator - generic power of the Laplacian - and show that it is given by the expression essentially different from the conventional exponential Wentzel-Kramers-Brillouin (WKB) ansatz. Rather it is represented by the generalized exponential function (GEF) directly related to what is known in mathematics as the Fox-Wright $\varPsi$-functions and Fox $H$-functions. The structure of its essential singularity in the proper time parameter is different from that of the usual exponential ansatz, which invalidated previous attempts to directly generalize the Schwinger-DeWitt heat kernel technique to higher-derivative operators. In particular, contrary to the conventional exponential decay of the heat kernel in space, we show the oscillatory behavior of GEF for higher-derivative operators. We give several integral representations for the generalized exponential function, find its asymptotics and semiclassical expansion, which turns out to be essentially different for local operators and nonlocal operators of noninteger order. Finally, we briefly discuss further applications of the GEF technique to generic higher-derivative and pseudodifferential operators in curved space-time, which might be critically important for applications of Horava-Lifshitz and other UV renormalizable quantum gravity models.


INTRODUCTION
Physical phenomena in higher derivative and nonlocal field theories are essentially different from conventional local quantum field theory (QFT) with the wave operators of second order in space-time derivatives. There are numerous manifestations of this difference including the problem with unitarity which arises due to higherderivative (Ostrogradsky) ghosts [1], violation of Lorentz invariance and peculiar causality properties of Hořava-Lifshitz gravity models (which are motivated by the attempts to solve this problem in renormalizable quantum gravity), noninteger conformal operator dimensions in conformal field theories and so on. These peculiarities are deeply rooted in mathematical formalism of higher derivative models and one of its fundamental ingredients-the heat kernel-the building block underlying the propagator of the theory in Feynman diagrammatic technique. Now the heat kernel method is one of the most powerful tools in mathematical physics, that has a wide range of applications extending from pure mathematics (spectral geometry) to the analysis of financial markets. Being combined with the background field method in QFT it provides directly in the coordinate space a calculational technique for the quantum effective action, studying renormalizability of field models, their quantum anomalies, critical phenomena, etc. This makes it indispensable for computations in the presence of external * barvin@td.lpi.ru † petr@phys.msu.ru ‡ vladvakh@gmail.com fields or in curved space-time, which is crucially important for gauge theories and quantization of gravity [2][3][4][5][6][7][8][9][10]. See also [11][12][13] and references there.
Importance of the heat kernel approach was understood long ago both by mathematicians [14][15][16][17][18][19] and physicists [20][21][22]. But the efforts of mathematician were mainly focused on general estimations and theorems regarding the proper time expansion of the functional trace of the heat kernel with curvature invariant coefficients [17][18][19][23][24][25][26][27], whereas the physicists would consider the two-point heat kernel itself with the separate point arguments [6,21,22]. This would give essentially more flexibility and efficiency in obtaining these coefficientswith the ultimate goal of physical applications in UV renormalization and gradient expansion. This is where the difference between the expression for this kernel for second order and higher order operators starts explicitly showing up. Gilkey-Seeley approach, which is based on functorial methods [28][29][30][31], does not feel this difference, while the Schwinger-DeWitt technique, which explicitly generates recurrent equations for the two-point HaMiDeW (Hadamard-Minakshisundaram-DeWitt) coefficients and their coincidence limits, is very vulnerable to the choice of the leading order heat kernel ansatz and breaks down when it is inappropriately chosen.
Consider a generic minimal second-order operator F (∇) = −∆ + . . . whose covariant derivatives form a Laplacian ∆ = g ab ∇ a ∇ b acting in a curved d-dimensional space-time with the coordinates x = x a and the Riemannian metric g ab . Then the ansatz for its heat kernel K F (τ |x, y) = e −τ F (∇) δ(x, y) has the form τ n a n [F |x, y], (1.1) where σ(x, y) = l 2 (x, y)/2 is the Synge world function and l(x, y) is the geodetic distance between the points x and y. In fact this ansatz has a semiclassical nature. This is because its exponential coincides with the principal Hamilton function S(τ |x, y) = σ(x, y)/2τ of the particle moving in the x-space and fictitious imaginary time τ with the Hamiltonian F (∇), and the preexponential factor is just the square root of the Van Vleck-Morette determinant D(x, y)/(2τ ) d = det ∂ 2 S(τ |x, y)/∂x a ∂y b . So the expansion in powers of the proper time corresponds to the conventional semiclassical expansion in . Two-point coefficients of this expansion then satisfy simple recurrent equations which can be systematically solved for their coincidence limits a n [F |x, x]-local invariants of the space-time curvature and the coefficients of the operator F (∇). Note the property of the expansion (1.1) that it isolates essentially singular part of the heat kernel in the exponential, which vanishes in the coincidence limit y = x, whereas important physical information is contained in the HaMiDeW-coefficients of the regular expansion a n [F |x, x].
It is straightforward to formally extend this semiclassical ansatz to higher-derivative or pseudodifferential operators of the minimal form F (∇) = (−∆) ν + . . . with some integer or noninteger ν, but this extension fails to generate solvable recurrent relations for the generalized Schwinger-DeWitt coefficients. The origin of this difficulty is that this semiclassical approach fails to perform the resummation of all negative powers of τ in the exponential, and the infinite power series in τ turns out to include infinitely many of its negative powers. Essentially singular part of the τ -expansion at τ → 0 does not get isolated in the exponential and does not seem to vanish in the coincidence limit y = x as it happens in (1.1) for second-order operators.
Apparently due to this difficulty the heat kernel method was only indirectly used in physical applications with higher order operators. Numerous problems like regularization by higher order covariant derivatives or higher derivative theories, namely, R 2 -gravity [32,33], nonlocal and superrenormalizable models [34,35] and Hořava-Lifshitz models [36,37] were treated by means of the reduction to minimal second-order operators which allow one to use the expansion (1.1). Such a reduction technique for a wide class of theories with the generalized causality condition was suggested in [6] and actually allowed to circumvent the use of the proper shorttime expansion of the form (1.1). Discussion of the heat kernel method for higher-order operators within similar reduction, functorial or other methods can be found in [29,[38][39][40][41][42], see also a series of papers by Gusynin et al. [43][44][45][46][47][48].
Nevertheless, the heat kernels of higher-order differential operators are themselves important as explicit objects, because these kernels represent the building block of Green's functions of these operators, which are needed not only in the UV limit of their coinciding arguments. Moreover, a consistent version of the expansion (1.1) for higher-derivative operators and nonlocal operators of pseudodifferential type should be a source of recurrent relations for the generalized HaMiDeW-coefficients, and this generalization is expected to be a much more powerful tool than the reduction technique mentioned above. So possible applications of the standard method to higher-order operators does not make it less interesting to study their heat kernels directly. This is the goal of the present paper.
To understand the nature of the generalization of (1.1) for minimal higher-derivative operators it is enough to consider the case of a flat space-time of the Euclidean signature with the world function σ(x, y) = (x − y) 2 /2 and the operator F (∇) = (−∆) ν -the νth power of the Laplacian ∆ = δ ab ∂ a ∂ b , so that the standard heat kernel takes the translationally invariant form Then, the basic fact for a generic and not necessarily integer ν can be formulated as where E ν,d/2 (z) is what we call generalized exponential function (GEF). It is defined as a two-parameter family of functions represented by the Taylor series This function obviously reduces to the usual exponential function at ν = 1 and recovers the standard Gaussian behavior of the heat kernel. On the contrary, for other values of ν it performs resummation of negative powers of the proper time, which is impossible with the usual semiclassical ansatz. It should be emphasized that the functions (1.3) were originally utilized in higher derivative models in [38] and [45], the latter paper also including their series expansion (1.4). Later they were used in the context of anomalous diffusion theory [49] and in application to Hořava-Lifshitz models [50]. However, thus far no systematic studies of these functions were undertaken and their potential role for the extension of the HaMiDeW technique to modified field and gravity models was underestimated. The goal of this work is to fill up this omission.
The paper is organized as follows. In Sec. 2 we derive the heat kernel of the operator (−∆) ν and its associated GEF in the form of the Taylor expansion and present its integral representation in terms of Bessel functions. In Sec. 3 we discuss the properties of the generalized exponential functions E ν,α (z) and consider their Mellin-Barnes integral representation generating their asymptotic behavior at z → ∞, which is responsible for the short time, τ → 0, or large |x| → ∞ limit of the heat kernel (1.3). Interestingly, this limiting behavior turns out to be different for fractional and integer powers ν. Contrary to the second-order case this asymptotics is power-law for fractional ν and exponential for integer ν, and moreover features oscillations for growing |x|. For integer powers ν this property is demonstrated in Sec. 4 where the asymptotic behavior of GEF is compared with the semiclassical heat kernel ansatz and the saddle-point approximation for the momentum space integral representation. In concluding section we briefly discuss further application of GEF to generic minimal and nonminimal higher-derivative operators in curved space-time, which will allow us to build a solvable recurrent equations for the full set of generalized HaMiDeW-coefficients.

THE HEAT KERNEL OF THE POWER OF LAPLACIAN
For the operator F = (−∆) ν its heat kernel has an obvious momentum space representation where k = | k | = √ k 2 and kx = k a x a . For ν = 1 this integral defines the well-known fundamental solution (1.2).
Note that the heat kernel (2.2) is invariant with respect to O(d)-rotations and homogeneous where c is an arbitrary constant. Therefore, it should have the form (1.3), where E ν,d/2 (z) is some unknown function of the ratio −x 2 /4τ 1/ν . Since it stands in place of the exponential function in the usual expression for the heat kernel (1.2), we will call it the generalized exponential function (GEF).
Let us find the expansion of the generalized exponential function E ν,d/2 (z) in its Taylor series. Using the relations where σ = x 2 /2, it is easy to verify by induction that for an arbitrary function f (cσ), where c is a constant, the following differentiation rule holds ! are the binomial coefficients and Γ(s) is the Euler gamma function. Then for f (z) = E ν,d/2 (z) and c = −1/2τ 1/ν we obtain at x = 0 On the other hand, these quantities can easily be calculated directly Comparing the expressions (2.7) and (2.8), we find all the derivatives of E ν,α (z) at z = 0 and thus get its Taylor expansion (1.4).

Bessel functions representation
Another integral representation expresses the generalized exponential function E ν,α (z) in terms of the Bessel function J α (z) or the Bessel-Clifford function C α (z). The latter is determined by the series and related to the Bessel function J α (z) by the equation which removes its branching point at z = 0 Therefore, the Bessel-Clifford functions have no singularities and represent single-valued entire functions on the whole complex z plane. Bessel function representation of K ν,d (τ, x) follows from integration over angles in the momentum space integral (2.2), which reads as is the volume of (d − 2)dimensional unit sphere, θ is the angle between the vectors x and k in (2.2) and x = | x |. Expanding exp(ikx cos θ) and integrating over θ on account of As a result the heat kernel has the following integral representation with the Bessel-Clifford function 15) Under the change of integration variable µ = k 2ν τ comparison of this expression with (1.3) gives the relevant integral representations for the GEF E ν,α (z) Note that substitution of the expansion (2.9) for the Bessel-Clifford function into (2.16) directly leads to the expansion (1.4), which confirms this representation. Another remark is that the functions E ν,α (z) and K ν,d (τ ; x) are directly related to Bessel-Clifford (2.9) and Bessel functions (2.10) in the limit ν → ∞. Indeed, replacing in this limit Γ (α + m)/ν by ν/(α + m) in the expansion (1.4), one has Interestingly, the heat kernel of (−∆) ∞ becomes independent of the proper time parameter τ because of the obvious limit τ 1/ν − −−− → ν→∞ 1.

GENERALIZED EXPONENTIAL FUNCTIONS AND THEIR PROPERTIES
Various properties of GEF follow from the Mellin-Barnes integral representation of this function. This representation can be obtained by converting the series (1.4) into the contour integral in the complex plane of an auxiliary parameter s, such that the residues at simple poles of the integrand generate various terms of this series. It is easy to guess that this integral reads Then the inverse Mellin transform obviously gives The location of the poles of ε ν,α (s) and the contour C on the complex s plane is schematically shown in Fig. 1. The function Γ(s) has a sequence of poles at l k = −k, k = 0, 1, 2, . . . running to the left with residues (−1) k /k!. And the function Γ((α − s)/ν) has poles at r k = α + kν, running to the right. The contour C begins at −∞ − i , runs under the real axis, bends around 0 and returns to −∞+i . The integral (3.1) is equal to the sum of residues at the poles l k , which exactly gives the series (1.4).
It is possible, however, that not all the points r k = α + kν are poles of the function ε ν,α (s), since they can be canceled by the poles of the function Γ(α − s) in the denominator. In particular, the point r 0 = α is never a pole of ε ν,α (s). Three cases are possible: if ν is irrational, then all other points r k are poles; if ν = p/q is an irreducible fraction, then the poles r qk are also canceled; and, finally, if ν is integer, then all the poles r k are canceled.
To determine the asymptotic behavior of the function E ν,α (−z), we deform the integration contour C into the contour C w , coming from w − i∞ vertically to w + i∞, where w = Re r k , w > 0 (see Fig. 1). Then we have For z → ∞ the residues decrease as powers z −r k , whereas the integral term is obviously O(z −w ) because of the constant factor z −w in the integrand. Then for noninteger ν by pushing w → +∞ we obtain the sum of all residues as the power-like asymptotic expansion of GEF (3.4) Note that the cancellation of the pole in Mellin image (3.2) results in the vanishing of the corresponding term in this expansion (3.4). In the case of integer ν all residue terms vanish. This means that in this case the function E ν,α (−z) decreases faster than any power of z, i.e. in an exponential manner. Then the function E ν,α (−z) is exactly equal to the integral over the contour C w for any positive w, and its exponential asymptotic behavior is entirely determined by this integral. This case is very special, and we will consider it in detail later in Sec. 4.
The properties of GEF of the above type are strongly associated with the properties of ratios of gamma functions and their products, which are in fact a part of the theory of the Fox-Wright Ψ -functions [51,52] and more general Fox H-functions [53][54][55][56][57]. Namely, our generalized exponential function, defined by the series (1.4), is a special case of the Fox-Wright function p Ψ q [ (a, A); (b, B); z ] labeled by numerous parameters For completeness we briefly give the definition of these functions in Appendix A, in particular because their theory gives subtle details of the domain of parameters ν and α, where the function E ν,α (z) is well defined. For the generalized exponential function E ν,α (z) we have the Mellin-Barnes integral (3.1) and two power series: the first series (1.4) in z near 0 and the second series (3.4) in 1/z near ∞. As it follows from the general theory of H-functions (see Eqs. (A.5), (A.6) and the statement below them), there are three possible cases: for ν > 1/2 the first series converges absolutely on the whole complex plane z, and E ν,α (z) is thus an entire function and has the essential singularity at z = ∞, while the second series diverges and is asymptotic for z → ∞. When ν < 1/2 the situation is the opposite: the series (3.4) absolutely converges on the whole complex plane z, except for the essential singularity at z = 0, and the series (1.4) diverges and is asymptotic for z → 0. Finally, for the critical value ν = 1/2 the series (1.4) converges inside the circle |z| < 1/4 and diverges outside it, while conversely the series (3.4) converges outside this circle and diverges inside it.
It turns out that for the critical value ν = 1/2 the series (1.4) can be summed up analytically 1 . Indeed, using the Legendre duplication formula (2.13) and the well-known expansion we can easily find that Thus, for ν = 1/2 GEF not only have power-law asymptotic behavior, but they really are power functions. The series (1.4) diverges for |z| > 1/4 due to the existence of a pole at the point z = 1/4. It is not difficult to verify that in this case even terms of the series (3.4) vanish and odd terms converge to the function (3.7) in the circle |z| > 1/4. As a result, the operator √ −∆ in a flat d-dimensional space has the following heat kernel This expression is related to holographic [58,59] and brane world [60][61][62] applications of effective action because it represents the massless limit of the simplest brane-to-bulk propagator e −τ √ M 2 − δ(x), M → 0, [61,62] and may be interesting in the context of the discussion of fractional powers of generalized Laplacians in [25].
The terms of the series (1.4) for E ν,α (z) are well defined for complex We note, however, that the singularities at the points α = −n (and hence at all the points for integers ν) are removable, since the poles of the gamma functions in the numerator and denominator cancel each other. Expanding Γ(−n+ ) ∼ (−1) n /n! we can redefine the coefficients in (1.4) otherwise.
(3.10) Thus, for integer ν the function E ν,α (z) is an entire function of z for any values of α, and for noninteger ν > 1/2 it is an entire function for all values of α except α = −n − kν = −m with positive integer k, n and m. Consequently, the functions K ν,d (τ, x) are well defined not only for all integer ν and d, but also for fractional ν and d satisfying these conditions.
The graphs of the functions E ν,α (z) and K ν,d (τ, x) for various values of the parameters, obtained by numerical summation of the series (1.4) in MATLAB, are shown in Figs. 2-4. Important distinction from the case of a monotonic exponential falloff for ν = 1 is that the heat kernel for ν = 1 oscillates as a function of Other interesting properties of E ν,α (z) include the following simple differentiation rule For integer β, it can be verified directly by differentiating the definition (3.5). However, this relation makes sense also for all β such that E ν,α+β (z) is well defined for noninteger ν (and for any complex β if ν is integer). For negative integer β it will give the principal primitive of the function E ν,α (z). For noninteger β, the symbol d β /dz β should be understood as a certain operator of fractional integrodifferentiation. Thus for each ν in the range 1/2 < ν ≤ ∞ the family of functions E ν,α (z) is closed under the operation of integrodifferentiation.
For noninteger ν the expression (1.3) is the solution of the heat equation in which the nonlocal operator (−∆) ν should be understood as a pseudodifferential operator defined by the Fourier transform [63]. The corresponding equations are called fractional diffusion equations and have been widely discussed in [49,[64][65][66][67]. However, in these papers fractional equations are usually considered in (1 + 1)-dimensional (τ, x)-space , i.e. the case of d = 1 in our notations.

INTEGER POWER OF LAPLACIAN AND SEMICLASSICAL EXPANSION
As we see, the asymptotic behavior of GEF E ν,α (−z) at z → ∞ is critically different for noninteger and integer values of ν. It is power-law for noninteger ν corresponding to the nonlocal operator (−∆) ν and quasiexponential O(z −∞ ) for integer ν corresponding to local differential operators of order 2ν. But for the heat kernel (1.3) this limit is associated with the semiclassical limit x 2 /4τ 1/ν → ∞, the proper time τ → 0 playing the role of . On the other hand, semiclassical approximation for the solution of the heat equation (or the Schrödinger equation in the imaginary time) begins with the Pauli-Van Vleck (or WKB) ansatz [68] det − 1 2π where S(τ |x, y) is the principal Hamilton-Jacobi function, i.e. the solution of the Hamilton-Jacobi equation The operator F (∇) with the space-time gradients replaced by the canonical momenta p. We work in Euclidean time, which explains the absence of imaginary i-factor. The derivation for a generic Hamiltonian not necessarily polynomial in derivatives (momenta) can be found in Appendix C of [69].
For the power of Laplacian this Hamiltonian equals F (p) = (−p 2 ) ν , and the solution of this equation readily expresses as This leads to the absolute value of the preexponential factor Regarding its phase factor, which should be determined here by the correct choice of the branch for the fractional powers, the Pauli-Van Vleck algorithm does not give this information in contrast to the trivial case of ν = 1. Neither does it prescribe a definite linear combination of these branches in the heat kernel asymptotics. Below all this will be attained by two different methods -application of the general technique of Fox H-functions and by the steepest descent approximation. Note now that in view of the discussion in the previous section this quasiexponential behavior is completely impossible for noninteger ν with a power-law behavior, so that standard semiclassical expansion seems to break down for nonlocal operators (−∆) ν . Therefore in what follows we consider only the case of positive integer powers. To underline this, we will further denote the power of the Laplacian by N , when it is integer, and by ν, when it can be either integer or noninteger. The goal of this section is to find the heat kernel asymptotic expansion for this higher-derivative case by an alternative method which allows to get correct complex branches and to compare them with the semiclassical result of the above type.
6) which for µ > 0 and growing K > 0 will generate decreasing terms of the exponential expansion at z → ∞. However, what is integrated in (3.1) is not just a gamma function of this type, but rather a nontrivial ratio of those. This ratio ε N,α (s), which is given by Eq. (3.2), can nevertheless be converted to the series of gamma function terms of the above type, so that the s-integration can be successfully done. For this purpose we use, first of all, the Euler reflection formula to provide positive coefficients of the integration parameter s in the arguments of all gamma functions (just like in (4.6)), This makes the sequence of gamma function poles running to the left of the complex plane of s. For integer N the ratio of sines reduces to the sum of complex exponential functions After the substitution of (4.8) into the Mellin transform (3.1) this leads to the sum with the phase factors both in the coefficients of this sum and in the arguments of the new functionsẼ ν,α (−e iωj z) which can be called the generalized exponential functions of the second kind. They read Eq. (A.7) generates the asymptotic expansioñ with the coefficients E m independent of the integration variable s. These coefficients start with E 0 = 1 and they are systematically calculable by the procedure of Appendix B. The essence of this expansion is that for large s it runs over ever decreasing terms, each term being smaller than the preceding one in view of the obvious relation Γ(µs−a−m−1) = Γ(µs−a−m)/(µs−a−m−1). The Mellin transform (4.13) of (4.16) on account of (4.6) yields the following asymptotic expansion for GEF of the second kind In contrast to the GEF of the first kind E ν,α (z), which has no singularities, the functionẼ ν,α (z) is singular at z = 0. However, it always has a simple exponential (i.e., not a power-law) asymptotic behavior (4.17) for z → ∞. Second, it is monotonic for −∞ < z < 0 without oscillations characterizing GEF of the first kind. Moreover, one can say that the source of these oscillations is in fact the set of sines in Eq. (4.8) and phase factors in Eq. (4.12).
Thus in view of the decomposition (4.12) the asymptotic expansion of the generalized exponential function finally reads as a sum of N series of terms where both the amplitudes and phases depend on the phases ϕ j of the complex factors (4.12) (4.19) Consequently, the expression in (1.3) with z = x 2 /4τ 1/N gives the heat kernel asymptotics for fixed τ and | x | → ∞ or for fixed | x | and τ → 0, (4.20) The phase factors e iϕj here coincide with the set of fractional powers (−1) (N −1)/(2N −1) in the coefficient of the Hamilton-Jacobi function (4.4). This confirms the semiclassical ansatz (4.2) along with establishing a concrete choice of the linear combination of its complex branches. The leading order term of this expansion consists of two complex conjugated branches corresponding to the maximal real part of the exponential with j = 0 and j = N − 1 and the phase factors exp(∓iϕ),

Steepest descent approximation
Alternatively this result can be reproduced by the steepest descent method which is the basis of the semiclassical approximation with the small parameter τ → 0. The change of variables p = τ 1/N k, y = τ (N −1)/N x, converts the momentum integral (2.2) to the form where S(p) = (p 2 ) N − ipy. Its short time τ → 0 asymptotics follows from the saddle points of the complex action S(p) at which this action is stationary, ∂S(p)/∂p a = 0. These 2N − 1 points read  27) where the set of phase factors is determined by the relations j = 0, 1, . . . , 2N − 2, and in fact extends the range (4.19) up to 2N − 1.
The most complicated part of the saddle point method is the proof of the existence of a correct steepest descent integration contour and the choice of relevant saddle points which should contribute to the asymptotics in question [70]. In a particular case of the integral (2.2) there is an obvious hint on their choice, confirmed by a rigorous analysis in [70], that the points with cos ϕ j < 0, j ≥ N , contribute the terms exponentially growing with | x |, which contradicts the known | x | −∞ falloff at | x | → ∞. Therefore only the remaining N points with j ≤ N − 1 can contribute to the asymptotics of the integral. Note that their respective phase factors coincide with those defined by the Eqs. (4.12) and (4.19), and thus correspond to the set of terms in the decomposition of the generalized exponential function into the sum of GEF of the second kind. Thus the steepest descent method leads to the same result as the Mellin transform within the formalism of Fox H-functions.

Nonuniformity of the semiclassical expansion
Semiclassical expansion (4.22) clearly shows a principal difference of higher-order operators from the secondorder case N = 1. For N > 1 the short time expansion (4.22), τ → 0, does not stand the limit | x | → 0 because it involves negative powers of | x |. In particular, it does not maintain the initial condition K ν,d (0, x) = δ(x). The exact heat kernel of course satisfies this condition, because for any smooth test function f (x) where we used the relation (3.2). Interestingly, in the recovery of this result for integer N each of the N terms of the decomposition (4.12) contributes one and the same 1/N th part of it, because their dependence on ω j drops out in the limit α → d/2, At the same time, if we try to reproduce the same result by using N branches (4.27) of the semiclassical expansion, corresponding to different terms of the above decomposition then the result will be critically different.
Each K (j) N,d (τ, x) is singular at x = 0, but this singularity is integrable. However, the result of this integration is different from (4.31) Even more striking discrepancy between the exact heat kernel and its asymptotics is that, while all the terms of the latter are singular in the limit x → 0, the GEF (1.4) and the exact K ν,d (τ, x) are both well defined in this limit K ν,d (τ, 0) = 1 The short time expansion of this coincidence limit which is a main goal of the Seeley-Gilkey technique [17][18][19]24] runs in powers of τ 1/ν , whereas the expansion (4.22) goes in the other fractional powers τ 1/ (2N −1) . The source of all these discrepancies 3 is the fact that, contrary to the second order Laplacian, the heat kernel asymptotic expansion is not uniform in x → 0. While in this limit the expansion (1.1) for a minimal second order operator F (∇) = −∆ + . . . is just a regularization of the delta function, for N > 1 the expansion (4.22) fails for x → 0 and does not have a chance to reproduce correct initial conditions with a pointlike support at x = 0. Obviously, there is no such a discrepancy in the case of N = 1 with a single j = 0 branch of the heat kernel expansion, so that the coincidence limit y = x can be directly taken in the asymptotic expansion (1.1). 3 There is additional controversy with the result (4.33)-while all K (j) N,d (τ, x) with j = 0 and j = N − 1 are exponentially subdominant and should be discarded according to asymptotic expansion theory, all their integrals (4.33) are of the same order of magnitude.

CONCLUSIONS
Thus we obtained the expression (1.3) for the heat kernel K ν,d (τ, x) of the operator (−∆) ν in the d-dimensional flat space, which is a direct generalization of the wellknown heat kernel (1.2) to local higher derivative and nonlocal (pseudodifferential) operators. This generalization is represented in terms of the newly introduced two-parameter family of generalized exponential functions E ν,α (z), α = d/2, z = −x 2 /τ 1/ν , defined by the Taylor series (1.4) and related to Fox-Wright Ψ -and Fox H-functions. We studied various properties of these functions and their integral representations. They include, in particular, the Mellin-Barnes representation which allows one to find, by the technique of Fox H-functions, their asymptotic expansion in the limit of z → ∞ corresponding to two equivalent asymptotics of the heat kernel as τ → 0 or as |x| → ∞.
This expansion turns out to be critically different for integer and noninteger values of ν. In contrast to the exponential behavior, anticipated on the ground of semiclassical considerations with τ → 0 playing the role of , for noninteger ν, that is for nonlocal operators (−∆) ν , this is a power-law falloff. For integer ν > 1 this asymptotic expansion matches with the exponential Pauli-Van Vleck ansatz or steepest descent approximation, but it essentially differs from the pure Laplacian case of ν = 1. In particular, this asymptotic expansion is not uniform for all values of | x | → 0 and does not stand a singular coincidence limit x = 0 which is, on the other hand, easily accessible directly from the short time expansion for ν = 1. In addition, the heat kernel for higher-derivative and nonlocal operators with ν = 1 features oscillatory behavior in x-space contrary to the monotonic exponential falloff for the pure Laplacian.
The coincidence limit of the heat kernel K F (τ |x, x) and its functional trace are the objects of major interest in the Schwinger-DeWitt technique for quantum effective action. For generic minimal second order operators with xdependent coefficients it is based on the expansion (1.1) for K F (τ |x, y) with split arguments (generically lacking translation invariance). This coincidence limit is also the subject of mathematical Seeley-Gilkey treatise by various functorial methods not directly appealing to offdiagonal elements of K F (τ |x, y). The absence of a nonsingular coincidence limit in the heat kernel asymptotic expansion for ν = 1 considered above seems to invalidate any attempt to use it for some generalization of the Schwinger-DeWitt method. So these expansions are not very physically interesting in field theory applications. However, a systematic way of the short time expansion of the heat kernel, which underlies UV properties of field models, is provided by the recurrent equations for the coefficients of the expansion (1.1), these equations heav-ily relying on the off-diagonal K F (τ |x, y). The Seeley-Gilkey functorial methods are not so universal and powerful enough to yield everything what physicists need in quantum gravity and other applications. For example, Hořava-Lifshitz gravity [36,37] are encumbered with the necessity of working with higher order and nonminimal operators whose leading symbol does not reduce to the power of Laplacian and, therefore, go outside of the scope of functorial methods. Derivation of recurrent equations for the two-point coefficients of (1.1) generalized to such operators then becomes indispensable. The generalized exponential functions introduced above provide fundamental building blocks of such recurrent equations, and the lack of uniformity of their asymptotics does not make them less efficient. Note that the characteristic feature of the Schwinger-DeWitt expansion (1.1) is a single overall exponential factor absorbing all essentially singular dependence on τ → 0. The attempt to directly generalize this expansion to ν = 1 with a single semiclassical exponential factor fails because it generates infinitely many negative powers of τ .
On the contrary, resummation of these singular terms can be performed with the aid of the generalized exponential functions, but in contrast to (1.1) these functions will not form a single overall factor, but rather comprise the series of terms with different α-parameters. As we are going to show in the coming paper [71], for a minimal differential operator F of an (integer) order 2N in a curved Riemannian space [(x − y) 2 /2 → σ(x, y)] the following generalization of the expansion (1.1) holds The generalized HaMiDeW-coefficients a j [F |x, y] satisfy the manageable chain of recurrent equations, which can be solved for the coincidence limit x = y. Note also that the coincidence limit of (5.2) is well defined, even though the asymptotics of the underlying E N,α (−z) are not uniform for z → 0. In fact, these asymptotics are not needed for this limit. Since E N,α (0) = Γ(α/N )/N Γ(α) we have the following expansion for the coincidence limit So GEF should be treated as entire building blocks of the formalism, the operations with them being based on their simple differentiation rule (3.11) and the value at z = 0. In our next paper [71] we will consider various properties of the generalized HaMiDeW-coefficients in (5.2). In particular, we will prove the generalized "functorial property" for an arbitrary power λ of a differential operator F , (previously known only in the coincidence limit x = y [28]), and also easily reproduce and extend the results of [28] for higher order operators. These coefficients and the computational methods based on them promise to be very efficient and are likely to simplify the calculation of beta functions for theories with higher derivatives and Hořava-Lifshitz type models [36,37]. All this makes the use of GEF and the associated heat kernel coefficients very prospective. The above formalism seems equally applicable to the case of generic noninteger ν. This case is, in particular, important in superrenormalizable quantum gravity models [34,35,72,73], within analytical regularization of Feynman graphs [74] or for the calculation of UV counterterms in Hořava-Lifshitz gravity models. For example, in (3+1)-dimensional Hořava gravity cubic in spatial curvature counterterms follow from the heat kernel of the operator which is a square root of the sixth order nonminimal differential operator [37]. However, as we saw above there is a number of new features (like the power-law heat kernel asymptotics confronting their exponential analogue for integer N or the mismatch with the semiclassical expansion) which might backfire under indiscriminate extension of this method. Here we only briefly comment on possible modifications due to these subtleties.
One modification follows from the recovery of the heat kernel diagonal elements by inverse Mellin transform from the operator zeta function F −s δ(x, y) | x=y used in [25]. For the operators of the form F = H ν , where H is a Laplace type (minimal second order) operator, and noninteger ν this method leads to additional terms [25] (5.6) While the coefficients A j (x) are in one-to-one correspondence with the HaMiDeW-coefficients a j [H|x, x] of the operator H and are local quantities, the coefficients B k (x) are determined through the values of the zeta function at certain values of s. Rather than being expressed in terms of a j [F |x, x], they turn out to be nonlocal and irrelevant to UV renormalization because they do not contribute to UV divergences in view of analytic expansion in τ starting with a linear term. However, according to our method, just in the case of operators of the form F = H ν , such additional terms do not arise.
Another modification known from mathematical studies [23,25] is the origin of logarithmic terms in the proper time expansion of the heat kernel. For special values of noninteger ν leading to −α = n + kν = m [with positive integer k, n and m-see the discussion after (3.9)] the functions E ν,α (z) in the expansion (1.2) are not defined because of gamma function singularities. This exceptional case can occur, in particular, for even order roots of the Laplace type operator in odd spacetime dimensions. In this case the logarithmic terms appear [23,25] (5.7) which are again unrelated to renormalization of UV divergences.
All these modifications can apparently be attributed to the fact that zeta-function approach of [23,25] actually represents a regularization which for nonlocal theories (corresponding to noninteger values of ν) leads to different results 4 . Absence of uniformity of the asymptotic small time expansion in the vicinity of the heat kernel diagonal discussed above shows up for noninteger values of ν. The search for an asymptotic expansion of GEF and heat kernel that would be uniform for all x and y (analogous, for example, to the uniform WKB asymptotic expansion of Legendre functions [75]) apparently could have resolved the problem of these discrepancies. This however goes beyond the scope of this paper, in particular, because nonuniformity of the asymptotic expansion of GEF is harmless in the calculation of the functional trace (5.1) if one uses generalized exponential functions as building blocks of the expansion and takes their exact values at zero argument.
To summarize, the generalized exponential functions can serve as a very efficient tool in quantum field theory and quantum gravity. Moreover, their connection with fractional calculus opens the prospect of applying the obtained heat kernels far beyond the area of QFT. This includes the theory of fractional differential equations which can be effectively used to construct phenomenological models of fractal media, systems with memory and nonlocal interaction. Numerous applications of fractional calculus to physical problems are discussed, for example, in [76] and references therein. The Fox-Wright Ψ -functions p Ψ q [(a, A); (b, B); z] are labeled by two sets of parameters A k , a k , r = 1, . . . , p, and B j , b j , j = 1, . . . , q, among which A k and B j are real and positive. These functions are defined by their Taylor series also with real and positive A i and B j . The poles l i,k of the gamma functions Γ(b i + B i s), i = 1, . . . , m, enumerating the poles index k being integer, run to the left of the complex plane of s, whereas the poles r j,k of Γ(1 − a j − A j s), j = 1, . . . , n, run to the right. It is assumed that the parameters A j , a j , B i and b i are such that these poles do not match, l i,k = r j,l . Then the contour of integration C is chosen to pass from −i∞ to i∞ and to separate the poles l i,k and r j,k . The Fox H-functions are related to Fox-Wright Ψfunctions in exactly the same way as the well-known Meyer G-functions to generalized hypergeometric functions. Obviously The general theory of H-functions and H-transforms can be found in [53][54][55][56][57]. Here we briefly sketch their main properties and the way of handling their asymptotic behavior. This behavior is characterized by the following three combinations of their parameters Note that the structure of the expression (A.3) allows one to relocate gamma functions between the numerator and the denominator using the Euler reflection formula (4.7). Under this operation only the parameters m and n change, while the parameters p, q, µ, β and a, as it is easy to see, remain intact.
The main result, based on the use of the Stirling formula, is as follows: for µ > 0 the contour C in (A.2) can be closed on the left of the complex plane, then m series obtained by summing the residues at the poles l i,k will converge absolutely on the whole complex z plane, defining, generally speaking, a multivalued function with an essentially singular point at z = ∞. If in this case we formally close the contour C on the right, then the sum of the residues at the poles r j,k will determine the asymptotic (divergent) power series as z → ∞. For µ < 0 the situation will be exactly the opposite: the sum of residues at the poles r j,k will absolutely converge at z = 0, and the divergent series of residues at the poles l i,k will determine the asymptotic behavior of the function at z → 0. Finally, in the case of the critical value µ = 0 the series obtained by closing the contour C on the left will converge inside the circle |z| < β −1 , and the series obtained by closing the contour C on the right will converge outside of it.
Exponential asymptotic behavior for z → ∞ appears when h m,n p,q [s] does not have any rightgoing poles r j,k , i.e., when in the expression (A.3) either the functions Γ(1 − a j − A j s) in the numerator are completely absent (n = 0), or all their poles are canceled with the poles of the functions Γ(1−b i −B i s) in the denominator (as it happens in the case of the functions E ν,α (z) for an integer ν). The general recipe for finding exponential asymptotics, which is explained in detail in [53], is the following: one first needs to use the Euler reflection formula to relocate gamma functions so that they have only leftgoing poles, i.e. to convert the expression (A.3) to the form, when all gamma functions with the coefficients A j are in the denominator, and all gamma functions with the coefficients B i are in the numerator. After that, one needs to use the asymptotic expansion for the ratio of gamma function products,