A proper fixed functional for four-dimensional Quantum Einstein Gravity

Realizing a quantum theory for gravity based on Asymptotic Safety hinges on the existence of a non-Gaussian fixed point of the theory's renormalization group flow. In this work, we use the functional renormalization group equation for the effective average action to study the fixed point underlying Quantum Einstein Gravity at the functional level including an infinite number of scale-dependent coupling constants. We formulate a list of guiding principles underlying the construction of a partial differential equation encoding the scale-dependence of $f(R)$-gravity. We show that this equation admits a unique, globally well-defined fixed functional describing the non-Gaussian fixed point at the level of functions of the scalar curvature. This solution is constructed explicitly via a numerical double-shooting method. In the UV, this solution is in good agreement with results from polynomial expansions including a finite number of coupling constants, while it scales proportional to $R^2$, dressed up with non-analytic terms, in the IR. We demonstrate that its structure is mainly governed by the conformal sector of the flow equation. The relation of our work to previous, partial constructions of similar scaling solutions is discussed.


Introduction
Renormalization group (RG) fixed points are crucial for determining the properties of the underlying quantum field theory. They are lighthouses on theory space because in their vicinity one can explore the corresponding quantum field theory via scaling techniques and, eventually, via conformal field theory methods. Fixed points also play a crucial role in providing a general notion of renormalizability and ultraviolet completion of a quantum theory. In fact, in the vicinity of an ultraviolet (UV) fixed point the RG flow can be linearized and characterized by directions along which it is attracted to or repelled from the fixed point for increasing energies, thus providing a generalized notion of relevant and irrelevant deformations of the theory. The UV completeness of the theory is then ensured by finding a given RG trajectory that is attracted towards the fixed point at high energies, implying that any UV cutoff can be removed in such a way that all essential dimensionless coupling constants remain finite. Such an UV complete theory is also predictive if the dimension of the critical hypersurface that is attracted towards the fixed point in the UV is finite. In this case, only a finite number of experiments is required to pinpoint the specific RG trajectory that completes gravity in the deep UV regime. Depending on whether the fixed point corresponds to a free (Gaussian fixed point) or interacting theory (non-Gaussian fixed point), the RG trajectories terminating at the fixed point in the UV are termed asymptotically free or asymptotically safe, respectively. While it is well known that the Gaussian fixed point (GFP) of gravity does not provide a suitable mechanism to renormalize the theory, the possibility that the theory possesses at least one suitable non-Gaussian fixed point (NGFP) is not speculative and has been explored under various convincing approximations. In this work we investigate the Asymptotic Safety mechanism for gravity relaxing some of these approximations. In particular we construct the first consistent and globally well-defined fixed functional for Quantum Einstein Gravity in four spacetime dimensions.

Quantum Gravity via Asymptotic Safety
The RG flow of any theory admits a GFP at which the theory itself becomes non-interacting. In the neighborhood of the Gaussian point it is thus possible to use the methods of perturbation theory to investigate the UV completion of the theory under consideration. These methods have been applied very successfully to Yang-Mills theories thanks to the property of asymptotic freedom. However, it has been known for many years that the GFP cannot provide a UV completion mechanism of General Relativity [1][2][3]. The question whether gravity can be formulated as a consistent and predictive (asymptotically safe) quantum field theory then becomes closely intertwined with the question whether there are other fixed points of the RG on the gravitational theory space which could provide such a UV completion. Indeed it has been conjectured already in the late seventies that there exists a suitable NGFP on the gravitational theory space which could render gravity asymptotically safe [4].
The primary tool for investigating the Asymptotic Safety mechanism has been a Wilsonian-type functional renormalization group equation (FRGE) [5][6][7] k∂ k Γ k = 1 2 Tr Γ (2) in which Γ k denotes an effective "average" action, Γ (2) k is its second variation with respect to the fields, and R k is a suitable infrared-regulator which ensures that the rhs of the equation is finite and peaked around the renormalization group scale k. The effective average action Γ k represents an effective action of the system in which the UV modes have been integrated out in the path-integral using the scale k as a reference, thus providing a natural effective action for processes occurring at energies E ≃ k. In the limit k → 0 the IR cutoff R k vanishes and the effective average action coincides with the full effective action of the system Γ = Γ k=0 . The flow (1.1) is exact in the sense that no approximation is used in its derivation.
Notably, the FRGE allows to compute approximate RG flows of a theory without the need of a small expansion parameter. For gravity this approach has been pioneered in [8] which, by now, has evolved in a very successful research program, see [9][10][11][12] for reviews. The key strategy for obtaining non-perturbative information on the RG flow is to project the exact flow on a subspace spanned by the "most relevant interactions" for the physics at hand. Practically, this strategy is implemented by making a suitable ansatz for Γ k limiting the set of interaction terms to those which are argued to describe the relevant physical information. The beta functions, which describe the scale-dependence of the coupling constants contained in the ansatz and drive the approximated RG flow, are constructed by substituting the ansatz for Γ k into the FRGE and retaining the terms contained in the ansatz itself only. While these beta functions lack the control parameter of the loop expansion of perturbation theory [13], they are fully non-perturbative in the couplings and thus, for sufficiently big truncations, reliable also away from the GFP.
The projection of Γ k to finite dimensional subspaces has led to a substantial body of evidence supporting the existence of a gravitational NGFP suitable for Asymptotic Safety. Starting from the Einstein-Hilbert truncation, which contains scale-dependent Newton's and cosmological constants only [8,14,15], the existence of the NGFP has successively been established on truncations of increasing complexity and field content, see [16][17][18][19] for selected references. For the case in which the gravitational degrees of freedom are carried by a (Euclidean) spacetime metric, this line of research has already established striking results. Foremost, all computations corroborated the existence of a NGFP with suitable properties for Asymptotic Safety. In particular, the NGFP has been put to test and proved robust to the inclusion of the problematic operators appearing in the perturbative counter-terms in Γ k [20,21], thus directly showing the feasibility of the program beyond perturbation theory. Despite this fixed point being associated with an interacting quantum field theory in which the (background) Newton's constant scales with an anomalous dimension η N = −2, classical power-counting may still constitute a good ordering principle for identifying the relevant deformations of the fixed point as shown in [22][23][24][25]. In particular the number of relevant parameters of the theory could be as low as three, implying that, compared to many other quantum field theory models, an asymptotically safe theory of gravity can be very predictive. Moreover, stringent bounds on the numbers of free matter fields compatible with the Asymptotic Safety mechanism have been obtained in Refs. [26,27]. At the level of the formalism, there has also been progress on computing properties and expectation values of the fluctuation fields around a given background [28][29][30][31][32]. By now it is appreciated that this type of computations, known as bi-metric computations, may be essential for precision computations of the properties of NGFP such as its critical exponents, addressing global aspects of gravitational RG flow by, e.g., establishing a suitable c-theorem [33], or establishing background-independence of the formalism [34].
While the case where the gravitational degrees of freedom are carried by metric variables is understood the best, this scenario is not the only way to realize Asymptotic Safety in gravity: a suitable NGFP has been found to occur rather independently of the type of gravitational degrees of freedom. Investigations of the gravitational RG flows based on the ADM-formalism [35,36], the Einstein-Cartan formalism [37][38][39][40] and the "tetrad only" formalism [41,42] also identify a NGFP whose properties are strikingly similar to the one encountered in the metric formalism. This is rather remarkable since, while these formulations are equivalent at the classical level, it is far from clear that this equivalence also holds at the quantum level. Moreover, there are also first indications that Asymptotic Safety may be realized in unimodular gravity [43][44][45]. This leads to the perhaps surprising conclusion that, at the level of finite-dimensional projections, the appearance of a NGFP in gravitational theories is rather commonplace instead of exceptional. 1

