Renormalization Group Equation and scaling solutions for f(R) gravity in exponential parametrization

We employ the exponential parametrization of the metric and a"physical"gauge fixing procedure to write a functional flow equation for the gravitational effective average action in an $f(R)$ truncation. The background metric is a four-sphere and the coarse-graining procedure contains three free parameters. We look for scaling solutions, i.e. non-Gaussian fixed points for the function $f$. For a discrete set of values of the parameters, we find simple global solutions of quadratic polynomial form. For other values, global solutions can be found numerically. Such solutions can be extended in certain regions of parameter space and have two relevant directions. We discuss the merits and the shortcomings of this procedure.


Introduction
Since decades the ultraviolet (UV) completion of gravity is one of the major open problems in theoretical physics. Several frameworks have been proposed, such as string theory, loop quantum gravity, causal dynamical triangulation, matrix models and their generalizations such as tensor models and group field theory. Perturbative analyses in quantum field theory have shown that the gravitational interactions, considered at quantum level, are not renormalizable starting from 2 loops [1] (one loop in the presence of quantum matter [2]). In this approach a weak coupling expansion is therefore unable to give a meaningful description of the quantum gravitational interactions up to UV scales, even if, due to the smallness of the Newton constant, it can be considered a low-energy effective theory including some quantum corrections in a consistent bottom-up approach [3]. On the other hand, the possibility of having an interacting UVcomplete quantum field theory of gravitation, originated from a non-Gaussian UV fixed point (FP) in the theory space, has been proposed by Weinberg [4] and fits naturally in the nonperturbative renormalization group framework. The existence of such a FP (scaling solution) would permit to have a renormalization group (RG) trajectory in the theory space characterized by all the dimensionless couplings remaining finite when the UV cutoff is removed.
In order to make the critical theory at the FP physically meaningful, the corresponding conformal field theory should have a finite number of relevant deformations so that a finite number of measurements at low energy would be sufficient to completely fix the theory, i.e. all the (infinitely many) couplings at any scale. To decide if this latter fundamental property may be achieved, any kind of investigation in this direction should start from a sufficiently large (possibly infinite dimensional) theory space, show the existence of one or more critical theories and eventually should be able to identify which one has the required properties under deformations that may lead to a definition of a UV complete and predictive theory.
The search of a gravitational FP has been conducted mostly using some approximation of the functional renormalization group equation (FRGE). In most cases, finite truncations (i.e. finitely many couplings) were studied and the general strategy was to establish that the addition of new couplings would not affect too strongly the results [5,6,7,8,9]. While the number of couplings considered simultaneously has become quite high, these analyses still fall short of exploiting the full power of the FRGE machinery, which, as the name suggests, is designed to deal with the renormalization of whole functions, or equivalently infinitely many couplings.
The study of truly functional truncations has been slowly gathering momentum in the last few years. In the case of f (R) truncations functional flow equations for the functions f had already been written in [7,8], but the first serious attempt to solve the functional FP equations was made in [10]. Negative results concerning the equations written in [7,8] have been reported in [11,12]. It has been argued that the solution may require truncations that go beyond functionals of the background field alone, which in this framework one cannot avoid to introduce, and taking into account the Ward identities of the quantum-background split symmetry. This is a general issue that goes beyond the f (R) truncations and progress in this direction has been made in [13,14]. In the meantime solutions were found in simplified (lower-dimensional and/or conformally reduced) settings [15,16] and recently also in the full four-dimensional case [17]. A different class of functional truncations of a scalar-tensor theory consists of actions of the form In this case, the system is closer to familiar models of scalar theories and one may hope to be able to use the experience gained there. In particular, given that a FP is known to exist in the Einstein-Hilbert truncation of three dimensional gravity, and a non trivial functional FP exists for pure scalar theory, it would seem reasonable to expect a functional FP for the combined system. However, nothing resembling such a solution was found in [18]. The reason for this has been discussed in [19] and a different flow equation was proposed, based on the use of an exponential parametrization for the metric. In this case, solutions of the corresponding FP equations have been found [19,20,21]. The use of the exponential parametrization can be motivated by the nonlinear nature of the metric itself since, for example in the case of euclidean signature, the set of all metrics at a point can be identified with the coset space GL(d)/O(d).
It is natural to use for this nonlinear space an exponential map [19,22] and, as we shall show, in such a case the off-shell effective action expanded in perturbations appears to be simpler than the one obtained using a linear parametrization. The same method has been used also in f (R) theory [23] and again some scaling solutions could be found in closed form, at least for specific cutoff procedures. In the present paper, we expand the analysis of [23] and also give details for the derivation of the equations. We then discuss the relation between these exact solutions and numerical solutions for generic values of the parameters, for which we give a constructive procedure.
In Sect. 2, we begin by deriving the Hessian for the f (R) gravity in d dimensions (Sect. 2.1), and discuss gauge fixing and the resulting determinants in Sect. 2.2. In Sect. 3, we give the details of the derivation of the flow equation. In particular in Sect. 3.1 we first set up the general form of the flow equation based on a coarse-graining procedure, dependent, for each irreducible spin component, on an endomorphism parameter, and use heat kernel to evaluate it in Sect. 3.2. Up to this point, our discussions are for arbitrary dimensions, but they are restricted to four dimensions from Sect. 3.3 onwards. We then discuss the alternative flow equation derived by the spectral sum in Sect. 3.4. In Sect. 4, we present global solutions of the spectral-sum based flow equation, first quadratic ones in Sect. 4.1 and then more general numerical solutions in Sect. 4.2. In Sect. 5, we also present and discuss global solutions of the flow equation derived by the heat kernel. Since the structure of the flow equations are basically the same, we only describe the solutions briefly. Finally in Sect. 6, we summarize our results and discuss the limitations of the present results and possible ways of overcoming them.

