Solar System Tests in Brans-Dicke and Palatini f(R)-theories

We compare Mercury's precession test in standard General Relativity (GR), Brans-Dicke theories (BD), and Palatini f(R)-theories. We avoid post Newtonian (PN) approximation and compute exact precession in these theories. We show that the well-known mathematical equivalence between Palatini f(R)-theories and a specific subset of BD theories does not extend to a really physical equivalence among theories since equivalent models still allow different incompatible precession for Mercury depending on the solution one chooses. As a result one cannot use BD equivalence to rule out Palatini f(R)-theories. On the contrary, we directly discuss that Palatini f(R)-theories can (and specific models do) easily pass Solar System tests as Mercury's precession.


Introduction
Standard general relativity (GR) has proven to be a reliable theory of gravitation in many an instance (see also [1]). All predictions of the theory in vacuum have been confirmed up to the experimental precision, both locally in the Solar System as in astrophysical regimes such as in binary systems.
However, when considering non-vacuum solutions, as one does in galaxies and cosmology, standard GR is still successful, yet at the price of introducing dark sources (dark matter and dark energy) to fit observations [see 2,3,4,5,6,7]. While there are overwhelming pieces of evidence of the gravitational effects of such sources, there is still no direct evidence of their fundamental constituents.
For these reasons, researchers have been considering the possibility that the effects which are currently ascribed to dark sources may be in fact purely gravitational effects due to modifications of the gravitational interaction itself. There are a number of candidate models which collectively are called modified gravitational theories or extended theories of gravitation [see 8,9,10,11,12,13,14,15,16]. In these models there are no fundamental dark source fields or particles. It is the gravitational interaction that is modified with respect to standard GR, and the modification is designed so that it preserves the local predictions in vacuum, while deviations at different scales justify observations in terms of effective sources.
In this paper we consider two specific classes of modified models, Brans-Dicke models (BD) [17] and Palatini f (R)-theories [see 18,19,20,21], especially in relation with classical Solar System tests, in particular Mercury's precession. The discussion aims also at highlighting that specifically in gravitational theories, observables are model dependent and when data are available, one needs to make predictions for each test within each extended model, making explicit all choices about observational protocols, since such choices may be different from what one does in standard GR.
The two modified models, Brans-Dicke and Palatini f (R)-theories, have been chosen because Palatini f (R)-theories are known to be dynamically equivalent via a conformal transformation, to a subset of Brans-Dicke models. Since Brans-Dicke theories have historically been used as a benchmark for Solar System tests, its parameters have been experimentally constrained and it has been shown that the preferred values of parameters correspond to standard GR. Now it happens that, by dynamical equivalence, Palatini f (R)-theories correspond to a subset of Brans-Dicke theories (with a specific potential as well as) with a value of parameter which is not compatible with Solar System tests. Hence this equivalence has been used to rule out all Palatini f (R)-theories. This is not the only argument for ruling out Palatini f (R)-theories. It has also been argued (and confuted) that Palatini f (R)-theories leads to singularities in polytropic stars [see 22,23,24,25].
In this paper we shall show in details how this is wrong due to multiple reasons interacting with each other. First, the dynamical equivalence requires a potential which was not assumed in the original Solar System tests analysis. Secondly, the value of the parameter for which the dynamical equivalence occurs is a singular value for Brans-Dicke models. Since the original analysis of Solar System tests in Brans-Dicke theories has been performed for generic regular values, one cannot even say that the singular value of the parameter has been ruled out (and in fact we shall show it is allowed). Lastly, we highlight how, even in view of the dynamical equivalence, test particles (e.g. Mercury itself) in the two theories are expected to go along different worldlines and in the specific example they do not only because one is considering a vacuum solution (so that the conformal factor is constant).
In general, we argue that dynamical equivalence may or may not extend to a complete physical equivalence, in which case the equivalence should preserve the action principle, as well as all the independent choices which define observational protocols [see also 26].
Test particles are an independent choice: when you use the eikonal approximation, equations for test particles are obtained, though they are not invariant with respect to redefinition of fields [e.g. 27,28]. Accordingly, the choice of test particle equations is just transformed into the choice of which field corresponds to test particle. Moreover, often one does not have a clear Lagrangian description of test particles in terms of fields and still one uses test particles.
Another choice, is space-time decomposition. In a relativistic theory there is no time and no space, just spacetime. Each observer may split spacetime into space-and-time, though each in a different way. Space and time lengths are thus relative to the choice and conventional. We use them extensively in astrophysics, just because a relativistic theory has no Dirac observables [29]. If we fix a space-time decomposition, that partially breaks general covariance, it conventionally reduces the symmetry group, so that non-trivial relative observables may be allowed.
One can show that defining atomic clocks then a space-time decomposition follows [see 30,31], and one can define space and time lengths out of each specific atomic clock. That is very well known in standard GR, though it extends to a Weyl geometry [32], as required in Palatini f (R)-theories.
In a previous paper (see [33]; see also [34]) we discussed cosmology in a particular Palatini f (R)theory, based on the function We discussed SNIa fit and showed that the system is quite strongly degenerate. If we provide α and β, then the SNIa dataset allows to fix γ, by the way to a value of about γ 10 −104 m −4 . That calls for independent measurements to fix α and β. Here we used Solar System tests to reduce degeneracy. It has also been argued (by an anonymous referee) that the best fit value α 0.1 we found there would fail in Solar System tests, we show here that this is not the case considering Mercury test.
Material is organised as follows. In Section 2, we fix notation in BD theories and Palatini f (R)theories. In Section 3, we review the dynamical equivalence. In Section 4, we consider static, spherically symmetric solutions which will be used to model the Solar System. In Section 5, we consider geodesic equation, first for generic static, spherically symmetric metric, then for a solution.

Brans-Dicke and Palatini f (R)-theories
A Brans-Diche (BD) theory is a gravitational theory for a metric g and a scalar field ϕ. The action principle in dimension m = dim(M ) = 4 is where R is the scalar curvature of g and U (ϕ) is a potential. The parameter ω is called the BD parameter.
From this action, one has vacuum field equations We shall eventually be interested also in the special case ω = − 3 2 in which field equations become which is of course a degenerate value since for ω = − 3 2 the field equation for ϕ drops order and becomes an algebraic equation depending on the potential U (ϕ). When no potential is assumed in the degenerate case ω = − 3 2 , the scalar field ϕ is left undetermined. The original analysis of Solar System tests [see 35] was carried over with no potential and for a generic regular value of ω.
In BD models, test particles go along timelike geodesics of g, which determines the geometry of spacetime as well as its metric structure. The scalar field ϕ is non-minimally coupled and it modifies the law in which gravitational field is mediated.

Palatini f (R)-theory
For a Palatini f (R)-theory, we start from fields (gµν ,Γ µν ), whereΓ is a (here torsionless) generic connection on the spacetime M , a priori independent of g, and a Lagrangian for some (regular enough) function f (R) of the scalar curvature R = g µνR µν , whereRµν is the Ricci tensor of the connectionΓ alone. Field equations read By tracing the first field equation by g µν , we obtain the so-called master equation, which must be identically satisfied along solutions. The function f (R) is called regular enough when the zeros of the master equation are simple and they form a discrete set.
To solve the second equation, one has to set ϕ ∝ f (R). That can be locally inverted as R ∝ r(ϕ).
For the model based on (1), we have the master equation With these values of the curvature we get and At that point the (vacuum) field equations are where the last approximations have been done for small values of βγ with respect to α √ αγ (which both have units of an inverse squared lenght) and small values of β √ αγ with respect to α 2 (which are both adimensional) Let us notice that, as a consequence of the assumption of being in vacuum, the conformal factor ϕ is constant and field equations reduce to Einstein equations with a comslogical constant, which is proven in general by universality theorem [see 36]. Universality theorem guarantees that vacuum solutions of Palatini f (R)-theories maintain the successes shown by standard GR solutions provided that the cosmological constant is small enough, as small as it is expected by observations.
In the specific example (1), the (effective) cosmological constant Λ ± is small enough iff γ is small enough and γ is the parameter which is better constrained by SNIa. We can consider the best value of γ γ 2.46 +3.84 −2.24 · 10 −104 m −4 (12) obtained in (1) setting α = 0.095, β = 0.25 m −2 . As a matter of fact, the predicted best fit value of Λ = 1.27 +0.58 −0.98 · 10 −52 m −2 is not far away from the value observed, e.g. by the one found by Planck survey (Λ Planck(2018) = (1.106 ± 0.023) · 10 −52 m −2 ) [see 37] This value for the (effective) cosmological constant will not be observable (or falsifiable) with experiments in the Solar System. If it were, then ΛCDM would be falsified as well [see 38,39].
Thus let us review the dynamical equivalence and then investigate how this rough though clear result can be compatible with exclusion of the values ω < 4 · 10 4 [see 40,41], hence including ω = − 3 2 , by BD and Mercury's precession.
3 Equivalence between Palatini f (R)-theories and BD Before going to field equations, one can prove dynamical equivalence at the level of the action. The proof of equivalence is in two steps. In the first step we show equivalence between Palatini f (R)-theory and a theory with an extra scalar field governed by the Helmholtz Lagrangian. In the second step we recast Helmholtz Lagrangian as a BD model by a suitable field transformation [see 16,27]. Let us consider the definition of the conformal factor ϕ = f (R) and solve it for the curvature R, Let us fix ϕ * = α + 3 3β 2 γ 4 . For ϕ < ϕ * , the equations has only one positive solution R = + R(ϕ).  ii) R * ≤ R < 0 and ϕ ≥ ϕ * ; R = r2(ϕ).
The Helmholtz Lagrangians are dynamically equivalent to f (R)-theory. By varying Helmholtz Lagrangian with respect to (g,Γ, ϕ), one gets field equations the last equation being equivalent to a branch of ϕ = f (R). Then we can defineg = ϕg and solve the second to getΓ = {g}. Finally, from the first one, one can trace to get the master equation (and get that R and hence ϕ are constant on shell). Since ϕ is constant {g} = {g} andR (µν) = R µν . Then, consequently, the first equation becomes Einstein with cosmological constant, equivalent to (11). As a matter of fact, this gives us the chance to test the Palatini f (R) model directly without passing through BD equivalence. Moreover, the other way around, this gives us also the chance to use equivalence to test the degenerate BD models which are equivalent to Palatini f (R) models. In both cases, we have that the value of the cosmological constant one gets from f (R) given by (1), with parameters (12), is compatible with what found by Planck survey in 2018. Accordingly, the models are not rules out by solar system tests (as well as by other test which are insensitive to the observed value of the cosmological constant).
If we select a specific potential U = f (r(ϕ)) − ϕr(ϕ) in a Brans-Dicke theory we get equations: With that potential, one has which in fact is the left hand side of the master equation.
Ricci tensor conformal transformations are given by Thus the first equation reads as By tracing the first equation, we get For There are solutions of the second equation with a non-constant ϕ. For The second equation implies that ϕ is constant, thenR µν = R µν . By tracing the first, one gets the master equation, from which R is constant. Then the first equation becomes Einstein with cosmological constant.
That shows that there is a dynamical equivalence between Palatini f (R)-theories and BD theories with ω = − 3 2 and a potential U (ϕ) = f (r(ϕ)) − ϕr(ϕ). For the function f (R) given by (1), the function r(ϕ) is quite complicated, but in fact we do not really need it. Anyway, currently we are on the positive branch (the first one) and as long as observational cosmology is concerned we can restrict to that branch. For a proof at the level of action see [16,27].