From fixed points to fixed functions
Currently, one of the key conceptual challenges in the Asymptotic Safety program is the leap from approximations for Γ k containing a finite number of scale-dependent couplings to truncations containing an infinite number of coupling constants. In the latter case, Γ k contains scale-dependent functions. Substituting such an ansatz into the exact functional renormalization group equation leads to a (system of) partial differential equations (PDE) which govern the scale-dependence of these functions. The fixed points appearing in finitedimesional approximations of Γ k are then promoted to fixed functions, namely global kstationary solutions of the PDE. (Since both fixed points and fixed functions induce fixed functionals, we use the notion of fixed functions in order to stress that we work in the realm of infinitely many scale-dependent coupling constants.) The simplest setting where these generalizations can be implemented for gravity includes an arbitrary function of the scalar curvature R where g denotes the (Euclidean) spacetime metric. This so-called f (R)-truncation provides the simplest example for a truncation including a scale-dependent function, since the PDE governing the k-dependence of f k (R) can be derived by using a maximally symmetric background with covariantly constant curvature. Starting from the flow equation derived in [23], the existence of fixed functions f * (R) have been investigated by various groups in d = 3 [52][53][54] and d = 4 [55][56][57][58][59] spacetime dimensions. Quite unsettling, the verification of a suitable NGFP at the level of fixed functions turned out to be extremely challenging: while the finite-dimensional computations always produced a suitable UV fixed point regardless of the computational details and setting, up to now the fixed function completing the four-dimensional NGFP in the ansatz (1.2) is still elusive. This deficit acted as catalyst to readdress some of the fundamental questions in the functional renormalization group approach to Quantum Einstein Gravity (QEG).
Substituting the ansatz (1.2) into the FRGE and projecting on the corresponding subspace leads to a non-linear PDE which is first order in k and third order in the Rderivatives. Imposing the condition that the fixed function is k-stationary, this equation reduces to a non-linear ordinary differential equation (ODE) of third order. Thus, locally, there is a three-parameter family of solutions to this equation. It was then noticed earlyon that the condition for the solution to exist globally may reduce the number of free parameters and may produce isolated fixed functions if there is a balance between the order of the ODE and its singularity structure [56] degree of ODE − number of singularities = 0 . (1.3) Thus any singularity reduces the number of free parameters by one. While the number of fixed singularities can essentially be read off from the ODE, determining the occurrence of movable singularities often requires studying explicit solutions. Thus, so far, only fixed singularities played a role in restricting the number of free parameters in the fixed function in the f (R)-truncation. As a important consequence of (1.3) it was then understood that the initial flow equation contains too many fixed singularities in order to admit global solutions [55,56]. Based both on physical arguments and on the experience maturated during the course of several f (R)-type analysis, we propose the following guiding principles towards the construction of a suitable fixed function in four-dimensional QEG, which will turn out to be a set of sufficient conditions for the existence of said fixed function: 1) The background should admit only one topology. The initial flow equation has been derived on a maximally symmetric, compact background sphere S d with positive scalar curvature. A detailed comparison of the flow equations for f k (R) in the three-dimensional conformally reduced setting established that the topology of background (non-compact with negative scalar curvature or compact with positive scalar curvature) manifestly influences the singularity structure of the PDE. Additionally, the flow equation on a hyperbolic background with negative scalar curvature does not follow from the analytic continuation of the flow equation on a sphere. We thus argue that it is potentially inconsistent to include more than one topological sector in the flow equation avoiding topology changes in the background. In the sequel, we restrict the construction to the domain R > 0 and do not demand that the solutions admit an continuation to the R ≤ 0 region, thereby avoiding that the fix functional has a suitable continuation into a domain with different topology.
2) The path-integral measure should not introduce spurious singularities.
As it turns out, the measure in the gravitational path integral has a crucial effect on the singularity structure of the ODE. Using a linear background field decomposition together with a standard measure for the metric fluctuations h µν and evaluating the flow equation based on the transverse-traceless decomposition of h µν generally introduces fixed singularities from the contributions of the gauge degrees of freedom, ghost terms, and Jacobians arising in the field decomposition. We believe that these singularities have to be considered as spurious, since their number can be changed at will by varying unphysical parameters. We will resort to the flow equation derived in [60] in order to avoid the spurious poles of the gauge-fixing sector and of the Jacobians of the field decompositions. This framework can be motivated by geometrical considerations, since it corresponds to a trivialization of the gauge-bundle in field space, and is equivalent to a particular simultaneous choice of regulator function R k and of gauge-fixing, which allows for a precise cancellation between the functional traces capturing the gauge-degrees of freedom and ghost contributions.
3) The asymptotic behavior is captured only by the exact heat-kernel.
The rhs of the FRGE (1.1) is a trace involving a (generally complicated) differential operator. In finite-dimensional truncations, this functional trace is typically evaluated by applying the early-time expansion of the corresponding heat-kernel. This is the natural approach if Γ k is organized via a derivative expansion, for example in powers of the background curvatures. It turns out, however, that at the level of full functional truncations the use of the early-time expansion and the evaluation of the trace based on the exact heat-kernel on S d leads to a different asymptotic behavior when k 2 ≪ R. This difference is caused by non-analytic contributions to the heat-kernel, related to the effect that on S d the diffusing particle may return to its origin multiple times. This effect is not captured by the (analytic) early-time expansion. These terms are, however, essential for reproducing the correct late-time behavior of the heat-kernel (and are yet another manifestation of the compactness of the background manifold). Therefore, we base the evaluation of the operator trace (1.1) on the exact heat-kernel including non-analytic contributions to correctly account its asymptotic behavior.

4) Scalar and TT-fluctuations should be integrated out simultaneously.
The initial flow equation has been derived using the Laplacian of the background sphere for the construction of a covariant regulator. In the nomenclature of [24], this corresponds to a type I cutoff. Our flow equation generalizes this construction by allowing a non-zero endomorphism in the regulator function, using a so-called type II cutoff. Besides making the solutions of the ODE more stable, the freedom of the endomorphism allows to satisfy both the singularity count (1.3) and the correct asymptotic behavior of the solution. Moreover, choosing different endomorphisms in the different sectors of the flow equation allows adjusting the scale at which the lowest eigenmode of the fluctuation field is integrated out. We thus need a physically motivated principle to fix this additional freedom. We will impose the principle of equal lowest eigenvalues [60], enforcing that the lowest energy mode in the transverse-traceless and scalar sectors are integrated out at the same value of k. Since the order of the highest derivative of f k (R) differs in the transverse-traceless and scalar sector of ODE, this choice in particular ensures that the order of the ODE does not change on the integration interval.

