On avoiding Ostrogradski instabilities within Asymptotic Safety

We study the renormalization group flow of gravity coupled to scalar matter using functional renormalization group techniques. The novel feature is the inclusion of higher-derivative terms in the scalar propagator. Such terms give rise to Ostrogradski ghosts which signal an instability of the system and are therefore dangerous for the consistency of the theory. Since it is expected that such terms are generated dynamically by the renormalization group flow they provide a potential threat when constructing a theory of quantum gravity based on Asymptotic Safety. Our work then establishes the following picture: upon incorporating higher-derivative terms in the scalar propagator the flow of the gravity-matter system possesses a fixed point structure suitable for Asymptotic Safety. This structure includes an interacting renormalization group fixed point where the Ostrogradski ghosts acquire an infinite mass and decouple from the system. Tracing the flow towards the infrared it is found that there is a subset of complete renormalization group trajectories which lead to stable renormalized propagators. This subset is in one-to-one correspondence to the complete renormalization group trajectories obtained in computations which do not track of the higher-derivative terms. Thus our asymptotically safe gravity-matter systems are not haunted by Ostrogradski ghosts.


Introduction
Constructing phenomenologically interesting quantum field theories which are valid at all length scales is one of the central topics in theoretical high-energy physics to date. For gravity, such theories may be realized through the Asymptotic Safety mechanism, see [1][2][3][4][5][6][7][8] for reviews. As first suggested by Weinberg [9,10] this mechanism could provide a consistent and predictive quantum theory of gravity within the well-established framework of quantum field theory. The key idea of this program is that the gravitational renormalization group (RG) flow possesses a non-trivial RG fixed point. At this fixed point (some of) the dimensionless couplings take non-zero values, so that the resulting theory is interacting. This is in contrast to the free (or Gaussian) fixed points underlying asymptotic freedom where the high-energy completion is provided by a free theory. Starting from the seminal work [11], there is, by now, substantial evidence that gravity in four spacetime dimensions actually possesses a non-Gaussian fixed point suitable for Asymptotic Safety. In particular, it has been shown that this fixed point is robust under the inclusion of the two-loop counterterm [12] and is connected to a classical regime through a crossover [13]. Besides ensuring the absence of unphysical divergences, this fixed point also comes with predictive power: any UV-repulsive direction of the fixed point allows to express the associated coupling as a function of the relevant parameters. The resulting relations may be tested experimentally, at least in principle.
While the prospects of obtaining a quantum description of the gravitational force valid at all length scales is already intriguing, it is also clear that a realistic description of our world also requires the inclusion of matter degrees of freedom. While there has already been significant effort geared towards understanding the role of the Asymptotic Safety mechanism for gravity-matter systems, the picture is still far from complete. In order to discuss potential UV-completions of gravity-matter systems it is useful to distinguish between the two cases where the matter sector of the underlying fixed point is Gaussian or non-Gaussian in the sense that matter self-interactions are either absent or turned on. On general grounds, one may expect though that non-trivial interactions in the gravitational sector also induce non-trivial matter self-couplings, see e.g. [14] for a discussion. Depending on the details of the approximation used to investigate the fixed point structure of the gravity-matter system, it is conceivable that a matter fixed point which is actually non-Gaussian may be projected onto a Gaussian one if the approximation used to probe it does not include self-interactions. Conversely, a fixed point identified as Gaussian may split into a Gaussian and non-Gaussian one once additional couplings are probed.
In order to get an idea which matter sectors could actually be compatible with Asymptotic Safety, Refs. [15][16][17][18] studied projections of the full RG flow where the matter sector contained an arbitrary number of minimally coupled scalars N s , vectors N v , and Dirac fermions N D . Complementary results for the case where spacetime carries a foliation structure have been reported in [19]. While all studies agree on the statement that the matter content of the standard model of particle physics leads to a fixed point structure suitable for realizing the Asymptotic Safety mechanism, the precise values for N s , N v , and N D supporting a NGFP are different. Restricting to the cases where the matter sector contains scalar fields only, [15][16][17][18] report an upper bound N s 16 − 20, while in [19] no such bound is present in agreement with the initial works [20,21]. This difference can be traced back to different choices for the coarse graining operators and definitions of Newton's constant employed in these works. In particular, Refs. [16][17][18] define Newton's constant based on the flat space graviton propagator while [15,19] resort to a background Newton's constant. As argued in [17] matter degrees of freedom contribute differently in these settings. The two pictures are in qualitative agreement if N s is small but start to deviate once the matter contribution becomes significant.
In a complementary approach, the fixed point structure arising within scalartensor theory has been studied in [22][23][24][25][26][27][28][29]. 1 This setup includes two arbitrary functions of the scalar field φ, a scale-dependent scalar potential V k (φ) and a function F k (φ) encoding the coupling of the scalar field to the Ricci scalar. In d = 3 this setting gives rise to a Wilson-Fisher type RG fixed point which can be understood as a gravitational-dressed version of the Wilson-Fisher fixed point known in a nondynamical flat background. In d = 4 the analogous analysis indentifies a fixed point with a Gaussian matter sector. In particular the scalar mass and φ 4 -coupling vanish at this fixed point. Ref. [29] supplements this setting by a third scale-dependent function K k (φ) dressing the scalar kinetic term. In this generalization also a non-Gaussian matter fixed point has been identified.
The influence of gravity on the flow of gauge-couplings has extensively been discussed in both perturbative [34][35][36][37] and non-perturbative [38][39][40][41][42] settings. Fundamental aspects related to the inclusion of fermions have been discussed in [43,44] and the compatibility of light chiral fermions with asymptotic safety has been argued in [45][46][47]. Starting from the prediction of the Higgs mass based on Asymptotic Safety [48], mass hierarchies in the standard model and its extensions have been studied in [49][50][51] while the influence of gravitational interactions on the flow of Yukawacouplings has been studied in [52][53][54][55][56][57]. 2 Based on these works there have been several key insights related to asymptotically safe gravity-matter systems. Firstly, non-Gaussian fixed points in the matter sector may come with a higher predictive power than their Gaussian counterparts. In Ref. [39] this property has been used to predict the value of the fine-structure constants based on the Asymptotic Safety mechanism. Secondly, a non-vanishing fixed point value for the U (1) hypercharge may provide a solution to the triviality problem of the standard model [41]. Thirdly, the Higgs mass can be predicted correctly based on the beta functions of the standard model completed by the Asymptotic Safety mechanism above the Planck scale [48].
These salient features are, however, also accompanied by the lurking danger that the non-vanishing gravitational interactions may induce potentially dangerous terms in the fixed point action. Typical candidates are higher-derivative terms contributing to propagators of matter fields, which are typically associated with Ostrogradski instabilities or the violation of unitarity, see [65,66] for reviews. In this work we initiate the study of this class of interaction terms for gravity-matter flows. For transparency we focus on the simplest possible model comprising the Einstein-Hilbert action supplemented by minimally coupled scalar fields including a higher-derivative term in the scalar propagator. We show that, as expected, the higher-derivative term is generated along the RG flow. Quite remarkably, the flow admits RG trajectories for which the ghost degrees of freedom decouple in the renormalized propagator. These findings constitute a highly non-trivial consistency test concerning the structure of asymptotically safe gravity-matter systems. From a complementary viewpoint they also provide the initial step towards extending the classical stability analysis of Horndeski [62] and "beyond Horndeski" theories [63,64] to the quantum level.
The remaining work is organized as follows. The Ostrogradski construction and its loop-holes are reviewed in Sect. 2. Sect. 3 introduces the setup of our RG computation incorporating a higher-derivative kinetic term in the scalar sector, and reports the resulting beta functions in Sect. 3.2. The properties of the RG flow are investigated in Sect. 4 and we discuss the consequences of our findings in Sect. 5. Technical details related to the evaluation of the flow equation using a non-smooth regulator are relegated to Appendix A.

