Towards apparent convergence in asymptotically safe quantum gravity

The asymptotic safety scenario in gravity is accessed within the systematic vertex expansion scheme for functional renormalisation group flows put forward in Christiansen et al. (Phys Lett B 728:114, 2014), Christiansen et al. (Phy Rev D 93:044036, 2016), and implemented in Christiansen et al. (Phys Rev D 92:121501, 2015) for propagators and three-point functions. In the present work this expansion scheme is extended to the dynamical graviton four-point function. For the first time, this provides us with a closed flow equation for the graviton propagator: all vertices and propagators involved are computed from their own flows. In terms of a covariant operator expansion the current approximation gives access to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Lambda $$\end{document}Λ, R, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R^2$$\end{document}R2 as well as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_{\mu \nu }^2$$\end{document}Rμν2 and higher derivative operators. We find a UV fixed point with three attractive and two repulsive directions, thus confirming previous studies on the relevance of the first three operators. In the infrared we find trajectories that correspond to classical general relativity and further show non-classical behaviour in some fluctuation couplings. We also find signatures for the apparent convergence of the systematic vertex expansion. This opens a promising path towards establishing asymptotically safe gravity in terms of apparent convergence.

In [3] the systematic vertex expansion in quantum gravity initiated in [1,2], see also [26,33], has been pushed to the graviton three-point function, for the first time including a dynamical graviton-scattering in the asymptotic safety analysis. Here we extend the vertex expansion to the graviton four-point function. Apart from the significant technical challenge such an upgrade of the approximation has posed, we think that this constitutes a necessary and significant progress towards asymptotically safe gravity: • As such it is an important step towards apparent convergence of the vertex expansion in quantum gravity: apparent convergence aims at the convergence of vertices as well as observables in the order of a given systematic expansion scheme; here we use the vertex expansion scheme. Together with the investigation of the regulator (in-)dependence of observables this provides a systematic error estimate in the present approach, and should be compared with apparent continuum scaling and extrapolation on the lattice. • The present approximation allows for the identification of diffeomorphism-invariant structures in the vertex expansion, i.e. R 2 and R 2 μν tensor structures as well as those of higher derivative invariants. This is not only important for getting access to the number of relevant directions at the asymptotically safe UV fixed point in quantum gravity, but can also provide non-trivial support and additional information for computations within the standard background field approximation.
• It is the first approximation in the asymptotic safety approach to gravity where the flow of the pivotal building block, the two-point correlation function or (inverse) propagator, is closed: The flows of all involved vertex functions are computed within given approximations. As the propagator is the core object in the present approach, we consider this an important milestone on the way towards asymptotically safe quantum gravity. As a main result of this paper, we find further significant evidence for a non-trivial UV fixed point in quantum gravity. This fixed point has three relevant directions and two repulsive ones. The three relevant directions can be associated with the cosmological constant (graviton mass parameter), Newton's coupling and the R 2 -coupling, see Sect. 5. We also investigate the stability of this UV fixed point and observe that the system is significantly less sensitive to the closure of the flow equations than previous truncations. In addition, we observe that the critical exponents also become less sensitive to the details of the approximations. These are two necessary signatures of apparent convergence. Furthermore, we investigate the infrared (IR) behaviour of the system and find trajectories connecting the UV fixed point with classical general relativity in the IR.
This work is structured in the following way: In Sect. 2 we elaborate on the important property of background independence of quantum gravity and its consequences for the dynamical correlation functions. In Sect. 3, we briefly introduce the functional renormalisation group and the covariant vertex expansion used in the present work. Next, Sect. 4 presents our setup and the derivations of the flow equations for the couplings. We especially focus on how the tensor structures of higher curvature invariants are embedded in our vertex expansion. In Sect. 5, we present our results, in particular the non-trivial UV fixed point with three attractive and two repulsive directions. Furthermore, we test the stability of the UV fixed point with respect to change of the precise truncation, and find that in almost all approximations we obtain a similar UV fixed point. In Sect. 6, we investigate the IR-behaviour of the UV-finite trajectories. Our analysis confirms the classical IR fixed points found before in the vertex expansion. In Sect. 7 we discuss the convergence of the present vertex expansion scheme with an increasing number of flowing n-point correlation functions. Finally, in Sect. 8 we summarise our results.

