Brans-Dicke theory in the local potential approximation

We study the Brans-Dicke theory with arbitrary potential within a functional renormalization group framework. Motivated by the asymptotic safety scenario of quantum gravity and by the well-known relation between f(R) gravity and Brans-Dicke theory at the classical level, we concentrate our analysis on the fixed-point equation for the potential in four dimensions and with Brans-Dicke parameter omega equal to zero. For two different choices of gauge, we study the resulting equations by examining both local and global properties of the solutions, by means of analytical and numerical methods. As a result of our analysis we do not find any nontrivial fixed point in one gauge, but we find a continuum of fixed points in the other one. We interpret such inconsistency as a result of the restriction to omega equal to zero, and thus we suggest that it indicates a failure of the equivalence between f(R) gravity and Brans-Dicke theory at the quantum level.


Introduction
Many models of modified gravity have been proposed and studied over time in an attempt to address the problems encountered in quantum gravity and cosmology [1]. New models are most commonly postulated as the starting point of a new research direction, however, one instance in which the logic is partially reversed, and a new model of gravity is hoped to emerge as the final result, is the asymptotic safety scenario [2,3,4,5,6]. The general idea behind such scenario is that the theory of quantum gravity should be sought within a large class of theories (e.g. all possible theories described by an action functional of a single metric field), out of which one single theory (or few isolated ones) should emerge with the peculiar characteristic of being a fixed-point of the renormalization group (RG) flow. IR-unstable trajectories emanating from such fixed point(s) would then define nonperturbatively renormalized theories of gravity. The use of functional renormalization group equations (FRGEs) has provided considerable evidence in support of the existence of a nontrivial fixed point theory for gravity, in a large number of formulations and approximations (see list of references in [5,6]).
Even within a given class of models, specified by a choice of variables and symmetries, it is obviously impossible to explore the entire space of theories and one has to resort to approximations. A common approximation in the asymptotic safety literature consists in truncating the theory space to a finite-dimensional subspace by making a guess of what might be the most important terms to keep track of in the effective action. The guess is then supposed to be repeatedly refined until little improvement of the results is gained by new refinements. In practice, even this is quite hard, and only recently such program has been implemented to a high order of refinement for truncations that only retain polynomials of the Ricci scalar [7,8,9,10,11]. At the same time, the functional nature of the renormalization group methods being used has just started being exploited further, opening up the possibility of exploring infinite-dimensional subspaces of the theory space. The main idea is to study the RG flow of gravity in a spirit similar to the derivative expansion of scalar field theory [12,13,14,15]. There, the leading order approximation is known as local potential approximation (LPA), and it consists in projecting the flow equation on a constant scalar field, thus allowing to study the flow of a generic effective potential V k (φ) = Γ k [φ = const.]/ d d x. The next to leading order includes a term with two derivatives and any field dependence, and so on at higher orders. At each order of the derivative expansion one is lead to study partial differential equations for functions of the field φ and of the running scale k. In gravity there is much more structure, and there are probably many options in organizing an expansion of this sort. A very natural option is to organize the action as if it was an expansion around a maximally symmetric background. For the latter the only non-zero component of the Riemann tensor is the Ricci scalar R, which is constant: we have ∇ µ R = S µν = C µνρσ = 0, where S µν is the traceless Ricci tensor, and C µνρσ is the Weyl tensor. The analogue of the derivative expansion can then be an expansion in S µν , C µνρσ and their derivatives (by the Bianchi identity ∇ µ R = 2d d−2 ∇ ν S µν ), with arbitrary dependence on R at each order. In the leading order of such an expansion we are left with an f (R) theory, whose study in such spirit was begun in [16,17,18,19,20].
As compared to the LPA for scalar field theory, in the f (R) approximation for gravity we face a number of additional technical complications, in particular a larger number of contributions to the FRGE, with a more complicated dependence on the unknown function, and the challenge of evaluating functional traces on a curved background. The latter in particular introduces some subtleties related to the presence of zero modes in compact backgrounds and to the staircase nature of the results obtained for the traces when using cutoffs with step functions [16]. Also for these reasons progress has been slow in this direction, and it is desirable to find alternative ways to study the same problem.
One possible simplification, which we will explore in this paper, is suggested by a well known classical property of the f (R) theory [21,22]. The classical action for f (R) gravity, 1 can be traded for an equivalent action, describing a scalar-tensor theory, The relation between the two Lagrangians is given by a Legendre transform, where R(φ) is found by solving the equation φ = −f ′ (R) for R, and as usual, regularity of the transform is guaranteed if f ′′ (R) = 0. From a RG perspective the advantage of such formulation is that we can study the running of the potential by projecting the FRGE on a flat background, thus sidestepping all the complications of curved backgrounds. In fact, we will construct flow equations for a more general version of (1.2), that is, a generic Brans-Dicke theory with a potential (see (2.1)). 2 Projection on a flat background will allow us to study such theory without truncating the potential to a polynomial form, thus performing an analysis similar to that of pure scalar theory [28,29,30,31,32,33].
Of course, at a quantum level the two theories might be inequivalent. They are both perturbatively nonrenormalizable, and standard perturbative reasonings could only apply at an effective field theory level. When looking for a UV completion in the form of a nontrivial fixed point, we study the RG equations in two different theory spaces, and in the full fixed point theory the scalar field might couple to other geometric invariants, or acquire its own dynamical term. Moreover, since the functions f (R) or V (φ) are not chosen a priori, but have to be found such that they correspond to an RG fixed point, it could happen that the regularity of the transform fails at one or several field values (or even in the full range of definition if for example f (R) is linear at the fixed point). As a consequence, if fixed points were to be found in both formulations, they might describe different physics. It might also happen that one formulation admits an asymptotic safety scenario and the other does not. 3 However, we also do not know a priori whether the (nonper-1 As in most asymptotic safety research (for a rare exception see [23]), we will be working in Euclidean signature. 2 Note that in the context of asymptotic safety, Brans-Dicke theory with ω = 0 was considered in [24] as a RG improvement of the Einstein-Hilbert truncation, in which the running gravitational and cosmological constants were promoted to fields as a result of an identification of scales with spacetime points. Clearly our work differs substantially from [24], as we study the RG equations directly for the Brans-Dicke theory. In a sense our work relates to [24] like the general f (R) studies [16,17,18,19,20] relate to the f (R) actions obtained by improvement of the Einstein-Hilbert truncation [25,26,27]. 3 In addition, we should also notice that often in the cosmology literature other "frames" are considered, in which a new metric field is defined via a conformal map, often together with a redefinition of the scalar field as well. Again, at the classical level these are all equivalent theories (although there has been some confusion on the issue in the past [34]), but at the nonperturbative quantum level this is not guaranteed (perturbatively there is on-shell equivalence, as shown for example in two dimensions [35], but at the nonperturbative level this has not been shown, although it might be possible following the developments of [36,37]). We will not study here other frames, having always in mind the original pure metric theory, whose metric we assume to define the coupling to ordinary matter.
turbatively renormalized) quantum theories are equivalent or not, and only a direct comparison (which we can at least do at the level of truncations) could allow us to settle the question. In any case, given that in asymptotic safety we are in principle allowing for extra degrees of freedom, there seems to be no reason to consider only pure metric theories of gravity, and the study of scalar-tensor theories is of interest in its own. Brans-Dicke theory is one of the oldest modifications of general relativity [38], and together with its variations and generalizations it finds plenty of applications in cosmology [1], and in quantum gravity (e.g. [39,40,41]). Other versions of scalar-tensor theories have already been studied in the context of asymptotic safety [42,43,44], but to the best of our knowledge, no study of this sort has been done before for the Brans-Dicke theory in the formulation we consider here (sometimes referred to as Helmholtz-Jordan frame). We will introduce more precisely our ansatz and setup in Sec. 2, together with the two choices of gauge fixing that we are going to employ. In Sec. 3 we will derive the FRGEs for both gauges for general dimension and Brans-Dicke parameter, while in Sec. 4 we will proceed to analyze their local properties in d = 4 and ω = 0. Finally, we will present the results of numerical integrations in Sec. 5, concluding in Sec. 6 with a discussion of our findings and of future prospects. In App. A we analytically solve the special case d = 2, which helps illustrating some of the points discussed in the conclusions.

The general setup: ansatz, variations and gauge fixing
The action (1.2) is a particular case (ω = 0) of the more general Brans-Dicke theory with a potential,Γ which in turn is a special case of dilaton gravity. The f (R) theory in the Palatini formalism is related to the same theory but with ω = −3/2 [45]. We have introduced a subscript k which stands for the running scale at which the effective average actionΓ k is defined [14]. To a large extent we will keep ω general, only to concentrate on the specific case ω = 0 for our numerical analysis (studying the running of ω would require using a non-constant background, or looking at the 2-point function, which we leave for future work). Note that (2.1) differs from other scalartensor theories studied in the asymptotic safety literature [42,43,44] in two important aspects: it is not invariant under φ → −φ (and of course φ is not restricted to be positive), and the kinetic term (when present, that is, when ω = 0) contains an inverse of the field. The point of view we wish to adopt in this paper is that the action (2.1) is a LPA approximation for the effective action, and that only to next order we would promote ω and the function coupled to R to general functions of φ. As we explained in the introduction, we will project the flow equation for (2.1) on a flat background and study only the running of the potential. However, for future reference, we will present in this section the results of variations and gauge fixing for a general maximally symmetric background metric and constant background scalar field.
Introducing the background splitting we make the usual approximation for the effective average action [46] and neglect the running of the gauge-fixing and ghost actions, S gf and S gh . For the FRGE we will need the second variation of the effective average action, therefore we expandΓ and find (omitting from now on the field dependencies of the action functionals) (2.5) We can exploit the gauge-fixing freedom to simplify the Hessian operator, adding to the original action the gauge-fixing term for some choice of gauge-fixing constraint F µ and of non-degenerate operator G µν . Physical results should be independent of the gauge choice, however, it is well known that the off-shell effective action is not gauge independent, and furthermore, the approximations we employ in the FRGE lead to additional gauge dependences. It is then important to test our analysis against different choices of gauge. We present in the following the two types of gauge which we will use in the forthcoming sections.

Feynman gauge
First we consider a Feynman gauge (α = 1) with and G (F )µν = φ g µν . (2.8) The total quadratic action becomes (2.9) Decomposing h µν =ĥ µν + 1 d g µν h, with g µνĥ µν = 0, we finally obtain (2.10) We note that via the gauge-fixing procedure we have introduced a kinetic term for the auxiliary field ϕ even in the case ω = 0. The kinetic term disappears for ω = −1/2, which is a special value for the Brans-Dicke theory in this gauge.
For the gauge sector we employ a standard Fadeev-Popov determinant which we rewrite in terms of a quadratic integral over complex Grassmann fields C µ andC µ . For constant background scalar field, the ghost action reads (2.11)