Higher-derivative terms and Ostrogradski instability
We start by briefly reviewing the classical Ostrogradski instability and its loopholes, mainly following the expositions [65,66].

The instability . . .
It was shown by Ostrogradski in the 1850's that non-degenerate classical systems containing time derivatives of finite degree larger than two give rise to Hamiltonians whose kinetic term is not bounded from below [67]. Irrespective of the exact form of the action, the unbounded Hamiltonian will yield several unwanted phenomena, related to the instability of the system. At the classical level, the presence of degrees of freedom coming with a wrong sign kinetic term allows to accelerate particles to infinite velocity while keeping the total energy of the system constant.
This type of instability also appears in the corresponding quantum system. While the presence of higher-derivative terms in the propagators lowers the degrees of divergencies arising in loop computations, the presence of positive and negative energy states may trigger an instantaneous decay of the vacuum. Naively, a way out may be to reinterpret the negative-energy creation and annihilation operators as positiveenergy annihilation and creation operators, respectively. Although this seems to cure the instability of the vacuum state, this procedure yields states with negative norm. Removing these states from the physical spectrum, however, yields a non-unitary S-matrix.
In the case of a non-interacting scalar field theory, the Ostrogradski instability can be nicely illustrated by the Källén-Lehmann representation [68]. This representation expresses the dressed propagator G(x − y) as a superposition of freely propagating particles with mass µ ≥ 0 and propagator For a unitary theory, the spectral density ρ(µ 2 ) is a sum over norm-states with positive coefficients, thus ρ(µ 2 ) ≥ 0. If ρ(µ 2 ) < 0 for some µ 2 in the physical sector of the theory, then unitarity issues arise.
In this article, we will study a system containing scalar fields φ where the propagator contains a fourth order kinetic term (see Sect. 3.1) 3 where Z denotes a wave-function renormalization and Y is the coupling associated with the higher-derivative term. This has a propagator expanded in a Fourier basis given by Using partial fraction decomposition, we can expand this in terms of free propagators: (2.5) We see that the Källén-Lehmann spectrum contains a massless state with positive density, and a state of mass with negative density. The latter state is called a (Ostrogradski) ghost. It is easy to see that the spectral density is not positive. Therefore the theory will generically be unstable.

. . . and its loop-holes
Although higher derivatives generically introduce severe fundamental flaws in a theory, there a number of ways to bypass this problem. This can be done at both the classical and the quantum level. One way for curing the Ostrogradski instability at the classical level is to lift the condition of non-degeneracy. In this case the higher-order time derivatives are removed by either combining them into total derivatives or using a gauge symmetry. In the former case, the total derivatives in the Lagrangian do not contribute to the dynamics. Provided that this procedure removes all higher-derivative terms, this results in a healthy theory. 4 In the latter case, gauge symmetry can be used to impose an extra condition to the equations of motion. If these constraints remove the higher derivatives, the instability is cured as well.
A second option consists of replacing the terms appearing in the straight bracket of eq. (2.3) by an entire function of the momentum possessing a single pole of first order. This strategy results in a non-local theory which contains time-derivatives of infinite order. In this case the propagator does not admit a partial fraction decomposition and the absence of poles in the physical spectrum implies that the theory is still stable. However, the question if the resulting non-local theory is well-posed is subtle. An exposition on the treatment of this class of theories is given in [71,72]. 5 When assessing the stability of a higher-derivative theory at the quantum level, the situation becomes even more involved. In this case the dressed propagator of the theory can be obtained from the effective action Γ and one expects that for a stable theory this propagator does not give rise to Ostrogradski ghosts. Following the discussion of the classical case above, this may be realized in two ways: a) pushing the mass of the Ostrogradski ghost to infinity. b) completing the dressed propagator into an entire function.
The first case can be illustrated by considering the action (2.3). At the quantum level the coupling Y will depend on the renormalization group scale k, which we indicate by Y k . The requirement that the higher order derivative term does not contribute to the dressed propagator corresponds to demanding that lim k→0 Y k → 0. At the level of the decomposition (2.5), sending Y → 0 means that the ghost mass goes to infinity. The ghost then decouples from the spectrum of the theory and does not entail an instability. 6 This scenario may be realized in two ways. Firstly, the system may exhibit a fixed point located at Y * = 0. The theory at the fixed point is scale invariant and ghost-free. Secondly, an RG trajectory may be attracted to the Y k = 0 hyperplane as k → 0. The ghost will drop out of the effective propagator rendering the renormalized theory effectively ghost-free.
When investigating case a), gravity plays an essential role. In its absence, the action (2.3) describes a one-parameter family of non-interacting theories parameterized by Y . The only ghost-free theory in this set is Y = 0. This picture changes once 4 The point that a healthy theory has to remove the entire tower of higher-derivative terms has been stressed in [70]. We are greatful to H. Motohashi for bringing this work to our attention. 5 For a more detailed discussion of infinite-order theories in the context of gravity we refer to [73][74][75][76]. 6 For a similar discussion in the context of higher-derivative gravity see [81].
a minimal coupling to the gravitational field is included. In this case the gravitational interactions induce a non-trivial flow of Y k , opening the door to the nontrivial scenarios described above. At this stage the following remarks are in order. Firstly, we stress that the condition that the theory should be ghost-free applies to the dressed propagator (obtained at k = 0) only. At finite values of k it is expected that the process of integrating out quantum fluctuations mode-by-mode will generate higher-order derivative terms in the intermediate description. This does not signal the sickness of the theory, as its degrees of freedom should be read off from the dressed propagator. Secondly, investigating the case b) will require generalizing the simple ansatz (2.3) to a scale-dependent function of the momentum. In [82] it has been shown that this class of models suffices to obtain the Polyakov effective action from a renormalization group computation. This generalization is beyond the present work though, so we will not discuss this case in detail.

