Towards the renormalization group flow of Horava gravity in $(3+1)$ dimensions

We compute the renormalization group running of the Newton constant and the parameter $\lambda$ in $(3+1)$-dimensional projectable Horava gravity. We use the background field method expanding around configurations with flat spatial metric, but non-vanishing shift. This allows us to reduce the number of interaction vertices and thereby drastically simplify the calculations. The gauge invariant $\beta$-function of $\lambda$ has two families of zeros, attractive in the infrared and ultraviolet respectively. They are candidates for the fixed points of the full renormalization group flow of the theory, once the $\beta$-functions for the rest of the couplings are added.


I. INTRODUCTION
Since the formulation of general relativity (GR) a century ago and of quantum mechanics a few years later the quest for a theory unifying the two -quantum gravity (QG) -has been one of the biggest endeavors in theoretical physics. While other known fundamental interactions are successfully described within the formalism of local perturbative quantum field theory (QFT), most approaches to QG involve departures from this framework. Despite a number of impressive achievement of these research programs, the quest for QG remains largely unaccomplished. Thus, it is reasonable to wonder if there is still a way to formulate QG within the realm of QFT.
However, a straightforward approach to this goal fails due to the dimensionful character of Newton's constant G in dimensions greater than two, in particular in (3 + 1) dimensions. This renders the theory non-renormalizable, with the number of divergences increasing at every loop order and requiring an infinite set of counterterms suppressed by higher powers of the Planck mass M P ∼ √ G −1 ∼ 10 19 GeV. Thus, perturbative GR can be thought at most as an effective low-energy theory of gravitational interactions, see e.g. [1][2][3], unable to describe physics above the ultraviolet (UV) cutoff given by M P .
A simple way out of this problem in four space-time dimensions is to extend the gravitational Lagrangian with terms quadratic in curvature [4][5][6] so that UV divergences are softened by higher powers of four-momenta in the propagators. Although this theory is thereby renormalizable and even asymptotically free for some regions of the parameter space, unitarity is jeopardized by the presence of four time derivatives in the action, which gives rise to a propagating Ostrogradsky ghost with negative norm in the Hilbert space [7,8].
It was pointed out by P. Hořava [9] that unitarity can be preserved in this kind of models if we sacrifice Lorentz invariance (LI) at very high energies. This allows one to introduce into the Lagrangian extra higher spatial derivatives which will soften UV divergences while at the same time keeping only two time derivatives, so that the ghosts are avoided. For this proposal to work, one must supplement the spacetime manifold with a preferred foliation into spatial slices and a privileged time direction. Then it is possible to formulate a theory containing only marginal operators and relevant deformations with respect to an anisotropic (Lifshitz) scaling in (d + 1) dimensions 1 , restricted symmetry in terms of foliation preserving diffeomorphisms (FDiff), where t (t) must be a monotonic function. That is, we are left with invariance under time reparameterizations and spatial time-dependent diffeomorphisms. Hořava's proposal gave rise to numerous follow-up publications studying different aspects of Hořava gravity (HG) (see [10][11][12][13] for reviews). These include an exhaustive study of the low energy regime of the theory, which led to the identification of a version of HG -the healthy non-projectable model [14] -which can reproduce the known phenomenology of GR at the scales in which the latter has been tested [15]. While its parameter space has been strongly constrained by the observation of a gravitational wave signal from a binary neutron star merger [16], it still remains phenomenologically viable [17]. In this theory, a certain amount of LI violation persists at all scales [18]. This can lead to interesting cosmological applications [19,20]. However, it also demands the presence of a mechanism 2 suppressing the percolation of LI violation to the visible matter sector, where LI has been tested to a very high precision [29].
Regarding renormalizability, the situation is less clear. Although the Lagrangian of the non-projectable model is power-counting renormalizable, the lapse element of the metric induces an instantaneous interaction [18,30] which can potentially lead to non-local loop divergences [31]. Whether this indeed happens or not is an open question.
There has been a significant progress in understanding the quantum dynamics of a reduced version of HG known as the projectable model, where the lapse N is restricted to be only time-dependent. Then it can be algebraically fixed to take a given value, say N = 1. In this way one gets rid of time reparameterization invariance as a symmetry of the theory, being left only with spatial diffeomorphisms. This version of the theory does not have a stable homogeneous vacuum with flat Minkowski metric and thus is unlikely to reproduce realistic phenomenology, at least in the controllable regime of weak coupling [18,32,33] (see, however, [10,34,35]). Nevertheless, it captures several salient features of QG, such as propagating spin-2 gravitons (in spacetime dimensions greater than 3) and a large gauge group of coordinate transformations precluding existence of local gauge-invariant observables. The projectable model has been proven to be renormalizable in the strict sense, i.e. all divergences are proportional to gauge invariant operators already present in the bare action, in any space-time dimension [31,36].
Furthermore, the computation of the full renormalization group (RG) flow of the theory in (2 + 1) spacetime dimensions reveals the presence of an asymptotically free UV fixed point [37]. Unlike (2 + 1)-dimensional GR, this theory possesses a propagating scalar degree of freedom and thus provides an example of a perturbatively UV complete gravitational field theory with non-trivial local dynamics. Other results along these lines include the study of a conformal reduction of the (2+1)-dimensional theory [38], the computation of the anomalous dimension of the cosmological constant [39], also in (2 + 1) dimensions, and the contribution of scalar matter into the renormalization of the gravitational couplings [40].
In this paper we make a step towards uncovering the RG flow of HG in (3+1) dimensions. This case differs from the (2 + 1)-dimensional model in two respects. First, the physical spectrum of the theory now includes a spin-2 mode -graviton. Second, on the technical side, the structure of the propagators and vertices is significantly more complicated than in (2 + 1)-dimensional HG, which makes a direct application of the approach used in [37] infeasible. We find a way to modify this approach that allows us to compute the 1-loop β-functions of the two couplings entering the kinetic part of the action. Renormalization of the potential part requires different techniques and will be reported elsewhere.
The paper is organized as follows. In Sec. II we review the structure of projectable HG and its relevant properties. In Sec. III we describe our calculation. In Sec. IV we report the results for the β-functions and discuss their implications. Section V is devoted to conclusions. Appendices provide additional details.