5) All quantum fluctuations need to integrated out.
A crucial question related to the existence of globally well-defined solutions is the integration range for which the solution ought to exist. The ODE describing the fixed function is typically formulated in terms of the dimensionless curvature r ≡R/k 2 . For fixed (background) curvatureR the interval k ∈ [0, ∞[ is then mapped to r ∈ [∞, 0], so it seems natural to demand that solutions are regular on the entire positive real axis. This is in agreement with the intuition from scalar field theory in a flat background and the requirement that the regulator R k ∝ k 2 vanishes for k 2 → 0, so that all quantum fluctuations are integrated out. Based on this logic the domain of definition for the fixed function may be inferred by imposing that a) all quantum fluctuations are integrated out, b) there are no residual regulator effects in the flow equation.
Depending on the particular choice of regulator, and in particular for Type II cutoffs in combination with the Litim regulator [61], these requirements may already be satisfied at a finite r max . In particular this construction removes the regulator effects not via the "canonical mechanism" sending k → 0 but by the stepfunction containing the eigenvalues of the differential operator becoming zero once the last fluctuation mode is integrated out.
In such a case, after the last eigenfluctuation is integrated out, the flow will be driven by the canonical scaling of the function f (R). Based on this mechanism the r-interval, on which the global solution is required to be regular, may be restricted to a finite interval r ∈ [0, r max ], where r max depends on the choice of endomorphism included in the regulator and on the regulator itself. If the flow equation is derived using a smooth regulator or the behavior of a smooth regulator is mimicked by applying a smoothing procedure the contribution of the quantum fluctuations are smeared over the entire half-axis r ∈ [0, ∞].
Since we apply the smoothing procedure [55] in the construction of our flow equation, we consequently require that a proper fixed function is globally well-defined on this domain.
In this work, we establish that these guiding principles lead to a PDE governing the scale-dependence of the ansatz (1.2) which gives rise to an isolated, unique and globally well-defined fixed function, thus reconciling previous results based on finite-dimensional polynomial truncations with a functional ansatz. While it could have been argued via the counting theorem that only a finite number of solutions was present, it is in fact amazing to discover that the solution resulting from the application of the above principles is unique. Let us stress again, however, that our guiding principles constitute a set of sufficient conditions for the existence of said fixed function. At this stage it is unclear whether all principles are strictly necessary for establishing a suitable fixed function. In particular, the list does not incorporate the recent observation [59,62] that the combination of a flow equation with the equation for the modified split Ward identities may lead to substantial simplifications of the PDEs prescribing the fixed functions of the theory.
The remaining paper is organized as follows. We start by illustrating the effects of the singularity theorem (1.3) together with the analytic and numerical techniques used to analyze the existence of fixed functions in a simple toy model in Sect. 2. Implementing the guiding principles discussed above, the ODE governing the existence of fixed functions in the f (R)-truncation is constructed in Sect. 3. The properties of the resulting ODE are discussed in Sect. 4 and a numerical study of its solutions is carried out in Sect. 5. The central result is the isolated, globally well-defined fixed function shown in Fig. 4. We end with a discussion in Sect. 6. The analogue of the full fixed function arising from the conformally reduced approximation of our flow equation is constructed in Appendix A.

Isolated global solutions: a toy model
Before embarking on the quest of constructing our flow equation for f (R)-gravity and its fixed functions, it is useful to illustrate the techniques used in the analysis by applying them to a simple, completely tractable example. For this purpose, we analyze the second order ODE , at x = 0, the equation admits a two-parameter family of solutions which we parameterize by (y 0 , y 1 ). Inspecting the rhs of eq. (2.1) one easily sees that the denominator has a root at x sing = 1 which constitutes a fixed singularity. A solution with generic initial conditions (y 0 , y 1 ) imposed at x = 0 will have a diverging second derivative at x sing . A necessary condition for regularity of the rhs is i.e., a globally well-defined y(x) has at least a simple root at x sing . Thus it is expected that not all solutions can be extended to global solutions which are well-defined on the entire positive half-axis. The fixed singularity provides a non-trivial boundary condition which puts additional constraints on the parameters (y 0 , y 1 ) characterizing the local solutions. We now discuss various techniques which allow to analyze this constraint.

The exact analytic solution
The linear ODE (2.1) is still simple enough so that its general solution can be given in terms of modified Bessel functions Evaluating the solution in the limit x → 1 − , the term containing I 1 vanishes while the  (2)) . Imposing the condition for a regular, global solution, c 2 = 0 and eliminating c 1 from the system (2.4) shows that the global solutions are given by a line in the (y 0 , y 1 )-plane, satisfying This line is displayed as the bold blue line in the right panel of Fig. 1. Thus the interplay between demanding that a solution to exist globally and the fixed singularity of the ODE reduces the number of free parameters by one.
Polynomial expansion at x = 0 Concerning the ODE describing fixed functions in quantum gravity, it is not likely that exact, analytic solutions can be found. Thus other techniques to derive the condition for a global solution (2.5) are required. A tractable method approximates the exact solution by a polynomial of finite order Substituting this ansatz into (2.1) and matching the coefficients multiplying powers of x up to x N −2 yields N − 1 equations relating the coefficients y n : This system can be used in two ways. Firstly, it can be complemented by the two equations appearing at orders x N −1 and x N : Combining eqs. (2.7) and (2.8) then shows that this completion scheme only leads to the trivial solution y n = 0, ∀ n ∈ 0, . . . , N . Alternatively, one can follow the strategy of fixing free coefficients by demanding the global existence of the solution. For this case, it is convenient to cast the system (2.7) into the form from which the y n can be determined as functions of the parameters y 0 , y 1 in a recursive manner. Fixing N = 12 and evaluating the regularity condition (2.2) for the polynomial (2.6) then yields the regularity condition This condition is illustrated by the thin green line in the right panel of Fig. 1 and is essentially identical to (2.5). Notably, fixing all boundary conditions at x = 0 (first strategy) or demanding the existence of a globally well-defined solution (second strategy) does not necessarily lead to the same result.

Numerical shooting method
The numerical shooting method approximates the space of initial conditions spanned by y 0 , y 1 by a lattice with a suitable lattice spacing. For illustrative purposes we use ∆y i = 0.1, corresponding to 4941 pairs of initial conditions on the (y 0 , y 1 )-plane shown in Fig. 1. At these lattice points, the ODE (2.1) is then solved numerically on the interval x ∈ [0, x sing −ǫ] employing a standard ODE-solver. Practically, we chose ǫ = 10 −3 . The numerical solutions provide a map (y 0 , y 1 ) → e(x; y 0 , y 1 )| x=x sing −ǫ . (2.11) The regularity condition (2.2) corresponds to e(x sing − ǫ, y 0 , y 1 ) = 0 and defines a submersion on the space of initial conditions (y 0 , y 1 ). Technically, this condition is implemented by looking for sign changes in e(y 0 , y 1 ) along the strips of constant y 0 . The resulting (interpolated) values for y 1 indicating a zero of (2.11) are given by the red squares displayed in the right panel of Fig. 1. Fitting a linear function to the data points yields This matches the exact result (2.5) with three digit precision. Thus all three methods, exact solutions, polynomial approximations and the numerical shooting method are suitable to quantitatively analyze the constraints implied by the existence of a global solution on the free parameters of the ODE.