RG flows including higher-derivative propagators
Following up on the general discussion of Sect. 2, we now perform a RG computation determining the scale-dependence of the higher-derivative coupling Y in a gravitymatter setting. The key results of this section are the beta functions (3.14), (3.15), (3.19) and (3.21) which govern the RG flow of our projection.

The functional renormalization group equation and its projection
Currently, the predominant tool for investigating the fixed point structure and RG flows of gravity and gravity-matter systems is the functional renormalization group equation (FRGE) for the effective average action Γ k [11,[83][84][85] Here t ≡ ln(k/k 0 ) denotes the logarithmic RG scale, Γ (2) k is the second variation of Γ k with respect to the fluctuation fields and Str contains an integral over loop momenta and a sum over component fields. The regulator R k provides a mass-term for fluctuation modes with momenta p 2 k 2 and vanishes for p 2 k 2 . The interplay of the R k -terms in the numerator and denominator then ensures that the RG flow of Γ k is actually driven by quantum fluctuations with momentum scale p 2 ≈ k 2 . In this way the FRGE realizes Wilson's picture of renormalization where the RG flow is generated by integrating out fluctuations shell-by-shell in momentum space.
The FRGE comes with some highly desirable properties. Firstly, it allows the computation of RG flows without specifying a fundamental action a priori. This feature makes the equation tailor-made for identifying interacting renormalization group fixed points. Moreover, the regulator R k vanishes for k = 0 so that all quantum fluctuations are integrated out as k → 0. As a consequence the effective average action agrees with the standard effective action in this limit, lim k→0 Γ k ≡ Γ. Finally, the framework turns out to be sufficiently flexible to probe settings where different classes of metric fluctuations are admitted by either implementing a linear split [11], an exponential split [86,87], or an ADM split [88][89][90] of the gravitational degrees of freedom. Throughout this work, we will implement a linear split, decomposing the physical metric g µν into a fixed background metricḡ µν and fluctuations h µν according to Covariant objects carrying a bar are then constructed from the background metric while unbarred ones are constructed from g µν . Furthermore, we will set d = 4 throughout. While the generalization to general dimension d is straightforward the rather lengthy nature of the beta functions in the general case obscures the relevant structures, so that we make this choice for clarity. A common technique for finding non-perturbative approximate solutions of the FRGE consists of making an ansatz for Γ k , including the operators of interest, and subsequently projecting the full flow onto the subspace spanned by the ansatz. The beta functions governing the scale-dependence of the couplings contained in the ansatz are then read off from the coefficients multiplying the interaction terms contained in the ansatz. In order to study the effects of higher-derivative terms appearing in the scalar propagators of gravity-matter systems, we make the following ansatz for the effective average action The gravitational part of this ansatz is taken of Einstein-Hilbert form It includes a scale-dependent Newton's constant G k and cosmological constant Λ k . The gravitational sector is supplemented by a gauge-fixing action Γ gf k and a ghost term S ghost [g,c, c;ḡ]. In order to facilitate the comparison with the results reported in [15], we implement the harmonic gauge This gauge-fixing is accompanied by a standard ghost-term 7 The gravitational part of Γ k is supplemented by N s scalar fields, where ∆ ≡ −g µν D µ D ν is the Laplacian constructed from the full metric. Besides a wave-function renormalization Z k , this ansatz contains a scale-dependent coupling Y k associated with a higher-derivative contribution to the scalar propagator.

Evaluating the flow equation
Starting from the ansatz (3.3), the goal is to find the beta functions determining the scale-dependence of G k , Λ k and Y k as well as the scalar anomalous dimension η s = −∂ t ln Z k . This information is obtained by substituting the ansatz into the FRGE and extracting the relevant interaction terms from the trace appearing on the right-hand-side. The explicit evaluation of this operator trace requires specifying the regulator function R k . Throughout this work, we will resort to a Litim-type profile function [92,93], The matrix-valued wave function renormalization Z k is obtained from the substitution rule → P k ≡ + k 2 r( /k 2 ). Following the nomenclature introduced in [2], the coarse graining operator is chosen either as Type I : = ∆ , where the endomorphism E ≡ qR is chosen such that all curvature terms appearing in Γ (2) k become part of the coarse-graining operator. Using the Litim-profile in the regulating procedure has the advantage that all operator traces relevant in this work can be performed analytically. This comes at the price that the regulator is not smooth and the extraction of external momenta from the traces is non-trivial. In particular, contributions arising at the boundary of the momentum integrals have to be taken into account carefully. Our strategy for incorporating such terms is explained in detail in Appendix A.
The projection of the operator trace is then done as follows. The flow of G k and Λ k can be read off from the terms proportional to´d 4 x √ḡR and´d 4 x √ḡ , respectively. These contributions are conveniently found by selectingḡ µν as the metric on a 4-sphere and taking the background value of the scalar fieldφ = 0. The resulting operator traces can then be evaluated using standard heat-kernel techniques [2,7,11]. In this way, one arrives at the beta functions for the gravitational couplings given in eq. (3.14). The flow in the scalar sector is efficiently computed on an Euclidean background geometryḡ µν = δ µν and by expanding the background scalar fieldφ(x) in terms of Fourier modes. Setting the fluctuation fields to zero, the scalar sector appearing on the left-hand side of the flow equation is Thus the scale-dependence of Z k and Y k is encoded in terms coming with two powers of the background scalar field and two and four powers of the momentum q, respectively. The Feynman diagrams generating these structures are depicted in Fig. 1. They consist of a pure graviton tadpole, and two diagrams with scalar-graviton loop formed by connecting two three-point vertices. The projection of the flow equation then requires extracting the contributions proportional to q 2 and q 4 from these diagrams. Following the procedure described in Appendix A, this results in eqs. (3.19) and (3.21). The result of these computations is conveniently expressed in terms of the dimensionless couplings 11) and the anomalous dimension of Newton's constant and of the scalar field The scale-dependence of the dimensionless couplings (3.11) is encoded in the beta functions which we define according to For the dimensionless variables, the system of differential equations is autonomous in the sense that the beta functions are independent of k.
The explicit expressions for the beta functions in the gravitational sector are . (3.14) The anomalous dimension of Newton's constant is y and N s dependent. Inspired by [11], it can be cast into the following form: The functions B 1 and B 2 encode the contribution of the gravitational sector. For a Type I regulator, these functions have been determined in the seminal paper [11]. For a Litim-type regulator, they read For the Type II regulator, cf. eq. (3.9), these functions become Besides the gravitational self-interaction, there is a contribution of the scalar sector to the running of λ and g. For the latter, the additional scalar part is captured by In absence of higher derivative terms in the action, i.e. y = 0 and β y = 0 and setting the relevant ghost contributions to zero, this result agrees with [15]. Note that the choice of regulator, eq. (3.9), enters into B 1 and B 2 only. Next, we turn to the beta functions of the scalar sector. The anomalous dimension for the scalar field can be expressed as where the λ and y dependent coefficients are given by (3.20) The system is completed by the beta function for the higher-derivative coupling y. Its general structure follows a similar pattern as η s : The functions S 5 to S 8 depend on λ and y and are found to be 1+y − 25 + 85 y , Eqs. (3.14), (3.15), (3.19) and (3.21) form an implicit system which can be solved for the beta functions β λ , β y and anomalous dimensions η N and η s . In absence of the higher-derivative terms in the scalar propagator, which can be switched off by setting y = 0 and β y = 0, the beta functions agree with the ones reported in [15]. This provides a non-trivial crosscheck of our derivation.