Landau gauge
As an alternative choice of gauge, we consider a Landau gauge (α = 0) with and The interesting aspect of such gauge is that it does not modify the kinetic term of ϕ, in particular it does not introduce one for ω = 0.
In this case, in order to simplify the non-minimal operators that appear in the second variation, we use the transverse-traceless decomposition of the metric fluctuations, given by with the component fields satisfying In the α → 0 limit, the ξ µ and σ field components decouple completely from the rest of the Hessian, and their contribution to the FRGE cancels exactly with the ghost contribution, when properly implemented [47]. We thus write the second variation of the action directly omitting the contribution of the longitudinal components: (2.16) Because of the change of variables (2.14), in this case there is also a Jacobian to keep track of, which we do by introducing auxiliary fields. The Jacobian for the gravitational sector leads to the auxiliary action [47] where the χ T µ and χ are complex Grassmann fields, while ζ T µ and ζ are real bosonic fields. The Jacobian for the transverse decomposition of the ghost action is given by with ψ a real scalar field.

The flow equation
The flow equation can be evaluated by means of the Functional Renormalization Group Equation (FRGE), which takes the generic form [12,13,14,15] being Γ (2) k (x, y) = and where Φ is a superfield collecting all the fields involved in the quantum action, i.e. Φ ≡ {ϕ, h µν , · · · }. R k is a generic cutoff operator, t ≡ log(k) is the RG running scale and STr identifies a functional supertrace, carrying a factor 2 for complex fields and a factor −1 for Grassmann fields.
We will construct the cutoff operator in such a way to implement the substitution rule being r(z) a dimensionless smearing function. That is, we choose a cutoff of the form R k = Γ k . A convenient choice of smearing function, leading to a considerable simplification of the functional traces, and which we will therefore use here, is the so-called "optimized" cutoff [48] which reads where Θ(x) is a Heaviside step function.

Feynman gauge
The Hessian of the effective action is mostly diagonal in field space, with the only exception of the {h, ϕ} sector, thus the supertrace in (3.1) can be easily decomposed into standard functional traces. In the Feynman gauge we obtain where H k is the modified inverse propagator, namely H k = Γ (2) k + R k . The evaluation of the first trace requires to invert the h-ϕ matrix, which is trivial since the matrix elements commute. The ghost term takes a factor of minus two with respect to the other terms, because of the complex Grassmannian nature of the ghost fields.
The trace over a generic Riemannian manifold can be evaluated by means of a heat kernel expansion, but since we are interested in projecting the flow equation on a flat background we can evaluate the trace over modes as a simple integral over momenta. The derivative of the cutoff operator with respect to the RG time returns which reduces to the sole Heaviside step function using the property that the distributional product of the delta function with its argument is zero. Because of the step function, the trace reduces to a momentum integral between 0 and k, thus automatically rendering the functional traces UV finite, a well-known feature of the FRGE. Performing the trace we obtain The trace over the tensor structure gives the factor d(d + 1)/2 − 1 for the h T µν contribution and a factor d for the ghosts, counting the number of their independent components. Since we are working on a flat manifold and constant background field both the Ricci scalar and the kinetic operator vanish, so that equation (3.5) reduces to an RG flow equation for the dimensionful potential. We cast the equation in an autonomous form, i.e. with no explicit dependence on k, by introducing the dimensionless quantities is the classical part of the equation, which is linear in the potential, and 14) is the quantum part, which contains all the loop contributions, and which is responsible for the nonlinear character of the equation.

Landau gauge
Working in the Landau gauge the supertrace in (3.1) reads 15) where the contributions of ghosts and longitudinal modes have been omitted, since they exactly cancel each other as explained before.
After performing the integral over momenta we obtain The RG flow equation for the dimensionless potential in such a gauge reads then where the classical part T tree is the same as in (3.12), and the quantum part reads

Analytical study of the equation
We want to search for fixed point solutions of equation (3.11) and (3.19), i.e. search for scale invariant solutions such that ∂ tṼk = 0, requiring them to be globally analytic [28,49,31]. The latter requirement has a well-understood physical and mathematical justification, being a necessary condition for the existence of the average effective action at all values of k, and hence of the full effective action in the limit k → 0 (which in d > 2 requires the existence of the solution forφ → ±∞, see (3.10)). In addition, the condition of global analyticity is expected to reduce the continuous set of solutions to a discrete subset of acceptable ones. For ∂ tṼk = 0, both (3.11) and (3.19) reduce to second order ordinary differential equations, thus we expect 2-parameter families of local solutions, parametrized by the initial value conditions. Extending such local solutions to global ones, we generally have to impose constraints coming from the analyticity requirement and from the symmetries of the problem. In our case we do not have any constraints originating from symmetries (e.g. we have noφ → −φ symmetry, hencẽ V ′ (0) = 0 in general), and we will have to study the equation on the full real line imposing asymptotic boundary conditions atφ ∼ ±∞. The latter, due to the non-linear nature of the equations, could contain less than two free parameters, implying that global solutions would also necessarily be parametrized by less than two degrees of freedom. Other explicit constraints can originate from fixed singularities of the equation, requiring analyticity conditions (e.g. [50]), and it is hoped that the equation does not have too many such fixed singularity, which would require an over constraining of the solutions [16,18].
We will apply the following strategy to select solutions: i) we look for singularities of the equations, either fixed or movable, and study the behavior of the solution in a neighborhood of the singularity, ii) we study the large field asymptotic solutions of the equation and count the degrees of freedom of each class of solutions, iii) we numerically look for global solutions satisfying all the constraints.
The study of the large field asymptotic solutions is important also for other two reasons, namely, the derivation of the full effective action at the fixed point [16], and the relation to the f (R) theory, as we will explain in the concluding section.
We will present most of the analysis for the case ω = 0, although occasionally we will refer to other values. As in the Landau gauge the ω = 0 value is a critical value, analogous to the ω = −1/2 value for the Feynman gauge, we will treat separately the two gauges, starting with the Feynman gauge. Most of our considerations apply to generic dimension d > 2, although we will most often specialize to d = 4. In Appendix A we will treat the special case d = 2.