Diffeomorphism invariance and background independence
The quantum field theoretical formulation of quantum gravity in terms of metric correlation functions necessitates the introduction of a background metricḡ, and quantum fluctuations are taken to be fluctuations about this background metric. This begs the question of whether diffeomorphism invariance and background independence of observables are guaranteed in such a framework. While this is an important question, its answer is not directly relevant for the computations presented here. Hence, the present chapter may be skipped in a first reading.
In the present work we perform computations in the linear split, where the full metric g is given by g =ḡ+h. More general splits, g =ḡ + f (ḡ, h) have been considered for example within the geometrical or Vilkovisky-deWitt approach, e.g. [26,[83][84][85], or the exponential split, e.g. [40,44,68,[86][87][88]. Observables, on the other hand, are background independent. This property is encoded in the Nielsen (NI) or split-Ward identities (SWI) that relate derivatives of the effective action [ḡ, φ] w.r.t. the background metricḡ to those w.r.t. the graviton fluctuations h. Here we have introduced the fluctuation superfield φ = (h, c,c). (1) The (anti-) ghost fields, c andc, stem from the Faddeev-Popov gauge fixing procedure, see next section. The effective action generates all one-particle-irreducible correlation functions and as such encodes the symmetries of the theory. Schematically, these identities read [6,20,26,84,[89][90][91][92][93][94][95] where we have suppressed the ghost fields. In the linear split we have C[ḡ, h](x, y) = δ(x − y) and the second term N [ḡ, h] carries the information about the nontrivial behaviour under diffeomorphism transformations of the gauge fixing sector and the regularisation. In turn, in the geometrical approach diffeomorphism invariance of the effective action is achieved by a non-linear split with f (ḡ, h) leading to the non-trivial prefactor C [ḡ, h] in (2). The term N [ḡ, h] then carries the deformation of the Nielsen identity in the presence of a regularisation but does not spoil diffeomorphism invariance. In both cases the Nielsen identity is a combination of a quantum equation of motion, the Dyson-Schwinger equation, and the Slavnov-Taylor identity (STI) or diffeomorphism constraint. The setup also entails that correlations of the fluctuation fields are necessarily background-dependent. This is easily seen by iterating (2). Moreover, in the linear split, diffeomorphism invariance of the observables is encoded in non-trivial STIs for the fluctuation correlation functions, while in the geometrical formulation, the nontrivial STIs are encoded in expectation values of f (ḡ, h) and its derivatives.
Due to (2) we have to deal with the peculiarity that background independence and physical diffeomorphism invariance of observables necessitate background-dependence and non-trivial STIs for the correlation functions of the fluctuation fields. This leads to seemingly self-contradictory statements: in particular, for the quantum effective action [ḡ, h] it entails that physical diffeomorphism invariance of observables is not achieved by diffeomorphism invariance w.r.t. diffeomorphism transformations of the fluctuation fields. The latter does not do justice to either diffeomorphism invariance or background independence.
This peculiarity can easily be checked in a non-Abelian gauge theory within the background field formulation: in a fluctuation gauge invariant approximation to the effective action, even two-loop universal observables such as the twoloop β-function cannot be computed correctly. Indeed, in this case it is well-known that only the non-trivial STIs for the fluctuation gauge field elevate the auxiliary background gauge invariance to the physical one holding for observables, see e.g. [96,97].
The above considerations underline the importance of a direct computation of correlation functions of the fluctuation field h. Indeed, the corresponding set of flow equations for (n) = (0,n) is closed in the sense that the flow diagrams only depend on (n) with n ≥ 2. Here, (n,m) stands for the n-th background field derivative and m-th fluctuation field derivative of the effective action, In turn, the flows for pure background, (n,0) , or mixed background-fluctuation functions, (n,m) with m = 0, necessitate the fluctuation correlation functions as an input: the background correlation functions can be iteratively computed in powers of the background metric. In other words, the dynamics of the system is solely determined by the pure fluctuation correlation functions. For this reason, several approaches for computing correlation functions of the fluctuation field have been put forward in the last years. Some of these approaches are set up for computing both correlation functions of the fluctuation field as well as those of the background field: vertex expansion [1][2][3]33,71,73,74] and bimetric approach [20,24,37,90]. Another set of approaches relies on utilising the Nielsen or split Ward identities explicitly or implicitly, [26,84,[92][93][94][95][98][99][100][101][102].
So far, the (n) for n ≥ 2 have been computed directly in only one approach, the vertex expansion, see [1][2][3]33] for pure gravity and [71,73,74] for gravity-matter systems. A mixture of vertex expansion and background approximation has been used in [33,67,72]. Present results include (n) for n = 0, 1, 2, 3, where higher vertices have been estimated by lower ones, relying on approximate covariance of the correlation functions. Such a structure has already been confirmed in the perturbative and semi-perturbative regime of QCD, see [103].

Effective action and functional renormalisation group
The set of (covariant) correlation functions of the metric, g(x 1 ) · · · g(x n ) , defines a given theory of quantum gravity. All observables can be constructed from these basic building blocks. The correlation functions are generated from the single metric effective action, [g] = [g, h = 0], which is the free energy in a given metric background g =ḡ + h at h = 0. Here we have restricted ourselves to a linear split. The underlying classical action is the Einstein-Hilbert action, We employ a linear, de-Donder type gauge-fixing, In particular we use the harmonic gauge given by β = 1. This choice simplifies computations considerably due to the fact that the poles of all modes of the classical graviton propagator coincide. We furthermore work in the limit of a vanishing gauge parameter, α → 0. This is favourable because then the gauge does not change during the flow since α = 0 is a RG fixed point [104]. Moreover, it allows for a clear separation of propagating and non-propagating degrees of freedom. The ghost part of the action reads wherec and c denote the (anti-) ghost field and M is the Faddeev-Popov operator deduced from (6). For β = 1 it is given by The gauge fixing and ghost term in (5) and (7) introduce the separate dependence onḡ and h leading to the non-trivial Nielsen identities in (2).

Flow equation
An efficient way of computing non-perturbative correlation functions is the functional renormalisation group. In its form for the effective action, see [5,105,106], it has been applied to quantum gravity [6]. For reviews on the FRG approach to gauge theories and gravity see e.g. [53][54][55][56]89,107]. The RG flow of the effective action for pure quantum gravity is given by Here, ∂ t ≡ k∂ k denotes the scale derivative, where k is the infrared cutoff scale.
(2) k = (0,2) k stands for the second fluctuation field derivative of the effective action, while R k is the regulator, which suppresses momenta below k. The trace in (9) sums over internal indices and integrates over space-time.
The introduction of cutoff terms leads to regulatordependent modifications of STIs and NIs that vanish for R k → 0. The respective symmetry identities have hence been named modified Slavnov-Taylor identities (mSTIs) and modified Nielsen-or split Ward identities (mNIs/mSWIs). The modification entails the breaking of the physical or quantum diffeomorphism invariance in the presence of a background covariant momentum cutoff. Still, background diffeomorphism invariance is maintained in the presence of the cutoff term.

Covariant expansion
The effective action k [ḡ, φ] depends on the background metricḡ and the fluctuation superfield φ = (h, c,c), see (1), separately. The functional flow equation (9) is accompanied by the functional mSTIs & mNIs for the effective action that monitor the breaking of quantum diffeomorphism invariance, see (2) in Sect. 2.
In order to solve (9), we employ a vertex expansion around a given backgroundḡ, to wit where the superscript fields in parentheses are a short-hand notation for field derivatives, and where contracting over super-indices i j occurring twice is implied. In this work, we include the full flow of the vertex functions up to the graviton four-point function.
As discussed in Sect. 2, the expansion coefficients satisfy mSTIs as well as mNIs with (n,m) k being defined in (3). For the sake of simplicity we now restrict ourselves to the gauge fixing used in the present work, (6) with α = 0. Then the fluctuation graviton propagator is transverse: it is annihilated by the gauge fixing condition.
An important feature of the functional RG equations is that for α = 0 the flow equations for the transverse vertices In other words, the system of transverse fluctuation correlation functions is closed and determines the dynamics of the system. On the other hand, the mSTIs are non-trivial relations for the longitudinal parts of vertices in terms of transverse vertices and longitudinal ones. This leads us to the schematic relation see [108] for non-Abelian gauge theories. In consequence, the mSTIs provide no direct information about the transverse correlation functions without further constraint. In the perturbative regime this additional constraint is given by the uniformity of the vertices, for a detailed discussion in non-Abelian gauge theories see [109]. Accordingly, our task reduces to the evaluation of the coupled set of flow equations for the transverse vertices The square of the Weyl tensor C 2 is eliminated via the Gauß-Bonnet term, which is a topological invariant. Higherderivative terms, such as are also taken into account. Without the constraint f (0) = 0, Eq. (14) also includes R 2 and R 2 μν , more details on this basis can be found in Sect. 4. Note that also non-diffeomorphisminvariant terms are generated by the flow. In Sect. 4.3 we discuss all invariants which are included in the parameterisation of our vertices.
For the background vertices (n,0) k we use the following: the NIs become trivial in the IR as we approach classical gravity, as shown in Sect. 6. Moreover, for one of the two IR fixed points this implies that the derivative with respect to a background field is the same as a derivative with respect to a fluctuation field. This allows us to impose the trivial NIs in the IR, and all couplings are related. Then, the couplings at k > 0 follow from the flow equation. However, for the fluctuation couplings this amounts to solving a fine-tuning problem in the UV, for more details see Sect. 6. The latter is deferred to future work.

Flows of correlation functions
In this chapter we discuss the technical details of the covariant expansion scheme used in the present work, including the approximations used and their legitimisation. In our opinion, a careful reading of this chapter is essential for a full understanding of the results obtained in the present work. This applies in particular to Sect. 4.1.

Covariant tensors and uniformity
The flows of the n-point correlation functions are generated from the FRG Eq. (9) by taking n-th order fluctuation field derivatives in a backgroundḡ, (see Fig. 1). In order to solve the flow equation, we employ a vertex ansatz [1,110] including the flow of all relevant vertices up to the graviton fourpoint function. This vertex ansatz disentangles the couplings of background and fluctuation fields by introducing individual couplings n and G n for each n-point function. These individual couplings are introduced at the level of the npoint correlators and replace the cosmological constant and Newton's coupling G N of the classical Einstein-Hilbert action after performing the respective field derivatives. In summary, for the flat backgroundḡ = δ our vertex ansatz reads where denote the tensor structures extracted from the classical gauge-fixed Einstein-Hilbert action (4). The only flowing parameter in these tensors T (φ 1 ...φ n ) is n , while G n ( p) carries the global scale-and momentum dependence of the ver- tex. In the above equations, p = ( p φ 1 , ..., p φ n ) denotes the momenta of the external fields φ i of the vertex. Apart from their flow equations, the n-point functions in (15) also satisfy standard RG-equations, see e.g. [89]. These RG-equations entail the reparameterisation invariance of the theory under a complete rescaling of all scales including k. With the parameterisation given in (15), this RG-running is completely carried by the wave function renormalisations Z φ i ( p 2 i ) of the fields φ i , see e.g. [2,71,110]. Consequently, the G n and n are RG-invariant, and hence are more directly related to observables such as S-matrix elements. This parameterisation of the vertices also ensures that the wave function renormalisations never appear directly in the flow equations, but only via the anomalous dimensions G n ( p) is the gravitational coupling of the n-point function, while n denotes the momentum-independent part of the correlation function. In particular, 2 is related to the graviton mass parameter M 2 := −2 2 . Finally, all the parameters Z φ i , G n , and n are scale-dependent, but we have dropped the subscript k in order to improve readability.
In principle, all tensor structures, including nondiffeomorphism-invariant ones, are generated by the flow, but for our vertex functions we choose to concentrate on the classical Einstein-Hilbert tensor structures in the presence of a non-vanishing cosmological constant. Despite the restriction to these tensor structures, the n-point functions have an overlap with higher curvature invariants via the momentum dependence of the gravitational couplings. For example, the complete set of invariants that span the graviton wave function renormalisation is given by where the superscript indicates that is a covariant tensor contributing to the two-point correlation function. Note also that we now drop the restriction on f present in (14). Then, this invariant naturally includes R 2 and R 2 μν as the lowest order local terms. If we also allow for general momentumdependencies, the corresponding covariant functions f are given by given by The lowest order local terms, R 2 and R 2 μν , are given by P (2) R 2 = 1 and P (2) R 2 μν = 1, respectively. Note that (19) also allows for non-local terms in the IR, i.e. anomaly-driven terms with P (2) R 2 = 1/∇ 2 , see e.g. [111]. In turn, higher curvature invariants do not belong to the set of the graviton wave function renormalisation since they are at least cubic in the graviton fluctuation field.
In the present work we resort to a uniform graviton propagator in order to limit the already large computer-algebraic effort involved. The uniform wave function renormalisation is then set to be that of the combinatorially dominant tensor structure, the transverse-traceless graviton wave function renormalisation, thereby estimating the wave function renormalisations of the other modes by the transverse-traceless one. Such uniform approximations have been very successfully used in thermal field theory. There, usually the tensor structures transverse to the heat-bath are used as the uniform tensor structure, for a detailed discussion see e.g. [112] and references therein. This approximation is typically supported by combinatorial dominance of this tensor structure in the flow diagrams. Indeed, as already indicated above, the transverse-traceless mode gives the combinatorially largest contribution to the flow of the vertices computed here. Note that such an approximation would get further support if the R tensor structures dominate the flows, which indeed happens in the present computation.
Within this approximation the R 2 tensor structures drop out on the left-hand side of the graviton flow, since R is already quadratic in the transverse-traceless graviton fluctuation field: in other words, the tensors defined by f (2) R 2 in (19) have no overlap with the transverse-traceless graviton.
The set of invariants that span the gravitational coupling Again, the invariants R 2 and R 3 can be excluded from this set due to their order in transverse-traceless graviton fluctuation fields. In consequence, G 4 ( p) is the only coupling in our setup that has overlap with R 2 contributions and higher terms in f (4) R 2 . Furthermore, in Sect. 4.3 we show that the by far dominant contribution to G 3 ( p) in the momentum range 0 ≤ p 2 ≤ k 2 stems from the invariant R. All higher momentum dependencies of the graviton three-point function are covered by the momentum dependence of the graviton wave function renormalisation. This was already observed in [3]. As already briefly mentioned above, it gives further support to the current uniform approximation: the assumption of uniformity allows us to restrict ourselves to computing the Einstein-Hilbert tensor structure for the transverse-traceless graviton as the combinatorially dominating tensor structure. The striking momentum-independence of the actual numerical flows supports a momentum-independent approximation of G 3 . In terms of (19) it implies that the dominant tensor structure for the transverse-traceless mode is given by f μν tensor structure vanishes approximately, see (27). In contrast to the situation for the two-and three point function, the R 2 invariant overlaps with our transversetraceless projection for the graviton four-point function. Indeed, its flow receives significant contributions from the invariant R 2 . It follows that for the graviton four-point function R is not the only dominant invariant in the momentum range from p = 0 to p = k, as we show in Sect. 4.3. In consequence we either have to disentangle contributions from R and R 2 tensor structures in terms of an additional tensor structure or we resolve the momentum dependence of G 4 ( p). In the present work we follow the latter procedure, see Sect. 4.5 for details.

Projection onto n-point functions
The flow equations for the couplings n and G n are obtained by the following projection onto the flow of the graviton npoint functions ∂ t (n) k . We use the classical Einstein-Hilbert tensor structures T (n) ( p; n ) as a basis for our projec-tion operators. Furthermore, we project onto the spin-two transverse-traceless part of the flow, which is numerically dominant. Moreover, it does not depend on the gauge. This transverse-traceless projection operator is then applied to all external graviton legs. The flow of the couplings n is then extracted with the help of the momentum-independent part of said tensor structures, namely n := T (n) (0; n )/ n . For the couplings G n we use G n := T (n) ( p; 0)/ p 2 . Dividing by n and p 2 ensures that the projection operators are dimensionless and scale-independent.
In principle, the flow of any n-point function depends on all external momenta p i , i ∈ {1, ..., n}, where e.g. p n can be eliminated due to momentum conservation. For the two-point function, the momentum configuration is trivial, and only one momentum squared, p 2 , needs to be taken into account. In contrast, this dependence becomes increasingly complex for the higher n-point functions: The three-point function depends on three parameters (two momenta squared and one angle), the four-point function already depends on six parameters, and so on. To simplify the computations, we use a maximally symmetric (n − 1)-simplex configuration for all n-point-functions, thereby reducing the momentum dependence to a single scalar parameter. This symmetric momentum configuration was already used for the graviton threepoint function in [3]. In the context of Yang-Mills theories, this approximation has been shown to be in good agreement with lattice computations on the level of the flow of the propagator [109]. Notably, in the symmetric momentum configuration all external momenta have the same absolute value p, and the same angles between each other. The scalar product of any two momenta in this momentum configuration then reads where δ i j denotes the Kronecker delta. Note that such a symmetric momentum configuration only exists up to the (d +1)point function, where d is the dimension of spacetime. In the following, the expressions Flow (n) stand for the dimensionless right-hand sides of the flow equations divided by appropriate powers of the wave function renormalisations. More explicitly, we define where the index i represents the projection on some tensor structure. In this work, we use the transverse-traceless projection operator TT , the projection operators G n and n mentioned earlier for the graviton n-point functions, as well as the transverse projection operator T for the ghost propagator. Note that the objects Flow Last but not least, we choose to model the regulator functions R φ i on the corresponding two-point functions at vanishing mass, i.e.
Here, r φ i ( p 2 i /k 2 ) denotes the regulator shape function. For all fields in this work, we choose the Litim regulator [113], to wit This choice allows for analytic flow equations for all couplings that are evaluated at vanishing external momenta. Furthermore, we introduce the dimensionless couplings At the UV and IR fixed points, the flow of these dimensionless couplings vanishes.

Momentum dependence of the graviton n-point functions
We now investigate the momentum dependence of the flow of the graviton n-point functions as defined in (22). We restrict ourselves to the momentum range 0 ≤ p 2 ≤ k 2 as well as to the transverse-traceless part of the graviton n-point functions.
The first non-trivial result is that the flows of the graviton three-and four-point functions projected on the tensor structure of the gravitational coupling and divided by (− n 2 η h ( p 2 ) − n + 2) are well described by a polynomial in p 2 , provided that the couplings λ n are small, i.e.
with some constants a i and b i that depend on the evaluation point in theory space. This momentum dependence is displayed in Fig. 2. We emphasise that these equations only hold in the momentum range 0 ≤ p 2 ≤ k 2 , if the flow is generated by Einstein-Hilbert vertices, and if the constant parts of the vertices are small, i.e. |λ n | 1. If the condition of small λ n is violated, then the flow as in (26) is non-polynomial. We did not compute the flow generated by an action including higher curvature terms, however, we suspect that the flow will still be polynomial but possibly of a higher degree.   (26). The flows are evaluated at (μ, λ 3 , λ 4 , g 3 , g 4 ) = (−0.4, 0.1, −0.1, 0.7, 0.5) and λ 6 = λ 5 = λ 3 as well as g 6 = g 5 = g 4 . The flows have such a simple polynomial structure as long as all couplings λ n remain small, i.e. |λ n | 1. Impor-tantly, the inclusion of a p 4 term in the left panel offers no significant improvement. Note that the constant parts of the functions are irrelevant for the beta functions since they are extracted from a different tensor projection. For p 2 > k 2 the momentum dependence of the flows is not polynomial anymore It is important to note that the graviton three-and fourpoint functions have a different highest power in p 2 . This is a second non-trivial result for the following reasons: as already mentioned before, the coupling g 3 ( p 2 ) has an overlap with R and R 2 μν , and higher derivative terms in f ifest itself in a p 4 -contribution to the flow of the graviton three-point function. Equation (26) and Fig. 2 show that such a p 4 -contribution as well as higher ones are approximately vanishing. This demonstrates in particular that the generation of R 2 μν is non-trivially suppressed. In other words, where the superscript indicates the three-graviton vertex.
On the other hand, the projection on g 4 ( p 2 ) overlaps with R, R 2 μν , R 2 tensor structures, and the related higher derivatives terms in f It also overlaps with curvature invariants to the third power with covariant tensors such as f (4) R 3 μν and similar ones. Note that it has no overlap with f (4) R 3 . Similarly to possible p 4 -contributions for the threegraviton vertex, p 6 -contributions and even higher powers in p 2 could be generated but are non-trivially suppressed. The p 4 -contribution to the flow, which is described in (26) and displayed in Fig. 2, could stem from either R 2 or R 2 μν tensor structures. Now we use (27). It entails that the graviton threepoint vertex does not generate the diffeomorphism invariant term R 2 μν although it has an overlap with it. This excludes R 2 μν as a relevant UV direction, which would otherwise be generated in all vertices. This statement only holds if we exclude non-trivial cancellations of which we have not seen any signature. Accordingly we set and conclude that this p 4 -contribution or at least its UVrelevant part stems solely from R 2 . It may be used to determine f (4) R 2 . In summary, the above statements about the momentumdependencies are highly non-trivial and show that R 2contributions are generated while R 2 μν and other higher derivative terms are strongly suppressed. These non-trivial findings also allow us to determine the most efficient way to project precisely onto the couplings of different invariants. This is discussed in Sect. 4.5.
We close this section with a brief discussion of the effect of higher derivative terms on perturbative renormalisability and the potential generation of massive ghost states. As already discussed in [114] in a perturbative setup, it is precisely the R 2 μν term which makes the theory perturbatively renormalisable. However, in this setup it gives rise to negative norm states. On the other hand, the R 2 term neither ensures perturbative renormalisability, nor does it generate negative norm states. This is linked to the fact that the R 2 term does not contribute to the transverse traceless part of the graviton propagator. Consequently, the non-trivial suppression of R 2 μν tensor structures might be interpreted as a hint that we do not suffer from massive ghost states. However, a fully conclusive investigation requires the access to the pole structure of the graviton propagator, and hence a Wick rotation. Progress in the direction of real-time flows in general theories and gravity has been made e.g. in [25,50,[115][116][117][118][119][120][121].

Higher order vertices and the background effective action
The results in the last section immediately lead to the question about the importance of the higher-order covariant ten-sor structures like e.g. f R n which have no overlap with the graviton n-point functions computed in this work. These are potentially relevant for the flows of G 5 and G 6 . These tensors have been dropped in the current work, thus closing our vertex expansion. However, we may utilise previous results obtained within the background field approximation for estimating their importance: first we note that R 2 gives rise to a new relevant direction, as we will show in Sect. 5.1. This has also been observed for the background field approximation [10,[14][15][16]31]. There it has also been shown that the critical dimensions of the R n -terms approximately follow their canonical counting [31]. Furthermore, our results so far have sustained the qualitative reliability of the background field approximation for all but the most relevant couplings. Indeed, it is the background field-dependence of the regulator that dominates the deviation of the background approximation from the full analysis for the low order vertices, and in particular the mass parameter μ of the graviton. This fielddependence is less relevant for the higher order terms. Thus, we may qualitatively trust the background field approximation for higher curvature terms. This means that they are of sub-leading importance and can be dropped accordingly. Finally, the above findings together with those from the literature suggest that an Einstein-Hilbert action is generating a diffeomorphism-invariant R 2 -term but not an R 2 μν term in the diffeomorphism-invariant background effective action . Moreover, no higher derivative terms are generated if a non-trivial wave function renormalisation Z h ( p 2 ) and graviton mass parameter μ = −2λ 2 are taken into account. Note that this only applies for an expansion with p 2 < k 2 . This is a very interesting finding as it provides strong non-trivial support for the semi-quantitative reliability of the background approximation in terms of an expansion in R for spectral values smaller than k 2 subject to a resolution of the fluctuating graviton propagator: μ and Z h have to be determined from the flows of the fluctuation fields or in terms of the mNIs.

Flow equations for the couplings
In this section we derive the flow equations for the couplings from the projected n-point functions.
The flow equations for μ and η h ( p 2 ) are extracted from the transverse-traceless part of the flow of the graviton twopoint function. We evaluate this two-point function at p 2 = 0 for ∂ t μ, and bilocally at −μk 2 and p 2 for η h ( p 2 ). The algebraic equation for η c ( p 2 ) can be obtained directly from the transverse part of the flow of the ghost two-point function. These equations are derived in the same fashion as in [2,71]. For details see also App. 1.
In the case of the couplings λ n and g n ( p 2 ), we project onto the flow of the graviton n-point functions. The flow equations for the couplings λ n are always obtained at p 2 = 0, since λ n describes the momentum-independent part of the graviton n-point functions.
In the case of the couplings g n ( p 2 ) it is technically challenging to resolve the full momentum dependence in the flow. Thus, we resort to a further approximation of the momentumdependence. We have checked that this approximation holds quantitatively. First we note that typically FRG-flows are strongly peaked at q ≈ k due to the factor q 3 from the loop integration and the decay for momenta q k due to ∂ t R k (q 2 ). This certainly holds for all the flows considered here. From this we can infer that we extract the leading contribution to the flow diagrams if we feed g n (k 2 ) back into the diagrams. In consequence we compute only the flow equations for g n (k 2 ), as they form a closed system of equations within the given approximation.
Conveniently, the momentum dependence of the flow for g 3 ( p 2 ) is trivial, see Fig. 2 in Sect. 4.3. Hence the approximation discussed above is of quantitative nature, and we obtain precisely the same equation as in [3].
In contrast, the flow of the graviton four-point function exhibits a p 4 contribution, implying a non-trivial g 4 ( p 2 ). Here our approximation is necessary to simplify the computation significantly. The flow equation for g 4 (k 2 ) is obtained from a bilocal momentum projection at p 2 = 0 and p 2 = k 2 , and furthermore uses an approximation that relies on the fact that the coupling λ 4 remains small. We refer to this equation as a bilocal equation. It is explicitly displayed in Appendix 1, see (E5). Within our setup this equation gives the best approximation of the vertex flows since it feeds back the most important momentum information into the flow. This further entails that the coupling g 4 (k 2 ) includes information about the invariants R and R 2 .

Disentangling R and R 2 tensor structures
In this section we present projection operators that disentangle contributions from R and R 2 tensor structures to the flows of the couplings g n ( p 2 ). In the present setup this only allows us to switch off the R 2 coupling and thus to check the importance of the R 2 coupling.
For the disentanglement, we have to pay attention to two things: First of all, a local momentum projection at p 2 = 0 is very sensitive to small fluctuations and in consequence not very precise with regard to the whole momentum range 0 ≤ p 2 ≤ k 2 . This was already discussed in [2,3] and is explicitly shown in Appendix 1. Hence, we have to rely on non-local momentum projections. Here the highest polynomial power of p 2 , as indicated in (26), dictates the simplest way of projecting on the p 2 -coefficient. The graviton threepoint function is at most quadratic in the external momentum, and consequently it is enough to use a bilocal projection at p 2 = 0 and p 2 = k 2 . The resulting equation is displayed in Appendix 1, see (E6).
The graviton four-point function, on the other hand, has p 4 as its highest momentum power, i.e. it is of the form see also (26). Thus a bilocal momentum projection would not extract the p 2 coefficient b 1 alone. Instead, we use a trilocal momentum projection at p 2 = 0, p 2 = k 2 /2, and p 2 = k 2 in order to solve the above equation for b 1 . I.e., we solve a system of linear equations and obtain The resulting flow equation is again presented in Appendix 1, see (E7). For even higher order momentum contributions we would have to use even more points of evaluation. These momentum projections together with the observation of (26) guarantee that we project precisely on the p 2 coefficient in the whole momentum range 0 ≤ p 2 ≤ k 2 .
A natural upgrade of the current approximations amounts to the introduction of a second tensor structure that is orthogonal to the Einstein-Hilbert one in terms of these projections. Within our uniformity assumption this is considered to be sub-leading, and the momentum-dependence of g 4 ( p 2 ) takes care of the contribution of the R 2 tensor structure f (4) R 2 . While the orthogonal projection on the respective flow is simple, its back-feeding demands a two tensor structure approximation of the three-and four-graviton vertex in the flow, the implementation of which is deferred to future work.
Here, we only perform a further check of the relevance of the R 2 tensor structure. This sustains the fact that the inclusion of the four-graviton vertex with its contribution of the R 2 tensor structure leads to an additional UV-relevant direction. To that end we generalise our ansatz for the graviton fourpoint function such that we can extract a flow equation for both the Einstein-Hilbert tensor structure as well as for the R 2 tensor structure. As already mentioned above, we cannot feed the generated coupling back into the flows, since they are given by vertices with Einstein-Hilbert tensor structures. Instead we compute the fixed point value that arises only from the Einstein-Hilbert tensor structures.
As the corresponding ansatz for the transverse-traceless graviton four-point function we choose (4) which is precisely the vertex that emerges from the sum of Einstein-Hilbert tensor structure and R 2 tensor structure. The related generating diffeomorphism-invariant action for this four-graviton vertex is where S EH is defined as in (4). The flow of 4 is then obtained by the trilocal momentum projection described below (29).
The explicit form of the resulting flow equation for the dimensionless coupling ω 4 := 4 k 2 is given in Appendix 1, see (E8). Note that in the present approximation, the flows do not depend on the coupling ω 4 since it does not feed back into the vertices.

Computational details
The  [122,123] and Mathematica [124]. For individual tasks, we employ specialised and in part self-developed Mathematica packages. In particular, we use VertEXpand [125] and xPert [126] for the generation of vertex functions, DoFun [127] to obtain symbolic flow equations, and the FormTracer [103,109,128] to create optimised FORM scripts to trace diagrams. Furthermore, we make use of the self-developed framework TARDIS [125] facilitating an automated and seamless usage of the aforementioned tools.

Asymptotic safety
In this section, we discuss the UV fixed point structure of our system. We first present our best result, which includes the tensor structures as presented in Sect. 4.1 and in particular in (18) and (20). The underlying UV-relevant diffeomorphism invariants turn out to be , R, and R 2 . The R 2 coupling is included via the momentum dependence of the gravitational coupling g 4 ( p 2 ), see Sect. 4.5. As a main result we find an attractive UV fixed point with three attractive directions. The third attractive direction is related to the inclusion of the R 2 coupling.
We further analyse the stability of this UV fixed point with respect to the identification of the higher couplings. We also analyse the previous truncation [3] and compare the stability of both truncations. Here we find that the improvement of the truncation increases the stability of the system. In particular, we find a rather large area in the theory space of higher couplings where the UV fixed point exists with three attractive directions throughout.
Lastly, we discuss the importance of the R 2 coupling. In Sect. 4.6 we have constructed projection operators that disentangle the contributions from R and R 2 tensor structures. This allows us to switch off the R 2 coupling and compare the stability of the reduced system to that of the full system. We find that the reduced system is significantly less stable, and that the area in the theory space of higher couplings where the fixed point exists is rather small. This highlights the importance of the R 2 coupling.

UV fixed point
In this section we display the UV fixed point structure of our full system. This means that we feed back the generated R 2 coupling via the momentum dependence of the gravitational coupling g 4 ( p 2 ), as discussed in Sect. 4.5. Fixed points are by definition points where the flows of the dimensionless couplings vanish. In consequence, we look for the roots of the Eqs. (E1), (E4), (E6), and (E5). We use the identification scheme g 6 = g 5 = g 4 and λ 6 = λ 5 = λ 3 . We find a UV fixed point at the values μ * , λ * 3 , λ * 4 , g * 3 , g * 4 = (−0.45, 0.12, 0.028, 0.83, 0.57) .
The fixed point values are similar to those of the previous truncation [3]. The biggest change concerns the graviton mass parameter, which is now less negative and thus further away from its pole. Moreover, it is remarkable that the new couplings λ 4 and g 4 are close to their lower counterparts λ 3 and g 3 , but not at precisely the same values. Since we use the difference between these couplings to parameterise the breaking of diffeomorphism invariance, this is more or less what we expected. This issue is further discussed in the next section. We do not have access to the full stability matrix of the UV fixed point due to the unknown flow equations of the higher couplings. For this reason, we discuss two different approximations of the stability matrix. The main difference between these two approximations concerns the order of taking the derivatives and identifying the higher couplings, which is explained in more detail in Appendix 1. We argue that in a well converged approximation scheme the most relevant critical exponents should not depend on the approximation of the stability matrix. Thus, we can use the two different approximations to judge the quality of the current level of truncation. In this work, we define the critical exponents as the eigenvalues of the stability matrix without a minus sign. We call the critical exponents of the first approximationθ i , and the ones of the second approximationθ i . The critical exponents using the first approximation are given bȳ while the critical exponents using the second approximation arẽ θ i = (−5.0, −0.37 ± 2.4i, 5.6, 7.9).
Hence this fixed point has three attractive directions in both approximations of the stability matrix. The third attractive direction compared to the system of the graviton three-point function [3] is related to the fact that the graviton four-point function has an overlap with R 2 , which we feed back via the momentum dependence of the gravitational coupling g 4 ( p 2 ). The R 2 coupling has also been relevant in earlier computations with the background field approximation [10,[14][15][16]31]. In addition, note that the most attractive eigenvalue is almost identical in both approximations of the stability matrix. This is a positive sign towards convergence since it is expected that the lowest eigenvalue is the first that converges, c.f. Appendix 1. Furthermore, the anomalous dimensions at the UV fixed point read where we have chosen to display the anomalous dimensions at the momenta that feed back into the flow. All anomalous dimensions stay well below the reliability bound η φ i ( p 2 ) < 2, as introduced in [71].

Stability
In the following we investigate the UV fixed point from the previous section by varying the identification of the higher couplings. Again we look for the roots of the Eqs. (E1), (E4), (E6), and (E5). These equations however still depend on the higher couplings g 5 , g 6 , λ 5 , and λ 6 . We have to identify these couplings with the lower ones or set them to constants in order to close the flow equations. It is a natural choice to simply set these higher couplings equal to lower ones, e.g. g 6 = g 5 = g 3 and λ 6 = λ 5 = λ 3 , as done in the previous section. The couplings would fulfil this relation exactly in a fully diffeomorphism invariant setup. However, such a diffeomorphism invariant setup is [3] dependent on the higher couplings g 4 and g 5 . Right: The system (μ, λ 3 , λ 4 , g 3 , g 4 ) dependent on the higher couplings g 5 and g 6 . The higher couplings λ n>nmax are always identified with λ 3 . The blue area marks the region where an attractive UV fixed point was found. At the border of this area the fixed point either vanishes into the complex plane or loses its attractiveness. In both systems the area where the fixed point exists is rather large and contains the identification g n>nmax = g 3 . Conveniently, the area increases for the better truncation, indicating that the system becomes more stable with an improvement of the truncation. The number of attractive directions is uniformly two in the left panel and three in the right panel not at hand. In fact, we can parameterise the breaking of diffeomorphism invariance via these couplings, e.g. by writing g n = g 3 + g n . Here we have designated g 3 as a reference coupling since it is the lowest genuine gravitational coupling. For this reason, it is also the most converged gravitational coupling within this vertex expansion, thus justifying this choice. In general we expect g n to be small and in consequence we vary the identification of the higher couplings only in this part of the theory space of higher couplings. The quantity g 4 is indeed small at the UV fixed point presented in the last section, see (34). More precisely, it takes the value | g 4 /g 3 | ≈ 0.3 at this UV fixed point.
In this analysis we choose to identify and λ 6 = λ 5 = λ 3 for simplicity, and investigate the existence of the UV fixed point as a function of the parameters α 1 and α 2 . In Fig. 3 the area where an attractive UV fixed point exists is displayed in blue. In the left panel, this is done for the previous truncation (μ, λ 3 , g 3 ) [3], and in the right panel for the current truncation (μ, λ 3 , λ 4 , g 3 , g 4 ). At the border of the blue area the UV fixed point either vanishes into the complex plane or loses its attractiveness. Remarkably, both areas are rather large, suggesting that the existence of the UV fixed point is quite stable. Even more conveniently, the area increases with the improved truncation, suggesting that the system is heading towards a converging limit. Note that the number of attractive directions of the UV fixed point is constant throughout the blue areas, namely two in the left panel and three in the right panel.
We further analyse the fixed point values that occur within the blue area in the right panel of Fig. 3 Hence, in particular the fixed point value of λ 3 is already confined to a very small interval, and also a very small number. The latter is important since some of our approximations rely on the fact that the λ n are small, see Sect. 4.5. The fact that λ * 4 is varying more strongly than λ * 3 is not surprising since we expect λ 3 to be better converged, being a lower coupling. The fixed point values of g 3 and g 4 seem to try to compensate the change induced by the identification. Thus, g * 3 and g * 4 become larger towards the identification g 6 = g 5 = 0 and smaller towards g 6 = g 5 = 2g 3 . The shape of the area in the left panel in particular suggests the relation g * 4 < g * 3 , which is fulfilled by the improved truncation almost throughout the whole area where the fixed point exists. This is indeed a nontrivial prediction that has been fulfilled by our approximation scheme. Fig. 4 Plot of the existence of an attractive non-trivial UV fixed point (blue) dependent on the higher couplings g 5 and g 6 . Here, the trilocal equation for the gravitational coupling g 4 was used, which allows us to switch off the R 2 coupling. We found two different fixed points with rather similar fixed point values. Each fixed point has its own area of existence in the theory space of the higher couplings. The blue area marks the unified area of both fixed points. Nevertheless, the area is significantly smaller than the areas displayed in Fig. 3. This reflects the importance of the R 2 coupling A further study of the dependence of the UV fixed point properties on the choice of identification is given in Appendix 1.

Importance of the R 2 tensor structure
In the previous subsection we have fed back the R 2 contributions to the flow via the momentum-dependent gravitational coupling g 4 ( p 2 ). In order to check the quality of our approximation and to investigate the influence of the R 2 tensor structure on the fixed point structure of the system, we switch off the R 2 contribution in this section. We do the latter by projecting onto the p 2 part of the flow via a trilocal momentum projection scheme, cf. Sects. 4.3 and 4.6. This is both an examination of the influence of R 2 on the results presented in the previous subsections, as well as a proof of concept for disentangling the tensor structures of different invariants. Our analysis in this subsection suggests that leaving out the contribution of R 2 leads to significantly less stable results.
In Fig. 4 we display the result for the same analysis as in the previous section, but with the trilocal equation (E7) for g 4 instead. We find two fixed points with rather similar fixed point values. However, we are only interested in identifying the area in the theory space of the higher couplings where at least one UV fixed point exits. Thus, we unify both areas and obtain the blue area displayed in Fig. 4. This area forms a rather narrow band whose total area is significantly smaller than for the momentum dependent gravitational coupling g 4 ( p 2 ), c.f. Fig. 3. The identification g 6 = g 5 = g 3 also does not lie within these regions, but just outside of them. Since we switched off the R 2 contribution, a less stable fixed point structure was to be expected, and consequently these results highlight the importance of the R 2 coupling.

IR behaviour
In this section, we discuss the IR behaviour of the present theory of quantum gravity. We only consider trajectories that lie within the UV critical hypersurface, i.e. trajectories that are UV finite, and which end at the UV fixed point presented in (34) for k → ∞. In this section we use the analytic flow equations given in Appendix 1 for simplicity, and set the anomalous dimensions to zero, i.e. η φ = 0. This approximation gives qualitatively similar results, as discussed in Appendix 1.
In the IR, it is particularly interesting to examine the background couplingsḡ andλ. In the limit k → 0 the regulator vanishes by construction and the diffeomorphism invariance of the background couplings is restored. Hence they become observables of the theory. The flow equations for the background couplings are displayed in Appendix 1.
In general we look for trajectories that correspond to classical general relativity in the IR. This implies that the quantum contributions to the background couplings vanish and in consequence that they scale classically according to their mass dimension. The classical scaling is described bȳ We use the classical scaling in the flow from the UV fixed point to the IR in order to set the scale k in units of the Planck mass M Pl . We need to find a large enough regime wherē g ∼ k −2 . This entails that Newton's coupling is a constant in this regime and sets the scale k via G N = M −2 Pl =ḡk −2 . In Fig. 5, two exemplary trajectories are displayed. In the left panel all couplings scale classically below the Planck scale and reach their UV fixed point values shortly above the Planck scale. All quantum contributions are suppressed simply by the fact that μ → ∞. In the right panel on the other hand some couplings exhibit a non-classical behaviour even below the Planck scale, which is triggered by the graviton mass parameter μ flowing towards the pole of the graviton propagator at μ = −1. This entails that the dimensionful graviton mass parameter M 2 = μk 2 is vanishing in the IR. This IR behaviour is analogous to the one observed in [2], and recently also [50]. Remarkably, not only μ is behaving non-classically but also λ 3 , even though it is not restricted by any pole. However, in this scenario the numerics break propagator at μ = −1. However, in this case the numerics break down at k ≈ 0.02M Pl due to competing orders of the factor (1 + μ) close to the singularity at μ = −1. The trajectories in both panels correspond to theories that behave like classical general relativity in the IR. Note that some couplings are plotted shifted or with a minus sign in order to keep them positive over the whole range Table 1 Properties of the non-trivial UV fixed point for different orders of the vertex expansion scheme, computed for momentum dependent anomalous dimensions η φ i ( p 2 ) and bilocally projected Newton's couplings g n (k 2 ). The critical exponentsθ i andθ i stem from two different approximation of the stability matrix as discussed in Appendix 1. The fixed points are computed with the identifications g 6 = g 5 = g max and λ 6 = λ 5 = λ 3 . We observe that the fixed point values are only varying mildly between the different orders of the vertex expansion. Notably, if we compare the critical exponents of the two approximations of the stability matrix, we observe that the difference becomes smaller with an increasing order of the vertex expansion. This is precisely what one would expect of a systematic approximation scheme that is approaching a converging limit down at k ≈ 0.02M Pl due to competing orders of the factor (1 + μ) close to the singularity at μ = −1.
In the left panel we have tuned the background couplings g andλ so that they are equal to the lowest corresponding fluctuation coupling in the IR, i.e.ḡ = g 3 andλ = λ 2 = −μ/2 for k M Pl . This is equivalent to solving a trivial version of the Nielsen identities (NIs). Since all quantum contributions are suppressed by the graviton mass parameter going to infinity in the IR, μ → ∞, the NI in (2) reduces to In consequence, we should see that all couplings coincide in this limit,ḡ = g n andλ = λ n . This is not the case in the left panel of Fig. 5 since we have only fine-tuned the background couplings, and thus we have two further degrees of freedom that could be used for fine-tuning, stemming from the three dimensional UV critical hypersurface. We postpone this fine-tuning problem to future work. In summary, we find different types of trajectories that correspond to classical general relativity in the IR. The main difference lies in the behaviour of the graviton mass parameter μ, which flows to infinity in one case and to minus one in the other case. Both scenarios are equivalent to general relativity in the end, in particular since only the background couplings become observables in the limit k → 0.

Towards apparent convergence
In this section we discuss and summarise the findings of this work concerning apparent convergence. On the one hand, the order of our vertex expansion is not yet high enough to fully judge whether the system approaches a converging limit. Nevertheless, we have collected several promising first hints that we want to present in the following.
In this work we have introduced two different approximations to the stability matrix, as presented in Appendix 1. We have argued that in a well converged approximation scheme the most relevant critical exponents should not depend on the approximation of the stability matrix. In Table 1 we display the UV fixed point properties for different orders of the vertex expansion. The first system is without the graviton four-point function and exactly the same as in [3]. Then we look at systems where we add either only an equation for g 4 (k 2 ) (c.f. (E5)), or only an equation for λ 4 (c.f. (E4)). Lastly, we display our best truncation including all couplings up to the graviton four-point function, see Sect. 5.1. We observe that the fixed point values of the couplings vary only mildly with an improving truncation, although there is no clear pattern to those variations. The most important piece of information is the difference between the critical exponents from the two different approximations of the stability matrix. While the difference is rather large in the truncation of the graviton three-point function, it is becoming smaller with each improvement of the truncation. At the level of the graviton four-point function, the critical exponents show only a small difference. This is precisely what we expect, and thus we interpret this as a sign that the system is approaching a converging limit.
Another important piece of information comes from the stability of the UV fixed point under different closures of the flow equation. In a well converged expansion scheme, the properties of the UV fixed point should be completely insensitive to the details of the closure of the flow equation. We have performed this analysis in Sect. 5.2. We observed that the area in which the UV fixed point exists in the theory space of higher couplings is indeed increasing with the improvement of the truncation. Furthermore, we saw that the UV fixed point values are confined to small intervals. We again interpret this as a sign that the system is approaching a converging limit.
In summary, we have already seen several signatures of apparent convergence although we are only at the level of the graviton four-point function within the present systematic expansion scheme. This suggests that we are on a promising path and that the present setup will eventually lead to a converging limit.

Summary
We have investigated quantum gravity with a vertex expansion and included propagator and vertex flows up to the graviton four-point function. The setup properly disentangles background and fluctuation fields and, for the first time, allows to compare two genuine Newton's couplings stemming from different vertex flows. Moreover, with the current truncation we have closed the flow of the graviton propagator: all vertices and propagators involved are computed from their own flows.
As a first non-trivial result we have observed that the vertex flows of the graviton three-point and four-point functions, in the sense of (26), are well described by a polynomial in p 2 within the whole momentum range 0 ≤ p 2 ≤ k 2 . The projection used for the flows takes into account the R, R 2 and R 2 μν tensor structures as well as higher order invariants with covariant momentum dependencies. Importantly, it is orthogonal to the R 2 tensor structure for the graviton three-point function, but includes it for the graviton four-point function. We have shown that the highest momentum power contributing to the graviton three-point function is p 2 . Therefore, R 2 μν and higher derivative terms do not contribute to the graviton three-point function. Thus, in particular R 2 μν is excluded as a UV-relevant direction. On the other hand, the flow of the graviton four-point function shows p 4 as its highest momentum power. Together with the three-point function result we infer that R 2 is UV-relevant and contributes to the graviton four-point function. This is a very interesting and highly nontrivial result.
At the moment, we cannot make final statements about higher R n terms directly from our analysis. Nonetheless, predictions can be made with a combination of the results presented here and previous ones obtained within the background field approximation as well as the vertex expansion: Firstly, our work sustains the qualitative reliability of background field or mixed approximations for all but the most relevant couplings. We have seen that the range of allowed Newton's couplings stemming from n-graviton vertices is growing with the level of the approximation. Moreover, in [71,74] it has been shown that already the substitution of the most relevant operator, the mass parameter μ, in a mixed computation with that in the full vertex expansion stabilises the results in a particular matter gravity system. Hence, this gives us some trust in the qualitative results for higher R n terms in the background field approximation. In [31] the f (R)potential has been computed polynomially up to R 34 , and the relevance of these operators follows the perturbative counting closely. Accordingly it is quite probable that the higher R n will turn out to be irrelevant in the full vertex expansion as well.
Based on the above observations we have also constructed projection operators that properly disentangle the contributions of different diffeomorphism-invariant tensor structures. This allowed us to switch off the R 2 coupling in order to analyse its importance for the system. In this case, we are led to an unstable system, which highlights the importance of the R 2 coupling for the asymptotic safety scenario. In the present work we include the R 2 contributions via the momentum dependence of the gravitational coupling g 4 ( p 2 ), leading to a very stable system in the UV.
In the full system with R 2 contributions we found an attractive UV fixed point with three attractive directions and two repulsive directions. The third attractive direction can be explained due to the overlap with R 2 , and is in agreement with previous R 2 studies in the background field approximation [10,[14][15][16]31]. We investigated the stability of this UV fixed point with respect to changes of the identification of the higher couplings and compared it to the stability of the previous truncation without the graviton four-point function. We characterised the stability via the area of existence in the theory space of higher couplings, and remarkably this area increased with the improved truncation. We interpret this as a sign that the systematic approximation scheme is approaching a converging limit.
Furthermore, we investigated the IR behaviour and found trajectories that connect the UV fixed point with classical general relativity. In particular, we found two different types of such trajectories. In the first category all couplings, including background and fluctuation couplings, scale classically according to their mass dimension below the Planck scale. In consequence the Nielsen identities become trivial in this regime and we can solve them in the IR. In the second category, the graviton mass parameter and the coupling λ 3 scale non-classically below the Planck scale, which is triggered by the graviton mass parameter flowing towards the pole of the graviton propagator μ → −1. In summary, the IR behaviour was found to be very similar to [2], and recently also [50].
Lastly, we discussed signs of apparent convergence in the present system by comparing the results to previous truncations. As mentioned before, we observed that the present system is more stable and less sensitive to the closure of the flow equation, which is expected from a converging system. We furthermore used two different approximations of the stability matrix and argued that the critical exponents belonging to the most attractive directions should not differ in a well converged expansion. Indeed we found that the difference of the critical exponents is decreasing with an improvement of the truncation. We interpret this as a sign towards convergence.
In the present approximation we have taken the , R and R 2 tensor structures into account. Furthermore, we have shown that the higher derivative tensor structures and the R 2 μν tensor structure are suppressed. There are also very strong indications for the irrelevance of the higher orders in R n . Altogether this suggests that the natural extension of the present work towards apparent convergence consists primarily of the inclusion of all tensor structures of the vertices and propagators on the given level n = 4 of the vertex expansion. In particular, this concerns the inclusion of the R 2 tensor structure with the orthogonal projection devised in the present work. Moreover, the different graviton modes should be furnished with their separate dispersion or wave function renormalisation. This is well in reach with the current technical status of the programming code and its implementation. Then, selected tensor structures of higher vertices could be used for further tests of apparent convergence. We hope to report on this in the near future.
The critical exponents that correspond to this approximation represent the critical exponents that belong to the computed phase diagram of the theory.
In the second approximation, the identification is performed after taking the derivatives, to wit This approximation is more closely related to the full stability matrix in the sense that it respects the fact that the higher couplings in the full system do not coincide with the lower ones. Note that these two different approximations only differ if we choose a non-trivial identification scheme, i.e. if the higher couplings are functions of the lower ones.
If the expansion scheme is well converged, then the contributions of the higher couplings to the flow of the lower couplings are small, e.g. (∂ λ nmax+2 λ n max )| FP ≈ 0. In this case, the most relevant eigenvalues of the stability matricesB and B coincide approximately. The stabilisation of the most relevant eigenvalues was also observed in an expansion in R n in [31]. In consequence, a huge deviation in the most relevant eigenvalues of both approximations would clearly indicate a lack of convergence. For this reason, we use the comparison of the different approximations as a first check of convergence.

Appendix B: Background couplings
In this section we present the flow equations for the background couplingsḡ andλ. They are in particular interesting in the limit k → 0 where they become observables. In this limit, the regulator term vanishes by construction and diffeomorphism invariance is restored, which implies that these couplings can be interpreted as physical observables only for vanishing k.
For notational convenience we reintroduce the coupling λ 2 = −μ/2. Following [44], we compute the flow of the background couplings with a curvature expansion on an Einstein space. We use a York-decomposition [129,130] and field redefinitions [9,57] to cancel the non-trivial Jacobians. The resulting flow equations are given by where the functions f R 0 and f R 1 read In consequence the fixed point equations for the background couplings are given bȳ , . (B3) Note that the background couplings are non-dynamical, i.e. they do not influence any other coupling. Furthermore, the background couplings only depend on the couplings of the two-point function. Hence only the graviton mass parameter μ (or equivalently λ 2 ) and the anomalous dimensions η h and η c directly affect them.

Appendix C: Dependence on the identification scheme
The flow of each n-point function depends on the couplings of the (n + 1)-point function and the (n + 2)-point function, see Fig. 1. For the highest couplings, we consequently do not have a flow equation at hand. In our setup, these are the couplings of the five-and six-point function, i.e. λ 5 , λ 6 , g 5 , and g 6 . In order to close the flow of our system, we need to make an ansatz for these higher order couplings. A natural choice is one that is close to diffeomorphism invariance, i.e. to identify these couplings with a lower order coupling. In our setup, there are two lower order couplings that correspond to a (partly) diffeomorphism invariant identification scheme, e.g. λ 5 can be identified with λ 3 or λ 4 . In a well converged expansion scheme, the details of the identification should not matter and lead to similar results. In this section, we compare the results for different identification schemes in order to evaluate the stability of our expansion scheme.
The properties of the non-trivial UV fixed point for different identifications schemes are displayed in Table 2. In all identification schemes except for the identification g n>4 → g 4 and λ n>4 → λ 4 we find an attractive UV fixed point. In this case we can see that the fixed point has just vanished in the complex plane. For all other identifications we observe that the fixed point values and the critical exponents vary only mildly. Especially the number of attractive directions is consistently three with the first approximation to the stability matrix, c.f. Appendix 1. With the second approximation the Table 2 Properties of the non-trivial UV fixed point for different identification schemes, i.e. different closures of the flow equations, as discussed in Appendix 1. The flow equations are computed with momentum dependent anomalous dimensions η φ i and bilocally projected Newton's couplings g n (k 2 ). The critical exponentsθ i andθ i stem from two different approximation of the stability matrix as discussed in Appendix 1. An attractive UV fixed point is found in most identification schemes with mildly varying fixed point values. In the first approximation of the stability matrix we always find three attractive directions, while in the second approximation of the stability matrix we find one or three attractive directions, since the real part of one complex pair of eigenvalues is quite close to zero. These results suggest that the present system is rather stable under change of the closure of the flow equations. In the case of the single identification without a physical UV fixed point we found that it had in fact just vanished in the complex plane number of attractive directions varies from one to three since the real part of one complex conjugated pair is close to zero.
In conclusion, this analysis suggests that our system is rather stable with respect to different identification schemes. Only one particular identification scheme has led to the disappearance of the attractive UV fixed point. This constitutes further support for our results in Sect. 5.2, where we found that our full system is very stable with respect to the identification of g 5 and g 6 .

Appendix D: Possible issues of a local momentum projection
In this section we want to point out some possible issues of a local momentum projection. A local momentum projection is for example a derivative expansion about a certain momentum, usually p = 0.
The full solution of a flow equation includes a full resolution of the momentum dependence of all vertex flows. For higher n-point functions this task is computationally extremely challenging due to the high number of momentum variables. We have already argued in Sect. 4.2 that this task can be tremendously simplified with a symmetric momentum configuration. We have further shown in Sect. 4.3 that the quantity Flow is polynomial in p 2 , at least for n = 3, 4. Thus it is possible to consistently project on each coefficient of this polynomial in the whole momentum range 0 ≤ p 2 ≤ k 2 by employing a non-local momentum projection. In contrast, a local momentum projection scheme does not capture the correct momentum dependence over the whole momentum range 0 ≤ p 2 ≤ k 2 in general since it is sensitive to local momentum fluctuations. Furthermore, it is very challenging to project on the p 4 coefficient or even higher momentum order coefficients due to IR singularities. All these statements are explicitly exemplified in Fig. 6.
Projection at p 2 = 0: 0.66p 2 Fig. 6 Fit of a local momentum projection at p 2 = 0 to the momentum dependence of the flow of the graviton three-point function (left) and the graviton four-point function (right) divided by (− n 2 η h ( p 2 ) − n + 2) as defined in (26). The flows are evaluated at (μ, λ 3 , λ 4 , g 3 , g 4 ) = (−0.4, 0.1, −0.1, 0.7, 0.5) and λ 6 = λ 5 = λ 3 as well as g 6 = g 5 = g 4 , i.e. the same values as in Fig. 2 where a non-local momentum projection was used. In comparison to the non-local momentum projection, the local momentum projection does not capture the correct momentum dependence in the whole momentum range 0 ≤ p 2 ≤ k 2 since it is sensitive to local momentum fluctuations. Furthermore, it is technically very challenging to project on the p 4 -term due to IR singularities. Note again that the constant parts of the flows are irrelevant for the beta functions since they are extracted from a different tensor projection On the other hand, the local momentum projection at p = 0 has the advantage that it allows for analytic flow equations, as discussed in Appendix 1. Analytic flow equations are more easily evaluated in the whole theory space, but, as the discussion above suggests, one should be be mindful of the fact that they easily introduce a large error.
We use the analytic flow equations in Sect. 6 precisely for the reason that they can easily be evaluated in the whole theory space. Thus we show now that the fixed point properties in this analytic system are qualitatively similar to the full system, despite the error that is introduced by the analytic equations.
The properties of the UV fixed point for different approximations are displayed in Table 3. Truncation 1 corresponds to our full system, i.e. with momentum dependent anomalous dimensions and bilocally evaluated gravitational couplings. Truncation 4 corresponds to the system used in Sect. 6, i.e. without anomalous dimensions and with gravitational couplings from a derivative expansion. Truncation 2 and 3 are in between those truncations, i.e. with anomalous dimension but gravitational couplings from a derivative expansion and without anomalous dimensions but with bilocally evaluated gravitational couplings, respectively.
We observe that the UV fixed point exists in all truncations, and that the properties of this fixed point vary only mildly. The fixed point values are all located within a small region, with the exception of our simplest truncation. There, the couplings λ 3 and λ 4 have a different sign compared to the other truncations. Considering the critical exponents we always find three attractive directions with the first approximation of the stability matrix and in three out of four cases three attractive directions with the second approximation of the stability matrix as well. These results suggest that it is an acceptable approximation to use the analytic flow equations if one is only interested in the qualitative behaviour of the system.

Appendix E: Derivation of flow equations
We obtain the flow equations for the individual coupling constants by projecting onto the flow of the graviton n-point functions, as explained in Sect. 4.5.
The equations for the graviton mass parameter μ and for the graviton anomalous dimension η h are extracted from the transverse-traceless part of the flow of the graviton two-point function. For p 2 = 0, we obtain for the flow of the graviton mass parameter.
We obtain an equation for the graviton anomalous dimension η h ( p 2 ) by evaluating the flow of the graviton two-point function bilocally at p 2 and −μk 2 , to wit The ghost anomalous dimension η c ( p 2 ) can be directly obtained from the transverse flow of the ghost two-point function, to wit In case of the higher order couplings, we employ the projection operators described in Sect. 4.2. For the couplings λ n , this leads to where the projection dependent constant C n is implicitly defined via n • n TT • T (n) (0; 1) = C n . Here, • denotes the pairwise contraction of indices.
As discussed in Sect. 4.5, the gravitational couplings g n ( p 2 ) are momentum dependent. In order to simplify the computation we make an approximation of the full momentum dependence. This approximation exploits the fact that the flows are peaked at p 2 = k 2 and consequently we set the feed back on the right-hand side of the flow equation to g n ( p 2 ) ≈ g n (k 2 ). This closes the flow equation for g n (k 2 ) and thus we only solve this equation. The easiest way to obtain the flow equation for g n (k 2 ) is a bilocal projection at p = 0 and p = k. The resulting equation for g 3 (k 2 ) is precisely the same as in [3]. For g 4 (k 2 ) we obtain The derivation of this equation is based on the assumption that λ 4 is small. In Sect. 4.6 we have laid out a strategy to disentangle contributions from different tensor structures, in particular those of R and R 2 . The flow equations for the g n are obtained by a projection onto the p 2 part of Flow (n) G divided by − n 2 η h p 2 − n + 2 , see Sect. 4.3 and Sect. 4.6. The graviton three-point function is at most quadratic in the external momentum, and consequently it is again enough to use Table 3 Properties of the non-trivial UV fixed points for different approximations within the full system. In the truncations 3 and 4, we set η φ i = 0, while in the truncations 1 and 2 we use momentum dependent anomalous dimensions. In truncations 2 and 4, the couplings g 3,4 are computed via a derivative expansion at p = 0, while in the truncations 1 and 3 the couplings g 3,4 (k 2 ) are evaluated with a bilocal projection between p = 0 and p = k. The quality of the truncation decreases from 1 to 4. The fixed point values are obtained with the identification scheme λ n>4 = λ 3 and g n>4 = g 4 . The critical exponentsθ i andθ i stem from two different approximation of the stability matrix as discussed in Appendix 1. The fixed point properties from different approximations are qualitatively very similar. In particular, all fixed points exhibit three relevant directions when the first approximation of the stability matrix is used. Using the second approximation of the stability matrix also results in three relevant directions in three out of four cases a bilocal projection at p 2 = 0 and p 2 = k 2 . Consequently, the flow equation for g 3 is quantitatively equivalent to the previous one if λ 3 is small. The graviton four-point function, on the other hand, has p 4 as its highest momentum power, and thus we use a trilocal momentum projection at p 2 = 0, p 2 = k 2 /2, and p 2 = k 2 . The flow equations of g 3 and g 4 are then given by (1 + η 4 ) ∂ t g 4 = 2g 4 − g 4 C 4 (∂ t λ 4 + 2λ 4 ) − 1 η h (k 2 ) + 1 where The constants C are implicitly defined via G n • n TT • T (n) ( p 2 ; n ) = C G n n n + C G n p 2 p 2 , and we use the abbreviation C n := C G n n /C G n p 2 . Again, • denotes the pairwise con-traction of indices. Note that the constants η n are chosen in such a way that η n = 0 for vanishing anomalous dimensions. Analogously, we can obtain a flow equation for the R 2 coupling of the graviton four-point function ω 4 by using a trilocal momentum projection, as explained in Sect. 4.6. We evaluate the flows at the same momenta as for the trilocal flow equation of g 4 . The equation for ω 4 then reads where The constants C are again defined via the contraction G n • n TT •T (n) ( p 2 ; n ) = C G n n n +C G n p 2 p 2 +C G n ω n 4 p 4 . Again, η ω is defined such that η ω = 0 for vanishing anomalous dimensions.
In the previous paragraphs we introduced abbreviations for constants that arise from the projection scheme. The explicit values of these constants are: