Asymptotic safety of quantum gravity beyond Ricci scalars

We investigate the asymptotic safety conjecture for quantum gravity including curvature invariants beyond Ricci scalars. Our strategy is put to work for families of gravitational actions which depend on functions of the Ricci scalar, the Ricci tensor, and products thereof. Combining functional renormalisation with high order polynomial approximations and full numerical integration we derive the renormalisation group flow for all couplings and analyse their fixed points, scaling exponents, and the fixed point effective action as a function of the background Ricci curvature. The theory is characterised by three relevant couplings. Higher-dimensional couplings show near-Gaussian scaling with increasing canonical mass dimension. We find that Ricci tensor invariants stabilise the UV fixed point and lead to a rapid convergence of polynomial approximations. We apply our results to models for cosmology and establish that the gravitational fixed point admits inflationary solutions. We also compare findings with those from $f(R)$-type theories in the same approximation and pin-point the key new effects due to Ricci tensor interactions. Implications for the asymptotic safety conjecture of gravity are indicated.

It is widely acknowledged that for quantum field theories to be fundamental and predictive up to highest energies, their short-distance behaviour should be controlled by an ultraviolet (UV) fixed point under the renormalisation group [1,2]. In particle physics, it is known since long that ultraviolet fixed points can be free, such as in asymptotic freedom of non-abelian gauge theories [3,4]. Ultraviolet fixed points can also be interacting, such as in asymptotic safety [5]. In the absence of gravity, general theorems for asymptotic safety [6] and proofs of existence in four dimensional theories are available, covering simple [7,8], semi-simple [9], supersymmetric gauge theories [10], and Standard Model extensions [11]. It has also been speculated that a quantum theory for gravity with matter may exist with the help of an interacting UV fixed point [12][13][14][15][16][17][18].
In gravity, and much unlike in non-Abelian gauge theories, proofs or theorems for asymptotic safety are presently not at hand. This is largely due to the dimensional nature of Newton's coupling whose canonical mass dimension must be compensated by large anomalous dimensions and strong quantum effects. Still, at strong coupling, the asymptotic safety conjecture can be tested systematically with the help of a bootstrap search strategy [19]. By now, an extensive body of evidence for asymptotic safety has grown over the years using the renormalisation group within increasingly sophisticated approximations . Invariably, the fluctuations of the metric field are found to be strong, strong enough for gravity to become anti-screening [12][13][14]. Canonical power counting continues to be a "good" ordering principle [5,19]: canonically relevant couplings remain relevant at interacting fixed points, canonically irrelevant couplings remain irrelevant, and the spectrum becomes near-Gaussian with increasing mass dimension [19,57]. It would then seem that asymptotically safe gravity becomes "as Gaussian as it gets".
Much of the above picture has been established for f (R) quantum gravity and high-order polynomial approximations thereof. In this work, we extend investigations of the asymptotic safety conjecture towards actions beyond Ricci scalars. We use the method of functional renormalisation to derive renormalisation group flow equation for gravity including at strong coupling. Our primary goal is to generalise the successful flow equations for f (R) type approximations to much more general classes of actions involving Ricci and Riemann tensor invariants, or products and combinations thereof. A main challenge for such a programme is that explicit flow equations for general actions of the form f (R, R µν , R µνρσ , · · · ) are not available. The main new idea is the observation that for certain subsets of actions, flow equations can nevertheless be derived and analysed. For want of example, we consider the UV fixed point of quantum gravity based on actions of the form f (R, R µν ) = F (R µν R µν ) + R · Z(R µν R µν ) . (1.1) where F and Z are unspecified functions, ultimately determined by the UV fixed point itself. For small curvature, the action reduces to the standard Einstein-Hilbert action. For large curvature, higher order terms modify the action and the dynamics of the theory. The virtue of actions of the form (1.1) is that, when Taylor-expanded in curvature invariants, only a single term arises per canonical mass dimension. In other words, evaluating the renormalisation group equations on spaces with constant curvature is sufficient to find the running for all polynomial couplings of the theory, as well as the non-perturbative renormalisation group flow of the functions F and Z. By and large, these ideas are similar to what had been done previously for f (R) type actions. Generalising those for actions involving more general curvature invariants will enable us to investigate the impact on the gravitational UV fixed point of curvature invariants beyond Ricci scalars. Unlike f (R) theories, our approximations also incorporate the full kinematics of fourth order gravity alongside higher order curvature invariants. As a first application of our results, we investigate the availability of scaling solutions for cosmology such as inflation, and compare findings with f (R)-type theories within the exact same approximation.
The outline of the paper is as follows. In Sec. II, we introduce the main new idea in more detail, and discuss functional renormalisation group techniques and approximations to be used in the sequel. In Sec. III we explain the technical algorithm used to compute the operator traces. Sec. IV introduces our models and the derivation of the RG flow for all couplings. In Sec. V, we search for fixed points within polynomial approximations and beyond. This is followed by an overview of scaling exponents in Sec. VI. We also show that our findings satisfy the boostrap test for asymptotic safety, Sec. VII. In Sec. VIII, we apply our findings to quantum cosmology and establish the existence of de Sitter solutions for small Ricci curvature. In Sec. IX results are compared contrasted with those based on f (R) type actions. A summary of findings and an outlook is given in Sec. X. Technicalities are summarised in three appendices.

II. RENORMALISATION GROUP
In this section we will put forward our main new idea, also recalling the general setup for the renormalisation group for gravity that we shall adopt.

A. Main idea and approximations
Our starting point is the scale-dependent effective action for gravity Γ k [φ,ḡ µν ] where momentum modes down to the wilsonian RG scale k have been integrated out. Here,ḡ µν denotes the background metric and φ the set of dynamical fields. In the case of interest φ includes the metric fluctuation h µν ≡ g µν −ḡ µν as well as ghosts and auxiliary fields arising from the functional measure. The effective average action obeys the flow equation [83,84] k∂ k Γ k [φ,ḡ µν ] = 1 2 STr Γ (2) where Γ (2) k denotes the second variation of the action with respect to the dynamical fields φ and R k is the regulator function which vanishes for large momenta p 2 k 2 while suppressing low momenta p 2 < k 2 in the path integral. The RG scale k therefore splits modes into high and low momentum and the left hand side of (2.1) gives the change in the effective action as k is varied. Evaluation of the super-trace on the right hand side of the equation allows one to obtain the non-perturbative beta functions.
To find an approximate solution to (2.1) we will specify an ansatz for Γ k [φ,ḡ µν ] while making simplifying assumptions in order that the trace may be evaluated in closed form. To this end we write the action as whereΓ k [g µν ] only depends on the metric fluctuation h µν through the full dynamical metric g µν = g µν +h µν and we demand thatΓ k [0,ḡ µν ] = 0. Our first simplification will be to specify thatΓ k [φ,ḡ µν ] is given by the sum of terms in the bare action S which arise from the functional measure namely the gauge fixing action S gf , the ghost action S gh and the auxiliary field action S aux arising from field redefinitions. Since these terms come from the bare action we will neglect their k dependence and therefore they only appear explicitly in the right hand side of the flow equation as part of the hessian operator Γ (2) k . The left hand side of the flow equation will therefore only contain the beta functions of couplings arising from the actionΓ k [g µν ] which is a diffeomorphism functional of the metric tensor g µν . It then suffices to evaluate the flow equation at vanishing field φ = 0 in order to obtain the beta functions for the couplings insideΓ k [g µν ]. This approximation is known as the single field or background field approximation. In general, we may introduce an ansatz forΓ k [g µν ] by writing it asΓ Here, O n,i are invariants of mass dimension 2n built out of the metric field and derivatives thereof, with i differentiating terms of the same dimensionality, and λ n,i the corresponding couplings. As a second approximation, we choose to evaluate the flow equation (2.1) on a spherical background with metricḡ µν obeying the relations whereR,R µν andR µνρσ are the Ricci scalar, Ricci curvature, and Riemann curvature tensors respectively. This choice then limits our ability to distinguish curvature invariants with the same canonical mass dimension. In order to obtain unique flows for the couplings in(2.3), we will therefore include in our approximations for the effective action operators O n such that for every n the corresponding operator is proportional to the nth power of the scalar curvature when evaluated on the sphere Given a choice for the operators O n we can close the approximation by first identifying g µν =ḡ µν and then expanding both sides of the flow equation in powers of the scalar curvatureR = R.
Our main point is that (2.6) is only one out of many possible choices: flows for the running of couplings can be found for any set of O n with unequal canonical mass dimensions. In this light, we put forward a strategy to derive RG flows for much more gravitational actions by replacing (2.6) withΓ where F k and Z k are functions ofX with a, b, c free parameters. By construction, (2.8) is a linear combination of the three independent quadratic curvature invariants with mass dimension four, which are the squares of the scalar curvature R, the Ricci tensor Ric and the Riemann tensor Riem. We also require thatX, when evaluated on spheres with constant scalar curvature r, givesX| spheres = r 2 /4. This normalisation condition implies that only two of the three parameters a, b and c are independent. The family of theories described by (2.7) contain the scale-dependent cosmological constant Λ k and Newton's coupling G k as the leading order coefficients The action involves new curvature invariants which have not been studied previously within nonperturbative RG for gravity. We stress that our choice (2.11) is ad hoc inasmuch as (2.10). The main purpose of this study is to establish whether or not the additional dynamics of the metric field is going to affect the UV behaviour in any significant manner.
In the remainder of this section we derive the generic form of the flow equation in accordance with the above approximations. Many of the results of this section are commonly used in RG studies of gravity and can be found in previous works [20,24,30,31,86]. We generalise and extend these results for the purpose of our study, also noting that they do not depend on the choice for the gravitational actionΓ k [g µν ]. We investigate the RG flow for the model (2.12) in Sect. IV.

