Asymptotic freedom and safety in quantum gravity

We compute non-perturbative flow equations for the couplings of quantum gravity in fourth order of a derivative expansion. The gauge invariant functional flow equation for arbitrary metrics allows us to extract $\beta$-functions for all couplings. In our truncation we find two fixed points. One corresponds to asymptotically free higher derivative gravity, the other is an extension of the asymptotically safe fixed point in the Einstein-Hilbert truncation or extensions thereof. The infrared limit of the flow equations entails only unobservably small modifications of Einstein gravity coupled to a scalar field. Quantum gravity can be asymptotically free, based on a flow trajectory from the corresponding ultraviolet fixed point to the infrared region. This flow can also be realized by a scaling solution for varying values of a scalar field. As an alternative possibility, quantum gravity can be realized by asymptotic safety at the other fixed point. There may exist a critical trajectory between the two fixed points, starting in the extreme ultraviolet from asymptotic freedom. We compute critical exponents and determine the number of relevant parameters for the two fixed points. Evaluating the flow equation for constant scalar fields yields the universal gravitational contribution to the effective potential for the scalars.

Quantum gravity can be formulated as a consistent quantum field theory for the metric if a fixed point for the flow of (generalized) couplings exists. If this fixed point is approached in the extreme ultraviolet, the quantum field theory is complete in the sense that it can be extrapolated to arbitrary short distances [1]. In short, one defines the microscopic theory by such a fixed point. The relevant parameters for small deviations from the fixed point correspond to the free parameters of the model, reflected by "renormalizable couplings". Their number is typically finite, such that the model is predictive. The fixed point can be either a free theory-this is the case for asymptotic freedom. In contrast, the case of non-vanishing interactions at the fixed point is called asymptotic safety [1][2][3]. See Refs. [4][5][6][7][8][9][10][11][12][13][14][15][16] for reviews. For asymptotic freedom the model is perturbatively renormalizable, while the case of asymptotic safety corresponds to a non-perturbatively renormalizable theory unless all dimensionless couplings at the fixed point are small.
For pure higher derivative gravity, based on terms quadratic in the curvature tensor R µν and the curvature scalar R, Stelle has shown perturbative renormalizability for the couplings C −1 and D −1 [17]. (In addition, the coefficient of a term linear in the curvature scalar and a cosmological constant are renormalizable couplings.) The perturbative renormalization flow for the couplings C −1 and D −1 has been computed by Fradkin and Tseytlin [18][19][20], and reads (3) These equations show a fixed point with asymptotic freedom, to which we refer as the SFT-fixed point. Quantum corrections will lead to a quantum effective action Γ which involves additional terms. In particular, a term linear in the curvature scalar is needed for any realistic theory of gravity. Expanding up to fourth order in derivatives (and omitting the Gauss-Bonnet topological invariant) a diffeomorphism invariant action for pure gravity takes the form where the (reduced) Planck mass M p and the cosmological constant V are additional relevant parameters beyond C −1 and D −1 . At the fixed point the role of V and M 2 p is negligible if these couplings are finite. Flowing away from the fixed point, however, C and D become finite. In this case M 2 p is no longer negligible. For both signs of M 2 p and arbitrary V the effective action (5) leads to tachyons and/or ghosts, such that Minkowski space is no longer a stable approximate solution of the field equations for V M 4 p . As a consequence, the effective action (5) seems not to be compatible with observation. A priori, it is not known if this is a shortcoming of the approximation to the effective action (truncation) or a basic flaw of quantum gravity based on the SFT-fixed point. Indeed, ghosts and tachyons can be artifacts of insufficient truncations [21,22].
There has been recently an intense discussion of perturbative quantum gravity based on the SFT-fixed point, dubbed "agravity" [23,24]. See also Refs. [25][26][27]. The conclusions are severely limited, however, by the simple fact that the couplings C −1 and D −1 flow outside the perturbative domain as the renormalization scale is lowered. For a judgement of the fate of higher derivative gravity it seems compulsory to understand the flow of couplings in the non-perturbative domain. Only this can give an answer to the question of stability for the effective action. Non-perturbative flow equations based on functional renormalization [2,28,29] seem the appropriate tool for such an investigation.
A simple "Einstein-Hilbert truncation" of the flowing effective action or average action Γ k omits the quartic terms ∼ C, D in Eq. (5). One finds a fixed point in the flow of the dimensionless couplings w = M 2 p /2k 2 and u = V /k 4 , often called Reuter fixed point and named here R-fixed point. The R-fixed point is associated to asymptotic safety and the gravitational interactions do not vanish at the fixed point. Its existence corresponds to non-perturbative renormalizability of quantum gravity as a quantum field theory for the metric. The R-fixed point remains present for large classes of extended truncations . It is also present if matter couples to gravity , as for the minimal standard model coupled to gravity. The non-perturbative character of the fixed point implies that the number of relevant parameters is no longer determined by the canonical dimension of couplings. This leads to an enhanced predictivity, as demonstrated by the successful prediction of the mass of the Higgs boson [97]. It is well possible that other parameters of the standard model may become predictable [98][99][100][101][102][103][104][105][106][107][108][109][110][111][112][113][114][115][116][117].
In the present paper we ask what is the relation between the SFT-and R-fixed points. We will use a truncation for which both fixed points are found. This opens the terrain for many interesting questions. Is the SFT-fixed point viable for a definition of quantum gravity? Is there a critical trajectory from the SFT-fixed point to the R-fixed point? What are the implications for the predictivity of quantum gravity coupled to matter? The present paper will not yet fully answer all these questions. It focusses on the derivation of the relevant flow equations for higher derivative gravity and provides for partial answers to these questions within the given truncation.
The minimum truncation needed for these questions takes for the effective average action the from (5), with four running couplings C, D, M 2 p and V depending on the renormalization scale k. Deriving these flow equations within functional renormalization is a technical challenge. There has been previous work for reproducing the perturbative β-functions for C and D [18][19][20][118][119][120][121][122][123], circumventing the technical issues by a rather special field-dependent gauge [55] or using directly expansions in small C −1 and D −1 [30,36]. We aim here for an understanding of the full non-perturbative flow equations for arbitrary values of the couplings C, D, M 2 p and V . This is needed in order to see both the SFT-and R-fixed points. The use of the gauge invariant flow equation [124][125][126] constitutes an important advantage, since the contributions from the physical fluctuations in the metric can be separated from the universal "measure contribution" of the gauge fluctuations. Still, a major technical issue is related to mode mixing. For particular classes of metrics, as Einstein spaces, the transverse-traceless (TT) fluctuations (t-fluctuations) cannot mix with the physical scalar degree in the metric (σ-fluctuations) due to symmetry. Deriving flow equations for Einstein spaces permits a rather straightforward application of known heat kernel expansions. One can, however, only obtain a flow equation for the linear combination D + 6C in this way. The required flow equations for C and D separately require geometries beyond Einstein spaces for which the t-and σ-modes mix.
We display the more technical parts in various appendices and concentrate in the main part on summaries of results and a discussion of crucial features for the derivation of the flow equations. In Section II we present the setup, the characteristic features of the flow equations and the fixed points. Section III summarizes the heat kernel method and the flow contributions from matter fluctuations. In Section IV we exhibit our central result on the flow contributions from the metric fluctuations. Section V discusses asymptotic safety for the R-fixed point. In Section VI we address the infrared region and conclusions are presented in Section VII.

II. SUMMARY: SETUP, FLOW EQUATIONS AND FIXED POINT STRUCTURE
In this section, we summarize our setup, the resulting flow equations and the corresponding fixed points or scaling solutions. The degrees of freedom are the metric and scalar fields. For the effective action we choose a truncation with five coupling functions U (ρ), F (ρ), C(ρ), D(ρ) and E(ρ), where ρ is an invariant formed from the scalar fields without derivatives. The function U (ρ) constitutes the effective potential for scalar fields, F (ρ) is a field-dependent effective squared Planck mass, and C(ρ), D(ρ), E(ρ) are coupling functions multiplying the gravitational invariants involving four derivatives of the metric. For a setting without a scalar field background one simply omits the dependence of U , F , C, D and E on ρ. Besides the metric and scalar fields we also include the fluctuations of fermions and gauge bosons. For the matter fields (scalars, fermions, gauge bosons) we do not compute the flow of derivative terms (kinetic terms), or interaction terms as Yukawa couplings.