Solutions with point-like sources
Let us here consider the static spherically symmetric solutions in Palatini f (R)-theory and BD theory, both in the case of generic parameter and no potential and for ω = − 3 2 and the potential U (ϕ) = f (r(ϕ)) − ϕr(ϕ) induced by dynamical equivalence.

Solution in Palatini f (R)-theory
In a Palatini f (R)-theory, in view of universality theorem, we get a static spherically symmetric solution which isg where we set dΩ 2 := dθ 2 + sin 2 (θ)dφ 2 for the volume element on the sphere. Thus In view of the fact that the conformal factor is constant, we can change coordinates by ( √ ϕr = r, √ ϕt = t) and obtain If we want a specific asymptotic behaviour, for example a metric which is asymptotically antide-Sitter, we can set a = 1 in the function A(r).
Since the solution of BD theory will be given in isotropic coordinates, let us first recast this solution in isotropic coordinates (t, ρ, θ, φ). Any static, spherically symmetric metric can be recast in isotropic form by a change of radial coordinate ρ = ρ(r) (hence dρ = ρ dr).
One has simply and comparing with the expression in pseudo-spherical coordinates one gets the conditions Hence one can integrate the last condition to get ρ(r) and then set A(ρ) := A(r(ρ)).
For example, the Schwarzschild metric in pseudo-spherical coordinates (t, r, θ, φ) is while in isotropic coordinates it reads as where the integration constant has been fixed to have lim r→+∞ ρ r = 1. Let us remark that for ρ → +∞ we get Minkowski metric.

Solution in BD theory
We can find a static and isotropic solution of BD equations by an ansatz in isotropic coordinates (t, ρ, θ, φ), namely If we fix the potential to be zero and ω = − 3 2 , we can solve the BD equation [see 42,43] as where we have defined That corresponds to what Weinberg does [35], though with no approximations. It is a solution for any (α 0 , β 0 , ϕ 0 , C) and λ is computed with the identity above. If we want to get the Schwarzschild metric at infinity, we need to set (α 0 = β 0 = 1). Even considering α 0 = β 0 = 1, we see that there is a 1-parameter family of static and spherically symmetric solutions, parameterized by C, unlike what happens in standard GR where the solution is unique. We see that there are two families of solutions: one with C = 0, and consequently λ = ±1, where the conformal factor ϕ is constant and indeed, when λ = 1, reduces to the Schwarzschild metric.
The second family is for C = 0 which has a non-constant conformal factor whenever m = 0. Accordingly, we can say that Schwarzschild solution is always there, although as a somehow isolated solution, while BD theory allows a whole family of solutions with a non-constant conformal factor.
For ω = − 3 2 and no potential, one has an arbitrary conformal factor and the Schwarzschild metrics only. If the potential is added as in view of the dynamical equivalence, then the conformal factor is frozen to be constant and the Schwarzschild-de Sitter solution is obtained.
In what follows we shall consider a sub-family of the general solution with C = 1/n, and consequently λ 2 = [2(n 2 + n + 1) + ω]/(2n 2 ) (which for ω → − 3 2 gives λ → 2n+2 2n , while for n → ∞ gives λ = ±1 which is the value associated to the Schwarzschild solution) and we shall compute an observable, such as the precession rate of Mercury. We shall show that even in the limit ω → − 3 2 such an observable is discontinuous and it does not reproduce the result for the Schwarzschild-de Sitter solution obtained by setting ω = − 3 2 . That eventually justifies the claim that the value ω = − 3 2 is degenerate and one cannot use a limit procedure to infer the result in the degenerate case, and then neither in the dynamically equivalent Palatini f (R)-theory (which in fact will pass the Mercury test).

Geodetics and exact relativistic Kepler laws
We shall predict precession of Mercury by solving a more general problem, i.e. writing exact Kepler laws for an arbitrary static spherically symmetric metric in isotropic coordinates restricted to the equatorial plane θ = π 2 , g = −c 2 A(r)dt 2 + B(r)dr 2 + r 2 dφ 2 .
Precession comes from failure of first Kepler laws, since it measures by how much the orbits fail to be closed. Orbits are g-timelike g-geodesics, in an allowed region where r is bounded from below and above. We could use two Lagrangians. The first Lagrangian is for geodesics parameterised by proper This Lagrangian has 3 d.o.f. (t, r, φ) and 3 first integrals (P, H, K), namely which we can solve as dr dτ The second Lagrangian is invariant with respect to reparameterisations This Lagrangian has 2 d.o.f. (r, φ) and 2 first integrals (E, J).
We have Now we can map the values of the first integrals in the two frameworks From the invariant Lagrangian we have where the former equation gives the time evolution and the latter the orbit trajectory.
Having orbits is related to having an allowed region [r − , r + ] in Φ (as well as in Ψ). By integrating the latter in the allowed region [r − , r + ], one gets 2φ = 2π +δ, where δ is the precession per orbit. When the precession is zero, the orbit is closed, i.e. the classic first Kepler laws. Accordingly, the precession δ is expected to be smaller and smaller getting away from the source, though not in the limit unless the solution is asymptotically flat (namely, Λ = 0). In the limit to standard GR, Mercury's precession should approach ∼ 43 arcsec/century. For the Earth precession should be negligible.
The second Kepler law is related to conservation of angular momentum, which indeed is conserved exactly also in the relativistic regime. Finally, we can get the period T (r) from the integral in the allowed region. In the Newtonian approximation one has T 2 ∝ (r − + r + ) 3 . The function T (r) contains the exact law, which in the limit must reproduce the Newtonian prediction.
Let us now make this explicit in some metrics which are relevant for the theories we are considering.