f (R) truncation
We use the background field method and split the metric into background and quantum parts. We perform the calculations by different parametrizations of the metric fluctuation. One is the linear split g µν =ḡ µν + h µν , (2.1) and the other is the nonlinear exponential type: Henceforth we assume that the background space is four-dimensional Einstein space with

Hessian
For the linear split (2.1), the Hessian is given by Here and in what follows, we suppress the overall factor √ḡ and the covariant derivative ∇ µ is constructed with the background metric. In the exponential parametrization (2.2), the metric has an infinite series expansion Therefore the linear term in the expansion of the action, which is generates an additional contribution to the Hessian by the replacement Substituting the York decomposition (B.1) into (2.4), we find where ∆ 2 , ∆ 1 and ∆ 0 are the Lichnerowicz Laplacians on the symmetric tensor, vector and scalar respectively, defined in Appendix A. This formula agrees with [7]. When the exponential parametrization (2.2) is used, we find Remarkably all terms containing ξ µ cancel out. For f (R) =R 2 , this agrees with our previous result [24,25]. If we use the gauge-invariant variable s = h + ∆ 0 σ, Eq. (2.9) can be rewritten as Note that the last term, which is the only non-gauge invariant one, is proportional to the field equation. It vanishes for F (R) =R d/2 , which isR 2 for four dimensions, because the action is scale invariant in this case.

Gauge fixing
Let us consider a standard gauge fixing term Using the York decomposition, this reduces to (2.14) We see that a specific combination of scalar degrees of freedom appears in this formula. Following [26], it is convenient to reparametrize the scalar sector in terms of the gauge-invariant variable s and a new degree of freedom χ defined as The gauge fixing function then reads Thus the gauge fixing action becomes From (2.15), we see that χ transforms in the same way as σ. On shell, the last term in (2.11) is zero, so the quadratic part of the action is written entirely in terms of the physical degrees of freedom h T T and s, and the gauge fixing entirely in terms of the gauge degrees of freedom ξ and χ. The ghost action for this gauge fixing contains a non-minimal operator Let us decompose the ghost into transverse and longitudinal parts 19) and the same forC. (This change of variables has unit Jacobian). The ghost action splits in two terms On shell, the h-h term in (2.11) goes away. The one-loop partition function is the product of several determinants. The determinants associated to the gauge invariant variables (h T T , s) are manifestly gauge independent. The fields (ξ, χ) in the gauge fixing term contribute (2.21) Then there are the ghost determinants and finally we have the Jacobians. The York decomposition has Jacobian Det ∆ 1 − 2 dR whereas the subsequent transformation (σ, h) → (s, χ) has unit Jacobian. Altogether, on shell, the gauge-dependent factors cancel out, and the result is gauge independent. The one-loop determinants are Note that in General Relativity, the second determinant in the denominator becomes a constant. The result then agrees with [27]. The trivial scalar determinant is a sign that the scalar degree of freedom s is non-propagating. In general, in f (R) gravity this degree of freedom propagates and contributes to the one-loop effective action. Off shell, things are more complicated. One can solve Then, the last term in (2.11) contributes additional terms to the s-and χ-operators (as well as a s-χ mixing term). There are two choices that simplify the situation. For β = 0, in which case χ = σ, we have h = s − ∆ 0 χ. In the case β → ∞ we have simply h = −∆ 0 χ. In this case the last term of (2.11) leaves the gauge-invariant modes h T T and s untouched and only produces a term proportional to χ 2 . However, for β → ∞ the gauge fixing strongly enforces the condition χ = 0, independently of α, so in this gauge the quadratic action becomes the same off-shell as on-shell. We can further kill the ξ µ gauge degree of freedom by choosing α = 0. This gauge should then be equivalent to what was called the "unimodular physical gauge" in [19].
In this gauge, we have h = 0, so s = ∆ 0 σ and (2.10) reduces to the first two lines in (2.9). As discussed in [19], the gauge conditions ξ µ ≡ ∆ 1 − 2R d ξ µ = 0 and h = 0 produce two ghost determinants In this case the spin-one determinant in the Jacobian of the York decomposition (2.23) is cancelled by the Jacobian of the transformation ξ µ → ξ µ . Collecting the rest, we see that this gauge reproduces the on-shell determinants (2.24). The only difference is that off shell the background curvature is generic, whereas on shell it becomes a function of the couplings in f . We take this equivalence of the one-loop effective action to the on-shell result be a distinct advantage of the exponential parametrization and of our gauge choice, in the sense that using this procedure the off-shell results are less sensitive to the contributions of unphysical degrees of freedom.
3 Flow equations