A. Effective action for gravity
Our ansatz for the truncated effective average action consists of a gravity part and a matter part For the gravity part, we consider the following truncated effective action, where R is the curvature scalar, C µνρσ is the Weyl tensor whose squared form is given by C µνρσ C µνρσ = R µνρσ R µνρσ − 2R µν R µν + R 2 /3. For the computation of the universal measure contribution to the gauge invariant flow equation we use the equivalent physical gauge fixing for which Γ gf and Γ gh are the gauge fixing and the ghost action for diffeomorphisms. These are given in Appendix A, Eqs. (A5) and (A6). The last term in the square brackets in Eq. (7) is the Gauss-Bonnet term which reads For constant E this is a topological invariant, while for a dynamical scalar field it contributes to the field equations. The effective action (7) contains the most general diffeomorphisms invariant terms for the metric with up to four derivatives. In Appendix B 1 we present the same action in terms of different linear combinations of invariants. The coefficients U (ρ), F (ρ), C(ρ), D(ρ) and E(ρ) are functions of real singlet-fields ϕ, ρ = ϕ 2 , N -component real fields φ a , ρ = φ a φ a /2, or complex scalar fields ϕ a , ρ = ϕ † a ϕ a . One can expand these coefficient functions into polynomials of ρ: Here V is the cosmological constant, m 2 is the scalar mass parameter and λ is the quartic scalar coupling. In the gravitational sector M 2 p is the Planck mass squared for ρ = 0 and ξ is the non-minimal coupling between the scalar field and the curvature scalar. The latter plays a crucial role for the realization of Higgs inflation [127][128][129][130][131].
The matter part Γ matter k consists of canonical kinetic terms for all matter fields. We also include gauge and Yukawa couplings, but set the effects of interactions and masses to zero in many parts of this work. We consider N S scalar bosons, N V vector bosons and N F Weyl fermions. The explicit form of the action for the matter part is given in Section III.