Fixed singularities
In order to look for fixed singularities, we search for poles of the denominator of the scale invariant flow equation ∂ tṼk = 0 written in normal form, i.e.
where N and D are polynomial functions obtained from (3.11). For d > 2 the only zero we find is atφ = 0, while for d = 2 the equation reduces to a first order equation with no fixed singularities. To test the consequences of such singularity in d > 2 we impose analyticity, and study the equation in a Laurent expansion. Locally, imposing analyticity means requiring the existence of a Taylor expansion of the solution, in other words we make the ansatzṼ (φ) = n≥0 v nφ n , and after plugging it into the equation we expand the latter in a Laurent series centered at the origin. At leading order, the equation in the Feynman gauge reduces to which vanishes either restricting to ω = −1/2 (the analogous case in Landau gauge will be ω = 0, see 4.2.1), or fixing the potential in the origin to As a consequence for d > 2 and ω = −1/2 we have one constraint, thus reducing the number of degrees of freedom at the origin to one. For technical reasons, when integrating the equation numerically, we need to start from an arbitrary small value of the field ǫ. The boundary condition at ǫ can then be parametrized in terms of the derivative of the field in zerõ being τ =Ṽ ′ (0) the free parameter, and evaluated by means of a MacLaurin series and higher order coefficients are likewise obtained solving the equation order by order in ǫ.

Movable singularities
The constrained differential equation admits now a one parameter family of local solutions parametrized by τ . Still, because of the non linearity of the equation, we expect most of the solutions to end at a movable singularity, i.e. at a singularity whose location depends on the initial condition. We want to study the behavior of solutions in the neighborhood of such singularities, in order to confirm analytically the existence of such singularities and be able to recognize them in the numerical integrations, as well as to discuss possible interpretations in the terms of the f (R) theory. We will present in the next section the results of our search for a set of values of τ for which the singularity goes to infinity. Letφ c be the value of the field at which the singularity occurs, and suppose that the singular behavior is such that there exists an n 0 ≥ 0 such thatṼ (n) (φ c ) ∼ ∞ for every n ≥ n 0 . In order to understand what values of n 0 can occur for our equation, it is convenient to recast the equation (3.11) in the following form where the P i are two polynomials containing the same monomials but with different coefficients. As the polynomials P i have the same structure we deduce that forφ →φ c their ratio will in general go to a constant for any value n 0 . Special situations can arise when some cancellation occurs in P 2 which does not happen in P 1 , and such cases will have to be discussed separately. As a consequence, in the general case the linear part of the equation cannot diverge, otherwise it could not be balanced by the rational part, i.e. both the potential and its first derivative do not diverge at the singularity, restricting the possible value of n 0 to n 0 > 1. At this stage, we can assume that in the neighborhood ofφ c the potential can be written as and that γ > 1 (so that n 0 > 1), and we can try to determine the value of γ by means of the method of dominant balance. In order to do so we can start with the guess that the second derivative is divergent atφ c , that is 1 < γ < 2. In such case, by studying the balance of terms we arrive at the equation γ − 1 = −γ + 2, leading to γ = 3/2, in accordance with our guess. Plugging (4.7) with γ = 3/2 into (4.6), we can iteratively work out all the coefficients in the expansion as functions of the parameter u 0 and of the singular field valueφ c . For example, in d = 4 we find for the leading order terms. The subleading corrections can be computed iteratively, but their expression is very long, and not particularly enlightening.
Other singular behaviors are possible if P 2 has a zero. Such situations are more easily uncovered by studying the equation written in normal form, (4.1). Assuming that the first derivative of the potential is divergent (or more divergent than the potential itself) atφ ∼φ c , we obtain the equationṼ leading to a simple pole solutionṼ (φ) ∼ (φ −φ c ) −1 , which is consistent with the assumption. Subleading corrections can be worked out, confirming the possibility that such type of singular behavior can appear in a solution of the fixed point equation.

