Renormalization Group in Six-derivative Quantum Gravity

The exact one-loop beta functions for the four-derivative terms (Weyl tensor squared, Ricci scalar squared and the Gauss-Bonnet) are derived for the minimal six-derivative quantum gravity (QG) theory in four spacetime dimensions. The calculation is performed by means of the Barvinsky and Vilkovisky generalized Schwinger-DeWitt technique. With this result we gain, for the first time, the full set of the relevant beta functions in a super-renormalizable model of QG. The complete set of renormalization group (RG) equations, including also those for the Newton and the cosmological constant, is solved explicitly in the general case and for the six-derivative Lee-Wick (LW) quantum gravity proposed in a previous paper by two of the authors. In the ultraviolet regime, the minimal theory is shown to be asymptotically free and describes free gravitons in Minkowski or (anti-) de Sitter ((A)dS) backgrounds, depending on the initial conditions for the RG equations. The ghostlike states appear in complex conjugate pairs at any energy scale consistently with the LW prescription. However, owing to the running, these ghosts may become tachyons. We argue that an extension of the theory that involves operators cubic in Riemann tensor may change the beta functions and hence be capable of overcoming this problem.


I. INTRODUCTION
Higher derivative terms play an important role in quantum gravity (QG), as they are required to provide renormalizability (see, e.g., [1] for introduction). However, including higher derivatives leads to the presence of ghosts or ghostlike states, representing a fundamental and, generally, unsolved problem. In the effective field theory approach, it is usually assumed that the higher derivative terms in both the classical action and quantum corrections are taken as small perturbations (see, e.g., Ref. [2] and references therein), regardless of the related inconsistencies, as it was recently discussed, e.g., in [3]. One of the main problems is that the effective approach is applicable only at low energies, while QG is supposed to cover also the UV region. An alternative way is trying to provide a consistent treatment of higher derivative terms in QG at the fundamental level 1 . One of the important aspects, in this case, is the analysis of loop corrections and of their possible implications on the nature of physical degrees of freedom of the theory [5,6].
In the present paper, we consider the higher derivative operators in the classical action as a viable nonperturbative completion of Einstein's gravity to achieve a super-renormalizable theory. In other words, the higher derivative operators here are not small in comparison to the Einstein-Hilbert action, but actually dominant in the ultraviolet regime, contrary to the effective approach. A usual proposal to overcome the nonrenormalizability of Einstein's theory is to introduce higher derivative terms, with four or more derivatives of the metric, improving the UV convergence of the theory at quantum level in the perturbative loop expansion.
properties of the theory. As it is usually done in gravitational contexts here we mean the stability of small perturbations of equations of motion (EOM) analyzed around a given classical background, so here we quadratize the action and for example look for spacetime evolution of perturbations and can find Lyapunov exponents of the stability (both for the perturbative and nonperturbative levels of changes of the background metric [49][50][51]). This notion of stability is, of course, different than the issue of stability of ground states of the quantum system describing vacua of these quantum field theory (QFT) models. As it is well known the latter concept is investigated by truly nonperturbative methods, for example, of vacua tunneling and spontaneous symmetry breaking in the ground states of these QFTs.
The unitarity and stability issues have been considered also in the polynomial higher derivative QG. In particular, in the recent papers [15] and [16], it was shown that the QG model with six derivatives is tree-level unitary if the massive poles of the classical propagator form complex conjugate pairs. The power counting in these models guarantees superrenormalizability and, therefore, the polynomial QG models also have universal running of all parameters, in the same sense as the nonlocal models have.
The difference between nonlocal and polynomial theories of QG is that, in the nonlocal case, there is no available reliable covariant technique for the one-loop calculations. On the other hand, in the polynomial models, the one-loop divergences can be obtained using Barvinsky and Vilkovisky universal traces [52], albeit such calculations may be rather cumbersome and technically nontrivial. Until now, the list of the known beta functions included the ones for Λ and G, derived in [10] and [53], respectively. In the present work, we report on the quantum calculations of the remaining beta functions in a minimal six-derivative model of QG, which includes the versions with real and complex poles corresponding to both normal particles and ghosts. The last case is the one which was used in [15] for the Lee-Wick analysis of the unitarity. Indeed, the beta functions obtained in what follows form, together with the previous ones [53], the full set and, therefore, enable one to explore how the running affects the positions of the poles of the propagator of metric perturbations on the complex plane. From this perspective, the derivation of the full set of equations describing the running is the necessary step in describing the theory at the quantum level.
From a technical point of view, the one-loop calculations in super-renormalizable QG present more difficulties compared to the ones in the four-derivative renormalizable models [35][36][37]. The level of complexity of such a calculation depends on the type of one-loop counterterms. The counterterm for the cosmological constant is actually easy to obtain [10]. However, the derivation of the divergence linear in the scalar curvature requires really big efforts and was done only recently in [53]. In the present work, we make the next step, calculating the simple one-loop divergences for the four-derivative sector. As a result, we arrive at the beta functions for the Weyl-squared, R 2 and Gauss-Bonnet terms. The calculation is really cumbersome and it is done for the simplest possible six-derivative theory without "killer" terms in the classical action, which here are third powers of the generalized curvature tensor R. Even in this simplest case, the intermediate expressions are too large for the presentation, hence they will be mostly omitted. Similar computations in four-, six-and general higher derivative gauge theory were performed in [54][55][56][57][58].
The work is organized as follows. In Sec. II we briefly review the renormalization in the six-derivative model of QG and show that the divergences in this theory do not depend on the gauge-fixing. Section III is devoted to the derivation of one-loop divergences, while in Sec. IV we introduce the renormalized Lagrangian, the beta functions for the four-derivative operators and the related renormalization group equations. In Sec. V we apply our results to the six-derivative Lee-Wick quantum gravity solving exactly the renormalization group equations for all running coupling parameters. Later in Sec. VI we discuss the issue of asymptotic freedom in this model. Finally, in Sec. VII we draw our conclusions and give a little discussion of additional related issues. In the special appendix A we present a detailed analysis of the Wick rotation in local higher derivative quantum field theories.