II. PROJECTABLE HOŘAVA GRAVITY
Hořava gravity [9] is formulated by considering the spacetime metric together with a preferred foliation along the time direction. This induces an Arnowitt-Deser-Misner (ADM) decomposition of the metric, The lapse N , the shift N i and the spatial metric γ ij on the foliation slices transform in the standard way under FDiff transformations, Their dimensions under the anisotropic scaling (1) are, 3 The FDiff symmetry, together with time-reversal invariance, parity and the requirement of power-counting renormalizability fix the action of HG to be Here G and λ are dimensionless coupling constants and the extrinsic curvature of the foliation is given by where dot denotes time derivative and ∇ i is the covariant derivative associated to the spatial metric; K ≡ K ij γ ij denotes the trace of K ij . FDiff transformations (4) are compatible with an assumption of the lapse function being constant on the foliation slices, N = N (t). This specifies the projectable version of the model [9], to which we restrict in what follows.
In this case, we can set N = 1 as an algebraic gauge fixing for time reparameterization invariance, eliminating it from the set of dynamical fields and leaving us only with spatial diffeomorphisms as a gauge symmetry. The potential V contains all possible marginal and relevant operators 4 invariant under FDiff and without time derivatives. In d = 3 and after using Bianchi identities and integration by parts, these reduce to [41], where R ij is the Ricci tensor corresponding to the spatial metric and we have used the fact that in three dimensions the Riemann tensor is expressed in terms of R ij . The spectrum of perturbations propagated by this action contains a transverse-traceless (tt) graviton and an additional scalar mode [31,32]. Both modes have positive-definite kinetic terms as long as G > 0 and λ ∈ −∞, 1 3 ∪ (1, ∞) and the theory admits unitary quantization. Their dispersion relations around a flat background 5 are 3 We say that a field Φ has dimension [Φ] if it transforms under (1) as Φ → b [Φ] Φ. Throughout the paper we will understand the notion of dimension in this sense. 4 I.e. operators with dimension less or equal than the absolute value of the dimension of the spacetime integration measure dtd d x, which is equal to 2d. 5 By this we mean the configuration N i = 0, γ ij = δ ij . It is a solution of equations following from (6), (8) provided the cosmological constant Λ is set to zero.
where k is the spatial momentum and we have defined A problem arises in the low-energy limit, where the dispersions relations are dominated by the terms proportional to k 2 . Due to the negative sign in front of this term, the scalar mode behaves as a tachyon at low energies, implying that flat space is not a stable vacuum of the theory. Attempts to suppress the instability by choosing λ close to 1 lead to the loss of perturbative control [18,32,33] (see, however, [10,34,35] for suggestions to restore the control by rearranging the perturbation theory at the classical level). Alternatively, the instability can be eliminated by tuning η = 0 or by expanding around a curved vacuum. In both cases the theory does not appear to reproduce GR in the low-energy limit, as there is no regime where the dispersion relation of the tt-graviton would have the relativistic form ω 2 tt ∝ k 2 . On the other hand, both dispersion relations (9) are perfectly regular at high momenta for ν 5 , ν s > 0. This is the relevant region for the structure of UV divergences which are the focus of the present paper.
For the rest of this work we restrict to the sector of marginal operators that dominate the UV behavior of the theory by setting η = µ 1 = µ 2 = 0. We also perform a Wick rotation to "Euclidean" time τ , which amounts to flipping the sign of the potential term. Thus we consider the action, We are interested in the 1-loop correction to this action and the corresponding RG flow of the couplings {G, λ, ν 1 , . . . , ν 5 }.