Schwarzschild
The Schwarzschild solution is defined as One has [b] = L and b is called the Schwarzschild radius.
Since the mass of the Sun is M = 1.9885 · 10 30 kg, its Schwarzschild radius is b = 2953.29 m.
The function Φ is then given by This function to the power − 1 2 will be integrated. In order to have an analytical expression of that integral it will be better to factorize the polynomial, i.e. expressing it as a function of perihelion and aphelion (r ± ) instead of as a function of (E, J). If we want a bounded orbit, the function Φ must be negative at big r, thus c 2 E 2 < 1, i.e. −1 < cE < 0. The polynomial p(r) = (c 2 E 2 − 1)r 3 + br 2 − J 2 r + bJ 2 has a zero since it is of odd degree. We have Hence, we can set initial conditions in the zero r = r 0 , whereṙ 0 = 0, and set φ 0 = 0, φ 0 = v 0 /r 0 . For these initial conditions, we have where we set A 0 = A(r 0 ). For having a bounded orbit we must have Let us also notice that (ξ 0 +b) 2 c 2 < c 2 . Accordingly, if we start from zero initial speed v0 = 0 and we increase it, at some point the orbit will become unbounded.
With these values of course we can express Φ(r, r 0 , v 0 ) and factorise out of it (r − r 0 ): If Φ has only one zero r = r 0 > b, the allowed region is [b, r 0 ] and the geodesic is falling towards the asymptotic goal at r = b, i.e. the horizon. To have a bounded periodic orbit the allowed region must be [b, r 0 ] ∪ [r − , r + ]. Necessary condition for that to happen is Let us remark that 0 < v 2 m < v 2 M < c 2 . In fact, v 2 m is obviously positive. It is also Hence, starting from zero speed and increasing it, the test particle will fall in until it reaches a limit speed vm and it will stay bounded until a new limit v M . Accordingly, one has bounded orbits for v 2 Then for v 0 = v m we have a new zero of Φ appearing in r 1 = 2br0 r0−b , which is less that r 0 , for any r 0 > 2b.
One has that it vanishes for ξ0 = b and it is negative for ξ0 > b since it goes to −∞ as ξ → ∞. Accordingly, unless b < r0 < 2b, for r0 > 2b one has the first stable orbit (for v0 = vm) in ξ ∈ [ξ1, ξ0]. Then, increasing v0 the perihelion grows to infinity. At some point il will reach and pass ξ0, then it will grow to infinity which is reached at v0 = v M .
When the perhelion becomes equal to the aphelion (circular orbit) and it keeps growing they get exchanged with each other. When we originally choose a zero r 0 of Φ we can always choose the biggest one. In this way we are interested to find the speed v c = v 0 for which we have circular motion. That happens when Φ(r 0 ) = 0 twice, i.e. for v 2 c = c 2 b 2r0 .
When r0 > 3b, we have vm < vc < v M . Once again for r0 < 3b we are too near and Newtonian dynamics is not a good approximation.
Let us summarise. We choose initial position at the aphelion, for r = r + > 3b. 1 i) For 0 < v 0 < v m (thus low angular momentum), the test particle falls in, end of the story.
ii) For v 0 = v m , we have the first "orbit" with a perihelion of r − = 2br+ r+−b < r + (which is not an orbit yet, since r − is an asymptotic goal). In what follows we are interested in iii) and iv) only since they capture all bounded orbits r ∈ [r − , r + ], parameterised by r + and r − . By the way, if we fix r − and r + the initial velocity we need for that is v 2 0 = c 2 br 2 − (r+−b) r 2 + (r++r−)(r−−b) . The first integrals for the bounded orbits are Accordingly, we can express the function Φ in terms of r ± as Ψ(r, r ± ) = − r 2 + r 2 − r(r + − r)(r − r − ) ((br − − r − r + + br + )r + br + r − ) .
Alright then, we have expressed everything in terms of r ± .
1 What happens for r + < 3b is not our concern now.
This is computed, theoretical quantities, not measured. We are not meaning we can observe the Earth period up to 10 −7 s. We mean that we can predict its value with arbitrary prediction so that we can compare observed value with it.
The Earth precession per orbit is The Mercury's precession per orbit is implying a precession per century (i.e. the cumulative precession for 100 Earth's periods T ♁ , which corresponds to about 415.21 Mercury's periods) of δ = 42.98 arcsec. Let us stress that comparing with Earth's period to define century allows us to avoid any direct reference to an external clocks. In some sense we are using relational time within the system itself.
To find third law we can consider Φ and Ψ, make the substitutions r → ρ/x, r ± → (1 ± e)/x and expand at x = 0 + (i.e. far away from the central mass). The first term in the series reproduces the corresponding Kepler function. The second term in the series gives corrections. This is a very convenient technique to define what is the Newtonian limit, since it is done before integration. Of course, it requires one does Kepler case first.

Post Newtonian approximation
One usually does is expanding the metric coefficients in isotropic coordinates in series of M G/ρ and assume it is not very different from Minkowski, i.e. the weak field approximation, namely This approximation is good enough for Mercury, it would fail too close to a black hole. In this approximation the Schwarzschild solution is recovered for α = β = γ = 1, while the BD solution is obtained for It is clear that in the very same moment one expands in series, the ability to spot isolated singular solutions is lost. Of course, for ω → ∞, one gets γ = 1 for the Schwarzschild solution, but for ω = −3/2 one gets γ = −1. As a matter of fact, one is assuming a priori to be on the main regular sequence of solution, which is incorrect (or a partial viewpoint) as we shall see. For this reason it is much better to stick to exact results rather than starting expanding in series. As an extra bonus, by doing it exactly, we also can test how close we need to be to see the strong field effects which, in principle, is a solid prediction of the theory.