Cutoff and functional renormalization group equation
For the definition of the coarse-graining, we have to choose some reference operator. In the Hessian on the four-sphere (2.11), the operator ∆ = −∇ 2 appears everywhere and is a natural choice. However, in order to gain some additional freedom, we follow [16,17] and add to ∆ terms proportional to the scalar curvature, with coefficients −α, −γ and −β for spin two, one and zero, respectively. These parameters should not be confused with the gauge fixing parameters which do not appear in the following. By the standard procedure, we then get the FRGĖ where the dot denotes the logarithmic derivative with respect to the scale k and P k (z) = z + R k (z), with the cutoff function R k (z). Because of the implicit dependence on f of the coarse-graining scheme, we refer to it as spectrally adjusted. The subscripts on the traces represent contributions from different spin sectors.
In order to guarantee that the operators ∆ − αR, ∆ − βR, ∆ − γR have positive spectrum, the parameters α, β, γ should satisfy certain bounds. The spectrum of ∆ is given in Appendix C. We recall that the = 1 modes of ∆ acting on spin one fields are Killing vectors, so that they do not contribute to the spectrum of h µν and have to be left out. For the same reason, the modes = 0 and = 1 of ∆ acting on scalars also have to be left out. Thus all spectra begin with = 2.
In four dimensions, requiring that the modes = 2 have positive eigenvalues leads to the conditions

Heat kernel evaluation
The evaluation of the traces is done as follows: First, for some differential operator z, consider for the spin j sector, whereW (s) is the inverse Laplace transform of W (z): Using the heat kernel expansion We choose the optimized cutoff profile [28] . For the contribution of the spin-two modes in (3.1), we find After moving to the dimensionless quantities, Similarly we find for spin 1 while for spin 0, we have Finally we recall that on the sphere the constant mode of h cannot be considered a gauge degree of freedom and therefore cannot be gauge fixed. One can choose a coarse-graining scheme where this infrared mode never appears in the flow equation. Its contributions would have to be added directly to the effective action only in the deep infrared at k = 0. Alternatively, we can consider the scheme where we add to the r.h.s. of Eq. (3.1) the following term: . (3.12)

Four dimensions
The heat kernel coefficients b 2n for ∆ acting on spin-two, one and zero are given in [8,17] for type I cutoff. We extend the calculation to our case and give the results in Appendix C. Substituting these heat kernel coefficients and equations (3.9), (3.10) and (3.11) in (3.1), we obtain where (3.14) If we include the contribution of the constant mode of trace h, we have an additional term to the r.h.s. of Eq. (3.13) 8 3 which has been obtained dividing the expression in Eq. (3.12) for d = 4 by the corresponding volume of the sphere and passing to dimensionless quantities.

Spectral sum approach for d = 4
Alternatively, we can compute the traces in Eq. (3.1) by summing directly the corresponding functions of the eigenvalues of the Laplacian on the sphere [29] TrW where the eigenvalues and the corresponding multiplicities are given in Table 3 in Appendix C and we recall that all the mode sums have to start from = 2. For simplicity, we restrict our analysis to four dimensions. We shall use the same optimized cutoff so that the support of R k (λ l (d, s) + E) will be restricted to the modes ≤¯ , where the upper bound is determined by the condition λ + E ≤ k 2 . In particular, for the different spins we have the followinḡ R-dependent restrictions: The sums extend up to the integer part of these upper bounds. In this way one would obtain a discontinuous structure in the flow equation. Since the sum can be written as an exact function (a polynomial) of these upper bounds, we keep them as real variables. To reduce the error, we follow [10] and perform the average of the sums taken up to¯ (s) and¯ (s) −1. This average has also the nice property of removing the square roots from all the spectral sums. We notice that the lower limit in the sums requires¯ (s) ≥ 2 in order to have a contribution from the integration of the quantum fluctuations. We shall discuss the consequences of this in Sect. 6. Using Eq. (3.16), we then obtain the following flow equation in terms of dimensionless variableṡ One can notice that the structure of the equations (3.13) and (3.18) is the same. In fact, (3.18) can be rewritten exactly in the form (3.13), with the coefficients 4 Solutions of the spectral-sum based equation We analyse the FP solutions of the flow for d = 4. In Ref. [23], we have presented exact quadratic solutions for both (3.13) and (3.18), so we first discuss this kind of solutions.