II. POWER COUNTING AND GAUGE-FIXING-INDEPENDENCE OF DIVERGENCES
The action of the theory that we consider in this paper has the form where the density of Lagrange function reads 4 Here C µνρσ is the Weyl tensor and E 4 is the Gauss-Bonnet term, which is the Euler characteristic in d = 4, In what follows we will use the pseudo-Euclidean notations and |g| denotes the absolute value of the square root of the metric determinant. The coupling ω κ is related to the Newton constant G, while ω Λ to the cosmological constant Λ. The theory with the Lagrangian (2) is the simplest one that describes the general form of the graviton propagator around flat spacetime in four spacetime dimensions with six derivatives. It is also the simplest QG model admitting, besides the massless graviton, only Lee-Wick complex conjugate poles at the classical level [15]. Let us note, from the very beginning, that there are two remarkable special cases in the theory (2). In order to have a nondegenerate classical action one needs that both coefficients ω C and ω R should be nonzero and the quantum calculations reported in the next section correspond only to this kind of model. Correspondingly, in the special case of ω C = 0 and ω R = 0, the theory has the propagating spin-two mode with four derivatives and the propagating spin-zero mode with six derivatives, while interaction vertices have six derivatives in both special and generic theories. For another special version, with ω C = 0 and ω R = 0, the situation is opposite, but according to the power counting arguments [10] in both special cases the theory is nonrenormalizable.
In the case of both ω C and ω R nonzero, the theory (2) is the simplest example of a super-renormalizable model of QG. The power counting procedure results [10] in the expression: where p is the number of loops of the diagram, d is the number of derivatives of the metric on the external lines or, equivalently, the mass dimension of the divergences, D is the superficial degree of divergences of the diagram. For the logarithmic divergences, D = 0 and we can see that the one-loop divergences have d = 4, while higher loops give less derivatives in the counterterms. This result has deep consequences for the general structure of renormalization in the model (2). Let us list here these aspects of the theory.
i. There are no UV divergences with d = 4 at higher loop orders than the first one. Thus, the beta functions for the four-derivative terms, which can be derived at one-loop order, are exact and hold to all perturbative orders in the loop expansion. The situation is different for the Einstein term with d = 2, which gains divergent contributions (and hence contributions to the beta functions) for p = 2 and for the cosmological constant term which gains such contributions for both p = 2 and p = 3. Starting from p = 4 (four-loop order) the loop corrections are finite, i.e., the theory is super-renormalizable.
ii. The terms with four, two and zero derivatives in the action (1) do not affect the counterterms with d = 4. This means that the practical calculations of the d = 4 counterterms in this theory can be performed for the reduced model with the action On the other hand, the d = 4 counterterms may be affected by inclusion of extra O R 3 -type terms into the action. These terms are not necessary, e.g., for renormalizability and should be regarded as nonminimal. The practical calculations in the next section were performed without these terms.
iii. The power counting (5) shows that the classical action and the divergences have different numbers of derivatives of the metric. This feature implies that the theory is super-renormalizable, but also results in the universality of quantum corrections that we mentioned in the Introduction. According to the general statements about gauge-fixing and parametrization dependence of the effective action (see, e.g., [60] for a general proof and [61] for the reduced simplified version valid for the one-loop approximation), both kinds of dependencies disappear on the mass shell. If we describe the ambiguity in gauge and parametrization by the parameters α i and use the logic of [36] (and, more explicitly, in [29,40]), then the difference between the effective actions Γ for the two choices of these parameters, namely α i and α 0 i , can be written as follows: where f µν = f µν (g κλ , α i , α 0 i ) is an unknown function of the metric and parameters α i and α 0 i and If we are interested in the divergent part of the effective action, f µν is a local expression which is at least of the first order in the loop expansion parameter . Thus the dimension of the right-hand side of (7) is at least equal to the dimension of ε µν . Furthermore, the equations of motion ε µν can be expanded into a power series in . At the lowest level of this expansion, one meets classical equations of motion, hence ε µν is at least of the sixth order in the metric derivatives. At the same time, according to the power counting (5), the divergent part of the left-hand side of (7) contains only four derivatives of the metric. Therefore, the divergent parts of the both sides of Eq. (7) vanish implying that f div µν = 0. The conclusion is that the right-hand side of (7) is a finite expression and, hence, the divergent part of the effective action Γ does not depend on the gauge-fixing or parametrization ambiguities.
All in all, we have shown that the beta functions of the theory (2) are defined in a universal way, that means they are free of ambiguities which are typical for QG based on Einsteinian general relativity or on the renormalizable fourderivative theory of gravity. As it was already mentioned above, the derivation of zero-and two-derivative divergences has been previously done in Refs. [10] and [53]. In what follows we will complete the set of beta functions for the theory (2) by deriving the exact beta functions for the four-derivative coefficients θ C , θ R and θ GB . Without loss of generality the calculation will be performed for the reduced model (6).

III. ONE-LOOP CALCULATION IN SIX-DERIVATIVE MODEL
Except for a larger volume of algebra, the one-loop calculations are performed using the same general algorithm developed in the super-renormalizable models [53]. Thus we can skip a great part of the explanations and give only specific details of the calculation of the UV-divergent four-derivative terms.
The background field method is defined by splitting the metric Since the result does not depend on the gauge-fixing, the parameter β in a linear gauge-fixing condition χ µ , can be chosen in the most simple "minimal" way, as we will show below. The same concerns the parameters α and γ in the weight operatorĈ = C µν , This operator is already in a self-adjoint form and it is originally defined when acting on a covariant vector χ ν . Together with the gauge-fixing χ µ the operatorĈ defines the gauge-fixing action [36], The action of the Faddeev-Popov complex ghosts has the form where the bilinear part depends on χ µ and on the generator of gauge transformations, Let us give a few more details concerning the choice of the gauge-fixing parameters. The bilinear form of the action is defined from the second variational derivative, where the first term includes six-derivative terms and O(R) is the rest of the bilinear form with four of less derivatives and nonvanishing powers of curvature R. The dimension is compensated by the powers of curvature tensor and its covariant derivatives; we can also commute some covariant derivatives and trade them to curvature in this case 5 . The corresponding terms are very bulky in the present case and it is not reasonable to include them in the paper. The six-derivative part of theĤ operator, after adding the gauge-fixing term, has the form For the sake of brevity, in this expression we left implicit symmetrization inside the pairs of indices (µ, ν) and (ρ, σ).
It is easy to see that to make the operatorĤ minimal, one has to choose the following values of the gauge-fixing parameters [53]: As we have explained in the previous section, this choice does not affect the one-loop divergences in superrenormalizable QG. Thus, we assume (17) as the most simple option.
With these values of the gauge-fixing, Eq. (16) becomes H µν,ρσ lead = ω C δ µν,ρσ + kg µν g ρσ 3 with δ µν,ρσ = 1 2 (g µρ g νσ + g µσ g νρ ) , where the coefficient k is In the expressions for β in (17) and for k in (19), the denominator blows up near the point where ω C = 3 2 ω R . This point corresponds to the leading six-derivative term of the gravitational action in (2) of the form R µν R µν . However, the final results show no singular dependence on the 4ω C − 6ω R . Therefore, the results in this special case can be regarded as a smooth analytic continuation from the generic case when 4ω C − 6ω R = 0.
The DeWitt metric in the internal space of metric fluctuations h µν has the form where The inverse DeWitt metric can be found in the form G −1 κλ,µν = y 1 2 g κµ g λν + y 1 2 g κν g λµ + y 2 g κλ g µν = y 1 δ κλ,µν + y 2 g κλ g µν , such that The coefficients can be easily found to be In the expression for y 2 we explicitly see that limits ω C → 0 and ω R → 0 are singular. We find that new operator H ′ is in the standard minimal form For the divergent part we have Tr logĤ ′ = Tr logĤ. Now we have all necessary elements to write down the general formula for the one-loop contribution to the effective action of the theory, The calculation of the divergent parts of the first two expressions in (26) is pretty much standard [53], with the use of the generalized Schwinger-DeWitt technique [52]. For this reason we shall skip most of the technical details and will just comment on the derivation of the last term.
In order to derive the O R 2 terms in the divergent part of effective action (26) coming from the operatorĈ in (11), we remind the reader that both degenerate γ = 0 and nondegenerate γ = 0 versions of it are possible, according to Eq. (17). These two cases should be considered separately. For the sake of simplicity, we focus only on the nondegenerate case. It turns out to be useful to introduce another conjugate operatorB, whose contribution can be found in [52], One can choose the parameter b in such a way that the product of two operators is a minimal sixth-order operator, The value of the parameter b is obtained from the last condition (28), namely To cancel the nonminimal terms we need Starting from (28), we find where U and V tensors are with nonzero powers of curvature. The last expression, in the right-hand side, can be elaborated using the universal traces in [52], while the formula for the divergent part of Tr log B νλ can be found in the same paper. The ellipsis in (31) denotes terms cubic in curvature or with background dimensions 6, hence not contributing to divergences. As we have already mentioned above, the intermediate expressions for our calculations are too cumbersome for the LaTeX file, so they will not be presented here.
Let us note that the computation of Tr logĈ was checked using three methods. Since theĈ operator is a nonminimal four-derivative differential operator with vector indices, then the computation of its trace of the logarithm is a bit tricky and one has to be careful. This is why we performed here additional verifications of our partial results for this trace. All the methods consist basically of creating new operators, with higher number of derivatives. By selecting adjustable parameters, those could be put into a minimal form. This reduction was achieved by an operatorial multiplication by some two-derivative spin-one operatorB containing one free parameter, e.g., b as in (31).
In the first method, we multipliedĈ by another operator from the right side. In the second case, the multiplication was from the left. And, in the third method, we used the symmetric form of multiplicationBĈB, which is important for the self-adjointness property of the resulting operator. In this case,B has the form of a two-derivative operator, whose trace of the logarithm is known (and was also checked below). We remark that in the first two methods the resulting operators were six-derivative ones, while in the last one it was an eight-derivative operator. All three methods of computation of Tr logĈ agree for the divergent terms quadratic in curvatures, appearing in the form of UV divergences of the theory (R 2 , R 2 µν and R 2 µνρσ ). Similarly, the computation of Tr logB was checked using three analogous methods. We used multiplication from both sides by the operatorÂ and also the symmetric formÂBÂ, whereB is a two-derivative operator, whose trace of the logarithm is searched for. Here,Â is a two-derivative nonminimal spin-one vector gauge operator, whose trace of the logarithm is well known [36,52]. Again, all three methods agree for terms quadratic in curvatures.
A. The results for the traces and effective action Skipping the intermediate details, we present the results for the functional traces. It proves useful to use the covariant cutoff regulator L [52], related to the dimensional regularization parameter ǫ by [52,62] log where n is the regularized dimension of spacetime. The contribution of the main operatorĤ to the divergent part of the expression of our interest (26) has the form where we introduced the useful notation Some observation is in order at this stage. From the quantum field theory arguments discussed above we cannot expect that the limit ω C → 0 (x → 0) will be continuous because when ω C = 0 exactly, the number of degrees of freedom of the theory changes. But still for every δ > 0 one might expect that the divergent part of the effective action is continuous and bounded for x > δ. Therefore, all the terms proportional to 1/x in (33) should cancel out in the final expression. This can be seen in (37) below as an effective test of the correctness of calculations 6 . The contribution of the bilinear ghost sector is related to the expression And finally, the contribution of the weight operatorĈ from (11) has the form Summing up the expressions (33), (35) and (36) according to (26), we arrive at the final result for the divergent contribution to the quantum effective action for the operators C 2 , R 2 and E 4 , where we used a more physical basis defined by Eq. (4). It is remarkable that in this expression there are no singular terms of the type 1/x, as it was expected. Note that these terms are present in the intermediate expressions for (33), (35) and (36), hence we conclude that their cancellation in the sum is a partial verification of the result. On the other hand, the smooth limit x → 0 does not mean that the result (37) remains valid for x = 0. As we explained above, in this case the form of the operatorĤ changes dramatically and the form of the divergences is supposed to change too. The calculation in this case can be done, e.g., using the method described in Chapter 8 of [61]. At the same time, the theory with ω C = 0 is nonrenormalizable and, hence, such calculation does not look sufficiently interesting.
The expression (37) confirms the expectations based on power counting, general covariance and locality. The result is covariant, local and has four derivatives of the metric, in agreement with Eq. (5). Furthermore, all the coefficients of the three divergent terms of (37) depend only on the dimensionless ratio (34) of the terms with six derivatives of the classical action, but not on the lower derivative terms, as it was in the case of the cosmological constant and the linear in R divergences [53].
Let us say that regardless of the simplicity of the formulas (33), (35) and (36) and of the final result (37), the intermediate calculations were huge and not only because of the size of the expressions for the blocks (where Mathematica [63] and its special package xTensor [64] for symbolic algebra manipulations provided a useful assistance), but also due to the complexity of all the steps starting from the quadratic expansions. The correctness of the calculations has been checked in several ways, as briefly described below.
First, the expression for the operator of the second variational derivative of the action was verified. The divergence of the second variational derivative operator (Hessian) with respect to gravitational fluctuations h µν , from each covariant term in the gravitational action is separately zero, namely This formula was explicitly checked to the order quadratic in curvatures and up to total of four derivatives. Another verification is that, in the total divergent effective action, we checked a complete cancellation of terms with singularities in y = 4ω C − 6ω R variable, namely of terms proportional to 1 y and 1 y 2 . This is a nontrivial cancellation, as the corresponding terms emerge in both Tr logĈ and Tr logĤ.
Finally, in order to control our Mathematica [63] code, the technically similar computation in the four-derivative QG was performed. We easily were able to reproduce results about one-loop UV divergences there [37,39]. We found a complete agreement for all the coefficients and the same dependence on the parameter x as in the original papers. This calculation was the last of the tests.
There is one interesting detail concerning the expression for the divergences (37) in the minimal six-derivative theory. It is easy to note that the coefficient of the C 2 term has a general form with x defined in (34) and A 0 and A 1 constants. On the other hand, the other counterterms (i.e., R 2 and E 4 ) come with constant coefficients, independent of x. This pattern is different (one can even say opposite) from the situation in the four-derivative QG [36,37,39], where only the divergence proportional to R 2 depends on the analog of (34), i.e., on x 4−der = θC θR , according to the general form At the same time, e.g., the C 2 divergence in the four-derivative QG comes with the universal (and gauge-fixing independent [40]) coefficient. The mentioned difference is noticeable, regardless we cannot explain it using some general principles.
One could also analyze a special value of the fundamental ratio x of the minimal six-derivative theory (2), which makes the C 2 sector of UV divergences completely finite. This value is We remind the reader for comparison that in the case of quadratic gravity, the special values for x 4−der , which made the coefficient of the R 2 term in divergences vanish, is For the sake of completeness we also remind the reader here of the divergent contributions to the quantum effective action with the Einstein-Hilbert [53] and the cosmological constant operators [10], The sum of (41) and (37) is the full set of the relevant one-loop divergences in the six-derivative model (2). The remaining counterterm R being a total derivative is, in fact, the most difficult to calculate. Also, this makes not much sense to derive this term in the nonconformal theory, because it is equivalent to a finite R 2 -type contribution and the overall R 2 term is the subject of renormalization conditions.