Structural properties of the beta functions
The system of beta functions (3.14), (3.15), (3.19) and (3.21) possesses several interesting properties. Firstly, η s and β y depend on the number of scalar fields N s only implicitly. This feature is readily deduced from the Feynman diagrams in Fig.  1 which do not contain closed scalar loops that could give rise to terms proportional to N s . The number of scalars then enters the flow in the scalar sector only indirectly through the value of the cosmological constant and the anomalous dimension of Newton's constant. This suggests that the fixed point structure and flow pattern obtained from the beta functions will be rather stable under a change of the number of scalar fields.
Moreover, the beta functions possess several singular loci where either a beta function or an anomalous dimension diverges. The projection of these singular lines onto the y = 0-plane is shown in Fig. 2. Inspecting β λ and β y one encounters two singular lines λ sing = 1 2 and y sing = −1 , where the denominators in the beta functions vanish. 8 In addition one obtains singular lines when the anomalous dimensions η N or η s develop a singularity. For η N this locus is independent of y and N s and implicitly parameterized by the relation Since B 2 (λ) depends on the choice of coarse-graining operator, there are two distinguished structures entailed by this relation. As illustrated in Fig. 2 the Type I choice leads to a singular locus which screens the line λ sing = 1 2 for positive Newton's constant while the Type II coarse graining screens λ sing = 1 2 for g < 0. This observation may actually become important when "quenching the cosmological constant" along the lines proposed in [91] which presupposes that an RG trajectory emanating from the classical regime can actually reach the singular locus λ sing = 1 2 . The hypersurface on which the scalar anomalous dimension η s diverges is given by a quadratic polynomial in g with λ and y-dependent coefficients η sing For y = 0 the resulting line is depicted as the purple line in Fig. 2. The hypersurface also screens the line λ sing = 1/2 for g > 0. In the Type I coarse graining procedure η sing s is sandwiched between η sing N and λ sing = 1/2, while for the Type II procedure, it actually provides the screening of the λ sing = 1/2-line. Thus we see that the inclusion of scalar matter actually alters the singularity structure of the beta functions. At the same time, we expect that the system is rather insensitive to the inclusion of matter fields. The later point will be confirmed in more detail by the analysis of the next section.

Properties of the renormalization group flow
We now discuss the properties of the RG flow entailed by the system (3.14), (3.15), (3.19) and (3.21). In Sect. 4.1 we study the flow of the subsystem where the effects of the higher-derivative terms are switched off. The results provide the basis for analyzing the effects related to the presence of higher-derivative terms in the scalar propagator in Sects. 4.2 and 4.4. Throughout the section we focus on the flow generated by the choice (3.16), restricting ourselves to the discussion of a Type I coarse-graining operator only.