III. THE BACKGROUND-FIELD APPROACH
In our calculation we use the background field method [42]. We expand the fields about arbitrary background values, and expand the action (11) up to second order in the fluctuations h ij , n i . One-loop divergences will be captured as contributions to operators carrying background fields after integrating out the fluctuations. This procedure should work in general and, similarly to the case in (2 + 1) dimensions, we should be able to obtain the contributions to the whole set of coupling constants {G, λ, ν a } by an appropriate choice of background. However, the complexity of computations quickly grows with the number of space-time dimensions. This stems mainly from the increased number of terms in the potential part of the action and their more complicated structure. Indeed, in (2 + 1) dimensions there is a single term R 2 which generates about twenty different contributions when expanded about an arbitrary background (12). In (3 + 1) dimensions we have five operators shown in (11) which generate hundreds of structures. This makes a complete calculation of the RG flow in any dimension higher than (2 + 1) highly challenging. In particular, it does not appear feasible to apply the diagrammatic approach of [37] to renormalization of the operators in the potential. Rather, one will have to use the covariant techniques developed in [43][44][45] or other advanced tools. Fortunately, the situation is different for the operators in the kinetic part of the action, i.e. the two operators involving the extrinsic curvature. We observe that the latter contains the shift vector N i in a specific combination. Then, exploiting the known gauge invariance of divergences [31,36], we can extract the renormalization of G and λ not only from the terms with time derivatives of the spatial metric as in [37], but also from terms with spatial gradients of N i . Since the shift only enters into the kinetic term, all interaction vertices contributing to the renormalization of these operators will also come exclusively from the expansion of the kinetic term. The potential, as we will see shortly, only contributes into the propagators. Furthermore, we can restrict to backgrounds with flat spatial metricγ ij = δ ij , as long as we allow the background shiftN i to take arbitrary space-dependent values. All this together means that we can derive the renormalization of G and λ avoiding the complexity mentioned above and allows us to carry out the computation in a reasonably simple way.
According to the above discussion, we consider a family of backgrounds characterized by flat spatial metric and an arbitrary space-dependent shift. Thus, we specify (12) to The Lagrangian is then expanded to quadratic order in h ij and n i . Note that we define the perturbation of the shift with lower indices: this choice turns out to be more convenient as it reduces the number of terms in the perturbed Lagrangian. The perturbation of the component with the upper indices is then given by where on the r.h.s. the indices are raised using the flat background metric δ ij . From now on we will not distinguish upper and lower indices in the expressions expanded around the background (13); all repeated spatial indices will be simply summed over.
We are interested in those operators in the effective action for the background fields that contain the background shift. With our choice of flat background spatial metric, these are By computing the one-loop contribution to these operators we will be able to capture the renormalization of G and λ.