IV. BETA FUNCTIONS AND RENORMALIZATION GROUP
Starting from the divergences derived in the previous section, we can now derive the beta functions of the theory. From (37) and (41), the total expression for the divergent part of the effective action can be presented in the form The derivation of the beta functions is a standard operation, explained, e.g., in the recent textbook [1]. Thus, we give only the final result. The running of all the relevant parameters in the six-derivative QG is described by the following renormalization group equations: It is worth noting that the one-loop beta functions for four-derivative terms (C 2 , R 2 , E 4 ) listed above are exact and universal in the theory (2), as we have demonstrated in Sec. II, while the beta functions for ω κ and ω Λ receive also corrections at two-loop and two-and three-loop orders respectively. However, all these beta functions (except the β Λ ) can be modified as the result of inclusion of additional terms into the initial action of the theory. For example, this happens if we add terms that are cubic in Riemann tensor (so called killers). These terms do not spoil the superrenormalizability of the theory. In particular, adding the term (which can be shown to be a sum of total derivatives and O R 3 terms [10]) is supposed to affect the renormalization group equations (43), (44), (45) and (46). For example, as computed in [53], the beta function for the parameter ω κ with the additional contribution of the ω GB · GB 1 term in the Lagrangian (2) reads Moreover, in this case, it is expected that beta functions β C , β R and β GB exhibit up to quadratic dependence on ω GB . In this section, we do not introduce this or other similar structures and consider the quantum theory based on the minimal model (2).
In the remaining part of this section, we solve the renormalization group equations (43), (44) and (45) for the four-derivative terms, while the equations for the Newton constant and the cosmological constant will be considered in the next section for the special case of Lee-Wick quantum gravity.
Let us start with (43). In this case the solution is where the condensed notation with the logarithmic RG-time t = log(µ/µ 0 ) was used. In the four-derivative theory, the perturbative coupling constant is λ C (see, e.g., [1]), where λ C = − 1 2θC . In this case, the stability of the theory requires λ C > 0, so θ C (t) < 0. Moreover, to have asymptotic freedom in the UV regime, the signs of θ C (0) and β C should be the same.
It is worth noting that this situation is common for the higher derivative models since in this case there are, typically, different degrees of freedom with different masses [36]. When these particles are separated by using auxiliary fields (see, e.g., [3] or [1] for a simplified example), the constants of rescaling for the propagators tend to zero in the case of asymptotic freedom for mutual interactions between these different degrees of freedom.
The solution to Eq. (44) reads Similarly to the Weyl-squared term, in a purely four-derivative theory this would be related to the coupling ξ by θ R = − 1 ξ and the stability requires ξ > 0. In this case, asymptotic freedom requires θ R (0) to be negative, something that cannot be easily established in the present case.
Finally, the solution of (45) is In the same sense as described above, for the asymptotic freedom we need θ GB (0) > 0.
Let us note that one can attribute a physical meaning to the running described above only in the UV (high energy regime) where the MS-based renormalization group reflects the physical running of effective charges that depend on the energy scale. At lower energies the logarithmic running does not take place due to the Appelquist-Carazzone decoupling theorem [65]. For higher derivative QG this theorem has not yet been tested, but at the semiclassical level there are solid results in this respect [66,67] obtained by means of diagrams and by the extended Schwinger-DeWitt technique [68,69]. Taking this into account, it is reasonable to assume that the MS-scheme-based running in (50), (51) and (52) takes place for the energy scales above a threshold defined by the magnitudes of the masses of the higher derivatives modes, which are actual degrees of freedom (i.e., beyond the massless graviton) that are present in the theory.
Indeed, a sufficiently intensive running can change the value of a relevant parameter even on a restricted interval of energy scales. In this respect, the most interesting is the running (51) of the coefficient θ R of the term R 2 in the action. However, it is easy to see that, according to the solution (51), the running of the coefficient θ R cannot provide its value of the order 10 8 − 10 9 , required for the phenomenologically successful Starobinsky model of inflation [45,70]. The perturbative running in super-renormalizable QG is logarithmic and the beta function is a parameterindependent constant. Thus, exactly as in the case of other quantum corrections (see, e.g., [46] for the review and further references), this running cannot change the value of the parameter θ R by many orders of magnitude starting, e.g., from the value of order one.