B. Flow equations
The coupling functions U , F , C, D, E depend on the renormalization scale k. Their k-dependence is determined by a truncation of the exact functional flow equation [28,29,[132][133][134][135]. The result of our computation can be written in the form Here the last two terms are contributions from metric fluctuations, where π (t,σ) denotes contributions from the TT tensor (t-mode) and the physical scalar metric fluctuation (σ-mode), while contributions from the gauge modes (the longitudinal modes in metric fluctuations and the ghost fields) are included in δ (g) . We employ dimensionless quantities The flow equations at fixedρ read for ∂ t = k∂ k Here M i denote the contributions from metric fluctuations. Their explicit form is displayed in Section IV B. The structure of the gravitational contributions to the flow equations can be understood by writing them in the form The interpolating functions L mix (m 2 t ,m 2 σ ) represent contributions from the t-mode, the σ-mode and the t − σ mixing, respectively. They are functions of the dimensionless mass terms for the t-and σ-modes where The interpolating functions are linear combinations of the "threshold functions", with whose explicit form can be read off from Section IV B, Eqs. (96)- (111). The explicit forms of L The validity of the flow equations is restricted tom 2 t > −1, while the issue form 2 σ is more complex. For Einstein spaces, one has R µν = (R/4)g µν such that C µνρσ C µνρσ = R µνρσ R µνρσ − R 2 /6 and G 4 = R µνρσ R µνρσ . For these geometries one has only two independent invariants In turn, an evaluation of the flow equation on Einstein spaces, which is done in most work in the literature, only yields flow equations for D + 6C and D + 2E. One finds that for these linear combinations, the mixing terms L (t−σ) mix (m 2 t ,m 2 σ ) cancel out. A computation of separate flow equations for C and D is not possible in this way. For their extraction more general geometries have to be included, and the mixing term plays a role.
The coupling functions depend on two variablesρ and k. We will restrict the discussion here to the case where either the k-dependence or theρ-dependence is neglected, such that the coupling functions depend either onρ or on k. The k-independent functions u(ρ), w(ρ) etc. define a "scaling solution", that generalizes the notion of a fixed point for a finite number of couplings. On the other hand, we may evaluate the flow with k atρ = 0, such that the terms ∼ρ ∂ρu etc. can be omitted in Eqs. (16)- (20). In this case, which we discuss in the following, one is left with the β-functions for five couplings. Fixed points correspond to zeros of these β-functions. We emphasize that the flow with k atρ = 0 can be directly mapped to theρ-dependence of the scaling solution. One simply replaces ∂ t by −2ρ ∂ρ. Thus the flow away from fixed points atρ = 0 translates directly to theρ-dependence of the scaling solution. Keeping this connection in mind, we omit in the following the term ∼ρ ∂ρ etc. in the flow equations. Then the coupling E does not appear in the flow equations for the other couplings, as appropriate for a topological invariant.

C. Asymptotic freedom
The flow equations (16)- (20) admit the SFT-fixed point Since w remains finite at the fixed point, this also implies Finite values of w and v play no role precisely at the fixed point since their relative contribution tom 2 t andm 2 σ vanishes. Close to the SFT-fixed point we can evaluate the flow equations in the limit w → 0, u → 0 orm 2 t → ∞, m 2 σ → ±∞. In this limit, the beta functions for the higher derivative couplings in Eqs. (23)- (25) read These results agree with the perturbative computation [123] in higher derivative gravity. They are universal, i.e. independent of regularization and gauge parameter choice. One infers the flow equation for the perturbative coupling For vanishing matter effects (N S = N V = N F = 0), Eq. (3) is reproduced. For arbitrary numbers of matter particles this is of the asymptotically free form with positive D −1 increasing as k is lowered. There is unavoidably a range where D −1 becomes large and the perturbative result is no longer valid. For the flow equation for C −1 one has Setting N S = 0 yields Eq. (2). For C, D ≥ 0 all terms are positive, driving C −1 towards smaller values as k-decreases while C increases according to Eq. (32). In the range of negative C −1 the ratio ω = 3C 2D (37) has two fixed points according to the flow equation Both fixed points occur for negative ω and therefore negative C for positive D. For pure gravity the numerical values are At the fixed point the ratio ω is not determined since the r.h.s of Eq. (38) vanishes ∼ D −1 for arbitrary values of ω. For the flow trajectories away from the fixed point ω increases with decreasing k for all trajectories with "initial values" ω > ω We conclude that the fixed points (39) concern the behavior of trajectories away from the fixed point, but do not fix the ratio C/D at the fixed point. Starting close to the fixed point with C < 0, D > 0 seems to be rather natural since in this case the Euclidean action (6) is bounded from below. We conclude that both C −1 and D −1 are independent (marginally) relevant parameters. For the Gauss-Bonnet coupling one has at the fixed point These values agree with that found in Ref. [123].
At the fixed point the ratios M 2 p /k 2 and V /k 4 take finite non-zero values, given for pure gravity by At the fixed point one finds Thus w * and v * depend on ω and therefore on the particular trajectory away from the fixed point. For ω = ω UV * one finds the numerical values where for ω = ω Both u and w correspond to relevant parameters, with critical exponents 4 and 2 following from Eqs. (16) and (17) for M U and M F not depending on u and w.
Taking things together the SFT-fixed point has a free undetermined parameter ω. For pure gravity it has four relevant or marginal couplings, with critical exponents given by The Gaussian fixed point characterizes asymptotic freedom of higher derivative gravity [18][19][20][118][119][120][121][122][123]. Finally, we mention the limit of the Weyl invariance which is realized for U → 0, F → 0 and C → 0 in the action (7). For that limit in the pure gravity system, one has Therefore, finite values of u, w and C are induced by quantum effects and the Weyl invariance is not conserved in our current setting.

D. Asymptotic safety
Solving β U = β F = β C = β D = 0 simultaneously, we find a further non-trivial fixed point resulting inm The critical exponents are given as Hence, there are three relevant directions and one irrelevant one. This result agrees with the cases of higher derivative truncations in maximally symmetric spaces, e.g. [53], in Einstein spaces [31], in an arbitrary space within a strong gravity expansion [55] and in a flat space within the vertex expansion scheme [42]. We will argue in Section V that this fixed point is the extension of the R-fixed point for one truncation. The huge absolute values of the critical exponents θ 3 and θ 4 are presumably artifacts of the truncation. Similar high values are observed in a truncation which omits the term ∼ D [48,49]. Including higher order curvature invariants ∼ R n the high value of θ 3 is reduced to a quantity of the order one, and further critical exponents become negative, corresponding to irrelevant parameters [50][51][52][53][54]. The small negative value of D * is not necessarily a cause of worry either. Taking the four-derivative truncation at face value it would imply a tachyon in the t-sector and therefore an instability. The derivative expansion yields, however, only a Taylor expansion of the inverse graviton propagator for small momenta. In the momentum range of the possible instability higher order terms can cure this issue. Furthermore, extended truncations could shift D * to positive values.  i (m 2 t ) we need to fix the two other parameters that we take as ω and v. We show two sets, one for the values of v and ω at the SFT-fixed point with ω = ω (UV) * , and the other for the R-fixed point. The same procedure is used for L Fig. 2. The curves are shown with v taken from the SFT-or R-fixed point, and ω = 0.5 for both curves. Figures for the other interpolating functions can be found in Section IV.
We have chosen normalizations for the interpolating functions such that their overall size is of the order one, except for L (t−σ) mix for the R-parameter set. This allows for a judgement of the size of the different contributions from the prefactors in Eqs. (21)- (25). In particular, for the asymptotically free SFT-fixed point the t-fluctuations dominate, and the t − σ-mixing gives only a small contribution of around 15 percent. In this region an omission of the t − σmixing contributions, which are technically the hardest part, does only lead to rather minor errors. For pure gravity (N S = N V = N F = 0) and atρ = 0 one obtains at the SFT-fixed point for fixed constantm 2 σ /m 2 t = 3C/D, Here the first, second and third square bracket are contributions from the t and σ modes and their mixing, respectively, while the last terms are the measure contributions. The flow of D and C is dominated by the t-fluctuations. The coupling D decreases until the negative contributions to β D from the mixing and σ-fluctuations get large enough to compensate the positive contribution from the t-fluctuations. In our truncation this happens for negative d and negativem 2 t . For the R-fixed pointm 2 t is indeed negative, cf. Eq. (47). At the R-fixed point the interpolating functions remain of the order one, with the exception of the mixing contribution, A fixed point with finite D can only occur for negative D since the decrease of D with decreasing k has to be stopped. The rather large value (54), together with a rather large negativem 2 t , may cast same doubts on the robustness of the R-fixed point with respect to extended truncations.   The flow equations (16)- (20) are valid for arbitrary constant scalar fields, or arbitrary constantρ. In particular, Eq. (16) describes the flow of the scalar potential. In principle, N S , N V and N F depend onρ, reflecting the flowcontributions from matter loops [84]. We focus here on the gravitational contributions encoded in M U (ρ). The central quantity is the gravity induced scalar anomalous dimension A. For A > 0 the quartic scalar couplings become irrelevant parameters. Together with the assumption that the flow of couplings below the Planck scale does not deviate much from the one of the standard model this leads to the successful prediction of the mass of the Higgs boson [97], or more precisely for the mass ratio between top quark and Higgs boson [11,108]. For A > 2 also the scalar mass terms become irrelevant couplings, driving all scalar masses rapidly to zero. This is the basis for a possible solution of the gauge hierarchy problem by the running of couplings [136,137].
The metric-fluctuation induced scalar anomalous dimension is given by the derivative of M U with respect to u for ρ → 0, Here the first and second term correspond to contributions from the t-mode (graviton) and σ-mode in the metric, respectively. In our previous papers [83,84] within the Einstein-Hilbert truncation (C → 0, D → 0 and E → 0), we have found The dominant graviton (the t-mode) contribution agrees with Eq. (55) for D → 0, while the part of the scalar mode has a difference of a factor 9/20 for C → 0. This difference arises from the different regularization scheme for the propagator of the scalar mode. At the SFT-fixed point the gravity induced scalar anomalous dimension vanishes according to On the other hand, for the R-fixed point one finds This is positive and below two. Adding contributions from matter fluctuations can change the value of A. In particular, a decrease of the fixed point value for w leads to an increase of A. Our investigation gives further support to the prediction for the mass of the Higgs boson, while the question of the gauge hierarchy will depend on the precise matter content.

III. FLOW GENERATOR AND HEAT KERNEL METHOD
We compute the flow equations (16)-(20) by the heat kernel method. For this purpose we consider arbitrary geometries close to flat space. Geometries beyond Einstein spaces are needed for the extraction of the gravitational contributions to the individual coupling functions for the higher derivative terms in the gravitational effective action. The evaluation of the heat kernels for the differential operators needed in this context is technically new terrain and will be described in the next section. In this section we introduce the method and present the flow-contributions of matter fluctuations that can be obtained by more standard techniques.
We want to evaluate the flow generator as a general functional of the ("background")-metric field. The sum i is over fluctuation contributions of degrees of freedom that do not mix. For our application the metric is closed to flat space, but not restricted otherwise. For the heat kernel method ζ i is represented as a trace over suitable differential operators, ζ i = tr W (z = ∆ i ) where ∆ i is an appropriate Laplacian acting on the degree of freedom i. The flow generator can be expanded as Here tr (i) acts on the space for the degree of freedom i, and the coefficients c where the values of the heat kernel coefficients b (i) n depend on the degrees of freedom of a field on which the Laplacian acts. The coefficients of (1, − 1 2 R, − 1 2 R 2 , 1 2 C µνρσ C µνρσ , G 4 ) yield (∂ t U, ∂ t F, ∂ t C, ∂ t D, ∂ t E) at fixed ρ. Switching to dimensionless couplings yields the terms −4u and −2w, and the change to the dimensionless scalar invariantρ induces the terms ∼ 2ρ∂ρ in the flow equations (16)- (20).
The detailed steps of these calculations are displayed in the Appendices A-E. In the main text we summarize the most important points. For the evaluation of Eq. (60) the explicit heat kernel coefficients are summarized in Appendix C. The threshold functions Q n are given by The contributions to the flow generator take typically the form where a, b, c are dimensionless couplings and R k constitutes the infrared cutoff. The quantity depends on the normalization of the fields. For contributions from the t-mode and matter modes, one has = 0, while = 1 is used for the σ-mode contributions. For the cutoff function R k (z) we employ a generalization of the Litim cutoff [138] such that z in the denominator of W (z) is replaced by where we neglected the anomalous dimensions, i.e. ∂ t a = ∂ t b = ∂ t c = 0. The full expressions for Q n [W ] with the anomalous dimensions can be found in Eq. (C9) in Appendix C. For more details on the heat kernel technique we refer to Refs. [7,84]. Our general setting will become more explicit once we evaluate next the individual contributions from the fluctuations of free and massless scalar fields, fermions and gauge bosons.

A. Scalar bosons
We start by evaluating contributions from a scalar field whose effective action is given by The second functional derivative with respect to ϕ yields where ∆ S = −D 2 is the Laplacian acting on a spin-0 scalar field and we define the mass term, In order to compute the flow-contributions, an appropriate IR cutoff function R k (∆ S ) is added so that the Laplacian ∆ S in Eq. (66) is replaced to P k = ∆ S + R k (∆ S ). The flow equation for the scalar contribution reads where the flow kernel is Using the heat kernel expansion (60) with the corresponding heat kernel coefficients to ∆ S (see Table I in Section C 3), one obtains For the Litim cutoff the threshold functions are given by Eq. (28). For N S massless scalars (m 2 ϕ = 0) this yields in Eqs. (16)-(20) the contributions ∼ N S .

B. Gauge bosons
We employ the effective action for a gauge theory as where F a µν is the field strength of the gauge field A a µ . Here, Γ gh are the actions of the gauge fixing and ghost fields (c a ,c a ) associated with the gauge field A a µ and are given respectively by with α V the gauge fixing parameter. The Hessian of A a µ , i.e. the second-order functional derivative of Eq. (71) with respect to A a µ , is computed as In the Landau gauge α V → 0, the flow equation for Eq. (71) reads For vanishing gauge couplings the r.h.s. of Eq. (74) becomes the sum of contributions from individual gauge bosons. For a single gauge boson the flow kernel is given by where the regulator R k (z) replaces the Lichnerowicz Laplacian The use of the Litim cutoff and the heat kernel expansion yields the contribution from the physical mode in A a µ as The contribution of the gauge mode and ghost yields for the "physical" Landau gauge the universal measure contribution, The two terms sum up to For N V gauge bosons this yields in Eqs. (16)-(20) the contributions ∼ N V .

C. Weyl fermions
We next consider the contributions from Weyl fermions to the flow of the gravitational interactions. Spinor fields in curved spacetime involve the vierbein fields e m µ . The covariant derivative D µ involves the spin connection constructed from e m µ . With / D = γ µ D µ and e = det(e m µ ) the effective action with a Yukawa coupling to a scalar field ϕ reads One obtains the Hessian as The squared covariant derivative becomes − / D 2 = −D 2 + R/4 = ∆ L 1 2 , which is the Lichnerowicz Laplacian acting on a spinor field. Thus, we regulate z = ∆ L 1 2 by employing the regular R k (z) to obtain the flow generator From the heat kernel technique, one finds with m 2 ψ = y 2 ρ/k 2 . The terms ∼ N F in Eqs. (16)-(20) account for the contributions form N F massless Weyl fermions (m 2 ψ = 0).

IV. FLOW GENERATORS FROM METRIC FLUCTUATIONS
This section summarizes the main technical achievement of this work, namely the computation of the metric fluctuations to the functional flow of all coupling functions of higher derivative gravity in fourth order. To this end, we use the heat kernel method again. The two-point function (Hessian) of the ghost fields can be written in terms of the form in the so-called non-minimal operator for which we use the heat kernel coefficients obtained in Refs. [139][140][141]. On the other hand, the Hessian for the metric field becomes complicated and cannot be simplified to be in such a non-minimal operator. Therefore, we use the off-diagonal heat kernel expansion introduced in Ref. [33] to evaluate the heat kernel coefficients for e.g. tr [e −s∆ R µν D µ D ν ] appearing in the flow kernel for metric fluctuations.

A. Physical metric decomposition and gauge invariant flow
Our starting point for the derivation of the flow equations is to employ the "physical decomposition" [124][125][126] of the metric fluctuations.
withḡ µν the argument of the effective action (often referred to as "background field"). Hereafter we omit the bar on the background field. The physical metric fluctuation f µν satisfies the transverse condition, i.e. D µ f µν = 0. The "gauge fluctuations" a µν denote the directions in field space generated by infinitesimal gauge transformations (diffeomorphisms). In a physical gauge they decouple from the physical fluctuations.
In second order in f µν the expansion of the effective action (7) yields where the covariant derivative operator acting on the physical metric fluctuations f is given by The "interaction piece" M(R,∆ T ) is a tensor depending on curvature tensors and the covariant derivatives. With one sees that T is orthogonal to I, namely T µν ρτ I ρτ αβ = 0. The part ∼ T is the kinetic part of the inverse graviton propagator. The flow equation for the system (with the physical gauge fixing action (A5) for α → 0 and β = −1) takes the form where π (f ) k is contributions from the physical modes f µν , whereas δ (g) k contains contributions from the gauge modes a µν and the ghost fields.

B. Physical metric fluctuations
The physical metric fluctuations f µν can be further decomposed as Here t µν is the TT tensor, i.e. satisfies D µ t µν = 0 and g µν t µν = 0, and σ is the scalar physical fluctuation of the metric defined by The tensorŜ µν obeysŜ such that D µŜ µν = 0 and g µνŜ µν = 1. With σŜ µν t µν = −σ(3∆ S − R) −1 R µν t µν one has, in a general spacetime, mixing effects between t µν and σ. Only in an Einstein spacetime, R µν = (R/4)g µν , the TT tensor mode (t-mode) decouples from the scalar mode (σ-mode).
The Hessian for the physical metric fluctuations takes the following structure: where ∆ T = −D 2 is the Laplacian acting on tensor fields. The kinetic parts are defined by They correspond to the inverse propagators of the t-mode and σ-mode. The interaction parts M are lengthy. Their explicit forms are shown in Appendix B: See Eqs. (B20)-(B24). The last terms involve the mixing between the tand σ-mode. This vanishes for Einstein spaces due to the different transformation of these modes under rotation symmetry.
For the particular case of flat spacetime we can perform a Fourier transformation. The propagators of the graviton and σ-mode are given respectively by with M 2 t = F/(2D) and M 2 σ = F/(6C). The last terms are taken for U → 0, where we observe the usual ghost for the t-mode of fourth-order gravity. We also note the negative sign of the kinetic term for the σ-mode in Eq. (91). Functional renormalization deals with these well known problematic feature by introducing an infrared cutoff.
For the two-point functions, we employ the regulator such that ∆ i is replaced by P k (∆ i ) in the kinetic terms K (t) and K (σ) in Eq. (92). The contribution to the flow generator from the physical metric fluctuations reads Due to the regulator functions for the t-mode and σ-mode the Laplacians (∆ T and ∆ S ) are replaced to P k (∆ i ) = ∆ i + R k (∆ i ). The last term in Eq. (95) arises from the regulated Jacobian which accounts in the functional integral for the decomposition of the metric fluctuations into physical modes and gauge modes. Computations are given in Appendix E. Here, we list the explicit forms of contributions from the t-and σ-modes to the beta functions (16)- (20): and from π We also define the interpolating functions for the mixing effects as In Eq. (102), the Euler characteristic and the number of Killing vectors and conformal ones are denoted by χ E and N , respectively. For geometries with a topology of Euclidean flat space they are given by (26), correspond to the regularized propagators of the t-and σ-modes, multiplied with Dk 4 and 3Ck 4 , respectively. For the particular Litim-type regulator employed here the combination P k (∆ i ) is effectively replaced by k 2 . By virtue of the gauge invariant formulation, or equivalent physical gauge fixing, the contribution of the different physical modes is well visible and separated from the sector of the gauge modes. Despite the somewhat lengthy expressions the structure of the contributions to the flow generator is well visible. The terms involving both factors of (1 +m 2 t ) −1 and (1 +m 2 σ ) −1 are due to the mixing of the t-and σ-modes for geometries that do not exhibit rotation symmetry.
The contributions from the t − σ-mixing correspond to L (t−σ) mix in Eqs. (23)- (25), which can be read off directly from the explicit expressions in Eqs. (98)-(102) and (106)-(111). The remaining parts in these equations define the interpolating functions L  21) and (22). We show these interpolating functions in Fig. 3 for the same parameter sets as for Fig. 1. The flow generator for the effective scalar potential ∼ M U is particularly simple and has a simple expression in terms of one loop diagrams in Fig. 4. It can be equivalently evaluated in flat space, with propagators taken for constant scalar fields. For the other expressions the t-and σ-propagators in Fig. 4 have to be evaluated in a curved background, and there is propagator mixing in the absence of rotation symmetry. A perturbative computation would expand the propagators in deviations from flat space, introducing external legs for metric fluctuations in Fig. 4. The increasing number of legs needed for the flow of higher derivative couplings is directly related to the increasing complexity of the expressions in Eqs. (96)- (111). A Taylor expansion of the functional flow equation in the number of external fluctuation fields produces the vertex expansion for the flow equations, which has been investigated extensively for quantum gravity [38][39][40][41][42][43][44][45][46][47]. In comparison to this work the gauge invariant flow equation has additional contributions from the field dependence of the infrared regulator and requires a particular physical gauge fixing. A more detailed discussion of these issues can be found in Ref. [142,143].