Schwarzschild-de Sitter
Let us discuss the Schwarzschild-de Sitter spacetime with A = 1 − b/r − λr 2 . That is a solution of Einstein equation in vacuum with cosmological constant Λ = 3λ. We need A > 0 at least in a region r ∈ [r 1 , r 2 ]. Since we have if λ > 0, we have lim r→+∞ A = −∞ and, if b > 0, we have lim r→0 A = −∞. Thus we would like a region in the middle where A > 0, which is easy to find since and in r = r * we have (71) If λ is too big there is no such a region where A > 0. That means that as λ grows, sooner or later the cosmological horizon will touch the Schwarzschild one.
Thus the cosmological constant has to be small enough. In this case the function A has at least two zeros, hence three (since it is odd degree).
If we want a bounded orbit, p(r) must have 4 zeros, {r m , r + , r − , r M }, in the region [r 1 , r 2 ] so one negative zero r = −(r m + r + + r − + r M ). Hence it is If one knows the mass of the star and cosmological constant, then r 1 and r 2 can be computed. Then one gives a planet with its r ± and can compute out of (r 1 , r 2 , r ± ) the value of (E, J, r m , r M ), i.e. the initial conditions. In other words, while (r 1 , r 2 ) are a convenient way of parameterizing the parameters of the system, namely (m, λ), (r + , r − ) are a convenient way of parameterizing initial conditions of a specific timelike geodesic.
Although it is good to discuss it once from scratch, one can also cut the discussion short by saying that we want to have an orbit, i.e. time-like geodesic that is bounded from above and below. That means we need an allowed region [r−, r+] for Φ and Ψ, so that both r± are simple zeros. Knowing that, at infinity, Φ is definite negative and it is positive around r1 < rm < r− ≤ r+ < r M < r2, we directly get (81) Note that we also have two allowed regions [r1, rm] and [r M , r2] corresponding to the test particle falling in and escaping to infinity, respectively. Anyway, here we are interested in solutions in [r−, r+].
Let us finally remark that the factorised form of Φ is very convenient for analytical computation. Most of the the possibility of treating the problem analytically relies on this factorisation, i.e. in using (r1, r2, r±) instead (λ, m, E, J) and express the rest in terms of them. Now that Φ(r) and Ψ(r) are fixed, one can compute the planet period and the precession per orbit as we did for Schwazschild. That can be done for Earth and for Mercury so to have the ratio between periods. Then we can compute the precession of Mercury in a century, obtaining which shows that the effect of the cosmological constant grows approximately linearly in the (log, log)graph and that with a cosmological constant Λ = Λ+ we are well within the limit in which we cannot observe it in Mercury perihelion. We put a bar to highlight on its left hand side the digits which agrees with the value computed with no cosmological constant. That bar also highlights on the right hand side the digits which are affected by the cosmological constant. Let us notice that we are not simply saying that as long as we consider Mercury's precession we can neglect the cosmological constant. We are in fact computing the relative error one does by neglecting it.

Solution in f (R)
We have considered the Schwarzschild-de Sitter metric and computed precession of Mercury in them. We found that if the cosmological constant is small enough the theoretical precession is compatible with the observed one.
In view of universality theorem we know that (vacuum) Palatini f (R)-theories in fact are equivalent to Einstein with a cosmological constant. However, extra care is needed in this case. Universality theorem claims thatg is a solution of Einstein equations with a cosmological constant the value of which is dictated by the function f (R) via the master equation. Now from the viewpoint of Ehlers Pirani and Schild (EPS) framework (see [44], [45]), one should expectg to govern geodesic motions, while g is related to distances, while in Schwarzschild-de Sitter model (as above) one has a single metric, namely g above, which dictates both geodesic equations and distances. Although in principle one should discuss whether this aspect plays a relevant role when applying the discussion above as it does in general, we have to remark that Solar System is modelled by a vacuum solution of Palatini f (R)-theory, in which hence the conformal factor is constant. Accordingly, {g} = {g}, i.e. the two metrics g andg actually define the same geodesics trajectories, the same timelike directions, and one hasR µν = R µν , i.e.
Accordingly, also g obeys the same field equations asg. We can use g only in vacuum, and apply the result above.
If the function f (R) determines a small enough cosmological constant (as it happens with the function (1) and parameters (12)) then it actually predicts precession of Mercury compatible with the observed one. Let s remark that Mercury test is passed despite the value α 0.095. That directly shows that α 1 is not required to pass this test. In view of dynamical equivalence between Palatini f (R) theory and BD (with a potential and ω = −3/2) this shows also that such a BD theory passes the test as well.
Also in this case, one should pay attention to the fact that in Palatini f (R)-theories geodesics and distances are related to two different metrics while in BD both are related to the same metric g. However, in vacuum, this is not an issue since g-timelikeg-geodesics are also (and the only) g-timelike g-geodesics. However, in non-vacuum solutions (as galactic dynamics or cosmology) the models would be actually different.
Since classical tests have been performed in BD theories (with no potential and generic ω, hence different from ω = −3/2, which is degenerate) and they show that ω must be ω > 10 4 , this was used to try, erroneously as we discussed above, to rule out all Palatini f (R)-theories at once. Besides, we saw directly that this argument is spurious, we shall review the test in BD model on an exact formulation. We have two reasons to do it: first, since we will not use PN approximation for it we can apply it close to the horizon, in the strong field regime. Secondly, we show that the value ω = − 3 2 is degenerate also with respect to the test and the role of the potential cannot be neglected.