Global quadratic solutions
There are solutions which are quadratic polynomial in the curvature and exist for a finite, discrete set of values of α, β and γ. They are obtained by plugging into the FP equation the ansatz and writing the equation as N D = 0. Here N is a polynomial of fifth order in r and N = 0 can be solved for the six unknowns α, β, γ, g 0 , g 1 and g 2 . We then find the following distinct solutions (Table 1), where the unknowns were evaluated numerically:  Table 1: Quadratic solutions of the spectral sum FP equation. In the last column, we report the results for the positive critical exponents, evaluated in an ninth-order polynomial expansion. The critical exponent 4 is present in all solutions and is related to the cosmological term. Those in the second line converge slowly and may not be accurate estimates. The fixed points in the last two lines have eigenvalues greater than 4 and are not reliable.
For these values of the couplings and endomorphisms, one has then to check the behavior of the denominator D. Since D has some zeros for at least one positive value of r, one may worry that the solutions might not exist at these possible singular points. However N vanishes identically as a function of r, so the residues of N/D at the zeros of D are zero, implying that the solutions are valid through these points.
Actually we notice that the solutions (α, β, γ, g 0 , g 1 , g 2 ) which are given numerically in Table 1, have simple non-numerical forms Indeed using these expressions, one may also work directly with Eq. (3.18) and check explicitly that upon substitution of any of the five solutions, the residues at all the possible singular points are zero after cancellations of several terms. Therefore these solutions are defined globally, and even extensible on the full real line.
It is instructive to further examine the possible singularities of Eq. (3.18) and their cancellation. This will be useful in the numerical search of global solutions for more generic values of α, β and γ to be discussed in the next subsection. Since in all the above solutions the function ϕ(r) has a minimum, one expects that the term in (3.18) containing d 2 /ϕ may be a possible source of singularities. We can check that these possible singularities are actually resolved by the zeros in the polynomial d 2 . Let us see this more concretely. The position of the minima for the 5 solutions are On the other hand, the zeros of d 2 are located atr(α) such that It is easy to check that one of these two zeros is located precisely at the minimum for the first four solutions, providing a simple mechanism of cancellation of the singularity. For the last solution, this is not the case and the cancellation of the singularity is more tricky and involves other three terms of Eq. (3.18). Actually there are also other possible fixed singularities, one which depends on α at r 1 = − 6 6α+1 and a second one, coming from the vector ghost term which leads to a γ-dependent fixed singularity at r 2 = − 4 4γ−1 . We also note that for the quadratic solutions, there is no contributions from the term ϕ = 0, which otherwise could have given rise to other fixed singularities.
An analysis of the eigenperturbations of the solution, based on a polynomial expansion around the origin, leads to the following conclusions. The first solution has two relevant directions with critical exponents θ (1) i = (4, 2.03) while the second one has two relevant direction θ (2) i = (4, 3.2) and possibly a third one with critical exponent close to zero. If the perturbation is at most quadratic in r, then its existence (at least for positive r) is subject to the same constraint as the solution itself. This is the case for δϕ(r) = 1, which is the eigenperturbation with eigenvalue λ = −4. We do not find other simple eigenperturbations of a finite polynomial form. Since the other eigenperturbations cannot be purely quadratic, the fixed singularities corresponding to the zeros of the coefficient of ϕ will appear. We find that for the first, fourth and fifth solutions, the values of α, β and γ lead to a linearized third order equation with three fixed singularities. Therefore there may be globally defined eigenperturbations in the domain r ≥ 0. On the other hand, the second and third solutions have values of α, β and γ corresponding to four fixed singularities in the differential equation for the eigenperturbations, which therefore cannot exist globally, but only up to some maximum value of r.
It follows from this discussion that for any infinitesimal change in the parameters α, β and γ, the equation would have at least three fixed singularities and another singularity at the minimum, which precludes the existence of a global solution. As we shall discuss in Sect. 6, the physical significance of the solution beyond the minimum is questionable, so this result may not be so negative after all.

Global numerical solutions
We now analyse the FP of the flow given in Eq. (3.18) at fixed endomorphisms, and search for global numerical solutions compatible with its analytical structure.
First of all, it is convenient to cast the FP equation in the normal form. Evidently the equation becomes singular at the zeros of the coefficient of ϕ . The location of these zeros varies with β but otherwise is fixed: If these zeros happen to be at positive r, the possible singularities there must be cancelled.
We have already seen that the potential singularities arise if there is an extremum in ϕ, but a regular solution can exist if the minimum point r min is located at a zero of the coefficient d 2 , which depends on α. There are also two possible fixed singularities at r = − 6 6α+1 and at r = − 4 4γ−1 . Since the FP equation is of third order, its solutions contain three parameters. Regularity condition at a fixed singularity requires that the residue at the singularity should vanish, and this fixes one of the three parameters. Eventually there might also be moving singularities, as for example happens to be in the local potential approximation analysis for an interacting 3-dimensional scalar field theory with Z 2 symmetry, where a unique solution is selected, the Wilson-Fisher FP of the Ising model. We show in Fig. 1 the positions of the singularities and possible minimum as a function of α, β and γ. If the solution has no minimum and there is no moving singularities, we can obtain a finite set of solutions by requiring the presence of just two singularities at strictly positive r in addition to the singularity at r = 0. This leads to three constraints which maximally restrict the initial value problem. This situation could be realized in several ways, provided the three singularities are kept different: either In the following, we search for a global solution with a shape similar to the analytic quadratic solution previously found, which are characterized by a minimum. We have already noted that if there is a minimum in ϕ, the singularity there must be cancelled by a zero in d 2 . Therefore the parameter α must be chosen such that at least one of the zerosr(α) of d 2 is positive. Thus, avoiding the presence of the fixed singularity shown as a dotted curve in Fig. 1 In the absence of moving singularity, we can then allow for other two fixed singularities which would fix the solution uniquely. In order to have this, we should require 1 3 ≤ β ≤ 5 6 and γ ≥ 1 4 or β ≥ 5 6 and γ < 1 4 , as it is evident from Fig. 1 (a) and (c).
With these considerations, in what follows, we analyse in detail a specific example belonging to the first region, This is chosen just for the purpose of illustration of our study, and our following construction of the solution should go through for other choices of the parameters if we choose the endomorphism parameters within the ranges specified above. For the choice of (4.6), there are fixed singularities at r = 0 and r = 2 and the position of the minimum can take the values r min = 6/5 or r min = 3. It is convenient to employ first a simple polynomial truncation around the origin to find an approximated form of the solution in the vicinity of the origin ϕ(r) = N m=0 g m r m . (4.7) We have performed a scan of the solutions of the FP equation at different values of N up to N = 16, investigated also the linearized evolution around the FP, and then diagonalized the corresponding stability matrix. The scaling exponents are extracted as its eigenvalues with the opposite sign. For orders N > 10, we find a high degree of stability. Here we give the first coefficients of the polynomial solution of order N = 16: ϕ p (r) = 0.00410949−0.00809798 r+0.00392097 r 2 −0.000254811 r 3 −3.35023×10 −6 r 4 −8.34593×10 −6 r 5 +· · · (4.8) The plot of the polynomial solutions at different order is shown in Fig. 2.
The pattern of the eigenvalues is pretty stable. We find that there are just two relevant directions with eigenvalues −4, −1.83, +1.38, 4.41, · · · and increasing the order of the truncation leads to the appearance of new irrelevant directions.
We also observe that the polynomial solutions have a minimum converging exactly to the point 6/5, which is one of the values necessary to remove the corresponding singularity. Therefore in our numerical search for the global solution we impose this property.
Our strategy to construct a global numerical solution is then the following: 1. We construct three polynomial expansions as functions of two free parameters at each of the three possible singular points r = 0, 6/5, 2 by imposing locally the regularity condition there. In particular we choose the following set of free parameters: (ϕ (0), ϕ (0)) , (ϕ(6/5), ϕ (6/5)) and (ϕ (2), ϕ (2)), for the polynomial expansion at r = 0, r = 6/5 and r = 2, respectively. 2. We then evolve ϕ(r) numerically from r = 0+ to r = 6/5 and using the shooting method, we search for a subset of the parameter plane (ϕ (0), ϕ (0)) which gives solution with a minimum at 6/5, where the possible singularity is cancelled by the zero of d 2 . This leads to a one-dimensional curve in the parameter plane. We repeat the same procedure for the evolution of ϕ(r) from r = 2− to r = 6/5. We then find also a one-dimensional curve, this time in the parameter space (ϕ (2), ϕ (2)).
We show the results of this process in Fig. 3 (a) and (b). We note that in the plane (ϕ (2), ϕ (2)) there are regions where one encounters moving singularities, but these regions do not overlap with the ones satisfying the regularity conditions. 3. We map the two curves in (ϕ (0), ϕ (0)) and in (ϕ (2), ϕ (2)) obtained as conditions to cancel the singularities into the plane of parameters (ϕ(6/5), ϕ (6/5)). This can be done by numerical evolution.
5. Then we use the regular polynomial expressions to fix the boundary conditions around the singular points, and construct the solution by evolving ϕ(r) numerically from 6/5− to 0, from 6/5+ to 2− and from 2+ to the right. The result is presented in Fig. 4 (a) and (b).
6. Finally we consider the leading asymptotic behavior as a function of a free real parameter A which we tune to obtain a match with the numerical solution. The result is shown in Fig. 4 (c).
We give here the first few terms of the asymptotic expansion: We have repeated the analysis for the same values of α and β but with γ = 1 4 + 10 −2 and found that the results are very similar, with just a mild deformation. Actually we find that there is a small difficulty in finding an acceptable asymptotic match in this case. This is due to the fact that at γ = 1 4 , the asymptotic behavior changes to ∼ r 2 log(r) with a negative coefficient, and close to this value of γ the radius of convergence of the asymptotic expansion seems to become much smaller.
In general there are connected regions in the space of the parameters α, β, γ inside which the global solutions can be continuously deformed. We have checked that similar solutions can be obtained for α = − 1 6 , β = 1 3 and any γ > 1 4 . They look really as members of a single family of solutions which are connected by continuous deformation. Also the estimated values of the critical exponents show a mild dependence on γ.
One has also to check that the eigenperturbations should not be redundant [12]. In particular in an f (R) truncation based on a background sphere, an operator is redundant if it can be written in the following form i.e. it is proportional to the equation of motion at the fixed point with a well defined a(r) in the full domain of interest. If the equation of motion a the fixed point have a nontrivial solution within the physical domain, then a regular eigenperturbation cannot be redundant. We find that the equation of motion are always satisfied by a positive value of r such that r < r min , implying that the eigenperturbations are essential. For both cutoff types the fixed points have 2 relevant directions. This is in sharp contrast to the polynomial solutions found previously [6,7,8,9] that had three relevant directions.

