The derivative expansion in asymptotically safe quantum gravity: general setup and quartic order

We present a general framework to systematically study the derivative expansion of asymptotically safe quantum gravity. It is based on an exact decoupling and cancellation of different modes in the Landau limit, and implements a correct mode count as well as a regularisation based on geometrical considerations. It is applicable independent of the truncation order. To illustrate the power of the framework, we discuss the quartic order of the derivative expansion and its fixed point structure as well as physical implications.


Introduction
One of the major open problems in fundamental physics is the formulation of a consistent quantum theory of gravity. Despite of several decades of research effort that went into it, no completely consistent and experimentally verified solution exists so far. A big contributor to this status of the field is that quantum gravity effects are expected to be extremely small, and even one-loop effects appear to be unmeasurably tiny at the present resolution of experiments. Thus a lot of guidance must come from theoretical considerations and consistency with Standard Model physics.
A conservative approach to construct a theory of quantum gravity was theorised by Weinberg [1], and goes under the name of Asymptotic Safety. It tries to achieve the quantisation of gravity via postulating quantum scale invariance at high energies induced by a second order phase transition. On the technical level, this translates to an interacting fixed point of the renormalisation group flow. Results obtained in 2 + dimensions indeed suggest the existence of such a fixed point, at least near two dimensions [2][3][4][5][6][7], but the extrapolation to the physical case of four dimensions remained difficult for some time.
One systematic way to study the stability of these results is the derivative expansion. In this, interactions which include up to a set amount of derivatives acting on the fundamental field are taken into account. In the context of gravity, this corresponds to powers of curvature tensors and their covariant derivatives. Surprisingly, to this point a complete non-perturbative discussion of the fourth order approximation has not been carried out in the context of Asymptotic Safety. This has both conceptual and technical reasons. Conceptually, the functional renormalisation group relies on the choice of a regulator. The regularisation of operators in a curved spacetime is much more involved than in a flat spacetime. Only partial results have been obtained regarding this problem, and most rely on the particular structure of the theory at second order in derivatives. On the technical side, the computational complexity increases tremendously with the approximation order.
In this paper, we provide a solution to the problem of regularisation in asymptotically safe quantum gravity motivated by geometrical arguments. It is based on the decomposition into gauge variant and invariant components of the field, and can be applied to any order of the derivative expansion, including resummations in terms of form factors. We also provide new technical insights that allow for a generic implementation of the computation of the renormalisation group running. To illustrate the framework, we compute the complete non-perturbative renormalisation group running in quantum gravity to fourth order in the derivative expansion.
The main results that we obtain in the concrete computation are as follows: • Fourth order gravity admits a non-perturbative fixed point.
• Spacetimes with a negative (positive) Euler characteristic dominate (are suppressed in) the Euclidean path integral of Asymptotic Safety, whereas spacetimes with vanishing Euler characteristic contribute with unit strength.
• While the inclusion of the square of the Weyl tensor provides a well-controlled extension of the Reuter fixed point, the squared Ricci scalar introduces some kind of instability into the results. This has been observed before, and the inclusion of higher order terms seems to stabilise the system [19,21,23,29,32,34,36,39,43,47,[51][52][53]110].
This paper is structured as follows. In section 2 we briefly discuss the functional renormalisation group (FRG) which lies at the heart of our investigations of Asymptotic Safety. In particular, we define a set of criteria that well-behaved regulators and flows should satisfy. Section 3 discussed the decomposition of fields into gauge variant and invariant components. We first provide a simpler discussion in the context of an Abelian gauge field, and then investigate how much we can transfer the structure to the gravitational case. This discussion leads to our proposal of a well-motivated regularisation scheme in gravity. Sections 4, 5 and 6 as well as appendix A collect some technical machinery that allow us to decrease the technical complexity to a manageable level. The setup is illustrated in section 7, where we carry out the computation to fourth order in derivatives. We then conclude and provide an outlook in section 8.

