We compute the heat kernel coefficients that are needed for the regularization and renormalization of massive gravity. Starting from the Stueckelberg action for massive gravity, we determine the propagators of the different fields (massive tensor, vector and scalar) in a general linear covariant gauge depending on four free gauge parameters. We then compute the non-minimal heat kernel coefficients for all the components of the scalar, vector and tensor sector, and employ these coefficients to regularize the propagators of all the different fields of massive gravity. We also study the massless limit and discuss the appearance of the van Dam–Veltman–Zakharov discontinuity. In the course of the computation, we derive new identities relating the heat kernel coefficients of different field sectors, both massive and massless.

In both mathematics and physics the heat kernel technique is a well-established method for the computation of traces of differential operators.1–6 Its area of application ranges from fluctuations of quantum fields on curved spacetimes, the determination of ultraviolet divergences, effective actions and quantum anomalies to quantum gravity and further. Concretely, it allows to define the regularized propagator and using this the one-loop effective action, and to compute counterterms and anomalies in a way that is naturally extended to field theories on curved space. Furthermore, it constitutes an essential ingredient in solving the gravitational functional renormalization group equation.7,8

The interest in the heat kernel can be traced back to the beginning of the theory of quantum fields in curved spacetimes. In the 1950s, Morette and DeWitt established the path integral approach to propagators in curved spacetime, and in the 1960s Schwinger and DeWitt developed the proper time formalism and the associated heat kernel expansion. The main idea is to express a given Green’s function as an integral of the heat kernel (which satisfies the heat equation with time τ, the proper time) over the proper time. The heat kernel depends on the background geometry, especially the background metric, and possibly other background fields such as gauge fields. Except for some special manifolds, such as maximally symmetric spacetimes, it is not possible to compute it exactly in general. However, it is possible to compute it approximately, and up to now two main expansion schemes have been developed: the small-time or local expansion (also known as the Seeley–DeWitt technique), and a non-local expansion for small curvature. The first scheme is employed in quantum field theory to compute ultraviolet divergences and anomalies, since these are directly determined by the lowest-order heat kernel coefficients in the small-time expansion. The non-local expansion is used instead to compute the finite part of the effective action in a covariant manner. In the framework of quantum gravity, non-local heat kernels are crucial in renormalization group studies, such as tests of the Asymptotic Safety conjecture, and heat kernels on maximally symmetric spaces are used to find β functions for gravitational couplings.9–12 A recent development in this framework is the consideration of Lorentzian renormalization group (RG) flows,13–18 which is important since a generic metric does not admit a Wick rotation. In general, including couplings with increasing scaling dimension such as higher curvature terms, an increasing number of heat kernel coefficients is needed, whose computation becomes increasingly difficult.

For gauge theories of Yang–Mills type it is straightforward to determine classical observables using the Becchi–Rouet–Stora–Tyutin (BRST) formalism.19–21 Since these generally are not linear in the elementary fields, in the quantum theory they become composite operators which need additional renormalization, beyond the usual renormalization of couplings.22 This can also be done using the heat kernel technique, point-splitting the classical expression and subtracting as many terms of the small-time expansion of the heat kernel as are needed to obtain a finite coincidence limit (depending on the scaling dimension of the operator). When using an effective action, one has to introduce an additional term coupling the composite operator to an external source, renormalize the extended effective action, and then take functional derivatives with respect to the source to obtain correlation functions including insertions of this operator. In the exact renormalization group framework, this has been recently used to study flows of volume and length operators in gravity,23–25 see also Refs. 26 and 27. For example, when studying the correlator of metric fluctuations in the vielbein formalism, the metric itself is a composite operator which needs to be regularized and renormalized in a suitable way, such as the point-splitting method28–30 or the heat kernel renormalization method explained above.

However, when defining observables in gravity another complication arises. Namely, in contrast to Yang–Mills gauge theories where the gauge symmetry is an internal symmetry and acts locally, i.e., it transforms fields at the same point, the diffeomorphism symmetry of gravity moves points around. It follows immediately that a local field cannot be gauge-invariant, and that observables must be necessarily nonlocal. Various approaches to construct such nonlocal observables have been considered, for example defining correlation functions involving the geodesic distance between points31–42 or relational observables. Nevertheless, once one has suitable observables at hand, also for gravity the usual framework of regularization and renormalization applies, and a natural question that arises is the gauge dependence of those regularized and renormalized quantities in the quantum theory. While the gauge independence of the observables is clear in the classical theory, in the quantum theory anomalies may arise which then might render the final result gauge-dependent. To check the gauge independence of the renormalized observables, it is useful to work in a general family of gauges that depend on a certain number of gauge parameters, and then verify that the end result is independent of these parameters. For this, it is necessary to know the Green’s functions and the heat kernel coefficients in such a general family of gauges. In Ref. 43, the authors constructed Green’s functions in a globally hyperbolic spacetime in general linear covariant gauges, both for vector gauge bosons and linearized Einstein gravity in the presence of a cosmological constant. In this work, we will generalize their result to the massive gravity case, relating the propagators of the massive gravity field sectors in a general linearized gauge to their expressions in the corresponding Feynman gauge.

The choice of the massive gravity theory is motivated from both a phenomelogical and cosmological point of view, and for theoretical interest. Modifications of Einstein gravity, of which massive gravity theories are a subset, give rise to an effective cosmological fluid which can explain dark matter and/or dark energy.44–46 On the other hand, while there are stringent upper bounds on the possible mass of the graviton,47,48 a tiny mass is not ruled out completely, even in light of the van Dam–Veltman–Zakharov (vDVZ) discontinuity.49,50 On the theoretical side, to the best of our knowledge the propagators for a massive spin-2 field in a general gauge have not been determined, and consequently also the corresponding heat kernel expansion is unknown. To avoid breaking diffeomorphism invariance, which a simple Fierz–Pauli mass term51–53 would do, we use the Stueckelberg formalism54,55 and add additional spin-1 and spin-0 fields that transform in such a way to keep diffeomorphism invariance of the complete theory. That is, we consider in this work only the free massive gravity theory, which determines the propagators. Moreover, we work with an on-shell background satisfying the background Einstein equations, which means that the propagators that we derive could, for example, later on be used to determine an on-shell effective action for gravity.56 While this is not strictly necessary, since one could also work with the off-shell (Vilkovisky–deWitt) invariant effective action,57 the resulting expressions are already complicated enough for an on-shell background, and would become unmanageable off-shell.