Type I and type II cutoff
We also note that [30] finds, in unimodular gravity, polynomial solutions that have two relevant directions. Since the cosmological term is absent in unimodular gravity, this would presumably correspond to three relevant directions in the full theory. The source of this difference cannot be in the parametrization, which is exponential in both cases.
According to the count of singularities discussed earlier, the type-I cutoff has four fixed singularities and therefore is not expected to admit global solutions. Most likely, the polynomial solution (4.11) given above is valid only in a neighborhood of the origin. The type-II cutoff has two fixed singularities due to β (in r = 0 and r = 2), so if one demands the existence of a minimum the solution would be completely determined. The solution (4.13) given above has a minimum and could be an approximation of a global solution. The solution (4.12) does not have a minimum, but it seems to develop a stationary inflection point with increasing order of the polynomial. For the sake of comparison we note that the solutions reported in [30] have a minimum. In the preceding numerical analysis we have studied cutoffs that differ from type II only in having γ > 1/4. This is because the exceptional value γ = 1/4 leads to a more complicated asymptotic behavior for large r.

Some results for the case of the Heat kernel equation
A similar analysis can be carried out for the flow equation Eq. (3.13) obtained from a heat kernel expansion, and for the linear variations around it. Again there exist a discrete set of endomorphisms parametrized by α, β and γ for which purely quadratic solutions exist. They were reported already in [23] and are given here in Table 2. The critical exponents have been obtained by a polynomial truncation. These numbers should only be taken as rough estimates.   In the last column, we report the results for the positive critical exponents, evaluated in an seventh-order polynomial expansion. The critical exponent 4 is present in all solutions and is related to the cosmological term. Those in lines 3 and 4 converge slowly and may not be accurate estimates. The fixed point in the last line has eigenvalues greater than 4 and is not reliable.
The situation is the same as for the equation obtained from a spectral sum. The only difference is that the values of the singular points have more involved expressions, and this is the reason why we presented detailed analysis for the spectral-sum equation. We find that these solutions are indeed globally defined, all residues at any singularity of the equation being zero.
As before, all the FP solutions have one trivial relevant direction associated to a constant eigenperturbation with eigenvalue −4. No other deformations which are polynomials in r of order no greater than two exist. Then we find that for the last solution the third order differential equation for the eigenperturbations is characterized by three fixed singularities while for all the others there are four or more fixed singularities so that no global eigenperturbations can exist. As already mentioned in the case of the spectral sum equation, this precludes the existence of global solutions in an infinitesimal neighborhood of the exact quadratic solutions. The continuous interpolations between the exact solutions that we have discussed in [23] are therefore necessarily restricted to a finite range of r. See Sect. 6 for a discussion of this point.
Let us now address the question how to find global numerical solutions using the same reasoning developed in Sect. 4. The zeros of coefficient ϕ of Eq. (3.13) are not simple analytically, so we give some numerical bounds. It is easy to see that for 0.394474 < β < 0.678204, one has only two fixed singularities at r = 0 and r = r 0 (β). Then the coefficient c 1 of 1/ϕ has at least one positive zero for α < 0.47552. As we have already seen, the zeros of c 1 in (3.14) naturally provide a simple mechanism for the cancellation of a possible singularity located at the minimum of ϕ(r), if present. In such a scenario with the presence of a minimum, which we have seen to be favoured by the dynamics in several cases, we must require that no other singularities be present. These may arise from the poles at r = − 6 6α+1 and at r = − 4 4γ−1 , as for the spectral sum case. Therefore we must require also α ≥ − 1 6 and γ ≥ 1 4 . As an example, we consider a case close to the one of Sect. 4.2 with the choice α = − 1 6 , β = 1 2 and γ = 1 2 . For such values of the parameters, we expect a minimum at r min 1.55718. We have studied the polynomial solutions around the origin, which represent reasonable approximation of the global solutions. The polynomial solutions have a minimum converging to the point r min 1.55718. We give here the first ones for the polynomial solution of order N = 16.
(5.1) In Fig. 5 (a), we show how the position of the local minimum of the polynomial solutions at different order N rapidly converge to the expected value of r min given above, while in Fig. 5 (b) we give the leading five eigenvalues (which include the two negative ones associated to the relevant deformations) for each polynomial solution available up to N = 11. In particular the nontrivial relevant direction for a polynomial approximation of order N = 11 is associated to λ = −1.799. Again when we increase the order N , we find new irrelevant deformations. The pattern at lower order is not much altered and numerical values of the lower eigenvalues are little sensitive to the increase of the order of the analysis. In this way we can construct global solutions.