Functional renormalisation group
The tool that we will use to investigate the non-perturbative renormalisation group flow is the functional renormalisation group. Its central object is the effective average action Γ k , which interpolates between a microscopic action S in the limit k → ∞, and the standard quantum effective action Γ at k = 0. Decreasing the fiducial scale k then corresponds to integrating out modes in the Wilsonian sense. The dependence of Γ k on k is governed by the following functional integro-differential equation [8][9][10]: In this, ∂ t = k∂ k is the logarithmic scale derivative, Γ k denotes the second functional derivative of Γ k , R k is a regulator term which acts as a momentum-dependent mass, and the super-trace STr indicates a sum over discrete (e.g. spacetime or gauge bundle), and an integral over continuous (e.g. momentum) indices. For reviews of the FRG see [97,[111][112][113][114][115].
From the flow equation (1) we can extract the beta functions of a theory, i.e., (integro-)differential equations that govern the scale dependence of couplings. Typically, one discusses them for the dimensionless versions of couplings, where we multiply the coupling with a power of the scale k to make it dimensionless. Combined zeros of all beta functions are called fixed points. We will indicate fixed point values of couplings by an asterisk. If all couplings vanish at a given fixed point, it is called Gaussian, otherwise we call it interacting or non-Gaussian. A fixed point is then characterised by its critical exponents, which describe the linearised flow around it. They are defined as minus the eigenvalues of the matrix of partial derivatives of the beta functions with respect to the couplings, evaluated at the fixed point. Positive (negative) critical exponents indicate relevant (irrelevant) directions, so that the flow is attracted towards (repelled away from) the fixed point when increasing the scale k. To end up with a predictive theory, a fixed point can only have finitely many positive critical exponents, since they correspond to the number of independent measurements that one has to perform to fix the theory uniquely. An asymptotically safe fixed point is then defined as an interacting fixed point with a finite number of positive critical exponents.
In quantum gravity, we have to add a gauge fixing term to the action to make the propagator that appears in the flow equation well-defined. On the practical level, this is implemented with the help of the background field method. We split the metric g into an arbitrary background metricḡ and (not necessarily small) fluctuations h around it, Other (non-linear) ways to perform this separation have been investigated, see e.g. [28,30,32,33,39,42,51,57,[116][117][118][119][120][121][122] for examples in the context of Asymptotic Safety. Such a split is also necessary to define the regulator term. The advantage of this method is that invariance with respect to background diffeomorphisms can be maintained in every step. However, as a downside, the regulator and the gauge fixing term break certain Ward identities, so that in principle one has to deal with modified Ward identities that have to be fulfilled together with the flow equation. In this work, we will focus on a background field approximation, which neglects these issues and corresponds to setting the fluctuation field h to zero after taking the second variation. For an in-depth discussion of these issues, see e.g. [15, 18, 23-25, 27, 38, 40, 41, 46, 49, 50, 56, 59, 67, 71, 98, 122-140].
In practice, we generally have to make approximations to solve the flow equation (1). Only in special cases, an exact solution is possible, see [141]. For recent progress in constructing exact solutions to the flow equation, see also [142]. In the following, we will discuss a (covariant) derivative expansion of the effective average action. The order of the expansion is then the maximal number of derivatives that act on the metric.
A generic problem in this setup is the systematic choice of a regulator. Partial results have been obtained in the literature, notably [12], but they typically rely on technical assumptions that potentially do not carry over to higher orders, or only at considerable technical cost. Two of the goals of this paper are to establish general criteria that a good flow should have, and a generic way to choose a regulator which gives rise to a good flow. Before we enter this discussion, we briefly note that in this paper we consider a Euclidean flow. Results with Lorentzian signature can be found in [13,37,100,117,[143][144][145][146][147][148]].

Criteria for a successful flow
Having introduced the machinery, we will now discuss some criteria that we expect a flow to have. The first condition is what we call a correct mode count. The idea is that the flow of the cosmological constant should have a very generic form, and effectively counts the number of physical degrees of freedom. For example, for a free scalar field, the contribution to the flow of the cosmological constant is In this, the factor of a half comes directly from the flow equation (1), and the rest of the prefactor of the integral, as well as its measure y d 2 −1 , come from the heat kernel. The integrand is then just the product of the regularised propagator and the scale derivative of the regulator. We thus demand that the contribution of any physical degree of freedom comes in this form: an integral over the full regularised propagator with the prefactors as above.
A second criterion that we want to implement is that the flow should admit a finite Landau gauge limit. It corresponds to a strict implementation of the gauge condition, and is a fixed point of the flow [41,149]. This puts indirect constraints on the regulator as well, which are however hard to spell out concretely. A well-motivated setup with an a priori guaranteed finite Landau limit will be presented below. In any case, this can be easily checked a posteriori, once the flow is computed.
A third criterion concerns the choice of operator that is regularised, and the tensor structure of the regulator. As a convention, we will always assume that the operator is normalised such that the (principal) Laplace part of the operator comes with a unit prefactor. We then require that the regularised operator, and the tensor structure of the regulator, have a physical or mathematical motivation. We will discuss this more concretely in the next section.

Field decompositions in curved spacetime
To solve the problem of finding a suitable regularisation, we will now discuss the decomposition of fields into components. A guiding principle will be that we rescale the modes such that no non-trivial Jacobians are introduced into the path integral. As a simple example, we will start with an Abelian gauge field, for which we can easily derive all necessary ingredients. We then make a short digression to discuss more general vector fields, and discuss their regularisation. This includes the gravitational Faddeev-Popov ghost, which has a slightly different structure than an Abelian gauge field. Finally, we will discuss the mode decomposition of the graviton and the regularisation strategy that this suggests.
Let us mention that in practice, it is easier to avoid working with decomposed fields when it comes to computing functional traces like in (1). We will provide an argument for this in subsection 3.4. The aim of this section is to motivate our choice of regularisation from a geometric perspective. In concrete computations, we then implement the regularisation on the level of the full fields with suitable projectors so that we can use standard heat kernel techniques. We will nevertheless provide all details in the hope that some readers with different applications might find them useful.

Transverse decomposition of an Abelian gauge field
We will now discuss the decomposition of an Abelian gauge field into modes, and their respective regularisation. The starting point will be the gauge fixing condition, from which we derive the decomposition. We then rescale some of the fields to eliminate the need of Jacobians. After a discussion of the projectors onto the different modes, we discuss the natural operator that arises from the decomposition, including its heat kernel properties. This will suggest a particular way to regularise the theory. Finally, we will make some comments on the origin of the natural operator, and briefly discuss how to regularise higher derivative Abelian theories.
Gauge fixing condition The standard choice for a linear, covariant gauge fixing is where D is the covariant derivative. The sought-after decomposition should be such that only the gauge mode appears in this expression, i.e., the physical mode of the decomposition is annihilated by F. In this case, it means that the physical mode is (covariantly) transverse.
Transverse decomposition By this argument, the decomposition into gauge invariant and gauge variant modes is the well-known decomposition into transverse and longitudinal components, Here, A T µ is the gauge invariant transverse mode andφ a scalar field underlying the longitudinal gauge variant component. In consequence, acting with the gauge fixing operator on the field gives where we introduced the Laplacian At this point, let us mention that there is still a potential gauge redundancy in the transverse mode. It can be shifted by the derivative of a solution to the Laplace equation, We will not address this problem in this work.
Jacobian When calculating an integral, every variable transformation gives rise to a Jacobian. This is also the case when we want to calculate a path integral and perform the decomposition (5). The arising Jacobian can be calculated by a standard trick [150]. We can consider the exponential integral Since the overall normalisation of the path integral is inessential, we neglected overall factors. We also dropped boundary terms upon integration by parts. The Jacobian has to be chosen to cancel the determinant, so that Avoiding the Jacobian We would like to avoid the introduction of such determinants, and formulate a decomposition which has field components of the same mass dimension. This is not the case for (5) -the scalarφ has a relative mass dimension of one unit less in comparison to A and A T . In this case, the solution is straightforward: we define a new scalar field This is well-defined as long as we exclude potential negative or zero modes of the Laplacian. With this definition, the new decomposition of the vector field does not give rise to a Jacobian. A path integral over A is thus the same as a path integral over A T and φ, up to the aforementioned subtleties of individual modes.
Projectors Let us now write down the projectors onto the transverse and longitudinal components. It is easy to see that project onto the longitudinal and transverse components, Anticipating the discussion for the graviton case, note that the projector onto the gauge mode is entirely built up from the gauge fixing operator, as we can write This might seem trivial in the Abelian case, but the structure is true more generally, and a consequence of demanding that different modes are orthogonal to each other.
Natural operator In connection with the projection operators, we will introduce the concept of the "natural" operator associated with the field, ∆ A . We define it as an operator of Laplace type which commutes with the projectors, and has a compatible index structure such that it maps a given field to a field with the same index structure. The principal part of the operator is then normalised to one. For a vector field the operator ∆ A can be constructed easily. Observe that so that is the sought-after operator. When calculating the renormalisation group running of couplings, using this operator simplifies calculations. The above operator is the unique local, Laplace-type operator whose set of eigenfunctions splits into transverse and longitudinal eigenfunctions.
Heat kernel coefficients of ∆ A Let us illustrate the special role of the operator ∆ A by considering its heat kernel coefficients. While it has been found that the heat kernel coefficients of a pure Laplace operator in the space of transverse functions is singular in even dimensions [151], we will illustrate now that this is not the case for the natural operator (17). To that extent, we consider the (traced) heat kernel coefficients in the longitudinal sector, The trace can be calculated with standard off-diagonal heat kernel techniques [151,152]. One finds that the heat kernel coefficients agree precisely with the heat kernel coefficients corresponding to a pure Laplacian acting on a scalar, This should not come as a surprise -the gauge fixing operator equips the longitudinal scalar φ with a plain Laplacian, One different way to interpret (19) is that the vector version of the Laplacian acting on a scalar is the operator ∆ A . From this, we can calculate the transverse heat kernel, The total contribution of a free Abelian gauge field in curved spacetime is thus This seems like a trivial statement -the contribution of a vector is the sum of the contributions of the individual modes. In a quantum field theory setting, where regularisation is necessary, this becomes a guiding principle. Only those regularisations that preserve this additive structure are proper. In particular, if one were to regularise only the flat part of the operator, this sum rule is violated, and heat kernel coefficients diverge in even dimensions [151]. This finding has particular relevance for the flow equation. For a classification of different regulator types see [153]. We thus can finally formulate our regularisation strategy for an Abelian vector field. The transverse part is regularised using the operator ∆ A , so that 1 In the longitudinal sector, we can either also adapt this operator if in practice we work without the explicit decomposition, or equivalently, we use a pure Laplacian if we use the decomposition, In practice, it is useful to expand the regulator in powers of curvature. This can be done easily, see e.g. [154]. This is possibile even within a form factor setup, see e.g. appendix C of [139].
General actions and identifying propagators and interactions Coming back to a more general scope, the precise shape of the operator ∆ A is no coincidence. The action of a free Abelian gauge field is proportional to the square of the field strength The corresponding two-point function of a free Abelian gauge field reads which is precisely the natural operator found above, together with a projector onto the transverse mode. This observation has a profound consequence for the definition of the nonperturbative Abelian gauge field propagator in curved spacetime. The most general term quadratic in the Abelian gauge field strength can be written as where E is some operator which is independent of A. If we want to preserve the above structure, namely that this action gives rise to a two-point function of the transverse Abelian gauge field of the form 2 with an arbitrary function e which defines Abelian gauge field propagation in curved spacetime, we have to chose The operator in brackets can be derived by demanding together with the assumption that it is a local Laplace-type operator. The above gives a unique prescription of how to split the action of an Abelian gauge field in an arbitrarily curved spacetime into propagator and interaction terms: first collect the pieces that survive in the flat spacetime limit, then complete the operator to E to lift it into curved space. This represents a minimally coupled Abelian gauge field with a non-trivially momentum-dependent propagator. Any term with two Abelian gauge field strengths and some power of curvature that is not of this form is then a genuine interaction term. Note that the regularisation prescription that we outlined above is still applicable.
2 Formally, we define a function of an operator by either its Taylor series, or an inverse integral transform of exponential type, e.g., an inverse Laplace transform. This covers most interesting functions. In particular, it includes the logarithm via To prove some formulas, we will work with inverse Laplace transforms. All manipulations that we perform here and later in the paper will however also go through without significant changes for functions like the logarithm.

General vector field
Before we continue with the case of the graviton, let us briefly discuss the regularisation of more general vector fields. This will also cover the gravitational Faddeev-Popov ghost. A general second order two-point function for a vector looks like Here b is a number andẼ is a multiplicative operator (often referred to as endomorphism).
For the gravitational Faddeev-Popov ghost, and with the gauge fixing (41) defined below, we have For later reference, we will call this ghost operator ∆ c , Let us rewrite the operator ∆ in terms of the operator ∆ A : 3 where we introduced the shifted endomorphism Due to the central role that is played by ∆ A , we will still use it as the operator in our regulator choice. This means that we treat E as the curvature interaction term. Concretely, the regularised version of (35) reads For simplicity, we chose the same regulator shape function in both sectors. This regularisation satisfies the mode count requirement: the contribution of a vector trace in the flat limit counts as d modes, independent of the value of b. As a matter of fact, for E = 0 the dependence on b drops out. This is reasonable since when we decompose the vector into transverse and longitudinal parts, we still have to canonically normalise the field φ. The rescaling then eliminates all occurrences of b. The above regulator choice implements this idea in the presence of a finite endomorphism. In turn, the flow equation automatically takes care of the abovementioned rescaling if we regularise like in (37). With standard off-diagonal heat kernel methods, one can then derive the contribution of the trace of a vector regularised in such a way within the FRG in a derivative expansion. In general dimension d, it reads In the last line we neglected terms with more than four derivatives, and we introduced the integrals Higher orders can be calculated systematically. As expected from the general form of the trace, terms with m powers of E come with (m + 1) powers of the propagator in the integrals. From the prefactor of the integral of the volume term, we can explicitly see that the mode count of d modes for a vector is implemented correctly. Since later we are interested in d = 4, we have to be careful in the evaluation of (d − 4)I 3 1 . One can show that lim

Decomposition of the graviton
We will now turn our attention to the decomposition of the graviton. In doing so, we will try to follow the same steps as for the Abelian gauge field. As it turns out, much of the construction can be done in a similar way, but there are some key differences. From the general theory of irreducible representations, we anticipate a rank two transverse traceless tensor, a transverse vector and two mixing scalars. The scalars can be diagonalised with respect to the gauge fixing, so that one linear combination is gauge invariant, while the other is gauge variant, and will be recombined with the pure gauge transverse vector. In the construction of the decomposition, we took inspiration from [150], but with view on our goal of defining a useful regularisation and avoiding Jacobians, our implementation differs in some details.
In the continuum approach to quantum gravity, the use of the background field method is hard to avoid. Thus, as indicated earlier, in the discussion below we will make use of it and construct the decomposition with respect to background quantities, indicated by an overbar. For an alternative approach towards defining a decomposition with respect to the full metric, see [23].
Gauge fixing condition Again we start by specifying a gauge condition. For gravity one typically considers the one-parameter family of linear covariant gauges where β is a gauge parameter that determines the way of how the two scalar modes mix to give the gauge invariant and the gauge variant scalar mode.
Transverse traceless decomposition We make the following ansatz for the decomposition of the metric fluctuation: Here h TT is the transverse traceless mode, which is the gauge invariant spin two mode. The pure gauge vector ζ is introduced by means of the gauge fixing operator, similar to the Abelian case. We do not decompose it further into transverse and longitudinal components. The scalar θ is then the gauge invariant scalar. By construction, h TT is annihilated by the gauge operator, We also require that the gauge condition also annihilates θ, Let us construct the operator Q. Observe that for the choice β = 0, the gauge fixing operator is traceless, so that the gauge invariant scalar mode is the trace. This motivates the ansatz where Q is symmetric. In the following we will assume that there is no term proportional to the metric in Q µν . If there were such a term, we could pull it out and rescale the field θ to enforce a unit coefficient as in the above equation. Acting with the gauge fixing operator on this gives Let us now assume that β = 0 so that we can fix the operator Q. We can rewrite the equation as This equation tells us that acting with a derivative from the left on the operator Q and contracting should essentially give back the derivative, i.e. it is some kind of longitudinal projector, but with an unusual index structure. For that reason, let us make the ansatz and derive the form of X so that the above equation is fulfilled. Note again that we could add a term proportional to the background metric to Q, but this would only yield a total rescaling of Q, as discussed above. With (48), we get We conclude that X must be the inverse of the operator in the brackets, In fact, this is precisely the kinetic operator of the Fadeev-Popov ghosts (34) associated to the gauge fixing operator (41), This means that, up to the condition that which is needed for positivity and invertibility, see (35), the operator and its inverse should exist inside the first Gribov region.
Jacobian Having derived the decomposition into gauge invariant and gauge variant modes, let us calculate the Jacobian that arises from this variable transformation. We consider the same integral as for the case of the Abelian gauge field. Before we do that, we first define the operator upon neglecting boundary terms. By this it is clear that Q † annihilates the gauge condition, Note that F † = −F since it is a linear differential operator, so we will omit the dagger symbol for it. Combining all properties, we see that in the calculation of the Gaussian integral all off-diagonal terms, that is those that mix the different modes, vanish. We also see immediately that the transverse traceless sector does not give rise to a Jacobian. The gauge vector integral reads so that the corresponding Jacobian is The choice of normalisation will be made clear below. In the gauge invariant scalar sector, so that the Jacobian reads Avoiding the Jacobians Once again, we would like to avoid the introduction of these Jacobians. For that matter, we rescale the fields by again assuming that all involved operators exist. The decomposition of the metric fluctuation into the set (h TT , ξ, σ), then gives rise to no Jacobians, and the decomposed fields have all the same mass dimension.
Projectors I We can now construct the projectors onto each of the individual components.
In doing so, we make use of the properties of the gauge fixing operator F and of Q. Let us start with the gauge invariant scalar. Acting with Q † on (63) gives From this it is immediately clear that the projector onto this mode reads In a similar fashion, acting with the gauge fixing operator onto h gives so that the gauge projector is Finally, we define the projector onto the TT mode by subtracting the two other projectors from the symmetric identity, The symmetric identity is defined as and maps symmetric rank two tensors to themselves. Inserting explicit expressions into these projectors seems to indicate that Π 2 depends on β. We will now rewrite everything to show that this is actually not the case.
Rewriting the operators In the above expressions, we have two different inverse operators, one constructed from the square of the gauge fixing operator, whose explicit form clarifies the choice of prefactor, the other is X which appears in the operator Q, The two operators agree for the gauge parameter choices β = 0, −1. It will be convenient to formally expand the operators in a Taylor series in β around zero, and resum the full series once the inverse is calculated. The central operator then is Once again we assume that the inverse of ∆ 1 exists. Using a geometric series, we can write X as From the first to the second line, we rewrote the terms in the sum into a form of another geometric series, which is performed in the next step. We also defined the scalar operator The inverse of the squared gauge fixing operator (70) can evidently be obtained from that result by the replacement β → −β 2 , so that Before we go back to the explicit form of the projectors, we can re-express the operator Q as This also entails the compact form of the expression Projectors II We will now present the explicit form of all projectors. The scalar projector reads For the gauge projector we find, after a short calculation, (79) Combining the two into the TT projector, we get Here we used the traceless projector to bring the expression into a compact form, where in the second equation we also introduced the trace projector Π Tr . As promised above, the TT-projector is indeed independent of the gauge parameter β.
Natural operator An obvious question is whether we can define a natural operator for the graviton. This would be a local Laplace-type operator which commutes with the projectors. As it turns out, such an operator does not exist. One can show this in the following way. Assume that there is an operator ∆ 2 which commutes with the spin two projector (80). In that case, we would have that since the projector is transverse. We can easily write down the most general form that this local operator can take, HereC is the background Weyl tensor, see (108) below, and the c i are numerical coefficients. All other potential tensor structures vanish when they act on h TT . Inserting this ansatz into (82), one finds that there is no choice of c i to make this equation true. One can find a nonlocal solution to (82), but due to the inherent difficulties in handling such operators, we will avoid that path in this work, and rather look for alternatives for the operator that we want to regularise.

Regularisation and decoupling in gravity
Having discussed the decomposition of the graviton into gauge variant and invariant modes, but not having found a natural operator, we now have to construct a regularisation scheme by other means. Let us first discuss the gauge variant vector mode. By means of the decoupling theorem [12], the vector mode decouples from the gauge invariant modes completely in the Landau gauge limit, which implements the gauge fixing condition strictly. This means that the functional trace (1) splits into the gauge invariant sector which involves the correlation functions derived from the given action, and a simple vector trace of the form (38) with the operator (70). At the same time, the trace over the Faddeev-Popov ghosts is the same trace but with the operator (34). As noted earlier, the two operators agree if either β = −1 or β = 0. These gauge choice thus implement an exact partial cancellation of these traces. The cancellation is only partial due to the extra factor of two for the ghosts. The regularisation of these two vector fields then follows the discussion in subsection 3.2. Now we will discuss the gauge invariant part. From the explicit form of the spin zero projector (78), we see that the gauge choice β = −1 leaves the gauge invariant scalar nonlocal. We would like to avoid such non-localities, so we will settle for β = 0 as our preferred gauge choice. In this case, the gauge invariant scalar mode is just the trace of the fluctuation field h.
To finally fix the regularisation, we take inspiration from the Abelian case and consider the two-point correlation function of the simplest gravitational action -the Einstein-Hilbert action (without cosmological constant): We then decompose the field via (63) where we can neglect the gauge variant vector since it decouples. In d = 4, the two-point function turns out to be diagonal, and up to overall prefactors containing G N , it has the two parts Here h = h µ µ is the trace of h µν , which is the field σ for β = 0. We thus propose a regularisation involving ∆ 2 and ∆, with traceless and trace projectors, respectively: The numerical prefactors are fixed uniquely by requiring that the regularised version of (85) reads Let us explain why we do not have to use the transverse traceless projectors for the regulator. First, the regulator clearly does not add a regularisation to the trace mode by construction, since we chose β = 0. Second, the regulator has some overlap with the gauge variant vector mode. However, since we implement the strict Landau gauge limit, the contribution of this regulator to the vector mode drops out of any final trace. It is thus unnecessary to employ the spin two projector, and we can rather use the much easier traceless projector. Note that this regulator also fulfils our mode count requirement.
Let us note that in dimensions other than four, the two-point function obtained from the Einstein-Hilbert action is not diagonal -there are off-diagonal terms of the form We will not discuss the regularisation with such off-diagonal terms here, since we are exclusively interested in the physical case d = 4 in this work.
(Not) Using the decomposition Having discussed the regularisation, we will briefly clarify why using a decomposition in practice does not necessarily simplify computations except in special cases. The central reason is that the inverse of a projected operator (in the projected subspace) is in general not the projected inverse of the unprojected operator. As a concrete example (neglecting the regularisation), assume that we would want to invert a plain Laplacian in the transverse sector of a vector field. That is, we are looking for the inverse where the inversion is to be understood in the transverse subspace. To compute this, we can complete the Laplacian to the natural vector operator ∆ A and invert the operator via a geometric series: Here Ric indicates the Ricci tensor. By contrast, the projected inverse reads These two expressions do not agree for an arbitrary manifold. In general, the two expressions only agree if the operator commutes with the projector. This once again highlights the central role of a natural operator. Coming back to gravity, in the previous subsection we found that there is no local natural operator for the transverse traceless part of the graviton. Rather, the Einstein-Hilbert action suggests to consider the spin two operator ∆ 2 , which does not commute with the spin two projector Π 2 . As a consequence, The inversion on the left hand side is to be understood on the space of transverse traceless tensors. This means that even if we use the decomposition, the computation of the transverse traceless propagator is complicated by the presence of the projectors.

Inversion of the graviton two-point function
To compute the non-perturbative renormalisation group flow with the FRG, we have to compute the regularised propagator, which is the inverse of the regularised two-point function. In this section we briefly illustrate how to perform this inversion for the graviton in a derivative expansion, without specifying a particular background metric. The idea is based on the observation that if we were to compute the propagator in flat spacetime, that is to zeroth order in the derivative expansion, we can simply go to momentum space and perform the inversion with standard techniques. Concretely, both the flat two-point function and the flat propagator have the form Here, p µ is the momentum vector. The determination of the flat part of the propagator from an arbitrary flat two-point function has been carried out in full detail e.g. in [41]. The key observation is that any difference between the flat propagator and the full propagator is, by definition, at least linear in curvature. We can thus define the covariantised flat propagator Crucially, we can choose where we put the scalar propagator functions g i , since any different choice only differs by terms at least linear in curvature. The functions g i can be computed in flat space. To obtain the full propagator, we notice that the product of the inverse of the full propagator G and G 0 is where X is, by definition, of linear and higher order in curvature. Crucially, since the inverse full propagator is just the regularised two-point function, which is known for a given action, X can be computed to the necessary order. Knowing X, we then can calculate the full propagator via where we suppressed the indices. For a fixed order of the derivative expansion, only finitely many terms of this sum contribute to the full propagator. Let us mention that one can of course also choose a different operator ordering. In complete analogy to the above, we can define a tensor Y via so that the full propagator reads In general the tensors X and Y do not agree. Which of the two orderings is more efficient is hard to predict generally, and has to be tested in practice.
The algorithm also works on more general backgrounds, for which one can derive the exact propagator. A prime example is the sphere -the Ricci scalar is finite and covariantly constant, so that (94) holds if we let the propagator functions also depend onR. This has been used in the context of affine gravity in [155].

Commutator rules
The complexity of computations beyond actions linear in the curvature increases extremely quickly. It is thus advantageous to employ tensor algebra packages to perform the necessary calculations to achieve reliable results. A key ingredient for a reliable code is the generic implementation of simplification rules like commutators to a given order. In this section, we will derive such recursive formulas for the commutator of a function of the Laplacian with either a curvature tensor or a covariant derivative. Our focus lies on formulas applicable in a finite order derivative expansion -an extension to the curvature expansion will be presented elsewhere.
Before we specify a detailed commutator, let us consider the following general case. Let f be a suitable function, X is an arbitrary operator and Y is a tensor of arbitrary rank. We are interested in a formula of the form where we want to find an explicit expression of the term indicated by the dots on the righthand side of the equation. To derive the expression, we will formally use an inverse Laplace transform and the Baker-Campbell-Hausdorff formula, In this equation, we use the multicommutator The reason why (100) is useful in a derivative expansion is that the multicommutators increase the order of the expression by at least one. In that way, in a finite order computation, only finitely many terms in this sum contribute.
With the help of the formulas that we prove in appendix A, we can rewrite (100) into the form The usefulness of this formula lies in the fact that it only involves the simple commutator, which is straightforward to implement.
Let us now specify the two cases of commutators that are needed in the derivative expansion. The first case is whenever X is a multiplication operator. The relevant example is that of a curvature tensor, potentially with a number of derivatives acting on it. In that case, we have This can be inserted into (102) and produces (104) Since the commutator increases the order of the expression by at least one derivative, the recursive application of (102) produces only finitely many terms for a fixed order of the derivative expansion.
The second case is when X =D µ is the covariant derivative. In that case, The commutator of two covariant derivatives is related to the Riemann tensor via We thus find This time, the commutator increases the order of the expression by at least two units, so that once again only finitely many terms contribute to any fixed order computation. Repeatedly applying the formulas (104) and (107) then gives the commutator of a function of the Laplace operator to the needed order.

Simplification of tensor expressions of maximal order
Before we finally discuss the application of our framework, we shall point out a way to simplify the calculation of truncated RG flows significantly. This simplification concerns operators which are already of the derivative order that one truncates at, and before such operators are traced. The key observation is that eventually these terms will be contracted with metrics only, and by assumption no higher order terms arise. Because of this, we can replace these operators by combinations of metrics and scalar curvature invariants that respect the symmetries of the term.
To make this concrete, let us first discuss the Einstein-Hilbert case, that is we truncate at the second order in derivatives. At this order, the most convenient way to see the simplification is to transition to the traceless basis, that is all occurrences of the Riemann tensor are replaced by the Weyl tensorC, the Ricci tensor and the Ricci scalar viā and then all Ricci tensors are replaced by traceless Ricci tensorsS and Ricci scalars viā Now since eventually all these tensors must be contracted with metrics only, it is immediately clear that we can drop the terms with Weyl and traceless Ricci tensors since their traces vanish. In other words, terms linear inC andS can only contribute to quartic and higher orders in derivatives, and we can setC Indeed this has been used implicitly in much of the literature on Asymptotic Safety, however with a different view, namely that a special background was chosen. Here we see that the procedure is indeed general and not related to a specific choice of background. At quartic order in derivatives the structure is slightly more complicated. In this case we can replace terms quartic in the curvature, but not those that are linear. Once again it is useful to employ the traceless basis. Since there are only three scalar curvature monomials at this order, in this basis only three combinations of curvatures do not vanish. In particular we find directly thatC since all complete contractions of these terms with metrics vanish. On the other hand, we find thatC µνρσCτ ωκλ → T µνρστ ωκλC αβγδC αβγδ , Here T is a rank 8 tensor constructed from the metric alone, which is too long to be displayed here. These equations can be derived by making the most general ansatz of the appropriate number of metrics, imposing the relevant symmetries, and finally computing some particular contractions to fix any remaining free coefficients. Note that both at quadratic and quartic order, if we neglect boundary terms (which we do in this work), all curvature tensors can be assumed to be covariantly constant. Only at sextic and higher orders, monomials with covariant derivatives appear.
Clearly one can also formulate similar equations in a Riemann basis, however the traceless basis disentangles the invariants maximally. The generalisation to higher orders is also straightforward, although increasingly lengthy. Nevertheless it is also clear that by using these relations the computational complexity can be decreased considerably, since a lot fewer tensor structures arise at any intermediate step of the calculation.

Application to quartic order
We will now apply the machinery introduced in the preceding sections to quantum gravity at the quartic order in the derivative expansion. This entails a theory space with a total of five coupling constants. This section contains a discussion of the action, the flow equations, the fixed point search strategy, the actual fixed point structure of the theory, and a discussion of the topological term.

Action
To fix our conventions, we will first discuss the ansatz for the effective average action. Concretely, this ansatz reads In this, G N is the Newton's constant, Λ is the cosmological constant, and G C 2 , G R 2 and G E are the quartic couplings. All these couplings depend on the FRG scale k, which for better readability we do not indicate explicitly. Moreover, E stands for the integrand of the four-dimensional Euler characteristic, For the discussion of fixed points, we introduce dimensionless coupling constants by a rescaling with the appropriate power of k, so that We will also use an overdot to indicate the derivative with respect to the RG time t, e.g., The action (113) is amended by a gauge fixing term of the form and a corresponding Faddeev-Popov ghost term Note that we deviate from standard convention by introducing the coupling G N also in the ghost action. This is necessary for an exact cancellation of traces as discussed in subsection 3.4. The gauge parameter α will be sent to zero to implement the Landau limit.
The final ingredient to specify is the regulator. Since this has been discussed in detail in section 3, we will not repeat it here.

Flow equations
In this section we present the flow equations for the dimensionless couplings of our system. For convenience, we introduce the dimensionless propagators We also introduce the notation for all regulator shape functions. It is related to the scale derivative of the regulator with some factors pulled out for convenience. The additional dependence onġ comes from the fact that the regulator tensor comes with a prefactor of 1/g, see (86), on which the scale derivative acts non-trivially. The additional factor of 1/g is then cancelled by a factor of g coming from the propagator. This cancellation is already taken into account in the above notations.
For the cosmological constant, we finḋ We did not simplify the prefactor of the integral to illustrate where the individual factors come from: the factor 8πg comes from the left-hand side, the factor of a half comes from the flow equation itself, and the 1/(16π 2 ) comes from the heat kernel. We can also see how the mode count is satisfied: we have a prefactor of 5 for the transverse traceless sector, a prefactor of 1 for the trace sector, and a prefactor of −4 for the combination of gauge variant vector (+4) and ghosts (−8). For identical regulators and at vanishing cosmological constant and higher order couplings, this gives a total of 2 modes that contribute to the flow of the cosmological constant, which is the correct number of physical polarisations of the graviton in d = 4. The flow of the dimensionless Newton's coupling readṡ The terms with more than one power of the propagator come from genuine interaction terms.
In the graviton sector, they are proportional to the higher order couplings due to our regulator choice (86). Next, we will present the flow equations for the fourth order couplings. A general feature of higher order couplings is that if the dimension is at or below the order, their flow equations feature non-integral terms. These stem from positive powers of the heat kernel expansion parameter, which map to derivatives of the function that is traced over, evaluated at zero argument. The flow of the R 2 coupling iṡ (123) For the flow of the C 2 coupling, we finḋ Finally, the flow of the coupling to the Euler term readṡ Since the Euler characteristic is a topological invariant in d = 4, its coupling does not appear in any of the flow equations, except in the scaling term in (125). Intriguingly, the flow of the Euler coupling can be written in terms of the other flows and a term without an integral. The complete set of flow equations (123), (124) and (125) has been obtained for the first time without using a special background [17,82,[162][163][164][165] (thereby neglecting some of the couplings) or expanding in some of the couplings [20,[166][167][168].
We note in passing that the above equations do not reduce to the standard one-loop result of perturbative Stelle gravity if the Einstein-Hilbert part of the action is neglected. There are several reasons for this. First, both the gauge fixing and the regulator are constructed with a non-perturbative setting in mind (meaning that all terms in the action are assumed to be non-vanishing), as they include an explicit factor of 1/G N . Thus to be able to probe the perturbative Stelle limit where G N → ∞, both would need a very careful rescaling. Second, and more importantly, we have used the background field approximation. It is well-known that this can introduce artefacts even into universal one-loop beta functions [149]. To resolve this issue, either the corresponding Ward identities have to be solved, or the limit k → 0 has to be taken. Both options go beyond the scope of the present work. It is however noteworthy that with the standard higher derivative gauge fixing which makes the fourth order kinetic term minimal, the universal one-loop beta functions come out in the usual way, see e.g. [166][167][168]. This leads to the conjecture that such minimal gauge fixings might generally not need Ward identities to obtain such a result. It would be interesting to understand which classes of gauges share this property.

Intermezzo: fixed point search strategy
Before we move on to discuss the fixed point structure of the theory, we will present a strategy to search for fixed points in an extended truncation that are continuously connected to a fixed point found in a smaller truncation. We illustrate the strategy going from one to two couplings, but the method is applicable generically.
Let us assume that we start with a truncation with a single coupling g 1 , and that we have found a fixed point for it, say g * 1 , so thatġ Now we enhance the truncation by adding a coupling g 2 which was set to zero before. That is, the above condition for the beta function of the coupling g 1 now readṡ while the beta function of the new coupling at this point in theory space is in general non-zero, The search strategy is then to change the new coupling g 2 by a small amount, say , and search for a fixed point in g 1 for this new value of g 2 . Let us assume that we find a zero ofġ 1 at the value g 1 , so thatġ The beta function of g 2 will now also have changed. We then repeatedly change g 2 by a small amount, find a fixed point for g 1 , and computeġ 2 . The result is that we getġ 2 as a function of g 2 on a partial fixed point. Clearly, if we find a g 2 such thatġ 2 = 0, we have found a fixed g * λ * g  While this strategy will in general not find all fixed points, one can extend its applicability by starting at any value of the new coupling, search for a fixed point in the old couplings, and start the procedure from there. This gives an efficient search strategy for fixed points in all of theory space.

Fixed point analysis
We now study the fixed point structure of the quartic order of the derivative expansion. To set the stage, we briefly present the fixed point at the quadratic order for our setup. We then add either of the couplings individually, and finally discuss the complete system. In all of the following discussion, we choose the regulator shape function that is a simple exponential regulator.
Einstein-Hilbert truncation The system with g C 2 = g R 2 = 0 has a single fixed point at the coordinates g * = 0.534 , λ * = 0.121 .
The critical exponents at this fixed point are Since the real part is positive, the fixed point is fully attractive. This is the well-known Reuter fixed point, and the results for the critical exponents are in reasonable agreement with results published in the literature [14,15,23,26,28,153,[169][170][171][172][173] when factoring in different regularisation schemes and choices of gauge fixing. C 2 truncation We now add the C 2 term to the system, and use the strategy outlined in subsection 7.3 to search for fixed points. Figure 1 depicts the result of this strategy. It shows the partial fixed point values for g and λ at a range of values of g C 2 , and the value ofġ C 2 at these coordinates. The horizontal dashed line indicates zero, whereas the vertical dashed line is at the value of g C 2 where its beta function vanishes. We find a single fixed point at with critical exponents The value of g * C 2 is positive, so that no new poles arise for positive squared Euclidean momenta in the spin two sector. The inclusion of the C 2 term thus has two effects: first, it makes the formerly complex conjugate critical exponents real, and second, it adds an irrelevant direction, so that the value of one of the couplings can be predicted. Considering the magnitude of the critical exponents, we find a slight to moderate reduction compared to the canonical scaling dimensions, which are 4, 2 and 0, respectively. This is in line with the conjecture of "near-Gaussian scaling exponents" [19,21,43,47,53]. R 2 truncation As a next step, we add the R 2 term to the Einstein-Hilbert system. Applying the fixed point search strategy for positive g R 2 does not yield a fixed point, see figure 2. The slope near g R 2 = 0 indicates that there might be a fixed point for negative g R 2 though. For this case, we however have to flip the sign of the trace regulator, since the coupling of the highest order term in a derivative expansion dictates the sign of the regulator. As a consequence, we introduce a new singularity into the flow, which sits at We are thus confined to the region λ ∈ (−3/4, 1/2), so that all propagators have the correct sign and integrals over the loop momentum exist. This also causes the fixed point search for negative g R 2 to not continuously connect to the case of vanishing coupling. We thus have to start our fixed point search strategy at a finite, negative value of g R 2 , search a fixed point in g and λ, and apply the recursive strategy from that point. The result of this procedure is shown in figure 3, and we indeed find a fixed point at g * = 0.739 , λ * = 0.382 , g * with critical exponents θ 1,2 = 4.44 ± 7.06i , This fixed point is much closer to the singular line λ = 1/2, and the deviation of the critical exponents from the canonical scaling is large. Moreover, the complex conjugate pair has a large imaginary part. Similar signs of instability have been observed previously in f (R)-type truncations [19,21,23,29,32,34,36,39,43,47,[51][52][53]110]. Including higher orders tends to tame these stronger variations.
Complete quartic order Having studied the two quartic terms individually, we now set out for the full system. There are two starting points to initialise our search strategy. Starting from the fixed point with finite g C 2 , (133), we do not find a continuously connected fixed point for positive g R 2 , mimicking the case for vanishing g C 2 . On the other hand, we do find a fixed point when starting from (137). It sits at and has the critical exponents θ 1 = 9.47 , θ 2,3 = −89.6 ± 120i , Judging from the very non-canonical values of the critical exponents, one might conclude that this fixed point is either an artefact of the truncation, or that higher order terms should have a large impact to stabilise the fixed point. It is likely that these extreme values arise due to the fixed point of the cosmological constant lying so extremely close to the singular line λ = 1/2. It is conceivable that higher order terms indeed shift it to smaller values, yielding more realistic critical exponents in the process, but in the end only an actual computation can give certainty about this.

Comparison: regulator without endomorphism
To investigate the stability of the results under changes of the regularisation, we briefly discuss the fixed point structure of the same system where we leave out the endomorphism in the traceless regulator. This entails using∆ instead of ∆ 2 as an argument in (86). Stable fixed points should only be quantitatively affected by this. At the Einstein-Hilbert level, and with only the C 2 term, this is indeed the case. In the former case, we find a fixed point at with critical exponents The critical exponents in this case are very close to the ones found above, see (132). With g C 2 = 0, the fixed is at with critical exponents This is in qualitative agreement with (134), and serves as an estimate of the truncation error. As soon as we include the coupling of the R 2 term, we do not find a fixed point, neither without nor with finite g C 2 . This indicates again the instability mentioned earlier, and we expect more stable results when higher order terms are included.
Let us finally note that the flow of the Euler coupling within this regularisation scheme is not of the simple form (125). It rather has an additional integral, which in particular depends on g C 2 . Taking the simplicity of (125) as a guiding principle, this might be taken as an a posteriori argument for the regularisation choice (86).

The flow of the Gauss-Bonnet term
As mentioned earlier, the Gauss-Bonnet term is a topological invariant. Consequently, the right-hand side of the flow equation is independent of the coupling constant g E in front of it. This implies, absent some unlikely cancellations, that this coupling will never have a fixed point, since its beta function, evaluated on the fixed point of all other couplings, is a constant: This can be seen very easily from the explicit form of the beta function, (125). Depending on the sign of the constant e, g E will go to either a positive or negative infinite value at high energies. This has an interesting consequence for the weight of different topologies in the (Euclidean) path integral. Disregarding the subtleties of the reconstruction problem [174][175][176], if we find that e is positive (negative), spacetimes with a negative (positive) Euler characteristic are enhanced, while spacetimes with a positive (negative) Euler characteristic will be suppressed. This mimics the idea of finite action [177], but the origin is the coupling instead of a divergent curvature invariant.
Due to the extremely simple structure of the beta function, we can make analytical statements about the sign of it at a fixed point of the other couplings. As a matter of fact, at our level of truncation, and with the normalisation condition that the regulators at vanishing argument are ±1 (depending on the sign of the quartic couplings), we find From this it follows that generically, at this level of the truncation, the beta function for the Euler coupling is positive at the fixed point. This suggests that manifolds with a complicated topology contribute most to the Euclidean formulation of Asymptotic Safety. Incidentally, the same conclusion holds for Stelle gravity at the asymptotically free fixed point. We illustrate this with numerical results. At the level of the Einstein-Hilbert truncation, we find e EH = 0.523 , for the fixed point with finite C 2 coupling, we have for the fixed point with negative R 2 coupling, the value is whereas finally, in the full system, we find For a Lorentzian path integral, one might argue that either sign suppresses spacetimes with non-vanishing Euler characteristic, since they "oscillate away". This would give a dynamical mechanism whereby only spacetimes with vanishing Euler characteristic will contribute to the path integral.

Conclusion
In this paper we have set up a systematic framework to study the derivative expansion of nonperturbative renormalisation group flows in quantum gravity. We proposed a set of criteria that well-defined flows and regulators should satisfy. Then we set out to construct a suitable regulator that fulfils these criteria. Geometric considerations guided this search. As a simple example, we discussed the case of vector fields first, and found a natural way to regularise them. Moving on to gravity, we encountered some difficulties which obstruct a similar regularisation. We then used input from the action of General Relativity to nevertheless set up a well-motivated regularisation scheme. From the discussion, it is clear that our approach is applicable to any order in the derivative expansion, including an expansion in form factors.
Having set up the formal structure, we then discussed some techniques that help in practical computations and allow for an efficient evaluation of renormalisation group flows via tensor algebra software. In particular, we presented a general algorithm to obtain the propagator from an arbitrary two-point function. On the more technical side, we provided commutator rules that can be implemented generically in computer code, and discussed simplifications for tensorial terms at the maximum considered order.
We finally applied all these methods to derive and analyse the non-perturbative flow equations at the quartic order of the derivative expansion in quantum gravity. The complete set of flow equations (121) - (125) has been presented for the first time without further assumptions on top of the truncation itself. After a brief generic discussion of our fixed point search strategy, we discussed the resulting fixed point structure at different levels of sophistication. In all approximations, we find an interacting fixed point. The inclusion of the R 2 term introduces previously observed instabilities which are expected to be resolved by the inclusion of higher order terms. Lastly, we discussed the flow of the Euler term, which we found to flow to positive infinity at the fixed point. This indicates that Euclidean Asymptotic Safety is dominated by manifolds with negative Euler characteristic. We speculated that in Lorentzian signature, only spacetimes with vanishing Euler characteristic would contribute. This includes the flat spacetime, and gives a dynamical principle to discard more exotic structures.
Several future directions are available from here. Clearly, to resolve whether the instability of the R 2 term persists upon improving the truncation, the derivative expansion should be extended to sextic order, see [178] for results in the conformally reduced case. We will report results on the full case elsewhere. Another road is the inclusion of form factors along the lines of [139,179], and discuss aspects of unitarity and causality of the theory e.g. in the context of scattering amplitudes [180,181]. This would also allow to compute spectral functions from the fully momentum-dependent background propagators [140].
In any improved truncation, it will moreover be interesting whether the simple form of the flow equation (125) for the Euler coupling remains to be valid. If this would be the case, this would give a constructive proof of our observation on which topologies contribute to the path integral.

A Some commutator formulas
In this appendix we prove some useful commutator formulas. The aim is to rewrite the multicommutator into a form slightly more useful for an implementation in a computer code. In the following, X shall be any operator.
First, we have which we prove by induction. The base case l = 1 is trivially seen to be true. For the induction step, assume that the formula is correct for some l ≥ 1, and calculate In the first line we used the product formula for the commutator, in the second line we used the induction hypothesis, in the third line we relabelled the sum index, and finally we combined all terms into a single sum. This establishes the claim (151). The second formula that we will prove is that The equality between the two sums follows by a relabelling of the summation index. Once again we will prove this formula by induction. The base case l = 0 is true by the definition of the multicommutator. Assume now that the formula holds for some l ≥ 0 and calculate (154) The first step uses the definition of the multicommutator. In the second step, we use the induction hypothesis. Afterwards, we use the definition of the commutator to split the sum into two. The next two steps implement a relabelling of the second sum. Then, we combine the overlapping parts of the two sums, then use a standard identity for the binomial coefficients. In the final step, we recombine the individual summands into the final sum. This proves the formula (153).
Finally, we will combine (151) and (153) and show that We will prove this formula by direct calculation, starting from (153) We started by using (153). We then have split off the first term of the sum, introduced a commutator and calculated one of the sums which cancelled with the first term. We also assumed l ≥ 1 for this in order for the splitting to make sense. We can now use (151) since that formula applies for the involved commutator as the sum starts at one: Here we have exchanged the order of the sums to be able to perform one of them. This completes the proof of the formula (155).