Minimally coupled scalar fields
The system (3.13) constitutes a set of autonomous coupled first order differential equations capturing the scale-dependence of {g k , λ k , y k }. Before delving into the analysis of the full system, it is useful to first analyze the subsystem obtained from setting y k = 0, β y = 0. In this approximation the contributions of the higher-derivative terms in the scalar sector are switched off and the projection of the flow equation is given by the Einstein-Hilbert action supplemented by an arbitrary number N s of minimally coupled scalar fields. The RG flow resulting from similar projections has been studied in [15,[17][18][19]21]. The analysis of this subsection then facilitates the comparison with these works.
Fixed point structure. The reduced system possesses two fixed points, a Gaussian and a non-Gaussian one. The Gaussian fixed point (GFP) is situated in the origin and its stability coefficients are determined by the mass-dimension of the coupling constants, (λ * , g * ) = (0, 0) , θ 1 = 2 , θ 2 = −2 . In the interval N s ∈ [−4, 16] the NGFP discussed above is the only non-trivial fixed point solution. Outside this window the simplified system possesses additional NGFPs. These are, however, located outside the physically interesting region located at g > 0 and to the left of the singular lines depicted in Fig. 2. Therefore, these fixed points will not be discussed in detail. Flows away from the NGFP. Beyond the vicinity of the NGFP, where the linearized approximation of the flow is valid, the RG trajectories can be constructed by integrating the beta functions of the reduced system numerically. In the case where the critical exponents of the NGFP are complex (N s ∈ [−6, 93]) the resulting phase diagram follows the same classification as in the case of pure gravity [13]. For the case N s = 1 three prototypical RG trajectories are shown in the left diagram of Fig. 4. The trajectories undergo a crossover from the NGFP, controlling the highenergy regime, to the GFP, controlling the classical regime of the theory. The RG trajectory connecting the two fixed points is called "Type IIa" and leads to a vanishing value of the renormalized cosmological constant lim k→0 Λ k = 0. Trajectories flowing to the left (right) to this line are called Type Ia (Type IIIa) and give rise to a negative (positive) value of the cosmological constant in the classical regime. The present set of flow equations do not allow to continue the Type IIIa solutions to k = 0: they terminate in the line η sing N shown in Fig. 2 at a finite value of k. The scalar anomalous dimension obtained along these sample RG trajectories is shown in the right panel of Fig. 4. Notably η s (k) ≤ 0 along the entire flow: at the NGFP one has η * s = −0.771 and the scalar anomalous dimension approaches zero when the flow enters the classical regime governed by the GFP. Thus the anomalous dimension induced by the gravitational quantum corrections suppress the propagation of scalar modes on all scales. The rapid increase of |η s | for the Type IIIa trajectory close to its termination point is a clear indication that the present approximation is insufficient in this regime and should thus not be given too much significance.

Fixed point structure including higher-derivative terms
We now focus on the fixed point structure of the full system (3.13) including the higher-derivative coupling y k . Following the structure of the last subsection, we first discuss the fixed point structure of the system.  Figure 4: Three prototypical RG trajectories obtained from numerically integrating the reduced system of beta functions for N s = 1 (left). The flow is governed by the interplay of the NGFP and GFP. The scalar anomalous dimension η s along the trajectories is shown in the right diagram. The initial scales k 0 are tuned such that the trajectories in the right diagram are disentangled. Notice that η s is negative semi-definite along the entire RG flow. In the UV (k → ∞) the anomalous dimension η s approaches its fixed point value η * s = −0.771 independently of the specific initial conditions. In the IR η s remains negative and vanishes asymptotically for the solutions of Type Ia and Type IIa. Trajectories of Type IIIa terminate in the singular line η sing N triggering the divergence of η s at a finite value of k.
Inspecting the beta functions, one finds that the GFP (4.2) has the following extension (λ * , g * , y * ) = (0, 0, 0) , θ 1 = 2 , θ 2 = −2 , θ 3 = −2 . Again there is a GFP for all values of N s and the anomalous dimensions vanish at this fixed point. The stability coefficients indicate that the GFP is a saddle point in the λ-g-y-plane exhibiting one UV-attractive and two UV-repulsive eigendirections.
In particular, it may serve as an IR attractor for RG flows starting at g k > 0 which subsequently leave the GFP regime along the unstable direction. The analysis of possible NGFPs starts with the following, intriguing observation: when restricted to y = 0, the beta function β y , given in eq. (3.21), simplifies to Thus β y supports a fixed point at y * = 0 if η * N = −2. From β g one finds that the latter condition is precisely the anomalous dimension of Newton's coupling at any NGFP. This shows that there is an extension of the NGFP discussed in the previous section to the full system, i.e., for all values of N s we obtain a NGFP with y * = 0. This family of NGFPs will be called NGFP 0 in the sequel. Remarkably, the balancing between the anomalous dimension η * N and the other contributions to β y works for d = 4 only. In any other spacetime dimension the fixed point is shifted away from the y = 0-plane.
A numerical investigation of the fixed point structure for N s ∈ [−200, 350] reveals the existence of 3 families of NGFPs, parameterized by N s , and located in the physically interesting region. The three families are conveniently labeled by the sign of the fixed point value y * which is either negative (NGFP − branch), zero (NGFP 0 branch), or positive (NGFP + branch). The positions and stability coefficients of these fixed  points are shown in Fig. 5. In addition the characteristics for the NGFPs found for N s = 1 are collected in Table 2. The detailed properties of the fixed point solutions are the following.  At this point the following remark is in order. Combining eqs. (2.6) and (3.11), the mass of the Ostrogradski ghost is Thus µ 2 will become infinite for any RG trajectory approaching a NGFP as k → ∞. This is just a consequence of the fact that a fixed point can not support a dimensionful scale. The relation (4.5) also reveals that the fixed points NGFP 0 are very special. Owed to their position at y * = 0 the mass of the Ostrogradski ghost is infinite for all values k. In this way, the NGFP 0 realize the first class of loopholes discussed in sect. 2.2. Thus the extra degree of freedom is not present and one expects that the resulting theory does not suffer from an Ostrogradski instability albeit living in a theory space which permits the presence of higher-derivative terms in the propagator a priori.