This article is structured as follows: In Sec. II we give some more background on massive gravity, determine its linearized action including the Stueckelberg fields that are necessary to preserve diffeomorphism invariance, and using the BRST formalism add the gauge-fixing and ghost terms that determine a quite general linear covariant gauge for the fields. Section III is devoted to the determination of the propagators in this gauge, which comprises three free gauge parameters, as well as their massless limit, and the verification of the Ward identities in the free theory. In Sec. IV, we give an overview of the small-time expansion of the heat kernel, compute the leading three heat kernel coefficients, and determine relations between the heat kernel coefficients of different spins. Section V comprises our main results, the heat kernel expansion of the propagators for all fields in massive gravity, including the massless limit. Finally, in Sec. VI we discuss our result and present an outlook for future work.

Conventions: Our conventions are a mostly plus metric, the Riemann tensor defined by ∇μνvρ − ∇νμvρ = Rμνρσvσ, and the Ricci tensor R μ ν = R ρ μ ρ ν . We work in n dimensions, use geometrical units = c = 1, and set κ 16 π G N Pl with Newton’s constant GN and the Planck length Pl.

Even though General Relativity (GR) has been very well tested and successfully describes gravity both a small and large scales,58 there are also some tensions. Of those, let us mention the observed galaxy rotation curves59,60 and primodial inflation,61–64 which cannot be explained using GR alone. Various models have been proposed to explain these effects and alleviate the tensions, of which dark matter–matter that only interacts gravitationally, or only very weakly with Standard Model particles — can explain the galaxy rotation curves and dark energy — an effective fluid with negative pressure — can explain inflation.64 While dark matter can be easily added as a new fundamental particle (or several) such as axions, and the current expansion can be explained by a cosmological constant, primordial dark energy is usually modelled by a single fundamental scalar field (the inflaton) with a suitable potential, see Refs. 65 and 67 for recent reviews. Other fundamental explanations of dark matter and dark energy come in the form of modified gravity models,68 of which massive gravity in its various incarnations45,69–72 is a popular one.

In massive gravity, the mediator of gravitational interactions, the spin-2 field usually called graviton, is given a small mass. From a particle physics perspective, this is a natural modification since GR can be thought of as the unique theory of a massless spin-2 particle;73–75 it also seems natural since we know that the carrier particles of the electroweak forces acquire a mass through the Higgs mechanism. Since the force mediated by massive particles falls off exponentially at large distances, a graviton with a small mass of the order of the Hubble rate mH could also mimic the effects of the current accelerated expansion of our universe without invoking a cosmological constant or other forms of dark energy. Another motivation for infrared modifications of gravity comes from brane world models,76,77 where extra dimensions can be of large or infinite size and the effective four-dimensional graviton propagator on our brane can behave like a massive one. We note that giving a mass to a free spin-2 field is easy, and was already done in 1939 by Fierz and Pauli;51–53 the construction of interacting theories of massive gravity is much harder. One possibility for a non-linear generalization of the Fierz–Pauli mass term was explored in Ref. 78, where an additional fully dynamical “reference” metric was introduced, rendering the model invariant under general coordinate transformations. Such theories are called bi-metric theories of gravity, and have been of particular interest due to the presence of accelerating cosmological solutions. Observational data from both the early evolution of the universe and solar system tests of gravity can be used to constrain the parameters of these theories.79–82 Another possibility is topologically massive gravity83,84 and its generalizations,85 which was the first interacting model of massive gravity in three dimensions. Its β functions can be computed in the asymptotic safety approach, and it turns out that as Einstein gravity in four dimensions,9,10 topologically massive gravity is asymptotically safe.86 Since Lorentz-invariant models of massive gravity often suffer from the presence of ghosts (fields with a wrong sign of the kinetic term), a possibility which recently also attracted attention is to allow for a violation of Lorentz invariance.87–91 From the cosmological point of view, this is particularly interesting because one can construct a consistent model of gravity where the tensor graviton mode is massive, while the linearized equations for scalar and vector metric perturbations are not modified. The Friedmann equation then acquires an effective dark-energy component, leading to an accelerated expansion, while gravity in the solar system is not modified.

However, in this work we are interested in the linearized Fierz–Pauli theory that preserves Lorentz invariance, and which can be obtained as the first approximation of a general massive gravity theory around a given fixed background. It turns out that even in the massless limit, the extra polarizations of a massive spin-2 field do not completely decouple in general, unlike what happens for an Abelian gauge theory. Instead, a scalar mode remains coupled to the tensor modes, and the resulting theory differs from the completely massless one, a effect that is known as the van Dam–Veltman–Zakharov (vDVZ) discontinuity.49,50,92,94 Nevertheless, whether decoupling occurs depends strongly on the concrete theory, which in our case means the background around which we study the theory. Let us mention here that the discontinuity is absent when one considers GR with a non-vanishing cosmological constant Λ. While for an anti-de Sitter background with Λ < 0 the resulting theory is ghost-free for all masses,94,95 in n-dimensional de Sitter spacetime with Λ = (n − 1)(n − 2)/2H2 > 0 the theory has a scalar ghost if 0 < m2 < (n − 2)H2; the upper bound (n − 2)H2 is known as the Higuchi bound.96 For a modern account of these issues we refer the reader to Refs. 97–99 and references therein.

Even though Fierz–Pauli theory preserves Lorentz invariance, the mass term breaks linearized diffeomorphisms, which are a symmetry of the massless theory that descends from the full Einstein–Hilbert action for gravity. To remedy this, we employ the Stueckelberg trick54,55 and add auxiliary fields that transform in such a way as to keep the full action invariant under linearized diffeomorphisms. The original massive Fierz–Pauli theory can then be recovered in a special gauge, the so-called unitary gauge. While both perturbative and non-perturbative computations in gravity have used heat kernel expansions (see for example Ref. 100 for a recent study), and also massive vector fields have been studied in great detail,101–103 to our knowledge the present work represents the first computation of the heat kernel coefficients for a massive theory of gravity in a general gauge.

We consider the Einstein–Hilbert action for gravity including a cosmological constant, and add Stueckelberg fields to make the theory massive while keeping gauge invariance.54,55 Working in a perturbative approach, we separate the full metric g ̃ μ ν into a background metric gμν and perturbations hμν according to
(1)
where κ 16 π G N is our perturbation parameter. Expanding the Einstein–Hilbert action to second order in hμν, we obtain
(2)
where Λ is the cosmological constant. The background equations of motion are obtained by varying SEH with respect to hμν and setting h to zero, which results in the well-known Einstein equations
(3)
When these are satisfied, we say that the theory has an on-shell background or simply is on-shell; the terms quadratic in the perturbation hμν in the action (2) then simplify to
(4)
A straightforward computation using the background Einstein equations (3) shows that this action is invariant under the gauge transformation
(5)
which arises from the linearization of the full diffeomorphism symmetry δ ξ g ̃ μ ν = L ξ g ̃ μ ν around the background (1) after the rescaling of the gauge parameter ξμκξμ.
To obtain the (free) massive gravity theory, we add the Fierz–Pauli mass term51–53 
(6)
which is however not invariant under the gauge symmetry (5). To preserve the symmetry, we add a compensating Stueckelberg field Aμ that transforms as δξAμ = ξμ, and replace hμνhμν − ∇μAν − ∇νAμ in the Fierz–Pauli action (6). To obtain a smooth limit m → 0, we further introduce another (secondary) Stueckelberg field ϕ and a gauge symmetry
(7)
and use the combination Aμμϕ instead of Aμ. Note that ϕ is inert under linearized diffeomorphisms. Taking all together, we obtain the extended Fierz–Pauli–Stueckelberg action
(8)
which is invariant under both the gravitational gauge symmetry (5) and the new gauge symmetry (7). The unitary gauge, in which we recover the original theory, is achieved by making a gauge transformation with parameters ξμ = −Aμμϕ and f = −ϕ, which sets both the new Aμ and ϕ to zero.
For vanishing cosmological constant Λ = 0 and in the limit m → 0, the different sectors (tensor hμν, vector Aμ and scalar ϕ) actually decouple. To see this, we follow Ref. 69 and rescale Aμm−1Aμ, ϕm−2ϕ. This results in
(9)
and the rescaled gauge transformations
(10a)
(10b)
Using now the on-shell equations of motion (EOM) (3) and adding the second-order Einstein–Hilbert action (4), we obtain the full action that we consider in this work:
(11)
and for Λ = 0 the last two terms (which would be problematic as m → 0) vanish, as well as the mass term for the vector Aμ. If we then take m → 0, the last two gauge transformations in (10a) vanish and we are only left with the original gravitational transformation (5) and the usual electromagnetic U(1) gauge transformation. Moreover, the full action (11) can be rewritten as
(12)
with the new field h ̃ μ ν h μ ν 2 / ( n 2 ) g μ ν ϕ , which is just a sum of the free actions and thus the tensor, vector and scalar sectors decouple in this limit. However, we will keep both Λ and m2 non-zero throughout the remainder of this work, and in fact define
(13)
to shorten the expressions. Note that λ is a dimensionless parameter, such that the mass m is the only dimensionful parameter that is left.

To quantize the Stueckelberg gauge theory of (free) massive gravity with action (11), we use the standard BRST formalism,19–21 where we introduce ghost, antighost and auxiliary fields for every symmetry. The gauge transformations (10a) are replaced by the action of the fermionic BRST operator s whose action on fields is obtained by replacing the gauge parameter by the corresponding ghost field, and the action on the ghost, antighost and auxiliary fields must be determined such that the action is nilpotent: s 2 = 0 .

We therefore need a diffeomorphism ghost cμ and a Stueckelberg ghost c. We also define their antighosts c ̄ μ and c ̄ , and introduce for each an auxiliary field Bμ and B. Then the BRST transformations on fields and ghosts are obtained from Eq. (10a) and read
(14a)
(14b)
(14c)
(14d)
In our case, the ghosts have trivial BRST transformations since we consider a free theory (but on a non-trivial background) such that the gauge transformations (10a) are linear. The antighosts and auxiliary fields correspondingly form trivial pairs:
(15a)
(15b)
To obtain the gauge-fixing and ghost terms, we add a BRST-exact term to the action. We use the gauge conditions
(16)
for the graviton and Stueckelberg fields, which are the most general ones compatible with the scaling dimension of the fields and a regular limit m → 0. We thus set
(17)
which after performing the BRST transformation reads
(18)
and
(19)
We see that the shifted auxiliary fields Bμ + ξ−1Fμ and B + α−1G decouple, and that the remainder of the gauge-fixing action (18) gives the required quadratic gauge-fixing terms. Since the gauge conditions (16) couple the different sectors, also the ghost action (19) couples the diffeomorphism and Stueckelberg ghosts. To determine the propagators of the fields and ghosts with reasonable effort, we thus make two simplifications: we work with the on-shell background (3) and the on-shell action (11), and make a choice for some of the gauge parameters, namely
(20)
We recall that λ = 4/(n − 2)Λ/m2 (13), such that the on-shell background has Rμν = λ/2m2gμν. The full action for the fields including gauge fixing now reads
(21)
where we defined
(22a)
(22b)
and the ghost action is
(23)
We see that the various sectors are still coupled. For a full decoupling, we would have to make a different gauge choice (the analogue of Feynman gauge), for which
(24)
and then use the new field h ̃ μ ν as for the limit (12).104 However, to verify the gauge dependence of the on-shell heat kernel coefficients, we have to keep at least some of the gauge parameters arbitrary, and the choice we made above is the one that is actually feasible.

To determine the propagators of the different fields in a general gauge, we have to relate them to the propagators in Feynman-type gauges, where the corresponding equation of motion is normally hyperbolic. In normally hyperbolic equations of motion, second-order derivatives only appear via the d’Alembertian ∇2, and their (retarded and advanced) propagators or Green’s functions can be constructed in an arbitrary globally hyperbolic spacetime using successive approximations.105–107 As we will see below, in a general gauge where second-order derivatives appear also in a different form the propagators can be obtained as a linear combination of the Feynman-gauge propagator and derivatives of propagators for fields of lower spin. For this, we will first derive various identities for the Feynman-gauge propagators, then treat the ghost sector where a spin-1 propagator in a general gauge is needed, and afterwards consider the fields sector where in addition also the spin-2 propagator in a general gauge appears. To derive these identities, we will make repeated use of the fact that the solution of a normally hyperbolic equation with retarded or advanced boundary conditions is unique, in particular that such an equation with vanishing source only has a vanishing solution. Lastly, we verify that our propagators satisfy the relevant Ward(–Takahashi–Slavnov–Taylor) identities.