The flow equation for f k (R)
In this section, we construct the PDE encoding the scale-dependence of f k (R) based on the guiding principles laid down in the introduction. The key results are eqs. (3.21) and (3.24) which constitute the starting points in our search for fixed functions in the subsequent sections.

Projecting the RG flow
The construction of the flow equation (1.1) employs the background field method, splitting the physical metric g µν into a fixed backgroundḡ µν and a fluctuation field h µν When projecting the RG flow onto the ansatz (1.2), it is convenient to chooseḡ µν as the metric on the compact four-sphere S 4 with background curvatureR and preform a transverse-traceless decomposition of the metric fluctuations with respect to this background The component fields consist of a transverse-traceless tensor h T µν (spin 2), a transverse vector ξ µ , and the two scalars σ, h (spin 0) satisfying the constraints Following [60], the projected flow for the ansatz (1.2) takes the form with the explicit expressions for the traces given by Here t = ln k/k 0 is the dimensionless "renormalization group time", the primes denote derivatives of f k (R) with respect toR and P k ≡ + R k ( ) denotes the regularized kinetic term constructed from the spin-dependent coarse-graining operator (s) = −D 2 + E (s) . In the nomenclature of [24], the inclusion of non-trivial endomorphisms E (s) in the regulator function gives rise to a type II regulator. On a spherical background E (s) is proportional toR and we set The traces T (2) and T  configuration space of the gravitational path integral [60]. For vanishing endomorphism the traces coincide with earlier constructions [22,23]. The next step consists in evaluating the traces (3.5). In the case where the flow equation is expanded in powers of the background curvatureR one typically resorts to the local (or early-time) expansion of the heat-kernel which is truncated at the desired order inR. The resulting expansion is a good approximation in the regime k 2 ≫R, but does not give rise to the correct IR behavior, k 2 ≪R, since, by definition, it does not include the non-analytic terms reflecting the compact nature of the background. In order to implement construction principle 3 to obtain a flow equation with the correct asymptotic for both k 2 ≫R and k 2 ≪R, we evaluate the functional traces as the sum over the eigenvalues of the differential operators. The spectrum λ l (d, s) and multiplicities M l (d, s) of the Laplacian on a spherical background are well known and summarized in Tab. 1. Introducing the functions W (i) the functional traces can be written as where ∆ (s) ≡ −D 2 denotes the standard Laplacian acting on fields with spin s.
In order to evaluate the infinite sums explicitly, it is convenient to switch to dimensionless quantities Moreover, we specify the regulator R k to be of Litim type [61] R (s) where θ is the Heaviside step function. Substituting the traces (3.7) into (3.4) and expressing the result in terms of dimensionless quantities then yields (3.10) Here the dots and primes denote derivatives with respect to t and r, respectively and we omitted the arguments of ϕ k (r) to ease readability. The coefficients in the transversetraceless sum arẽ while the coefficients in the scalar trace are given by (3.12) The explicit form of the step functions is given by Each eigenmode thereby gives a specific contribution to the trace whose weight is determined through the Hessian Γ (2) k . From the step functions (3.13) one explicitly sees that lowering k (which, for fixed background curvatureR, corresponds to increasing r) the eigenmodes are integrated out consecutively starting from fluctuations with large values l: once k 2 crosses the corresponding eigenvalues the step-functions ensure that the contribution of the associated spherical harmonic drops out from the rhs of the flow equation and does no longer contribute to the running of ϕ k .
The structure of the step functions (3.13) also clarifies the role of the endomorphisms E (s) . The parameters α and β allow to adjust the value k (s),term at which the last fluctuation mode of the corresponding trace is integrated out. Substituting the l-value of the lowest eigenmode and equating the argument of θ (s) to zero yields (3.14) For r > r (s),term the corresponding trace vanishes identically. The condition that both the spin-two and scalar trace are completely integrated out at the same value r term can be implemented by requiring that the coarse-graining operators (s) have equal lowest eigenvalues (ELE). Equating r (2),term and r (0),term shows that the ELE condition is satisfied along the line In the following, we will implement this condition thereby realizing the guiding principle 4. This choice is motivated on physical grounds, since it implies that all fluctuations are integrated out at one fixed value of k. It turns out that the ELE-condition can consistently be imposed with all other guiding principles.
In principle the endomorphisms can also be used to integrate out all fluctuation modes on a compact interval by selecting r (2),term = r (0),term < ∞. For r > r (0),term the rhs of the flow, encoding the quantum effects, becomes trivial and the flow equation reduces to the classical scaling relation. This was the strategy implemented in [54], which constructed fixed functions on the compact interval 0 ≤ r ≤ r (0),term = 6. While this compactification simplifies the numerical analysis, it will not be implemented here and we will consider fixed functions which are well-defined on the entire positive half-axis 0 ≤ r < ∞.