C. Measure contribution
The measure contribution includes the spin-1 gauge modes a µ in the metric fluctuations (a µν = D µ a ν + D ν a µ ) and the ghost modes C µ ,C µ , as well as the Jacobian arising from the decomposition (83). The gauge invariant formulation tells us that there is a simple relation [125,126] between the metric fluctuations δ (g) k and the ghost ones k . Thanks to these relations, the total measure contribution takes the simple form with the differential operator acting on a vector field, Here ∆ V = −D 2 is the Laplacian acting on vector fields. In Eq. (115), we employ a regulator which replaces D 1 by P k . More explicit calculations using the heat kernel method are presented in Appendix E 2. Here we show only the result: As an important advantage of the gauge invariant formulation the gauge sector decouples from the physical sector. As a consequence, this measure contribution is universal in the sense that it does not depend on the couplings in the physical sector. The expression (117) does not involve the coupling function u, w, c, d.

V. R-FIXED POINT AND CRITICAL EXPONENTS
The fixed points of the flow equations (16)- (20) obtain by setting the terms ∼ ∂ t and ∼ρ∂ρ to zero. For pure gravity, with N S = N F = N V = 0, the fixed point conditions are These are four non-linear equations for the four couplings u, v, c and d that we have solved numerically. The numerical investigation finds several fixed points. Part of them are outside the validity of our truncation, occurring for example atm 2 t < −1. We have selected the R-fixed point (47). We will next argue that this fixed point corresponds to the fixed point of asymptotic safety found earlier for other truncations of the effective average action. For this purpose we compute these truncations in our setting for the gauge invariant flow equation and the particular infrared regulator employed in this paper. This allows for comparison with earlier results, and shows how the R-fixed point depends on the truncation.
at which the critical exponents read These values are in a range found for the R-fixed point within earlier truncations. The value for w * is similar to Eq. (47), while u * is substantially larger. This is due to the difference of the values form 2 t andm 2 σ between Eqs. (47) and (119). For the values in Eq. (119) the β-function for D remains positive, "destabilizing" the fixed point (119) once the higher derivative terms are included. As we have argued before, substantially more negative values for m 2 t are needed in order to find a fixed point for D. This is the main reason for the shift in the fixed point values between Eqs. (47) and (119). The fixed point w * in Eq. (119) corresponds to the dimensionless Newton constant g N * = 1/(16πw * ) = 0.92.
Let us finally look at the flow equation for E evaluated for the fixed point (119). Since the Gauss-Bonnet term is topological, its coefficient E does not contribute to any beta functions, whereas the beta function of E depends on other couplings. It takes into account the number of (normal and conformal) Killing vectors, and therefore monitors topological features of the background geometry. We find a vanishing β E for N G = 3.1χ E , but we will not impose β E = 0 for the search of fixed points.