Mercury test in BD theory
Let us here consider a BD theory with no potential (and ω = −3/2, as this is used to determine the solution as long as the scalar field is concerned). Test particles (and the planets Earth and Mercury) are assumed to go along geodesics of g. On the equatorial plane, in isotropic coordinates, one has the Lagrangian again with total energy and angular momentum as first integrals. In isotropic coordinates, the coefficients are a bit different to what described so far, so we need to repeat the discussion from scratch. First integrals are which can be inverted for velocities as dρ dt Considering Φ, if ρ → +∞, given that A → 1 and B → 1, we have in order to have bounded orbits. Now that we know how it works, we can either study the function Φ for all parameters and then select parameters which describe bounded orbits, or simply require two solutions Φ(ρ ± ) = 0 that bound an allowed region [ρ − , ρ + ]. Indeed, if we consider the equations Φ(ρ ± ; E, J) = 0, we can solve for (E(ρ ± ), J(ρ ± )) and then replace them back into the Weierstrass functions to obtain dρ dt The second issue to discuss is how to use orbital parameters. Since we are in isotropic coordinates, the coordinate ρ is not endowed with a direct meaning of a distance and it is different from the coordinate r in pseudo-spherical coordinates [see 31]. In view of the form of the metric in spherical coordinates, spacetime is foliated into Euclidean spheres (we mean metric spheres, not only topological spheres) arameterised by (t, r), with the coordinate r being the radius each sphere would have if it were embedded into a Euclidean space. That in fact corresponds to the observation protocol for measuring astronomical distances, which in fact uses Newtonian approximation and classical Kepler laws. Accordingly, the observed orbital parameters have to be related to r, not to ρ. However, we know that and we can determine the values of ρ ± which correspond to the observed r ± for Mercury or the Earth.
In order to consider the solution (36) in BD theory, we know that standard GR corresponds to a big value of ω, C = 0 and λ = 1. In view of the constraint among (ω, λ, C) we can consider a sub-family of solutions given by setting λ = 1, C = − 1 n and computing ω = 2(n − 1) and hence table the exact precession of Mercury for each value of n.
This gives us a (partial) insight on how precession of Mercury depends on (ω, C) in BD theory. We expect to obtain a good agreement with observed values for n → ∞ and observe rather incompatible values for small values of n. We can also consider the limit n → 1 4 , which corresponds to ω → − 3 2 , i.e. the degenerate parameter although with no potential. The whole computation of the precession for a specific value of n involves exact calculations up to a finite number of numerical integrations of converging improper integrals. This is not much different from what one does when studying the graph of a transcendent function, in which a finite number of evaluations of the function (e.g. to determine zeros or critical points) are eventually performed numerically. Accordingly, we can say our computation is an analytical exact result or, if you prefer, call it semi-analytical.
For each value of n we computed the orbital period of the Earth, of Mercury, the precession per orbit of Mercury and finally obtain the precession P (n) of Mercury during 100 orbits of the Earth.
These values are, theoretical predictions in different solutions of BD theories. They can be computed at arbitrary precision. Here we checked that the shown digits are not affected when the overall precision required is increased.
We see that in fact one can distinguish among different values of n (and ω) by means of P which is definitely observable. As the observed value of about 42.98 arcsec century −1 , observations are able to exclude small values of ω, including ω = − 3 2 , as expected. For the degenerate value ω = − 3 2 one has two possible solutions, one for C = 1 4 which sits in the sequence we analyzed and one for C = 0 which does not and it is corresponds to ordinary Schwarzschild. In fact, for the solution with C = 0 the value of ω is undetermined and the corresponding solution is somehow isolated in the solution space. It passes the tests obviously, since it is the same solution of standard GR.
Accordingly, it is not completely correct to say that BD theories with small value of ω, even with no potential, are ruled out by observations. The solutions with C = 0 are, while of course C = 0 is not.

Constraints in parameter space
We consider BD solutions for ω and C. We sample parameter spaces computing precession of Mercury (in arcsec century −1 ) and computing ∆P (ω, C) the relative error (·100) with respect to the observed value of about 42.98045118132.
The values of ω and C which have interestingly low errors are too spread in the (ω, C) plane. It is expected low errors in the limit ω → +∞ and C → 0 that corresponds to standard GR in the non-degenerate sequence. Thus we plot ∆P (ω, C) on the axes x = log(ω), y = − log(C), so that standard GR corresponds to the limit x → +∞ and y → +∞. Indeed, we see in Fig. 2 that in the region x → +∞ and y → +∞ one has the smallest errors, as expected.

Conformal factor and Weyl transformations
We can still explore one possibility. When we used aphelion and perihelion in isotropic coordinates we used the metric g as EPS framework dictates. However, we have two metrics and we should check the difference with what is predicted if we usedg instead. This is relatively easy, all we have to do is determining ρ ± by using the equation instead of (94). In the case n = 1000, we obtain the precession prediction to bẽ P (1000) = 42.95 arcsec century −1 to be compared with the value computed in (95), namely P (1000) = 42.95 arcsec century −1 . Thus we see a tiny difference, a difference yet, as tiny as expected since the conformal factor is very close to 1 at the orbit of Mercury or the Earth. Still, the difference is there and is expected to become bigger in stronger regimes which, by the way, says that the difference between using g org to describe distances is in principle observable.