We recall that the scalar propagator of (squared) mass M2 fulfills the EOM
(25)
For the massive vector propagator, the Feynman gauge action reads
(26)
and the Feynman gauge massive vector propagator G FG , M 2 thus fulfills the equation
(27)
Here and later on, primed indices refer to the point x′ and unprimed ones to x, and we have introduced the bitensor of parallel transport g ν μ ,108,109 which is defined as the unique solution to the equation
(28)
For the on-shell background with Rμν = λ/2m2gμν that we are using, Eq. (27) reduces to
(29)
and we see that there is a mismatch between the explicit mass M2 of the propagator and the effective one M2 + λ/2m2 appearing in the equation.110 
Taking the divergence of the equation of motion (29) for the Feynman gauge vector propagator and commuting covariant derivatives, we obtain
(30)
where the last equality follows from properties of the parallel transport bitensor [see Eq. (13.3) in Ref. 111], and combining with the scalar equation of motion (25) we get
(31)
Since ∇2M2 is a normally hyperbolic differential operator, which has unique (retarded and advanced) solutions, we obtain the identity109 
(32)
and upon exchanging x and x′ also
(33)
For the massive tensor propagator, the Feynman gauge action is the tensor part of the total massive gravity action (21) with gauge parameters ξ = ζ = α = 1 and reads
(34)
with
(35)
However, this EOM operator is not normally hyperbolic, since the coefficient of ∇2 is not (a multiple of) the identity on symmetric rank-2 tensors. Instead, we consider the operator
(36)
which is normally hyperbolic. Therefore, it has unique Green’s functions which can be constructed using standard methods,105–107 and we can assert that there exists a (rescaled) trace-reversed Feynman-type gauge propagator G μ ν ρ σ TG , m 2 satisfying
(37)
The Feynman-type gauge tensor propagator is then obtained by reversing this operation:
(38)
and satisfies
(39)
We now need to derive the trace and divergence identities for this tensor propagator. The easier one is the trace identity, for which we contract the propagator Eq. (39) with gμν to obtain
(40)
and from this
(41)
The differential operator is again normally hyperbolic, such that we have the unique solution
(42)
and upon exchanging x and x′ also
(43)
These trace identities generalize the known ones for the massless spin-2 field.112 
Using this relation, we can rewrite the propagator Eq. (39) as
(44)
Contracting with ∇μ and commuting covariant derivatives, from this formula we obtain
(45)
where we also used the contracted second Bianchi identity
(46)
whose right-hand side vanishes for our on-shell background. On the other hand, using the equation of motion (29) for the Feynman gauge vector propagator of mass M2, we compute
(47)
such that we infer
(48)
Finally choosing M2 = (1 − λ)m2, it follows that
(49)
which is the massive generalization of the known identity in the massless case.112 Exchanging x and x′, we also obtain
(50)
For the trace-reversed Feynman-type gauge propagator G μ ν ρ σ TG , m 2 (38), the trace and divergence identities (42), (43), (49), and (50) simplify to
(51a)
(51b)
(51c)
(51d)
To determine the propagators of all fields in massive gravity, we start with the simpler ghost sector. If we introduce the composite index A, we can define cA = (cμ, c) and c ̄ A = ( c ̄ μ , c ̄ ) and write the ghost action (23) as
(52)
with the components of the equation-of-motion operator PAB given by
(53a)
(53b)
(53c)
(53d)
The propagator
(54)
satisfies
(55)
where a sum over the composite index B is implicit, or in components
(56a)
(56b)
(56c)
(56d)
We make the general ansatz
(57a)
(57b)
(57c)
(57d)
where GA/B/C/D are scalar propagators of various masses, or linear combinations thereof, and MA is an unspecified mass. Plugging the ansatz (57) into the equations of motion (56) and commuting covariant derivatives (using that on-shell we have Rμν = λ/2m2gμν), we obtain the system of equations
(58a)
(58b)
(58c)
(58d)
The third equation is solved by GC(x, x′) = 0, the last equation tells us that GE = −G0(x, x′), and the second equation is solved by
(59)
with M B 2 = ζ ( 2 ζ 1 ) 1 ( ξ λ ) m 2 . Using the divergence identity (32) and the equation of motion (29) for the Feynman gauge vector propagator, we obtain M A 2 = ( ξ λ ) m 2 and the equation
(60)
which is solved by
(61)
So we obtain the ghost propagators
(62a)
(62b)
(62c)
(62d)
with the masses
(63)
We can also easily compute the massless limit. For this, let us define the mass derivatives
(64)
of the scalar propagator, which fulfill the equation of motion
(65)
We then have
(66)
and, keeping λ fixed, obtain the massless ghost propagators
(67a)
(67b)
(67c)
(67d)
Note that λ fixed means that both Λ and m2 tend to zero, with their ratio fixed, and we see that the massless limit of the ghost sector does not depend on this ratio. However, we can also keep Λ fixed and send m2 → 0. For this, we simply replace again λ = 4/(n − 2)Λm−2, whence the limit m2 → 0 results in
(68a)
(68b)
(68c)
(68d)
and we see that the propagator G μ ( x , x ) = i c ̄ μ ( x ) c ( x ) is actually divergent. This is related to the rescaling that we did at the beginning, which only has a sensible limit for vanishing cosmological constant.69 
On the other hand, if we first take the limit Λ → 0 and then m → 0, using the mass expansion (66) we obtain the finite result
(69a)
(69b)
(69c)
(69d)
We thus see that the limits Λ → 0 and m → 0 do not commute, which is an example of the van Dam–Veltman–Zakharov (vDVZ) discontinuity.49,50
Analogously to the ghost sector, we write the full gauge-fixed action (21) as
(70)
with the composite field HA = (hμν, Aμ, ϕ) and the Hermitean equation-of-motion operator PAB. Since the shifted auxiliary fields B ̃ μ and B ̃ decouple, we can treat them separately at the end. The propagator
(71)
satisfies
(72)
where a sum over the composite index B is implicit, or in components
(73a)
(73b)
(73c)
(73d)
(73e)
(73f)
(73g)
(73h)
(73i)
The equation-of-motion operators PAB can be read off from the full action (21), and are given by
(74a)
(74b)
(74c)
(74d)
(74e)
(74f)
As for the ghost fields, we make a general ansatz
(75a)
(75b)
(75c)
(75d)
(75e)
(75f)
(75g)
(75h)
(75i)
where the propagators with two indices are (linear combinations of) Feynman-gauge vector propagators of different masses, and the ones without indices (linear combinations of) scalar propagators of different masses. This system is solved as for the ghost sector, using the divergence and trace identities for the Feynman gauge propagators. We do not give any details of the lengthy computation (which we performed using the xAct suite of tensor algebra packages114,115,148), but only the result:
(76a)
(76b)
(76c)
(76d)
(76e)
(76f)
(76g)
(76h)
(76i)
where we recall the masses (63)
(77)
and defined the abbreviations
(78)
Also for the field propagators, we can take the massless limit, and again we have the choice of keeping λ fixed (and thus sending Λ → 0 together with m), or holding Λ fixed. In the first case, we obtain
(79a)
(79b)
(79c)
(79d)
(79e)
(79f)
(79g)
(79h)
(79i)
where we also defined the mass derivative of the Feynman-gauge vector propagator
(80)
which fulfills the equation of motion
(81)
and the second mass derivative of the scalar propagator
(82)
which fulfills
(83)
We see that while the vector Aμ decouples in this limit, such that the mixed propagators of Aμ and other fields vanish, the scalar ϕ is still coupled to the graviton hμν. This is again related to the vDVZ discontinuity:49,50,69 the scalar mode of the graviton does not decouple in the massless limit, and instead still interacts with other scalar fields.
On the other hand, we can also take the massless limit while holding Λ fixed. Defining Λ ̃ 4 / ( n 2 ) Λ , this results in
(84a)
(84b)
(84c)
(84d)
(84e)
(84f)
(84g)
(84h)
(84i)
As in the ghost sector, we have a divergence for small masses, which comes from the fact that the limit m → 0 is not continuous for Λ ≠ 0, i.e., ultimately from the fact that we have rescaled some fields. However, since the divergent terms are pure gauge, namely symmetrized derivatives, they will not contribute to correlators of invariant observables.
On the other hand, if we first take the limit Λ → 0 and then m → 0, we again obtain a finite result
(85a)
(85b)
(85c)
(85d)
(85e)
(85f)
(85g)
(85h)
(85i)
and also in this limit the scalar ϕ is still coupled to the graviton hμν, a manifestation of the vDZV discontinuity.
Finally, we consider the shifted auxiliary fields B ̃ μ and B ̃ defined by (22), which decoupled in the action (21). Since they have no kinetic term, their propagator simply reads
(86)
Using the propagators (76a), we can then determine the propagators involving the original auxiliary fields Bμ and B. First we determine propagators involving a single Bμ or B, which follow from the fact that propagators involving B ̃ μ or B ̃ and a dynamical field ϕ, Aμ or hμν vanish. This results in
(87a)
(87b)
(87c)
(87d)
(87e)
(87f)
These results were again computed with xAct, using the trace (42) and (43) and divergence identities (49) and (50) for the Feynman-type gauge tensor propagator, as well as the divergence identities (32) and (33) for the Feynman gauge vector propagator and its equation of motion (29). We now take the massless limit, holding either λ or Λ fixed. In the first case, we obtain
(88a)
(88b)
(88c)
(88d)
(88e)
(88f)
and we see that most propagators vanish, except the one between the diffeomorphism auxiliary field Bμ and the graviton hμν and the one between the Stueckelberg auxiliary field B and the Stueckelberg vector Aμ. That is, in the massless limit the auxiliary field sector decouples into separate tensor and vector sectors, which is of course very reasonable. In the second case, we obtain instead
(89a)
(89b)
(89c)
(89d)
(89e)
(89f)
where we set as before Λ ̃ 4 / ( n 2 ) Λ = λ m 2 . We see once more that this is not a sensible limit, since some propagators diverge. Taking thus first the limit Λ → 0 and afterwards the massless one, we obtain the same limit as in Eq. (88a).
Using the propagators (86) and (87a), we can then also determine the propagators involving two auxiliary fields Bμ or B. We obtain
(90a)
(90b)
(90c)
We see that the auxiliary fields do not propagate, such that the mixed propagators (87a) only become relevant when the auxiliary fields appear in external states, i.e., explicitly in the correlation function. This is relevant for the verification of Ward identities, which we do in Subsection III E.