Phase diagram including higher-derivative terms
We now extend the local analysis of the RG flow, based on its fixed point structure and stability coefficients, to a global picture. For concreteness, we focus on the case N s = 1. The details of the fixed point structure arising in this setting is summarized in Table 2. Since the essential features of the flow are set by its fixed point structure, it is clear that the analysis applies to an entire window −6 N s 12 where the fixed point structure and stability coefficients exhibit the same qualitative behavior.
The global structure of the RG flow is obtained by integrating the beta functions (3.13) numerically. A characteristic set of trajectories obtained this way is shown in Figs. 7 and 8. Fig. 7 then shows the RG trajectories connecting the 3 NGFPs (gray lines) and the NGFPs with the GFP (blue lines). Since both NGFP ± act as UV-attractors in the λ-y-g-plane and the NGFP 0 possesses one IR-attractive eigendirection there is a single RG trajectory emanating from either NGFP ± for k → ∞ and ending at the NGFP 0 as k → 0. The GFP possesses 2 IR-attractive eigendirections. As a result, one finds a unique trajectories which starts from NGFP 0 and connects to the GFP k → 0 (light blue line). This trajectory is the intersection of the two-dimensional UV-critical hypersurface of NGFP 0 with the two-dimensional IR-critical hypersurface of the GFP. In addition there are two families of solutions which originate from NGFP ± and end at the GFP, again coming from the intersection of the 3-dimensional UV-critical hypersurfaces of the NGFPs with the IR-critical hypersurface of the GFP. These flows are exemplified by the dark blue lines. All together this set constitutes the generalization of the Type IIa trajectory displayed in Fig. 4. Fig. 8 then illustrates the generalization of the trajectories of Type Ia and Type IIIa to the λ-y-g-plane. These trajectories may emanate from all three NGFPs and subsequently cross over to the GFP. From the vicinity of the GFP they either flow to large negative values λ k (Type Ia) or positive λ k (Type IIIa) such that their projection to the λ-g-plane resembles the left diagram of Fig. 4. The latter class again terminates in the hypersurface η sing N at a finite value k. Notably, for all physically interesting trajectories which exhibit a crossover to the GFP, y k flows to zero in the IR, provided that the underlying trajectories do not terminate at a finite value k. When evaluating the scalar anomalous dimension η s along the RG trajectories shown in Figs. 7 and 8 one again obtains the qualitative behavior shown in the right diagram of Fig. 4: for large values of k η s is determined by its fixed point value η * s . Once the RG trajectory enters the vicinity of the GFP quantum effects become small, η s 1 asymptotically.

Ghost-free RG flows in the infrared
In order to determine the stability of the theory in the presence of higher-derivative terms one has to study the renormalized scalar propagator obtained from the effective average action Γ k in the limit k → 0. Defining Y 0 ≡ lim k→0 Y k the (squared) mass of the Ostrogradski ghost is (cf. eq. (2.6)) Hence instability will disappear from the spectrum if Y 0 = 0. Thus the focus of the investigation is on the IR behavior of y k . Fig. 7 demonstrates that all physically interesting RG trajectories have the property that the dimensionless coupling y k goes to zero in the IR. This leaves three potential scenarios for the dimensionful coupling Y k = y k k −2 : 1. The dimensionless coupling y k approaches zero slower than quadratically. The canonical scaling of Y k will dominate the flow and Y 0 diverges. In this case the ghost becomes massless and eats up the scalar degree of freedom, see eq. (2.5).
2. The dimensionless coupling falls off faster than k 2 . The anomalous scaling dominates the flow, and Y k → 0. The Ostrogradski ghost decouples and the theory is stable.
3. The dimensionless coupling converges exactly quadratically. The dimensionful coupling Y k approaches a constant, which can be either zero or nonzero. The theory is stable only if this constant is zero.
We will now discuss the IR behavior of the several classes of trajectories. Most of the physically interesting trajectories fall into the classes Type Ia, Type IIa, or Type IIIa introduced in Fig. 4. The only trajectories which are not captured by this classification are the trajectories connecting the NGFPs which will be discussed separately. Our investigation reveals that the phase diagrams shown in Figs. 7 and 8 realize all of the three cases described above.
Trajectories ending at the GFP (Type IIa). We start our analysis by considering Type IIa trajectories for which the cosmological constant Λ k flows to zero for k → 0. In this case the IR completion of the trajectory is provided by the GFP (4.3). The IR attractive hypersurface of the GFP is spanned by the two eigenvectors associated with the negative stability coefficients θ 2 = θ 3 = −2. The explicit expression for these eigenvectors are e 1 =ŷ and e 2 = 2+Ns 16πλ +ĝ, whereŷ,λ andĝ are the unit vectors along the y, λ and g-axis, respectively. By linearizing the flow at the GFP one finds that along these scaling directions Hence, there is a single RG trajectory, specified by Y k 0 = 0, for which Y 0 = 0 and the mass of the Ostrogradski ghost becomes infinite. This is the trajectory that has no initial component in theŷ-direction, i.e. the one that approaches the GFP along e 2 . Integrating the beta functions numerically one finds that this trajectory belongs to the UV-critical hypersurface of NGFP − .
Trajectories of Type Ia and IIIa. Fig. 4 illustrates the existence of RG trajectories where λ k flows towards negative or positive infinity as k → 0. The corresponding solutions are then classified as trajectories of Types Ia and IIIa, respectively. In order to determine the IR behavior of these trajectories, we numerically integrate the beta functions. Trajectories of Type IIIa terminate at η sing N at a finite value of k and can not be completed to k = 0 in the present approximation. Therefore, we limit our analysis to trajectories of Type Ia which extend up to k = 0. The IR values Y 0 ≡ lim k→0 Y k arising within this class of solutions are conveniently illustrated by studying the behavior of RG trajectories piercing the y-g-plane located at λ = −0.1 since the flow is essentially perpendicular to this plane. The resulting structure is illustrated in Fig. 9. The plot shows that Type Ia trajectories can emanate from all three NGFPs: trajectories coming from NGFP 0 pass the plane at the blue line while trajectories above (below) this line lie in the UV-critical surface of NGFP + (NGFP − ). Trajectories where Y 0 = 0 span the black line in this diagram. Thus there is a 1-dimensional surface of solutions where the renormalized squared mass of the Ostrogradski ghost, (4.6), is infinite such that the resulting degree of freedom does not propagate. Imposing the physical requirement that the renormalized scalar propagator does not give rise to an Ostrogradski ghost may then be used to fix one of the free parameters of the theory from stability considerations.
Trajectories flowing to NGFP 0 . The final option for taking an IR limit consists in approaching NGFP 0 along its IR-attractive eigendirection. From Fig. 7  which end at NGFP 0 as k → 0. Linearizing the RG flow at the NGFP 0 and using the stability coefficient along the IR attractive eigendirection listed in Table 2 yields the RG evolution of y k for these trajectories: Since the scaling of the dimensionless y is significantly smaller than k 2 , the dimensionful Y diverges as k → 0 for all initial values y = 0. As a consequence the IR value of the ghost mass vanishes and the two terms describing the propagation of the scalar field in eq. (2.5) mutually cancel. Loosely speaking, the physical degree of freedom is eaten by the ghost so that the scalar does not propagate anymore. Verifying the robustness of this cancellation-mechanism requires the inclusion of further powers p 6 , p 8 , . . . in the scalar propagator. This analysis is beyond the scope of the present work, however, and will be addressed in a forthcoming publication [94].