Smoothing the staircase
A feature that makes (3.10) particularly difficult to work with is that the PDE possesses an infinite number of discontinuities. At the points , l ≥ 0, (3.16) where k 2 crosses an eigenvalue of (s) , the step functions (3.13) remove the corresponding eigenmode and the rhs of (3.10) changes discontinuously. This feature is particular to the use of the Litim regulator and would be absent in the case of a smooth (e.g. exponential) regulator. In order to arrive at a smooth PDE we follow the strategy of [55] and approximate this staircase-behavior by a smooth function. This is achieved as follows. Firstly, we observe that, for r > 0, the sum over eigenvalues is truncated to a finite sum. In order for a term to contribute to this series, the corresponding argument of the theta-function has to be positive. This gives an r-dependent upper boundary N r on the number of terms in the finite sum 12 − r(12α + l(l + 3) − 2) > 0 ⇐⇒ l < r((17 − 48α)r + 48) − 3r 2r =: N (2) (r) (3.17) for the tensor modes and for the scalar modes. Again, it can be seen that only for r = 0 infinitely many terms are contributing. Regarding the l-dependence, the spectral sum contains only a polynomial in l, as can be seen in eq. (3.12). For r > 0 each monomial can be summed using the relation with Bernoulli numbers B n . Substituting M → N (i) (r) yields a smooth approximation bounding the staircase-function from above. A smooth approximation that bounds the sum from below is obtained by a sum from n = 1 to n = N r − 1. In practice, we approximate the staircase by the average of this two approximations, setting Applying the summation formula (3.20) to the rhs of the flow equation (3.10) yields the desired smooth PDE govering the scale-dependence of ϕ k (r) The coefficients c i andc i are polynomials in r and depend on the endomorphisms and (3.23) The RG flow of the function ϕ k (r) is governed by the PDE (3.21). Within this type of functional truncations, fixed points are promoted to fixed functions ϕ * (r) given by globally well-defined, k-stationary solutions of the PDE. Dropping all t-derivatives from (3.21) then yields a third-order ODE for the fixed function ϕ * (r) with coefficients given in (3.22) and (3.23). Since there is no risk of confusion we will omit the asterisk and denote the fixed function only by ϕ(r).
One particular effect of the smoothing procedure applied in this subsection is that rhs of the flow equation is non-zero on the entire interval r ∈ [0, ∞[ even if the endormorphisms are chosen in such a way that the quantum fluctuations in the initial flow equation (3.10) are integrated out on a finite r-interval. Thus, implementing the guiding principle 5, we insist that a "globally well-defined" solution is required to be regular on the full interval r ∈ [0, ∞[.

Conformally reduced flow equation
Based on the full flow equation (3.21), it is straightforward to write down the "conformally reduced" approximation where the RG flow is solely driven by fluctuations of the scalar mode. Eliminating the contributions from the transverse-traceless fluctuations by setting c 1 andc 2 to zero yieldṡ with coefficients (3.23). Setting all t-derivatives equal to zero yields the fixed point equation which is the analogue of the ODE determining the fixed functions in the full system (3.24).
Comparing the ODE's of the conformally reduced and full system one finds that they share all the characteristic features: both settings give rise to a non-linear ODE of third order. Moreover, the coefficient multiplying the highest derivative of ϕ is identical, indicating that both equations possess the same fixed singularities. On this basis, one may expect that both equations also share the same number of globally well-defined solutions. This expectation will be confirmed in Appendix A where eq. (3.26) will be analyzed in detail.

Analytical properties of the stationary flow equation
The ODE (3.21) encoding the fixed functions of QEG in the f (R)-approximation is a complicated non-linear third order differential equation making the construction of exact analytic solutions a difficult task. Nevertheless some features of the ODE can be accessed by analytic methods. We start with analyzing the fixed singularities in Sect. 4.1 before zooming into the asymptotic IR behavior of solutions in Sect. 4.2.

Singularity structure
Being of third order, the ODE (3.21) gives rise to a three-parameter family of solutions which exist at least locally. The condition that such a local solution can be extended to a global solution, existing for all r ≥ 0, gives rise to constraints on these parameters. If the ODE contains fixed singularities, global solutions have to fulfill regularity conditions at these points in order to be able to pass through the singularity. This constrains the space of solutions and generically one expects that each fixed singularity reduces the number of free parameters by one. The existence of isolated, globally defined solutions generally requires a total number of three singularities. This condition is met, for example, if the ODE has three fixed singularities in its domain of validity. In order to determine its fixed singularities, the ODE (3.21) is cast into normal form by solving for ϕ ′′′ (r). One then sees that the coefficient c 4 multiplying rϕ ′′′ (r) determines the fixed singularities of the ODE: zeros of c 4 constitute points where the ODE develops a pole. Due to the factor r multiplying ϕ ′′′ (r), one pole is always real and located at r = 0. This pole originates from switching to dimensionless quantities and the fact that R is not dimensionless. The coefficient c 4 is a forth order polynomial in r, giving rise to four additional roots. Thus set of fixed singularities is given by r sing,1 = 0 , r sing,2 = 9 5 + 9β , r sing,3 = 1 β , r sing,4,5 = 6 6β − 1 . (4.1) The position of the roots r sing,i , i = 2, 3, 4, 5 depends on the endomorphism β and is shown in Fig. 2. Notably, r sing,2 diverges at β = − 5 9 and the double pole r sing,4,5 is shifted to infinity if β = 1 6 , so that these are distinguished choices for the endomorphism. The analysis in Fig. 2 shows that that for β = 0 the ODE has two fixed singularities at finite, positive values r. By choosing the endomorphism in the interval 0 < β ≤ 1 6 the number of fixed singularities is increased to three. For β > 1 6 the number of fixed singularities further increases to four. Thus the inclusion of the endomorphism allows to control the singularity structure of the ODE to a certain extend.
Based on the singularity counting argument presented in Sect. 1.2, we then select β such that the number of fixed singularities matches the order of the ODE. In anticipation of the numerical analysis, we will impose the condition of ELE and choose In this case, the coefficient c 4 reduces to a third order polynomial which considerably simplifies the numerical analysis of the ODE computations.

Asymptotic IR behavior of the solutions
The second property which can be accessed analytically is the scaling of the solutions ϕ(r) in the IR limit r → ∞. In order to zoom into this asymptotic region, we rescale r and ϕ by a constant scale parameter ǫ, and require that the new variablesr andφ are strictly of the order one. Since a > 0 the limit ǫ → 0 corresponds to r → ∞. Retaining the leading terms in ǫ only, eq. (3.24) simplifies to whereφ =φ(r) and the primes denote derivatives with respect to the argument. The lhs is independent of the scale parameter a, which can be traced back to the fact that the lhs of eq. (3.24) is invariant under rescalings of the form r → λr, ϕ → ϕ. Similarly, the rhs is independent of the scale parameter b reflecting that the rhs of eq. (3.24), is invariant with respect to the rescaling ϕ → λϕ, r → r. Notably, the rhs of (4.4) receives contributions from the transverse-traceless sector only. The choice β = 1/6 reduces the order of polynomials c i so that the contribution of the scalar sector is sub-leading in the large-r limit. Depending on the ratio of a and b, eq. (4.4) exhibits three different scaling behaviors: classical scaling for b > 2a, quantum scaling for b < 2a and balanced scaling b = 2a. These cases will be discussed below.  which entails thatφ with A 1 and A 2 free integration constants. Notably, the existence of theφ(r) = A 1r 3 2 solution is not intrinsic to the special choice of endomorphisms used to derive (4.4) but appears for generic values α and β where the scalar sector also contributes to the asymptotic IR scaling.

c) balanced scaling b = 2a
In this case, the classical and quantum terms balance and eq. (4.4) entails that  For this equation to be satisfied,r has to have the same power on both sides implying classical scaling a = 2. For this choice the lhs vanishes, however, so that a = 2 does not constitute to a consistent solution.
The structure of (4.9) suggests modifying the scaling ansatz including a logarithmic termφ (r) = Ar 2 (1 − η log(r)) . Also, in this case the appearance of a logarithmic contribution in the scaling solution constitutes a generic feature and is not related to the particular choice of endomorphisms. The scaling solutions obtained in the cases a), b), and c) constitute the three admissible asymptotic behaviors of potential solutions of the full ODE in the IR. A priori it is unclear which phase is realized by a potential fixed function. Nevertheless, this scaling analysis provides important information on the IR dynamics of potential fixed functions in QEG.

Numerical construction of the fixed function
In this section, we use a numerical shooting method to construct solutions of the ODE (3.24) which are regular at the fixed singularities. The central result of this construction is the unique, isolated and globally well-defined fixed function shown in Fig. 4.