Behavior at large field values
We apply here the method of dominant balance to study the large field regime of the differential equation (3.11). We have already seen in (4.6) that whatever is the leading term (forφ → ∞ in this case) the quantum part of the equation in general goes as a constant plus subleading corrections, hence we have two possibilities: either the potential diverges at infinity, and the classical part of the equation defines the leading order, or the potential goes to a constant, and there must be some balance between linear and nonlinear part. In the first case, in theφ → ∞ limit the solution goes asṼ (φ) ∼ Aφ d d−2 + subleading terms , (4.10) where A is a free parameter. Subleading terms can be calculated by solving iteratively the differential equation for an ansatz of the typẽ (4.11) For d = 4, for example, the first few coefficients a n (A) are The coefficients are all inversely proportional to the bare parameter A, so that this expansion cannot be continued to A = 0, and that case must be treated separately. The asymptotic solution so far constructed defines a one-parameter family of solutions parametrized by the variable A, but as the equation is second order, we can ask if the asymptotic solutions have more degrees of freedom.
In order to answer such question, we follow [31,18] and perturb the flow equation in the neighborhood of the solution we just found, i.e. we introduce a perturbation to the potential, substitute it into (3.11), and expand to linear order in ǫ. ReplacingṼ (φ) with (4.11), and keeping only the leading terms in the coefficients of the linear operator acting on the perturbation, in d = 4 we obtain the linear equation which allows a solution which goes asymptotically like where B 1 and B 2 are two integration constants. Note that eq. (4.14) seems to reduce to a first order equation for ω = −1/2, but as we will see for the Landau gauge for ω = 0 (which is the analogue of the case ω = −1/2 in the Feynman gauge) for that critical value of ω we simply need to include the subleading correction of the coefficient of δṼ ′′ (φ). Whereas the power-law solution in (4.15) merely shifts A in (4.11), the exponential solution would seem to be a new degree of freedom. However, for positiveφ (and ω > −1/2, otherwise the role of positive and negativeφ are interchanged) it grows faster than the solution it is perturbing, contradicting our asymptotic analysis, hence it must be discarded. On the other hand, for negativẽ φ it is an exponentially small perturbation, hence it is acceptable. As the perturbation is smaller than any power at largeφ, while the leading solution (4.11) contains only powers, it is not difficult to see that the full equation decomposes in a hierarchy of equations, according to powers of the exponential correction, that is, the exponential acts like an ǫ parameter and we can iteratively solve the equation to obtainṼ whereṼ [0] (φ) is the leading solution (4.11), while for ω = 0 we find Z(φ) = 768 π 2 A 2φ3 + 4224 π 2 Aφ 2 + 64 24 A + 769 π 2 φ , (4.17) and so on, leaving A and B as free parameters. The presence of a new degree of freedom atφ ∼ −∞ creates an interesting situation, as we already know that we have an analyticity constraint atφ = 0, hence if we had just one-parameter families of solutions at both plus and minus infinity it would be unlikely to have a global solution. 4 There remains to consider the special case A = 0, which we now proceed to examine for d = 4 and ω = 0. From the previous discussion of dominant balance we would expect in such case a solution that asymptotes to constant. Nevertheless, we should be careful as in that analysis we have excluded special cases leading to cancellations in the denominator of the quantum part of the equation. By plugging into the equation an ansatz of the typẽ we find at leading order the equation A 1 = 0, in accordance with the previous analysis. However, a careful look at the higher orders of the expansion reveals the presence of poles at A 1 = 1 and A 1 = 3/2, meaning that for those values the general expansion is not valid, and a separate treatment is needed. In fact, we find that such special values of A 1 also lead to solutions that are solvable with an iterative algorithm. 5 In all three cases (A 1 = 0, 1 and 3/2) we find no free parameter in the expansion (4.19), but by studying the linear perturbations we discover the presence of exponentially small corrections at negativeφ for A 1 = 0, exponentially small corrections at both positive and negativeφ for A 1 = 1, and a non-integer power correction at negativeφ for A 1 = 3/2. It is quite easy to see that exponentially small corrections always carry one new degree of freedom, while the analysis in the case of the non-integer power is slightly more tedious and we have not pushed it further (also because in our numerical analysis we saw no evidence of the A 1 = 3/2 asymptotic behavior for the Feynman gauge). Just as an example of the type of results, for A 1 = 0 we find that the coefficients in and so on, leaving B as the only free parameter.
In conclusion, we found four isolated sets of solutions atφ → ±∞. As we will explain in the concluding section, from the point of view of the f (R) theory the most interesting solutions are those in the first class, i.e. (4.10), for which we have found the presence of two degrees of freedom atφ → −∞ and one atφ → +∞ (or the opposite for ω < −1/2).

Fixed singularities
We repeat here the analysis of the analyticity of the differential equation for the Landau gauge, starting with the study of the fixed singularity inφ = 0. Following 4.1.1, we recast the differential equation in its normal form (4.1) and then we expand it in a Laurent series employing a Taylor expansion for the potential. In this gauge we find that at leading order the equation reduces to which vanishes constraining the potential at the origin as 25) or restricting to ω = 0, which is the case we are interested in. Comparing (4.24) with (4.2) we note once more that the case ω = 0 in the Landau gauge is analogous to the case ω = −1/2 in the Feynman, so that the analytic properties of the equation in the two gauges are the same for those two particular values.
For ω = 0 we have now an equation free of singularities. As a consequence, since the equation is unconstrained, we have (for d > 2) two degrees of freedom at the origin,Ṽ (0) andṼ ′ (0), and at least one atφ ± ∞, so that it seems more likely to find global solutions. On the technical side, the absence of a singularity atφ = 0 also means that in this case it is possible to integrate numerically from the origin without employing a MacLaurin expansion.

Movable singularities
As in the Feynman gauge we expect the non linearity of the equation to involve the presence of movable singularities. Since the polynomials P i in equation (4.6) contain the same monomials in both gauges, the analysis carried out in 4.1.2 with the method of the dominant balance still holds and we find in general the singular behavior (4.7) with γ = 3/2. However, because of the gauge dependence of the off-shell effective action, we end up with different coefficients for both the analytic and divergent part. For example, for d = 4 and generic ω we obtain etc. Also similar to the Feynman gauge is the presence of simple pole singularities, with (4.9) replaced byṼ (4.28)