Conclusions and outlook
In this work, we use the effective average action Γ k to study the renormalization group flow of gravity coupled to scalar matter. Our ansatz for Γ k is given by the Einstein-Hilbert action coupled to an arbitrary number of minimally coupled scalar fields. The novel feature of the setup is the inclusion of a higher-derivative term in the scalar propagator. At the classical level these types of actions suffer from the so-called Ostrogradski instability reviewed in Sect. 2: the appearance of degrees of freedom with a wrong-sign kinetic term, so-called Ostrogradski ghosts, renders the theory either unstable or non-unitary. At the same time it is clear that a generic RG flow will generate such potentially dangerous higher-derivative terms dynamically. This work initiates the systematic study of these types of terms in the RG framework with the goal of assessing their hazard potential for asymptotically safe theories. The quantity that actually encodes the relevant information on the spectrum of the theory is the renormalized propagator. Exploiting that the effective average action obeys lim k→0 Γ k = Γ with Γ being the standard effective action, this quantity can be accessed in the IR-limit of the flow. Within the present approximation the stability properties of the theory are captured by the IR-value of the (squared) Ostrogradski ghost mass µ 2 = Y −1 0 . The ghost decouples from the spectrum if Y 0 = 0, so that the setting may give rise to stable (or equivalently unitary) theories even though the generic actions include higher-derivative kinetic terms.
The detailed study of the RG flow then established the following picture. In absence of the higher-derivative term the setting gives rise to a unique non-Gaussian fixed point (NGFP) suitable for rendering the gravity-matter system asymptotically safe. Upon including the scale-dependent Ostrogradski ghost mass, this NGFP splits up into three NGFPs which are labeled by the sign of Y * . Notably there is one fixed point solution NGFP 0 for which Y * = 0 for all values of k.
When projected to the λ-y-g-plane (see Figs. 7 and 8) the system of NGFPs essentially possesses a UV-critical hypersurface with three relevant directions. Within this space we have identified a two-dimensional subspace of RG trajectories that have a ghost-free IR limit. Phrased differently, the Ostrogradski ghost mass corresponds to a relevant direction of the NGFPs coming with a new free parameter. This freedom can be fixed by the requirement that the theory should contain only physical degrees of freedom in the IR. In this way the construction elegantly circumvents the potential danger of Ostrogradski instabilities by introducing a new free parameter and a mechanism to fix its value simultaneously. The analysis in Sect. 4 shows that the set of complete, unitary RG trajectories obtained from the full λ-g-y-system (3.13) is in one-to-one correspondence with the one found in the reduced system excluding the higher-derivative coupling.
As a byproduct, our analysis also provided new insights on potential bounds on the number of scalar fields compatible with the asymptotic safety mechanism. Throughout the calculation, we used a coarse graining operator of Type I (see [2] for an extended discussion), and extracted the running of η N from the background Newton's constant. The resulting analysis indicates that there are NGFPs suitable for realizing asymptotic safety for all values N s . The characteristic fixed point properties shown in Fig. 3 are strikingly similar to the ones found for foliated gravity-matter systems [19]. Notably, our results also agree with the ones reported in [15], where an upper bound N s 17 has been obtained. The crucial difference between the two settings lies in the choice of coarse-graining operator in the gravitational sector: our analysis uses a Type I coarse-graining operator while [15] resorts to a coarse-graining operator of Type II. If the analysis of Sect. 4.1 is repeated for a coarse-graining operator of Type II, which effectively replaces eq. (3.16) by (3.17), the reduced system (3.13) gives rise to the same upper bound on the number of scalar fields N s 17. From Fig. 2 one then expects that the singular line η sing N plays a decisive role in stabilizing the NGFP for large values N s .
Our analysis demonstrates that the existence of unitary RG trajectories is a non-trivial feature. A priori, a kinetic function of polynomial type is bound to have multiple roots, yielding a ghost in the particle spectrum. Further investigations suggesting themselves include studying a) polynomial truncations including further powers of the momentum, p 6 , p 8 , . . ., or b) truncations of non-polynomial type. In the first case, a higher-order truncation allows to investigate whether RG properties in lower orders are stable. In the second case, non-polynomial kinetic functions open up the possibility to have analytic kinetic functions without multiple roots, giving a ghost-free spectrum. An example is a propagator of the type e −∆ (∆ + m 2 ) −1 studied, e.g. in the context of non-local gravity models [73][74][75][76]. 10 The unitarity conditions on such kinetic functions are studied in a separate paper [69] and we hope to come back to the other points in the near future as well.