B. Field decompositions
In order to calculate the trace we will adopt heat kernel methods. To do so we first bring the second variation into a form where it is a function of Laplacians and scalar curvature only. For this reason we use the transverse traceless decomposition of the original field h µν [87] which was first introduced in the context of the functional renormalisation group in [86]. The various new fields are subject to the constraints, (2.14) Here h =ḡ µν h µν is the trace of the fluctuation, h T µν denotes the transverse-traceless part of h µν , ξ µ is a transverse vector that together with the scalar σ gives the longitudinal-traceless part of h µν according to (2.13). Such a decomposition is advantageous because it leads to partial diagonalisation of the propagator on a spherical background (except between h and σ) and therefore it becomes a simple task to invert it.
Note that after decomposing our original field h µν into its transverse traceless decomposition it does not receive any contributions from the modes of the σ field that obey the conformal Killing equation∇ and similarly from the modes of the ξ µ field that obey the Killing equation These modes, when evaluated on the sphere, correspond to the lowest two modes of σ and the lowest mode of ξ µ respectively and therefore they should be excluded from the trace evaluation. For details of the heat kernel coefficients of the constrained fields as well as for the exclusion of lowest modes see appendix B.

C. Gauge fixing and ghosts
Because of the diffeomorphism invariance of the metric field the effective action has to be supplemented with a gauge fixing term so that only the physical modes of the field are taken into account. It is convenient to choose a gauge fixing term of the form with the gauge fixing condition F µ = 0. Here we will choose F µ to be given by with κ = (32πG N ) −1/2 and h = h µνḡ µν being the trace of the fluctuation. The gauge fixing condition that was originally used in [20] corresponds to δ = d 2 − 1. Here and from now on bared geometrical quantities, such as∇ above, mean that they are constructed by the background metric g µν . Then, by substituting F µ into the gauge fixing action we get For the above gauge fixing, the corresponding ghost action is with the Faddeev-Popov operator given by As with the metric fluctuations, it is convenient to decompose the ghost fields into transverse (C T µ andC T µ ) and longitudinal (η andη) parts In this decomposition, the modes that are unphysical and must be excluded from the trace evaluation are the constant modes of η andη corresponding to lowest modes of the Laplacian. In addition, it has also been argued that the lowest mode of C T andC T and the second lowest modes of η and η should be removed to ensure an exact cancellation with the fluctuations of the field ξ and σ from the gauge fixing term (see below) in the gauge δ = 0, [30][31][32]. Here, we point out an independent motivation for leaving out these additional modes by noting that they corresponds to zero modes of the ghost operator. As such, if retained, these modes would lead to a vanishing Fadeev-Popov determinant, which is why we will leave them out in this paper. We come back to this aspect in Sect. IV B. More technical details about this construction can be found in App. B.
The gauge fixing action (2.19) is already quadratic in the fields. Now we can substitute the metric decomposition (2.13) into (2.19) to express S gf in terms of the metric components as It follows that the contributions to the hessians coming from the gauge fixing action take the form where we have dropped the bars for notational simplicity, since after computing the Hessians we set g µν =ḡ µν and therefore it remains only one metric field. However, it should be kept in mind that all the geometric quantities are constructed with the background metric.
After substituting the ghost decomposition (2.22) into the ghost action (2.20) we have Now we substitute the Faddeev-Popov operator (2.21) into the above equation and we perform the second variation in order to get for the Hessians of each ghost component field where again we have dropped the bars after setting g µν =ḡ µν .

D. Auxiliary fields
The metric decomposition (2.13) is merely a coordinate transformation and as such it induces the Jacobian of the transformation J gr . Here we are going to determine this quantity and we will follow the Faddeev-Popov trick so that we transform the contributions from the determinants into contributions from auxiliary fields. We begin by writing the following relation between the original field h µν and the components of the transverse traceless decomposition which is valid when integrated over d d x √ḡ . Then, at the level of the path integral the Jacobian of the field transformation takes the form with the operators M (0) and M (1T ) coming from the contributions of the scalar σ and the transverse vector ξ µ respectively and having the form The terms containing the transverse-traceless field h T µν and the trace field h do not contribute to the transformation Jacobian since they do not involve operators and under the path integral they are simple Gaussian integrals contributing only a constant. Now, we would like to express these determinants as new contributions to the action in terms of auxiliary fields. For this reason, we follow the Faddeev-Popov trick and write them as Gaussian integrals. We start with the determinant of the scalar field where λ andλ are complex Grassmann fields coming from the denominator of the above expression and ω is a real field coming from the numerator. Thus the action for the scalar auxiliary fields reads S aux (0) = d d x √ḡ λ M 0 λ + 1 2 ωM 0 ω and the Hessians forλ, λ and ω are given, after dropping the bars, by Similarly for the determinant detM (1T ) can be written in terms of a path integral with the aux- ζ T µ coming from the contribution of the transverse vector with c T µ andc T µ being complex Grassmann transverse vector fields coming from the denominator of the above expression and ζ T µ a real transverse vector field coming from the numerator. The corresponding Hessians, after dropping the bars, are In the same way that the metric decomposition (2.13) induces the Jacobian of the transformation we get contributions J gh from the decomposition of the ghost fields (2.22). Now, the original ghost fields C µ andC ν obey the following identity with the components C T µ ,C T µ , η andη The fields C T µ andC T µ do not involve differential operators and thus they contribute only a constant. Now, at the level of the path integral the Jacobian of the transformation comes only from the scalars η andη and after performing the Gaussian integral it can be written as As before we rewrite this expression in terms of the contribution to the action of auxiliary fields, so that it takes the form where now the fieldss and s are complex conjugate scalars with the corresponding auxiliary field action S aux (gh) = − d d x √ḡs −¯ s. The Hessian for this auxiliary field, after dropping the bars, is given by Since we exclude the two lowest modes of the ghostsη and η we must also remove these modes from the spectrum ofs and s. This completes the introduction of the basic degrees of freedom and auxiliary fields.