Behavior at large field values
Since the method of the dominant balance leads to similar conclusions for both gauge choices, we expect also for the Landau gauge to find generically an asymptotic solutions of the form (4.29) We can iteratively solve the differential equation for this ansatz, obtaining in d = 4 and so on. As for the other gauge, we see that the coefficients are inversely proportional to A, so that also in this gauge we have to treat separately that case. Before studying those other solutions we focus on the number of free parameters of (4.29), by introducing a perturbation δṼ . We then linearize the equation for the perturbation and study the leading terms, obtaining the equation For ω = 0 the analysis is similar to the one we presented for the Feynman gauge. For ω = 0 we need to include the next order term in the coefficient of δṼ ′′ , and thus consider the equation which admits solutions with the asymptotic behavior The novelty here is that the leading power in the exponent is fourth rather than third order (a consequence of the different power ofφ in the coefficient of δṼ ′′ in (4.32) with respect to (4.31)), so that the solution does not discriminate positive from negativeφ, but rather leads to constraints on A. For A < 0, the solution (4.33) contains an exponential degree of freedom which grows faster then the perturbed function in both positive and negative field regimes, so that we must discard it. Interestingly such sector is the unphysical one, since negative A defines the asymptotic behavior of an unbounded potential. On the other hand, for A > 0 the perturbation is exponentially small both at positive and negativeφ, hence it is always acceptable, and we can work out the subleading corrections as done before for the Feynman case. The higher power in the exponent means that we have to solve more iteration steps before getting to the power-law corrections, but as we do not gain any qualitative insight from such analysis, we do not report further on that, the main message being that now we have two degrees of freedom at both plus and minus infinity. Regarding the case A = 0, making the ansatz (4.19) we find again (d = 4 and ω = 0) the same three special values A 1 = 0, 1 and 3/2, as in the Feynman gauge. The main difference appears in the case A 1 = 3/2, for which the expansion (4.19) now contains one degree of freedom, i.e. b 1 is a free parameter in terms of which all the other b n are expressed: By perturbing around such solution we find that in order to discover new solutions we have to include at least the next-to-leading order coefficients for largeφ in the linear equation, yielding (4.35) whose asymptotic solutions are a superposition of a solution that simply perturbs (4.34), and a series of logarithmic corrections, that carries a second degree of freedom, namely the free parameter c 1 .

Numerical results
In order to find global solutions we integrate out fromφ = 0 and search for a set of initial conditions τ such that the movable singularity goes to infinity in both the positive and negative field region. We present here our analysis for both gauges for ω = 0 and d = 4, starting with the Feynman gauge.

Feynman Gauge
We start a numerical integration at the origin (actually atφ = ±ǫ as explained in Sec. 4.1.1), and similarly to what done in [31], we plot the location at which we hit a singularity, as a function of the free parameter τ =Ṽ ′ (0). When we see a spike in such a plot, we interpret it as a hint of a possible global solution. Since spikes can occur as artifacts due to the scale of the plot, ending instead at a finite value, the next step is to show that such spike can be made arbitrarily long by increasing the numerical precision and by refining the mesh. In addition, in our case we have to produce such type of plots at both positive and negativeφ, looking for spikes that occur at the same value of τ in both ranges. At negativeφ the plot of the singularities looks like in Fig. 1. We apparently find a spike in the negative region for an initial condition τ ∼ 1.638, which however, when zooming in, reveals a richer fine structure, actually three peaks being present (only two of which are shown in the right panel of Fig. 1). Such triple peak can be understood in terms of transition between different types of singular behavior. The most clear explanation is obtained in terms of the numerator and denominator of the normal equation, N and D in (4.1), which we plot in Fig. 2 for four representative cases. We find that for τ 1.638534 and τ 1.638597 both N and D diverge, together with their ratio, at someφ c thus signaling the pole type of singularity found in (4.9). In the range between those two value we find that D vanishes at someφ c , reaching zero with an infinite slope; at the same N reaches a finite value, and we deduce that we are hitting a singularity of the type (4.7) with γ = 3/2. The transitions between γ = −1 and γ = 3/2 coincide with two of the peaks observed in the fine structure of Fig. 1. We interpret the remaining spike at τ ∼ 1.638591 as signaling a transition (as τ increases) from a regime in which N is always positive, to one in which it changes sign twice before hitting hittingφ c . As seen in the zoomed plot in Fig. 1, spikes can be pushed farther away from the origin, however, high precision is needed and we have not tried to reach much beyondφ c ∼ −0.1. In fact, it turns out that a more detailed investigation of the spikes is not worth, as the remaining part of the plot, for positiveφ, turns out to be quite disappointing. Integrating in the positive field region for any initial condition, including the neighborhood of τ ∼ 1.638, we encounter a singularity, as can be seen in Fig. 3, so that we would have not in any case a global solution. Only one type of singular behavior is found in the positive domain, a typical example of which is shown in Fig. 4, and from which we recognize a behavior consistent with (4.7) and γ = 3/2.
We did not find other spikes in both negative and positive region for other values of τ (outside the plot range in Fig. 3), so that in the end we conclude that there are no global solutions in d = 4 and ω = 0 in the Feynman gauge.