A. Covariant gauge fixing
The action (11) is invariant under the residual part of FDiff that preserves the condition N = 1, namely time-dependent spatial diffeomorphism of the form In order to proceed with the quantization of the theory we must thus append the bare action with a gauge fixing term for this symmetry. Within the background-field approach, the latter is chosen to be invariant under the background gauge transformations that in our case read, This ensures that the one-loop background effective action consists of gauge invariant operators and allows one to extract the renormalized couplings as the coefficients in front of such operators [42]. By comparing (17) to (4) it is straightforward to see that the shift perturbation with upper index δN i ≡ N i −N i and the metric perturbation h ij transform as a contravariant vector and a covariant tensor respectively 6 . They serve as the building blocks to construct the gauge fixing term. We choose the gauge fixing introduced in [31] which ensures that the propagators of all fluctuations decay uniformly in frequency and momentum in a way compatible with the anisotropic scaling. It reads, where Here h ≡ h i i is the trace of the metric fluctuations, ∆ ≡ ∂ i ∂ i is the flat-space Laplacian, σ and ξ are free parameters and is the covariant time-derivative of δN i . Note that we have added to the gauge-fixing function a contribution proportional to the background extrinsic curvatureK i j with an arbitrary coefficient ς. This additional freedom, which was not present in the original gauge fixing of [31], can be used to simplify the Lagrangian for perturbations and as an extra check of gauge invariance of our results.
The presence of the non-local operator O ij leads to appearance of non-local terms in the gauge-fixing action that are quadratic in δN i . Following [31,37] we eliminate this nonlocality by integrating in an extra field π i , so that we rewrite, Then the whole gauge-fixing action takes the local form, Finally, we have to add the anticommuting Faddeev-Popov ghost c i and antighostc i . The action for them is constructed out of the gauge fixing function as where the BRST operator s can be understood as implementing a gauge transformation with the gauge parameter replaced by the ghost, To sum up, the total action that we are going to work with reads, where S (2) is the result of expanding (11) to quadratic order in h ij , n i and S gf , S gh are given by (22), (23) respectively.

B. Propagators
To set up the computation along the lines of [37], we need to derive the propagators for the dynamical field {h ij , n i , π j , c i ,c j }. They are obtained by setting the background shift N i to zero in the total quadratic action for the fluctuations (25). The resulting Lagrangian splits into several pieces, In deriving these expressions we used integration by parts. Note that our choice of the gauge fixing ensures cancellation of the terms mixing h ij and n i in the expansion of the bare Lagrangian [31]. Inversion of the differential operators appearing in (26) yields the propagators, wherek i ≡ k i /k is the unit vector along the momentum k i , and the pole structures are given by, These expressions coincide with those derived in [31]. Note that the first two poles correspond to the dispersion relation of the transverse-traceless and scalar gravitons of the physical spectrum of the theory, while the rest are gauge dependent and their positions can be tuned by changing the gauge parameters σ and ξ.

C. Vertices
The second ingredient needed for the computation are vertices describing the interaction of externalN i -legs with the dynamical fields. The relevant vertices are depicted in Fig. 1 and are obtained by expanding the Lagrangian up to the given order in the number of fields. The corresponding formulas are rather lengthy and are relegated to the Appendix A.
The most important advantage of focusing on the renormalization of shift-dependent operators is that all vertices come exclusively from the kinetic terms and are completely independent of the structure of the potential, as well as of the gauge parameters σ and ξ. Interestingly, this implies that the vertices are universal for HG in any spacetime dimensions. Another implication is that they are independent of the couplings ν 1 , ν 2 , ν 3 . As the propagators (27) also do not contain these couplings, we conclude that the one-loop renormalization of G and λ is insensitive to them altogether.

D. One-loop diagrams
We are interested in capturing the renormalization of the operators ∂ iNj ∂ iNj and (∂ iNi ) 2 , which in momentum space read p i p iNjNj , p i p jNiNj , where p i denotes the momentum flowing through theN i -legs. The diagrams contributing to these operators are shown in Fig. 2. The prototypical form of the integrals appearing in the fish diagrams (upper row and first two diagrams in the lower row in the figure) is Here ω and k are the loop frequency and momentum, p is the external momentum, Π 1,2 are the propagators of the fields in the loop, and Q is a polynomial in momenta and frequency coming from the vertices. The propagators and vertices carry tensor indices that we have suppressed. We have also used that the external frequency vanishes as the background shiftN i was chosen to be time-independent. Similarly, the tadpole diagrams (the last two diagrams in Fig. 2) contain integrals of the form, The above integrals feature logarithmic divergences that we want to extract. To this end, we Taylor expand the integrands in (29), (30) and focus on the terms quadratic in the external momentum p i . Next we average over directions of the loop momentum using the standard formulas specified to d = 3 dimensions [46], In this way we are left with a sum of terms containing various tensor structures multiplied by the integrals with constant exponents a, b, A, B. Here P I,J , I, J ∈ {tt, s, 1, 2} are the pole structures listed in (28). Note that at most two different pole structures can appear in an individual integral. The exponents in (33) satisfy the relation so that the integrals are logarithmically divergent.
To regularize the UV divergences we use the Schwinger representation for each pole function, where Γ(z) is the gamma-function. The integration over ω and k is then performed using the formula, Finally, if the resulting integral involves two Schwinger parameters, we integrate over one of them. Namely, we write, where we made the change of variables z = s 2 /s 1 and used the relation (34). The integral over z is convergent and for all values of a, b, B appearing in our calculation can be taken in terms of elementary functions. The integral over the last remaining Schwinger parameter is kept to the end of the calculation and encodes the logarithmic divergence.
We have implemented the computational procedure described in this section in a Mathematica code. Manipulation of tensor objects was performed using the package xAct 7 . Some parts of the computation were repeated using an alternative code written in FORM 8 [47]. For a check of the procedure and the code we also performed the calculations in the case of (2 + 1) dimensions where we reproduced the results of [37].