E. Gravitational path integral
We now turn to the introduction of a Wilsonian momentum cutoff and a path integral formulation of the theory. After gauge fixing and employing the field decompositions we can collect all propagating fields into a "superfield" ϕ with components where we indicate the origin of each field. We also recall that C T µ ,C T ν η,η, c T µ ,c T ν ,λ and λ are Grassmann fields. In the conventions put forward here, the (gravitational) partition function can then be written as a functional integral over the set of fields (2.43), with J denoting external currents and the fields (2.43) the integration variables. The term S in (2.44) denotes the fundamental action for all fields of the theory. The term ∆S k introduces a cutoff at momentum scale k, in the sense of Wilson. Its task is to regularise the path integral (2.44) and to suppress the propagation of momentum modes with momenta q smaller than the RG scale k. Technically, this is achieved by taking ∆S k to be quadratic in the fluctuation fields, The cutoff functions R ij k , which take the role of momentum-dependent mass terms, are at our disposal. We demand that ∆S k → 0 for k → 0, to ensure that the partition function reduces to that of the full physical theory in the limit where the momentum cutoff is removed. In the UV limit k → ∞, we also require ∆S k to grow large to ensure that fluctuations are switched off. Coming back to our settings in the context of gravity, we have that The auxiliary field action denotes the sum S aux = S aux (0) + S aux (1) + S aux (gh) , which, together with the gauge fixing S gf and ghost action S gh , has been defined in the previous subsection. In turn, the gravitational action S grav has been kept arbitrary for now. It will be specified in Sect. IV for the purposes of the present paper. Similarily, the cutoff action (2.45) takes the form Notice that the momentum cutoff is "diagonal" in field space, (2.43 [17,88], where φ = ϕ J stands for the expectation value of the quantum field in the presence of the source J. 1 The renormalisation group scale-dependence of Γ k , (2.1), is then obtained from the exact flow of the partition function ∂ t Z k = − ∂ t ∆S k J by standard manipulations [20,83,84]. With the above ingrediences at hand, and following textbook procedures, we are led to a path integral representation for the scale-dependent gravitational effective action (2.48) Notice that in (2.48), the cutoff term (2.45) is now evaluated around the shifted fields ϕ → ϕ − φ. Technically, the representation (2.48) takes the form of an integro-differential equation for the quantum effective action of gravity. Moreover, it is compatible with the flow equation (2.1). Below, we will primarily be concerned with identifying gravitational effective actions which display interacting UV fixed points Γ * , and, as such, qualify as fundamental candidates for a quantum theory of gravity [17].

F. Wilsonian momentum cutoff
The last ingredient of the flow equation (2.1) consists in the choice of the regulator term R k . In general, for each component field, the second variation of the effective action 1 √ḡ δ 2 Γ k /δφ i δφ j , which we will often denote as Γ (2)φ i φ j k , will be a function of some differential operator ∆. The regularised inverse propagator is then given by (2.49) The general prescription we are going to adapt is that the Wilsonian regulator term R φ i φ j k for the inverse propagator (2.49) should be chosen in such a manner that it leads to the replacement for all component fields (for explicit examples, see Sect. III B). Here, R k (∆) denotes the profile function of the Wilsonian momentum cutoff. Hence, at this stage we are still left with some freedom for choosing the cutoff scheme by taking different definitions for the differential operators ∆, and different choices for the profile function R k . For the former, one may write the differential operator as where∇ 2 represents the covariant derivative for both diffeomorphism invariance and all other possible gauge symmetries, and E denotes a linear map acting on the field. Following the classification in [32] a type I cutoff is defined as E = 0 such that the Wilsonian regulator will only depend on −∇ 2 and not on any endomorphism E. Type II regulators include a non-vanishing endomor-phism E which is independent of the RG scale k. Finally, type III regulators are those where the endomorphisms E additionally depend on the RG momentum scale.
To proceed we also have to make a choice for the profile function. For our calculation we are going to choose an optimised cutoff [89][90][91][92] where Z denotes a suitably chosen differential operator in real space. Most of the time we will be using Z = ∆ (type I cutoff). In momentum space, Z is taken to be momentum squared q 2 . A key feature of the optimised cutoff is that it vanishes identically for momenta larger than the RG cutoff scale q 2 > k 2 . For smaller momenta, it acts like a momentum-dependent mass term. Its scale derivative equally vanishes for momenta larger than the cutoff scale, and is given by The second term is absent provided that Z carries no scale-dependence of its own, that is for type I and type II momentum cutoffs. It has been shown that this regulator leads to particularly good stability and convergence properties of the functional flows [89,91]. It also provides algebraic simplifications for the flows which will become of great importance when we will be concerned with gravity. Another key benefit of the optimised cutoff in gravity is that the heat kernel expansions truncates at a finite order, allowing for closed expressions for the flow equations (see Sect. III). For a more detailed analysis of the optimised regulator benefits see [93,94]. Since its introduction, the regulator (2.52) has been very popular in the context of quantum gravity and it was used in the previous studies of f (R) gravity [30,31], too.

III. TRACE COMPUTATION ALGORITHM
In this section we give a general algorithm for computing functional traces of a general operators which can be written as a Taylor series in terms of − = −∇ µ∇ µ up to a maximal power p, no longer contain derivative operators. The value of p can take any integer value 1 ≤ p ≤ ∞. The case p = ∞ can arise if we consider, for example, terms in the action of the form d d x √ gR e /M 2 R for some mass scale M . Here we shall consider type I regulators explicitly however the results are easily generalised to the case of type II and III regulators. In particular for a more general regulator one replaces − → ∆ = − + E which in turn leads to different values for the coefficients A

A. Heat kernels
For the computation of the functional trace in (2.1) we use the heat kernel techniques [95,96] which we will recall in this section. The RHS of the flow equation (2.1) consists of the functional trace (i.e. the sum over all indices and the integration over all momenta) over the quantity (Γ In general, the functional trace of a functional W (∆) of an operator ∆ is given by the sum where λ i are the eigenvalues of the operator. By introducing the Laplace anti-transform W (∆) = ∞ 0 dt e −t∆W (t) we can write (3.2) as where K(t) = e −t∆ is the heat kernel of the operator ∆. In equation (3.3) the subscript s denotes the spin of the field which the operator ∆ acts on and takes the values 0 for scalars, 1 for vectors and 2 for second rank tensors. The trace of the heat kernel has a well known early time (t → 0) asymptotic expansion given by [95] Tr (3.4) The coefficients b n are called the heat kernel coefficients and are sums of curvature invariants and their derivatives. The traces of the coefficients that we will need here evaluated on the spherical background are given in appendix B. Now we define Now the functional trace is expressed as the early time expansion for the heat kernel of the operator ∆ instead of the spectral sum (3.2). This becomes apparent when we make use of the Mellin transform and we express the functional Q n [W ] in terms of the original function W (z). Then for every n > 0 and similarly for Q −n with n ≥ 0 ∈ Z In the following, we develop a general algorithm for computing these functionals just with the knowledge of the second variation. It will then become clear that the usage of the optimised cutoff (2.52) truncates the series and we end up with a closed expression for the flow equation which is one of the great benefits of the optimised cutoff. Therefore with a finite number of heat kernel coefficients we can derive the full flow. That the series truncates beyond some order becomes apparent from the last formula. As we will see later, with the optimised cutoff, W will be a polynomial function of z with the highest power determined by the number of derivative operators in the second variation of the action.
As will become clear later, there are some cases where we have to exclude certain modes from the trace computation due to constraints that some fields obey. Here we are going to state the general mechanism of how this is taken into account. More details about the specific exclusions that we are going to use can be found in App. B. Since the trace of an arbitrary smooth functional W (O) can be represented as its spectral sum, the general rule for omitting the missing modes from the trace of an operator valued function is where the m primes on the LHS indicate the first m modes to be subtracted, l labels the eigenvalues Λ l (d, s) in the spectrum of the operator −∇ 2 , D l (d, s) gives the multiplicity of each eigenvalue, d is the spacetime dimension and s = 0, 1, 2 for scalars, vectors and tensor respectively. For the operator −∇ 2 acting on scalars, transverse vectors and transverse-traceless symmetric tensors D l (d, s) and Λ l (d, s) are given in Table 8 of App. B.
In what follows we present an algorithm for computing the trace having as input the general form of the second variation. We then determine the form that the regulator term should have, we move on by implementing the optimised cutoff and finally we evaluate the integrals for the functionals Q n [W ] defined in (3.6) and (3.7). We also evaluate these functionals for the case that we have off-diagonal terms in the second variation with the mixing of two components. Finally, we compute the two lowest excluded modes for a scalar field and the lowest excluded mode for a vector field.

B. Wilsonian momentum cutoff
We will now discuss the specifics of the regulator term R φ i φ j k , which needs to be added to (3.1). For convenience, we introduce the short-hand notation to denote the φ i φ j -component of the inverse regularised propagator. We will drop the φ i φ j indices when it is clear from the context. The coefficients A m (3.1) depend on the momentum scale k through the couplings of the theory. For our purposes it is enough to determine the regulator for type I cutoff. However, the process described here has a straightforward generalisation to the other types of cutoffs described in [32]. For type I cutoffs we define the Wilsonian regulator (see Sec. II F) by demanding that the addition of the regulator term R φ i φ j k to the inverse propagator of the corresponding field leads to the replacement of the operator − by − + R k (− ), where R k (− ) is a suitably chosen profile function characterising the Wilsonian momentum cutoff for the field modes. Using (3.9) together with (3.1) and the requirement on the regulator, the inverse regularised propagator takes the form Solving (3.10) with (3.9) for R k , we find It is readily confirmed that R k → 0 entails R k → 0, as it must, alongside R k → ∞ for R k → ∞. Now that we have an explicit form for the regulator term we can proceed to determine its scale derivative, which is an essential component of the flow equation (2.1). For type I cutoff the operator − does not contain any couplings and thus it does not carry any RG scale-dependence. Therefore the ∂ t derivative of the regulator reads (3.14) In this paper, we will use one and the same shape function R k for all fields. This choice can be omitted in more general settings.

C. Diagonal pieces
Now we have all the ingredients that we need in order to evaluate the trace of the flow equation. For this purpose we split the RG flow into two parts according to (3.12). Using (3.13), (3.14) the equation for the diagonal part reads where G is defined in (3.10) in terms of the coefficients A m and the profile function R k . With the optimised cutoff (2.52) and its scale derivative (2.53) we can evaluate the functionals (3.6) and (3.7) for the above expressions. We note that for a type III cutoff there will be an extra term in (2.53) proportional to ∂ t E. One observes that both the above integrals in (3.15) will have as a common factor the step function θ(k 2 − z) and therefore we can change the integral limits from ∞ 0 to and replace the step functions by 1 everywhere. After substituting − = z, we then have for the coefficients C 1 , C 2 and G in (3.15) the following results, (3.16) Both G and C 1 have become independent of z for the optimised cutoff. Now we turn our attention to the evaluation of the integrals Q n [W ] for positive n, (3.6). The integrand W (z) in (3.6) can be split in two parts according to (3.15). Defining W 1 = (C 1 ∂ t R k )/(2G) and W 2 = C 2 /(2G) we evaluate the corresponding integrals I i n = Q n [W i ] separately. Also these two integrals contain step functions which are treated as above. Then we have (3.18) By summing the two contributions (3.17) and (3.18) we get an expression for the coefficients Q n of the heat kernel expansion for positive n, For negative integer n, the z-dependence of W (z) in the Q-functionals (3.7) comes only from the term. It is then enough to compute the derivative of this term. In the previous sections we have disregarded the explicit form of the θ(k 2 −z)-functions, whose effect, as overall factors, amounts to a change of integration limits. Here, we use the properties of the step function as introduced above, and additionally exploit its idempotence, and that its derivative dθ( will not contribute at z = 0 and k 2 = 0. Altogether, we find After evaluating this expression at z = 0, the only term that survives is the one with m = n. Therefore, we conclude that the coefficients Q −n are given by The results (3.19) and (3.21) will be needed for the subsequent derivation of the RG flows of couplings.

D. Non-diagonal pieces
When there are off-diagonal components in the second variation (3.1), the inversion of (3.1) and (3.9) become more difficult. For the case where a mixing arises between two field components φ 1 and φ 2 , the inverse of the corresponding submatrix G φ i φ j involving the fields φ 1 and φ 2 is given by Performing the matrix multiplication in the field space φ 1 and φ 2 , the trace takes the form The expressions for G φ i φ j (3.10), for the regulator (3.11) and the coefficients of the expansion (3.12) remain the same as given previously, except that we restore the indices φ i φ j in all explicit expressions involving the momentum cutoff. For the remaining step we assume that the second variation G φ i φ j is symmetric. Adopting the optimised momentum cutoff (2.52), G becomes z-independent and no longer contributes to integration over z, see (3.16). The only remaining z−integration is the one over ∂ t R k , which is known from the previous section. Putting everything together, we find that the coefficients Q n for the non-diagonal piece for positive and negative index are given by respectively. Here, we have used the expressions and G φ i φ j is given as in (3.16).

E. Excluded modes
As explained in Sec. III A we often encounter the trace of a function where some eigenmodes of the operator − have to be excluded. For example we saw in Sec. II B that after performing the decomposition of the fluctuation h µν into its components through the transverse traceless decomposition, the lowest mode of ξ µ and the two lowest modes of σ have to be excluded. In order to incorporate this fact into the evaluation of the traces we make use of (3.8). For our calculation we will need the form of Tr (0) , Tr (0) and Tr (1T ) . We start with the case where the lowest mode of a scalar has to be excluded. According to Tab. 8 in App. B, for the lowest mode of a scalar field we have The part W (0) of the above equation which has to be excluded from the scalar trace is simply evaluated by setting z = 0 in the expression for W (z). Then we denote the lowest excluded mode of a scalar as X (0) and we have For the exclusion of the two lowest modes of a scalar field we have according to Tab. 8 In general, the result will involve theta functions coming from the choice of the profile function and more precisely from (2.52) and (2.53). For the specific expression under consideration the theta functions are of the form θ(k 2 − R d−1 ). Since in the following we are going to focus on an expansion in small R k 2 we can evaluate these theta functions in the approximation R k 2 1 in order to get for the exclusion of the two lowest scalar modes Finally we have to determine the exclusion for the lowest mode of a transverse vector. Again, by reading off the multiplicities and the eigenvalues from Tab. 8, we find As before we evaluate the theta functions coming from the profile function for the approximation R k 2 1 and we have for the lowest exclusion mode of the vector This completes the collection of the technical tools that we need to calculate the renormalisation group flow for certain general classes of gravitational actions projected on the sphere. In the next section we will apply this formalism for a concrete model of quantum gravity involving Ricci tensor invariants.

IV. RICCI TENSOR EXPANSION
In this Section, we derive RG flows for theories of gravity whose actions depend on certain Ricci tensor invariants, generalising the f (R) type actions towards more complicated curvature invariants. This will give us new information on gravitational renormalisation group flows complementary to existing results in the f (R)-approximation. A new qualitative feature is that now the propagator of the tensor modes will also contain higher-derivative contributions in contrast to the f (R)-approximation.

A. Hessian
Here we are going to use the methods developed in the previous Section to derive the Hessians for all the fields that contribute to our ansatz. As explained in Sect. II A the effective average action takes the form φ i φ j for these fields. The missing element is the computation of the second variation for the gravitational partΓ we are yet to fixΓ k [g], which we termed S grav in (2.46). To that end, we assume that the gravitational effective average action is formed of two unspecified functions of the square of the Ricci tensor F k (R µν R µν ) and Z k (R µν R µν ), with the latter one additionally multiplied by the Ricci scalar Once the equation is evaluated on a spherical background and expanded in power of the curvature we then have exactly one invariant O n ∝ R n (see (2.5)) for all integer values of n. The ansatz (4.2) provides a natural extension to preceding work on the effective average action in gravity including more complicated tensor structures with non-trivial dynamics. In order to compute the Hessians and to evaluate the flow equation we need to know Γ (2) k . To proceed we expand the gravitational actionΓ k [g] as to extract the part quadratic in h µν which takes the form where all the geometric quantities of the above equation are constructed from the background metric with δg µν = h µν . From now on we drop the bar notation since we identify the metric with the background. Moreover it is understood from now on that F k → F k (R µν R µν ) as well as Z k → Z k (R µν R µν ) and that the primes denote derivatives with respect to the argument. Using the expressions in App. A we find that the quadratic part takes the form We then use the metric field decomposition as defined in Sect. II B to find the Hessians in terms of each component field. These are given in App. A. In Tab. 1 we summarise the contribution from each individual component field after adding the contributions from the gauge fixing part, the ghost part and the auxiliary fields.  Table 2. Summary of second variations (continued) of the gravitational action, suitably decomposed into independent component fields.

B. Gravitational flow equation
Having evaluated the second variation for our ansatz we are ready to compute the flow equation using the machinery developed in Sect. III. We need to introduce dimensionless variables for the two functions that make up the effective average action. Note that after taking the second variation and before computing the trace we evaluate all the expressions on the spherical background metric. Thus the arguments of the functions F k and Z k become R 2 d and we will mostly drop the arguments. Then we define Here, and in order to avoid a proliferation of symbols, we denote the dimensionful Ricci scalar as R, and the dimensionless Ricci scalar curvature as R, with Hence, from now on, R denotes the Ricci scalar in units of the RG scale k. For the RG scale derivatives of the functions F k and Z k we have For the dimensionless functions f and z the arguments become R 2 d . Moreover, the notation f (n) denotes the nth derivative of f with respect to its argument. Gathering together all the above we can compute the LHS of the flow equation for the effective average action given by (4.1). Then we have The RHS of the equation is given by the sum of the traces for the individual components after using the algorithm described in Sec. III and subtracting the excluded modes, as indicated. Then, in terms of the components we find ss + Rs s (4.12) where the Hessians for each field component are given in Tab. 1 and 2. The primes at the traces denote the number of lowest modes to be excluded as described in Sect. III and the traces are to be computed using the algorithm for the Q-functionals and the relevant b n coefficients listed in App. B. Moreover, the auxiliary fields have inherited the primes from the fields from which they originate. We note that the prescription to additionally leave out the lowest two modes from the ghosts as detailed in Sect. II C has lead to an additional prime on each of the three traces in the last line of (4.12).
In order to simplify the computation of the traces we fix the gauge and focus on a specific spacetime dimension d = 4. For the gauge parameters, we choose This choice of gauge results in three simplifications of the flow equation. Firstly, since we take the Landau gauge limit α → 0 the gauge fixing terms diverge. This entails that the Landau gauge is a fixed point under the renormalisation group flow [97]. Secondly, since terms proportional to 1 α are also included in the regulator, only the terms proportional to 1 α survive in the limit α → 0. As a result the non-diagonal term σh vanishes since it has no dependence on α (for δ = 0) while the denominator which involves the components hh and σσ becomes very large. Finally, the limit (4.13) entails that the physical and the gauge degrees of freedom totally decouple. The gravitational degrees of freedom are encoded in h T h T and hh, while the gauge degrees of freedom are contained in ξξ and σσ. Together with several cancellations which occur in the gauge (4.13) the flow equation (4.12) simplifies substantially and takes the final form (4.14) After performing the trace computation on the RHS of (4.14) and combining with the corresponding LHS given in (4.11), we obtain the explicit flow for the two functions f and z introduced in (4.6) and (4.7) as Here, the expression on the RHS encodes the contributions from fluctuations and it splits in several parts as for any a = 0. For this reason, the RHS becomes negligible in the limit 1/f → 0 and 1/z → 0, corresponding to the classical limit where quantum fluctuations are absent. Thirdly, the scaling exponents of the theory are sensitive to small deviations form the fixed point. To determine the scaling exponents, we also need to know the functions I i , i = 1, · · · 5, in addition to the fixed point solution itself. The fixed point and the scaling exponents are determined in Sect. V and Sect. VI, respectively.

V. INTERACTING FIXED POINTS
In this section, we turn to a detailed analysis of fixed points of the gravitational flow (4.15) and (4.16) for gravitational actions of the form (4.1), (4.2), particularly in view of the asymptotic safety conjecture. Interacting ultraviolet fixed points of the gravitational effective action Γ * allow for a fundamental definition of quantum gravity, (2.46). We discuss the classical and quantum fixed points, using polynomial expansions and numerical integration techniques.

A. Classical fixed points
We begin by discussing the classical fixed points in the absence of fluctuations, following [57]. In particular, the RG flow (4.15) simplifies and becomes Most notably, we observe that while the functions f and z mix on the quantum level, the mixing is absent in the classical limit. Consequently, the RG flow (5.1) can be solved analytically by introducing the new functionf in terms of which (5.1) becomes Alternatively, we can write down a flow for the inverse, Solutions and fixed points of (5.3) and of (5.4) have been analysed in [57]. There, it has been shown that theories described by (5.3) develop a Gaussian fixed pointf * = 0 and an infinite Gaussian fixed point 1/f * = 0. The most general solution of (5.3) is given bȳ for arbitrary function H(x) which is solely determined by the boundary conditions at some reference energy scale t = 0.

B. Quantum fixed points
Next we turn our attention to quantum fixed points in the presence of fluctuations. At an interacting fixed point of the theory we observe ∂ t f = 0 = ∂ t z and the RG flow reduces to We then proceed to expand (5.6) in powers of the dimensionless curvature R. To this end we approximate the functions f and z polynomially as where the couplings λ n are dimensionless. We recall that the argument X of the functions f = f (X) and z = z(X) is given by X = R 2 /4. It arises from evaluating the square of the Ricci tensor R µν R µν on spheres with constant scalar curvature R = R k 2 , all in units of the RG scale k. The overall normalisation of the functions f and z is chosen such that the couplings λ 0 and λ 1 relate to the dimensionless Newton coupling g = G k k 2 and the dimensionless cosmological constant λ = Λ k /k 2 as λ = −λ 0 /(2λ 1 ) (5.9) g = −1/λ 1 .

C. Fixed point data
We solve the fixed point condition β i = 0 using (5.12) to identify fixed point candidates at each order of the approximation from N = 2 to N = 21 [85,98]. The analoguous problem in f (R) type theories of gravity allowed for an exact recursive solution of all fixed point couplings in terms of the two lowest ones [19,57].
More generally, exact recursive solutions are possible provided that the RG flows β i depend only linearly on the coupling λ n with the highest index n. In the present case, this is not the case as the dependence on the highest coupling comes out quadratic. This new feature arises owing to the fluctuations of the Ricci tensor as opposed to those of the Ricci scalar. For this reason, we adopt  Table 3. Summary of the fixed point couplings λ n (N ) for n = 0, ..., 11 and for selected approximation orders N . Overall, we observe a fast convergence. Of all couplings, only λ 8 (N ) (red-shaded) starts off with the wrong sign at first appearance (N = 9). a different strategy to solve for fixed points. We begin with the fixed point in the Einstein-Hilbert approximation (N = 2), where a unique interacting fixed point is found by solving β 0 = 0 = β 1 for (λ * 0 , λ * 1 ). As boundary conditions, we impose β i = 0 and λ i = 0 for all i > 1. Provided a fixed point is found, we then increase the approximation order and use the previously-found values (λ * 0 , λ * 1 ) as starting values to search for interacting fixed points in the couplings (λ * 0 , λ * 1 , λ * 2 ) in its vicinity. This is repeated from N = 2 up to N = 21. We find that the procedure converges rapidly. At each and every order, we observe an UV fixed point in accordance with the asymptotic safety conjecture. Starting from order N = 3 we also find spurious fixed points which either disappear in the next order or show inconsistencies in the number of relevant eigenvalues. These spurious fixed points are discarded.
In Fig. 1, we plot the decadic logarithm of the fixed point condition in field space, (5.6), for all approximation orders N . With increasing N , we observe that the accuracy of the fixed point solution and the domain of validity increases for all scalar curvatures within the radius of convergence R c . The latter can be estimated from Fig. 1 as the absolute value of the scalar curvature beyond which the accuracy in the fixed point solution starts decreasing with increasing N . We find The estimate (5.13) agrees with the location of the first singularity along the real axis of the fixed point equation, discussed in (5.22) below. As such, our findings indicate that the polynomial approximation exhausts the maximum radius of convergence set by a pole in the real field axis of the RG flow itself. Poles could have arisen in the complex field plane [68,[99][100][101][102][103], but they do not in the case at hand.
Our numerical results for the fixed points are shown in Tabs. 3 and 4, and in Fig. 2. From the data in Tab. 3 we notice that the convergence of the first few couplings is fast. Also, all couplings arise with the correct sign already at their first appearance, except λ 8 , see Tab. 3. At order N = 3, the inclusion of the dimensionless coupling λ 2 proportional to R µν R µν is expected to have a strong effect. It does lead to a significant shift on the couplings λ 0 and λ 1 , with λ 1 changing by about 25%. From approximation order N = 4 onwards, however, the inclusion of higher dimensional operators rapidly stabilises the fixed point values. We estimate the asymptotic values of the first  few couplings by taking the average over the four best approximation orders to find λ 0 = 0.3250067185 ±2 · 10 −10 λ 1 = −0.8306876250 ±3 · 10 −10 λ 2 = 0.1915686687 ±2 · 10 −10 λ 3 = −0.6291586052 ±9 · 10 −10 λ 4 = −0.229072717 ±6 · 10 −9 λ 5 = −0.00337078 ±2 · 10 −8 λ 6 = 0.11247802 ±8 · 10 −8 λ 7 = 0.0061012 ±2 · 10 −7 . (5.14) The indicated error corresponds to a standard deviation in the data from the mentioned average. The fast convergence can also be read off from Fig. 2, showing that the accuracy in the couplings increase by roughly a decimal place for approximation order N → N + 2. This result is quite the opposite of what has been observed previously in f (R) models of gravity. There, the coupling proportional to the square of the Ricci scalar leads to a strong impact on the fixed point coordinate. Also, the R 2 coupling itself showed a slow convergence with increasing approximation order. A brief remark on accuracy and error estimates is in order. The error estimates in (5.14) refer to the accuracy with which the relevant fixed point equations have been solved. In principle, results at low polynomial orders cannot be taken for granted due to the non-perturbative nature of the graviton [19]. However, the high numerical accuracy and the fast convergence make sure that the fixed point is a reliable approximation of the full, non-perturbative, fixed point of (4.15), (4.16). On the other hand, our error estimate in (5.14) does not account for proper systematic errors. For strategies to estimate the latter using functional renormalisation we refer to [94].  Figure 4. Shown is the non-perturbative "fixed point functional" f * (R 2 /4) + R z * (R 2 /4) given in (5.7), (5.8) as a function of the dimensionless scalar Ricci curvature R and the polynomial approximations from N = 3 up to N = 19 (in steps of 2, thin black lines). The highest approximation order (thick red line) also coincides with the numerical integration of (5.18), (5.19). Notice that the radius of convergence, indicated by the shaded area, is maximal, and given by R p in (5.22).

D. Universal coupling ratios
Fixed point couplings are non-universal. Still, some universal quantities of interest are given by specific products of couplings which remain invariant under global re-scalings of the metric field g µν → g µν . Note that the classically dimensionless coupling λ 2 remains invariant under the rescaling (5.15). All other couplings scale non-trivially. Consequently, various products of couplings can be formed which stay invariant under (5.15). Such invariants may serve as a measure for the relative strength of the gravitational interactions [104]. For couplings including up to λ 3 RR µν R µν we may construct three independent invariants with values λ 0 /λ 2 1 = 0.4709956080 ± 5 · 10 −10 λ 0 λ 2 3 = 0.1286508384 ± 4 · 10 −10 λ 1 λ 3 = 0.5226342675 ± 6 · 10 −10 .

(5.17)
Tab. 4 show the convergence of couplings, universal coupling ratios, and the leading scaling exponents. We postpone a detailed discussion of universality until Sect. VI.

E. Beyond polynomial expansions
In order to go beyond the polynomial approximations we can integrate the fixed point condition numerically. In a first step, we consider even and odd parts of the fixed point condition in R (5.6), leading to These two equations can then be taken as differential equations for the functions f (x) and z(x) which depend on R through the variable x = R 2 4 . Taking the initial conditions from our polynomial solutions we can solve these two equations to give f and z. Since the equations depend on up to the third derivatives of these functions we must give three initial conditions for both f and z.
Resolving them for the highest derivatives, the coupled differential equations take the form where For the numerical integration, we therefore will have to chose initial conditions away from (5.22). Initial conditions are provided through the polynomial fixed point in its domain of validity. This limits the integration of the fixed point equations (5.18), (5.19) to the regime 0 < x < 1.00649. We fix initial conditions at a sufficently small value for x, say x = 1/100. We have checked that results do not depend on this choice. In Fig. 4 we plot the function f + Rz as a function of R and compare it to the polynomial approximations from orders N = 3 to N = 19. We observe that the polynomial solutions are in good agreement with the numerical solution all the way to R = ±2.00648. The existence of the formal singularity at R p = 2.00648 and the fact that our polynomial approximation agrees with the numerical solution up to this point indicates that the radius of convergence for the polynomial approximation is maximal: there are no other convergence-limiting singularities in the complex plane of scalar Ricci curvature with |R sing. | < R p .

VI. UNIVERSALITY AND SCALING EXPONENTS
In this section, we turn to the computation of universal scaling exponents {ϑ n } characterising the UV fixed point.

A. Scaling exponents
Unlike fixed point values which are non-universal, scaling exponents are universal and, in principle, mesurable in an experiment. They are deduced as the eigenvalues of the stability matrix at the fixed point M ij = ∂β j /(∂λ j )| * , which in turn is obtained from (5.12). At each order N in the approximation we retain N operators in the effective action. We therefore find a set of N eigenvalues {ϑ n (N ), n = 0, · · · , N − 1} , (6.1) which we order according to the magnitude of their real parts, Re ϑ n ≤ Re ϑ n+1 . Negative (positive) eigenvalues correspond to relevant (irrelevant) operators in the UV. It is expected that a UV fixed point displays, at best, a finite number of relevant operators in order to ensure predictivity of the theory. At each order of the approximation from N = 2 to N = 21 we find that there are always three UV attractive directions, corresponding to three negative eigenvalues, exactly as was the case for the f (R) approximation up to order N = 35 [19,57]. Some scaling exponents are complex conjugate pairs. It is then convenient to denote their real and imaginary parts of the eigenvalues (6.1) as θ = −Re ϑ and θ = ±Im ϑ. Specifically, the three leading relevant exponents {ϑ 0 , ϑ 1 , ϑ 2 } are a complex conjugate pair together with a real eigenvalue, which we write as ϑ 0 = −θ 0 + iθ 0 , ϑ 1 = −θ 0 − iθ 0 and ϑ 2 = −θ 2 . All irrelevant eigenvalues ϑ n for n ≥ 3 settle as complex conjugate pairs, which we write pairwise as ϑ n = −θ n +iθ n and ϑ n+1 = −θ n −iθ n (n ≥ 3 and odd). Using these conventions, in Fig. 5 we show the real part of the eigenvalues ϑ n as functions of the polynomial approximation order. In Fig. 6, we show the fast convergence of {θ 0 , θ 2 , θ 3 , θ 5 , · · · , θ 15 }, again as functions of the polynomial approximation order.
The values of the critical exponents can be found in Tab are quite stable from order to order, except for the first occurrence of the perturbatively marginal coupling λ 2 with eigenvalue θ 2 at order N = 3. The reason for this has been given in [57]: the invariant √ gR µν R µν is classically marginal meaning that the RG flow for the corresponding coupling does not have a term linear in the coupling. Consequently, higher order interactions are needed to help the classically marginal coupling to develop a fixed point of its own. In other words, the operators with irrelevant scaling dimension are necessary to help the relevant and marginal invariants to develop stable fixed points. This structural pattern is generic for fixed points of quantum fields theories within polynomial or vertex expansions.
inclusion of more complicated tensor structures.

B. Degeneracy
We briefly comment on the complex eigenvalues. Complex conjugate pairs of eigenvalues are a consequence of interactions and indicate a degeneracy amongst the scaling of eigenoperators. Technically, they arise due to large off-diagonal entires in the stability matrix. The degeneracy indicates in that two linearly independent eigenperturbations decay with the exact same power law (associated to the real part of the eigenvalue), the sole difference being a relative phase (associated to the imaginary part of the eigenvalue) [57]. If complex conjugate exponents persist in the full physical theory, it would be important to fully understand the dynamical origin. However, the degeneracy may very well be an artifact of our inability to retain all gravitational couplings. If so, we expect that the degeneracy is lifted, or becoming less pronounced, provided additional gravitational couplings are retained in the effective action beyond those studied here. In the present case, we find that imaginary parts (if they arise) are small. In particular, we find that the ratio of real to imaginary part Imϑ n /Reϑ n and the angles φ n = arctan Im ϑ n Re ϑ n , (6.4) if nonzero, becomes increasingly small with increasing n. Our results for the angles (6.4) are shown in Fig. 7. In view of the discussion given above, small imaginary parts imply that the degeneracy amongst pairs of eigenperturbations is mild, and increasingly milder for increasingly higher order invariants. We conclude that our findings are consistent with the picture given above. x Figure 8. The bootstrap test for asymptotic safety: from left to right, each line D i shows the i th largest eigenvalues from all approximation orders N ≥ i, connected by a line, and x(N ) = N + 1 + i. At fixed i, and with increasing approximation order N corresponding to increasing x, we observe that all curves D i consistently grow with N . This establishes that invariants with a larger canonical mass dimension lead systematically to larger scaling exponents (see main text).

C. Gap
We now turn to the gap in the eigenvalue spectrum, following [57]. Here, our results establish three relevant eigendirections corresponding primarily to the cosmological constant √ g, the Ricci scalar √ gR and the square of the Ricci tensor √ gR µν R µν . The smallest negative eigenvalue is of the order ≈ −1. In turn, the most relevant of the irrelevant eigendirections, primarily induced by the invariant √ gRR µν R µν , has the eigenvalue ≈ −5.00. Classically, its Gaussian eigenvalue reads −2.
The large numerical value means that fluctuations have made the operator less relevant. Denoting the difference between the smallest relevant eigenvalue and the smallest irrelevant eigenvalue (or their real parts if complex) as the "gap" ∆, we then find from (6.2) that the gap in the eigenvalue spectrum is given by Notice that the gap is roughly of the same size as the first irrelevant eigenvalue. Furthermore, the gap is only slightly larger than the gap ∆ ≈ 5.5 observed in the f (R) approximation [19,57]. This comparison seems to suggest that the invariant √ g R R µν R µν is a less relevant operator than √ gR 3 .

VII. BOOTSTRAP AND NEAR-GAUSSIANITY
In this section, we apply the bootstrap test for asymptotic safety [19] and investigate the large order behaviour of scaling exponents.

A. Bootstrap hypothesis
A key challenge in the reliable determination of asymptoticaly safe UV fixed points consists in the fact that the set of relevant, marginal, and irrelevant operators at the fixed point are not known beforehand. This implies that the set of universal eigenvalues {ϑ n } (7.1) at the interacting fixed point is not known beforehand. This is in stark contrast to theories with asymptotic freedom, where the set of positive, vanishing, or negative eigenvalues is known a priori to coincide with the Gaussian exponents. Furthermore, fixed point searches in interacting theories generically require approximations, and even if stable fixed points are established, one still has to acertain that invariants neglected thus far are not spoiling the result. In contrast, in asymptotically free theories the finite set of relevant and marginal operators is known beforehand allowing for systematic approximations. For this reason, in [19], a bootstrap strategy has been put forward to help circumnavigate the lack of a priori information about interacting fixed points, applicable for UV and IR fixed points alike. As our working hypothesis we will assume that canonical power counting continues to be a good guiding principle at an interacting fixed point, or, in other words, we assume that • the relevancy of invariants at an interacting fixed point continues to be governed by the invariant s canonical mass dimension . 2) The motivation for this hypothesis is that all known interacting fixed points in quantum field theory share this property. The value of a working hypothesis consists in the fact that it can be tested and falsified, systematically, within given approximations. If successful, it establishes the existence of (asymtpotically safe) fixed points under the hypothesis that canonical power counting is applicable.

B. Testing asymptotic safety
To establish that our results conform with the bootstrap hypothesis, we display in Fig. 8 the ordered scaling exponents per approximation level. The quantities D n displayed there are defined as follows, and the scaling exponents ϑ n refer to (6.1). The meaning of (7.3) is that D n (N ) gives the n th largest eigenvalue at approximation order N . In Fig. 8, the line D 1 as a function of the approximation order N thus shows how the largest eigenvalue increases with increasing approximation order provided N has become sufficiently large. Similarly, we observe that all lines D i grow with increasing approximation order. This result establishes very clearly that the addition of higher order invariants N → N + 1 systematically adds a new eigenvalue to the spectrum which is larger, and thus less relevant, than those already present. We conclude that the asymptotically safe fixed point detected here is fully consistent with our working hypothesis.

C. Large order behaviour
In [19,57] is has been observed for f (R) type theories of gravity that higher order powers in the Ricci scalar lead to scaling exponents which become increasingly Gaussian. The result as such is quite intriguing, suggesting that the quantum dynamics of gravity makes the theory "as Gaussian as it gets": the scaling of a few invariants with a low mass dimension receives strong corrections with scaling exponents deviating strongly from classical values. On the other hand, the scaling of invariants with large canonical mass dimension receive only very mild corrections leading to near-Gaussian exponents. Exact Gaussian scaling is clearly not achievable owing to the perturbative non-renormalisability of the theory. Here, we want to check whether this picture persists in case the higher order interactions are substantially different from those investigated in [19,57].
To that end, we display the entire eigenvalue spectrum (6.1) for all approximation orders in Fig. 9. The colour coding (from blue to red) indicates growing approximation order N . We make two observations. Firstly, three negative eigenvalues are neatly visible for all approximation orders. They differ by order unity from their classical values. Secondly, with increasing n, we observe that the eigenvalues approach Gaussian values from above. Both of these features agree with the picture detected in f (R) type approximations [19,57]. The near-Gaussian behaviour thus appears to be a feature of quantum gravity rather than a feature of or limitation to specific technical approximations.

D. Origin of near-Gaussianity
In Fig. 10, we further investigate the origin for the near-Gaussian behaviour. Shown are the scaling exponentsθ m , which, for each approximation order N , correspond to the largest eigenvalue within the set of eigenvalues (6.1) or its real part if complex, We also indicate if the largest exponent is real (full dot) or a complex conjugate pair (crossed dot). We observe that with increasing m (or N ) the largest eigenvalue rapidly approaches Gaussian values. This pattern is qualitatively similar to the pattern observed for f (R) type approximations [57]. As opposed to the f (R) models, the imaginary part of high order exponents comes out smaller. For this reason, the approach to near Gaussianity is largely insensitive to wehther the largest exponent has a small imaginary component, or not.
In Fig. 11 we investigate the rate at which near-Gaussian behaviour is achieved in the scaling exponents. Shown is a semi-logarithmic plot for the relative shift of the eigenvalues away from Gaussian values, introducing v n (N ) = 1 − Re ϑ n (N ) ϑ G,n .
The colour-coding of the data shows the trend that |v n (N )| decreases with increasing N . Based on the data up to 1/N max ≈ 0.05, we conclude that (7.5) resides in the 10-20% range, decreasing with increasing n. In addition, we estimate the asymptotic behaviour of (7.5) by taking the average values for each n over all approximation orders N . These are indicated in Fig. 11 by the dashed black line to guide the eye. The mean values show a much smoother dependence on n, slowly decaying with increasing n. Their envelope is fitted by a simple exponential, v n ≈ v · exp − n n e . (7.7) All mean values from n > 5 onwards, and most entries from the high-order data sets, are below the The significance of (7.7) with (7.8) is as follows. The parameter v is a measure for the mean relative variation in (7.5) at low n, and the parameter n e states at which order the relative variation becomes reduced by a factor of e. With N max /n e ≈ 3, the reduction at N max is by a factor of 20, consistent with (7.6). The important piece of information here is that the data shows a consistent and fast asymptotic decay towards near-Gaussian values. Extrapolation of (7.7), (7.8) predicts that v n (N ) → 0 (7.9) for sufficiently large n, and 1/N → 0. We also note that the apporach towards near-Gaussianity is much faster than in f (R) type theories. In the later case, v ≈ 0.22 and n e ≈ 46.7 has been found [57]. We may conclude that the dynamics of Ricci tensor fluctuations strengthens the near-Gaussian behaviour of higher order invariants. As a last comment, we stress that near-Gaussian eigenvalues are not mandatory for the asymptotic safety conjecture to apply [19,57]. In fact, even more substantial modifications of eigenvalues up to v n (N ) < 1 (7.10) at large orders with intrinsically non-Gaussian corrections would still be compatible with the asymptotic safety conjecture. In this light, the quantum modifications of the high-order eigenvalues at the fixed point of our model are moderate.  Table 5. Shown are the three real anti-de Sitter and de Sitter solutions R − < 0 < R + < R ++ within the radius of convergence of the underlying expansion and their dependence on the approximation order N .

VIII. FROM ASYMPTOTIC SAFETY TO COSMOLOGY
Cosmology, thanks to the wealth of data available from observation [105][106][107], offers an important territory to test the asymptotic safety scenario for gravity. Provided that asymptotic safety is realised in nature, it is conceivable that the characteristic quantum gravitational modifications have impacted during the very early universe, including its phase of inflationary expansion and the phase of late time acceleration. A number of studies have explored these possibilities by exploiting characteristics of an asymptotically safe fixed point using renormalisation group improvements of the effective action or of the gravitational equations of motion including those of Friedmann-Robertson-Walker universes [48,[108][109][110][111][112][113][114][115][116][117][118][119][120][121][122][123][124][125][126].
In this section, we apply our findings to models of cosmology. Our main interest relates to the availability of cosmological scaling solutions such as inflation in the very early universe, or possibly accelerated late time expansion.

A. Stationarity of the quantum effective action
Given the fast convergence of the polynomial fixed point solution, we now ask the question whether the UV fixed point solves the equation of motion. Since we have evaluated the flow equation on the sphere, solutions to the equation of motion correspond to (Euclidean) de Sitter spaces which give a stationary point of the action. Specifically, varying the action (4.2) with respect to the metric and evaluating it on a sphere we obtain The requirement of stationarity, meaning that δΓ/δg µν vanishes, entails the "equation of motion" Alternatively one obtains the very same condition by first evaluating the entire action on a sphere with constant Ricci curvature, and then searching for stationary points of the "potential" Γ(R). This leads to where the prefactor 24π R 2 originates from the volume of the four-sphere. Requiring stationarity of the potential (8.4), ∂Γ(R)/∂R = 0 with R = 0 then leads to (8.2) with (8.3). Its solutions R dS identify stationary points for both (Euclidean) de Sitter and anti-de Sitter corresponding to positive and negative R dS respectively. The set of solutions with positive R dS are of relevance for cosmological scenarios with inflation and determine the Hubble parameter [68].

B. De Sitter and anti de Sitter solutions
For the polynomial fixed point solutions, the de Sitter conditions are polynomial as well and can possess several roots as we increase the order N . We can look for solutions of (8.3) at each order N of the fixed point solution by plotting Γ * (R) and looking for stationary points. Overall, we observe several real solutions, with The positive real solutions R + and R ++ correspond to a minimum and a maximum of Γ * , and the additional negative root R − to a minimum. In Tab. In Fig. 12 we demonstrate that the full numerical solution shares the same stationary points We have also searched for solutions of (8.3) in the complex field plane, see Fig. 13. At each and every order in the polynomial expansion a large number of complex roots are found. Concretely, once N > 3, the equation of motion (8.3) has N − 1 roots in the complex plane. All of these solutions are shown in Fig. 13 for all approximations N ≤ N max . We also indicate the radius of convergence set by the singularity, R c , The radius R s is estimated from the polynomial expansion of (8.3) using the same techniques as in [68]; the band between R c and R s thus serves as an error bar.
We emphasize that most of the roots are unphysical in that they arise outside the domain of validity where |R| > R s . Moreover, all stable and reliable solutions are on the real axis (highlighted by arrows in Fig. 13). Most importantly, both positive roots are firmly within the domain of validity, while the stable negative root resides at the boundary.
For f (R)-type theories within the same set-up, stable solutions in the complex field plane where found, close to the real axis [68].

IX. COMPARISON
We are now in a position to evaluate the impact of Ricci tensor invariants in the effective action. Our results are based on the effective action where f k (R, R µν R µν ) ≡ F k (R µν R µν ) + R Z k (R µν R µν ), see (1.1), (4.2). Findings are compared with closely related studies of quantum gravity for an f (R)-type action of the form see [19,57,68]. The main point is that the exact same regularisation schemes and momentum cutoffs have been used in either case. Therefore, qualitative and quantitative differences can now unambiguously be associated to the dynamical differences in the setup (9.1) vs (9.2), that is, the presence (or absence) of Ricci tensors in the effective action. Our results are summarised in Tab. 6.
A few comments are in order: a) UV critical surface and scaling. On the basis of either of (9.1) and (9.2), the UV critical surface is three-dimensional, and the UV fixed point satisfies the bootstrap test. Canonical mass dimension continues to offer a "good" ordering principle [57]: classically relevant couplings receive strong quantum correction and remain relevant at an interacting UV fixed point. Classically irrelevant couplings remain irrelevant, and, moreover, their scaling exponents become increasingly near-Gaussian with increasing canonical mass dimension. The main difference relates to the rate with which exponents approach near-Gaussian values, which is found to be substantially larger as soon as Ricci tensor invariants are retained. Finally, the classically marginal couplings, corresponding to the invariants √ gR 2 and √ gR µν R µν respectively, both become relevant in the quantum theory with scaling exponents of order unity.
b) Speed of convergence. Within polynomial expansions of either of (9.1) and (  for f (R) gravity with (9.2) the rate of convergence is slow, owing to an eigth-fold periodicity pattern hidden underneath the fixed point equation [19,57]. It has also been shown that the slow rate of convergence is ultimately related to the presence of a near-by pole in the complex field plane [68]. Poles in the complex field plane are regular encounters in nonperturbative fixed points [99][100][101][102][103]. Including Ricci tensor invariants, (9.1), we have observed a substantially faster convergence both for the fixed point couplings as well as for the scaling exponents, see Figs. 3, 5,and 6, and Tab. 4. We conclude that Ricci tensor fluctuations have had an impact on the singularity structure of the UV fixed point solution.
c) Radius of convergence. The renormalisation group flow for either of (9.1) and (9.2) displays a technical pole at R p ≈ 2.00648 · · · , see (5.22). However, the polynomial radius of convergence for f (R) gravity with (9.2) was found to be much smaller [68], R c ≈ 0.82 . . . 0.91. On the other hand, evaluating the flow for the Ricci tensor action (9.1) on spheres with constant curvature, we found that the radius of convergence is maximal, R c = R p ≈ 2.00648, and much larger than for f (R) gravity in an otherwise identical setup. The fact that the maximal range in R is exhausted indicates that the fixed point effective action for (9.1) does not display convergence-limiting singularities in the complex field plane within |R| < R p . We conclude that Ricci tensor interactions stabilise the theory.  highest order coupling. Unlike f (R) theories, which contain a scalar degree of freedom in addition to the polarisation of the graviton, the inclusion of terms such as R µν R µν leads to a fourth order inverse propagator for the graviton. Consequently, as soon as Ricci or Riemann tensor interactions are present, the corresponding equation becomes quadratic. The reason for the latter is that the contribution from the transverse traceless part h T µν contains the highest-order coupling itself. Ultimately, the appearance of the highest-order coupling in the contribution of h T µν comes from the fact that its second variation contains a fourth-order differential operator 2 , a term which is absent in the f (R) case. For further aspects of asymptotic safety in fourth order gravity models, see [33,34,37,70,80].
f) Recursive structure. In the f (R) case, it is straightforward to set up an iterative recursive fixed point solution for all polynomial couplings, as done in [19,57]. In principle, the same philosophy could be used here as well, except that roots arise at the leading order (see Sec. V C). The proliferation of nested roots, after iteration, makes a computer-algebraic solution of the recursive relations cumbersome. We have then limited ourselves to approximations up to the 21 st order [85,98]. Although the approximation order is much lower than the 35 orders achieved in the corresponding f (R) study [19,57], this is more than compensated by the substantially enhanced rate of convergence, see points a), b), c).
In conclusion, we have established that the dynamics of the Ricci tensor has a strong impact onto fixed point solutions. Most noticeably, the UV fixed point is stabilised and convergence is faster. Moreover, de Sitter solutions now arise very naturally within the domain of validity. Clearly, f (R) models for quantum gravity and the model put forward here are both "ad hoc" inasmuch as a full quantum theory of gravity is expected to contain additional invariants not account for here. It is encouraging that the broad features are not affected by these choices. For future work, it will be interesting to extend these investigations for actions involving functions of Riemann and Ricci tensor invariants, and to explore the extend to which the UV fixed point depends on it.