Ward(–Takahashi–Slavnov–Taylor) identities116–119 are relations between correlation functions that stem from the gauge invariance of the underlying classical theory. In the BRST formalism, they arise from the fact that any BRST-exact term has vanishing expectation value in a physical state, i.e., s ( anything ) = 0 . In our case, they relate the different propagators of the fields, ghosts and auxiliary fields, and we will display explicitly a subset of all possible identities.

A non-trivial group of identities is obtained by considering the expectation value of the BRST transformation of a field and an antighost. Using the explicit BRST transformations (14a) and (15a), we obtain
(91a)
(91b)
(91c)
(91d)
(91e)
(91f)
which connects the ghost propagators (62) with the ones containing an auxiliary field (87a). Inserting the explicit expressions for the propagators, one checks easily that the identities (91a) are indeed fulfilled.
Also the non-propagation of the auxiliary fields, i.e., that fact that their propagators vanish (90a), is the consequence of a Ward identity. Namely, we consider the expectation value of the BRST transformation of an auxiliary field and an antighost. Using the explicit BRST transformations for antighosts and auxiliary fields (15a), it follows that
(92a)
(92b)
(92c)
(92d)
Various other identities can also be derived as consequences of a Ward identity, for example the divergence identities (32) and (33) for the Feynman gauge vector propagator and the trace (42) and (43) and divergence identities (49) and (50) for the Feynman-type gauge tensor propagator. For details of the derivation and further identities, we refer the reader to Ref. 43.

The heat kernel technique was introduced in quantum field theory as a way to treat functional traces and determinants of local differential operators of Laplace type, i.e., second-order elliptic differential operators whose principal symbol is the metric (also called minimal operators). Later on, it was generalized to also include non-minimal operators, where a more general principal symbol is admissible, and which arise for example in gauge theories in a gauge different from Feynman gauge, as well as hyperbolic operators. In the following we will give a short overview of the technique, and refer to the reviews and books1,2,5,120–122 for more in-depth results.

Formally, the heat kernel associated to a Laplace-type operator P is given by
(93)
i.e., it is the kernel of the heat semigroup whose generator is the operator P, which we take to be negative definite (or more generally, −P should be bounded from below). More concretely, we consider a Laplace-type operator PAB acting on sections of a vector bundle whose elements are denoted by vA, which is a differential operator of the form
(94)
with 1 A B the identity on the bundle under consideration, and where EAB contains at most first-order derivatives. For example, for the standard scalar Laplacian acting on (spin-0) functions, we have 1 = 1 and E = −ξR with the scalar curvature R and a parameter ξ R , for spin-1 vector fields we have 1 μ ν = g μ ν and for a spin-2 symmetric tensor we have 1 ( μ ν ) ( ρ σ ) = g μ ( ρ g σ ) ν ; in general, 1 A B is a properly symmetrized product of metrics for a spin-s tensor field. The heat kernel KAB(x, x′, τ) is then the solution to the heat equation
(95)
with boundary condition
(96)
which is unique assuming some mild conditions on the growth at infinity.123,124 Furthermore, K is bounded for large x or x′ for any τ > 0.123 
For n-dimensional flat space with P A B = 1 A B 2 , the heat kernel is known explicitly:120,122
(97)
In Fourier space, it is clear that for τ = 0 we obtain a δ distribution and the correct normalisation. However, for a general Riemannian manifold, no explicit closed-form expression can be given. Instead, it is possible to determine successive approximations to the heat kernel, and two main schemes are known: the first one is the asymptotic local or early-time expansion in terms of powers of τ (see Ref. 123, Theorem 2.30)
(98)
where σ is the Synge world function,125 equal to half of the square of the geodesic distance between x and x′, and Δ is the van Vleck–Morette determinant.126 The quantities A k A B are called heat kernel coefficients, and their coincidence limits A k A B ( x , x ) are local functions of the curvature tensors and their covariant derivatives. The second scheme is based on a non-local expansion of the heat kernel around the flat-space limit in powers of the curvature tensors,3,127,128 and resums certain powers of τ in the expansion. Since we will only use the first scheme in this work, we do not enter into further details.