IV. THE β-FUNCTIONS OF G AND λ
Evaluation of the one-loop diagrams depicted in Fig. 2 produces a divergent correction to the effective Lagrangian forN i , Comparing this with (15) we extract the renormalized couplings, We interpret the logarithmic divergence as the result of integrating out high-momentum modes in the Wilsonian approach to renormalization. Thus, we identify 9 , where Λ UV is a UV cutoff for momentum and k * is the RG reference scale. β-functions are defined as sensitivities of the couplings to k * , which yields, The computation of the coefficients C 1,2 for general gauge parameter {σ, ξ, ς} is rather cumbersome due to the presence of four different pole structures in the propagators that generate a large number of terms for every diagram. The situation is simplified if one chooses the gauge parameters in such a way that some of the pole structures coincide. In this section we fix which collapses the two poles corresponding to the gauge modes, P 1 = P 2 . In order to keep track of the gauge (in)dependence of the result we allow the two remaining gauge parameters σ, ς to take arbitrary values. Results for two other choices of gauge are presented in Appendix B.
We observe that both expressions are independent of the gauge parameter ς. On the other hand, β G depends on σ. This is not unexpected. Indeed, a change of gauge is known to affect the background effective action by adding to it a contribution which is a linear combination of equations of motion [48,49]. Such terms vanish for on-shell background configurations, but not in general. The background used in our calculations is not on-shell, so the effective action evaluated on it is gauge-dependent. In Appendix C we show that a change of gauge induces the following shift of the renormalized couplings, where is an infinitesimal parameter. We see that G, as well as ν a , transform non-trivially and therefore their renormalization and β-functions are not gauge-invariant. On the other hand, λ remains untouched by (46). It represents an essential coupling of the theory and its β-function must be gauge-independent. This is indeed born out by our computation, see (44). Other essential couplings can be chosen, for example, as G/ √ ν 5 , α and ν a /ν 5 , a = 1, 2, 3.
Let us discuss the properties of the gauge-invariant function β λ , eq. (44). Of particular interest is the manifold in the parameter space of the theory where this function vanishes. It constitutes a potential location of the fixed points of the full RG flow. The contribution of G/ √ ν 5 factors out, so that the sign of β λ and the location of its zeros are determined only by two parameters λ and α. We have the following regions in the plane of these parameters 10 , see Fig. 3: In this region β λ is always positive. So it does not contain any fixed points.
• λ > 1 , α < 9/4 β λ is positive at λ < λ α and negative at λ > λ α , where λ α is again given by (47). Now the line (λ α , α) is attractive along the λ-direction in UV. When embedded into the full 6-dimensional space of the essential couplings of the theory, this line is promoted to a codimension-one surface. UV fixed points of the full RG flow, if any, should lie on this surface.
Presence of UV attractive zeros of β λ at finite values of λ > 1 is analogous to the situation in (2 + 1) dimensions, where such zero is known to correspond to an asymptotically free UV fixed point [37].
To carry the comparison with the (2+1)-dimensional case further, it is instructive to work out the behavior of β λ when λ approaches the boundaries of the allowed regions, λ → 1/3 − or λ → 1 + , with all other couplings in the Lagrangian held fixed. As implied by (9b), these limits correspond to degenerate dispersion relation of the scalar graviton. One writes (see eqs. (10), (45)), Substituting this into the expression (44) one obtains, In other words, the β-function of λ vanishes (diverges) as the square-root when the frequency of the scalar graviton diverges (vanishes). This is the same qualitative behavior as in (2 + 1) dimensions 11 , see eq. (22a) of [37].