B. R 2 truncation
Next, we analyze the R 2 truncation for pure gravity. The earlier works on that truncation have been done in the Einstein spacetime R µν = (R/4)g µν [31,32,81] or the maximally symmetric spacetime R µνρσ = (R/12)(g µρ g νσ − g µσ g νρ ) [49,[52][53][54]. In the present work, we do not assume such a special spacetime background. We first consider the R 2 truncation defined by setting D = E = 0 and then solving β U = β F = β C = 0. We find a fixed point for at which we have the critical exponents The critical exponents for F and U are close to those of the Einstein-Hilbert truncation, whereas the R 2 coupling has a huge value of the critical exponent. Such a situation has been actually seen in the previous works and indicates the insufficiency of the truncation. By improving the truncation, i.e. including higher order operators such as R 3 , R 4 , etc., the critical exponents converge to reasonable values [53,54].
As an alternative R 2 -truncation, we can look at the flow of the linear combinationC = C + D/6, setting C =C for the flow generators. This corresponds to the previous investigations for Einstein spaces. For this procedure we find a fixed point at This fixed point value gives the critical exponents Comparison with Eq. (121) demonstrates that the fixed point values depend substantially on the choice of the precise truncation, while the first two critical exponents are more robust. The extraction of the flow of C is ambiguous and needs a specification of assumptions on the ratio D/C.

C. Vanishing cosmological constant
Another interesting approximation or truncation neglects the scalar potential or cosmological constant. A numerical search for simultaneous zeros of the beta functions β F , β C , β D for u = v = 0 and N S = N V = N F = 0 shows three For these fixed points the critical exponents are found as FP 1 : θ 1 = 2.10 , θ 2 = 9.51 , The fixed point FP 1 is close to the R-fixed point in Eq. (47), while the critical exponents are sensitive to the omission of u.

D. Full system
For a numerical search of the fixed points for the full system, we look for solutions β U = β F = β C = β D = 0. This yields the R-fixed point in Eq. (47). Whilem 2 t takes substantial negative values as necessary for a stop of the flow of D, there is no sign of a problematic behavior. This is exemplified by the flow of D for which we show the behavior of β D = M D /960π 2 as a function of D in Fig. 5. For all couplings except for D, we use the fixed point value (47). The location of the fixed point value D * is shown by a fat blue dot. While the vicinity of the pole atm 2 t = −1 yields the needed substantial negative contributions to β D , the effect is not dramatic. In the current setup we have found no other fixed point withm 2 t > −1. The fixed points FP 2 and FP 3 in the truncation for u = 0 may therefore be considered as artifacts due to the truncation.

VI. INFRARED REGION
The coupling w corresponds to a relevant parameter, both at the SFT-and R-fixed points. Away from the fixed points w typically increases and one enters the infrared region of large w. In this region the Einstein-Hilbert action becomes a very good approximation. We will discuss this issue by investigating the scaling solution in the region of largeρ. Equivalently, the infrared region can be discussed at ρ = 0 for the k-dependence of the coupling constants.
The scaling solution for largeρ is characterized by w ∼ρ, i.e.
with constant ξ ∞ and w ∞ . If forρ → ∞ the coupling functions u, C, D remain finite or increase only slowly withρ, the leading behavior in the infrared region is given by In this limit the gravitational contributions to the flow generators simplify considerably sincem 2 t andm 2 σ vanish. One finds (with M such that M D and β D take constant values Hereβ i obtains from β i by omitting the term ∼ρ∂ρ. Similarly, one obtains and The scaling solution exists for a constant value of u = u ∞ , such that v = u/w ∼ρ −1 vanishes indeed forρ → ∞. The coupling C increases logarithmically, with constantC adjusted such that the asymptotic behavior (137) for largeρ matches the behavior at the SFT-fixed point forρ → 0. This increase is rather slow, according to the value A C ≈ 0.015 from Eq. (135). Similarly, the coupling function D(ρ) for the squared Weyl tensor decreases logarithmically (A D ≈ 0.031 > 0) Forρ → ∞ both c and d vanish ∼ lnρ/ρ. With F = ξ ∞ ϕ 2 + 2w ∞ k 2 the effective Planck mass depends on the scalar field ϕ. For the cosmology of this type of models it becomes dynamical. For typical solutions ϕ and F diverge for infinitely increasing time [144,145]. The observed Planck mass can be associated with the present value of the scalar field, M 2 = F (t 0 ) = ξ ∞ ϕ 2 (t 0 ). If we choose the renormalization scale k such that u ∞ k 4 coincides roughly with the present dark energy density, u ∞ k 4 ≈ (2 · 10 −3 eV) 4 , the increasing scalar field has reached today a large value ϕ 2 (t 0 )/k 2 =ρ(t 0 ) ≈ 10 60 . This large value is connected to the huge age of the universe in Planck units. The infrared limit forρ → ∞ becomes a very good approximation.
For processes involving finite length scales the associated non-zero momenta typically act as an additional infrared cutoff, replacing effectively in the logarithms of Eqs. (137) and (138), lnρ → ln(ϕ 2 /(k 2 + q 2 )), resulting in The corrections of higher derivative gravity to Einstein's gravity are suppressed by DC µνρσ C µνρσ /M 2 R or CR/M 2 . Associating R ∼ q 2 , C µνρσ C µνρσ ∼ q 4 , and using M 2 /q 2 ∼ wk 2 /q 2 , the suppression factor Dq 2 /M 2 ≈ (D/w)(q 2 /k 2 ) = dq 2 /k 2 ≈ 10 −30 Dq 2 /k 2 is tiny for all momenta sufficiently below the Planck mass. Modifications of the ghost pole for the graviton propagator in higher derivative gravity due to the logarithmic running of D in Eq. (139) have been discussed in Ref. [11]. Instead of a ghost pole on the real axis for q 2 there is a pair of poles in the complex plane. It is not clear if this is compatible with an acceptable graviton propagator. Leaving the issue of the graviton propagator open, a model of quantum gravity based on a scaling solution with a flow inρ from the SFT-fixed point forρ → 0 to the infrared limit forρ → ∞ seems acceptable. This would realize "fundamental scale invariance" [146]. The observable gravitational interactions on length scales smaller than the present horizon would agree with the predictions of Einstein gravity if the matter sector is also characterized by a scaling solution, with all particle masses proportional to the scalar field ϕ and ϕ-independent dimensionless couplings. Observable mass ratios are then time independent despite a cosmological time evolution of ϕ. On cosmological scales the scaling solution with u = u ∞ , w ∼ρ gives rise to dynamical dark energy [144]. By a Weyl scaling to the Einstein frame the effective Planck mass and particle masses become constants and all dependence on the renormalization scale k is eliminated [146].

VII. CONCLUSIONS
Within functional renormalization we have computed the flow of higher derivative gravity with up to four derivatives of the metric. We have coped with this technical challenge by use of the gauge invariant flow equation (or equivalent physical gauge fixing). This has permitted us a decomposition into fluctuation modes for which the individual contributions can be separated, helping for a qualitative understanding of different effects.
In the double limit C −1 → 0, D −1 → 0 we recover the perturbative β-functions and the associated asymptotic freedom at the SFT-fixed point. According to the perturbative β-functions the coupling D −1 increases as the renormalization scale k is lowered and runs outside the perturbative domain. The non-perturbative character of the functional flow equation allows us to follow the flow of D −1 and C −1 in the non-perturbative domain where D and C are no longer large. We find that for D > 0 nothing stops the decrease of D with decreasing k. As a consequence, D reaches zero at some k 0 .
For negative D the sign of the β-function for D can change. The zero of β D at D < 0 permits an additional fixed point with small finite D * < 0. This fixed point lies outside the perturbative domain of higher derivative gravity. At the fixed point also C, the dimensionless effective Planck mass w and the dimensionless effective cosmological constant u take fixed values. This additional R-fixed point can be associated with asymptotic safety. It is possible to formulate a consistent quantum field theory for gravity based on the R-fixed point.
The negative value of D, as well as the effective mass termm 2 t being not very far from the breakdown of the truncation atm 2 t = −1, may cast same doubts on the reliability of our truncation. An extended truncation would be much welcome for a test of robustness. On the other hand, investigations of the momentum dependence of the graviton propagator by different functional renormalization approaches [42,47] suggest that the R-fixed point persists for an arbitrary momentum dependence of the inverse graviton propagator beyond the terms in fourth order in momentum included in the present paper. It is an interesting question if a critical flow trajectory connects the SFT-and R-fixed points.
The existence of the R-fixed point seems not mandatory for a consistent quantum field theory of metric gravity. The microscopic theory may be defined as well at the SFT-fixed point and be characterized by asymptotic freedom. The trajectories away from the SFT-fixed point may join directly the "infrared region" for which an effective theory of gravity based on the Einstein-Hilbert action becomes a very good approximation.
It is possible that the effective action for quantum gravity is described by a scaling solution for which the functions u, w, D and C depend on the dimensionless scalar invariantρ, without explicit dependence on k. In this case, the limitρ → ∞ defines an infrared fixed point, and flow trajectories could connect directly the SFT-and infrared fixed points. The behavior of this scaling solution for largeρ is given by with constant ξ ∞ and u ∞ . While C and D depend logarithmically onρ for largeρ the ratios d = D/w, c = C/w vanish ∼ lnρ/ρ. These ratios are the relevant quantities to measure the deviations from an effective action based on the Einstein-Hilbert truncation.
We may replace for the logarithmic flow of C and D the renormalization scale k 2 by k 2 + q 2 , with q 2 a typical squared momentum q 2 or corresponding covariant (negative) Laplacian ∆ = −D 2 . Takingρ = ϕ 2 /k 2 with ϕ a scalar singlet field, the gravitational effective action according to the scaling solution is given for q 2 ϕ 2 and k 2 ϕ 2 by Together with a suitable kinetic term for ϕ, and particle masses ∼ ϕ as appropriate for the quantum scale symmetry at an IR-fixed point, this effective action seems compatible with observation. The scale k is the only overall scale and we may choose it as k ≈ 2 · 10 −3 eV. The present effective value of the Planck mass is ∼ ξ 1/2 ∞ ϕ. It is dynamical and reaches for present cosmology the value 2 · 10 18 GeV, such that k 2 /ϕ 2 ≈ 10 −60 explains the weakness of gravity. The scalar field ϕ can be associated with the cosmon, providing for dynamical dark energy.
Note added. After submitting our paper on arXiv, a recent paper [147] has computed the heat kernel coefficients for a more general non-minimal forth-order derivative operator, D 4 +Ω abc D a D b D c +D ab D a D b + H a D a +P . Because gravitational systems with higher derivative operators could contain such a type of derivative in the two-point functions, those results may make computations simple and then may be quite useful for investigations of asymptotically safe gravity and its extended systems.

Appendix A: Setup
Our main tool in this work is the functional renormalization group [28,29,133,135]. A central object is the scale-dependent 1PI effective action (or effective average action) Γ k . For k → 0, one obtains the full 1PI effective action Γ k=0 = Γ. The scale change of Γ k is described by the following functional differential equation [28]: Here, R k is a regulator which suppress low momentum modes with |p| < k such that high momentum modes with |p| > k are integrated out, ∂ t = k∂ k is the dimensionless scale derivative, and Γ (2) k is the Hessian, i.e. the full two-point function defined by the second-order functional derivative with respect to field variables. Tr denotes the functional trace acting on all internal spaces of fields such as momenta (eigenvalues of the covariant derivative), flavor and color. See Refs. [148][149][150][151][152][153][154][155][156][157][158][159] on the derivation and technical aspects of the flow equation (A1) and its applications to various systems.
Although the flow equation (A1) is exact, namely it is derived without any approximations, one needs to make approximations to solve Eq. (A1) even for a simple system. In general, the effective average action Γ k includes an infinite number of effective operators generated by quantum effects. Therefore, approximations can be made by restricting an infinite-dimensional theory space into a finite subspace.
In this appendix, we explain our setups and techniques to derive the flow generators or beta functions for the coupling functions in a derivative expansion for quantum gravity in fourth order in derivatives. Firstly, we give the specific form of the effective average action for the gravitational system in order to investigate the asymptotic freedom and asymptotic safety scenarios for quantum gravity. Secondly, properties of the decomposed metric fluctuation field are summarized. Thirdly, we list useful identities for the covariant derivatives acting on various spin fields which are used to evaluate the trace for spacetime indices in the flow equation (A1). Finally, we show the expanded projectors into polynomials of curvature invariants. The decomposition of the metric fluctuation field is realized by acting appropriate projectors on the metric fluctuation field. Those involve generally an infinite number of the covariant derivatives and curvature operators. For our purpose in this work, it is convenient to expand the projectors in polynomials of curvature invariants up to its squared forms.

Effective average action
The starting point to derive the flow equations for the gravitational system is to split the metric field g µν into a background fieldḡ µν and a fluctuation fields h µν : Hereafter, quantities, variables and operators contracted by the background field are presented by a bar on them. The effective average action as a truncated system reads Here the action for the gravity part is parametrized with where the Gauss-Bonnet term G 4 = R 2 − 4R µν R µν + R µνρσ R µνρσ is topological in four dimensional space, and C µνρσ C µνρσ = G 4 + 2R µν R µν − 2/3R 2 = R 2 µνρσ − 2R 2 µν + 1/3R 2 is the squared Weyl tensor. The gauge and ghost actions for diffeomorphisms are given respectively by where C µ andC µ are the ghost and anti-ghost fields and the gauge fixing function obeys with h =ḡ µν h µν the trace mode of the metric fluctuation. The constants α and β are the gauge fixing parameters.
In this work, we use the physical gauge fixing β = −1 and α → 0 for which the gauge fixing action (A5) deals with the path integral for the gauge field on the gauge orbit satisfyingD ν h νµ = 0. The physical gauge fixing acts only on the gauge modes which are generated by the action of a gauge transformation on the background metricḡ µν . For the matter part, we give the action for free N S -scalars, N F -Weyl fermions, N V -gauge bosons, where √ e denotes the determinant of vierbein e µ a . The gauge fixing and the ghost action for the gauge symmetries read where α V is the (dimensionless) gauge fixing parameter for the additional gauge symmetries in the matter sector beyond diffeomorphisms.

Physical metric decomposition
In general, the fluctuating metric h µν is a second rank symmetric tensor, namely it has 10 degrees of freedom in four dimensional spacetime. This tensor can be decomposed into physical and gauge degrees of freedom. In this work, we employ the physical metric decomposition [124] which is given by where f µν is the physical metric satisfying the transverse condition, i.e.D µ f µν = 0, and a µ is the spin-1 vector gauge mode. Thus the physical metric f µν has 6 degrees of freedom, while there are 4 degrees of freedom in the gauge mode a µ , corresponding to infinitesimal diffeomorphism transformations. We introduce the trace mode (the physical spin-0 scalar mode), σ :=ḡ µν f µν , and split the physical metric fluctuations, The two parts read where we have introduced the unit matrix acting on symmetric tensors and the projection tensors, The projection tensors satisfy the conditions as the orthogonal basis, Using these projections, one can define the projection operators P a and P f on the gauge and physical modes, where These projectors satisfy and, at the lowest order, tr P a = 4 , tr P f = 6 . (A19) In Eq. (A17) and the following we denote the Laplacian∆ = −D µD µ acting on scalars, vectors or tensors by∆ S , ∆ V and∆ T , respectively.
The physical metric fluctuations and the gauge modes are further decomposed into irreducible representations of the Lorentz group, where t µν is the transverse-traceless (TT) tensor and κ µ a transverse vectorD µ κ µ = 0, while σ(=ḡ µν f µν ) is the spin-0 scalar field defined above. The projector on the physical fluctuations is defined bŷ This operator satisfies trŜ =ḡ µνŜ µν = 1. The decomposition (A12) can be written as wheres µν satisfies the traceless condition, i.e..ḡ µνs µν = 0 and is explicitly given bỹ We will later need further identities for the physical scalar metric fluctuation σ, aŝ For an expansion linear in the curvature tensor one has and where we define the transverse projector, The projection operators which project out t µν andŜ µν σ from h µν are where Under infinitesimal diffeomorphism transformations for the metric fluctuations, the TT tensor t µν and the physical spin-0 scalar field σ are invariant, whereas the transverse vector field κ µ and the gauge spin-0 scalar field u are not.

Identities for covariant derivatives
We summarize some identities for covariant derivative which are needed to evaluate the flow generators. We start with the commutator of two covariant derivatives acting on arbitrary tensor φ α1α2...αn , From this, one obtains where A (µ B ν) = (A µ B ν − B ν A µ )/2. In this work we assume covariantly constant background curvature and drop the termD νR ν µ λ ρ .
From the identity The tensor Φ µν αβ satisfiesḡ Then, one can define useful Laplacians, called "Lichnerowicz Laplacians" [160] acting on a spin-2 tensor field, while in a general background one has∆ and then∆ Now we see that where we have usedD µŜ µν = 0.

Expansion of projectors
To calculate the traces in the flow generators, we expand the projectors into polynomials of curvature operators. In this work, we calculate contributions up to the quadratic order of the curvature invariants, i.e.R 2 ,R µνR µν and R µνρσR µνρσ . As one can see from Eq. (B20) and (B21), the interaction parts involve at least one curvature scalar. Therefore, it is enough to have the linear order of curvature operators.
We start by expanding the inverse ofD 1 defined in Eq. (A17): Then, from Eqs. (A15) and (A16), one has where superscripts, (0) and (1), on the projectors denote the order of curvature operators, and we have defined P a . Here and hereafter, it is supposed that all projectors act on h ρσ , i.e., indices ρ and σ are contracted and thus those are not open-indices. In particular, the projector P where P µ ν is defined in Eq. (A27). The physical scalar projector is expanded aŝ from which the projector for σ given in Eq. (A29), reads Thus one has the traceless-transverse projector so that where We should note here that although the lowest order projectors P To derive the beta functions using the flow equation (A1), we need the inverse two-point functions which are second functional derivatives with respect to metric fluctuation fields. We first summarize the relations between the three couplings in the higher derivative operators for different bases. Using the fact that the Gauss-Bonnet term is topological we will only need the second functional derivative for the squared Ricci scalar curvature and the squared Ricci tensor. For these invariants we will list the explicit forms of Hessians for the physical metric fluctuations f µν and the decomposed ones t µν , σ, a µν .

Basis of gravitational interactions
For the higher derivative operators, one can write different bases. As given in Eq. (7), one of the bases for the higher derivative terms is the Weyl basis which reads On the other hand, the calculations of the beta functions using the heat kernel method yield results in the basis spanned by R 2 , R 2 µν and R 2 µνρσ rather than the operator basis in Eq. (B1). Therefore, one recasts (B1) as In order to obtain the loop contributions for the squared Weyl tensor and the Gauss-Bonnet term (8), we compare between the actions (B1) and (B2) and the read the relations between these coupling constants such that or equivalently Once the beta function for C, D and E are computed, we can obtain the beta functions for C, D and E using these relations. Let us consider the variations for the effective action (A4). The term E(ρ) R µνρσ R µνρσ is written in term of the Gauss-Bonnet term such that where we used the fact that the Gauss-Bonnet term is topological, i.e. this term can be written as a total derivative. The second term on the right-hand side does not contribute to the Hessians for a constant E(ρ). The effective action (A4) is written as where we used the relations (B3) in the second equality, and the squared Weyl tensor is given as Therefore, we do not have to calculate the variation for R µρνσ R µρνσ . It is sufficient to evaluate variations for the effect action, where we define

Variations
We list the second variations for the effective action in terms of the physical metric fluctuations f µν around a general backgroundḡ µν , i.e. for the action (B8) we show Γ grav, (2) k, The second variations for each part are computed as where the tensor Φ µν ρσ is defined in Eq. (A36).

a. Metric fluctuations
The Hessians Γ (2) k for the action (B8) is defined by For the physical metric decomposition (83), it has a simple the structure, The two-point functions of the physical metric fluctuation read Here the regulated kinetic terms for the t-mode and the physical scalar σ-mode read, respectively, The interacting parts, denoted by M are computed up to the squared order of curvature operators. For the t-mode we find while for the σ-mode one finds In Eqs. (B20) and (B21) we have defined the shorthand and N is given by Eq. (A21). The off-diagonal parts read In this work we employ the physical gauge fixing β = −1 and α → 0. For this choice, the gauge fixing action (A5) with (A7) takes the form, so that the Hessian for the physical mode f µν does not involve the gauge parameter α. The Hessian for the gauge mode is given by The propagator appearing in the flow equation involves the inverse of the (regulated ) Hessian. The Landau gauge α → 0 decouples the gauge mode in the propagator matrix. For α → 0 we do not have to specify the explicit form of Γ gravity,(2) aa . Then the Hessian for the spin-1 gauge mode a µ takes a simple form The remaining Hessian for the physical modes is given by In the limit of constant scalar fields considered in the present work the mixing terms between ϕ and metric fluctuations (t µν and σ) vanish. The Hessian (B28) becomes block-diagonal, with a separate two-point function for the scalar ϕ given by Here we define the effective derivative dependent scalar mass term with We next consider the Hessians for the ghost fields. From the ghost action (A6) with β = −1, one finds This differential operator is the same as the gauge modes in the metric, see Eq. (B27).
In order to obtain the Jacobian arising from the metric decomposition, we consider the Gaussian path integral for the metric fluctuation where E (t) is the identical matrix (A13) acting on the space satisfying the TT condition. One finds the Jacobian arising from the decomposition of the metric fluctuation, The differential operator has already appeared in the Hessians for the gauge modes (B27) in the metric and the ghost field (B32). The Jacobian for spin 1 mode can be eliminated by the redefinitions of the vector fluctuatioñ In the basisã µ , the Jacobian does not arise, while its Hessian has to be appropriately modified. In this work, we do not use this redefined field (B36). In Section E 1 e, we evaluate the flow contributions from the regulated Jacobian explicitly.
where the heat kernel trace is expanded as Note that tr (i) denotes The threshold functions Q n [W ] are given by for which the Mellin transformation yields In particular, when the flow kernel takes the form, (∆ i ) m W (∆ i ), the threshold function reads The flow kernel in this work takes typically the form, with k-dependent constants A, B and C . Here P k is the regulated momentum, i.e. P k = z + R k (z) and the regulator function in the numerator of Eq. (C7) is introduced such that∆ i is replaced to P k , namely one gives R k = A P k (z) 1+2 − z 1+2 + B P k (z) 2+2 − z 2+2 + C P k (z) 2 − z 2 for which the derivative with respect to t yields Let us now calculate the threshold functions (C5). To this end, we employ the Litim-type cutoff function [138], i.e. R k (z) = (k 2 − z)θ(k 2 − z). Thus, for n ≥ 1, one finds (1 + 2 )(n + 2 + 2 ) n + 1 + 2 a −p 1 + 2 (1 + )(n + 1 + 2 ) (1 + 2 )(n + 2 + 2 ) b a + 2 (n + 2 + 2 ) (n + 1 + 2 ) (n + 2 )(1 + 2 )(n + 2 + 2 ) c a × 1 − η n 2(n + 2 + 2 ) 2k 2n−4p−4p 2n where a = Ak −2 , B = b and c = Ck −4 . Here the anomalous dimension is defined as We define the dimensionless threshold function as 2n where nΓ(n) = Γ(n + 1) = n! is used. Note that the threshold functions for p ≥ 1 can be obtained from d 0 so that For n = 0, one has with the anomalous dimension of field Indeed, this result agrees with setting n = 0 for Eq. (C9).

Projected heat kernel
We specify the heat kernel expansion (C2) for each spin mode, We start with the symmetric spin-2-mode case: with the heat kernel coefficients, b Let us now consider the case where the projection operator P a given in Eq. (A15) is inserted, namely, tr (2L) e −s∆ T = tr (2) e −s∆ T P a . Here, the evaluation of the heat kernel coefficients is performed as follows: (C17) Once obtaining the heat kernel coefficients b (2L) i , we obtain the heat kernel expansion for the physical 2T-mode from tr (2T) e −s∆ T = tr (2) e −s∆ T P f = tr (2) e −s∆ T − tr (2) (C18) Here the first term corresponds to Eq. (C15) and the second to Eq. (C17). In the same manner, we can evaluate the heat kernel expansion for the spin 0 and TT cases, In the next subsection, we exhibit the heat kernel coefficients b The heat kernel coefficients for the Laplacian∆ T acting on the TT tensor in a general background spacetime was computed in Ref. [33], and result in b (2TT) 0 where the superscript (2TT) stands for the TT spin-2 field, χ E is the Euler characteristic associated to zero modes involved in the metric decompositions, and N is the sum of the number of Killing vectors and conformal ones. In the Weyl bases (C1) and (C2) this yields For the Laplacian∆ S acting on spin-0 scalar field, heat kernel coefficients are well-known as Next, we give a formula in order to evaluate the flow contribution from the (spin-1) measure modes, i.e., a µ , C µ , C µ and the Jacobian. As we have derived in the last section, the measure mode in the physical gauge fixing (β = −1, α → 0) takes the following uniform differential operator: For differential operators taking the typical form D 1 ν µ = δ ν µ∆V + aD µD ν −R µ ν with a constant a, the heat kernel coefficients have been evaluated [139][140][141]  For free matter fields (vector, Weyl spinor and scalar), it is convenient to define the following Laplacians, called the "Lichnerowicz Laplacians" [160],∆ These Laplacians in an Einstein spacetime obeȳ The heat kernel coefficients for the Lichnerowicz Laplacians are summarized in Table I.

Off-diagonal heat kernel expansion
The Hessians include the following operator: for which the off-diagonal heat kernel method [33] is used to evaluate the flow generator, Here we have The first two lowest order terms read ḡ αµRβν +ḡ ανRβµ +ḡ βµRαν +ḡ βνRαµ +ḡ αβRµν +ḡ µνRαβ .
For instance, one can compute the flow generators as follows: One can see from this that the lowest order term is obtained by the replacementR µνD µDν → −(R∆ S )/4 (equivalentlȳ R µν →ḡ µνR /4 or −D µDν →ḡ µν∆S /4). Indeed, at the lowest order, one has where we have used Eq. (C6), and the values of the Gamma function, Γ(3) = 2 and Γ(2) = 1.
Using the off-diagonal heat kernel expansion for a curvature tensor O µνρσ of order ofR 2 , e.g.,R µνRρσ , one can calculate where we have performed the following computations: and Using Eq. (C37), the flow generator with the TT projected operator is given by We have the projected curvature tensors by P t and P σ (orŜ µνŜρσ ) given in Eq. (A52) and Eq. (A53). Here we show explicit forms of tensor products up to of order of R 2 which appear in the flow generators. The Hessian for the TT metric fluctuation involves the following terms: We have also the following tensor products of order ofR 2 : In the Hessian for the spin-0 mode, we havê These terms involve inverse Laplacians (∆ S ) −n (n > 0). These inverse Laplacians cause an IR divergence (corresponding to zero eigenvalue of∆ S ) when performing the integration for z ∼∆ S within the flow equations. To avoid it, we will employ a field redefinition for σ as given in Eq. (E4).
In the flow kernel of σ, we calculate the following tensor-product terms of order ofR 2 : Appendix D: Summary of matter contributions We briefly summarize the matter contributions. For massless N S -scalars, N V -vector bosons and N F -Weyl fermions (corresponding to N D = N F /2 Dirac fermions), one finds Projecting it on each curvature operator and using the definition of the threshold functions (28), one obtains For U and F , their matter dependence agrees with Ref. [74] and our last paper [84]. With the relations (B4) among coupling constants, one infers the beta functions for R 2 , the squared Weyl tensor and the Gauss-Bonnet term, For the higher derivative terms, we obtain the same result in Refs. [161].
Here, the first term on the right-hand side is the contribution from the physical metric fluctuations, while the remaining terms are the measure contributions coming from the measure modes of metric fluctuations with the Jacobian and the ghost fields. As one will see in Appendix E 2, the measure contribution takes a simple form with the differential operatorD 1 defined in Eq. (A17). Below, we evaluate the different contributions to Eq. (E1).

Physical metric contribution
We evaluate the physical metric contribution in Eq. (E1), whose flow generator reads π (f ) The explicit form of J grav0,k is presented in Eq. (E32). The Hessian for the TT mode, Γ k | tt is given in Eq. (B17).

a. Full propagator and cutoff function
As we will see later, the second term in Eq. (B19) causes an IR divergence due to∆ −1 S in the denominator of S µνŜ µν . To avoid this divergence, we adopt a redefinition with a positive parameter. In the truncation which we employ in this work, the choice = 1 is enough to remove the IR divergence. We now employ the regulator functions for the TT mode and the physical scalar mode such that Laplacians in Eqs. (B18) and (B19) are replaced to P k (

. Expansion of flow generator
Let us here calculate the propagator, i.e. the inverse of Eq. (B17). One finds the propagators for each field, Then the flow generator (E3) is evaluated so that Here the flow generators for the TT mode and the spin 0 physical scalar mode are given, respectively, by The each contribution from the TT mode is calculated as follows: The lowest order term is simply evaluated by using only the heat kernel technique such that Here we define the flow kernel of the TT mode, (with = 0 in Eq. (C7)) In the next order, we need to calculate the tensor products between the projection operator and the curvature tensors. One has The next-to-next order reads The mixing term is given by where the flow kernel is defined as Let us evaluate the contribution from the spin-0 physical scalar mode. To this end, we define the flow kernel as P 2 k in the denominator of Eq. (E17) comes from the redefinition of σ as given in Eq. (E4). In order to remove divergences arising from 1/∆ S in the Hessian Γ  We next calculate the next order contributions for which we need to use the formulae (C41)-(C48). The next-to-next order contribution reads Finally we evaluate the mixing effect which reads where the flow kernel is given by (E22) e. Jacobian We have seen in Section B 3 c that the decomposition of the physical metric fluctuation (B35) and the field redefinition (E4) with = 1 yield the Jacobian as One can write the Jacobian (E23) by using auxiliary fields as det (2) Γ Jac θ , whereχ, χ are Grassmannian variables, while θ is an ordinary real variable. These field can be decomposed into the TT mode and a scalar mode, i.e. χ T = (χ TT µν , χ) and θ T = (θ TT µν , θ). We insert regulators for these fields as Jac + R Jac Jac + R Jac k θ.
Here R Jac k replaces the squared Laplacian in Eq. (E24) to P k (∆ S ). More specifically, we give We evaluate the inverse matrix Γ Jac + R Jac As in Eq. (E28), the regulator matrix R Jac k has only a finite term in the (2,2)-component (i.e. scalar-mode component), so that it is enough to take the same component in the inverse matrix, i.e. up to the squared curvature operators, one has tr (2) ∂ t R Jac k Γ Jac + R Jac k = tr (0) ∂ t R Jac0 where the product between the off-diagonal parts and the unity matrix in the (1, 1)-component has to be understood as that between the off-diagonal parts and the TT-mode projector.
Then, the flow equation from these contributions reads where the flow kernel of the Jacobian is W Jac,p [∆ S ] = ∂ t R Jac k (∆ S ) (P 2 k ) p+1 . (E33)

Measure contribution
Let us first evaluate the measure contribution without decomposing into the transverse spin-1 vector mode and spin-0 scalar mode, i.e., The contribution from the spin-1 metric fluctuation (first term on the right-hand side) has the squared differential operator, (D 1 ) 2 , as given in Eq. (B27), so that one finds where the regularized Jacobian arising from the decomposition (B35) is evaluated as The spin-1 ghost contribution are given by the same form, Thus, we see the relation k , so that the measure contribution takes a simple form, Let us choose the regulator such that the differential operatorD 1 is replaced by k 2 . Using the heat kernel method with the coefficients (C24-I) or (C24-II), one can evaluate the measure contribution, We note a difference between the present result and Ref. [84] where we have considered the Einstein-Hilbert truncation in the maximally symmetric spacetime and have derived the beta functions of U and F . In Ref. [84] we have decomposed the gauge mode a µ and the ghost field C µ into the transverse modes and the scalar modes such that {κ µ , C T µ ,C T µ }, V {u, C,C}, andD µ V T µ = 0. For these decompositions, the differential operator (C23) acts on these fields as where we have usedR µν = (R/4)ḡ µν . The differential operator∆ S in the Hessian for the spin-0 field is subtracted by the Jacobians arising from the field decomposition and the overall factor 2 in the spin-0 measure mode does not contribute in the flow generators. In Ref. [84] we have individually regularizedD 1T =∆ V −R/4 andD 0 =∆ S −R/4 and obtained measure contributions to the flow generators of the spin-1 and spin-0 measure modes, separately Thus, the total measure contribution was obtained in Ref. [84] as This result differs from Eq. (E39). This disagreement arises from an overall factor 2 in the measure contribution of the spin-0 mode: In the present work we regularize the differential operator (C23), which corresponds to regularizing 2D 0 rather thanD 0 in Eq. (E42). Hence, the result in this study can be reproduced by replacing k 2 → k 2 /2 in the spin-0 contribution η 0 in Eq. (E42), namely, This result agrees now with Eq. (E39). We conclude that the difference in results arises from different regularization procedures. The regularization of the operator (C23) seems more universal and will be adopted here.