Landau Gauge
The search of global solutions is more complicated in the Landau gauge since we have two degrees of freedom at the origin. In order to search for fixed points we adopted the following strategy: i) we integrate numerically from the origin (since there is no fixed singularity we can directly impose initial conditions atφ = 0) for a fixed value ofṼ (0) varying the initial condition τ =Ṽ ′ (0), ii) we repeat the integration for a discrete set of positive and negative values ofṼ (0). As for the Feynman gauge we restrict our research to ω = 0 and d = 4.
We start withṼ (0) > 0, for which we illustrate a representative outcome at negativeφ in Fig. 5.
In this case we find a spike at τ = 1.5 and a continuum set of analytic solutions occurring for τ < τ c , where τ c is a critical value which depends on the initial conditionṼ (0), i.e. τ c ≡ τ c (Ṽ (0)). The peak at τ = 1.5 actually corresponds to an exact solution of the differential equation in normal form, which for generic d > 2 is given by the simple linear functioñ  v(φ) (implying also that (5.1) does not admit linear perturbations). We are thus led to deem (5.1) unacceptable. Regarding the continuum set at negativeφ, we find it for an initial conditions τ smaller then a critical value τ c which, as we already mentioned, depends on the value of the initial conditioñ V (0). VaryingṼ (0) we observed the value of τ c to oscillate between a minimum value τ min ∼ 0.96 and a maximum τ max ∼ 1.12. By increasing the numerical precision we were able to prolong at will the entire group of solutions and we found all of them to behave asymptotically as Aφ 2 , being A a function of the initial conditions.
A typical solution is illustrated in Fig. 6. The seemingly sharp edge in the second derivative is actually an optical artifact: working at high precision, and zooming around the edge one finds that the curve is smooth, as depicted in Fig. 7. We can understand the presence of such a short-scale transition as the rapid vanishing at largeφ of the exponential part of the solutions we discussed As it can be seen in Fig. 5 all the numerical integrations performed using initial conditions with τ > τ c lead (with the exception of τ = 3/2) to a singularity, which we found to be characterized by the exponent γ = 3/2. An acurate analysis reveals a transition in the way the solutions behave before reaching the movable singularity (i.e. the large field regime of the solution), from V (φ) ∼ A φ 2 at τ ∼ τ c , toṼ (φ) ∼ 3 2φ at τ ∼ 3/2. Such transition, together with the spurious solution (5.1), makes the equation particularly stiff around τ = 3/2, as it can be seen from the noise in Fig. 5. However, because of the presence of a singularity we did not put much effort on a more precise numerical integration of this group of solutions.
Integrating towards positiveφ we discover an interesting situation: forṼ (0) > 0 no solutions meet any singularity. We were able to push the integration to arbitrarily largeφ > 0 without encountering singularities for all values τ , and we found solutions with τ < 3/2 to behave asymptotically likeṼ (φ) ∼ 3 2φ , and solution with τ > 3/2 to go asṼ (φ) ∼ Aφ 2 . Combining our findings for positive and negativeφ we conclude that the solutions withṼ (0) > 0 and τ < τ c form a continuous set of global solutions. Atφ = 0 andṼ (0) = 0 the equation is singular. Imposing an analyticity condition at the origin we find that τ = (1 ± √ 19)/4. We did not study these special solutions in detail. ForṼ (0) < 0 the typical situation is depicted in Fig. 8. All the singular solutions we found, for both positive and negative field values, diverge with exponent γ = 3/2. We found in the positive field region a continuum of solutions which do not end on a movable singularity for τ > 3/2, while at negativeφ we met no singularity for τ < 3/2, in both cases with an asymptotic behavior V (φ) ∼ 3 2φ . The two sets have no overlap, hence there are no global solutions in this case. In conclusion, in the Landau gauge in d = 4 and ω = 0, we found a two parameter family of global solutions forṼ (0) > 0 and τ < τ c (Ṽ (0)). Such result could have been expected to some extent, as in the Landau gauge we have no fixed singularity at the origin, and we have at least two classes of asymptotic behavior with two degrees of freedom each at both positive and negativẽ φ. The global solutions we found behave asymptotically asṼ (φ) ∼ Aφ 2 forφ → −∞, and as V (φ) ∼ 3 2φ forφ → +∞. The latter is an indication of an unusual character of such solutions, as that type of asymptotic behavior is the result of a balance between the classical and quantum parts of the RG equation, to be contrasted to the usual situation, where for k → 0 (i.e. the large field regime) only the classical part survives.