Discussion
The present investigation is a step forward in the search of scaling solutions for f (R) gravity. It had been shown in [11,12] that the flow equations for f written in [7,8,10], either do not have scaling solutions, or the solutions are such that all their perturbations are redundant. It is therefore encouraging that the equations studied here do have some global scaling solutions. It is particularly striking that for some choices of cutoff, the scaling solutions have an extremely simple, quadratic form. This is probably due to the much simpler form of the equations. Similar equations have also been written in [16], and in [30] for unimodular gravity, and it is worthwhile to comment on the differences between them. Our equation, as well as the one for unimodular gravity, has the property that the function f does not appear in its beta function in undifferentiated form. It only appears undifferentiated in the flow equation for the dimensionless variable ϕ, through the classical scaling term. In the Einstein-Hilbert truncation, this fact corresponds to the absence of the cosmological constant from the beta function of Newton's coupling. This is related to the absence of trace fluctuations in both calculations -in unimodular gravity by definition of the theory and in our calculation by gauge choice. This is in turn due to the fact that in the exponential parametrization only the volume term contributes to the Hessian of the trace. The other difference is in the gauge and ghost sector. In our equation, there are contributions due to a spin one and a spin zero ghost. The spin zero ghost is absent in unimodular gravity because there is no need to impose the tracelessness of h µν as a gauge condition. On the other hand, no ghosts are present in [16]. In spite of these differences, the equation of [16] was shown in [17] to have a global scaling solution whose general shape is quite similar to ours. The equation for unimodular gravity has only been analyzed at polynomial level so far.
A potential problem of higher derivative gravity theories is the presence of physical ghosts. This is usually analyzed in a perturbative framework, but the significance of this analysis in the neighborhood of an interacting fixed point is questionable. Several mechanisms have been proposed that could get around this issue. In any case, we notice that in our analysis, which is restricted to an f (R) truncation in a de Sitter background, no physical ghosts are present.
The significance of the positive results we have found is reduced by several circumstances. The first is the restriction of the action to purely background-dependent terms, the so-called "single metric truncation". The effective action at finite cutoff cannot be a function of a single metric, so the classical invariance under the "shift symmetry"ḡ µν →ḡ µν + µν , h µν → h µν − µν (in the linear parametrization) is broken. 4 The dangers of the single-field truncations have been discussed in [31,32]. One should therefore consider truncations involving either two metrics [13] or the background metric and a fluctuation field [33,34,35] or else solve the flow equation together with the modified Ward identities of split symmetry [14]. In any case the scaling solutions found here can be at best an approximation of a genuine scaling solution.
A major difference between our results and those obtained earlier for the f (R) truncation [7,8,9] is the number of relevant directions. Previously there were three relevant directions (of which two had a complex pair of eigenvalues) involving a strong mixing of all operators, but mainly the cosmological term, the Hilbert term and the R 2 term. Due to the exponential parametrization, the cosmological constant does not appear in our beta functions and is therefore an isolated relevant direction with eigenvalue −4. Moreover, we only find one additional real relevant direction instead of two. This hints at some genuine difference between the results, but it is not clear at this point whether this is due to the different parametrization or we are looking at a different solution. In the section devoted to the search of numerical solutions we have discussed the regions in the parameter space, corresponding to different endomorphisms (coarse graining schemes), which may allow the existence of different global solutions. We note that in unimodular f (R) gravity, Ref. [30] found two relevant directions for a polynomial approximation of the fixed point. Given that the cosmological term is absent in unimodular gravity, this agrees with the "old" counting, in spite of the use of the exponential parametrization. It would be interesting to know the spectrum of eigenperturbations of the solution in [17], which was obtained in the linear parametrization.
In connection with the exponential parametrization, we also remark that while desirable, parametrization independence is in practice hard to maintain at the quantum level. In the context of asymptotically safe gravity the Vilkovisky-de Witt formalism has been used in [36]. Our philosophy here is that one could instead try to exploit the parametrization dependence to one's advantage by choosing a parametrization that reduces the dependence of the off-shell effective action on other choices. 5 It seems to us that the exponential parametrization has such virtues. Whether it defines the same quantum theory as the linear parametrization is a question that we cannot answer for the time being.
Another issue is the use of a "spectrally adjusted" cutoff, which means that the cutoff contains the function f . At a conceptual level, it is desirable to have a definition of coarse-graining that is independent of all couplings. With a spectrally adjusted cutoff, on the other hand, the definition of what one means by high and low momentum modes changes in the course of the flow. A practical consequence of this choice is the appearance of the third derivative in the flow equation. The reason for using these cutoffs 4 In the exponential parametrization, the split symmetry has a more complicated form. 5 For related remarks, see also [37,38].
in spite of these issues is that they lead to simpler equations. It would however be very important to consider a cutoff that does not have this feature [39].
A third issue had already been discussed earlier [15], and is related to the compactness of our background manifold. What is the meaning of coarse-graining on length scales that are larger than the size of the manifold? Otherwise said, if the spectrum of all operators has a finite gap δ, then for k < δ the flow equation does not integrate out any modes. If we imagine that the sphere has a fixed size ∼R −1/2 , then for r 1, one is in the regime where k 2 R . This puts into question the physical meaning of the behavior of the scaling solution for large r.
To make this more precise, recall that the mode sums begin at = 2 and consider again Eq. (3.17). If¯ becomes smaller than 2, then no modes are integrated out. This happens for respectively. Thus for fixed α, β, γ, the spin two, zero and one modes do not contribute to the flow of the function ϕ for respectively. There is no mode contributing to the flow equation for values of r such that all three bounds are satisfied.
Conversely we could demand that α, β and γ be chosen such that for all r there is some mode contributing to the flow equation. This leads to requiring that the inequalities (3.2) be violated. The meaning of this is clear: one is precisely demanding that even when k → 0 there is some mode contributing to the flow, and this is equivalent to requiring that there are some modes with zero or negative eigenvalue. Such procedure would clearly be artificial, because it amounts to a shift of the eigenvalues which does not change the eigenfunctions. Even if the constant mode has formally infinite wavelength, it would still be the case that the largest physical length scale that one can talk of in such a space is its diameter.
It is sometimes said that the equation obtained from the spectral sum method is superior to the one obtained from the use of the heat kernel expansion, because the heat kernel equation uses an asymptotic expansion that makes it valid only for small r. On the other hand, the equation obtained from the spectral sum, at least when used for a compact background and in conjunction with the optimized cutoff, requires a smoothing procedure that introduces an element of arbitrariness, which is especially strong at the low end of the spectrum. In view of the preceding considerations, the behavior at small k is in any case of dubious physical significance, so the difference between the two procedures is probably not so important.
An apparently unrelated issue is that the propagator of the spin-two mode has the wrong sign when f > 0, which occurs to the right of the minimum of our scaling solution. It is interesting to note that for each scaling solution that we have found, the spin-two bound in (6.2) is precisely saturated at the minimum. This suggests that the solution should only be considered physical in the region between zero and the minimum 0 ≤ r =R/k 2 ≤ r min . (6.3) In this region the issue of the sign of the propagator is not present. It appears from this discussion that one important source of ambiguities is the compactness of the background manifold. If one wants to understand better the existence of global scaling solutions in f (R) gravity within the realm of single-metric truncations, the best course seems to be to consider a non-compact background.

B York decomposition
The York decomposition is defined by When this is squared, we get Note that we can freely insert the covariant derivatives inside the above expression.
C Heat kernel coefficients on the d-sphere The curvature tensors satisfȳ The volume of the d-sphere is The heat kernel coefficients can be found by summing over eigenvalues λ (d, s) of the operator ∆ weighted by their multiplicity M (d, s)  Table 3. We use the Euler-MacLaurin formula  mode n = 1 for the spin one field ξ µ (Killing vectors) and for the field σ one has to leave out the modes n = 0 (constant) and n = 1 (related to the five conformal Killing vectors that are not Killing vectors), so the sum should start from n = 2. The integral parts in (C.4) are given by Note that these coefficients b 2n for n ≥ 2 disagree with [8,17] even if we set α = β = γ = 0, but this is because there the contributions from the Killing were included and the sum was taken from a = 0 (s = 0), a = 1 (s = 1), a = 2 (s = 2).