V. SIX-DERIVATIVE LEE-WICK QUANTUM GRAVITY
The theory proposed and studied by Stelle in [7] shows good quantum properties like renormalizability and asymptotic freedom [36,37], but the presence of a ghost instability at the classical level makes the theory nonunitary in its original quantization based on the Feynman prescription [7]. However, recently a new quantum prescription, based on the Cutkosky, Landshoff, Olive, and Polkinghorne [20] approach to the Lee-Wick theories [21,22], was introduced by Anselmi, and Piva [17,18]. This new prescription allows one to tame the ghost instability of Stelle's theory. In this way, the unitarity problem is definitely solved at any perturbative order in the loop expansion [19]. At the classical level the ghost (what Anselmi and Piva called "fakeon" because it can only appear as a virtual particle) is removed from the spectrum obtained by solving the equations of motion for the fake field by means of advanced and retarded Green's functions and by fixing to zero the homogeneous solution [71,72]. This is the classical equivalent of removal of ghosts in the quantum theory from the spectrum of allowed asymptotic states.
The described prescription is very general and can be applied to real as well as complex poles implying ghosts (fakeons) but also to normal particles. In particular, it can be applied to make perturbatively unitary the theory proposed in [15,16] and named "Lee-Wick quantum gravity". This class of theories is based on the general higherderivative action proposed in [10] (with more than four derivatives). Thus, we arrive at a large class of superrenormalizable or finite and unitary higher-derivatives theories of quantum gravity. In order to guarantee tree-level unitarity, the theory in [15,16] has been designed to possess only complex conjugate massive poles in the propagator, besides the massless graviton. It is worth mentioning that Stelle's theory [7] with the Anselmi-Piva prescription [17,18] is the only strictly renormalizable theory of gravity. On the other hand, the theories proposed in [10,15] represent an infinite class of super-renormalizable or finite models for QG. An important issue of uniqueness of these theories will be addressed in a possible future project.
The minimal six-derivative Lee-Wick model is obtained starting from (1) and (2), by fixing the parameters of the action to provide a suitable spectrum of the theory, i.e., one has to ensure the presence of the usual real graviton field and a pair of ghosts with complex conjugate poles, with the complex square of the mass. The action of the theory, including the cosmological constant term, can be written in the form where we also added for our convenience a generalized Gauss-Bonnet term GB 1 from (48), which will also be useful later in this section. Due to this addition the model ceases to be minimal. Let us start by considering the classical case without cosmological constant and show that the spectrum of such a theory contains only complex conjugate poles and the graviton. Since E 4 and GB 1 terms do not contribute to quadratized action for the graviton field in Minkowski spacetime (so we set ω Λ = 0), the propagator, up to gaugedependent terms, in Fourier space and in four spacetime dimensions, reads where P 2 and P 0 are spin projectors (see, e.g., [1]) and We remark that the above form of the coefficient in front of the spin-0 part, proportional to the projector P 0 , of the propagator, although originally it was obtained in the harmonic gauge, is valid generally for any gauge choice. Actually, functions H 2 (k 2 ) and H 0 (k 2 ) do not show any dependence on gauge-fixing parameters. Hence also masses of additional spin-0 or spin-2 modes are completely gauge-fixing-independent [1], both on the flat spacetime as well as around (A)dS backgrounds. If we define z = k 2 , the location of the zeros for both H 2 and H 0 is Thus, the conditions to have two pairs of complex conjugate poles are It is evident from the above constraints (59) and (60), that for ω κ > 0, we have to assume ω R < 0 (below we will use ω R = −|ω R |) and ω C > 0. At the quantum level, the relations (59) and (60) may change because of the running with the energy scale µ and some coupling constants like θ C , θ R , ω κ above, acquire nontrivial t-dependence. We remark that in an hypothetical ideal situation one should look for the full answer in the full quantum effective action and from this functional, one should read the positions and residues of the poles of the propagator. These would then correspond to the effects of true quantum propagation of dressed particles. The effective action functional Γ for this task should be understood on the tree-level. However, as it is well known the computation of the effective action is a very difficult task. Even at the one-loop level all its finite terms are much more complicated than just the technically demanding calculation we present in this paper. Therefore, our strategy is to look at the RG-improved form of the tree-level action of the theory. In this approach, we substitute tree-level classical values of couplings by their quantum running analogs obtained by integrating the beta functions that we have computed in the first part of the paper. This constitutes a first step towards taking the impact of quantum effects on the propagation of modes and on the spectral properties of the quantum theory since we do not have at the moment the knowledge of the full effective action and presently its full form is beyond our computational reach.
From such an RG-improved classical action of the theory we derive the propagator and analyze the structure of poles and their residues for all dynamical particles of the model by standard procedure (treating the RG-corrected action as given at the tree-level). We also want to remark that we expect that the quantum perturbative corrections (for example at the one-loop level) will not change too much the positions of the poles compared to the classical values. Hence, according to the perturbative approximation very likely a classical Lee-Wick theory will also remain such after the inclusion of quantum corrections.
Such an RG-improvement procedure for the action and its couplings is standard in QFT without gravity and around flat spacetime. We strongly think that the case of gravity, or of 6-derivative QG in particular, should not be more special in this respect. Gravitation shall not be treated like completely exceptional QFT, it is simply a quite peculiar QFT with a special particle content and a special set of symmetries. One can easily understand this by viewing gravity from the higher spin perspective. Models with higher spins are perhaps the most profoundly generalized QFTs, while theories of pure QG are special cases of former ones. Our philosophy is to use and apply the same methods which were already successfully applied in the case of QFT on a flat background for nongravitational interactions, but this time to quantum gravity. We think that gravity is not much different from other interactions as considered in the particle physics framework and could be treated on the same footing as they are and could be, for example, completely understood in the Feynman way of viewing gravity [73] as another spin-2 field (on flat background) with self-interactions dictated by a suitable gauge principle.
Furthermore, in order to verify the stability of the spectrum, we must also consider the running of the Newton constant G related to the ω κ coupling constant. Hence, we should check under what conditions on the free parameters ω C , ω R , initial conditions θ C (0), θ R (0) and ω κ (0), the ghosts remain forming complex conjugate pairs at any energy scale. Here we set ω GB = 0 in (53). Making explicit the dependence on t, (59) and (60) become Using the conditions discussed in the previous section, the requirement of having an asymptotically free theory implies θ C (0) > 0 (for β C > 0) and θ R (0) < 0 (for β R < 0) consistently with (50) and (51) respectively. Hence, for θ R (0) < 0 we can rephrase (50) and (51) as Furthermore, if we assume θ C (0) > 0, we need the following condition to be satisfied to ensure asymptotic freedom: Finally and most importantly, we should ensure that there are no tachyons in the spectrum, namely that the real parts of the two roots in (57) and (58) are negative according to our signature of the metric. Unfortunately, this is not possible, because in (57) θ C (t) and ω C are both positive as well as in (58) θ R (t) and ω R are both negative. In the spin-zero sector, we could flip the sign of ω R , but in this case due to (60), the ghost would become a real particle (see, e.g., [74] for a detailed discussion 7 ). This problem cannot be solved in the framework of the minimal model (2). Thus, the only way out is to extend the Lagrangian (2) by adding new operators, like GB 1 in (53), cubic in the Riemann tensor and/or its contractions, to make the beta function of the θ R coupling positive. In the latter case, we can achieve asymptotic freedom and condition of no tachyons at the same time, because (51) will be replaced by a new running with a positive beta functionβ R > 0 and with a positive initial condition θ R (0) > 0, namely with a positive overall sign. Much simpler is to eliminate the tachyon for the spin-two sector. As stated above, in (59) ω C cannot be negative if we want to have a complex pole. However, we can fix x to make zero the beta function β C or to flip the overall sign in (50) and choose, at the same time, θ C (0) to be negative. In this case, we guarantee that θ C (t) < 0, so the real part of the mass square parameter in (57) is negative. Also, under these conditions, the group velocity is smaller than the speed of light [74] and macrocausality is satisfied [75].
In the six-derivative theory, the cosmological constant ω Λ runs with the energy and the assumption of being in Minkowski spacetime is not valid anymore. Therefore, we must study the stability of the theory in (A)dS backgrounds. To this end, we shall follow the general analysis in [76][77][78][79]. At the quantum level, the conditions (61) and (62) are not sufficient to guarantee the absence of ghosts and we have to study the perturbative spectrum of the theory in (A)dS.
In order to expand the action (53) to the second order, in the graviton fluctuation h µν , around the onshell (A)dS background (ḡ µν ), we use the following York decomposition of the graviton field [76,80]: where the spin-two transverse and traceless fluctuation h ⊥ µν contains 5 degrees of freedom because it satisfies ∇ µ h ⊥ µν = g µν h ⊥ µν = 0. The transverse vector A ⊥ ν , satisfying ∇ ν A ⊥ ν = 0, accounts for 3 degrees of freedom. Finally, B and h are two real scalars. However, A ⊥ ν automatically drops out of the second variation of the action of the theory (53) due to the diffeomorphism invariance and out of the two scalars only the following combination φ = B − h appears there. Here, again, we set ω GB = 0 in (53). We end up with the following second order variation of the action around an (A)dS background parametrized by the cosmological constant Λ (according to the relationR µν = Λḡ µν ) [76][77][78][79], where we introduced the following two polynomials (being analogs of (55) and (56) in the case Λ = 0), and the canonically normalized fields In order to have the same stability properties (at the linear level) in the six-derivative theory as in Einstein's twoderivative gravity, we require the two polynomials H 2 and H 0 to have only complex conjugate roots. This requirement corresponds to the Lee-Wick prescription in the presence of the cosmological constant Λ. Introducing again the variable z = − = k 2 (here by =ḡ µν ∇ µ ∇ ν we mean the d'Alembert operator in the (A)dS spacetime), the roots of (69) and (70) are where Now we assume ω R < 0 and ω C > 0 according to the conditions in (59) and (60) in Minkowski spacetime. Moreover, we take θ R (0) < 0, in order to achieve asymptotic freedom in the UV regime. The running onshell cosmological constant parameter Λ(t) is defined by (when we use effective equations of motion from Einstein's gravity) In the full quantum theory, for consistency we end up with the situation that the scalar curvature of the onshell (A)dS background upon which quantum fluctuations are considered, is proportional to the running cosmological constant parameter Λ(t) according to the relationR = 4Λ(t). This signifies that the background curvature depends on the energy scale µ used to probe it by means of fluctuations. For consistency, the classical conditions for Lee-Wick particles in (A)dS ∆ 2 < 0 and ∆ 0 < 0 with (74) and (75) have now to be considered with the running parameter Λ(t) instead of Λ. Therefore, the conditions (61) and (62) are now replaced by In the case of the modified beta function for ω R as in (66), it proves useful to make the replacement θ R (0) → −θ R (0). Thus, θ R (0) must be chosen positive to achieve the asymptotic freedom in agreement with (66).
Integrating the renormalization group equations (46) and (47), we finally get the following general solutions: where θ C , θ R stand for the initial values θ C (0), θ R (0). The physical (running) cosmological constant given in (76) is proportional to the ratio of ω Λ (t) and ω κ (t), which are cubic and quadratic polynomials of t respectively, namely we have where the above coefficients a κ , b κ , a Λ , b Λ and c Λ are just nonrunning constants and ω κ and ω Λ are initial values (at t = 0) of running couplings ω κ (t) and ω Λ (t) respectively. Therefore, in general, the cosmological constant grows linearly in the ultraviolet regime, i.e.
However, the space of parameters for couplings in front of cubic operators, like ω GB , is certainly large enough to make vanish the linear contribution in t to the cosmological constant Λ(t) that now can approach a constant value in the ultraviolet regime. One can see that this is not possible just with playing with values of ω C and ω R alone since the trinomial does not have real solutions (its discriminant is negative). Note that the two parameters ω R and ω C do not scale with the energy, i.e. they do not run. One can envisage three different cases for the UV-limiting behavior of the physical cosmological constant. In the first (worst) scenario, the term proportional to t 3 in (80) is not canceled and the Λ(t) goes to infinity in the UV regime and hence the spacetime becomes fuzzy or with a foamlike structure [36,37]. This can be avoided by inclusion of cubic operators with tuned coefficients. In the second case, the cosmological constant goes to zero in the high energy regime and we end up with the Minkowski vacuum in the ultraviolet. In the last scenario, if we do not select the initial conditions entailing the second case, we have to face (A)dS vacuum at short distance with some finite radius of curvature.
By inclusion of the cubic in curvature operators, one can imagine a situation when the Lee-Wick conditions (77) and (78) are always satisfied, namely the poles of ghosts stay in a complex pair at any energy scale from IR to UV. In the second case considered above, with the initial condition Λ(0) = 0, running Λ(t) starts from zero in the IR, it increases/decreases and falls again to be zero in the UV. For the third choice the spacetime evolves from (anti)de Sitter in the IR to another (A)dS in the UV. Therefore, in the ultraviolet regime, we can have the following two sensible scenarios: free gravitons propagating on Minkowski or free gravitons propagating in (A)dS backgrounds.
It is worthwhile to make a comment about the generality of the results obtained here. As we mentioned in the previous section, the operator GB 1 and other cubic operators in the Riemann tensor and/or its contractions may affect the beta functions for the couplings θ C , θ R , θ GB and ω κ . In the present work, we did not compute these contribution of GB 1 . However, the contributions to β C , β R and β GB dependent on ω GB can affect the analysis in this section. Indeed, the ω GB -dependent contributions can flip the signs of all these beta functions and the corresponding conditions for asymptotic freedom in respective couplings θ C , θ R and θ GB will require to flip the signs of θ C (0), θ R (0) and θ GB (0) respectively as well.
The RG-improved version of the action is a reasonable approximation to solve the problem of the locations of poles of the propagator and their shifts due to quantum effects. It is definitely true that conditions for complex conjugate pairs of poles for LW particles here are rather restricting. However, since we expect that quantum corrections are small and perturbative, then we expect the same situation to be present on the full quantum level and we express the hope that it is probably possible to remain with the LW QG also on the one-loop level when quantum perturbative effects are taken fully into account.

VI. ASYMPTOTIC FREEDOM
In this section, we discuss the issue and necessary conditions for asymptotic freedom (AF) of the six-derivative QG models. First, we remark that usual textbook definitions of AF as typically used in two-derivative theories and for couplings like Yang-Mills charge g (that is requiring that β g < 0) could be a bit misleading. Here one has to exert special care when discussing AF in higher derivative theories. Actually, we want to point out that in the more general situation the condition β < 0 loses its validity. For example, even in the example of QCD, the beta function of the reparametrized coupling ω = g −2 is not negative in the deep UV regime. Actually near the trivial Gaussian UV fixed point (FP) (still describing AF pure Yang-Mills theory), it reaches a constant positive asymptotic value 11 3 C2(G) (4π) 2 . Based on this simple example, we need to allow for an extension of the definition of AF for the case of different coupling parametrizations and also for theories that in the kinetic terms for fluctuations contain higher than two-derivative terms. In a general setup, one has to be very careful in using the condition β < 0 for AF and assigning a potential physical meaning to it.
Actually, precisely such need we find in our six-derivative gravitational theory where we write our beta functions for ω-like (like in the QCD example discussed above) type of couplings (compare with formula (2)). Moreover, in higher derivative theories the couplings have not yet been measured and we can always consistently achieve AF by switching at our wish the signs of the classical coupling constants, which are actually the initial conditions for the RG flow equations. Therefore, we propose the following definition of physical asymptotic freedom (PAF). It is simply termed as the requirement that the interaction terms are suppressed compared to the kinetic terms describing free propagation of particles with quantum effects included. What would be a more physical definition of AF? The definition of PAF above is certainly more physical than the general theoretical requirement of β < 0 (holding only for some way of writing couplings and only in two-derivative theories). In more general theories with more general types of couplings, the system of beta functions and conjunctions of conditions and inequalities for AF are more intricate.
Our definition of PAF is one of the simplest and is very physical (and this implies that quantum scattering amplitudes, of course due to interactions, should go to zero in the deep UV regime). The reference for such a definition of PAF is the seminal paper by Fradkin and Tseytlin [36] and [37], where it is shown that our accepted definition is likely the physically reasonable one. We here follow exactly the same rescaling procedure as done by the authors above (see below where we will perform it explicitly and with all the details). The procedure is based on the rescaling of the graviton field by the same overall factor and showing that the physical interactions are suppressed compared to kinetic operators and that this suppression factor really tends to zero asymptotically at infinite energy.
In analogy with Stelle's quadratic gravitational theory [7], ω κ in our paper plays the role of the coupling f −2 in Stelle's theory and all the other couplings, which can run to a constant value in the UV, run logarithmically with the energy scale, or they do not run at all because they are in front of terms in the action containing more than four derivatives of the metric, playing the role of general ω couplings in Stelle's theory (see [36] about the definitions of f and ω there). We are here assuming that in the six-derivative theory all the running couplings are AF in the UV according to the common definition, namely there are no Landau poles for them at any finite energy scale.
Additional conditions for the RG flow towards an AF theory in the UV regime are that the initial values for the couplings, which we set conventionally at t = 0 are selected in order to avoid tachyons and real ghosts. This choice must make physical sense for the theory, that is the theory must be well defined and for example in gauge theories we cannot have a negative energy of perturbative excitations around vacuum (this enforces that the electric or YM charge is real, not imaginary). Moreover, we demand that the flow of all running couplings does not meet any singularity, namely there are no Landau poles for any finite value of the logarithmic energy scale t. The zero point values are also excluded from the flow, because there the coupling could change the sign and this could for example invalidate the positivity of energy as discussed above. These two conditions translate into inequalities for the initial values of the couplings and their first derivatives at t = 0. For example, assuming that ω C (t = 0) > 0 for AF we must also require that β C (t = 0) > 0. (For ω-type couplings it is natural to assume β ω > 0 for AF in UV, contrary to g-type couplings, like the electric or YM charge, for which we shall assume β e < 0 or β g < 0 for AF in the UV regime.) To clarify we point out that for the ω-type of couplings, the situation with ω = 0 corresponds to what is commonly called as Landau pole and ω → ∞ is a trivial Gaussian fixed point of the RG flow. In terms of g-type couplings, the situation is opposite and g = 0 is a free asymptotically noninteracting theory and g → ∞ is a sign of loss of perturbativity and of meeting a hypothetical Landau pole.
Assuming that all the couplings flow in the regular (described above) way and reach AF in the UV regime, we can eventually perform a rescaling analysis of the graviton field in order to compare the kinetic term with the interactions. For this propose we rescale by the coupling ω κ to have the two-derivative kinetic term of the graviton field in the canonical form.
Therefore one can easily prove that PAF is the feature of our model in the UV, provided the specific signs of beta functions that we have computed in earlier sections of this paper. By such an elegant analysis with field rescalings we do much more that just analysis of signs of beta functions. And we show below that PAF works in a more complicated setup of higher derivative kinetic (this part is usually not considered in standard textbooks on QFTs) and interaction terms.
In this way we convince the reader that PAF is a more physical criterion than theoretical β < 0 and that it indeed describes the physical features of quantum interactions compared to quantum propagation of free modes. Moreover, this is also the definition used in Stelle's theory and higher derivative gauge theory in [36] (see the appendix of that paper). Therefore we confidently can say that under specific conditions the theory is free in the UV since only propagation of free modes matters in the UV regime. In short, we find a free graviton particle in the UV.
In order to complete the whole story, here we explicitly show the asymptotic freedom of the theory (53), assuming the cosmological constant to be zero or constant in the ultraviolet regime realizing the second case discussed above. Since the Newton constant proportional to ω −1 κ in the UV falls off faster than all other running couplings (cf. (50), (51), (52), (81) and (82)), we rescale the graviton field as follows: The effective (RG-improved) action from (53), at the second or higher order in the perturbationh µν and around flat spacetime, reads where we explicitly compared the terms with the same number of derivatives and we neglect writing terms linear in the perturbationh µν due to onshell condition. The above higher order terms are symbolically indicated by O(. . .) and (0) means the d'Alembert operator constructed with the unperturbed metric. Moreover, with the labels (1) and (2) we pointed the expansion at the first and second order in the perturbationh µν . The dominant kinetic terms (quadratic inh µν ) around flat background were all indicated in (87), while interaction terms, possibly containing also derivatives, are of orderh 3 µν or higher. Due to the schematic form of the running in the UV, one sees that the ratios in O(. . .) brackets are suppressed for large t (e.g., Similarly 1/ ω κ (t) ∼ t −1 and we have also assumed, consistently with the previous discussion, that Λ(t) ∼ t 0reaches a constant value (or zero) in the UV regime. One could also instead of looking at asymptotic relations in (88) valid for t → +∞, consider the exact form of the RG flows for couplings as found in (79) and in (63), (64). This is why it is clear from the expansion (87) that in the ultraviolet regime, t ≫ 1, all interaction terms involving three or more gravitons are suppressed with respect to the kinetic terms.
One can also see the universal feature of the rescaling procedure that we have applied above. Indeed, the theory turns out to be AF regardless of which running coupling parameter is used to do the rescaling of the gravitational fluctuation field h µν . Above we used ω κ (t), provided it has a regular form of RG running and the corresponding inverse coupling ω −1 κ (t) does not meet any zero (ω −1 κ → 0) (giving AF theory), nor a Landau pole (ω −1 κ → ∞) at any finite t value. But for this purpose we could equally well use a running parameter θ C (t) or θ R (t), provided that the same conditions for the regular RG flows are satisfied. We still should assume that the cosmological constant parameter Λ does not exhibit any RG running. This universality of the rescaling analysis is another confirmation that our results about PAF in the UV regime in our theory are very robust.
So far so good, but unfortunately both ghosts in (72) and in (73) have positive real parts and, therefore, are tachyons in our signature of the metric [74]. One of the ways to overcome this obstacle consists in adding at least two operators cubic in the curvature to provide the vanishing beta functions for θ C and θ R . Alternatively, the beta function β C can be made zero also by selecting x in (40). However, this choice is inconsistent with a zero beta function for the cosmological constant. Let us check this statement. The new modified beta functions will have the following form: Here c i = (c C , c R , c GB , c κ ) are contributions coming from the cubic in curvature nonminimal terms. All such c i can be, in principle, made arbitrary by choosing the corresponding coefficients in the classical action. It is clear that proper choices of c C , c R can make zero the beta functions β C and β R . We analyze such a situation here. Therefore, we can consistently take the initial conditions θ C (0) and θ R (0) to be zero because, in this case, the theory is finite for what concerns the operators C 2 and R 2 in the divergent action (37). This also means that these two terms can be removed completely from the classical action (53) and consistently they will not be recreated by quantum corrections (θ C (t) = θ R (t) = 0). Moreover, we can replace in (93) to finally get β Λ = 0. Hence, we can fix Λ = 0 consistently at any energy scale, at the one-loop level. Under these conditions, at the end, only ω κ (t) > 0 can run and there we have asymptotic freedom as elucidated above. By adjusting to zero the value of the c κ coefficient, one can make absent also the running of ω κ and the theory is completely UV-finite. Finally, the relations (77) and (78) turn into From (94), it is obvious that we must consider only two cases: either ω C , ω R > 0 or ω C , ω R < 0. In the first case, the spectrum consists of the graviton, two LW-like spin-two fields (from ∆ 2 < 0) with purely imaginary mass square parameter, one real scalar and a scalar tachyon (from ∆ 0 > 0). Thus, the theory is unitary according to Anselmi-Piva's prescription, but it is unstable because of the presence of the tachyon. On the other hand, if we take ω C , ω R < 0 the situation is reversed. We end up with the poles of the ghosts in the spin-two sector becoming real: one ghost field and a ghost tachyon. In the spin-zero sector we meet two LW-like fields with purely imaginary mass square parameter. One concludes that the conditions (95) and (96), together with (94) are inconsistent, mainly due to the presence of ghosts and tachyons. When one allows for the running of the Newton constant ω κ → ω κ (t), the situation does not get any better. Therefore, in order to avoid tachyons, one may try to make the cosmological constant nonzero and make it to run with the energy scale as well as the Newton constant, according to the following formulas: We still assume that θ C (t) = θ R (t) = 0 both in the classical action and in its RG-improved quantum version. The ghost poles are now located in the following points (in accordance with (72) and (73)): spin-2 : z In order to avoid tachyons in the spin-0 sector case we must take Λ < 0 to have negative real contribution to the mass squared, while the spin-2 sector is without this problem. The two fields of spin-2 have both zero real parts of the mass square parameter, which is then a purely imaginary quantity. We assume that ω C > 0, ω R < 0 (or ω R > 1/15 ω C ) and also that ω Λ (0) > 0. Moreover, ω κ (0) > 0 and c κ > 0, so ω κ (t) > 0. Note that now the running parameter ω Λ (t) > 0 and hence the cosmological constant Λ(t) is negative according to (76), while the constraints (77) and (78) simplify to Both constraints cannot be satisfied since both ∆ 2 < 0 in (101) and ∆ 0 < 0 in (102) effectively lead in the large t leading asymptotics to the expression which obviously cannot be negative for any ω C and ω R . This is a simple consequence of the fact that in (101) and in (102) the Λ(t) 2 term has the leading t 2 asymptotics higher than that of ω κ (t) ∼ t (cf. (98)). One can also consider the situation where only one of the couplings ω Λ or ω κ shows quantum renormalization running. In the first case, when c κ = 0, the running ω Λ (t), which scales like t according to (97), again dominates the relations (101) and (102) in the UV regime of large t due to the fact that Λ(t) ∼ t. Now repeating verbatim the same argumentation as above, one sees that it is impossible to satisfy both ∆ 2 < 0 in (101) and ∆ 0 < 0 in (102) for all values of t. When ω C = 15ω R and c κ = 0, the coupling ω κ is the only one which runs and it is linear in t. In this situation, ω κ (t) is the dominant quantity for large t in relations (101) and (102), so effectively we can forget about Λ(t) there (which falls off here like 1/t). This means that in the large t asymptotics we are back to relations (95) and (96) obtained for Λ(t) = 0 that we showed above cannot be satisfied simultaneously together with ω C = 15ω R .
In the final comment, we can return to the case when the full theory is completely UV-finite, but when we keep a nonzero value of the ω Λ coupling. The choice of signs explained above implies that Λ < 0. A careful analysis of the inequalities (101) and (102) with all the couplings taking constant values and satisfying the relation ω C = 15ω R shows that again it is impossible to find a nonzero allowed interval for values of the ω R coupling. Simply, restricting to ω R > 0, the ∆ 2 < 0 implies ω R < 3 160 ω κ Λ −2 , while on the other side ∆ 0 < 0 implies that ω R > 3 8 ω κ Λ −2 . One therefore concludes that some of the couplings θ C or θ R must run with energy too to provide enough room for keeping the pair of Lee-Wick particles with complex conjugate poles at any energy scale.
Eventually, the full theory consistent with the analysis in the last paragraph will include all or some of the following cubic curvature invariants, written in general dimensions but due to various identities holding in four spacetime dimensions [4] only three out of the five operators containing a Riemann tensor in (105) are independent there. Therefore we can add the following six nonminimal terms below, which are cubic in gravitational curvatures with arbitrary coefficients s 1 , . . . , s 6 . The explicit form of the addition to the Lagrangian in (2) reads 8 With these additional operators at our disposal, we will likely get the modified beta functions in the general form (89)-(92). One can notice the following schematic relations of how the running couplings depend on the constant parameters s i with i = 1, . . . , 6 (where we keep only the highest powers). First, the couplings θ C (t), θ R (t) and θ GB (t) show up to quadratic dependence on s i (so they are proportional to s i s j in the numerators). The denominators of such ratios can be of three different types: ω 2 C , ω 2 R or ω C ω R . The running Newton coupling ω κ (t) will generally depend on a cubic form in s i coefficients in the numerators, while the running cosmological constant ω Λ (t) on up to a quartic one.
One last comment concerns our taken for granted trust in the energy dependence of the spectrum. Consistently with the perturbative expansion, the classical spectrum cannot be largely modified in the loop expansion. However, asymptotic freedom allows us to infer about the ultraviolet regime of the theory (mathematically this is the infinite energy regime). Therefore, we are forced to consider the spectrum at any energy scale with particular attention to the UV regime.

VII. CONCLUSIONS AND DISCUSSIONS
For the first time, we derived the exact beta functions for the four-derivative coefficients in the minimal, superrenormalizable model of the six-derivative QG. This minimal theory is the simplest super-renormalizable model, with only seven free parameters in the classical action (2).
The result is applicable to the minimal Lee-Wick quantum gravity proposed in [15,16], which has been shown to be unitary at any perturbative order and also admits the procedure suggested in [17][18][19]. For this model, we solved the renormalization group equations, showing asymptotic freedom and two possible perturbative scenarios in the ultraviolet regime: in one case the cosmological constant runs towards zero, while in the other case it has a negativevalue fixed point. In the former scenario a free graviton propagates in Minkowski spacetime, while in the latter one the graviton propagates in the (anti-) de Sitter spacetime, which may be unstable in the minimal six-derivative model because of the presence of tachyons in the spectrum. Therefore, we proposed an extension of the theory that very likely will be free of tachyons, while keeping asymptotic freedom on the (A)dS space in the ultraviolet regime.
Our novel results for UV divergences (proportional to C 2 , R 2 and the E 4 term) obtained in this work are completely free of any problem of gauge dependence. We emphasize that these results are gauge-independent and also gaugefixing-independent. This is another remarkable and very positive and physically welcomed feature of our model. And it is in opposition for example to the case of the beta function of the G N coupling in Stelle's quadratic theory in d = 4 spacetime dimensions, which still shows various ambiguities [53]. The UV divergences that we have found are fully unambiguous, gauge-independent and they represent true physical quantities since they have the meaning not only in the perturbative approach as one-loop expressions, but also they are perturbatively exact (as also explained below). The mathematical fact that these UV divergences are gauge-independent is proven by the powerful theorem (originally due to Kallosh, Tarasov, and Tyutin in [26]) that we have already explained in the point iii) of Sec. II. Hence there is not any problem with gauge dependence of our results.
Similarly, we just remark for completeness that the locations of the poles in both spin-2 and spin-0 parts of the gravitational propagator of the theory both around flat and (A)dS backgrounds are completely independent of the gauge choices used to compute the propagator. Hence also the masses of spin-0 and spin-2 particles present in the spectrum are completely gauge-fixing independent. This fact was crucial for the consistency and the physical sense of the analysis that we have presented in Sec. V and in Sec. VI. This is in accordance with the existing literature on the topic.
We also put special attention to the fact that the condition for AF depends crucially on the initial conditions of the RG flow and also on the signs of the beta functions. This point we explained at length also before. In many places above in Sec. VI, we gave reasonable physical inequalities for the signs of the initial values of the couplings. For example, we required that ω R (0) < 0, ω C (0) > 0, ω GB (0) > 0, etc. We also demanded for consistency and AF in the UV the corresponding matching signs of the beta functions such that the RG flow for all finite values of t is regular and does not touch neither zero nor infinite values.
We would like to remind the reader that exactness of our beta functions has to do with the absence of perturbative loop divergences at loop levels higher than the first one. Simply the theory is super-renormalizable. This conclusion is based on power counting analysis of UV divergences.
In this gravitational model with six derivatives, one is sure that besides the one-loop level there are no new perturbative UV divergences proportional to the C 2 , R 2 and the E 4 terms. Since one reads the beta functions of the couplings of the theory from these divergences, then one concludes that expressions for the beta functions at the one-loop level take into account all loop contributions and hence corresponds to full loop resummation. Assuming that there are no nonperturbative contributions to the beta functions of couplings, then this means that the beta functions that we have obtained are exact and there will be no additional contribution to them due to any other quantum effects. Just one-loop quantum effects encoded in UV divergences are completely enough to settle the exact form of the RG running of couplings θ C , θ R and θ GB in this model. Moreover, the structure of the beta functions for ω κ and for ω Λ as written in (46) and (47) respectively must be repeated also when two-, two-and three-loop contributions are taken additionally into account there. This is due to an expected form of divergences at two and three loop orders, which compared to formulas in (46) and (47) from dimensional reasons can be only multiplied by the powers (both positive and negative) of the fundamental ratio x of the theory. Hence the integration of these formulas leads to the running the same as encoded in formulas (81) and (82) respectively, where only the constant numerical coefficients may be different after inclusion of higher loop effects. But the leading t 2 and t 3 behaviors respectively for the running coupling parameters ω κ (t) and for ω Λ (t) in the UV regime are a universal feature here. Together with exactness of one-loop expressions for running θ C (t), θ R (t) and θ GB (t), this is one of the most beautiful and powerful features of the super-renormalizable model we have considered in this work. Therefore this model constitutes a good and promising theoretical laboratory for studying RG flows.
Let us mention two possible extensions of the results described in this paper. The beta functions (43), (44), (45), as well as the beta functions for the linear in curvature term and the complete one-loop beta function for the cosmological constant, can be used for testing alternative approaches to QG. For instance, there are interesting approaches, that are regarded nonperturbative, such as functional renormalization group or dynamical triangulations. It would be interesting to see whether these methods can reproduce the beta functions that were derived here. Such a comparison represents one of the possibilities to extend the present work.
Moreover, this paper lays the foundations for future computations with more involved (one can say nonminimal) models for a super-renormalizable QG. As we have discussed above, there are many possible extensions of the minimal Lagrangian (2). In this respect, the situation is very similar to the minimal Standard Model of particle physics. However, the important difference is that for the QG models there are no experimental or phenomenological restrictions on the higher derivative parameters of the simplest model (2), or on its nonminimal extensions. On the other hand, it would be interesting to explore the IR limit of the theory with complex poles without or with loop corrections, e.g., in the way it was done recently in [3]. In principle, a detailed analysis of this aspect of the models under discussion could be useful for many phenomenological applications.

APPENDIX A: WICK ROTATION IN A GENERAL LOCAL HIGHER DERIVATIVE THEORY
Here we address the issue of Wick rotation in higher derivative theory. The RG equations are by definition obtained in the deep UV regime of the Euclidean theory. And the question how to rotate them back to Minkowskian/Lorentzian signature has been addressed in Appendix A of Ref. [81]. Let us here expand the latter proof.
First of all let us remind the reader of the main feature of dimensional regularization. We can consider the following prototype of a one-loop integral suitable for a general higher derivative theory, where s and r are integers. We can consider the case in which r = 0 and s = −N and define the following integral in dimensional regularization: In order to compute the above integral we replace r = 0 and s = −N on the right side of (107), namely For N + D/2 > 0 (strictly) or N > −D/2 the Gamma functions in (109) are singular, but the limit in (108) is well defined and it gives zero for N > −D/2, i.e.
Other more rigorous proofs can be found in [82,83]. In short, following [83], the integrals in (109) can be defined as sums of other integrals convergent in a proper domain of the complex variable D (the dimension of spacetime). The result of these integrals turns out to be analytic in D and, therefore, its uniqueness allows us to compare integrals initially defined in different domains. The values of such integrals can finally be combined to reconstruct the integrals (109), which turn out to be zero. We now apply the result (110) to the integrals on the arcs in the first and third quarter of the energy complex plane and we show that they vanish in dimensional regularization. Once again, this is a general result true in any regularization scheme, but easy to prove in dimensional regularization.
Following the discussion in [81], the general one-loop integral (107) can be written in the following form: where here C is a scalar function of the external energy p 0 and momenta p combined in a Lorentz-covariant vector p = (p 0 , p) and possibly of masses, a j (k) is a polynomial of k and b i (p, k) is polynomial of k and the external Ddimensional momenta p, while s 1 and r 0 are integers. (Similarly, we have the D-dimensional Lorentz vector k = (E, k) over which we integrate.) The energy integral in (111) is divergent for 2r + 1 2s, hence, we define 2r = 2s + n − 1 and take n 0. In order to compute the contribution to the integral (111) along the arcs, let us take |E| much larger than all the other scales in the integrand. Therefore, the divergent contributions are obtained expanding the denominator in (111) as follows: Notice that the integral in the energy E is actually a convergent angular integral multiplied by the radius |E|. Therefore, we have to focus only on integrals that diverge in the limit |E| → ∞ after the angular integration has been performed. Since the theory is local, the number of such divergent integrals is finite and, therefore, we can add a finite number of integrals in (D − 1) dimensions that are "zero" as a particular feature of the dimensional regularization in (D − 1) spacetime dimensions. Indeed, all such integrals are polynomial integrals of k that turn out to be zero when the dimensional regularization (110) is implemented in the (D − 1)-dimensional manifold. Therefore, whenever we need to close the integration path, in the lower or upper half plane, we actually add a finite sum of vanishing contributions.
Let us expand about this statement. In general, we add to the integral (111) a finite number of "zeros" by means of using the dimensional regularization (DIMREG) scheme. However, for the sake of simplicity, we here consider only the case r = s or n = 1 (for Einstein's gravity r = s = 1), hence, in the expansion (112) only the first volume contribution is divergent, Now, we observe that when we take the limit |E| → +∞ (large radius in the first and third sector in the energy complex plane) only the first integral in the sum (113) diverges. However, it is identically zero or, which is equivalent, it cancels with the introduced "zero", as a particular feature of dimensional regularization. Notice that in the last step we did not use any expansion, but only the limit for large |E| in the first of the two integrals. Of course, the sum of two integrals is the sum of the integrands. In the general case 2r + 1 2s and introducing the definition 2r = 2s + n − 1, the integral reads: where c 1 , . . . , c n are constants.
Coming back to our calculation in the main text of the paper, it was done exclusively and from the beginning in the Euclidean domain (because of using Euclidean-signature based Barvinsky-Vilkovisky trace technology).
Afterwards, in order to get the theory in Minkowski spacetime, we need the Anselmi-Piva (AP) prescription [17][18][19] for which it is technically crucial that the contributions on the arcs discussed in this appendix vanish in DIMREG in (D − 1) dimensions. However, such prescription does not affect the divergences that are computed in Euclidean signature. Using the AP prescription the divergences are the same in Euclidean as in the Lorentzian case without any need to do explicitly Wick rotation. One could even say that the UV divergences that we have obtained as considered in Lorentzian theory are independent of any prescription scheme of how to encircle/include/exclude poles of the propagator in the Euclidean theory and this issue of how to close the integration contour on the complex plane is not essential for performing formal Wick rotation to the Minkowskian domain. This is so since these are subleading effects (they do matter for finite terms, but not for the leading in the UV terms of the expansion of loop integrals). Hence the UV divergences in the form we have found are the same in Euclidean as well as in Lorentzian signature. Here we assume that there exists a well defined calculational procedure in the Euclidean regime, or that the same calculations can be equivalently performed using Feynman diagrams, but for the UV divergences the results will be the same. The situation is a verbatim repetition of the case of beta functions of higher derivative Stelle's quadratic theory [7] in D = 4 spacetime dimensions. Also there the beta functions are valid universally both in Euclidean and Minkowskian theory. Of course, in the latter case, one has to understand the various curvature tensors (Riemann, Weyl, etc.) as computed using the spacetime metric with Lorentzian signature as opposed to the former case of