Conclusions
In this article we have presented a study of the Brans-Dicke theory (2.1) for an arbitrary potential V (φ) in the framework of the functional renormalization group. We have derived a differential equation in the local potential approximation for a generic parameter ω and dimension d, subsequently focusing our analysis on the case ω = 0 and d = 4, because of its classical equivalence with the f (R) theory. The main motivation for this work came from the asymptotic safety scenario of quantum gravity, which in the literature has been investigated mostly in the pure metric formulation, by means of truncations of the exact renormalization group equations. An important test for such approximation methods would be to show that at least some subclasses of truncations correspond to a series expansion of a functional approximation that explores an infinite dimensional subspace of the full theory space, a classical successful example of that being the local potential approximation in scalar field theory. However, for the case of gravity, such direction is progressing slowly because of the notorious difficulties related to working on curved backgrounds. The idea we have explored in this paper was to exploit the classical equivalence between Brans-Dicke theory and f (R) gravity in order to be able to study a functional approximation of gravity on a flat background. Besides such motivation, a study of alternative formulations of gravity or modified theories is interesting in its own, and studies like the one we presented here can help address the question of quantum equivalence of different formulations.
In order to achieve our goals we have evaluated the renormalization group equation for a generic potentialṼ (φ) in two different gauge choices, namely a Feynman and a Landau gauge, allowing us to discern possible gauge artifacts. As a result of our study, we found a number of important differences between the two gauges. In particular, we found no global solutions of the fixed point equation in the Feynman gauge, whereas we found a two-parameter family of global solutions in the Landau gauge. While some gauge dependence was expected (due to the approximations employed and to the fact of working off-shell, see for example [51,47]), we would have expected that at least qualitative features like the number of fixed points, and of associated relevant directions would be gauge independent (in principle together with any observable quantity, but in practice this property is expected to hold only approximately due to the approximations used). Being the results in our two gauges so different even at a qualitative level, we are led to infer some inconsistency of the model under consideration in the present approximation. Motivated by the relation to f (R) gravity, we did not analyze the case ω = 0 in detail, but we can identify the freezing of the Brans-Dicke parameter to ω = 0 as the culprit of the inconsistent scenario we uncovered. We expect the strong gauge dependence to be lifted once the Brans-Dicke parameter is promoted to a running coupling ω k , in the sense that in any gauge there will be some critical value ω c where something special happens (e.g. a discrete or continuous set of fixed points appears), the value of ω c being gauge dependent, but not so the overall picture. 6 For example, we already know that in the Feynman gauge the value ω = −1/2 gives very similar results to the Landau gauge at ω = 0, and it would be interesting to test whether such critical values correspond to fixed points of ω k for the two gauges, reached either in the UV or in the IR. In view of our results and of the possible solution we just outlined, we can draw an important conclusion: due to its renormalization group flow, the Brans-Dicke theory at the quantum level needs a running coupling ω k = 0 in order to be consistent. Since for ω = 0 the equivalence with the f (R) theory is broken, we are led to suggest that Brans-Dicke theory and f (R) gravity are inequivalent at the quantum level. Needless to say, this should not be intended as a proof of inequivalence, but rather as a logical interpretation of the results we found.
We should point out another aspect which also hints to a non-equivalence of Brans-Dicke theory and f (R) gravity at the quantum level. As we explained, the condition for a solution of the FRGE to be a valid fixed point is that it should be a global solution. While it is quite clear from our analysis that, at least within the present approximation, no nontrivial fixed point can be found for the Brans-Dicke theory at ω = 0 in the Feynman gauge, we should be careful in translating such statement back into f (R) gravity. Due to the nonlinearity of the Legendre transform it could happen that a problematic singularity in one theory would turn into a harmless one in the other, or vice versa. We should indeed remember that the following relations hold (here in dimensionless variables): As a consequence, if a singular point |φ c | < ∞ is such that the first derivative of the potential is divergent, then in the f (R) theory it simply means thatφ c is mapped toR c = ±∞, depending on the sign ofṼ ′ (φ c ). Although that would correspond to a strange situation in whichf ′ (R) does not diverge at infinity (usually the asymptotic behavior is a power law dictated by the tree level part of the equation [16,18,19], implying that at infinityf ′ (R) diverges for any d > 2), that would not be something we can discard as unacceptable. This is precisely what happens in reverse for the Landau gauge: we found global solutions forṼ (φ), but their first derivative is such that asymptoticallyṼ ′ (φ) ∼ 3/2 forφ → +∞, and thus their transform would lead to an f (R) theory valid only up toR c = 3/2. On the other hand, if the potential is such that only its derivatives of order greater or equal to two are divergent, then the singular point is mapped to |R c | < ∞, and thus also the transform of the potential is not a global function. The latter is precisely the case for the Feynman gauge, for which we saw that the singularities at positiveφ are characterized by an exponent γ = 3/2, that is, they have a finite first derivative at the singular point.
Regardless of its connection to the f (R) approximation, the study of Brans-Dicke theory is interesting in its own, and, being a nonrenormalizable theory, it is natural to wonder whether an asymptotic safety scenario applies to it. From such point of view, we should emphasize that what we have presented here is the result of the leading order in an approximation which should be systematically improved. The local potential approximation we employed can be considered, in fact, as a "double LPA" since we neglected both the renormalization of the coupling Z of the operator φ R (having set from the start Z = 1) and of the parameter ω. Both could be promoted to functions Z(φ) and ω(φ), thus leading to a next-to-leading order approximation which could uncover an anomalous scaling of φ and the existence of nontrivial fixed points.

A The two-dimensional case
In two dimensions the fixed point equations in both gauges reduce to ω-independent first order equations. The analysis is thus quite different in this case, it is actually much easier, and we can proceed mostly by analytical means.
Explicitly, the equations in d = 2 reduce to 7 for the Feynman gauge, and toṼ for the Landau gauge. Both equations can be easily integrated, leading to algebraic equations implicitly defining the solutionṼ (φ). As equation (A.1) is slightly more complicated to study than equation (A.2), but at the end it leads to very similar results, we will present the explicit analysis only for the Landau gauge. The fact that the two gauges lead to similar results in this ωindependent case supports our hypothesis that in higher dimensions the strong gauge dependence we found is an artifact of the restriction to ω = 0. The constant of integration C = 4 π (v 0 − φ 0 ) + log y 0 parametrizes a one-parameter family of global solutions, which hence are all acceptable fixed points. The asymptotic behavior of the Lambert function is such thatṼ (φ) ∼ φ for φ → +∞, andṼ (φ) ∼ e 4πφ+C for φ → −∞.
We can study the linear perturbations around such fixed points, by writing as usual 7 The field is dimensionless in two dimensions, hence we omit the tilde.
A being an arbitrary normalization constant, which we can fix to one. Given the exponential fall-off at φ ∼ −∞ of the fixed point solutionṼ (φ), we see that we must impose the constraint λ ≤ 2 in order to avoid exponentially growing perturbations. Indeed the asymptotic behavior of the eigenperturbations is v(φ) ∼ 1 − 2−λ 2 (4πφ) −1 for φ → +∞, and v(φ) ∼ (4π) for φ → −∞. Apart from the upper bound on λ, we do not have other restrictions, hence the perturbations form a continuous spectrum. However, for λ < 2 all the perturbations are redundant, corresponding to a field redefinition φ → φ + (Ṽ ′ (φ)) −λ/2 . We are left with only one essential perturbation, the constant one, v(φ) = A. One special solution of the fixed point equation isṼ (φ) = 0, whose linear perturbations satisfy with solutions v(φ) = A e (2−λ)2πφ . In order to avoid exponentially growing solutions in this case we have to restrict to λ = 2, that is, the only allowed perturbation is again a constant potential, which is a relevant perturbation, and which actually is an exact solution of the full flow equation.
We conclude noting that all the solutions in d = 2 do not admit an f (R) interpretation, as (6.1) together with the asymptotic behavior of the fixed point solutions imply thatR ∈ (0, 1).
The departure from f (R) is of course most evident in theṼ (φ) = 0 case, where the equation of motion obtained by varying φ is simply R = 0.