Non-perturbative corrections to the one-loop free energy induced by a massive scalar field on a stationary slowly varying in space gravitational background

The explicit expressions for the one-loop non-perturbative corrections to the gravitational effective action induced by a scalar field on a stationary gravitational background are obtained both at zero and finite temperatures. The perturbative and non-perturbative contributions to the one-loop effective action are explicitly separated. It is proved that, after a suitable renormalization, the perturbative part of the effective action at zero temperature can be expressed in a covariant form solely in terms of the metric and its derivatives. This part coincides with the known large mass expansion of the one-loop effective action. The non-perturbative part of the renormalized one-loop effective action at zero temperature is proved to depend explicitly on the Killing vector defining the vacuum state of quantum fields. This part cannot be expressed in a covariant way through the metric and its derivatives alone. The implications of this result for the structure and symmetries of the effective action for gravity are discussed.


Introduction
A construction of self-consistent quantum theory of gravity remains elusive due to several conceptual and technical problems. The main technical problem is, of course, a high non-linearity of classical theory of gravity (general relativity) and, as a result, nonrenormalizability of its quantum counterpart. The major conceptual issue is related to the problem of time in quantum gravity (see for a review [1] and also [2][3][4][5]) or, put another way, it regards the problem of definition of a natural vacuum state and a representation of the algebra of observables. On the other hand, there is a widespread belief resting on perturbative calculations on a flat background that the second problem does not actually exist. According to this point of view, quantum gravity is a mere another one effective quantum field theory similar to the Fermi theory of weak interactions. The aim of the present paper is to show that such a viewpoint is somewhat naive as long as it does not take into account non-perturbative corrections to the effective action.
One of the most powerful methods to find the non-perturbative corrections to the effective action (the generating functional of the one-particle irreducible Green functions) is the celebrated background field method (see, e.g., [6,7]) with its renormalization group improvements (see, e.g., [7]). It is, perhaps, the only one for many particle quantum systems in dimension D ≥ 3 without extra symmetries. We shall use this method to obtain the one-loop non-perturbative (in the gravitational constant) contributions to the effective action of gravity from a massive scalar field at zero and finite temperatures (for the renormalization group improvement of the perturbative contributions see, e.g., [8]). We shall derive the explicit expressions for these corrections in the case of a stationary slowly varying in space gravitational background (the stationary infrared limit). As far as we know, this problem has not been solved yet for such a general formulation.
The fact that the effective action has to possess the non-perturbative terms of the form that we shall derive in this paper was repeatedly noticed in the literature (see, e.g., [3,9,10]). However, neither the explicit form of these corrections for a sufficiently wide class of background metrics nor even their expression terms of the metric and coincides with the standard expression for the conformal anomaly [16,20,21,[24][25][26][27].
Section 4 is the heart of the paper, where all the general results of the preceding sections are collected together in order to obtain the high-temperature expansion of the free energy of a scalar field and describe its properties. Also, in this section, we heavily rely on the results of the papers [16] and [28]. In fact, it is the combination of the results of these papers that made it possible to derive a complete high-temperature expansion with the non-perturbative contributions. In Sec. 4.1, we shall develop a perturbative procedure for the heat kernel that allows us to deduce systematically non-perturbative corrections to the effective action. Using the procedure elaborated, we shall find the first correction to the leading (Gaussian) contribution to the heat kernel obtained in [28] (for the Euclidean case see [29,30]). The counting scheme, which sorts an infinite set of Feynman diagrams for the heat kernel and distinguishes the most relevant ones, will be also presented there. Section 4.2 is devoted to evaluation of the explicit expressions for non-perturbative contributions to the high-temperature expansion. From the technical point of view, this problem is reduced to a calculation of certain double integrals in the complex planes: one integral is over the proper-time and another one is over the energy of a mode. This issue is solved for both massive and massless cases in the weak field limit.
In conclusion we shall outline some implications of the results obtained and the possible further directions of research. In Appendix A, we prove a certain property for a variant of the zeta function of the Laplace type operator coming from the hyperbolic type operator Fourier transformed over the time variable. This property is used in Sec. 2. In Appendix B, the general formulas of perturbation theory are gathered and the first correction to the leading contribution to the heat kernel is explicitly calculated.
Since, in the present text, we shall try to reduce the duplication of formulas from [16] and [28] to a minimum, the reader is strongly encouraged to have these papers close at hand. In the course of the discussion, several misprints made in [16] and [28] will be also corrected. Besides, in Sec. 2, we shall extensively employ the analytic regularization technique for singular integrals [31]. The acquaintance with Sec. 3 of [31] will be required. The knowledge of the notion of heat kernel and the heat kernel expansion technique [22] is also desired in Sections 2 and 3. In order to realize better the main idea of the procedure developed in Sec. 4 and the structure of density of states considered in Sections 2 and 4, it is recommended to know the results of [32,33] and especially [33]. As for the construction of quantum field theory (QFT) on a curved stationary background, we shall assume that the reader is familiar with the paper [18] and Sections 17, 18 of [6]. Of course, the other references we provide in the paper are also advised for reading.
We shall use the conventions adopted in [6] R for the curvatures and the other structures appearing in the heat kernel expansion. The square and round brackets at a pair of indices denote antisymmetrization and symmetrization without 1/2, respectively. The Greek indices are raised and lowered by the metric g µν which has the signature −2. Also we assume that the metric possesses the timelike Killing vector ξ µ : that allows us to make the decomposition ( [34], Sec. 84; [13,[35][36][37]) where g µ = ξ µ /ξ 2 is a one-form dual to the Killing vector (the Tolman temperature one-form). Notice that this decomposition is not the Arnowitt-Deser-Misner (ADM) one, but in some sense dual to it. The decomposition (3) is constructed by the use of the vector field, while the ADM one is associated with the system of hypersurfaces or, equivalently, with the integrable one-from. In case of a static spacetime, these decompositions coincide so long as the family of hypersurfaces is identified with the integral manifolds of the one-form g µ . In the system of coordinates, where ξ µ = (1, 0, 0, 0), we have the relations g ik g kj = −δ j i , ξ 2 det(−ḡ ij ) = g, g 00 +ḡ ij g 0i g 0j = (g 00 ) −1 , g i =ḡ ij g j0 , g 00 , −ḡ µν = g µν − ξ 2 g µ g ν = g 00 − (g 00 ) −1 g 0j g i0 g ij .
We have changed the sign of the metricḡ ij in comparison with [16]. The Latin indices corresponding to the space are raised and lowered by the positive-definite metricḡ ij . The curvatures associated with this metric will be distinguished by the overbars, e.g.,R. Note that we consider a general stationary spacetime, i.e., the Tolman temperature one-form is supposed to be non-integrable ( [38], App. C; [34], Sec. 88; [13,[35][36][37]), in general. The system of units is chosen such that c = = 1.
2 High-temperature expansion

General formulas
Let H(ω) be a Fourier transform of a kernel of the wave operator on a stationary background (see, for instance, (62)). Suppose the corresponding operator is of a Laplacian type, depends analytically on ω, has the spectrum bounded from above at fixed ω in the Hilbert space of square-integrable functions, and there is no accumulation points in the discrete spectrum. Consider the operator (cf. [39], Sec. 1.10) where the contour C runs along the imaginary axis from top to bottom and encircles the origin from the left. For the special case, ν = 0, H 0 + (ω) = θ(H(ω)) is the projector to the subspace of the total Hilbert space. This subspace is spanned on the eigenvectors of H(ω) corresponding to the positive eigenvalues.
If the system is placed in a sufficiently large "box" in space so that the spectrum of H(ω) is discrete, then H −ν + (ω) is trace-class for any ν ∈ C. Consequently, there exists the entire function of ν, ζ + (ν, ω) = Tr H −ν + (ω).
The investigation of passing to an infinitely large "box" can be found, for example, in [40]. If the operator H(ω) also possesses a continuous spectrum, then ζ + (ν, ω) = k θ(ε k (ω))ε −ν k (ω) + εc(ω) 0 dεn(ε, ω)ε −ν , Re ν < 1, where ε k (ω) are the eigenvalues of H(ω) and n(ε, ω)dε is the number of states of the continuous spectrum in the interval [ε, ε + dε] (usually, it is proportional to the volume of a system). The quantity ε c (ω) specifies the boundary of the continuous spectrum. The generalization of formula (8) and the following ones to the case of several zones is quite obvious. As an example, we present the density of states for a free scalar particle in a (d + 1) dimensional Minkowski space where D := d + 1. The quasiclassical approximation for the Laplacian type operators we study gives for ζ + (ν, ω) at large ω (see [16,[19][20][21][24][25][26][27] and below), In fact, to derive this asymptotics, one needs to substitute (9) to (8) and evaluate the integral. For those H(ω), which depend nonquadratically on ω, for example, for the wave operator in a dispersive media, it is also reasonable to expect the asymptotic behavior (10) so long as a media becomes transparent in the ultraviolet regime. Consider, in general, that ε ′ c (ω) > 0 for ω > 0 and there exists ω c > 0 such that ε c (ω c ) = 0. Also suppose that 1. n(ε, ω) is a smooth function of ω and ε for ω > 0, ε ∈ (0, ε c (ω));
The first condition is rather technical and may be weaken. It makes it possible not to take care of the existence of derivatives. The second condition says that the integral over ε in (8) converges on the upper integration limit. The last condition is necessary for the Gelfand-Shilov analytical regularization (in its standard form) [31] of the integral over ε in the neighborhood of ε = 0. If these conditions are met, the function (8) can be analytically continued to the region Re ν ≥ 1. So, the integral in (8) is understood as analytically regularized [31] in the case when zero belongs to the continuous spectrum of H(ω). As follows from the general procedure [31], the function (8) has simple poles at ν ∈ N under the above restrictions on the density of states n(ε, ω).
Using the function ζ + (ν, ω), it is easy to obtain the expression for the one-loop correction to the Ωpotential. Assuming the condition of the vacuum stability is fulfilled for the eigenvalues ε k (ω) of H(ω) at the points, where ε k (ω) = 0, (see, e.g., [41]) we have (see for details, e.g., [16]) The contributions from both particles and antiparticles are taken into account in this expression. If the vacuum is stable (11) then ζ + (ν, 0) = 0, i.e., H(0) has no positive eigenvalues. In Appendix A, we shall prove this statement deforming the wave operator of free fields, which possesses this property, into the wave operator with interaction H(ω). The proof given in Appendix A requires, additionally, the condition of a "smooth deformability" of the eigenvalues of a family of operators under consideration. In what follows, we put ζ + (ν, 0) = 0.
Integrating by parts in (12) and taking into account that ζ + (ν, 0) = 0, we get It is convenient to introduce the function On substituting (8) into (14), the integral I ν (µ) falls into two pieces corresponding to two summands in (8). Suppose µ does not lie on the real positive semiaxis ω, where ε k (ω) = 0 for some k. Then, taking into account (11), we have for the first contribution where ω k (ε) is the inverse function to ε k (ω). According to [31], each integral in the sum over k as an analytical function of ν has singularities in the form of simple poles at the points ν ∈ N. It is convenient to write the second contribution as Suppose the function, is finite and has finite derivatives with respect to ε at ε = 0 or, put another way, recalling the dependence of H(ω) on m 2 , we suppose that the average number of particles has finite derivatives with respect to m 2 . This is a rather strong restriction on the class of operators under consideration. In particular, this condition fails to hold for massless particles in the Minkowski space. That is why we shall calculate the partition function of massless particles proceeding to the limit m 2 → 0. If the condition mentioned fulfills, the contribution (16) as a function of ν has simple poles at the points ν ∈ N. In this case, the function I ν (µ)/Γ(1 − ν) is an entire function of ν.

Function
It is natural to associate with ζ + (ν, ω) another function σ −1 ν (ω), which appears in the high-temperature expansion. Introduce the function Henceforward, the symbol reg denotes the analytical regularization of the integral [31] with respect to the parameter ν. This function is analytic in the ω plane and has a branch cut discontinuity on the positive real semiaxis, where At large values of ω, from (10) and (18) the asymptotics follows , where the exponentiation is defined as x α := |x| α e iα arg x , arg α ∈ [0, 2π), i.e., it possesses the cut along the real positive semiaxis in the ω plane. This asymptotic behavior holds true for all Re ν < (d + 1)/2. Indeed, for Re ν ∈ ((d − 1)/2, d/2), differentiating (18) with respect to ω, we come to .
On integrating this expression and taking into account that Re ν < d/2, we arrive at (20) for the strip Re ν ∈ ((d − 1)/2, d/2). Proceeding further, we can extend the domain of applicability of the asymptotics (20) up to Re ν < (d + 1)/2. Taking into consideration the asymptotic behavior (10), one can see that the function σ −1 ν (ω) has the singularities for Re ν > d/2 only in the form of simple poles at ν ∈ N in the complex ν plane (the singularities of ζ + (ν, ω)) since the integral over ω ′ converges in that case. It also follows from the asymptotics (10) that the integral (18) (after the substitution ω ′ → 1/ω ′ ) as the function of ν possesses the poles at (d−2ν +1) ∈ N apart from the poles at ν ∈ N. If these new poles do not coincide with the poles at ν ∈ N then they are simple. Otherwise, the second order poles emerge. The asymptotics (20) confirms this observation. If the ω expansion of the stepless part of ζ + (ν, ω) in the vicinity of the infinite point has the form (25) and ζ + (ν, 0) = 0 then there are no additional singularities of the function σ −1 ν (ω) in the ν plane. By the use of the function σ −1 ν (ω), the integral (14) can be cast into the form where H is the Hankel contour. If one knows the expansion of σ −1 ν (ω) in the neighborhood of the infinite point of the ω plane convergent outside of the disk with the radius ω c 1 then the representation (22) allows one to derive the high-temperature expansion of I ν (µ) easily. For the sake of definiteness, suppose Im µ ∈ Figure 1: A schematic pattern of the contours H and H ′ on the complex ω plane for the bosonic case. Here crosses denote the poles of the Bose-Einstein distribution function (22). There is a cut along ω > ωc making the functions σ −1 ν (ω) and ζ 0 + (ν, ω) single-valued. As a matter of fact, the function σ −1 ν (ω) may have singularities in the interval [0, ωc), and the function ζ 0 + (ν, ω) may have singularities inside the disk |ω| < ωc. These singularities are not depicted.
(−π/β, π/β). Then we deform the contour H into H ′ so that |ω| > ω c on the contour and the contour itself lies beyond the disk of radius |µ| centered at µ (see Fig. 1). The former condition allows us to expand σ −1 ν (ω) in the vicinity of the infinite point, while the latter one permits us to expand each term of the series (see, for example, the asymptotics (20)) in terms of decreasing powers of (ω − µ).
Upon deformation of the contour H into H ′ , we get The last term in the expression for bosons is the contribution from the pole of the integrand (22) after the deformation of H into H ′ has been performed. Other poles do not contribute because of this deformation in the case when |µ| < π/β and |µ ± πi/β| > ω c , for fermions; |µ| < 2π/β and |µ ± 2πi/β| > ω c , for bosons.
The first inequality is the requirement that the two poles nearest to the point ω = µ do not lie inside the disk |ω − µ| < |µ|, while the second one is the requirement that both of them are outside of the disk |ω| < ω c . Substituting the development of σ −1 (ω) in the neighborhood of the infinite point to (23), expanding each term of the series in the inverse powers of (ω − µ), and then integrating termwise the series (see the similar integral (27) below), we arrive at the expansion in the increasing powers of β. The integrals determining the coefficients of this expansion are reduced to zeta functions. The high-temperature expansion converges provided the inequalities mentioned above hold. In case of need, the convergence radius of the high-temperature expansion can be increased by "blowing up" the contour H ′ and taking into account the additional poles (the leading Matsubara frequencies) of the integrand (23).

High-temperature expansion
In many cases, it is more convenient to use another representation of the function I ν (µ) for calculation of the high-temperature expansion. It is that representation which we shall use. Suppose that (cf. (10)) where ω is real and ω c is the boundary of the series convergence domain. Usually, it coincides with the boundary of the continuous spectrum, but this will be irrelevant to the derivation of a general formula for the high-temperature expansion. The function ζ 0 + (ν, ω) is the so-called stepless part of ζ + (ν, ω) (see [19,32,33,[42][43][44][45] and others) or the Thomas-Fermi approximation for ζ + (ν, ω). It is that part of ζ + (ν, ω) which is given by a naive heat kernel expansion (see, e.g., [16,32,33]). For the Laplacian type operators, the expansion of ζ 0 + (ν, ω) takes the form (25). The discreteness of quantum numbers is completely ignored in this contribution. The second contribution to ζ q + (ν, ω) is essentially quantum one. Generally, the infinite point of the ω plane is the transcendental branch point for the function ζ q + (ν, ω) at noninteger ν, i.e., the development of ζ q + (ν, ω) in the inverse powers of ω contains infinitely many terms at positive powers of ω. A typical example of the function with such a singularity is ω ν exp(−aω). Further, in Sec. 4.2, we shall see in a concrete example that exactly this type of functions arises.
At first, consider the contribution from the stepless part of ζ + (ν, ω) to I ν (µ). Suppose that the functions ω d−2ν−k in (25) are continued from the real axis to the complex plane as ω α = |ω| α e iα arg ω , arg ω ∈ [0, 2π). Then, for bosons, where the contour C c coincides with the contour H ′ for Re ω ≤ ω c . In the second equality, it is supposed that the restrictions (24) on µ are fulfilled, i.e., we can draw the contour H ′ as depicted in Fig. 1 such that the additional contributions from the poles do not appear. Also, in this equality, using the standard trick we have passed from the integration over the ray ω ≥ ω c to the integration over the contour H ′ − C c taking into account the branch cut discontinuity of the function ω d−2ν−k on this ray. Having performed the transformation (26), the high-temperature expansion of the first contribution in I 0 ν (µ) is readily evaluated. One can develop the function ζ 0 + (ν, ω) on the contour H ′ as series (25). Then, substituting this development to the integral, we obtain As regards the second and the third contributions in I 0 ν (µ), their high-temperature expansion can be obtained provided that |ω − µ| < 2π/β on the interval [0, ω c ] and on the contour C c . Subject to this condition, the function defining the Bose-Einstein distribution can be expanded in the Laurent series in β convergent on the integration contour. Keeping in mind that µ has been already constrained by (24), the condition stated results in one additional inequality |µ| + ω c < 2π/β, where it is implied that the contour C c can be deformed into a part of the circle |ω| = ω c . Obviously, all of the mentioned inequalities can be satisfied since β → 0. In this case, The expression forσ l ν (µ) can be rewritten in a more simple form, when Re ν > (d + 1 + l)/2. We deform the contour C c so that it runs from above and below the part of the real axis [ω c , +∞) and is closed by an arch of an infinite radius. Due to the condition on ν, the contribution from the arch vanishes and the integral over the contour adjoining the semiaxis [ω c , +∞) can be expressed in terms of the branch cut discontinuity of the function. Thenσ Similarly to how it was done for the function σ −1 ν (µ), one can prove that for Re ν > (d + 1 + l)/2 the function σ l ν (µ) has singularities in the ν plane in the form of simple poles at ν ∈ N. In case when Re ν ≤ (d + 1 + l)/2, there appears the poles at (d + l − 2ν + 2) ∈ N in addition to the poles at ν ∈ N. If the last condition is satisfied for some ν ∈ N then there are second order poles at these points. All the rest poles of the function σ l ν (µ) in the ν plane are simple. The analytic properties in the ν plane of the function coincide with the properties of the functionσ l ν (µ). It should be mentioned that if H(ω) depends on the mass squared such that ∂ m 2 H(ω) = −1 then the following recurrence relations hold In particular, the recurrence relation (91) of [19] can be deduced from these ones.
So, in the case of the Bose-Einstein statistics, we obtain the high-temperature expansion of I 0 ν (µ) (cf. [16,19]) If the series (25) converges for ω ≥ ω c , the convergence domain of the series (33) is determined by the restrictions on µ specified above. The analogous expansion for the fermions reads as where |µ| + ω c < π/β and the first condition in (24) should be met. It is instructive to compare the derived expansions with the ones presented in [16,19]. In fact, in these papers the analytical regularization (continuation) of the integral is considered instead of ζ + (ν, ω) and so is the functioñ instead of I ν (µ). This function is the entire function of ν subject to the restrictions imposed on the spectrum of H(ω). In this connection, it should be noted that an inaccuracy was made in formula (35) of [16]. The revised version of this formula is Nevertheless, this inaccuracy does not affect the rest formulas of the paper. The recurrence relation of the form (32) becomes and all the rest relations can be written in a similar fashion. For the stepless contribution from the antiparticles, we take into account that It is convenient to continue the functions ω d−2ν−k from the real axis to the complex plane just as it was done for the contribution from the particles. Then, for bosons, we have is an even function of µ. If ζ + (ν, −ω) = ζ + (ν, ω) then I ν (µ)+J ν (µ) is an even function of µ too. In particular, the Ω-potential (13), which includes the contributions from particles and antiparticles, is an even function of µ in this case. Now we are going to show how to take into account the essentially quantum corrections ζ q + (ν, ω) to the high-temperature expansion of I ν (µ). Let us assume that, for ω → +∞, the function ζ q + (ν, ω) is developed as a series with a typical term ω α(ν) e −aω , Re a ≥ 0.
Furthermore, α(ν) decreases as the term number increases and |a| > a 0 > 0. In this case, we have for bosons where one can put ω 0 to be equal to µ + 2π/β for real a and µ, |µ| < 2π/β. If a is a complex number then the integration contour has to be rotated so that it goes along the line of the steepest descent of the function e −aω . Then the last two terms in (45) are suppressed by the exponent e −aω 0 for β → 0. In what follows, such terms will be neglected. However, it should be noted that the expansion (45) becomes asymptotic after these exponentially suppressed terms have been cast out. The function I q ν (µ) possesses the same singularities in the ν plane as ζ q + (ν, ω) does and may have only simple poles at ν ∈ N. For fermions, the expansion can be obtained in a similar way, whereas, in this case, |µ| < π/β and ω 0 = µ + π/β for a real.
To sum up, if we neglect the exponentially suppressed terms at β → 0, the high-temperature expansion of I ν (µ) reads as (33), (34), whereσ l ν (µ) have to be replaced by σ l ν (µ). This holds true for the contribution from the antiparticles too.

Descent formulas
The main tool to derive the high-temperature expansions is the heat kernel expansion technique [24] and its various resummations [11,46]. As a rule, the coefficients of the high-temperature expansion depend explicitly on the Killing vector, which determines the stationarity of the space-time, singles out the privileged set of mode functions of the quantum fields, and, consequently, determines the Fock vacuum state. However, some coefficients turn out to be independent of the Killing vector [16,20,21,24,26,27]. One can convince oneself in this fact by a direct calculation, but it is easier to use the descent formulas [21] that relate the expansion coefficients of the heat kernel in (d+1) dimensional space to the expansion coefficients of the same heat kernel in d dimensional space. In deriving these formulas, it is supposed that the coefficients of the Laplacian type operator (the background fields) determining the heat kernel are independent of the x 0 variable. The method connecting the Green functions in (d + 1) dimensional space with the Green functions in d dimensional space is known in mathematical physics as the descent method (see, e.g., [47]). Thus, we shall call the formulas relating the heat kernel expansion coefficients in (d+ 1) dimensional space with the corresponding coefficients in d dimensional space as the descent formulas. Similar formulas can be found in [22,39,48].
Consider the operator of a Laplacian type, acting on the space of square-integrable functions depending on D := d + 1 arguments. Assume that the coefficients of this operator are independent of the variable x 0 , i.e., in particular, the Riemannian metric g µν associated with this operator possesses the Killing vector ξ µ ∂ µ = ∂ x 0 , Performing the Fourier transform over the variable x 0 and rewriting the resulting operator H(ω) in terms of the Killing vector, as done in (8) of [28] and (15) of [16], we arrive at where f (x) is some function independent of x 0 . Expanding the left-and the right-hand sides of this equation in τ and canceling the arbitrary function f (x), we come to We have resummed some terms of the heat kernel expansion into the exponent on the right-hand side as made in (18) of [16]. Writingã k (ω, x) in the form (see (38) of [16]) and integrating over ω in (49), we obtain the equality Whence, equating the coefficients at the same power of τ , we deduce the descent formula The left-hand side of this expression is independent of the Killing vector field and, consequently, the explicit dependence on the Killing vector field on the right-hand side cancels out. Such a cancelation happens for the metric of an arbitrary signature since it has a pure algebraic origin. The explicit form of the first three formulas read as Notice that the coefficient at the logarithm in formulas (44), (45) of [16] (see also (130) below) coincides with the coefficient a 2 [20,21,24,26,27]. Now we prove another descent formula [21]. Let us given Then, substituting these asymptotic expansions to (48) and equating the terms at the same powers of τ , we find [21] b whereā (j) k are determined in the same way asã where the left-hand side has to be expanded in a series in τ . The descent formula (55) follows from (52) at m = 0. The descent formulas (55) enable us to prove an interesting property of the high-temperature expansion. If we expand the terms in the free energy standing at the zeroth power of temperature and at the negative powers of the effective mass squaredm 2 in the inverse powers of m 2 then the coefficients of the expansion will be independent of the Killing vector field. The resulting asymptotic series in m −2 is the standard large mass expansion of the effective action at zero temperature. This may serve as an indirect verification of formulas (41), (42) of [16] for the contribution of the stepless part of ζ + (ν, ω) and, in general, of the correctness of the energy cutoff regularization scheme.
Indeed, according to formulas (41), (42) of [16], the terms described above are of the form where the ellipses denote the remaining terms. The same formula, but with an opposite sign, holds for fermions too. Being rewritten through m 2 , this expression reads Now, taking into account thatā (j) k = 0 for j odd, we can write this formula as where, in the last equality, we have exploited the descent formula (55). The last formula is the standard form of the expansion of the effective action at zero temperature in the inverse powers of a large mass. This expansion can be found, for example, in (6.42) of [23]. The difference in sign (−1) k results from the different definition ofā k (compare (54) with (6.39) of [23]). Note that formula (57) is correct for D odd. As far as even D is concerned, one ought to get rid of all the terms, which are singular due to the gamma function in the numerator, i.e., all the terms at the nonnegative powers of the mass. So, we see that if D is odd, the finite part of the effective action at zero temperature (without the non-perturbative corrections) being expanded in asymptotic series in the inverse powers of a large mass m does not depend on the Killing vector. If D is even, this is valid for the terms at the negative powers of m 2 and the logarithmic term only. Later on, when we derive the explicit expressions for the non-perturbative contributions to the high-temperature expansion, we shall turn back to the interpretation of this result.

Perturbation theory
Now we apply the general formulas obtained above to the concrete model. Consider a massive scalar field on a stationary gravitational background at a finite reciprocal temperature β. For simplicity, we restrict our considerations to the case of a vanishing chemical potential. In the adapted coordinates, where ξ µ = (1, 0, 0, 0), the free energy F takes the standard form where T µ ν is the energy-momentum tensor, Σ is the Cauchy surface, which we take to be x 0 = const. The operator H is the Hamiltonian of the scalar field expressed in terms of the creation-annihilation operators associated with the stationary mode functions (u α ,ū α ). These mode functions are the eigenvectors of the Lie derivative with respect to the Killing vector (see, e.g., [3,6,18]) viz., they depend on time as e −iωat in the adapted system of coordinates. The mode functions corresponding to the energy ω span the kernel of the Klein-Gordon operator, where all the time derivatives should be replaced by −iω. This operator, which we denote as H(ω), must be supplemented by the appropriate boundary conditions. To simplify further calculations, we assume that the system considered is large enough to neglect the boundary effects or the space represents a compact manifold without boundary. The operator H(ω) is of a Laplacian type, it is Hermitian with respect to the measure |g| on the square-integrable functions depending on x, and possesses the spectrum bounded from above at fixed ω.
The one-loop correction to the free energy can be cast into the form (12) with µ = 0. As we derived in the previous section, the high-temperature expansion is written as (33), (34) with the replacementσ l ν → σ l ν , when the exponentially suppressed contributions are neglected. The first terms in (33) and (34) (i.e., the terms in these formulas that are proportional to the product of the gamma and zeta functions) were found in [16,20,21,[24][25][26][27]32]. They are determined by a stepless part of ζ + (ν, ω) and can be obtained with the help of the asymptotic heat kernel expansion in the large mass m. In order to find the second terms of the high-temperature expansions (33) and (34), one needs a non-perturbative (that is not in the form of an asymptotic series in τ ) expression for the heat kernel exp(−τ H(ω)) taken on the diagonal.
Hamiltonian. The heat kernel is a mere evolution operator with the Hamiltonian H(ω) taken at the imaginary time −iτ . Therefore, in order to find the approximate, but non-perturbative (in the sense mentioned above), expression for it, we can employ the standard perturbation theory in quantum mechanics assuming that the coefficients of H(ω) are nearly constant. For the reader convenience and for the conformity of notation we describe such a perturbation theory in Appendix B in some detail. In order to use the ordinary formulas of quantum mechanics and to take completely into account the dependence of the heat kernel on the metric, we shall work with the operators self-adjoined with respect to the standard scalar product with a trivial measure (not |g|). To this end, one has to perform a similarity transform changing the measure and making it trivial. As a result, the operator H(ω) passes toH(ω) (see [28] for details, the Euclidean version see in [29]), where p i := −i∂ i , ξ 2 = g µν ξ µ ξ ν , g µ := ξ µ /ξ 2 , h µ := ∂ µ ln ξ 2 , and the relations (4) hold. The connection ∇ i is constructed by the use of the positive-definite metricḡ ij . Then the heat kernel and its trace become respectively. Henceforward, it will be convenient to change the sign of the Hamiltonian and regard −H as the generator of evolution instead ofH. In that case, the quantum mechanical evolution operator is exp(−is[−H]). Putting s = iτ , we obtain the heat kernel from the latter operator.
System of coordinates. According to the general procedure expounded in Appendix B, it is necessary to split the Hamiltonian (63) into the free part H 0 , quadratic in the variables x i (τ ) and p i (τ ), and the perturbation V . Besides, we demand that, for an every finite order of the perturbation series, the approximate evolution operator possesses all the symmetries of the exact evolution operator. In our case this means that the contributions of the perturbation theory are to be written in a general covariant form in terms of the metric g µν and the Killing vector ξ µ , which defines the vacuum state as described above. Since the perturbation theory will represent a certain expansion in derivatives of the coefficients of the operatorH, we need to define a covariant gradient expansion. With this end in view, we introduce the Riemann normal frame of the metricḡ ij with the origin at the middle of the geodesic of the metricḡ ij connecting the points x and y (the so-called midpoint prescription). This system of coordinates is not uniquely defined: apart from the global Euclidean rotations in space around the origin, one can change the time variable as t → t + ϕ(x), where ϕ(x) is some smooth function of x (see [34], Sec. 88, for details). Under the latter transform, the fields g i are changed by a gradient of the function ϕ(x) similarly to the electromagnetic potentials. Let us seize this opportunity to redefine t and impose the Fock gauge on g i : in the Riemann normal frame specified above. The tensor f µν is defined in (5). These conditions fix unambiguously (up to global space rotations) the system of coordinates in the spacetime. This allows us to restore in a unique way the general covariant expressions from the derivatives ofḡ ij and g i taken at the origin of the frame [49]. For instance, where we have employed the formulas for the developments ofḡ ij andḡ in the Riemann normal coordinates (see, e.g., [28,29,49,50]). The tensorR ijkl , its covariant derivatives, and contractions constructed with the help of the metricḡ ij are expressed through the curvature tensor R µνρσ of the metric g µν , the Killing vector ξ µ , and their covariant derivatives and contractions (see, e.g., Appendices in [13,16,28]).
Free Hamiltonian. Now we need to single out unambiguously the quadratic part H 0 of the Hamiltonian and define the power counting scheme for the diagrams. The quadratic part H 0 determines the base of the perturbation theory and its propagators, while the power counting scheme allows us to order the infinite set of diagrams and distinguish the most relevant ones under the assumption of smallness of the field derivatives. In many respects this procedure is analogous to the effective field theory approach [51,52], but slightly simpler as long as we consider quantum mechanics with a finite number of degrees of freedom. In order to introduce such a grading, we shall make, at first, some estimations of the typical magnitudes of the structures appearing in the expansion of the coefficients of the operatorH. Let L be a characteristic scale of variations of the gravitational field. For example, for a spherically symmetric metric, this quantity is of the order of the distance from the center of a gravitating object to the point where the derivative expansion is sought. In the weak field limit, at a large distance from the gravitating object, where ε : we can use the expressions given in [34], Sec. 105. As a result, we have for solutions to the Einstein equations. The estimation of f ij has been obtained for the maximal (critical) value of the angular momentum of a gravitating object J = M r g /2, where M is a mass of a body and r g is the Schwarzschild radius. In the strong field limit, where ε ≈ 1, i.e., near the ergosphere, we can use the Kerr solution to find the estimations. In this case, the most relevant contributions come from the terms containing the negative powers of ξ 2 and the derivatives acting on ξ −2 , viz., In the weak field limit, the condition of slow variation of the fields g µν and ξ µ turns into |ω|L ≫ 1.
As follows from (69), in the strong field limit, ε ≈ 1, this condition is substituted for i.e., the derivatives ofḡ ij and g i are made dimensionless with the aid of the appropriate power of ω. Further, we shall see that this is indeed the case. We specify the free Hamiltonian H 0 determining the base of a perturbation theory by imposing the following two requirements: 1. H 0 is no more than quadratic in x i (τ ) and p i (τ ); 2. H 0 is quadratic in ω and p i , and does not include the terms at lower powers of ω and p i (apart from the constant term m 2 ).
The first condition is necessary to construct the perturbation theory stated in Appendix B. It is concerned with the fact that, for the systems that do not possess some special symmetries, i.e., for a general background, a general solution to the Heisenberg equations can be constructed only for the quadratic Hamiltonians. The second condition is related to the requirement that the free Hamiltonian must include the most relevant terms in the short-wave approximation (71). As we shall see below (80), p i ∼ ω for the solutions to the Heisenberg equations. The term m 2 is taken into account non-perturbatively since usually m 2 L 2 ≫ 1. On the other hand, the inclusion of the terms at lower powers of ω and p i into H 0 is unjustified since these small corrections are overlapped by the succedent terms of the perturbation series 2 . The Hamiltonian H 0 containing all the terms satisfying the above two conditions is uniquely defined Recall that we changed the overall sign of the Hamiltonian. The first condition imposed on the Hamiltonian H 0 looks rather technical and requires some extra substantiations. First of all, recall that our chief goal is to obtain the generating functional of one-particle irreducible Green functions (the effective action) at a finite temperature and small external momenta. In constructing the effective action with the aid of the background field method, the background metric has not to be a solution of the classical equations of motion (the Einstein equations), in general. So as to find at small external momenta we only need a sufficiently wide class of metrics slowly varying in space. It is clear that the quadratic approximation described above works well for stationary metrics closely approximated by, for example, in the adapted system of coordinates. In the first case, the quadratic approximation gives the exact answer, whereas in the second case the fulfillment of the condition (71) is implied. It is also reasonable to expect that the quadratic approximation is good enough for the solutions to the Einstein equations, when the metric varies slowly (see the approximation [11] and its numerical verifications in [14,56] and others). In terms of the contributions of classical trajectories to the evolution operator (see [33] for details), the quadratic approximation accurately describes the two leading contributions: the contribution from the shortest geodesic of a given energy connecting the points x and y (this provides the Thomas-Fermi contribution to ζ + (ν, ω)) and the contribution from the geodesic of a fixed energy connecting the points x and y, reflected once from the turning point (this provides an oscillating contribution to ζ + (ν, ω)). For example, the approximation made should work well for the potential of two spherically symmetric gravitating bodies in the vicinity of the point where the attractive force is approximately zero. In the neighbourhood of this point, the gravitational potential along the line connecting the centers of the gravitating bodies has the form of an inverted parabola. Consequently, there are complex classical trajectories under the potential barrier, which are "wound" around this parabola. According to the general results of ( [33], see also [57]), the contributions of such trajectories are represented by oscillating exponentially suppressed terms in ζ + (ν, ω). Usually, in QFT, the particle creation process (in our case, the creation of particles by a gravitational field [3,4,58,59]) is related to such trajectories in the sense that the module of the matrix element between the vacuum state of a scalar field in a flat spacetime and the vacuum state defined with respect to the creation-annihilation operators associated with the stationary mode functions, which take into account the interaction with gravity, is less than unity. This fact can be revealed by the presence of imaginary terms in the effective action constructed in an appropriate way. Note that in the framework we develop, which is based on the representation of a free energy (6), (13), the imaginary contributions to the free energy are absent by virtue of the fact that (13) is real. However, if one uses the Schwinger representation where that integral is understood in the sense of analytical regularization over ν by analogy with the generalized function (x − i0) −ν [31], then the effective action at a finite temperature, formally constructed as (13) with the replacement of ζ + (ν, ω) by ζ(ν, ω), will possess imaginary terms. A more detailed investigation of this question will be given elsewhere.
Ingredients of the perturbation theory. The Hamiltonian H 0 determines the averages and propagators (177) in the interaction picture and the matrix element out|in , where the states |in and |out are specified in (170) and (171). The explicit expressions for these ingredients of the perturbation theory can be obtained for an arbitrary quadratic Hamiltonian. Nevertheless, to simplify the subsequent formulas we consider the case when The vectors υ 1 i ,ῡ 1 i , and υ 2 i are orthonormal with respect to the standard Hermitian scalar product, the overbar denotes complex conjugation, the vector υ 2 i having real components. Also, for definiteness, we assume that λ 1 < 0 and λ 2 > 0. In a weak gravitational field, this relations hold for a vacuum solutions to the Einstein equations (see (48) of [28]). In that case, we obtain [28] (for the Euclidean version see [29]) By construction, the above expression for the matrix element out|in is a bi-density with respect to the coordinates x and y. So as to obtain the bi-scalar, one has to multiply this expression by∆ where∆(x, y) is the covariant van Vleck determinant for the metricḡ ij . Notice that the formula (30) of [28] contains a misprint in the sign of the expression standing in the exponent: exp(−iω ± τ ) should be substituted for exp(iω ± τ ).
The averages and the propagators (177) are written as (see Appendix B) We shall depict the propagators D xx and the averagesx by solid lines on the Feynman graphs. The propagators D pp and the averagesp will be denoted by dashed lines, while for the propagators D xp the half solid half dashed lines will be used. On developing the Hamiltonian (63) as a covariant Taylor series and casting out the quadratic part H 0 , we can distinguish four types of vertices: where n is the number of x lines joined to the vertex and the dots denote possible additional solid lines. The vertices of the types V 2p and V 1p have no less than two x lines joining to the vertex, while the vertices V x possess no less than three such lines. The ordering of operators in the vertex is taken into account by the infinitesimal shifts of the time arguments of operators. The time arguments are shifted in such a way that the T -ordering places them in the proper order as they stand in the Hamiltonian (see, for instance, (201)). Let us introduce a grading on the set of diagrams. We attribute, formally, every vertex by the coupling constant α n , where n is the number of derivatives of the fieldsḡ ij and g i in the vertex. As an example, see (81) and the vertices of the order α 2 in (201). For the vertices V 2p , V 1p and V x , the power of α is equal to the number of the lines D xx and D xp joining by the x leg to the vertex. As for the V 0 type vertices, one should keep in mind the fact that the vertex without external lines has already the order α 2 . The order of the whole diagram in α is defined as the product of the "coupling constants" α n of all the vertices of the diagram. The order of the diagram in ω is defined as where E p is the number of the external p lines, while I pp and I xx are the numbers of the internal lines D pp and D xx , respectively. The numbers of the corresponding vertices are denoted by V x , V 2p , and V 0 . Formula (82) easily follows from the explicit expressions for the averages and propagators (80). This formula also allows for the fact that the integral over τ in the vertex produces extra ω −1 . It can be seen from that the proper-time τ enters the arguments of the trigonometric functions and exponents in the combination τ ω.
Having redefined τ → ω −1 τ , the arguments of the functions mentioned cease to depend on ω, and every integration over τ results in ω −1 after that. It follows from (82) that every closure of the external lines into a loop diminishes the order of a diagram in ω by one. Thus, in order to find the contributions to the connected part of the evolution operator matrix element (187), one ought to draw all the connected tree diagrams of a given order in α, which determines the order of the diagrams in derivatives. The tree diagrams give the leading contribution in ω. Then, the external lines are closed into loops, which results in the corrections to the tree contribution of the order of ω −L , i.e., the expressions will look like where f (x) is some function independent of s and ω specified by the Feynman rules, while L is the number of loops. The loop expansion is equivalent to the quasiclassical procedure for expansion in ω −1 or elaborated in [60,61].
In the paper [28], it was verified by the explicit calculation that the development in s of the evolution operator ensuing from the perturbation theory described above reproduces the standard asymptotic expansion of the heat kernel (see, e.g., [22]). As a matter of fact, in order to reproduce all the terms of the asymptotic expansion of a given order n in derivatives, one has to calculate all the contributions of the perturbation theory up to the order α n . The resultant perturbation theory is rather cumbersome even for the diagonal matrix elements of the evolution operator. In Appendix B, we provide all the terms of the perturbation series of the order α 2 for the logarithm of the heat kernel diagonal (see (209)-(218)). Nevertheless, such a perturbation theory, as it is formulated in Appendix B, admits a simple realization in a computer program. The only technical issue, which could arise, is the evaluation of the integrals in vertices. However, as seen from (80), in our case the integrals of the perturbation theory are reduced to the integrals of a product of exponents, which are easily calculated analytically (see formulas (204) and (206)).

Non-perturbative corrections
The purpose of our investigation is to derive the explicit expression for the divergent and finite parts of the high-temperature expansion of a free energy. For a four dimensional spacetime, it implies we need to find all the contributions of the perturbation theory for the heat kernel with fourth derivatives of the fields, i.e., up to the order α 4 . The prospects are rather ominous having in view the explicit expressions for the terms of the order of α 2 (see (209)-(218)). Fortunately, as we shall see below, the higher contributions of the perturbation theory are not necessary for our aim. In fact, a good approximation for the high-temperature expansion to the order we consider can be obtained by a mere combination of the results of the papers [16] and [28]. It is important at this point that the higher orders of perturbation theory do not change the arrangement of singularities of the evolution operator in the s plane (see (80)).
According to the general results of the previous sections (33), (34), and (45), so as to find the hightemperature expansion, it is sufficient to obtain the explicit expressions for the functions σ l ν defined in (31) provided we neglect exponentially suppressed terms at β → 0. The first terms in the expansion (33), (34) were found in [16] up to the required order in β. To simplify the calculations and to adjust the notation to [16], we shall use the normalization (36), (37). At coinciding arguments, the Gaussian contribution (78) to the diagonal of the heat kernel can be cast into the form where Notice that the misprint was made in formula (56) of [28]: one has to exchange th ↔ cth in this formula. The branch of the multivalued function (84) in the complex τ plane is specified by the conditions that the function G 0 (τ ) should be holomorphic in the strip Re τ ∈ (−π/ω √ 2λ 2 , 0), it should be real-valued on the segment of the real axis belonging to this strip, and the system of cuts is to be symmetric with respect to the reflection in the real axis. Consequently, the structure of singularities of the function G 0 (τ ) looks as depicted in Fig. 2 and the square root of the sine in (84) is defined as x −1/2 = |x| −1/2 e − arg(x)/2 , where arg x ∈ [0, 2π). Later on, it will be convenient to continue (84) analytically to the complex ω plane. Then we shall define ω 3/2 as an analytic function with the cut along ω < 0, i.e., the principal branch of the function ω 3/2 will be taken.

Massive case
Partition of the contour. Let us consider, at first, the massive case m 2 > 0, mL ≫ 1. We shall assume that the singularities of the functions G(τ ) and G 0 (τ ) are located at the same points of the τ plane (this is the case at any finite order of the perturbation theory) and have the same "structure". Also we shall suppose that the nearest to the origin singularity lies on the real axis of the τ plane and rather than on the imaginary axis, that is we assume Usually, this inequality is satisfied (see (48) of [28]) in the weak field limit. Then the integral determiningζ + (ν, ω) is rewritten as where C = C 1 ∪C 1 and the contour C 1 goes from the point i/(m ge g s ω) to the point −i/(m ge g s ω) passing the origin of the τ plane from the left. Henceforth, we use the convenient notation and also g s := (g µ g µ ) 1/2 . Notice that m 2 ge and m 2 gm are not expressible in a covariant way in terms of the metric g µν , its curvature, and their contractions. This can be easily checked by considering m 2 ge and m 2 gm at the origin of the system of coordinates we work. The expressions for m 2 ge and m 2 gm do not depend onḡ ij , whereas any scalar possessing the dimension m 2 constructed solely in terms of the metric does depend on g ij and its derivatives of the second and higher orders taken at the origin of this frame.
In virtue of the assumption on the structure of singularities of G(τ ), the heat kernel can be developed as a convergent series in τ on the contour C 1 : where we have saved out the exponential factor from the series. The expression standing in the exponent gives the major contribution in the short-wave approximation and is determined by this condition unambiguously. Also if we redefine the variable ω → mω and take into account that mL is a very large quantity then, at the leading order, we come to the expression standing in the exponent again. So it is that form which is suitable for the large mass or the short-wave expansions. As a result, the contribution from the contour C 1 becomes The contribution of the first term in the square brackets provides the stepless part ofζ + (ν, ω). It is this contribution which was found in [16]. The second term in the square brackets leads to an appearance of the essentially quantum oscillating contributions (44) toζ + (ν, ω). These contributions and the contributions coming from the second term in (86) were not taken into account in [16]. Notice that the series in (89) consisting of the contributions from the first term in the square brackets does not converge. It is an asymptotic series. Only does the inclusion of the exponentially suppressed at m 2 → ∞ terms (the contributions of the the second term in the square brackets) make the series (89) convergent. A wrong impression may also arise that the contribution of the second term in the square brackets in (89) cancels out the second term in (86). This is not the case as the order of the summation over k and integration over τ are not interchangeable in (89). Often, carrying out considerations on the so-called physical level of rigor, one disregards such mathematical "subtleties". However, they are relevant in our case. Further, we shall ascertain that these two contributions toζ + (ν, ω) and, correspondingly, to the free energy are represented by different expressions.
Contour C 1 . Thus, we can distinguish three types of terms inσ l ν : the two contributions (89) from the contour C 1 and the contribution (86) from the contourC 1 . The first contribution coming from the integration along the contour C was found in [16]. Denoting byτ l ν the contribution of the stepless part of ζ 0 + (ν, ω) tõ σ l ν , we have (see (33) of [16]) where j ≤ [4k/3]. In the last equality, we have expanded the expression in m 2 instead ofm 2 . The introduction of the effective mass allows us to simplify drastically the calculation of higher terms of the heat kernel expansion, but the assignment of a physical sense to it is not always warrantable (see the footnote on p. 15). So, the two other contributions toσ l ν are left to calculate. Let us start with the second term in (89). Recovering the explicit dependence of the coefficients of the heat kernel expansion on ω as in (50), we see that we need to take the integral, in order to find the contribution toσ l ν . The contourC + 1 goes to the point i/(m ge g s ω) and the contourC − 1 goes from the point −i/(m ge g s ω) so that the integrals over τ converge. Consider the first integral. Rotate the integration contours in the ω and τ planes in such a way that during this rotation the integral remains convergent and its value is left unchanged. In the ω plane we rotate the integration contour clockwise by the angle π/2, while in the τ plane we rotate the contour counterclockwise such that it runs along the real axis from −∞ to −1/(m ge g s |ω|) after the rotation. Recall that the function τ k+ν−1−d/2 has a branch cut discontinuity for τ > 0. The integration contour in the τ plane does not intersect this branch cut during the rotation. Thus we come to where Γ(α, x) is the incomplete gamma function. The last integral can be also obtained by a straightforward calculation. At that, it is convenient to add −i0 to the terms standing in the round brackets in the exponent. Analogously, for the second integral, we rotate counterclockwise the integration contour in the ω plane by the angle π/2, while in the τ plane we rotate the contour clockwise such that it goes along the real axis from −∞ to −1/(m ge g s |ω|) after the rotation. Then we have Adding the above two expressions and performing a change of the variable, we arrive at Remark that this expression is equal to zero for (l + j) odd. Moreover, sinceā (j) k = 0 for j odd, this contribution vanishes for l odd. Also we see from this formula that the contribution considered is essentially quantum one, because the asymptotics of the expression standing under the integral over ω (in fact, the contribution toζ + (ν, iω)) has the form (44). Now we employ the asymptotic expansion for the incomplete gamma function at large arguments (see, e.g., [62]), and make the substitution ω = e t . Then the integral is easily evaluated by the WKB method in the vicinity of the saddle point t = 0. As a result, we obtain for the leading contribution Collecting the factors remaining in (89), we write the contribution of the second term in (89) toσ l ν as σ l ν II ≈ −e iπν sin π(l + 1) 2 ∞ k,j=0 where j is an even number and j ≤ [4k/3]. The expression (97) is regular at ν = 0.
ContourC 1 and the singular points. Now we have to evaluate the contribution toσ l ν from the second term in (86). In calculating this contribution, we shall use the Gaussian approximation (84) to G(τ ). Usually, the contributions we are about to calculate are small in comparison with the main contribution toσ l ν from the stepless part ofζ + (ν, ω). So, the use of G 0 (τ ) for the evaluation of such contributions is justified. The higher corrections to G 0 (τ ) like those found in (209)-(218) give even smaller contributions. The analysis of the asymptotics of G 0 (τ ) at large τ shows up that for the integration contourC 1 in the τ plane can be rotated to the right, whereas for it can be rotated to the left. Inasmuch as we are interested in the high-temperature expansion at a large mass m and ω ∼ m after the stretching of the integration variable, it is convenient to rotate the contour to the right, when and to the left, otherwise. This condition is independent of m 2 after the integration variable has been stretched ω → mω. It is clear that the final answer,ζ + (ν, ω) andσ l ν , does not depend on a choice of the concrete condition on ω, when the contourC 1 is to be rotated to the right or left. It is only important that the integral over τ will be convergent.
When we rotate the integration contourC 1 to the right, the singular points τ = in/(m gm g s ω), where n ∈ Z, also contribute to the integral over τ . We shall consider the contribution of these points subsequently, but now we find the leading contribution toσ l ν from the integral along the contourC 1 after the rotation to the left or to the right. As seen in Fig. 2, this integral is saturated in the neighbourhood of the boundary points ±i/(m ge g s ω). At these points, we have Here S = S 0 + τ m 2 is the expression standing in the exponent in (84). The preexponential factor becomes at these points .
Now we ascertain that, in the weak field limit, the corrections of order α 2 found in (209)-(218) give smaller contributions to S than the terms written above. For the reader convenience, we provide here the orders in ε for the structures appearing in the development of S. In the weak field limit, it follows from (68) that At the boundary points ±i/(m ge g s ω), we have t ∼ ℓ ∼ 1 and so where we have depicted the types of the diagrams only. The contributions of the terms proportional to ∂ b f ca are even smaller. The corrections (104) should be compared with the terms in (101). It is evident the contributions (104) can be neglected. Further, in order to find the contribution toσ l ν , we act in the same way as we did in calculating the contribution (97). For the part of the integration contourC 1 lying in the upper half-plane of the complex τ plane (C + 1 ), we make a transform ω → e −iϕ/2 ω, putting eventually ϕ = π, and the contour in the τ plane is deformed to part of the negative real axis (−∞, −1/(m ge g s |ω|)]. For the part of the integration contour C 1 lying in the lower half-plane of the complex τ plane (C − 1 ), we perform the analogous rotation, but to the opposite direction. The integral over τ is evaluated by the WKB method for the boundary point. Therefore, we can safely forget about the singular points on the imaginary axis when evaluating this integral. Summing up the two expressions corresponding to the lower and upper parts of the contourC 1 , we obtain in the leading order Note that this contribution is essentially quantum one as long as the spectral density corresponding to it behaves as (44). The mention also should be made that this contribution vanishes for l odd and, in particular, for l = −1 (the contribution to the coefficient at β −1 in the high-temperature expansion for the bosonic case). The resulting integral over ω is saturated near the saddle point ω = m/ĝ provided mĝ/(m ge g s ) ≫ 1.
Supposing that this inequality is fulfilled, we arrive in the leading order at where the last equality has been obtained under the assumption that ε ≪ 1 and the estimations (68) are satisfied, i.e., in the weak field limit. In this limit, one should set g s = 1. Nevertheless, we have kept the explicit dependence on g s to save the covariance of the expression under dilatations of the Killing vector: ξ µ → λξ µ , λ = const. Obviously, the expression obtained is finite at ν = 0.
In the massive case, it remains to find the contribution of the singular points τ = in/(m gm g s ω), when the inequality (100) holds. First of all, we redefine the integration variable τ as τ → ω −1 τ . Then, the singular points pass into τ = in/(m gm g s ). The integrals in the vicinity of these singular points will be evaluated by the WKB method, which is justified in the case when the inequality (71) is satisfied (see (107) below). The integration contour is depicted in Fig. 2. Also, in order to simplify the calculations, we assume that ε ≪ 1 and the estimations (68) are fulfilled. Then the expression standing in the exponent (84) is well approximated in the neighbourhood of the singular points by the three leading terms of the Laurent series in the vicinity of these singular points (see Fig. 2). In particular, such a truncated expansion provides a good approximation for the integrand in the neighbourhood of the saddle points where the integral is saturated. The preexponential factor is replaced by the leading term of the Laurent series in the vicinity of the singular point. Thus we have where x := τ − in/(m gm g s ), and all the terms of the order of ε and higher are thrown away. As seen from this expansion, the variable x is of the order of unity in a neighbourhood of the saddle point. The preexponential factor inσ l ν expanded in a neighbourhood of the singular points reads as e iπν e ±iπ(3−2ν)/4 (−1) n 16π 2 i where the upper sign corresponds to positive n's, the lower sign is for negative ones, and the first term of the Laurent series is only retained. The integral over x is easily evaluated As we see, this contribution toζ + (ν, ω) is essentially quantum one. From the formulas (107), (109), and (110) we deduce that a pair of singular points corresponding to ±n, n > 0, makes a contribution toσ l ν of the form Usually b ⊥ and, consequently, B n are rather small (see Eq. (48) of [28]). Therefore, we shall consider separately two cases, when (a) B n = 0, and (b) Recall that λ 1 < 0. The origin of the conditions (112) becomes clear soon.
In the case (a), the major contribution to the integral over ω comes from the boundary point. The sine argument also possesses the stationary point, but its contribution is exponentially suppressed at m/m gm ≫ 1 as compared to the contribution of the boundary point. The derivative of the sine argument calculated at the boundary point takes the form As a result, making use of the standard formula for the leading contribution from a boundary point and neglecting the corrections of the order of ε in the preexponential factor, we havẽ In the case (b), let us employ the asymptotic formula for the Bessel function at large arguments [62]: Then, substituting ω = mg −1 s ch ψ, the integral over ω is reduced to ∞ ψ 0 dψ sh 1/2 ψ ch 3/2+l−ν ψ cos(2B n m sh ψ − π/4) sin mn m gm (a n ch ψ − 1/ ch ψ) where ψ 0 := arcch(g/g), a n :=g It is evident that the trigonometric functions in (116) are highly oscillating. To analyze the saddle points, we write them via the exponents. This shows that the saddle points are located at the extremum points of the function mn m gm (a n ch ψ − 1/ ch ψ) ± 2B n m sh ψ.
The plot of this function is given in Fig. 3. This function has the extrema at the points In the latter case the signs are consistent with (118). It follows from the assumptions (112) that the extremum at the point ψ = m ge B n /n falls into the integration interval over ψ (the first condition in (112)) and is sufficiently sharp (the second condition in (112)). In that case, the contributions from the boundary point and other saddle points are small in comparison with the contribution of the extremum point ψ = m ge B n /n (see the integration contour in Fig. 3). In particular, the contribution to the integral over ψ from the exponent with (118) (multiplied by an imaginary unit) taken with the plus sign is small as compared to the contribution of the analogous exponent, but with the minus sign in formula (118). Developing (118) as a series in ψ up to the second order and keeping only the terms of the order of ε 1/2 and lower, we come to (the minus sign is taken) Evaluating the Gaussian integral, the contribution toσ l ν becomes The total contribution of the singular points τ = in/(m gm g s ) toσ l ν is obtained by summing the above expressions (114) or (121) over all natural n. The arising series in n, converges for any ν ∈ C. At ν = 0 and m ge /m gm > 1, the series is rapidly convergent and well approximated by its first term. Notice that, in deriving the expressions (114) and (121), we neglected the exponentially suppressed terms as against the leading contribution. These disregarded terms may be of the order of or even larger than the exponentially suppressed terms in (97) and (106). We are only interested in the leading non-perturbative contributions. The corrections (209)-(218) can be estimated as follows. In the vicinity of the singular points τ = in/(m gm g s ω), the most relevant contributions come from the terms in the t and ℓ tensors that are the most singular at these points. For the tree and one-loop contributions containingR (acb)d , they are proportional to x −2 , while for the two-loop contributions with the V 2p vertex they are proportional to x −1 . Hence, we have Taking x at the extremum point of (107), substituting ω = mg −1 s ch ψ, and expanding the resulting expression in ψ, we find where the terms at most quadratic in ψ have been only kept. These contributions should be compared with the terms in (120) at the same powers of ψ. We see that the contributions (124) are negligible in the weak field limit for a large mass m. As for the rest contributions of the order α 2 with the vertices V 1p and V 0 , they are much smaller than those we have considered.
High-temperature expansion. Thus,σ l ν can be written as where the items are presented in (90), (97), (106), and (122). It is only the first term inτ l ν which possesses singularities in the complex ν plane. The mention should be made that the non-perturbative contributions (97), (106), and (122) inσ l ν cannot be developed as a Laurent series with a finite principal part in the inverse powers of m 2 . The point m 2 = ∞ is a transcendent branch point for these contributions. Having cast out the exponentially suppressed contributions at β → 0, the high-temperature expansion takes the form where β T := ξ 2 β is the Tolman reciprocal temperature, and the limit of a vanishing ν is implied. In accordance with the general properties of the high-temperature expansion investigated in Sec. 2 and the direct calculations [16], F b is an entire function of ν. Besides, as long as e iπν can be factored out in front of the expression regular at ν → 0, this factor can be omitted in proceeding to the limit ν → 0. Also note that the first and second terms in (126) taken separately possess the poles in the ν plane. However, these divergencies are mutually canceled out. For the "scalar fermions" the analogous expansion reads as In fact, one just needs to add the non-perturbative corrections we found to the high-temperature expansions given in Sec. 7 of [16]. At that, the formulas of [16] should be rewritten in terms of m 2 rather than the effective mass squaredm 2 . In particular, lnm 2 must be expanded in the inverse powers of the mass squared as Then the explicit expressions for the finite and divergent parts of the high-temperature expansion presented in [16] are left unchanged save lnm 2 , which should be replaced by ln m 2 .
One of the interesting consequences of the result we have obtained is that we explicitly separated the perturbative and non-perturbative contributions to the effective action at zero temperature, the convergent (in terms of the Feynman diagrams) perturbative corrections standing at the negative powers of the mass squared being independent of the Killing vector and expressed through the metric only. Indeed, the one-loop nonrenormalized contribution of one bosonic degree of freedom to the effective action at zero temperature can be obtained from (127) by the use of the formula [16,19] where β −1 plays the role of a cutoff parameter and T is the time interval tending to infinity. Upon cancelation of singularities in the ν plane, the divergent and finite parts of the high-temperature expansion of F f for D even have the form (see (42) of [16]) where ψ(x) := Γ ′ (x)/Γ(x) and γ is the Euler constant. The prime at the sum over n says that the singular terms are discarded, the second term is zero by definition when the argument of the factorial is negative, and the last term in curly brackets vanishes by definition when the gamma function entering the numerator tends to infinity. In Sec. 3, we proved that the coefficient in front of ln β in (130) and the terms at the negative powers of m 2 (the last term in the curly brackets in (130)) do not depend explicitly on the Killing vector and are expressed solely in terms of the metric. The formula (130) allows us to understand why we do not "see" the dependence on the Killing vector when evaluating the effective action using the standard perturbation theory on a flat background. The sum of all the diagrams presented on the left picture in Fig. 4 corresponds exactly to the terms in the curly brackets in (130). Despite the fact that the terms at the nonnegative powers of m 2 does explicitly depend on the Killing vector, they can be completely canceled out by the appropriate counterterms added to the initial action (see [16] for details). In the framework of the dimensional regularization, this cancelation happens automatically, whereas for the other regularization schemes this cancelation is achieved by requiring the fulfillment the Ward identities for the vertex functions. The terms standing at the negative powers of m 2 correspond to the convergent contributions of the perturbation theory, are unambiguously evaluated, and can be represented in a general covariant form in terms of the metric. On the other hand, the last three terms in (130) cannot be reproduced at any finite order of the perturbation theory on a flat background. In order to reveal these non-perturbative terms, one has to sum the asymptotic series of perturbation theory in some way. This, in essence, is equivalent to the introduction of the background field. In a certain sense, one may say that the effective action is nonanalytic in the vicinity of the point g µν = η µν .
One could think, keeping in mind the analogy with the vector gauge fields in the presence of their sources, that the components of the metric g 0µ can be expressed through the components of the energy-momentum tensor with the help of the constraints. Obviously, even if such a situation took place, the problem of the explicit dependence on the Killing vector would not be solved since, in this case, we would also need some extra structure distinguishing the time components of the energy-momentum tensor. One could also hope that the non-perturbative contribution found would be canceled out if we took into account the contributions from the graviton sector (see Fig. 4). However, this is not the case in the one-loop approximation. For a vanishing background scalar field, the non-perturbative corrections depend nontrivially on the mass m, whereas the one-loop contributions from the graviton sector are independent of the mass of a scalar field. The non-perturbative corrections obtained cannot be canceled by the counterterms either, inasmuch as these corrections describe the actual physical phenomena similar to the Landau oscillations at a nonzero chemical potential [42-45, 63, 64] or the vacuum polarization effects [65,66] in quantum electrodynamics. Also recall that, in the standard model, m 2 is expressed in terms of the SU (2) doublet φ of the Higgs field. The addition of the counterterm nonpolynomial in φ 2 violates the renormalizability of the Higgs sector (even with the inclusion of the nonminimal term R|φ| 2 ). Thus we conclude that the effective action at zero temperature (its non-perturbative part, at least) depends explicitly on the Killing vector and cannot be written solely in terms of the metric.

Massless case
Main contribution. Let us turn now to the massless case m 2 = 0. In this case, the inequality (100) is always satisfied and we can deform the integration contour as depicted in Fig. 2 Tr e −τH(ω) .
Then it is necessary to evaluate the contributions from the contours C 0 , C ge , and C gm . We begin with the main contribution C 0 . Substituting the heat kernel expansion, we obtaiñ The integral over τ is the incomplete gamma function γ(α, x). Therefore, where we have expandedā k (ω, x) in ω as in formula (50). The integral over ω converges in the strip of the complex ν plane This allows us to find unambiguously the functionσ l ν | C 0 in the form of an analytic function of ν. Making use of the formula [62], which is easy to prove integrating by parts, we arrive at Contour C ge . We shall calculate the contributions of the contours C ge and C gm by the steepest descent method supposing that the estimations (68) and the inequality (71) are fulfilled. Let us start with the contour C ge . Upon stretching the integration variable, τ → ω −1 τ , the expression S 0 standing in the exponent in (84) is written as in the neighbourhood of the singular points. Here, x = τ − (2n + 1)/(m ge g s ) and we have discarded all the terms of the order ε and higher. The extremum points are now readily found The value of x ± is of the order of ε 0 . The contribution of the two extremum points x ± becomes where and we have only kept the leading contribution in ε in the preexponential factor. As we see, the contribution of the contour C ge toζ + (ν, ω) is essentially quantum one and has the form (44). The integral over ω converges for Re ν < (l + 2) and is easily taken. Theñ where Im acts on the expression in the square brackets just as the operation of taking the imaginary part for ν real. The total contribution of the contour C ge is obtained by summing (141) over all natural numbers n: The mention should be made that, in contrast to the massive case, the essentially quantum contributionσ l ν possesses the "infrared" poles at ν = l + 2, l + 3, . . .. Proceeding to the limit of a vanishing ν and taking into account that the second term in the round brackets under the Im sign in (141) is small, we derive in the leading order dx |g| 4πg l s |b | λ 2 m ge 2n + 1 l+4 m gm (l + 2)! sh(π(2n + 1)m gm /m ge ) where we have introduced the function The special function Li ν (z) is the polylogarithm which is defined as Using the Poisson summation formula, it is no difficult to find the asymptotic expansion of φ µ (ν, z) for z small: where, in the first line, the summation is carried over all the collections p = {p k }, p k = 0, ∞, k = 1, ∞.
The integral over x, for those values of µ and ν where it diverges, should be understood in the sense of an analytical continuation in these parameters. In particular, for µ = 1, Contour C gm . It is left to evaluate the contribution toσ l ν from the contour C gm . From formula (111) with m 2 = 0, we have for the contribution of a pair of singular points This contribution is evidently quantum one. The integral over ω is reduced to the Gauss hypergeometric function ([62], 6.699) for 0 < b < a and − Re λ < Re ν < 3/2. In our case, the integral over ω on the right-hand side of (148) equals where a n is given in (117). We see that the integral vanishes at ν = 0 and l even. In particular, the contribution of the contour C gm to the finite part of the effective action (l = 0, ν = 0) is equal to zero. The vanishing of this contribution at even l is the exact property of the heat kernel (84) for m 2 = 0, i.e., it holds not only for the weak field limit that we are considering now. Further, we shall be interested only in the leading contribution in ε. In this limit, a n and the hypergeometric function equal to unity. Then, (148) is written asσ and the contribution of all the singular points on the imaginary axis of the τ plane takes the form Notice that the oscillations are not present inσ l ν . This is a characteristic feature of a massless case at zero chemical potential (see, e.g., [19]).
High-temperature expansion. In order to obtain the high-temperature expansion, we just need to substitute the expressions (136), (143), and (152) to the general formula (126) or (127) and cancel the poles at ν → 0, if they arise. For even D such singularities do appear in the first terms of (126) and (127) and in the contribution (136) of the contour C 0 toσ l ν . Consider separately these two contributions to the free energy. Bearing in mind thatā (j) k = 0 for j even and j ≤ [4k/3], we get for the divergent and finite parts of the high-temperature expansion in the bosonic case where the omission points denote the terms at β in higher powers. The primes at the sum signs recall that all the singular terms arising at certain values of the indices k, j, and s are thrown away. The term in the second line is deliberately written in such a form to make a comparison of (153) with (41) of [16] simpler (the values of the digamma functions are presented in (43) of [16]). Matching (153) with formula (41) of [16], it is easy to see that the first term in (153) is the same as the first term in (41) of [16] expressed through m 2 and then m 2 set to zero. The coefficient at the logarithm in (153) coincides with the same coefficient in (41) of [16] at m 2 = 0. The expression in the argument of the logarithm is obtained by the replacement m 2 → m 2 ge e −2γ /4. The second term in the second line of (153) equals two times the expression in the third line in (44) of [16] expressed via m 2 (see (B21) of [16]) at m 2 = 0, the terms at the negative powers of m 2 being absent here. The coefficient at β −1 T and the last term in (153) are the expansions in derivatives of the metric and the Killing vector field. These expansions do not coincide with the analogous terms in (41) of [16] on identifying m 2 ge → m 2 (it is not surprising, of course). In particular, the coefficients at m 2 ge in the last term of (153) are not expressible solely in terms of the metric. In the fermionic case, the formula analogous to (153) can be cast into the form Recall that the coefficientsā (j) k are expressed through the coefficientsã (j) k found in [16] (see (56)). The total high-temperature expansion without the exponentially suppressed terms at β → 0 becomes where the non-perturbative contributions are presented in (143) and (152).
In conclusion, we outline how the procedure above is changed when m ge < m gm . In that case, the expansion in τ of the heat kernel G(τ ) is convergent in the disc |τ | < 1/(m gm g s ω). Therefore, it is convenient for both massive and massless cases to part the integration contour C in (6) into two: C 1 passing from the singular point τ = i/(m gm g s ω) to τ = −i/(m gm g s ω), and the contourC 1 being the rest part of the contour C (see the structure of singularities in Fig. 2). The integral along the contour C 1 is evaluated in the same way as we did for the integrals along the contours C 1 in the massive case and C 0 for the massless one. The integration contourC 1 is rotated in an appropriate way subject to the sign of the expression standing on the left-hand side of (100). Then the integral over τ is calculated by the steepest descent method just as we did for the contourC 1 and the contributions of the singular points on the imaginary axis of the τ plane. A detailed study of the expressions resulting from this procedure will be given elsewhere.

Discussion
The results that we obtained are valid within the bounds of the approximations made, when the Klein-Gordon equation on a curved background holds and the basic principles of QFT are kept. Using the background field gauge [6,7], one may expect on the basis of the perturbative analysis on a flat background that, in the infrared limit, the effective action Γ[g µν ] is expressed in a covariant form in terms of the metric, its curvature, and covariant derivatives constructed with the help of this metric. However, as we saw, this can be done only for the perturbative contributions after a suitable renormalization procedure (see the details of this procedure in [16,18]). Non-perturbative corrections depend explicitly on the Killing vector, which defines the vacuum state and, consequently, the representation of the algebra of observables according to the GNS construction. This dependence is just a manifestation of the dependence of averages on the state with respect to which they are calculated. Of course, there could be a situation that the unique vacuum state is singled out by a certain natural local structure made of the metric g µν alone. In particular, for a stationary gravitational background this would mean that the Killing vector, or some its analog, should be constructed in terms of the spacetime curvature and its covariant derivatives. However, there is no such a natural construction (see, e.g., the discussion in [1]). The Killing vector field is not expressed through the metric in a local form. Furthermore, if one tries to construct QFT on the background metric, which is "produced" by some extended massive object (say, a star) that at the initial moment moved uniformly and rectilinearly, then was accelerated by some external force, and after that became moving uniformly and rectilinearly, but with other velocity, then one comes almost inevitably to the conclusion that, under the assumption of locality and causality of QFT, the field ξ µ must be dynamical.
The non-perturbative corrections we found are extremely small far from the ergosphere. For example, the leading non-perturbative contribution to the energy density in the massive case (121), (122) is of the order ofσ On the other hand, there is another suggestive example in QFT. It is quantum chromodynamics (QCD) with the spontaneously broken chiral symmetry. At the present moment, the well established proof (starting from the QCD Lagrangian) of the existence of spontaneous breaking of the chiral symmetry and the description of the quark condensate properties are not given. However, the very assumption that this symmetry is spontaneously broken leads to a great number of nontrivial phenomenological consequences formalized in chiral perturbation theory (see, e.g., [52]). In particular, as long as the non-perturbative corrections to the effective action for gravity depend explicitly on the Killing vector, this vector field (or g µ = ξ µ /ξ 2 ) and its derivatives must appear in the initial Lagrangian in accordance with the general prescriptions of the effective field theories [51,52] and renormalization theory (see, e.g., [67,68]). Therefore, a complete cancelation of the perturbative quantum contributions to the effective action depending on the Killing vector by the appropriate counterterms is unnatural, at least. At the same time, one should bear in mind that an immediate transfer of reasonings used in constructing chiral perturbation theory is illicit. For instance, the Goldstone theorem does not apply in the case we consider since it rests on the assumption that an open system tends to a minimum energy state (the ground state). It is clear that this principle does not work in general relativity so long as the energy operator is an illdefined concept till the symmetry is broken, i.e., until a certain external structure defining the representation of the algebra of observables and the energy operator is introduced.
In general relativity it is natural to demand the fulfillment of the principle of general covariance from the in-in effective action Γ[g ± µν , Φ ± , . . .] (in the background field gauge) for all the fields including those defining the vacuum state and the representation of the algebra of observables. The omission points in the argument of the effective action denote these fields. The calculations presented in this paper and many others [3, 4, 11-16, 18, 28, 56] suggest that the role of such fields is played by some quantum timelike vector field ξ µ , or the Tolman temperature one-form g µ constructed in terms of this vector field. For stationary spacetimes, the average of this field should coincide with (or be close to) the Killing vector field. The requirement of general covariance appears to be sufficient to determine the dynamics of this field [15]. In this case, the so-called background independence is nothing but the existence of the in-in effective action Γ[g ± µν , Φ ± , g ± µ ] satisfying the Ward identities following from the general covariance: where the approximate equality means the fulfilment of this equality on the solution to the equations of motion for the fields Φ + and g + µ . The same equality holds for the "minus" fields as well. The vacuum state of quantum fields and all their correlators are restored from Γ. The dynamics of the fields g ± µ are found from the Ward identities (157), i.e., in fact, from the self-consistency condition for the quantum Einstein equations. In particular, the equations of motion for the average field g µ look like (157), but with the "plus" and "minus" fields identified upon variation. Note that, for a stationary spacetime, this last equation possesses the solution coinciding with the Killing vector [15], although other solutions to (157) may also exist in this case. Therefore, the Killing vector field, which is usually taken to define the vacuum state, does not violate the Ward identities for the average field (see, e.g., [13]). The higher Ward identities, i.e., (157) without identification of the "plus" and "minus" fields, can be satisfied only if the field g µ is quantized.
As shown in [15], the equations (157) for the field g µ are the equations of a hydrodynamical type. It is interesting to note that the sound speed in this "fluid", i.e., the speed of propagation of small disturbances of the field g µ , is determined by the asymptotics of the energy-momentum tensor in the weak field limit [15,16,28] and independent of the details of the model. Of course, inasmuch as the field g µ interacts with the metric g µν , whose perturbations propagate with the velocity of light, the disturbances of the metric induce the perturbations of the field g µ that propagate with the velocity of light too, as for the usual fluid in a gravitational field. The quantization of the field g µ can be constructed employing the well-known relativistic Lagrangian for an Eulerian fluid (see, e.g., [69][70][71][72][73][74][75]).
It would be interesting to investigate the quantum dynamics of the field g µ in the weak field limit, when the dispersion law for its small perturbations is known, and obtain possible phenomenological restrictions on its vertices. The immediate generalization of the results of the paper is, of course, the development of a procedure, similar to the one presented in the paper, for higher spin fields at a nonzero chemical potential including the interaction with the background fields other than metric.
In the case 1, the system is unstable at a = 1 (since ω k (1) is complex) and N is an even number. Indeed, if the roots remain complex for a ∈ (a 0 , 1] then N = 0. If they returned to the real axis and then go out to the complex plane then, as follows from the properties of ω 1,2 k (a) discussed above, N must be even. In the case 2, the system is unstable since ω 1 k and ω 2 k lie on the one side from the point ω = 0 and ε ′ (ω 1 k )ε ′ (ω 2 k ) < 0, which contradicts (11). The case 3 is only left. Consequently, N is an odd number provided the vacuum is stable at a = 1, and so ε k (0, 1) < 0 for any k. Thus we have shown that ζ + (ν, 0) = 0 subject to the above mentioned assumption on a smooth deformability of H(ω, 0) into H(ω, 1) and the vacuum stability at a = 1.
Define the in-and out-states in the Heisenberg representation and the interaction picture as |in : = U 0,t in |in, t in , |out : = U 0,tout |out, t out , |in : = S t in ,0 |in = U 0 0,t in |in, t in , |out : = S tout,0 |out = U 0 0,tout |out, t out . (171) Then out, t out |U tout,t in |in, t in = out|S tout,t in |in = out|in , and the Heisenberg averages calculated in the standard in-out perturbation theory take the form out|T {φ(t n ) · · · φ(t 1 )}|in = out|T {φ I (t n ) · · · φ I (t 1 )S tout,t in }|in .
In order to construct the perturbation theory, let us define the generating functional of the free Green functions If H 0 (t, φ) is quadratic in the fields φ µ (this will be assumed henceforward), the evolution operator in the last line can be explicitly calculated. It is not difficult to show that Therefore, varying (174), we obtain where out|in φ µ := out|φ i I (t)|in , out|in D µν : One of the simplest ways to take the noncommutativity of the operators in the vertex V (φ) into account is to shift their arguments by infinitesimal quantities so that under the T -ordering they would stand in the proper order. This resolves the ambiguity of the expressions like D ij (t, t) arising in the perturbative computations. The Feynman diagram technique is just a graphic representation of the right-hand side of (178). The easiest way to get it is to use the last formula in (178). From now on, we follow [77]. The following relation holds This formula is obviously valid for the functions F k = exp(p k φ k ) (no summation). All the other functions that can be Taylor expanded are obtained by a differentiation of k exp(p k φ k ) with respect to the parameters p k , these parameters vanishing after a differentiation. The differentiations with respect to φ k and p k commute. The formula (179) implies that where Developing the exponent as a series, we come to where π = {π kl , k ≤ l}, π kl = 0, ∞, k, l = 1, ∞, and S int k ≡ S int (φ k ). The right-hand side of the last formula can be represented graphically as a sum of all possible graphs with n numbered vertices, where π kl defines the number of lines 3 coming from the vertex k to the vertex l, the vertices with L lines correspond to the expressions and the lines correspond to iD µν , the summation over all the identical Greek indices is implied. The coefficient at the diagram equals It is easy to see that this quantity does not depend on the way of numeration of the vertices. On identifying φ k = φ, the vertices (183) become independent of the number k. Collecting all the identical diagrams with unnumbered vertices, we see that the coefficient at every the unnumbered diagram is given by (185) multiplied by n!/s, where n! is the number of permutations of vertices in the diagram and s is the number of permutations of the numbered vertices in the graph that do not affect it (the symmetry factor). Thus, we have the following Feynman rules for the right-hand side of (178): where π kl , k ≤ l, is the number of lines coming from the vertex k to the vertex l for some fixed (arbitrary) numeration of the diagram vertices. The computation of higher orders of the perturbation theory can be simplified if we observe that where "conn.part" means that only the connected diagrams with their coefficients are left in the expression. The connected diagram is a graph where one can pass moving along the lines from any vertex to any other vertex. A graph consisting of a single vertex is connected by definition. Let us prove the formula (187). If where c a is the numeric coefficient (186) of the connected diagram W a , then where the sum is performed over all the possible collections m = {m a , a = 1, ∞}, m a = 0, ∞. On the other hand, let us gather all the identical disconnected diagrams in the expansion of the right-hand side of (178) containing the diagrams W 1 , . . ., W n as the connected parts that are replicated with the multiplicities m 1 , . . ., m n . The symmetry factor for such a diagram is where s Wa is the symmetry factor for the graph W a . The factorial in this formula corresponds to the number of permutations of the numbered vertices in m a copies of the graph W a (all the vertices of one copy are permuted with all the vertices of the other). Other factors in the coefficient (186) at the diagram fall into the product of multipliers corresponding to the connected components W a . Therefore, the coefficient at such a diagram is equal to where c a is the coefficient (186) of the diagram W a . Comparing the derived expression with (189), we arrive at (187). Sometimes it is useful to make certain components of the fields φ µ → φ i a (t) explicit, where the index a numbers the components. Then the Feynman rules can be cast into the form where π ab kl is the number of lines of the type ab connecting (directly) the vertices k and l. The symmetry factor s is calculated in the same way, but taking into account that the lines are of a different type now. The statement about the connectedness of the diagrams defining (187) remains valid. Also, in some cases, it is convenient to single out the "mean field"φ explicitly and expand all the expression in it. It follows from the first formula in the second line of (178) that, in this case, the additional ingredient in the Feynman rules arises external line =φ µ , and the coefficient at the diagram acquires the extra factor where π k is the number of external lines joining to the vertex k. With this definition, the free endpoints of the external lines must not be considered as vertices in calculating the symmetry factor s. If the free endpoints of the external lines are regarded as vertices, the factor (194) ought not to be placed. It is reproduced automatically in evaluating s. It is clear that the statement about the connectedness of (187) holds true even after the inclusion of external lines.
In order to simplify the subsequent formulas, it is convenient to use the basis specified by a tetrad where its vectors are defined in (77). It is useful to introduce the following tensors in this basis where σ(r) = ±, r = 1, n + k, the sum over σ denotes the sum over 2 n+k combinations of σ(r), the bar over σ(r) implies the sign change of σ(r), n is the number of indices of x, and k is the number of indices of p.
The quantities b a ± ,κ a , andf a are the eigenvalues of the matrices (199). The tensor (204) is recovered in the initial basis as t i 1 j 1 ···i n+k j n+k = a 1 ,...,a n+k t a 1 ···a n+k x···xp···p e i 1 a 1ē i 1 a 1 · · · e i n+k a n+kē i n+k a n+k .
The tensor t arises in calculations of the tree one-vertex diagrams. In fact, this tensor equals the half of the integral over σ from −1 to 1 of the product of a corresponding number of the "first" terms (containing exp(ib ± σ)) in the expressions forx andp in (197). For the one-vertex diagrams involving tadpoles, the integral over σ will contain the product of several sines and cosines coming from the propagators with the coinciding arguments (200), in addition to the factors fromx andp. Therefore, it is also useful to introduce the tensors ℓ a 1 ···a n+k b 1 ···b m+l x···xp···ps···sc···c := (−i) k+m 2 m+lκ a n+1 · · ·κ a n+k sinκ b 1 · · · sinκ b m+l × × σ,ρ σ(1) · · · σ(n)ρ(1) · · · ρ(m) sin where ρ(q) = ±, q = 1, m + l, the sums are over all 2 n+k+m+l combinations of (σ(r), ρ(q)), m is the number of indices s, and l is the number of indices c. Furthermore,κ ± := ±κ. In the initial basis, the tensor is constructed as given in (205). The obvious properties hold for the ℓ tensor The tensors t and ℓ are symmetric with respect to permutations of the indices in each of the group x, p, s, and c. Here are some specific values: With the help of these tensors it is not difficult to write all contributions of the perturbation theory of the order α 2 . Exploiting (197), (200), and (202), we arrive at The coefficients at the diagrams are calculated making use of the formulas (192). The summation is implied over all the tetrad indices. Note that four terms (see (202)) correspond to each of the diagram on the left apart from the last one. However, in virtue of the symmetry properties of the tensorsR acbd and ∂ b f ca , some terms are similar or vanish. For instance, the contraction of the fields with the indices a, c or b, d for V 2p and the indices a, c for V 1p vanishes in the loop diagrams. The total contribution to the logarithm of the diagonal of the evolution operator G(ω, s; x, x) of the order α 2 is given by the sum of the expressions (209)-(218).