X. SUMMARY
We have put forward a systematic study of the asymptotic safety conjecture for gravity using actions beyond Ricci scalars. The main novelty is that our models are sensitive to the dynamics of more complicated tensor invariants, while still benefitting from the simplicity of f (R) models for gravity. We have illustrated our approach for actions of the form (4.1), (4.2) involving functions of the squared Ricci tensor R µν R µν . The results are quite promising: the theory displays a stable interacting UV fixed point which, in high-order polynomial approximations, shows fast convergence with the order of approximation (Tab. 3). The radius of convergence comes out maximal (Tab. 6), and much larger than the one found in f (R)-type models of gravity. Moreover, for small curvature and within the radius of convergence, the UV fixed point is identified non-perturbatively without resorting to a polynomial approximation (Fig. 4).
The theory has three relevant parameters related to the cosmological constant, the Ricci scalar and the square of the Ricci tensor. Higher order invariants become increasingly irrelevant, and scaling exponents converge very fast (Tab. 4), much faster than in theories with Ricci scalars only. It is also found that the theory becomes near-Gaussian with increasing mass dimension of curvature invariants (Fig. 8, 9). Clearly, higher order curvature invariants stabilise the fixed point, and the theory becomes "as Gaussian as it gets" [19,57]. Our findings seem to confirm once more that quantum scale invariance of gravity can be tested self-consistently by means of a bootstrap, despite of the fact that no small parameter has been identified. Also, the expansion in curvature invariants according to canonical mass dimension is viable. On the other hand, we find that the precise nature of higher order curvature invariants might not be centrally important for a fixed point to arise.
We also investigated implications of our model for quantum cosmology and the existence of inflationary solutions in the early universe. Most notably, we find that the UV fixed point of quantum gravity admits stable de Sitter solutions (Tab. 5). Here, the Ricci tensor fluctuations play a decisive role: without them, such as in the f (R) approximation [68], de Sitter solutions do not arise (Fig. 14). We conclude that the existence of de Sitter solutions is strongly sensitive to the dynamics of metric fluctuations in the deep UV.
Our setup can be developed in several directions. It is straightforward to extend it towards more general theories and to explore the impact of Riemann or Weyl tensor invariants, also using general Einstein spaces for the projection. In the same spirit, the role of matter fields on the gravitational fixed point can be explored [81,98,127]. Further directions for improvement include alternative routes for the projection [69], ideas from the resummed -expansion [74], and the role of the background field [88,128]. It would be particularly interesting to clarify how findings persist in settings that distinguish between background and fluctuation fields (see [78,82] for recent progress), or which are manifestly background independent (see [129] for related advances in gauge theories). Ultimately, the bootstrap search strategy can be used with either of these to test, or even falsify, the asymptotic safety conjecture for gravity. These topics are left for future work.
with ∇ µ V T µ = 0. Similarly, any symmetric tensor h µν is decomposed according to subject to the constraints so that h T µν is the transverse-traceless part of h µν , ξ µ is a transverse vector, σ is a scalar and h is the trace part of h µν . From now on we use the notation (2T ) for a transverse-traceless symmetric tensor, (1T ) for a transverse vector.
In order to see how this affects the calculation, we need to know how the coefficients b n are modified when the operator is restricted to act on (2T ) tensors or on (1T ) vectors which can be related to the spectrum of the unconstrained fields. We start with the transverse vectors and we note that the spectrum of every vector can be expressed as the union of the spectrum of a (1T ) vector and of the longitudinal mode ∇ µ η. The spectrum of ∇ µ η can be related to that of the scalar field η through commutation of covariant derivatives We note however that for a constant scalar η, the vector V µ receives no contribution from the longitudinal mode. So we have to subtract from the scalar trace the constant mode. Thus we write for the trace of a transverse vector Tr (1) e t∇ 2 = Tr (1T ) e t∇ 2 + Tr (0) e t( where the last term corresponds to the zero mode of the scalar field. Thus we can relate the spectrum of the transverse vector to that of the unconstrained vector. Now we turn our attention to the constraints of the transverse traceless tensors. For the symmetric tensors we use the commutation relations As in the case of the transverse vector there are modes that do not contribute to the trace. These modes are (i) the d(d+1)

2
Killing vectors for which ∇ µ ξ ν + ∇ ν ξ µ = 0 so that they do not contribute to h µν ; (ii) a constant scalar σ as in the case of vectors and (iii) the d + 1 scalars which correspond to the second lowest eigenvalue of −∇ 2 − 2 d−1 R. Thus we have for the trace of a symmetric tensor Tr (2) Again we can relate the spectrum of the constrained transverse traceless tensor to that of rank-2 tensors, vectors, and scalars. In order to clarify how the contributions from the exponents play a role in our calculation we expand the exponential in powers of R such as ∞ m=0 c m R m . Taking into account that the volume of the sphere goes like V ∼ R −d/2 and that the heat kernel coefficients go like b n ∼ R n/2 we find that Since ultimately we are interested in comparing powers of R, the exponentials contribute when 2m = n − d, and the coefficients b n receive contributions only when n ≥ d. Another way to see where the excluded modes enter is from the expansion e −tz = 1 − tz + 1 2 t 2 z 2 + .... In order for the parameter t to be included in Q i [W ] = ∞ 0 dt t −iW (t) we see directly from (3.5) that the corresponding power m of the expansion ∞ m=0 c m R m is such that 2m = n − d.

Heat-kernel coefficients
Here we summarise the trace of the heat kernel coefficients tr s b n = b n | s for the fields that we will be interested in after taking into account their constraints. We write 0 for the scalars, 1 for the vectors and 2 for the tensors. For the scalars we have For the transverse vector fields we have  Table 8. Summary of the discrete eigenvalue spectrum Λ l of the operator −∇ 2 on scalars, transverse vectors and transverse-traceless symmetric tensors (s = 0, 1, 2 respectively) and their multiplicities D l for dimension d labeled by the parameter l. In Table 8 we summarise the eigenvalues Λ l of the Laplace operator −∇ 2 on scalars, transverse vectors and transverse-traceless symmetric tensors and their multiplicities D l . The flow terms appearing in (C2) arise through the Wilsonian momentum cutoff ∂ t R k , which we have chosen to depend on the background field. All the terms I 0 [f, z], ..., I 5 [f, z] arise from tracing over the fluctuations of the metric field for which we have adopted the transverse traceless decomposition. The term I 0 [f, z] also receives f and z-independent contributions from the ghosts and from the Jacobians originating from the split of the metric fluctuations into tensor, vector and scalar parts. To indicate the origin of the various contributions in the expressions below, we use superscipts T, V, and S to refer to the transverse traceless tensorial, vectorial, and scalar origin respectively. Then we have for the various components