Implementing regularity at the fixed singularities
Our starting point is the ODE (3.24) for the specific choice of endomorphisms (4.2). The analysis (4.1) shows that in this case, the ODE possesses three fixed singularities located at finite values r: r 1,sing = 0 , r 2,sing = 18 13 , r 3,sing = 6 .
Regularity of the solution at these points imposes non-trivial boundary conditions on ϕ(r i ) and its derivatives. Casting the ODE into normal form by solving for ϕ ′′′ , and expanding the resulting rhs at one of the fixed singularities yields a Laurent expansion of the form with the residues e i given by Here it is implicitly understood that ϕ(r) and its derivatives appearing in e i are evaluated at r i,sing . Keeping ϕ ′′′ (r)| r=r i finite then requires e i = 0. A globally well-defined solution has to satisfy the resulting boundary condition (BC) at all three singular points (5.1). On this basis, we expect that the BCs fix all three free parameters characterizing the solutions of the ODE locally, so that the set of globally well-defined solutions is discrete. The next step consists in characterizing the solutions of the ODE by free parameters. Imposing that ϕ(r) satisfies the regularity conditions (5.3) ensures that the corresponding solution can be expanded in a Taylor series at r i,sing since all its derivatives remain regular. In particular a fixed function has an analytic expansion in the UV (r = 0), so that the solution locally takes the form Substituting this ansatz into the ODE and imposing regularity at r 1,sing by requiring e 1 = 0 uniquely fixes the constants g n , n ≥ 2 in terms of (g 0 , g 1 ). In particular, Thus, implementing the regularity condition e 1 , we obtain a two-parameter family of (local) solutions specified by (g 0 , g 1 ).
In a similar fashion, one can perform an analytic expansion of ϕ(r) at the other singular points Combining the ODE with the corresponding regularity condition again fixes the the coefficients b n and c n with n ≥ 2 in terms of (b 0 , b 1 ) and (c 0 , c 1 ), respectively. For the c n this solution is again unique, while the occurrence of (ϕ ′′ ) 2 in e 2 makes the solution b 2 (b 0 , b 1 ) double-valued. In the following we will use the parameterization (g 0 , g 1 ) to characterize potential global solutions, which then already satisfy regularity at r 1,sing = 0.