For non-minimal operators, there is no direct heat kernel expansion. Unfortunately, many of the operators that are of physical interest are not of Laplace-type. A well-known example is Yang-Mills theory in a covariant gauge with parameter ξ different from Feynman gauge ξ = 1, where the bundle is composed from the tangent bundle and an internal Lie algebra bundle (depending on the gauge group). While for a semi-simple Lie algebra the internal part of PAB (the Cartan–Killing form) is proportional to the identity δab, the covariant spacetime derivatives are uncontracted, and P A B = δ a b g μ ν 2 ( 1 ξ ) μ ν is a non-minimal differential operator. One possibility to treat such operators is the use of covariant projectors which lead to Laplace-type differential operators acting on subspaces of the original field space, i.e., on a different vector bundle. Typical examples for such projections are the transverse decomposition of a vector field, or the transverse traceless decomposition of a symmetric tensor field, for which one can in particular use the second approximation scheme explained above.129–132 Alternatively, one can write the heat kernel for the non-minimal operator as the sum of the one for the minimal operator and covariant derivatives of the heat kernel of lower spin fields as in Ref. 2, Sec. 2. This is akin to the determination of the corresponding Green’s functions in terms of the Feynman-type gauge ones, and is based on the Ward identities that hold in the free theory analogously to the trace and divergence identities that we computed in Sec. III. In fact, we will derive below trace and divergence identities for the (minimal, Feynman-type gauge) heat kernels (and their coefficients), using which the corresponding non-minimal heat kernel could be computed. Unfortunately, the relation between the two approaches to non-minimal operators and their heat kernel expansions is not straightforward, mainly because the known expansions are in general only asymptotic. One can obtain a convergent expansion by multiplying the coefficients A k A B ( x , x ) in the local expansion (98) by cutoff functions whose support decreases quickly enough with k (see Ref. 133, Sec. 5.2). In this way, the same coincidence limits (as x′ → x) as for the unmodified expansion hold, but in any finite neighborhood the sum only has a finite number of terms; only for analytic spacetimes one can hope to have a convergent local expansion without cutoffs. Fortunately, for most purposes it is enough to know the coincidence limits of a finite number of coefficients and their derivatives, which are always well-defined.

Given the heat kernel, the unique solution of the differential equation
(99)
for the Green’s function of Pm2 is then obtained by integrating over the proper time τ:
(100)
Using integration by parts and the fact that K is bounded, one sees that this indeed fulfills the required differential equation:
(101)
For later convenience, we have separated the mass so that P is the massless differential operator. As we see, a positive mass ensures convergence for large τ.
We now assume that the differential operator has the form (94) where E does not involve derivatives, which is all that we need. Inserting the asymptotic expansion (98) into the Eq. (95) satisfied by the heat kernel and comparing powers of τ, we obtain the transport equations
(102)
with the initial condition
(103)
where we also used the following identities for the Synge world function:111 
(104)
Taking into account the boundary condition for the heat kernel, the initial condition is solved by
(105)
with the spin-s parallel propagator G . For integer spin, it is given by the symmetrized product of spin-1 parallel propagators gμρ:
(106)
where A = μ1μs, B = ρ 1 ρ s , and where gμρ satisfies111 
(107)
The solution for the higher-order coefficients A k A B can be given exactly in Riemann normal coordinates.134–136 In these coordinates, the geodesics from x′ to x are straight lines yμ(λ) = (x′)μ + λ[xμ − (x′)μ], such that σ ( x , x ) = 1 2 g μ ν ( x ) [ x μ ( x ) μ ] [ x ν ( x ) ν ] and ∇μσ(x, x′) = xμ − (x′)μ. For k > 0, we thus compute
(108)
and it follows that
(109)

The main application of the heat kernel in quantum field theory is in the renormalisation of the effective action or of composite operators, for which only the coincidence limit of a finite number of coefficients and their derivatives are needed. In particular, for the renormalisation of the stress tensor in n dimensions one needs the coincidence limit of the coefficients Ak with k = 0, …, n. To check that the renormalized stress tensor is covariantly conserved, one needs in addition the coincidence limit of their derivatives μ 1 μ A k with k + n + 1. To determine these coincidence limits, one can either use the explicit formula (109) in normal coordinates, or one can take covariant derivatives of the transport Eq. (102) and their coincidence limit. In turn, these can be obtained recursively using the coincidence limits of derivatives of the Synge world function σ and the van Vleck–Morette determinant Δ.

Several methods have been developed to find explicit expressions for the heat kernel coefficients.137 DeWitt28 determined the first two coefficients with a covariant recursive method. Sakai134 relied on the Riemann normal coordinates to find the third coefficient in the scalar case (i.e., for a single scalar field on a curved space). For the general case, this coefficient was found by Gilkey138 using a noncovariant pseudo-differential-operator technique (see also Ref. 139). The integrated and traced fourth and fifth coefficients for an arbitrary field theory in flat space were found in Ref. 140 through the evaluation of a noncovariant Feynman graph. Avramidi141,142 presented a new covariant non-recursive procedure and found the fourth coefficient for the general case; the coefficient a5 in flat space was computed in Refs. 140 and 143. Various resummation methods have also been developed, see Ref. 144 and references therein. Finally, the connection between the Riemannian heat kernel and the Lorentzian analogue, which we treat in the following, was elucidated for stationary spacetimes by Strohmaier and Zelditch.145 

In the Lorentzian case, the heat kernel expansion is only formal even for positive mass, since ∇2 is now an unbounded operator whose spectrum (in general) is the full real line. We thus consider instead the Wick-rotated version of the heat kernel, which fulfills121,146,147
(110)
with the boundary condition
(111)
In this case the heat kernel satisfies a Schrödinger-like equation with PAB playing the role of the Hamiltonian, since the problem has been reformulated into an equivalent “quantum mechanical” problem, solving the “evolution” Eq. (110) in the Schwinger proper time.122 We assume that there exists a unique solution of this equation which does not grow faster than exponentially. The Feynman propagator corresponding to Pm2 can then be defined as
(112)
where the limit ϵ → 0 is to be taken in a distributional sense, i.e., after integrating over τ and smearing with test functions. Since we assume that K does not grow too fast, the τ integral is convergent for ϵ > 0, and analogously to the Riemannian case we can easily check that the Feynman propagator fulfills the right equation:
(113)
For n-dimensional flat space with P A B = 1 A B 2 , we again have an explicit expression for the heat kernel120,122
(114)
In Fourier space, it is again clear that for τ = 0 we obtain a δ distribution and the correct normalisation. For a general pseudo-Riemannian manifold, we make an ansatz analogous to the Riemannian case (98) for the small-τ asymptotic expansion:1,121,122
(115)
and assume the same form (94) for the differential operator. Let us remark that this expansion, as the one in the Riemannian case, is only valid for minimal operators, such as a vector field in Feynman gauge. For non-minimal operators, such as the Proca operator for a transverse vector field, one has to use either covariant projectors or add the appropriate multiple of derivatives of the heat kernel of lower spin fields as explained above.
Inserting the expansion (115) into the Eq. (110) satisfied by the heat kernel and comparing powers of τ, we obtain the transport equations
(116)
and the initial condition
(117)
We see that we obtain the same transport equations as in the Riemannian case, which is due to our inclusion of an explicit factor of ik in the expansion (115).