V. CONCLUSIONS
In this paper we have computed the renormalization group flow of G and λ in (3 + 1)dimensional Hořava Gravity. We used the background-field method and exploited the gauge invariance of the kinetic term which allowed us to focus on the renormalization of operators constructed from the shift field. Setting the spatial background metric to be flat we avoided major complications arising due to the involved structure of the potential term. We have considered a variety of gauge choices and found that the β-function of G depends on the gauge, as expected given that this coupling is not essential, i.e. it cannot be defined using only on-shell quantities.
On the other hand, the coupling λ is essential and its β-function given by eq. (44) is gaugeindependent. β λ vanishes on two disjoint codimension-one surfaces in the 6-dimensional parameter space of essential couplings. One of these surfaces is UV attractive along the λ-direction, whereas the other is UV repulsive. The former surface represents a potential location for UV fixed points of the full RG flow. Of course, to make a definitive conclusion about the existence of such fixed points one has to calculate β-functions of the remaining couplings of the theory that enter in the potential part of the action. This calculation will require other methods than the ones used in this paper and we plan to report it elsewhere.
As a final remark, we note that the approach described in this paper can be straightforwardly generalized to compute the running of G and λ in any number of spacetime dimensions (d + 1). Indeed, all interaction vertices used in our calculation come from the kinetic term in the action and the part of the gauge-fixing term that contains the shift field. These are universal. The dimension-dependent potential entered in the calculation only via its quadratic expansion that determined the propagators of the dynamical fields. This quadratic part can be reconstructed for arbitrary d from the requirements of invariance under the Lifshitz scaling and the linearized spatial diffeomorphisms and will be fully characterized by the combination of couplings appearing in the dispersion relations of the physical modes. oneN i , with or without spatial derivatives, and two dynamical fields: whereN ij ≡ ∂ iNj + ∂ jNi and in deriving (A1a) we eliminated some terms that combine into a total time derivative 12 . Second, we need the quartic interactions involving twoN i and two dynamical fields. In this case divergent contributions come from tadpole diagrams where the two lines of the dynamical fields stemming from the same vertex connect with each other to form a loop. Thus, to find the divergent contributions of the form (15), we can restrict only to those terms in the quartic Lagrangian that contain exactly two spatial derivatives acting on theN i -legs. These have the form, In this appendix we report the results of the computation in two gauges outside the family considered in Sec. IV. In both cases the gauge parameters σ, ξ are chosen in the way to reduce the number of pole structures by making the gauge poles coincide with the poles of the physical modes. The third parameter ς is kept arbitrary.
i) The first choice is It leads to the following β-function for the coupling G, where α is defined in (45).
ii) The second choice, yields We observe that while β G is independent of ς, it is different in the two gauges, as expected from the gauge-dependence of the renormalized Lagrangian.
On the other hand, we have found that the β-function of λ is the same in both gauges (B1), (B3) and coincides with the expression (44) from the main text. Thus gauge-independence provides a further check of our calculational procedure.
where is an infinitesimal parameter and A is a functional of the background fields which is a linear combination of the equations of motion [48,49]. Our task is to see if such a functional exists in the projectable HG. To be in line with the main text, we work in Euclidean time and focus on the high-energy part of the action given by (11).
As we are interested only in the one-loop divergence of Γ eff which is an integral of a local Lagragian density, the functional A must also have this form. Moreover, it should be invariant under FDiff and contain only marginal operators with respect to the Lifshitz scaling. Variation of (11) with respect to N i and γ ij produces the equations, where for simplicity we set N i to zero by a choice of gauge. Wo do not write explicitly the result of the variation of the potential term which is too lengthy. The l.h.s. of the first equation is a vector operator of dimension 4. Any other local covariant vector operator in the theory has at least dimension 3 (this is the case for ∇ i R). Thus, we cannot use the l.h.s. of (C2a) to construct a local Lagrangian density of dimension 6. On the other hand, the l.h.s. of (C2b) already has dimension 6. To get a scalar density under FDiff we contract it with γ ij . This yields, The last term on the l.h.s. represents the result of the variation of the potential part with respect to spatial Weyl transformations, γ ij → Ωγ ij . By inspection of the expression (11), we find that for Weyl rescalings with a space-independent factor Ω the potential part transforms as √ γ V → Ω −3/2 √ γ V. This implies that for a general rescaling the variation of √ γ V has the form, Therefore, the l.h.s. of (C3) can be written as Integrating over the whole spacetime and neglecting total derivatives one obtains, Adding it to the effective action as in (C1) induces a shift of the couplings (46) discussed in the main text. Let us mention another way to arrive to the expression (C6). Namely, one can restore time reparameterizations by reintroducing in (11) the lapse function N (t). Then variation with respect to N produces a global Hamiltonian constraint 13 , Clearly, setting N back to 1 and integrating the l.h.s. of this equation over time we recover (C6).