Extending the local solutions away from r = 0
The next step consists in extending the two-parameter family of solutions (5.4) to the interval r ∈ [0, ∞[ implementing the boundary conditions (5.3). In contrast to the example discussed in Sect. 2, we do not have an analytic solution of the ODE at our disposal. Moreover, the polynomial expansion is difficult to implement to sufficiently high order due to the complexity of the ODE. Therefore, we rely on the numerical shooting method illustrated at the end of Sect. 2. Since we need to implement the boundary conditions at r 2,sing and r 3,sing , we divide the positive half-axis into three intervals The shooting algorithm is applied twice: firstly, the solution is constructed on the interval I 1 . Solutions meeting the boundary condition e 2 = 0 are subsequently extended into the interval I 2 where the condition e 3 = 0 is tested. Finally, solutions satisfying all three boundary conditions are numerically extended onto I 3 . Owing to the non-linearity of the ODE, a solution might also terminate in a movable singularity so that it does not extend up to the fixed singularities marking the upper boundaries of the integration intervals. Practically, we implement this search strategy for fixed functions as follows. Because of the fixed singularity at r 1,sing = 0, the initial conditions for the numerical integration need to be imposed at r = ǫ 1 rather than at r = 0. The initial values for ϕ(r)| r=ǫ 1 and its derivatives in terms of the free parameters (g 0 , g 1 ) are obtained from the expansion (5.4), yielding We then discretize the (g 0 , g 1 )-plane by a lattice and use the initial conditions obtained at each lattice point to numerically integrate the ODE from ǫ 1 to r 2,sing − ǫ 2 . This defines a two-parameter family of solutions ϕ(r; g 0 , g 1 ) , r ∈ [ǫ 1 , r 2,sing − ǫ 2 ] . (5.10) Evaluating the BC e 2 , eq. (5.3), based on these solutions provides a map (g 0 , g 1 ) → e 2 (r 2 − ǫ 2 ; g 0 , g 1 ) . (5.11) The values for e 2 obtained in this way are shown in the top diagram of Fig. 3. Here lattice points with e 2 > 0 are marked green while points with e 2 < 0 are indicated in red. The regular points satisfying e 2 ≈ 0 are given by the black line. This set of regular points constitutes a one dimensional subspace that will be referred to as regular line. Every point on the regular line can be identified by g 0 , which is the only remaining free parameter. Thus implementing regularity of the solution at the fixed singularities r 1,sing and r 2,sing reduces the dimension of space of local solutions by two.
In the next step, we cross the fixed singularity by utilizing that the solutions regular at r = r 2,sing admit the polynomial expansion (5.6). The two coefficients b 0 and b 1 are determined by matching the expansion with the numerical solution (5.10)  Figure 3. Summary of the results obtained from the numerical shooting algorithm. Solutions satisfying the boundary condition e 2 , characterized by the two free parameters (g 0 , g 1 ) are shown in the top diagram. Points with e 2 (g 0 , g 1 ) > 0 are marked green while points with e 2 (g 0 , g 1 ) < 0 are marked red. The black line indicates the parameters which lead to solutions which are regular at r 1,sing and r 2,sing . Along this line g 1 is determined as a function of g 0 . The extension of these solutions to the interval I 2 is depicted in the bottom diagrams. Notably, most of the solution do not extend to r 3,sing but terminate in a movable singularity. The point of termination r max of the numerical solution as a function of g 0 is shown in the left panel. For the solutions extending to r 3,sing , the regularity condition e 3 (g 0 (g 1 )) is evaluated in the right panel. The minimum shown by e 3 (g 0 ) indicates that there is one isolated solution which satisfies all three BCs, e i = 0.
Here we explicitly indicated that the numerical solution depends on the single remaining free parameter g 0 , stressing that we consider the regular line only. Given the map we evaluate the polynomial expansion at r 2,sing + ǫ 2 in order to obtain the initial conditions for the numerical extension of the regular solutions to the interval I 2 ϕ(r 2 + ǫ 2 ) = ϕ 2 (r 2 + ǫ 2 ; g 0 ) , ϕ ′ (r 2 + ǫ 2 ) = ϕ ′ 2 (r 2 + ǫ 2 ; g 0 ) , ϕ ′′ (r 2 + ǫ 2 ) = ϕ ′′ 2 (r 2 + ǫ 2 ; g 0 ) . (5.14) The results obtained from the second numerical integration are shown in the bottom diagrams of Fig. 3. Notably, only a small window of points on the regular line, 0.029 The values of e 3 (r 3 − ǫ 3 ; g 0 ) are shown in the bottom-right diagram of Fig. 3. The map has a pronounced minimum at g 0 ≈ 0.0299. At this minimum the BC e 3 ≈ 0 is satisfied within the numerical accuracy of the shooting algorithm. This indicates that there is a single, isolated solution which is also regular at r 3,sing . Finally, the isolated solution is extended to the interval I 3 . Following the strategy employed when crossing the singularity at r 2,sing , the numerical solution is matched to the polynomial expansion (5.7) This determines the values of c 0 and c 1 for the regular solution. Following the strategy implemented in eq. (5.14), the parameters c 0 , c 1 are converted into initial conditions for a numerical integration on I 3 . The solution ϕ * (r) is then completed by carrying out this numerical integration and matching to the scaling regimes discussed in Sect. 4.2. 2 The numerical values of all expansion coefficients for the regular solution are summarized in the first line of Tab. 2. The resulting isolated and globally well-defined fixed function ϕ * (r) is shown in Fig. 4 and its properties are discussed further in Sect. 5.3. This solution is the first proper fixed function obtained within four-dimensional Quantum Einstein Gravity and constitutes the central result of this work. We close this subsection with some technical comments on the implementation of the shooting algorithm. Notably, the position and width of the plateau displayed in the bottomleft diagram of Fig. 3 has a mild dependence on the parameters entering the numerical integration. It is conceivable that it reduces to a single point (spike) when implementing a higher numerical precision. Owing to the fact that we are using a polynomial expansion at the singular points to continue the solutions through a fixed singularity the entire construction is very robust under variations of the auxiliary parameters ǫ 1 , ǫ 2 and ǫ 3 . Changes in ǫ i induce small changes (within error-bars) of g 0 that propagate into the coefficients b 0 , b 1 and c 0 , c 1 . These quantitative changes do not affect the qualitative picture shown in Fig.  3, however.

Properties of the complete fixed function
The proper fixed function ϕ * (r) satisfying all BCs (5.3) is shown in Fig. 4. Here, we will summarize its main properties.
Asymptotic expansion in the UV At r = 0, which corresponds to k → ∞ if the background curvature is held fixed, the fixed function admits the analytic expansion (5.4). The first two expansion coefficients are given in the first line of Tab. 2. In order to facilitate the comparison with results obtained in finite-dimensional approximations, it is convenient to compute the corresponding values of the dimensionless Newtons constant g * ≡ −(16πg 1 ) −1 and cosmological constant λ * ≡ −g 0 /(2g 1 ) together with the universal (gauge-independent) product τ * = g * λ * λ * = 0.411 , g * = 0.547 , τ * = 0.224 . (5.17) These values are in good agreement with the results obtained from expansions of flow equations truncated at finite order R N [8,16,17,65]. In particular τ * fits well to results obtained within "geometric" flow equations [31,60]. We take this agreement as strong evidence that, despite the different boundary conditions imposed on the flow equations, the polynomial expansions and the fixed function ϕ * (r) actually describe the same universality class.

Consequences of the global minimum
The fixed function possesses a global minimum located at r min ≈ 2.500. At this minimum ϕ * (r min ) < 0. This feature together with the resulting non-monotonicity distinguishes our solution from earlier constructions [55,56], where ϕ * (r) turned out to be positive definite and monotonically increasing. As a consequence our solution actually passes the redundancy test [58], requiring that has a zero for some value r. From Fig. 4 we then deduce Since E 4 (r) is continuous this implies that E 4 (r) has a zero within the interval r ∈ [0, r min ]. This provides a strong indication that linear perturbations of the fixed function cannot be absorbed by a field redefinition and constitute genuine eigenperturbations of the fixed function.
Asymptotic behavior in the IR Finally, it is interesting to investigate the asymptotic behavior of ϕ * (r) for r → ∞. In order to determine the scaling regime, we first approximate the large-r behavior through the ansatz  At this stage it is also illuminating to substitute the ansatz (5.20) into the ODE (3.24) and determine the coefficients A, u 1 and u 2 from a large-r expansion. Surprisingly, the ansatz provides an asymptotic solution if and only if the endomorphisms satisfy This condition is not satisfied for the choice (4.2) and becomes compatible with the ELE condition (3.15) for two specific values β ± = (18 ± √ 285)/78 only. Therefore we conclude that the asymptotic behavior of ϕ * (r) is described by the scaling solution obtained in the balanced regime and involves also non-analytic terms not present in the ansatz (5.20). This is actually in agreement with earlier works [56,58,66] which also observed the presence of non-analytic terms in this asymptotic scaling regime.

Conclusions and discussion
The existence of a NGFP on the gravitational theory space is a crucial ingredient for constructing a consistent and predictive quantum theory of gravity along the lines of Weinberg's Asymptotic Safety scenario. In this work, we used the functional renormalization group equation for the effective average action Γ k , (1.1), to establish the existence of a suitable fixed function in the realm of f (R)-gravity in four spacetime dimensions. Despite the vast body of earlier works [55,56,58,67], this is the first time that a suitable fixed function has been found and its construction constitutes an important step towards generalizing finite-dimensional truncations of the gravitational RG flow to approximations containing an infinite number of scale-dependent coupling constants. Our solution is globally well-defined on the positive real axis and its UV-expansion at r = 0 yields a positive cosmological constant and Newton?s constant. The structure and properties of the fixed function is in good agreement with the fixed points solutions found at the level of polynomial f (R)-truncations [60]. This construction constitutes an important step towards generalizing finite-dimensional truncations of the gravitational RG flow to approximations containing an infinite number of scale-dependent coupling constants.
The fixed function is based on non-linear ODE of third order, which construction principles discussed in the introduction. Owing to the inclusion of non-trivial endomorphisms in the regulator the ODE possesses three fixed singularities at finite r ≡ R/k 2 . This feature restricts the space of globally well-defined solutions to one isolated fixed function. The condition of "equal lowest eigenvalues" thereby plays an important role for the existence of this solution: numerical analysis indicates that once the transverse-traceless and scalar sectors of the flow equation are out of tune, the numerical solutions typically terminate in a movable singularity. Combining a scaling analysis of the ODE together with the numerical solution shows that f * (R) = AR 2 (1 + η log R + . . . ), matching the IR behavior of earlier constructions [55,56,58,67].
In Appendix A we analyzed conformally reduced version of the full flow equation. The resulting ODE shares many properties with the full flow. In particular it is of the same order and exhibits the same fixed singularities as the full system. At the same time the conformally reduced system is significantly simpler which allows investigating the existence of fixed functions by employing complementary techniques. As the central result we established that the conformally reduced equation also exhibits one globally well-defined isolated fixed function whose properties are qualitatively similar to the ones found in the full system (cf. Figs. 4 and 7). The result establishes that conformally reduced approximation is a very useful setting when searching for fixed functions of the gravitational flow equation.
At this stage, we have not yet analyzed the structure of admissible, globally welldefined, infinitesimal deformations of the fixed function. Substituting the corresponding ansatz describing these deformations into the PDE (3.21) leads to a linear third-order ODE which also contains the (possibly complex) scaling exponent of the deformation. In principle the corresponding free parameters can again be determined by applying the combination of analytical and numerical techniques employed in the construction of the fixed function. This numerically very demanding, challenge is, however, beyond the scope of the present work. The match between the fixed function and its finite-dimensional sibling let us expect that the full functional also shares the stability properties of the polynomial approximation. Presupposing that all relevant eigenvalues can be extended to global deformations, this implies that there should be three relevant deformations. This expectation that there is a finite number of relevant directions is also supported by the general arguments [57]. As an important first check in this direction, we explicitly verified that our solution passes the redundancy test [58], indicating that the fixed functional might indeed have non-trivial deformations which cannot be absorbed by a field redefinition.
Let us again stress that carrying out the extension of the fixed points seen in finitedimensional truncations to the level of fixed functions is an important step in the Asymptotic Safety program: the functional approximations are much more sensitive to structural details of the theory, as, e.g., the functional measure used to construct the flow equations. These details have a crucial influence on the singularity structure of the ODE describing the fixed function and may provide important insights on conceptual properties required to achieve Asymptotic Safety. In this context, we find it encouraging that the geometrically motivated flow equation actually admits an interesting fixed function solution.
Notably, our findings are not in tension with the recent results [59] where no fixed function for four-dimensionally conformally reduced gravity has been found. The flow equation and search strategy implemented in this work differs from our setting in various ways, which can be made explicit based on the construction principles detailed in the introduction. Firstly, the flow equation constructed in [59] utilizes a conformally flat background while our equation has been derived for a compact background with positive curvature. Secondly, the analysis in [59] implemented the modified split ward identities of the effective average action to reduce the fixed point equation to a linear PDE. While this certainly simplifies the corresponding fixed function analysis the special choice of regulator or the implementation of split symmetry restoration at the level of the truncation may be overconstraining the system such that no globally well-defined solution may be found.
In this light, it would be very interesting to investigate to which extend the construction principles of our proposal are strictly necessary for the existence of a satisfactory fixed function. In particular it would be interesting to investigate if a flow equation based on the exponential parametrization of the metric fluctuations g µν =ḡ µρ e h ρ ν , (6.1) recently advocated in [68][69][70] also gives rise to fixed functions at the level of the f (R)approximation. Compared to previous constructions, this parameterization corresponds to a change in the path integral measure. Combined with a clever choice of gauge-fixing this may lead to a second order ODE describing the fixed function. This structure is much simpler to analyze since fixed functions may be identified based on a single shooting instead of the double-shooting method required in this work. In [69,70] it was already observed that this choice significantly improves the IR behavior of the gravitational RG flow.
(3.26) and the the fixed singularities are again given by (5.1). The boundary conditions that fixed functions have to satisfy are determined by Laurent expanding according to (5.2) and explicitly read Any regular solution satisfying the BCs (A.1) can be expanded in a Taylor series at r i,sing according to (5.4), (5.6) and (5.7). Substituting these expansions into (3.26) allows for determining all coefficients g n , b n and c n with n ≥ 2 in terms of the first two coefficients (g 0 , g 1 ), (b 0 , b 1 ) and (c 0 , c 1 ), respectively.
In order to construct globally well-defined solutions of (3.26) we repeat the numerical shooting algorithm used in Sect. 5.2. Beginning at r = 0 we have a two-parameter family of solutions characterized by (g 0 , g 1 ). To extend these solutions to [0, ∞[, we again split the positive half-axis into the intervals (5.8). The first shooting on I 1 is realized by numerically integrating from r = ǫ 1 to r 2,sing − ǫ 2 using the initial values (5.9). After discretizing the (g 0 , g 1 )-plane we perform for every lattice point the integration and evaluate the BC e 2 with it. This provides a map (g 0 , g 1 ) → e 2 (r 2 − ǫ 2 ; g 0 , g 1 ) . (A. 2) The resulting values e 2 are visualized in Fig. 8, where lattice points with e 2 > 0 are marked green while points with e 2 < 0 are tagged red. Regular points have to satisfy e 2 ≈ 0. Clearly, there is a line of points separating two regions of different signs (i.e. different colors in the diagram) giving rise to a line of regular solutions. The points on this line are again parametrized by g 0 . Next, we repeat the shooting procedure on I 2 . Utilizing the polynomial expansion at r 2,sing , initial values for the numerical integration are again determined by matching according to (5.12). Integrating from r 2,sing + ǫ 2 to r 3,sing − ǫ 3 again reveals a window of g 0 -values where the corresponding solutions extend up to r 3,sing . This regular window is qualitatively similar to the one displayed in Fig. 3 and located at 0.052 g 0 0.054. The unique value for g 0 minimizing the last BC e 3 is listed in Tab. 2 and the corresponding globally well-defined fixed function ϕ * (r) is displayed in Fig. 7. It shares all features of the result obtained from the full flow (3.24), including a minimum located at r min ≈ 2.153. This result underlines the importance of the conformal mode for the existence of fixed functions.

A.2 Fixed functions from bootstrapping method
We close this section by presenting an alternative bootstrapping strategy relying on polynomial expansions only [71]. In the conformally reduced setup all expansion coefficients g n , b n and c n are much smaller in size enabling us to go to much higher orders in the Taylor  Figure 6. Solutions regular at r = 0 are parametrized by g 0 and g 1 . Evaluating these solutions at the BC e 2 , eq. (A.1), yields the data points depicted above. Points with e 2 > 0 are marked green while points with e 2 < 0 are tagged red. There is a line of points separating two regions of different signs (i.e. different colors) giving rise to a line of regular solutions satisfying e 2 ≈ 0. The black line shows g 1 as a function of g 0 obtained by imposing the bootstrap condition g N +1 (g 0 , g 1 ) = 0 for N = 17 and solving for g 1 . The two approaches show very good agreement.
expansions, especially at r = 0. Thus, we start by truncating the Taylor series (5.4) to a finite polynomial of the order N ϕ 1 (r; g 0 , g 1 ) = g 0 + g 1 r + N n=2 g n (g 0 , g 1 ) r n .

(A.3)
If the Taylor series (5.4) has a finite radius of convergence R 1 and if additionally r 2,sing ≤ R 1 , then (A.3) can be substituted into the BC e 2 . Again, this provides a map of the form (A.2). The condition e 2 (g 0 , g 1 ) = 0 relates g 0 and g 1 analytically. Unfortunately, this relation grows quickly in complexity with increasing order N and thus, is of no real practical use. However, if the finite Taylor polynomial (A.3) serves a good approximation at the order N , the next coefficient should by very small, i.e. g N +1 (g 0 , g 1 ) ≈ 0. Thus, a much simpler condition that can be imposed is given by which constitutes the central idea of the bootstrapping strategy (or "minimal procedure" in the terminology of [71]). Being a rational function in g 0 and g 1 , at every order new solutions of (A.4) emerge, which are mostly complex or singular and thus, discarded as spurious.
Pushing the order up to N = 17, only one real and non-singular solution of g 18 = 0 has been found and is displayed in Fig. 6 (black line). Most remarkably, the analytic solution g 0 → g 1 is in perfect agreement with the regular line obtained from the shooting algorithm. where M is the order to which the expansion (5.6) has been truncated. Imposing (A.5) yields a line b 0 → b 1 in the (b 0 , b 1 )-plane (spurious solutions are discarded), c.f. Fig. 8 (black line). In contrast to the expansion at r = 0, no convergence of the line can be observed revealing the asymptotic nature of the expansion at r 2,sing . Thus, we keep the order M low (M = 3 at most) and "transport" the previously obtained regular line from r = 0 to r 2,sing into the (b 0 , b 1 )-plane. This is done by substituting the analytic solution g 0 → g 1 into the Taylor polynomial (A.3) and simultaneously setting r = r 2,sing . The transported line is defined as g 0 → (ϕ 1 (r 2,sing ; g 0 , g 1 (g 0 ), ϕ ′ 1 (r 2,sing ; g 0 , g 1 (g 0 )) ≡ (b 0 (g 0 ), b 1 (g 0 )) , (A. 6) and shown (blue line) in Fig. 8. Inspecting Fig. 8 shows Comparing with the values from the shooting method in Tab. 2 shows a good qualitative agreement. In summary, the bootstrap strategy provides a further argument in favor of the existence of a globally well-defined fixed function. However, due to the asymptotic nature of the expansion (5.6) a precise determination of the coefficients characterizing the fixed function is impossible. Thus the shooting method employed in the main part of this work constitutes the most reliable method for constructing fixed functions in this setting.