A Expanding trace arguments including step functions
In this appendix, we collect the technical details underlying the derivation of the beta functions in the scalar sector. In this case, it is most convenient to choose a flat background spacetime wereḡ µν = δ µν . This allows to use momentum space techniques to evaluate the diagrams shown in Fig. 1.

A.1 Explicit form of vertex functions and propagators
We start by deriving the relevant propagators and interaction vertices from the ansatz for the effective average action (3.3). The result is conveniently expressed in terms of the variations Γ (k,l;m) k where the number of derivatives with respect to the metric fluctuations and scalar fluctuations are denoted by k and l, respectively. The number m denotes the number of remaining background scalar fields. Moreover, we use the index w to specify whether the building block is associated with the graviton (w = hh) or scalar fluctuations (w = φφ).
By expanding the gravitational sector up to second order in h µν one finds that the (inverse) gravitational propagator is given by is the unit on the space of symmetric tensors and [P h ] µν αβ ≡ d −1 δ µν δ αβ the projector on the trace mode. The (inverse) scalar propagator is obtained from (3.7) and reads For later convenience, we introduce the following short-hand notations for the scaledependent coefficients α w n multiplying the p 2n terms in the (scalar part) of eqs. (A.1) and (A.2), and all coefficients α w n with n ≥ 3 vanishing. In addition to the propagators, one also needs the (momentum-dependent) threeand four-point vertices containing one and two derivatives with respect to the background scalar field. Denoting the momenta associated with the graviton fluctuations, scalar fluctuations, and background scalar field byp, p, and q, respectively the 3-point vertex obtained from (3.7) is Finally, the 4-point vertex is All vertices are understood to contain the appropriate symmetrizations in the external indices and are subject to momentum conservation. Moreover, we set Y k = 0 in order to keep the expressions for the vertices at a readable length. The contributions proportional to Y k are easily generated by a computer algebra program. Their precise form is irrelevant for the discussion of the general structures below. Applying the implicit regulator prescription p 2 → P k = p 2 + R k (p 2 ) to the propagators (A.1) and (A.2) For the Litim-type cutoff [92,93] the dimensionful profile function R k is given by The key advantage of this regulator is that it allows for an analytic evaluation of the loop integrals shown in Fig. 1. The distributional character of the regulator renders the expansion in the external momenta q non-trivial, however. The next subsection discusses how this expansion can be implemented consistently, also taking into account the non-trivial boundary terms arising in the expansion procedure.

A.2 Loop-integrations with a distributional regulator
The loop integrals entailed by Fig. 1 contain a trace over spacetime indices and an integration over loop momenta. 11 The spacetime indices are taken into account by stringing together the propagators and vertices contracting the corresponding index structures. This results in q-dependent scalar loop-integrals of the form l=0 α w 2 l p 2l + R w 2 k (p 2 ) n .
(A.8) Here m and n encode the number of propagators appearing in the diagram and, in a slight abuse of notation, the symbol R hh k (p 2 ) is used to refer to the scalar part of (A.6a). Diagrams containing 3-point vertices have (m, n) = (1, 2) while the tadpole diagram comes with (m, n) = (0, 2). The function F k (q, p, cos(ϑ)) captures the momentum dependence of the vertices and is polynomial in q and p. In particular it has a well-defined series expansion around q = 0. Noting that the vertices (A.4) and (A.5) come with one and two powers of the external momentum, respectively, it is easy to verify that this expansion starts at order q 2 .
For a general profile function R k the integrals eq. (A.8) cannot be computed analytically. Moreover, the presence of the external momentum q and the scaledependent couplings make their numerical evaluation computationally very expensive. The profile function (A.7) allows to bypass this problem by restricting the p-integration to a compact domain and giving rise to cancellations in the propagators. The former property can be verified by noting that the logarithmic k-derivative of (A.6), evaluated for a Litim profile, has the form (A.10) Inspecting (A.8) for the case m = 0 (tadpole diagram) reveals that the stepfunctions appearing in the numerator and denominator have the same support. As a result the integrals simplify significantly I w 1 w 2 (0,n) ≡ˆdΩˆ1 Here´dΩ denotes an angular integration and the spacetime indices on F k andb w k (p 2 ) are suppressed for readability. Owed to the simple structure of the denominator, which is independent of p and ϑ the evaluation of these integrals is rather straightforward.
The case where m = 0 is non-trivial, however. Owed to the step-function in the numerator the full integration domain is reduced to a d-dimensional ball of radius k, i.e., p ∈ [0, k] and cos(ϑ) ∈ [−1, 1]. In this domain the second set of propagators again undergoes the simplification (A.11). In the first set of propagators the regulator leads to terms proportional to Θ(k 2 −( p+ q) 2 ), however. As illustrated in Fig. 10, the value of the step function has a non-trivial dependence on the absolute value of q and the angle ϑ. Thus, unless q = 0, there is always a part of the integration domain on which the denominator does not become trivial. As a result performing the integral becomes very involved. In order to complete the evaluation of the flow equation we then expand the integrands around q = 0, taking the distributional character of the integrand into account. This allows us to obtain analytic expressions for the resulting integrals. This is achieved as follows. encounter terms in which the delta-distribution has to be evaluated on the boundary of the integral domain. Using Θ(0) ≡ 1 2 , these can be evaluated by noting that

A.3 Master integrals
We close the discussion by deriving a set of master integrals, which form the basis of our loop computations I w,m (q, cos(ϑ)) ≡ˆk 0 dp (2π) d f (p) Θ(( p+ q) 2 −k 2 ) This result completes the discussion on carrying out the momentum integrals entailed by Fig. 1. Note that the surface terms do not enter into the computation of the scalar anomalous dimension. They contribute to higher-order kinetic terms in the propagator only.