For xx′, the factor e i σ ( x , x ) 2 τ in the heat kernel expansion (115), or e σ ( x , x ) 2 τ in the expansion (98), ensures that the proper time integrals (100) and (112) converge at τ = 0. However, in the coincidence limit x′ = x which is needed for applications, this factor becomes 1 and these integrals are divergent for small τ. This is of course the well-known divergence of the propagators in the coincidence limit.148 In particular, it is clear from the expansions (98) and (115) that the coefficients A0, …, An/2−1 have divergent prefactors if one computes the coincidence limit of the propagator, and more if one is interested in the coincidence limit of derivatives acting on it.

To regulate this divergence and extract the finite part of the result, analytic regularisation schemes are well suited for the heat kernel expansion. A well-known such scheme is dimensional regularisation, where one continues the dimension n of spacetime into the complex plane, computes results in a region where all integrals are convergent, and then analytically continues the result back to the physical dimension. In particular, for R e n small enough (depending on the number of covariant derivatives acting on the heat kernel), the proper time integrals are convergent for small τ even in the coincidence limit. The original divergence of the propagator then manifests itself in poles as n approaches the physical dimension. While dimensional regularisation is versatile and usually the simplest one, it also has problems, namely when the objects under study cannot naturally be defined in arbitrary dimension n. In particular, this concerns chiral, conformal and supersymmetric theories, where the symmetries are dimension-dependent. These symmetries can then be broken in the quantum theory and one obtains chiral, conformal and supersymmetry anomalies.149,150 Of course anomalies are physical and cannot depend on the choice of a regularisation scheme, but only on renormalisation conditions (as can be seen in the dichotomy between covariant and consistent anomalies151). Therefore, for dimension-dependent symmetries where the advantages of dimensional regularisation are absent, it can be easier to employ a different regularisation scheme. For the heat kernel, such a scheme is obtained by inserting a factor τδ in the integrals (98) and (115) and performing the analytic continuation in δ.

We would like to compute the coincidence limit of up to three covariant derivatives of the propagator in n = 4 dimensions. For this, we need the coefficients Ak with k ≤ 4 and certain derivatives of them, whose coincidence limit we compute in arbitrary dimensions. As ingredients, we need the coincidence limit of derivatives of the Synge world function σ and the van Vleck–Morette determinant Δ, which enter the transport Eq. (102) (we recall that those are the same in the Riemannian and Lorentzian case). These can be computed recursively using Eq. (104), taking the required number of covariant derivatives and then the coincidence limit. The results are dimension-independent and have been determined many times, see for example;126,152,153 for completeness, we here reproduce the required ones:
(118a)
(118b)
(118c)
(118d)
(118e)
(118f)
(118g)
(118h)
Non-symmetrized derivatives can be obtained from this by expanding the symmetrized expressions and commuting the covariant derivatives into the required order. Such computations are best done with the help of a computer algebra package such as xAct,115–117,154 see also Ref. 155 for a specialist package. In particular, we have
(119a)
(119b)
where we used the second Bianchi identity to simplify the results.
The coefficients Ak can then be recursively computed using the transport Eq. (102) with the initial condition (105), taking into account the explicit form of the differential operator (94). We first rewrite the transport Eq. (102) as
(120)
and using the results (118a), take the coincidence limit of this equation, which gives
(121)
Taking a covariant derivative of Eq. (120) and then the coincidence limit, we obtain
(122)
where we also used the contracted second Bianchi identity ν R μ ν = 1 2 μ R . The pattern continues: to determine the coincidence limit of covariant derivatives acting on the coefficient Ak, we need to know the coincidence limit of up to + 2 covariant derivatives acting on the coefficient Ak−1. To determine the coincidence limit of A2 and its first derivative, which are needed to check conservation of the stress tensor in four dimensions, we thus need to know the coincidence limit of up to three derivatives acting on A1 and up to five derivatives acting on A0.
Generalizing the coincidence limit (107) for the spin-1 parallel propagator, the spin-s parallel propagator satisfies111,156
(123)
To obtain the coincidence limit of unsymmetrized derivatives, we have to commute them into the right order, using that
(124)
where RμνAB is the curvature tensor of the inner bundle described by the composite indices A. In particular, for scalars (spin 0) we have Rμν⋅⋅ = 0, while for vectors (spin 1) Rμνρσ is the usual Riemann curvature tensor, and for symmetric tensors (spin 2) with A = ρσ and B = αβ we have
(125)
A straightforward computation157,158 shows that the second Bianchi identity is fulfilled:
(126)
From the above equations, we thus obtain the coincidence limits of A1 and its first derivative:
(127a)
(127b)
To obtain the same for A2, we need in addition the coincidence limit of two and three derivatives acting on A1. Acting with two covariant derivatives on the transport Eq. (120) for k = 1, symmetrising them when they act on σ or Δ, using the explicit result (105) for A0 and taking the coincidence limit we obtain
(128)
To derive this result, we used the xAct tensor package151–154 and the FieldsX extension,159 the contracted second Bianchi identity and the identity
(129)
which follows from the first Bianchi identity. The same procedure for three derivatives gives
(130)
We now can insert these results into the Eqs. (121) and (122) for k = 2, which together with the previous results (127a) gives us
(131)
and
(132)
where we defined
(133)

In the following sections, we specialize these results to scalar, vector and tensor coefficients.