Conclusion and Perspectives
We considered the standard test of Mercury in different contexts. Our treatment is analytical and we do not resort to weak field approximations so that our framework is still valid for satellites orbiting a black hole at few Schwarzschild radiuses.
It has been argued that since BD theory is ruled out by observations (for small values of ω) and they are dynamically equivalent to Palatini f (R)-theories, then these are ruled out as well.
We showed that this is not the case for a number of reasons. First of all, dynamical equivalence is with a degenerate value of the parameter ω = − 3 2 , moreover with a potential which has not been considered in the original BD test. Moreover, in BD theory one has no Birkhoff theorem so one has a three parameters' family of static spherically symmetric solutions.
We considered standard GR, standard GR with a cosmological constant, Palatini f (R)theories and BD theories. These produce Schwarzschild solution with or without a cosmological constant Λ as well as a more general family (36) of solutions of BD theory.
We showed that Schwarzschild, as well as Schwarzschild-de Sitter solutions pass the test provided that the cosmological constant Λ is small enough. For Palatini f (R)-theories this imposes constraints on the function f (R) which for example are met in some of the models based on (1), among which one has the best fit values (12) found in [33]. Among solutions of BD theory we showed that, besides Schwarzschild solution which is also a solution of BD theory for any value of ω and which passes the test, also the other solutions of BD theory passes the tests provided that ω is big enough.
This shows directly how Mercury test does not technically rule out BD for small ω, it rules out some solutions of it. Of course, the Schwarzschild solution cannot be ruled out, rather one does not need BD theory to have a Schwarzschild solution.
It also shows directly that dynamical equivalence is irrelevant for ruling out Palatini f (R)theories which in fact are ruled out only if they produced too big values of the effective cosmological constant. In particular, in the family (1) considered in [33] where the best fit value (12) of (α, β, γ) were found to model Ia supernovae, the parameter fit were found to be strongly degenerate. In particular, the value of β was found to be poorly localised by SN Ia data (as it can be expecetd since the β parameter has very tiny effect in a universe where supernovae can occur) and also for any given value of α one could find a best value for γ which produces a good agreement with observations.
During the peer review of [33] it has been argued that the best fit value of α 0.1 would fail to model Solar System. Here we showed that this is not the case. The ruling out of theories has to be carried over at the level of observables not of actions. Despite α 0.1 produces a vacuum action which does not approximate the standard GR vacuum action, the constant factor has no effect on observables in vacuum, hence in the Solar System tests. This is also implied by universality theorem for Palatini f (R)-theories.
In fact, if SN Ia fix a best fit value for γ, then Mercury test produces constraints for the ratio γ α . The two set of data in fact remove the degeneracy for α and γ, leaving the one connected to β.
As for future perspectives, we need to extend the analysis to the other classical Solar System test (light deflection and radar delays) to check whether these add constraints. We need to extend our exact approach to these cases so that conclusions are robust against weak field approximations and stay valid in a strong field regime.
Further constraints may arise from tight binary systems in which the weak field approximations can be at stake, as well one can look for consequences in collapse events relevant to gravitational waves. Being our method viable for different theories this would open a way to use gravitational wave phenomenology to be used for reliably distinguish different gravitational theories, especially when much more precise data will be available with new experimental surveys; see e.g. [46], [47].
Under this viewpoint of distinguishing different theories in terms of observables only, this paper is in a series with the aim of discussing validity of a specific family of theories (namely, in this case (1)). Before even discussing the validity of a model, one needs to fix parameters. This is not different from what it is routinely done in QFT, when a general theory of electromagnetic field and its interaction with charged fields has to be calibrated by choosing one cross section (the Compton scattering) to fix the renormalised parameters e m . Only after this calibration one can predict other cross sections and validate or falsify the model.
Gravitational theories are not different even in the classical regime. They have parameters to be calibrated by choosing some conventional observations. What we are doing is choosing SNIa to fix γ, Mercury's precession to fix α. Further investigation must be devoted to fix β (e.g. by using elements formation). Only at that point, even when a model has survived calibration, one can use the calibrated model to predict expected values (e.g. power spectrum) to falsify the theory on the basis of observations.
In other words, that is the way to go in modified gravity. We are adding parameters (e.g. in f (R)-theories we have the (potentially infinite) parameters which fix the function f (R)). In order to get stuck to finitely many parameters, we consider a family of functions at a time. Then we need extra experiments to calibrate the theory (in standard GR we only have G which in fact can be fixed by lab experiments) before we are allowed to say we defined a model. The more parameters we add the harder the calibration, which is what is fair to pay for an extended model of gravitation.
All other heuristic arguments about validity of a theory are based on physical intuition which is often model dependent and moreover it has been developed in standard GR and sometimes uncritically applied to different gravitational theories. We provided above a number of such arguments which eventually have been shown to be inconsistent. Being stuck to observation is the only robust way in gravitational theories (although probably in general) to really falsify a theory.
If we take this approach seriously, the current situation is almost desperate for modified gravity models. Let us close by remark and stress that here and in the series of investigation to come we are still trying to falsify a specific family (1)) of Palatini f (R)-theories. We still have no clue at all about how to do it for a generic Palatini f (R)-theory, which depends on potentially infinitely many parameters, less then ever for a generic modified gravity theory. However, this is what needs to be done. INdAM-GNFM.