For scalars, the inner bundle is one-dimensional and trivial, and consequently its curvature tensor vanishes, Rμν⋅⋅ = 0. However, we consider a non-minimal coupling to curvature such that the differential operator (94) reads ∇2ξR, and thus E = −ξR. Our results (131) and (132) then specialize to
(134a)
(134b)
and we recall that this result is valid in any dimension.
In four dimensions, we can express the Riemann tensor in terms of the Weyl tensor to write A2 in terms of the square of the Weyl tensor
(135)
and the Euler density
(136)
This results in
(137)
and agrees with known results in the literature (see for example Ref. 128).
For vectors, the inner bundle is the tangent bundle of the manifold, and thus the bundle curvature tensor is the usual Riemann curvature tensor. We also consider a non-minimal coupling to curvature of the form Eμν = −ξRμνζgμνR; the vector propagator in Feynman gauge (29) has ξ = 1 and ζ = 0. Our results (131) and (132) then specialize to
(138)
and
(139)
where we used the contracted second Bianchi identity to simplify the result. We note that the coincidence limit of A2μν(x, x′) is symmetric in μ and ν, which follows from the fact that the heat kernel coefficients themselves are symmetric:147,160,161 A2μν(x, x′) = A2νμ(x′, x). However, the difference lim x x ρ A 2 μ ν ( x , x ) 1 2 ρ A 2 μ ν ( x , x ) is antisymmetric in μ and ν, which follows from this symmetry together with the Synge rule111 
(140)
where the derivative ∇μ acts on x′. This provides a useful check on the computation.
In four dimensions we can again simplify these expressions somewhat by using dimension-dependent identities.162,163 These follow from the fact that antisymmetrising over five or more indices gives a vanishing result in four dimensions, and in general antisymmetrising over more than n indices in n dimensions. In particular, the Weyl tensor satisfies
(141)
such that
(142)
in terms of the Weyl tensor and the Euler density (136). We have also separated the full trace, such that the first two lines are traceless.

The coefficients again agree with results in the literature (see for example Ref. 128).

For symmetric tensors with A = ρσ and B = αβ, the bundle curvature tensor reads
(143)
and we also have 1 A B = δ α ( ρ δ β σ ) . Again we consider a non-minimal coupling to curvature of the general form
(144)
which is the most general form compatible with the symmetries Eμνρσ = E(μν)(ρσ) = E(ρσ)(μν) that make the differential operator hermitean, and we recall that the trace-reversed and rescaled tensor propagator in Feynman gauge (37) has ξ = 1 and ζ = χ = η = θ = 0. Our results (131) and (132) then specialize to
(145)
with
(146)
and
(147)
where we used the contracted second Bianchi identity to simplify the result.
In four dimensions, we can again use dimension-dependent identities to simplify the result, in particular (141). This gives
(148)
with the constant
(149)
where we also separated the (full) trace. For ξ = 1, the trace-reversed and rescaled tensor propagator in Feynman gauge, this further simplifies to
(150)
Again this agrees with known results in the literature (see for example Ref. 128). Recently, heat kernel coefficients for integer spin were computed for anti-de Sitter spacetime164,165 and the sphere166 (as the Euclidean analogue of de Sitter space), which also agree with our results if we specialize them to the appropriate background spacetime.

From the trace and divergence identities for the vector and tensor propagators in Feynman-like gauges, we obtain relations between their corresponding heat kernel coefficients.

Inserting the heat kernel expansion (115) in the divergence identity (32) for the vector propagator, we obtain
(151)
We now employ the small time regularisation, inserting a factor τδ and then computing the integrals over τ using the integral representation
(152)
of the modified Bessel function given in Ref. 167, Eq. (10.32.10). This results in
(153)
with Σ 2 M 2 ( σ + i ϵ ) , and we see that the small-time asymptotic expansion of the heat kernel is a large-mass expansion if we hold Σ fixed. It follows that each power of τ in the asymptotic expansion (151) must vanish separately, such that we have
(154)
and
(155)
for k ≥ 0.
Repeating this derivation with the identity (33), we also obtain
(156)
and
(157)
for k ≥ 0. Moreover, taking covariant derivatives of these relations and the coincidence limit, we obtain relations between the coincidence limits of the coefficients. In particular, without covariant derivatives we have
(158a)
(158b)
and using one covariant derivative we obtain
(159a)
(159b)
(159c)
In a completely analogous way, we derive from the trace and divergence identities (87) for the trace-reversed Feynman gauge tensor propagator relations for the tensor heat kernel coefficients. Inserting the heat kernel expansion (115), we obtain
(160)
and
(161)
where we have reexpressed λ in terms of Λ = n 2 4 m 2 λ . Since the relations between coefficients are obtained from a large-mass expansion, it is important that we keep the cosmological constant Λ and the mass m separate to derive them, since otherwise we would obtain an inconsistent expansion and wrong relations. Afterwards, we can of course express Λ again by the rescaled λ. As in the vector case, we then obtain the relations for the coefficients by comparing powers of τ, for which we now have to expand the exponential in Λ as well. For the trace (160), that means that the sought relations can be obtained from
(162)
where we used the general Leibniz rule
(163)
We thus have the relations
(164a)
(164b)
following from the trace identity, and we see that the trace of each tensor coefficient is a polynomial in the cosmological constant. In the same way, from (161) (and the analogous equation with x and x′ exchanged) we obtain
(165a)
(165b)
(165c)
(165d)
As in the vector case, we obtain relations for coincidence limits of the coefficients by taking covariant derivatives and then the coincidence limit of these relations. In particular, without covariant derivatives we have
(166a)
(166b)
and using one covariant derivative we obtain
(167a)
(167b)
(167c)
where we also used the previous relations (155) and (157) for the vector coefficients, and took into account that the background is on-shell, R μ ν = λ 2 m 2 g μ ν .

Having obtained all the ingredients, the heat kernel expansion can be used to derive the regularized propagators of a massive theory of gravity. In Sec. III we obtained the propagators in the various fields sectors, which we now express using the heat kernel.

The propagators of the ghost sector are given in Eq. (62). Expressing them in terms of the heat kernel (100) and inserting the heat kernel expansion (115), we obtain for the time-ordered propagators the expansions
(168)
(169)
and
(170)
where we recall the masses (77)
(171)
It is now straightforward to take the massless limit. Since we have seen that the limit m2 → 0 with Λ fixed is divergent, while we obtain a sensible limit keeping λ = 4 n 2 Λ m 2 fixed, we only give the expressions for the latter one. Now G F ( x , x ) (62d) is already massless, and the limit m → 0 of G μ F ( x , x ) vanishes (67b), while for the diffeomorphism ghost propagator we obtain
(172)
where we used the limit
(173)
Compared to (168), we have shifted the summation index in the last term to absorb the extra factor of τ.
The propagators of the field sector are given in Eq. (76a). As for the ghost sector, we express them in terms of the heat kernel (100) and insert the heat kernel expansion (115). However, there are two differences: first, that the propagators in the field sector also contain mass derivatives which are defined in Eqs. (64) and (80). For them, we thus obtain instead of Eq. (112) the heat kernel representation
(174)
which is less singular for small τ, consistent with the fact that the most singular term of the propagator is mass independent. Second, the tensor propagator in Feynman gauge G FG , m 2 does not have a straightforward heat kernel expansion, but the trace-reversed and rescaled tensor propagator G TG , m 2 (38) has. For the Feynman gauge tensor propagator, we obtain from this
(175)
Taking these differences into account, we obtain for the time-ordered propagators the expansions