# Cosmography and Large Scale Structure byf(R)-Gravity: New Results

Authors: Salvatore Capozziello; Vincenzo Salzano
Journal: Advances in Astronomy
 (2009)

---

## Abstract

The so-calledf(R)-gravity has recently attracted a lot of interest since it could be, in principle, able to explain the accelerated expansion of the Universe without adding unknown forms of dark energy/dark matter but, more simply, extending the General Relativity by generic functions of the Ricci scalar. However, apart several phenomenological models, there is no final f(R)-theory capable of fitting all the observations and addressing all the issues related to the presence of dark energy and dark matter. An alternative approach could be to “reconstruct” the form of f(R) starting from data without imposing particular classes of model. Besides, adopting the same philosophy, we take into
account the possibility that galaxy cluster masses, estimated at X-ray wavelengths, could be explained, without dark matter, reconstructing the weak-field limit of analytic f(R) models. The
corrected gravitational potential, obtained in this approximation, is used to estimate the total mass of
a sample of 12 well-shaped clusters of galaxies.


---

## Body

## 1. Introduction

As soon as astrophysicists realized that Type Ia Supernovae (SNeIa) were standard candles, it appeared evident that their high luminosity should make it possible to build a Hubble diagram, that is, a plot of the distance-redshift relation, over cosmologically interesting distance ranges. Motivated by this attractive consideration, two independent teams started SNeIa surveys leading to the unexpected discovery that the Universe expansion is speeding up rather than decelerating as assumed by the Cosmological Standard Model [1–5]. This surprising result has now been strengthened by more recent data coming from SNeIa surveys [6–13], large scale structure [14–18] and cosmic microwave background (CMBR) anisotropy spectrum [19–25]. This large data set coherently points toward the picture of a spatially flat Universe undergoing an accelerated expansion driven by a dominant negative pressure fluid, typically referred to as dark energy [26].While there is a wide consensus on the above scenario depicted by such good quality data, there is a similarly wide range of contrasting proposals to solve the dark energy puzzle. Surprisingly, the simplest explanation, namely the cosmological constantΛ [27, 28], is also the best one from a statistical point of view [29–31]. Unfortunately, the well known coincidence and 120 orders of magnitude problems render Λ a rather unattractive solution from a theoretical point of view. Inspired by the analogy with inflation, a scalar field ϕ, dubbed quintessence [32, 33], has then been proposed to give a dynamical Λ term in order to both fit the data and avoid the above problems. However, such models are still plagued by difficulties on their own, such as the almost complete freedom in the choice of the scalar field potential and the fine tuning of the initial conditions. Needless to say, a plethora of alternative models are now on the market all sharing the main property to be in agreement with observations, but relying on completely different physics.Notwithstanding their differences, all dark energy models assume that the observed apparent acceleration is the outcome of some unknown ingredient, at fundamental level, to be added to the cosmic pie. In terms of the Einstein equations,Gμν=χTμν, the right hand side should include something more than the usual matter and radiation components in the stress-energy tensor.As a radically different approach, one can also try to leave unchanged the source side (actually “observed” since composed by radiation and baryonic matter), but rather modifying the left hand side. In a sense, one is therefore interpreting cosmic speed up as a first signal of the breakdown of the laws of physics as described by the standard General Relativity (GR). Since this theory has been experimentally tested only up to the Solar System scale, there is no a priori theoretical motivation to extend its validity to extraordinarily larger scales such as the extragalactic and cosmological ones (up to the last scattering surface!). Extending GR, not giving up to its positive results at local scales, opens the way to a large class of alternative theories of gravity ranging from extra-dimensions [34–38] to nonminimally coupled scalar fields [39–42]. In particular, we are interested here in fourth order theories [43–54] based on replacing the scalar curvature R in the Hilbert-Einstein action with a generic analytic function f(R) which should be reconstructed starting from data and physically motivated issues. Also referred to as f(R)-gravity, some of these models have been shown to be able to both fit the cosmological data and evade the Solar System constraints in several physically interesting cases [55–59].In this review paper, we will face two of the main problems directly related to the dark energy and dark matter issues: cosmography and clusters of galaxies. These are typical examples where the standard General Relativity and Newtonian potential schemes fail to describe dynamics since data present accelerated expansion and missing matter. Our goal is to address them byf(R)-gravity.


### 1.1. Cosmography: Why?

It is worth noting that both dark energy models and modified gravity theories seem to be in agreement with data. As a consequence, unless higher precision probes of the expansion rate and the growth of structure will be available, these two rival approaches could not be discriminated. This confusion about the theoretical background suggests that a more conservative approach to the problem of cosmic acceleration, relying on as less model dependent quantities as possible, is welcome. A possible solution could be to come back to the cosmography [60] rather than finding out solutions of the Friedmann equations and testing them. Being only related to the derivatives of the scale factor, the cosmographic parameters make it possible to fit the data on the distance-redshift relation without any a priori assumption on the underlying cosmological model: in this case, the only assumption is that the metric is the Robertson-Walker one (and hence not relying on the solution of cosmological equations). Almost eighty years after Hubble discovery of the expansion of the Universe, we can now extend, in principle, cosmography well beyond the search for the value of the only Hubble constant. The SNeIa Hubble diagram extends up to z=1.7 thus invoking the need for, at least, a fifth order Taylor expansion of the scale factor in order to give a reliable approximation of the distance-redshift relation. As a consequence, it could be, in principle, possible to estimate up to five cosmographic parameters, although the still too small data set available does not allow to get a precise and realistic determination of all of them.Once these quantities have been determined, one could use them to put constraints on the models. In a sense, we can revert the usual approach, consisting in deriving the cosmographic parameters as a sort of byproduct of an assumed theory. Here, we follow the other way around expressing the model characterizing quantities as a function of the cosmographic parameters. Such a program is particularly suited for the study of fourth order theories of gravity. As it is well known, the mathematical difficulties entering the solution of fourth order field equations make it quite problematic to find out analytical expressions for the scale factor and hence predict the values of the cosmographic parameters. A key role inf(R)-gravity is played by the choice of the f(R) function. Under quite general hypotheses, we will derive useful relations among the cosmographic parameters and the present day value of f(n)(R)=dnf/dRn, with n=0,…,3, whatever f(R) is. (As an important remark, we stress that our derivation will rely on the metric formulation of f(R) theories, while we refer the reader to [61, 62] for a similar work in the Palatini approach.) Once the cosmographic parameters will be determined, this method will allow us to investigate the cosmography of f(R) theories.It is worth stressing that the definition of the cosmographic parameters only relies on the assumption of the Robertson-Walker metric. As such, it is however difficult to state a priori to what extent the fifth order expansion provides an accurate enough description of the quantities of interest. Actually, the number of cosmographic parameters to be used depends on the problem one is interested in. As we will see later, we are here concerned only with the SNeIa Hubble diagram so that we have to check that the distance modulusμcp(z) obtained using the fifth order expansion of the scale factor is the same (within the errors) as the one μDE(z) of the underlying physical model. Being such a model of course unknown, one can adopt a phenomenological parameterization for the dark energy equation of state (EoS) and look at the percentage deviation Δμ/μDE as function of the EoS parameters. (Note that one can always use a phenomenological dark energy model to get a reliable estimate of the scale factor evolution even if the correct model is a fourth order one.) We have carried out such exercise using the CPL model, introduced below, and verified that Δμ/μDE is an increasing function of z (as expected), but still remains smaller than 2% up to z~2 over a wide range of the CPL parameter space. On the other hand, halting the Taylor expansion to a lower order may introduce significant deviation for z>1 that can potentially bias the analysis if the measurement errors are as small as those predicted by future SNeIa surveys. We are therefore confident that our fifth order expansion is both sufficient to get an accurate distance modulus over the redshift range probed by SNeIa and necessary to avoid dangerous biases.


### 1.2. Clusters of Galaxies: Why?

In the second part of this review we will apply thef(R)-gravity approach to cluster of galaxies. In fact, changing the gravity sector has consequences not only at cosmological scales, but also at galactic and cluster scales so that it is mandatory to investigate the low energy limit of such theories. A strong debate is open with different results arguing in favor [63–68] or against [69–71] such models at local scales. It is worth noting that, as a general result, higher order theories of gravity cause the gravitational potential to deviate from its Newtonian 1/r scaling [72–77] even if such deviations may be vanishing.In [78], the Newtonian limit of power law f(R)=f0Rn theories has been investigated, assuming that the metric in the low energy limit (Φ/c2≪1) may be taken as Schwarzschild-like. It turns out that a power law term (r/rc)β has to be added to the Newtonian 1/r term in order to get the correct gravitational potential. While the parameter β may be expressed analytically as a function of the slope n of the f(R) theory, rc sets the scale where the correction term starts being significant. A particular range of values of n has been investigated so that the corrective term is an increasing function of the radius r thus causing an increase of the rotation curve with respect to the Newtonian one and offering the possibility to fit the galaxy rotation curves without the need of further dark matter components.A set of low surface brightness (LSB) galaxies with extended and well measured rotation curves has been considered [79, 80]. These systems are supposed to be dark matter dominated, and successfully fitting data without dark matter is a strong evidence in favor of the approach (see also [81] for an independent analysis using another sample of galaxies). Combined with the hints coming from the cosmological applications, one should have, in principle, the possibility to address both the dark energy and dark matter problems resorting to the same well motivated fundamental theory [82–85]. Nevertheless, the simple power law f(R) gravity is nothing else but a toy-model which fail if one tries to achieve a comprehensive model for all the cosmological dynamics, ranging from the early Universe, to the large scale structure up to the late accelerated era [83, 84].A fundamental issue is related to clusters and superclusters of galaxies. Such structures, essentially, rule the large scale structure, and are the intermediate step between galaxies and cosmology. As the galaxies, they appear dark-matter dominated but the distribution of dark matter component seems clustered and organized in a very different way with respect to galaxies. It seems that dark matter is ruled by the scale and also its fundamental nature could depend on the scale. For a comprehensive review see [86].In the philosophy off(R)-gravity, the issue is to reconstruct the mass profile of clusters without dark matter, that is, to find out corrections to the Newton potential producing the same dynamics as dark matter but starting from a well motivated theory.In conclusion,f(R)-gravity, as the simplest approach to any extended or alternative gravity scheme, could be the paradigm to interpret dark energy and dark matter as curvature effects acting at scales larger than those where General Relativity has been actually investigated and probed.Let us discuss now how cosmography and then galaxy clusters could be two main examples to realize this program.


## 1.1. Cosmography: Why?

It is worth noting that both dark energy models and modified gravity theories seem to be in agreement with data. As a consequence, unless higher precision probes of the expansion rate and the growth of structure will be available, these two rival approaches could not be discriminated. This confusion about the theoretical background suggests that a more conservative approach to the problem of cosmic acceleration, relying on as less model dependent quantities as possible, is welcome. A possible solution could be to come back to the cosmography [60] rather than finding out solutions of the Friedmann equations and testing them. Being only related to the derivatives of the scale factor, the cosmographic parameters make it possible to fit the data on the distance-redshift relation without any a priori assumption on the underlying cosmological model: in this case, the only assumption is that the metric is the Robertson-Walker one (and hence not relying on the solution of cosmological equations). Almost eighty years after Hubble discovery of the expansion of the Universe, we can now extend, in principle, cosmography well beyond the search for the value of the only Hubble constant. The SNeIa Hubble diagram extends up to z=1.7 thus invoking the need for, at least, a fifth order Taylor expansion of the scale factor in order to give a reliable approximation of the distance-redshift relation. As a consequence, it could be, in principle, possible to estimate up to five cosmographic parameters, although the still too small data set available does not allow to get a precise and realistic determination of all of them.Once these quantities have been determined, one could use them to put constraints on the models. In a sense, we can revert the usual approach, consisting in deriving the cosmographic parameters as a sort of byproduct of an assumed theory. Here, we follow the other way around expressing the model characterizing quantities as a function of the cosmographic parameters. Such a program is particularly suited for the study of fourth order theories of gravity. As it is well known, the mathematical difficulties entering the solution of fourth order field equations make it quite problematic to find out analytical expressions for the scale factor and hence predict the values of the cosmographic parameters. A key role inf(R)-gravity is played by the choice of the f(R) function. Under quite general hypotheses, we will derive useful relations among the cosmographic parameters and the present day value of f(n)(R)=dnf/dRn, with n=0,…,3, whatever f(R) is. (As an important remark, we stress that our derivation will rely on the metric formulation of f(R) theories, while we refer the reader to [61, 62] for a similar work in the Palatini approach.) Once the cosmographic parameters will be determined, this method will allow us to investigate the cosmography of f(R) theories.It is worth stressing that the definition of the cosmographic parameters only relies on the assumption of the Robertson-Walker metric. As such, it is however difficult to state a priori to what extent the fifth order expansion provides an accurate enough description of the quantities of interest. Actually, the number of cosmographic parameters to be used depends on the problem one is interested in. As we will see later, we are here concerned only with the SNeIa Hubble diagram so that we have to check that the distance modulusμcp(z) obtained using the fifth order expansion of the scale factor is the same (within the errors) as the one μDE(z) of the underlying physical model. Being such a model of course unknown, one can adopt a phenomenological parameterization for the dark energy equation of state (EoS) and look at the percentage deviation Δμ/μDE as function of the EoS parameters. (Note that one can always use a phenomenological dark energy model to get a reliable estimate of the scale factor evolution even if the correct model is a fourth order one.) We have carried out such exercise using the CPL model, introduced below, and verified that Δμ/μDE is an increasing function of z (as expected), but still remains smaller than 2% up to z~2 over a wide range of the CPL parameter space. On the other hand, halting the Taylor expansion to a lower order may introduce significant deviation for z>1 that can potentially bias the analysis if the measurement errors are as small as those predicted by future SNeIa surveys. We are therefore confident that our fifth order expansion is both sufficient to get an accurate distance modulus over the redshift range probed by SNeIa and necessary to avoid dangerous biases.


## 1.2. Clusters of Galaxies: Why?

In the second part of this review we will apply thef(R)-gravity approach to cluster of galaxies. In fact, changing the gravity sector has consequences not only at cosmological scales, but also at galactic and cluster scales so that it is mandatory to investigate the low energy limit of such theories. A strong debate is open with different results arguing in favor [63–68] or against [69–71] such models at local scales. It is worth noting that, as a general result, higher order theories of gravity cause the gravitational potential to deviate from its Newtonian 1/r scaling [72–77] even if such deviations may be vanishing.In [78], the Newtonian limit of power law f(R)=f0Rn theories has been investigated, assuming that the metric in the low energy limit (Φ/c2≪1) may be taken as Schwarzschild-like. It turns out that a power law term (r/rc)β has to be added to the Newtonian 1/r term in order to get the correct gravitational potential. While the parameter β may be expressed analytically as a function of the slope n of the f(R) theory, rc sets the scale where the correction term starts being significant. A particular range of values of n has been investigated so that the corrective term is an increasing function of the radius r thus causing an increase of the rotation curve with respect to the Newtonian one and offering the possibility to fit the galaxy rotation curves without the need of further dark matter components.A set of low surface brightness (LSB) galaxies with extended and well measured rotation curves has been considered [79, 80]. These systems are supposed to be dark matter dominated, and successfully fitting data without dark matter is a strong evidence in favor of the approach (see also [81] for an independent analysis using another sample of galaxies). Combined with the hints coming from the cosmological applications, one should have, in principle, the possibility to address both the dark energy and dark matter problems resorting to the same well motivated fundamental theory [82–85]. Nevertheless, the simple power law f(R) gravity is nothing else but a toy-model which fail if one tries to achieve a comprehensive model for all the cosmological dynamics, ranging from the early Universe, to the large scale structure up to the late accelerated era [83, 84].A fundamental issue is related to clusters and superclusters of galaxies. Such structures, essentially, rule the large scale structure, and are the intermediate step between galaxies and cosmology. As the galaxies, they appear dark-matter dominated but the distribution of dark matter component seems clustered and organized in a very different way with respect to galaxies. It seems that dark matter is ruled by the scale and also its fundamental nature could depend on the scale. For a comprehensive review see [86].In the philosophy off(R)-gravity, the issue is to reconstruct the mass profile of clusters without dark matter, that is, to find out corrections to the Newton potential producing the same dynamics as dark matter but starting from a well motivated theory.In conclusion,f(R)-gravity, as the simplest approach to any extended or alternative gravity scheme, could be the paradigm to interpret dark energy and dark matter as curvature effects acting at scales larger than those where General Relativity has been actually investigated and probed.Let us discuss now how cosmography and then galaxy clusters could be two main examples to realize this program.


## 2. The Cosmographic Apparatus

The key rule in cosmography is the Taylor series expansion of the scale factor with respect to the cosmic time. To this aim, it is convenient to introduce the following functions:(1)H(t)≡+1adadt,q(t)≡-1ad2adt21H2,j(t)≡+1ad3adt31H3,s(t)≡+1ad4adt41H4,l(t)≡+1ad5adt51H5,
which are usually referred to as the Hubble, deceleration, jerk, snap and lerk parameters, respectively. It is then a matter of algebra to demonstrate the following useful relations: (2)Ḣ=-H2(1+q),Ḧ=H3(j+3q+2),d3Hdt3=H4[s-4j-3q(q+4)-6],d4Hdt4=H5[l-5s+10(q+2)j+30(q+2)q+24],
where a dot denotes derivative with respect to the cosmic time t. Equation (2) make it possible to relate the derivative of the Hubble parameter to the other cosmographic parameters. The distance-redshift relation may then be obtained starting from the Taylor expansion of a(t) along the lines described in [87–89].


### 2.1. The Scale-Factor Series

With these definitions the series expansion to the 5th order in time of the scale factor will be(3)a(t)=a(t0){H0(t-t0)-q02H02(t-t0)2+j03!H03(t-t0)3+s04!H04(t-t0)4+l05!H05(t-t0)5+O[(t-t0)6]},(4)a(t)a(t0)=1+H0(t-t0)-q02H02(t-t0)2+j03!H03(t-t0)3+s04!H04(t-t0)4+l05!H05(t-t0)5+O[(t-t0)6].
It is easy to see that (4) is the inverse of redshift z, being the redshift defined by (5)1+z=a(t0)a(t).
The physical distance travelled by a photon that is emitted at time t* and absorbed at the current epoch t0 is (6)D=c∫dt=c(t0-t*)
Assuming t*=t0-(D/c) and inserting in (4) we have: (7)1+z=a(t0)a(t0-D/c)=11-(H0/c)D-(q0/2)(H0/c)2D2-(j0/6)(H0/c)3D3+(s0/24)(H0/c)4D4-(l0/120)(H0/c)5D5+O[(H0D/c)6].The inverse of this expression will be(8)1+z=1+H0cD+(1+q02)(H0c)2D2+(1+q0+j06)(H0c)3D3+(1+32q0+q024+j03-s024)(H0c)4D4+(1+2q0+34q02+q0j06+j02-s12+l0)(H0c)5D5+O[(H0Dc)6].
Then we reverse the series z(D)→D(z) to have the physical distance D expressed as function of redshift z: (9)z(D)=𝒵D1(H0Dc)+𝒵D2(H0Dc)2+𝒵D3(H0Dc)3+𝒵D4(H0Dc)4+𝒵D5(H0Dc)5+O[(H0Dc)6]
with (10)𝒵D1=1,𝒵D2=1+q02,𝒵D3=1+q0+j06,𝒵D4=1+32q0+q024+j03-s024,𝒵D5=1+2q0+34q02+q0j06+j02-s12+l0.
From this we have(11)D(z)=czH0{𝒟z0+𝒟z1z+𝒟z2z2+𝒟z3z3+𝒟z4z4+O(z5)}
with (12)𝒟z0=1,𝒟z1=-(1+q02),𝒟z2=1+q0+q022-j06,𝒟z3=-(1+32q0+32q02+58q03-12j0-512q0j0-s024),𝒟z4=1+2q0+3q02+52q03+72q04-53q0j0-78q02j0,-18q0s0-j0+j0212-s06-l0120.
In typical applications, one is not interested in the physical distance D(z), but other definitions:  (i)
the luminosity distance:(13)dL=a(t0)a(t0-D/c)(a(t0)r0),(ii)
the angular-diameter distance:(14)dA=a(t0-D/c)a(t0)(a(t0)r0),wherer0(D) is(15)r0(D)={sin(∫t0-D/ct0cdta(t))k=+1;∫t0-D/ct0cdta(t)k=0;sinh(∫t0-D/ct0cdta(t))k=-1.
If we make the expansion for short distances, namely if we insert the series expansion of a(t) in r0(D), we have (16)r0(D)=∫t0-D/ct0cdta(t)=∫t0-D/ct0cdta0×{1+H0(t0-t)+(1+q02)H02(t0-t)2+(1+q0+j06)H03(t0-t)3+(1+32q0+q024+j03-s024)H04(t0-t)4+(1+2q0+34q02+q0j06+j02-s12+l0)×H05(t0-t)5+O[(t0-t)6](1+2q0+34q02+q0j06+j02-s12+l0)}=Da0{1+12H0Dc+[2+q06](H0Dc)2+[6+6q0+j024](H0Dc)3+[24+36q0+6q02+8j0-s0120](H0Dc)4+[12+24q0+9q02+2q0j0+6j0-s0+12l072]×(H0Dc)5+O[(H0Dc)6]}.
To convert from physical distance travelled to r coordinate traversed we have to consider that the Taylor series expansion of sin-sinh functions is(17)r0(D)=[∫t0-D/ct0cdta(t)]-k3![∫t0-D/ct0cdta(t)]3+O([∫t0-D/ct0cdta(t)]5)
so that (4) with curvature k term becomes (18)r0(D)=Da0{ℛD0+ℛD1H0Dc+ℛD2(H0Dc)2+ℛD3(H0Dc)3+ℛD4(H0Dc)4+ℛD5(H0Dc)5+O[(H0Dc)6]}
with(19)ℛD0=1,ℛD1=12,ℛD2=16[2+q0-kc2H02a02],ℛD3=124[6+6q0+j0-6kc2H02a02],ℛD4=1120[24+36q0+6q02+8j0-s0-5kc2(7+2q0)a02H02],ℛD5=1144[3kc2(15+10q0+j0)a02H0224+48q0+18q02+4q0j0+12j0-2s0+24l0-3kc2(15+10q0+j0)a02H02].
Using these one for luminosity distance we have(20)dL(z)=czH0{𝒟L0+𝒟L1z+𝒟L2z2+𝒟L3z3+𝒟L4z4+O(z5)}
with: (21)𝒟L0=1,𝒟L1=-12(-1+q0),𝒟L2=-16(1-q0-3q02+j0+kc2H02a02),𝒟L3=124(2kc2(1+3q0)H02a022-2q0-15q02-15q03+5j0+10q0j0+s0+2kc2(1+3q0)H02a02),𝒟L4=1120[5kc2(1+8q0+9q02-2j0)a02H02-6+6q0+81q02+165q03+105q04-110q0j0-105q02j0-15q0s0-27j0+10j2-11s0-l0-5kc2(1+8q0+9q02-2j0)a02H02].
While for the angular diameter distance it is
(22)dA(z)=czH0{𝒟A0+𝒟A1z+𝒟A2z2+𝒟A3z3+𝒟A4z4+O(z5)}
with(23)𝒟A0=1,𝒟A1=-12(3+q0),𝒟A2=16[11+7q0+3q02-j0-kc2H02a02],𝒟A3=-124[2kc2(5+3q0)H02a0250+46q0+39q02+15q03-13j0-10q0j0-s0-2kc2(5+3q0)H02a02],𝒟A4=1120[5kc2(17+20q0+9q02-2j0)a02H02274+326q0+411q02+315q03+105q04-210q0j0-105q02j0-15q0s0-137j0+10j2-21s0-l0-5kc2(17+20q0+9q02-2j0)a02H02].If we want to use the same notation of [87], we define Ω0=1+kc2/H02a02, which can be considered a purely cosmographic parameter, or Ω0=1-Ωk=Ωm,0+Ωr,0+ΩX,0 if we consider the dynamics of the Universe. With this parameter (12)–(14) become(24)𝒟L,y0=1,𝒟L,y1=-12(-3+q0),𝒟L,y2=-16(12-5q0+3q02-j0-Ω0),𝒟L,y3=124[52-20q0+21q02-15q03-7j0+10q0j0+s0-2Ω0(1+3q0)],𝒟L,y4=1120[359-184q0+186q02-135q03+105q04+90q0j0-105q02j0-15q0s0-57j0+10j2+9s0-l0-5Ω0(17-6q0+9q02-2j0)],𝒟A,y0=1,𝒟A,y1=-12(1+q0),𝒟A,y2=-16[-q0-3q02+j0+Ω0],𝒟A,y3=-124[-2q0+3q02+15q03-j0-10q0j0-s0+2Ω0],𝒟A,y4=-1120[1-6q0+9q02-15q03-105q04+10q0j0+105q02j0+15q0s0-3j0-10j2+s0+l0+5Ω0].Previous relations in this section have been derived for any value of the curvature parameter; but since in the following we will assume a flat Universe, we will used the simplified versions fork=0. Now, since we are going to use supernovae data, it will be useful to give as well the Taylor series of the expansion of the luminosity distance at it enters the modulus distance, which is the quantity about which those observational data inform. The final expression for the modulus distance based on the Hubble free luminosity distance, μ(z)=5log10dL(z), is(25)μ(z)=5log10·(logz+ℳ1z+ℳ2z2+ℳ3z3+ℳ4z4),
with (26)ℳ1=-12[-1+q0],ℳ2=-124[7-10q0-9q02+4j0],ℳ3=124[5-9q0-16q02-10q03+7j0+8q0j0+s0],ℳ4=12880[-469+1004q0+2654q02+3300q03+1575q04+200j02-1148j0-2620q0j0-1800q02j0-300q0s0-324s0-24l0].


## 2.1. The Scale-Factor Series

With these definitions the series expansion to the 5th order in time of the scale factor will be(3)a(t)=a(t0){H0(t-t0)-q02H02(t-t0)2+j03!H03(t-t0)3+s04!H04(t-t0)4+l05!H05(t-t0)5+O[(t-t0)6]},(4)a(t)a(t0)=1+H0(t-t0)-q02H02(t-t0)2+j03!H03(t-t0)3+s04!H04(t-t0)4+l05!H05(t-t0)5+O[(t-t0)6].
It is easy to see that (4) is the inverse of redshift z, being the redshift defined by (5)1+z=a(t0)a(t).
The physical distance travelled by a photon that is emitted at time t* and absorbed at the current epoch t0 is (6)D=c∫dt=c(t0-t*)
Assuming t*=t0-(D/c) and inserting in (4) we have: (7)1+z=a(t0)a(t0-D/c)=11-(H0/c)D-(q0/2)(H0/c)2D2-(j0/6)(H0/c)3D3+(s0/24)(H0/c)4D4-(l0/120)(H0/c)5D5+O[(H0D/c)6].The inverse of this expression will be(8)1+z=1+H0cD+(1+q02)(H0c)2D2+(1+q0+j06)(H0c)3D3+(1+32q0+q024+j03-s024)(H0c)4D4+(1+2q0+34q02+q0j06+j02-s12+l0)(H0c)5D5+O[(H0Dc)6].
Then we reverse the series z(D)→D(z) to have the physical distance D expressed as function of redshift z: (9)z(D)=𝒵D1(H0Dc)+𝒵D2(H0Dc)2+𝒵D3(H0Dc)3+𝒵D4(H0Dc)4+𝒵D5(H0Dc)5+O[(H0Dc)6]
with (10)𝒵D1=1,𝒵D2=1+q02,𝒵D3=1+q0+j06,𝒵D4=1+32q0+q024+j03-s024,𝒵D5=1+2q0+34q02+q0j06+j02-s12+l0.
From this we have(11)D(z)=czH0{𝒟z0+𝒟z1z+𝒟z2z2+𝒟z3z3+𝒟z4z4+O(z5)}
with (12)𝒟z0=1,𝒟z1=-(1+q02),𝒟z2=1+q0+q022-j06,𝒟z3=-(1+32q0+32q02+58q03-12j0-512q0j0-s024),𝒟z4=1+2q0+3q02+52q03+72q04-53q0j0-78q02j0,-18q0s0-j0+j0212-s06-l0120.
In typical applications, one is not interested in the physical distance D(z), but other definitions:  (i)
the luminosity distance:(13)dL=a(t0)a(t0-D/c)(a(t0)r0),(ii)
the angular-diameter distance:(14)dA=a(t0-D/c)a(t0)(a(t0)r0),wherer0(D) is(15)r0(D)={sin(∫t0-D/ct0cdta(t))k=+1;∫t0-D/ct0cdta(t)k=0;sinh(∫t0-D/ct0cdta(t))k=-1.
If we make the expansion for short distances, namely if we insert the series expansion of a(t) in r0(D), we have (16)r0(D)=∫t0-D/ct0cdta(t)=∫t0-D/ct0cdta0×{1+H0(t0-t)+(1+q02)H02(t0-t)2+(1+q0+j06)H03(t0-t)3+(1+32q0+q024+j03-s024)H04(t0-t)4+(1+2q0+34q02+q0j06+j02-s12+l0)×H05(t0-t)5+O[(t0-t)6](1+2q0+34q02+q0j06+j02-s12+l0)}=Da0{1+12H0Dc+[2+q06](H0Dc)2+[6+6q0+j024](H0Dc)3+[24+36q0+6q02+8j0-s0120](H0Dc)4+[12+24q0+9q02+2q0j0+6j0-s0+12l072]×(H0Dc)5+O[(H0Dc)6]}.
To convert from physical distance travelled to r coordinate traversed we have to consider that the Taylor series expansion of sin-sinh functions is(17)r0(D)=[∫t0-D/ct0cdta(t)]-k3![∫t0-D/ct0cdta(t)]3+O([∫t0-D/ct0cdta(t)]5)
so that (4) with curvature k term becomes (18)r0(D)=Da0{ℛD0+ℛD1H0Dc+ℛD2(H0Dc)2+ℛD3(H0Dc)3+ℛD4(H0Dc)4+ℛD5(H0Dc)5+O[(H0Dc)6]}
with(19)ℛD0=1,ℛD1=12,ℛD2=16[2+q0-kc2H02a02],ℛD3=124[6+6q0+j0-6kc2H02a02],ℛD4=1120[24+36q0+6q02+8j0-s0-5kc2(7+2q0)a02H02],ℛD5=1144[3kc2(15+10q0+j0)a02H0224+48q0+18q02+4q0j0+12j0-2s0+24l0-3kc2(15+10q0+j0)a02H02].
Using these one for luminosity distance we have(20)dL(z)=czH0{𝒟L0+𝒟L1z+𝒟L2z2+𝒟L3z3+𝒟L4z4+O(z5)}
with: (21)𝒟L0=1,𝒟L1=-12(-1+q0),𝒟L2=-16(1-q0-3q02+j0+kc2H02a02),𝒟L3=124(2kc2(1+3q0)H02a022-2q0-15q02-15q03+5j0+10q0j0+s0+2kc2(1+3q0)H02a02),𝒟L4=1120[5kc2(1+8q0+9q02-2j0)a02H02-6+6q0+81q02+165q03+105q04-110q0j0-105q02j0-15q0s0-27j0+10j2-11s0-l0-5kc2(1+8q0+9q02-2j0)a02H02].
While for the angular diameter distance it is
(22)dA(z)=czH0{𝒟A0+𝒟A1z+𝒟A2z2+𝒟A3z3+𝒟A4z4+O(z5)}
with(23)𝒟A0=1,𝒟A1=-12(3+q0),𝒟A2=16[11+7q0+3q02-j0-kc2H02a02],𝒟A3=-124[2kc2(5+3q0)H02a0250+46q0+39q02+15q03-13j0-10q0j0-s0-2kc2(5+3q0)H02a02],𝒟A4=1120[5kc2(17+20q0+9q02-2j0)a02H02274+326q0+411q02+315q03+105q04-210q0j0-105q02j0-15q0s0-137j0+10j2-21s0-l0-5kc2(17+20q0+9q02-2j0)a02H02].If we want to use the same notation of [87], we define Ω0=1+kc2/H02a02, which can be considered a purely cosmographic parameter, or Ω0=1-Ωk=Ωm,0+Ωr,0+ΩX,0 if we consider the dynamics of the Universe. With this parameter (12)–(14) become(24)𝒟L,y0=1,𝒟L,y1=-12(-3+q0),𝒟L,y2=-16(12-5q0+3q02-j0-Ω0),𝒟L,y3=124[52-20q0+21q02-15q03-7j0+10q0j0+s0-2Ω0(1+3q0)],𝒟L,y4=1120[359-184q0+186q02-135q03+105q04+90q0j0-105q02j0-15q0s0-57j0+10j2+9s0-l0-5Ω0(17-6q0+9q02-2j0)],𝒟A,y0=1,𝒟A,y1=-12(1+q0),𝒟A,y2=-16[-q0-3q02+j0+Ω0],𝒟A,y3=-124[-2q0+3q02+15q03-j0-10q0j0-s0+2Ω0],𝒟A,y4=-1120[1-6q0+9q02-15q03-105q04+10q0j0+105q02j0+15q0s0-3j0-10j2+s0+l0+5Ω0].Previous relations in this section have been derived for any value of the curvature parameter; but since in the following we will assume a flat Universe, we will used the simplified versions fork=0. Now, since we are going to use supernovae data, it will be useful to give as well the Taylor series of the expansion of the luminosity distance at it enters the modulus distance, which is the quantity about which those observational data inform. The final expression for the modulus distance based on the Hubble free luminosity distance, μ(z)=5log10dL(z), is(25)μ(z)=5log10·(logz+ℳ1z+ℳ2z2+ℳ3z3+ℳ4z4),
with (26)ℳ1=-12[-1+q0],ℳ2=-124[7-10q0-9q02+4j0],ℳ3=124[5-9q0-16q02-10q03+7j0+8q0j0+s0],ℳ4=12880[-469+1004q0+2654q02+3300q03+1575q04+200j02-1148j0-2620q0j0-1800q02j0-300q0s0-324s0-24l0].


## 3.f(R)-Gravity versus Cosmography

### 3.1.f(R) Preliminaries

As discussed in the introduction, much interest has been recently devoted to the possibility that dark energy could be nothing else but a curvature effect according to which the present Universe is filled by pressureless dust matter only and the acceleration is the result of modified Friedmann equations obtained by replacing the Ricci curvature scalarR with a generic function f(R) in the gravity action. Under the assumption of a flat Universe, the Hubble parameter is therefore determined by (we use here natural units such that 8πG=1):(27)H2=13[ρmf'(R)+ρcurv],
where the prime denotes derivative with respect to R and ρcurv is the energy density of an effective curvature fluid(note that the name curvature fluid does not refer to the FRW curvature parameter k, but only takes into account that such a term is a geometrical one related to the scalar curvature R):(28)ρcurv=1f'(R){12[f(R)-Rf'(R)]-3HṘf′′(R)}.
Assuming there is no interaction between the matter and the curvature terms (we are in the so-called Jordan frame), the matter continuity equation gives the usual scaling ρM=ρM(t=t0)a-3=3H02ΩMa-3, with ΩM the present day matter density parameter. The continuity equation for ρcurv then reads:(29)ρ̇curv+3H(1+wcurv)ρcurv=3H02ΩMṘf′′(R)[f′(R)]2a-3
with(30)wcurv=-1+R̈f′′(R)+Ṙ[Ṙf′′′(R)-Hf′′(R)][f(R)-Rf′(R)]/2-3HṘf′′(R)
the barotropic factor of the curvature fluid. It is worth noticing that the curvature fluid quantities ρcurv and wcurv only depends on f(R) and its derivatives up to the third order. As a consequence, considering only their present day values (which may be naively obtained by replacing R with R0 everywhere), two f(R) theories sharing the same values of f(R0), f'(R0), f′′(R0), f′′′(R0) will be degenerate from this point of view. (One can argue that this is not strictly true since different f(R) theories will lead to different expansion rate H(t) and hence different present day values of R and its derivatives. However, it is likely that two f(R) functions that exactly match each other up to the third order derivative today will give rise to the same H(t) at least for t≃t0 so that (R0,Ṙ0,R̈0) will be almost the same.)Combining (29) with (27), one finally gets the following master equation for the Hubble parameter:(31)Ḣ=-12f'(R){3H02ΩMa-3+R̈f′′(R)+Ṙ[Ṙf′′′(R)-Hf′′(R)]}.
Expressing the scalar curvature R as function of the Hubble parameter as:(32)R=-6(Ḣ+2H2)
and inserting the result into (31), one ends with a fourth order nonlinear differential equation for the scale factor a(t) that cannot be easily solved also for the simplest cases (for instance, f(R)∝Rn). Moreover, although technically feasible, a numerical solution of (31) is plagued by the large uncertainties on the boundary conditions (i.e., the present day values of the scale factor and its derivatives up to the third order) that have to be set to find out the scale factor.


### 3.2.f(R)-Derivatives and Cosmography

Motivated by these difficulties, we approach now the problem from a different viewpoint. Rather than choosing a parameterized expression forf(R) and then numerically solving (31) for given values of the boundary conditions, we try to relate the present day values of its derivatives to the cosmographic parameters (q0,j0,s0,l0) so that constraining them in a model independent way gives us a hint for what kind of f(R) theory could be able to fit the observed Hubble diagram. (Note that a similar analysis, but in the context of the energy conditions in f(R), has yet been presented in [90]. However, in that work, the author give an expression for f(R) and then compute the snap parameter to be compared to the observed one. On the contrary, our analysis does not depend on any assumed functional expression for f(R).)As a preliminary step, it is worth considering again the constraint equation (32). Differentiating with respect to t, we easily get the following relations(33)Ṙ=-6(Ḧ+4HḢ),R̈=-6(d3Hdt3+4HḦ+4Ḣ2),d3Rdt3R=-6(d4Hdt4+4Hd3Hdt3+12ḢḦ).
Evaluating these at the present time and using (2), one finally gets(34)R0=-6H02(1-q0),Ṙ0=-6H03(j0-q0-2),R̈0=-6H04(s0+q02+8q0+6),d3R0dt3=-6H05[l0-s0+2(q0+4)j0-6(3q0+8)q0-24],
which will turn out to be useful in the following.Let us now come back to the expansion rate and master equations (27) and (31). Since they have to hold along the full evolutionary history of the Universe, they naively hold also at the present day. As a consequence, we may evaluate them in t=t0 thus easily obtaining(35)H02=H02ΩMf'(R0)+f(R0)-R0f'(R0)-6H0Ṙ0f′′(R0)6f'(R0)-Ḣ0=3H02ΩM2f'(R0)+Ṙ02f′′′(R0)+(R̈0-H0Ṙ0)f′′(R0)2f'(R0).
Using (2) and (34), we can rearrange (35) as two relations among the Hubble constant H0 and the cosmographic parameters (q0,j0,s0), on one hand, and the present day values of f(R) and its derivatives up to third order. However, two further relations are needed in order to close the system and determine the four unknown quantities f(R0), f'(R0), f′′(R0), f′′′(R0). A first one may be easily obtained by noting that, inserting back the physical units, the rate expansion equation reads(36)H2=8πG3f'(R)[ρm+ρcurvf'(R)]
which clearly shows that, in f(R) gravity, the Newtonian gravitational constant G is replaced by an effective (time dependent) Geff=G/f'(R). On the other hand, it is reasonable to assume that the present day value of Geff is the same as the Newtonian one so that we get the simple constraint:(37)Geff(z=0)=G→f'(R0)=1.
In order to get the fourth relation we need to close the system, we first differentiate both sides of (31) with respect to t. We thus get(38)Ḧ=Ṙ2f′′′(R)+(R̈-HṘ)f′′(R)+3H02ΩMa-32[Ṙf′′′(R)]-1[f'(R)]2×Ṙ3f(iv)(R)+(3ṘR̈-HṘ2)f′′′(R)2f′(R)×(d3R/dt3-HṘ+ḢṘ)f′′(R)-9H02ΩMHa-32f′(R),
with f(iv)(R)=d4f/dR4. Let us now suppose that f(R) may be well approximated by its third order Taylor expansion in R-R0, that is, we set(39)f(R)=f(R0)+f′(R0)(R-R0)+12f′′(R0)(R-R0)2+16f′′′(R0)(R-R0)3.
In such an approximation, it is f(n)(R)=dnf/Rn=0 for n≥4 so that naively f(iv)(R0)=0. Evaluating then (38) at the present day, we get(40)Ḧ0=Ṙ2f′′′(R0)+(R̈0-H0Ṙ0)f′′(R0)+3H02ΩM2[Ṙ0f′′(R0)]-1[f'(R0)]2-(3Ṙ0R̈0-HṘ02)f′′′(R0)2f′(R0)-(d3R0/dt3-H0R̈0+Ḣ0Ṙ0)f′′(R0)-9H03ΩM2f′(R0).
We can now schematically proceed as follows. Evaluate (2) at z=0 and plug these relations into the left hand sides of (35), (40). Insert (34) into the right hand sides of these same equations so that only the cosmographic parameters (q0,j0,s0,l0) and the f(R) related quantities enter both sides of these relations. Finally, solve them under the constraint (37) with respect to the present day values of f(R) and its derivatives up to the third order. After some algebra, one ends up with the desired result:(41)f(R0)6H02=-𝒫0(q0,j0,s0,l0)ΩM+𝒬0(q0,j0,s0,l0)ℛ(q0,j0,s0,l0),(42)f'(R0)=1,(43)f′′(R0)(6H02)-1=-𝒫2(q0,j0,s0)ΩM+𝒬2(q0,j0,s0)ℛ(q0,j0,s0,l0),(44)f′′′(R0)(6H02)-2=-𝒫3(q0,j0,s0,l0)ΩM+𝒬3(q0,j0,s0,l0)(j0-q0-2)ℛ(q0,j0,s0,l0),
where we have defined(45)𝒫0=(j0-q0-2)l0-(3s0+7j0+6q02+41q0+22)s0-[(3q0+16)j0+20q02+64q0+12]j0-(3q04+25q03+96q02+72q0+20),(46)𝒬0=(q02-j0q0+2q0)l0+[3q0s0+(4q0+6)j0+6q03+44q02+22q0-12]s0+[2j02+(3q02+10q0-6)j0+17q03+52q02+54q0+36]j0+3q05+28q04+118q03+72q02-76q0-64,(47)𝒫2=9s0+6j0+9q02+66q0+42,(48)𝒬2=-{6(q0+1)s0+[2j0-2(1-q0)]j0+6q03+50q02+74q0+32},(49)𝒫3=3l0+3s0-9(q0+4)j0-(45q02+78q0+12),(50)𝒬3=-{2(1+q0)l0+2(j0+q0)s0-(2j0+4q02+12q0+6)j0-(30q03+84q02+78q0+24)},(51)ℛ=(j0-q0-2)l0-(3s0-2j0+6q02+50q0+40)s0+[(3q0+10)j0+11q02+4q0-18]j0-(3q04+34q03+246q0+104).
Equations (41)–(51) make it possible to estimate the present day values of f(R) and its first three derivatives as function of the Hubble constant H0 and the cosmographic parameters (q0,j0,s0,l0) provided a value for the matter density parameter ΩM is given. This is a somewhat problematic point. Indeed, while the cosmographic parameters may be estimated in a model independent way, the fiducial value for ΩM is usually the outcome of fitting a given dataset in the framework of an assumed dark energy scenario. However, it is worth noting that different models all converge towards the concordance value ΩM≃0.25 which is also in agreement with astrophysical (model independent) estimates from the gas mass fraction in galaxy clusters. On the other hand, it has been proposed that f(R) theories may avoid the need for dark matter in galaxies and galaxy clusters [74, 76, 78, 81–83, 91]. In such a case, the total matter content of the Universe is essentially equal to the baryonic one. According to the primordial elements abundance and the standard BBN scenario, we therefore get ΩM≃ωb/h2 with ωb=Ωbh2≃0.0214 [92] and h the Hubble constant in units of 100km/s/Mpc. Setting h=0.72 in agreement with the results of the HST Key project [93], we thus get ΩM=0.041 for a baryons only Universe. We will therefore consider in the following both cases when numerical estimates are needed.It is worth noticing thatH0 only plays the role of a scaling parameter giving the correct physical dimensions to f(R) and its derivatives. As such, it is not surprising that we need four cosmographic parameters, namely (q0,j0,s0,l0), to fix the four f(R) related quantities f(R0), f'(R0), f′′(R0), f′′′(R0). It is also worth stressing that (41)–(44) are linear in the f(R) quantities so that (q0,j0,s0,l0) uniquely determine the former ones. On the contrary, inverting them to get the cosmographic parameters as function of the f(R) ones, we do not get linear relations. Indeed, the field equations in f(R) theories are nonlinear fourth order differential equations in the scale factor a(t) so that fixing the derivatives of f(R) up to third order makes it possible to find out a class of solutions, not a single one. Each one of these solutions will be characterized by a different set of cosmographic parameters thus explaining why the inversion of (41)–(51) does not give a unique result for (q0,j0,s0,l0).As a final comment, we reconsider the underlying assumptions leading to the above derived relations. While (35) are exact relations deriving from a rigorous application of the field equations, (40) heavily relies on having approximated f(R) with its third order Taylor expansion (39). If this assumption fails, the system should not be closed since a fifth unknown parameter enters the game, namely f(iv)(R0). Actually, replacing f(R) with its Taylor expansion is not possible for all class of f(R) theories. As such, the above results only hold in those cases where such an expansion is possible. Moreover, by truncating the expansion to the third order, we are implicitly assuming that higher order terms are negligible over the redshift range probed by the data. That is to say, we are assuming that(52)f(n)(R0)(R-R0)n≪∑m=03f(m)(R0)m!(R-R0)mforn≥4
over the redshift range probed by the data. Checking the validity of this assumption is not possible without explicitly solving the field equations, but we can guess an order of magnitude estimate considering that, for all viable models, the background dynamics should not differ too much from the ΛCDM one at least up to z≃2. Using then the expression of H(z) for the ΛCDM model, it is easily to see that R/R0 is a quickly increasing function of the redshift so that, in order (52) holds, we have to assume that f(n)(R0)≪f′′′(R0) for n≥4. This condition is easier to check for many analytical f(R) models.Once such a relation is verified, we have still to worry about (37) relying on the assumption that the cosmological gravitational constant is exactly the same as the local one. Although reasonable, this requirement is not absolutely demonstrated. Actually, the numerical value usually adopted for the Newton constant GN is obtained from laboratory experiments in settings that can hardly be considered homogenous and isotropic. As such, the spacetime metric in such conditions has nothing to do with the cosmological one so that matching the two values of G is strictly speaking an extrapolation. Although commonly accepted and quite reasonable, the condition Glocal=Gcosmo could (at least, in principle) be violated so that (37) could be reconsidered. Indeed, as we will see, the condition f'(R0)=1 may not be verified for some popular f(R) models recently proposed in literature. However, it is reasonable to assume that Geff(z=0)=G(1+ɛ) with ɛ≪1. When this be the case, we should repeat the derivation of (41)–(44) now using the condition f'(R0)=(1+ɛ)-1. Taylor expanding the results in ɛ to the first order and comparing with the above derived equations, we can estimate the error induced by our assumption ɛ=0. The resulting expressions are too lengthy to be reported and depend in a complicated way on the values of the matter density parameter ΩM, the cosmographic parameters (q0,j0,s0,l0) and ɛ. However, we have numerically checked that the error induced on f(R0), f′′(R0), f′′′(R0) are much lower than 10% for value of ɛ as high as an unrealistic ɛ~0.1. We are confident that our results are reliable also for these cases.


## 3.1.f(R) Preliminaries

As discussed in the introduction, much interest has been recently devoted to the possibility that dark energy could be nothing else but a curvature effect according to which the present Universe is filled by pressureless dust matter only and the acceleration is the result of modified Friedmann equations obtained by replacing the Ricci curvature scalarR with a generic function f(R) in the gravity action. Under the assumption of a flat Universe, the Hubble parameter is therefore determined by (we use here natural units such that 8πG=1):(27)H2=13[ρmf'(R)+ρcurv],
where the prime denotes derivative with respect to R and ρcurv is the energy density of an effective curvature fluid(note that the name curvature fluid does not refer to the FRW curvature parameter k, but only takes into account that such a term is a geometrical one related to the scalar curvature R):(28)ρcurv=1f'(R){12[f(R)-Rf'(R)]-3HṘf′′(R)}.
Assuming there is no interaction between the matter and the curvature terms (we are in the so-called Jordan frame), the matter continuity equation gives the usual scaling ρM=ρM(t=t0)a-3=3H02ΩMa-3, with ΩM the present day matter density parameter. The continuity equation for ρcurv then reads:(29)ρ̇curv+3H(1+wcurv)ρcurv=3H02ΩMṘf′′(R)[f′(R)]2a-3
with(30)wcurv=-1+R̈f′′(R)+Ṙ[Ṙf′′′(R)-Hf′′(R)][f(R)-Rf′(R)]/2-3HṘf′′(R)
the barotropic factor of the curvature fluid. It is worth noticing that the curvature fluid quantities ρcurv and wcurv only depends on f(R) and its derivatives up to the third order. As a consequence, considering only their present day values (which may be naively obtained by replacing R with R0 everywhere), two f(R) theories sharing the same values of f(R0), f'(R0), f′′(R0), f′′′(R0) will be degenerate from this point of view. (One can argue that this is not strictly true since different f(R) theories will lead to different expansion rate H(t) and hence different present day values of R and its derivatives. However, it is likely that two f(R) functions that exactly match each other up to the third order derivative today will give rise to the same H(t) at least for t≃t0 so that (R0,Ṙ0,R̈0) will be almost the same.)Combining (29) with (27), one finally gets the following master equation for the Hubble parameter:(31)Ḣ=-12f'(R){3H02ΩMa-3+R̈f′′(R)+Ṙ[Ṙf′′′(R)-Hf′′(R)]}.
Expressing the scalar curvature R as function of the Hubble parameter as:(32)R=-6(Ḣ+2H2)
and inserting the result into (31), one ends with a fourth order nonlinear differential equation for the scale factor a(t) that cannot be easily solved also for the simplest cases (for instance, f(R)∝Rn). Moreover, although technically feasible, a numerical solution of (31) is plagued by the large uncertainties on the boundary conditions (i.e., the present day values of the scale factor and its derivatives up to the third order) that have to be set to find out the scale factor.


## 3.2.f(R)-Derivatives and Cosmography

Motivated by these difficulties, we approach now the problem from a different viewpoint. Rather than choosing a parameterized expression forf(R) and then numerically solving (31) for given values of the boundary conditions, we try to relate the present day values of its derivatives to the cosmographic parameters (q0,j0,s0,l0) so that constraining them in a model independent way gives us a hint for what kind of f(R) theory could be able to fit the observed Hubble diagram. (Note that a similar analysis, but in the context of the energy conditions in f(R), has yet been presented in [90]. However, in that work, the author give an expression for f(R) and then compute the snap parameter to be compared to the observed one. On the contrary, our analysis does not depend on any assumed functional expression for f(R).)As a preliminary step, it is worth considering again the constraint equation (32). Differentiating with respect to t, we easily get the following relations(33)Ṙ=-6(Ḧ+4HḢ),R̈=-6(d3Hdt3+4HḦ+4Ḣ2),d3Rdt3R=-6(d4Hdt4+4Hd3Hdt3+12ḢḦ).
Evaluating these at the present time and using (2), one finally gets(34)R0=-6H02(1-q0),Ṙ0=-6H03(j0-q0-2),R̈0=-6H04(s0+q02+8q0+6),d3R0dt3=-6H05[l0-s0+2(q0+4)j0-6(3q0+8)q0-24],
which will turn out to be useful in the following.Let us now come back to the expansion rate and master equations (27) and (31). Since they have to hold along the full evolutionary history of the Universe, they naively hold also at the present day. As a consequence, we may evaluate them in t=t0 thus easily obtaining(35)H02=H02ΩMf'(R0)+f(R0)-R0f'(R0)-6H0Ṙ0f′′(R0)6f'(R0)-Ḣ0=3H02ΩM2f'(R0)+Ṙ02f′′′(R0)+(R̈0-H0Ṙ0)f′′(R0)2f'(R0).
Using (2) and (34), we can rearrange (35) as two relations among the Hubble constant H0 and the cosmographic parameters (q0,j0,s0), on one hand, and the present day values of f(R) and its derivatives up to third order. However, two further relations are needed in order to close the system and determine the four unknown quantities f(R0), f'(R0), f′′(R0), f′′′(R0). A first one may be easily obtained by noting that, inserting back the physical units, the rate expansion equation reads(36)H2=8πG3f'(R)[ρm+ρcurvf'(R)]
which clearly shows that, in f(R) gravity, the Newtonian gravitational constant G is replaced by an effective (time dependent) Geff=G/f'(R). On the other hand, it is reasonable to assume that the present day value of Geff is the same as the Newtonian one so that we get the simple constraint:(37)Geff(z=0)=G→f'(R0)=1.
In order to get the fourth relation we need to close the system, we first differentiate both sides of (31) with respect to t. We thus get(38)Ḧ=Ṙ2f′′′(R)+(R̈-HṘ)f′′(R)+3H02ΩMa-32[Ṙf′′′(R)]-1[f'(R)]2×Ṙ3f(iv)(R)+(3ṘR̈-HṘ2)f′′′(R)2f′(R)×(d3R/dt3-HṘ+ḢṘ)f′′(R)-9H02ΩMHa-32f′(R),
with f(iv)(R)=d4f/dR4. Let us now suppose that f(R) may be well approximated by its third order Taylor expansion in R-R0, that is, we set(39)f(R)=f(R0)+f′(R0)(R-R0)+12f′′(R0)(R-R0)2+16f′′′(R0)(R-R0)3.
In such an approximation, it is f(n)(R)=dnf/Rn=0 for n≥4 so that naively f(iv)(R0)=0. Evaluating then (38) at the present day, we get(40)Ḧ0=Ṙ2f′′′(R0)+(R̈0-H0Ṙ0)f′′(R0)+3H02ΩM2[Ṙ0f′′(R0)]-1[f'(R0)]2-(3Ṙ0R̈0-HṘ02)f′′′(R0)2f′(R0)-(d3R0/dt3-H0R̈0+Ḣ0Ṙ0)f′′(R0)-9H03ΩM2f′(R0).
We can now schematically proceed as follows. Evaluate (2) at z=0 and plug these relations into the left hand sides of (35), (40). Insert (34) into the right hand sides of these same equations so that only the cosmographic parameters (q0,j0,s0,l0) and the f(R) related quantities enter both sides of these relations. Finally, solve them under the constraint (37) with respect to the present day values of f(R) and its derivatives up to the third order. After some algebra, one ends up with the desired result:(41)f(R0)6H02=-𝒫0(q0,j0,s0,l0)ΩM+𝒬0(q0,j0,s0,l0)ℛ(q0,j0,s0,l0),(42)f'(R0)=1,(43)f′′(R0)(6H02)-1=-𝒫2(q0,j0,s0)ΩM+𝒬2(q0,j0,s0)ℛ(q0,j0,s0,l0),(44)f′′′(R0)(6H02)-2=-𝒫3(q0,j0,s0,l0)ΩM+𝒬3(q0,j0,s0,l0)(j0-q0-2)ℛ(q0,j0,s0,l0),
where we have defined(45)𝒫0=(j0-q0-2)l0-(3s0+7j0+6q02+41q0+22)s0-[(3q0+16)j0+20q02+64q0+12]j0-(3q04+25q03+96q02+72q0+20),(46)𝒬0=(q02-j0q0+2q0)l0+[3q0s0+(4q0+6)j0+6q03+44q02+22q0-12]s0+[2j02+(3q02+10q0-6)j0+17q03+52q02+54q0+36]j0+3q05+28q04+118q03+72q02-76q0-64,(47)𝒫2=9s0+6j0+9q02+66q0+42,(48)𝒬2=-{6(q0+1)s0+[2j0-2(1-q0)]j0+6q03+50q02+74q0+32},(49)𝒫3=3l0+3s0-9(q0+4)j0-(45q02+78q0+12),(50)𝒬3=-{2(1+q0)l0+2(j0+q0)s0-(2j0+4q02+12q0+6)j0-(30q03+84q02+78q0+24)},(51)ℛ=(j0-q0-2)l0-(3s0-2j0+6q02+50q0+40)s0+[(3q0+10)j0+11q02+4q0-18]j0-(3q04+34q03+246q0+104).
Equations (41)–(51) make it possible to estimate the present day values of f(R) and its first three derivatives as function of the Hubble constant H0 and the cosmographic parameters (q0,j0,s0,l0) provided a value for the matter density parameter ΩM is given. This is a somewhat problematic point. Indeed, while the cosmographic parameters may be estimated in a model independent way, the fiducial value for ΩM is usually the outcome of fitting a given dataset in the framework of an assumed dark energy scenario. However, it is worth noting that different models all converge towards the concordance value ΩM≃0.25 which is also in agreement with astrophysical (model independent) estimates from the gas mass fraction in galaxy clusters. On the other hand, it has been proposed that f(R) theories may avoid the need for dark matter in galaxies and galaxy clusters [74, 76, 78, 81–83, 91]. In such a case, the total matter content of the Universe is essentially equal to the baryonic one. According to the primordial elements abundance and the standard BBN scenario, we therefore get ΩM≃ωb/h2 with ωb=Ωbh2≃0.0214 [92] and h the Hubble constant in units of 100km/s/Mpc. Setting h=0.72 in agreement with the results of the HST Key project [93], we thus get ΩM=0.041 for a baryons only Universe. We will therefore consider in the following both cases when numerical estimates are needed.It is worth noticing thatH0 only plays the role of a scaling parameter giving the correct physical dimensions to f(R) and its derivatives. As such, it is not surprising that we need four cosmographic parameters, namely (q0,j0,s0,l0), to fix the four f(R) related quantities f(R0), f'(R0), f′′(R0), f′′′(R0). It is also worth stressing that (41)–(44) are linear in the f(R) quantities so that (q0,j0,s0,l0) uniquely determine the former ones. On the contrary, inverting them to get the cosmographic parameters as function of the f(R) ones, we do not get linear relations. Indeed, the field equations in f(R) theories are nonlinear fourth order differential equations in the scale factor a(t) so that fixing the derivatives of f(R) up to third order makes it possible to find out a class of solutions, not a single one. Each one of these solutions will be characterized by a different set of cosmographic parameters thus explaining why the inversion of (41)–(51) does not give a unique result for (q0,j0,s0,l0).As a final comment, we reconsider the underlying assumptions leading to the above derived relations. While (35) are exact relations deriving from a rigorous application of the field equations, (40) heavily relies on having approximated f(R) with its third order Taylor expansion (39). If this assumption fails, the system should not be closed since a fifth unknown parameter enters the game, namely f(iv)(R0). Actually, replacing f(R) with its Taylor expansion is not possible for all class of f(R) theories. As such, the above results only hold in those cases where such an expansion is possible. Moreover, by truncating the expansion to the third order, we are implicitly assuming that higher order terms are negligible over the redshift range probed by the data. That is to say, we are assuming that(52)f(n)(R0)(R-R0)n≪∑m=03f(m)(R0)m!(R-R0)mforn≥4
over the redshift range probed by the data. Checking the validity of this assumption is not possible without explicitly solving the field equations, but we can guess an order of magnitude estimate considering that, for all viable models, the background dynamics should not differ too much from the ΛCDM one at least up to z≃2. Using then the expression of H(z) for the ΛCDM model, it is easily to see that R/R0 is a quickly increasing function of the redshift so that, in order (52) holds, we have to assume that f(n)(R0)≪f′′′(R0) for n≥4. This condition is easier to check for many analytical f(R) models.Once such a relation is verified, we have still to worry about (37) relying on the assumption that the cosmological gravitational constant is exactly the same as the local one. Although reasonable, this requirement is not absolutely demonstrated. Actually, the numerical value usually adopted for the Newton constant GN is obtained from laboratory experiments in settings that can hardly be considered homogenous and isotropic. As such, the spacetime metric in such conditions has nothing to do with the cosmological one so that matching the two values of G is strictly speaking an extrapolation. Although commonly accepted and quite reasonable, the condition Glocal=Gcosmo could (at least, in principle) be violated so that (37) could be reconsidered. Indeed, as we will see, the condition f'(R0)=1 may not be verified for some popular f(R) models recently proposed in literature. However, it is reasonable to assume that Geff(z=0)=G(1+ɛ) with ɛ≪1. When this be the case, we should repeat the derivation of (41)–(44) now using the condition f'(R0)=(1+ɛ)-1. Taylor expanding the results in ɛ to the first order and comparing with the above derived equations, we can estimate the error induced by our assumption ɛ=0. The resulting expressions are too lengthy to be reported and depend in a complicated way on the values of the matter density parameter ΩM, the cosmographic parameters (q0,j0,s0,l0) and ɛ. However, we have numerically checked that the error induced on f(R0), f′′(R0), f′′′(R0) are much lower than 10% for value of ɛ as high as an unrealistic ɛ~0.1. We are confident that our results are reliable also for these cases.


## 4.f(R)-Gravity and the CPL Model

A determination off(R) and its derivatives in terms of the cosmographic parameters need for an estimate of these latter from the data in a model independent way. Unfortunately, even in the nowadays era of precision cosmology, such a program is still too ambitious to give useful constraints on the f(R) derivatives, as we will see later. On the other hand, the cosmographic parameters may also be expressed in terms of the dark energy density and EoS parameters so that we can work out what are the present day values of f(R) and its derivatives giving the same (q0,j0,s0,l0) of the given dark energy model. To this aim, it is convenient to adopt a parameterized expression for the dark energy EoS in order to reduce the dependence of the results on any underlying theoretical scenario. Following the prescription of the Dark Energy Task Force [94], we will use the Chevallier-Polarski-Linder (CPL) parameterization for the EoS setting [95, 96]:(53)w=w0+wa(1-a)=w0+waz(1+z)-1
so that, in a flat Universe filled by dust matter and dark energy, the dimensionless Hubble parameter E(z)=H/H0 reads(54)E2(z)=ΩM(1+z)3+ΩX(1+z)3(1+w0+wa)e-3waz/(1+z)
with ΩX=1-ΩM because of the flatness assumption. In order to determine the cosmographic parameters for such a model, we avoid integrating H(z) to get a(t) by noting that d/dt=-(1+z)H(z)d/dz. We can use such a relation to evaluate (Ḣ,Ḧ,d3H/dt3,d4H/dt4) and then solve (2), evaluated in z=0, with respect to the parameters of interest. Some algebra finally gives(55)q0=12+32(1-ΩM)w0,(56)j0=1+32(1-ΩM)[3w0(1+w0)+wa],(57)s0=-72-334(1-ΩM)wa-94(1-ΩM)×[9+(7-ΩM)wa]w0-94(1-ΩM)×(16-3ΩM)w02-274(1-ΩM)(3-ΩM)w03(58)l0=352+1-ΩM4[213+(7-ΩM)wa]wa+1-ΩM4[489+9(82-21ΩM)wa]w0+92(1-ΩM)[67-21ΩM+32(23-11ΩM)wa]w02+274(1-ΩM)(47-24ΩM)w03+812(1-ΩM)(3-2ΩM)w04.
Inserting (55)–(58) into (41)–(51), we get lengthy expressions (which we do not report here) giving the present day values of f(R) and its first three derivatives as function of (ΩM,w0,wa). It is worth noting that the f(R) model thus obtained is not dynamically equivalent to the starting CPL one. Indeed, the two models have the same cosmographic parameters only today. As such, for instance, the scale factor is the same between the two theories only over the time period during which the fifth order Taylor expansion is a good approximation of the actual a(t). It is also worth stressing that such a procedure does not select a unique f(R) model, but rather a class of fourth order theories all sharing the same third order Taylor expansion of f(R).


### 4.1. TheΛCDM Case

With these caveats in mind, it is worth considering first theΛCDM model which is obtained by setting (w0,wa)=(-1,0) in the above expressions thus giving(59)q0=12-32(1-ΩM),j0=1,s0=1-92ΩM,l0=1+3ΩM+272ΩM2.
When inserted into the expressions for the f(R) quantities, these relations give the remarkable result:(60)f(R0)=R0+2Λ,f′′(R0)=f′′′(R0)=0
so that we obviously conclude that the only f(R) theory having exactly the same cosmographic parameters as the ΛCDM model is just f(R)∝R, that is, GR. It is worth noticing that such a result comes out as a consequence of the values of (q0,j0) in the ΛCDM model. Indeed, should we have left (s0,l0) undetermined and only fixed (q0,j0) to the values in (59), we should have got the same result in (60). Since the ΛCDM model fits well a large set of different data, we do expect that the actual values of (q0,j0,s0,l0) do not differ too much from the ΛCDM ones. Therefore, we plug into (41)–(51) the following expressions:(61)q0=q0Λ×(1+ɛq),j0=j0Λ×(1+ɛj),s0=s0Λ×(1+ɛs),l0=l0Λ×(1+ɛl),
with (q0Λ,j0Λ,s0Λ,l0Λ) given by (59) and (ɛq,ɛj,ɛs,ɛl) quantifying the deviations from the ΛCDM values allowed by the data. A numerical estimate of these quantities may be obtained, for example, from a Marko chain analysis, but this is outside our aims. Since we are here interested in a theoretical examination, we prefer to consider an idealized situation where the four quantities above all share the same value ɛ≪1. In such a case, we can easily investigate how much the corresponding f(R) deviates from the GR one considering the two ratios f′′(R0)/f(R0) and f′′′(R0)/f(R0). Inserting the above expressions for the cosmographic parameters into the exact (not reported) formulae for f(R0), f′′(R0) and f′′′(R0), taking their ratios and then expanding to first order in ɛ, we finally get(62)η20=64-6ΩM(9ΩM+8)[3(9ΩM+74)ΩM-556]ΩM2+16×ɛ27,η30=6[(81ΩM-110)ΩM+40]ΩM+16[3(9ΩM+74)ΩM-556]ΩM2+16×ɛ243ΩM2,
having defined η20=f′′(R0)/f(R0)×H04 and η30=f′′′(R0)/f(R0)×H06 which, being dimensionless quantities, are more suited to estimate the order of magnitudes of the different terms. Inserting our fiducial values for ΩM, we get:(63)η20≃0.15×ɛforΩM=0.041,η20≃-0.12×ɛforΩM=0.250,η30≃4×ɛforΩM=0.041,η30≃-0.18×ɛforΩM=0.250.
For values of ɛ up to 0.1, the above relations show that the second and third derivatives are at most two orders of magnitude smaller than the zeroth order term f(R0). Actually, the values of η30 for a baryon only model (first row) seems to argue in favor of a larger importance of the third order term. However, we have numerically checked that the above relations approximates very well the exact expressions up to ɛ≃0.1 with an accuracy depending on the value of ΩM, being smaller for smaller matter density parameters. Using the exact expressions for η20 and η30, our conclusion on the negligible effect of the second and third order derivatives are significantly strengthened.Such a result holds under the hypotheses that the narrower are the constraints on the validity of theΛCDM model, the smaller are the deviations of the cosmographic parameters from the ΛCDM ones. It is possible to show that this indeed the case for the CPL parametrization we are considering. On the other hand, we have also assumed that the deviations (ɛq,ɛj,ɛs,ɛl) take the same values. Although such hypothesis is somewhat ad hoc, we argue that the main results are not affected by giving it away. Indeed, although different from each other, we can still assume that all of them are very small so that Taylor expanding to the first order should lead to additional terms into (62) which are likely of the same order of magnitude. We may therefore conclude that, if the observations confirm that the values of the cosmographic parameters agree within ~10% with those predicted for the ΛCDM model, we must conclude that the deviations of f(R) from the GR case, f(R)∝R, should be vanishingly small.It is worth stressing, however, that such a conclusion only holds for thosef(R) models satisfying the constraint (52). It is indeed possible to work out a model having f(R0)∝R0, f′′(R0)=f′′′(R0)=0, but f(n)(R0)≠0 for some n. For such a (somewhat ad hoc) model, (52) is clearly not satisfied so that the cosmographic parameters have to be evaluated from the solution of the field equations. For such a model, the conclusion above does not hold so that one cannot exclude that the resulting (q0,j0,s0,l0) are within 10% of the ΛCDM ones.


### 4.2. The Constant EoS Model

Let us now take into account the conditionw=-1, but still retains wa=0 thus obtaining the so called quiessence models. In such a case, some problems arise because both the terms (j0-q0-2) and ℛ may vanish for some combinations of the two model parameters (ΩM,w0). For instance, we find that j0-q0-2=0 for w0=(w1,w2) with:(64)w1=11-ΩM+(1-ΩM)(4-ΩM),w2=-13[1+4-ΩM(1-ΩM)(4-ΩM)].
On the other hand, the equation ℛ(ΩM,w0)=0 may have different real roots for w depending on the adopted value of ΩM. Denoting collectively with wnull the values of w0 that, for a given ΩM, make (j0-q0-2)ℛ(ΩM,w0) taking the null value, we individuate a set of quiessence models whose cosmographic parameters give rise to divergent values of f(R0, f′′(R0) and f′′′(R0). For such models, f(R) is clearly not defined so that we have to exclude these cases from further consideration. We only note that it is still possible to work out a f(R) theory reproducing the same background dynamics of such models, but a different route has to be used.Since bothq0 and j0 now deviate from the ΛCDM values, it is not surprising that both f′′(R0) and f′′′(R0) take finite non null values. However, it is more interesting to study the two quantities η20 and η30 defined above to investigate the deviations of f(R) from the GR case. These are plotted in Figures 1 and 2 for the two fiducial ΩM values. Note that the range of w0 in these plots have been chosen in order to avoid divergences, but the lessons we will draw also hold for the other w0 values.Figure 1
The dimensionless ratio between the present day values off′′(R) and f(R) as function of the constant EoS w0 of the corresponding quiessence model. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Figure 2
The dimensionless ratio between the present day values off′′′(R) and f(R) as function of the constant EoS w0 of the corresponding quiessence model. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.As a general comment, it is clear that, even in this case,f′′(R0) and f′′′(R0) are from two to three orders of magnitude smaller that the zeroth order term f(R0). Such a result could be yet guessed from the previous discussion for the ΛCDM case. Actually, relaxing the hypothesis w0=-1 is the same as allowing the cosmographic parameters to deviate from the ΛCDM values. Although a direct mapping between the two cases cannot be established, it is nonetheless evident that such a relation can be argued thus making the outcome of the above plots not fully surprising. It is nevertheless worth noting that, while in the ΛCDM case, η20 and η30 always have opposite signs, this is not the case for quiessence models with w>-1. Indeed, depending on the value of ΩM, we can have f(R) theories with both η20 and η30 positive. Moreover, the lower is ΩM, the higher are the ratios η20 and η30 for a given value of w0. This can be explained qualitatively noticing that, for a lower ΩM, the density parameter of the curvature fluid (playing the role of an effective dark energy) must be larger thus claiming for higher values of the second and third derivatives (see also [97] for a different approach to the problem).


### 4.3. The General Case

Finally, we consider evolving dark energy models withwa≠0. Needless to say, varying three parameters allows to get a wide range of models that cannot be discussed in detail. Therefore, we only concentrate on evolving dark energy models with w0=-1 in agreement with some most recent analysis. The results on η20 and η30 are plotted in Figures 3 and 4 where these quantities as functions of wa. Note that we are considering models with positive wa so that w(z) tends to w0+wa>w0 for z→∞ so that the EoS dark energy can eventually approach the dust value w=0. Actually, this is also the range favored by the data. We have, however, excluded values where η20 or η30 diverge. Considering how they are defined, it is clear that these two quantities diverge when f(R0)=0 so that the values of (w0,wa) making (η20,η30) to diverge may be found solving(65)𝒫0(w0,wa)ΩM+𝒬0(w0,wa)=0,
where 𝒫0(w0,wa) and 𝒬0(w0,wa) are obtained by inserting (55)–(58) into the definitions (45)-(46). For such CPL models, there is no any f(R) model having the same cosmographic parameters and, at the same time, satisfying all the criteria needed for the validity of our procedure. Actually, if f(R0)=0, the condition (52) is likely to be violated so that higher than third order must be included in the Taylor expansion of f(R) thus invalidating the derivation of (41)–(44).Figure 3
The dimensionless ratio between the present day values off′′(R) and f(R) as function of the wa parameter for models with w0=-1. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Figure 4
The dimensionless ratio between the present day values off′′′(R) and f(R) as function of the wa parameter for models with w0=-1. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Under these caveats, Figures3 and 4 demonstrate that allowing the dark energy EoS to evolve does not change significantly our conclusions. Indeed, the second and third derivatives, although being not null, are nevertheless negligible with respect to the zeroth order term thus arguing in favour of a GR-like f(R) with only very small corrections. Such a result is, however, not fully unexpected. From (55) and (56), we see that, having set w0=-1, the q0 parameter is the same as for the ΛCDM model, while j0 reads j0Λ+(3/2)(1-ΩM)wa. As we have stressed above, the Hilbert-Einstein Lagrangian f(R)=R+2Λ is recovered when (q0,j0)=(q0Λ,j0Λ) whatever the values of (s0,l0) are. Introducing a wa≠0 makes (s0,l0) to differ from the ΛCDM values, but the first two cosmographic parameters are only mildly affected. Such deviations are then partially washed out by the complicated way they enter in the determination of the present day values of f(R) and its first three derivatives.


## 4.1. TheΛCDM Case

With these caveats in mind, it is worth considering first theΛCDM model which is obtained by setting (w0,wa)=(-1,0) in the above expressions thus giving(59)q0=12-32(1-ΩM),j0=1,s0=1-92ΩM,l0=1+3ΩM+272ΩM2.
When inserted into the expressions for the f(R) quantities, these relations give the remarkable result:(60)f(R0)=R0+2Λ,f′′(R0)=f′′′(R0)=0
so that we obviously conclude that the only f(R) theory having exactly the same cosmographic parameters as the ΛCDM model is just f(R)∝R, that is, GR. It is worth noticing that such a result comes out as a consequence of the values of (q0,j0) in the ΛCDM model. Indeed, should we have left (s0,l0) undetermined and only fixed (q0,j0) to the values in (59), we should have got the same result in (60). Since the ΛCDM model fits well a large set of different data, we do expect that the actual values of (q0,j0,s0,l0) do not differ too much from the ΛCDM ones. Therefore, we plug into (41)–(51) the following expressions:(61)q0=q0Λ×(1+ɛq),j0=j0Λ×(1+ɛj),s0=s0Λ×(1+ɛs),l0=l0Λ×(1+ɛl),
with (q0Λ,j0Λ,s0Λ,l0Λ) given by (59) and (ɛq,ɛj,ɛs,ɛl) quantifying the deviations from the ΛCDM values allowed by the data. A numerical estimate of these quantities may be obtained, for example, from a Marko chain analysis, but this is outside our aims. Since we are here interested in a theoretical examination, we prefer to consider an idealized situation where the four quantities above all share the same value ɛ≪1. In such a case, we can easily investigate how much the corresponding f(R) deviates from the GR one considering the two ratios f′′(R0)/f(R0) and f′′′(R0)/f(R0). Inserting the above expressions for the cosmographic parameters into the exact (not reported) formulae for f(R0), f′′(R0) and f′′′(R0), taking their ratios and then expanding to first order in ɛ, we finally get(62)η20=64-6ΩM(9ΩM+8)[3(9ΩM+74)ΩM-556]ΩM2+16×ɛ27,η30=6[(81ΩM-110)ΩM+40]ΩM+16[3(9ΩM+74)ΩM-556]ΩM2+16×ɛ243ΩM2,
having defined η20=f′′(R0)/f(R0)×H04 and η30=f′′′(R0)/f(R0)×H06 which, being dimensionless quantities, are more suited to estimate the order of magnitudes of the different terms. Inserting our fiducial values for ΩM, we get:(63)η20≃0.15×ɛforΩM=0.041,η20≃-0.12×ɛforΩM=0.250,η30≃4×ɛforΩM=0.041,η30≃-0.18×ɛforΩM=0.250.
For values of ɛ up to 0.1, the above relations show that the second and third derivatives are at most two orders of magnitude smaller than the zeroth order term f(R0). Actually, the values of η30 for a baryon only model (first row) seems to argue in favor of a larger importance of the third order term. However, we have numerically checked that the above relations approximates very well the exact expressions up to ɛ≃0.1 with an accuracy depending on the value of ΩM, being smaller for smaller matter density parameters. Using the exact expressions for η20 and η30, our conclusion on the negligible effect of the second and third order derivatives are significantly strengthened.Such a result holds under the hypotheses that the narrower are the constraints on the validity of theΛCDM model, the smaller are the deviations of the cosmographic parameters from the ΛCDM ones. It is possible to show that this indeed the case for the CPL parametrization we are considering. On the other hand, we have also assumed that the deviations (ɛq,ɛj,ɛs,ɛl) take the same values. Although such hypothesis is somewhat ad hoc, we argue that the main results are not affected by giving it away. Indeed, although different from each other, we can still assume that all of them are very small so that Taylor expanding to the first order should lead to additional terms into (62) which are likely of the same order of magnitude. We may therefore conclude that, if the observations confirm that the values of the cosmographic parameters agree within ~10% with those predicted for the ΛCDM model, we must conclude that the deviations of f(R) from the GR case, f(R)∝R, should be vanishingly small.It is worth stressing, however, that such a conclusion only holds for thosef(R) models satisfying the constraint (52). It is indeed possible to work out a model having f(R0)∝R0, f′′(R0)=f′′′(R0)=0, but f(n)(R0)≠0 for some n. For such a (somewhat ad hoc) model, (52) is clearly not satisfied so that the cosmographic parameters have to be evaluated from the solution of the field equations. For such a model, the conclusion above does not hold so that one cannot exclude that the resulting (q0,j0,s0,l0) are within 10% of the ΛCDM ones.


## 4.2. The Constant EoS Model

Let us now take into account the conditionw=-1, but still retains wa=0 thus obtaining the so called quiessence models. In such a case, some problems arise because both the terms (j0-q0-2) and ℛ may vanish for some combinations of the two model parameters (ΩM,w0). For instance, we find that j0-q0-2=0 for w0=(w1,w2) with:(64)w1=11-ΩM+(1-ΩM)(4-ΩM),w2=-13[1+4-ΩM(1-ΩM)(4-ΩM)].
On the other hand, the equation ℛ(ΩM,w0)=0 may have different real roots for w depending on the adopted value of ΩM. Denoting collectively with wnull the values of w0 that, for a given ΩM, make (j0-q0-2)ℛ(ΩM,w0) taking the null value, we individuate a set of quiessence models whose cosmographic parameters give rise to divergent values of f(R0, f′′(R0) and f′′′(R0). For such models, f(R) is clearly not defined so that we have to exclude these cases from further consideration. We only note that it is still possible to work out a f(R) theory reproducing the same background dynamics of such models, but a different route has to be used.Since bothq0 and j0 now deviate from the ΛCDM values, it is not surprising that both f′′(R0) and f′′′(R0) take finite non null values. However, it is more interesting to study the two quantities η20 and η30 defined above to investigate the deviations of f(R) from the GR case. These are plotted in Figures 1 and 2 for the two fiducial ΩM values. Note that the range of w0 in these plots have been chosen in order to avoid divergences, but the lessons we will draw also hold for the other w0 values.Figure 1
The dimensionless ratio between the present day values off′′(R) and f(R) as function of the constant EoS w0 of the corresponding quiessence model. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Figure 2
The dimensionless ratio between the present day values off′′′(R) and f(R) as function of the constant EoS w0 of the corresponding quiessence model. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.As a general comment, it is clear that, even in this case,f′′(R0) and f′′′(R0) are from two to three orders of magnitude smaller that the zeroth order term f(R0). Such a result could be yet guessed from the previous discussion for the ΛCDM case. Actually, relaxing the hypothesis w0=-1 is the same as allowing the cosmographic parameters to deviate from the ΛCDM values. Although a direct mapping between the two cases cannot be established, it is nonetheless evident that such a relation can be argued thus making the outcome of the above plots not fully surprising. It is nevertheless worth noting that, while in the ΛCDM case, η20 and η30 always have opposite signs, this is not the case for quiessence models with w>-1. Indeed, depending on the value of ΩM, we can have f(R) theories with both η20 and η30 positive. Moreover, the lower is ΩM, the higher are the ratios η20 and η30 for a given value of w0. This can be explained qualitatively noticing that, for a lower ΩM, the density parameter of the curvature fluid (playing the role of an effective dark energy) must be larger thus claiming for higher values of the second and third derivatives (see also [97] for a different approach to the problem).


## 4.3. The General Case

Finally, we consider evolving dark energy models withwa≠0. Needless to say, varying three parameters allows to get a wide range of models that cannot be discussed in detail. Therefore, we only concentrate on evolving dark energy models with w0=-1 in agreement with some most recent analysis. The results on η20 and η30 are plotted in Figures 3 and 4 where these quantities as functions of wa. Note that we are considering models with positive wa so that w(z) tends to w0+wa>w0 for z→∞ so that the EoS dark energy can eventually approach the dust value w=0. Actually, this is also the range favored by the data. We have, however, excluded values where η20 or η30 diverge. Considering how they are defined, it is clear that these two quantities diverge when f(R0)=0 so that the values of (w0,wa) making (η20,η30) to diverge may be found solving(65)𝒫0(w0,wa)ΩM+𝒬0(w0,wa)=0,
where 𝒫0(w0,wa) and 𝒬0(w0,wa) are obtained by inserting (55)–(58) into the definitions (45)-(46). For such CPL models, there is no any f(R) model having the same cosmographic parameters and, at the same time, satisfying all the criteria needed for the validity of our procedure. Actually, if f(R0)=0, the condition (52) is likely to be violated so that higher than third order must be included in the Taylor expansion of f(R) thus invalidating the derivation of (41)–(44).Figure 3
The dimensionless ratio between the present day values off′′(R) and f(R) as function of the wa parameter for models with w0=-1. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Figure 4
The dimensionless ratio between the present day values off′′′(R) and f(R) as function of the wa parameter for models with w0=-1. Short dashed and solid lines refer to models with ΩM=0.041 and 0.250, respectively.Under these caveats, Figures3 and 4 demonstrate that allowing the dark energy EoS to evolve does not change significantly our conclusions. Indeed, the second and third derivatives, although being not null, are nevertheless negligible with respect to the zeroth order term thus arguing in favour of a GR-like f(R) with only very small corrections. Such a result is, however, not fully unexpected. From (55) and (56), we see that, having set w0=-1, the q0 parameter is the same as for the ΛCDM model, while j0 reads j0Λ+(3/2)(1-ΩM)wa. As we have stressed above, the Hilbert-Einstein Lagrangian f(R)=R+2Λ is recovered when (q0,j0)=(q0Λ,j0Λ) whatever the values of (s0,l0) are. Introducing a wa≠0 makes (s0,l0) to differ from the ΛCDM values, but the first two cosmographic parameters are only mildly affected. Such deviations are then partially washed out by the complicated way they enter in the determination of the present day values of f(R) and its first three derivatives.


## 5. Constrainingf(R) Parameters

In the previous section, we have worked an alternative method to estimatef(R0), f′′(R0), f′′′(R0) resorting to a model independent parameterization of the dark energy EoS. However, in the ideal case, the cosmographic parameters are directly estimated from the data so that (41)–(51) can be used to infer the values of the f(R) related quantities. These latter can then be used to put constraints on the parameters entering an assumed fourth order theory assigned by a f(R) function characterized by a set of parameters p=(p1,…,pn) provided that the hypotheses underlying the derivation of (41)–(51) are indeed satisfied. We show below two interesting cases which clearly highlight the potentiality and the limitations of such an analysis.


### 5.1. Double Power Law Lagrangian

As a first interesting example, we set(66)f(R)=R(1+αRn+βR-m)
with n and m two positive real numbers (see, e.g., [98] for some physical motivations). The following expressions are immediately obtained:(67)f(R0)=R0(1+αR0n+βR0-m),f'(R0)=1+α(n+1)R0n-β(m-1)R0-m,f′′(R0)=αn(n+1)R0n-1+βm(m-1)R0-(1+m),f′′′(R0)=αn(n+1)(n-1)R0n-2,-βm(m+1)(m-1)R0-(2+m).
Denoting by ϕi (with i=0,…,3) the values of f(i)(R0) determined through (41)–(51), we can solve(68)f(R0)=ϕ0,f′(R0)=ϕ1,f′′(R0)=ϕ2,f′′′(R0)=ϕ3,
which is a system of four equations in the four unknowns (α,β,n,m) that can be analytically solved proceeding as follows. First, we solve the first and second equation with respect to (α,β) obtaining(69)α=1-mn+m(1-ϕ0R0)R0-n,β=-1+nn+m(1-ϕ0R0)R0m,
while, solving the third and fourth equations, we get(70)α=ϕ2R01-n[1+m+(ϕ3/ϕ2)R0]n(n+1)(n+m),β=ϕ2R01+n[1-n+(ϕ3/ϕ2)R0]m(1-m)(n+m).
Equating the two solutions, we get a systems of two equations in the two unknowns (n,m), namely,(71)n(n+1)(1-m)(1-ϕ0/R0)ϕ2R0[1+m+(ϕ3/ϕ2)R0]=1,m(n+1)(m-1)(1-ϕ0/R0)ϕ2R0[1-n+(ϕ3/ϕ2)R0]=1.
Solving with respect to m, we get two solutions, the first one being m=-n which has to be discarded since makes (α,β) goes to infinity. The only acceptable solution is(72)m=-[1-n+(ϕ3ϕ2)R0]
which, inserted back into the above system, leads to a second order polynomial equation for n with solutions(73)n=12[1+ϕ3ϕ2R0±𝒩(ϕ0,ϕ2,ϕ3)ϕ2R0(1+ϕ0/R0)],
where we have defined(74)𝒩(ϕ0,ϕ2,ϕ3)=(R02ϕ02-2R03ϕ0+R04)ϕ32+6(R0ϕ02-2R02ϕ0+R03)ϕ2ϕ3+9(ϕ02-2R0ϕ0+R02)ϕ22+4(R02ϕ0-R03)ϕ23.
Depending on the values of (q0,j0,s0,l0), (73) may lead to one, two or any acceptable solution, that is, real positive values of n. This solution has then to be inserted back into (72) to determine m and then into (69) or (70) to estimate (α,β). If the final values of (α,β,n,m) are physically viable, we can conclude that the model in (66) is in agreement with the data giving the same cosmographic parameters inferred from the data themselves. Exploring analytically what is the region of the (q0,j0,s0,l0) parameter space which leads to acceptable (α,β,n,m) solutions is a daunting task far outside the aim of the present work.


### 5.2. The Hu and Sawicki Model

One of the most pressing problems off(R) theories is the need to escape the severe constraints imposed by the Solar System tests. A successful model has been recently proposed by Hu and Sawicki [56] (HS) setting (note that such a model does not pass the matter instability test so that some viable generalizations [99–101] have been proposed):(75)f(R)=R-Rcα(R/Rc)n1+β(R/Rc)n.
As for the double power law model discussed above, there are four parameters which we can be expressed in terms of the cosmographic parameters (q0,j0,s0,l0).As a first step, it is trivial to get:(76)f(R0)=R0-RcαR0cn1+βR0cn,f′(R0)=1-αnRcR0cnR0(1+βR0cn)2,f′′(R0)=αnRcR0cn[(1-n)+β(1+n)R0cn]R02(1+βR0cn)3,f′′′(R0)=αnRcR0cn(An2+Bn+C)R03(1+βR0cn)4.
with R0c=R0/Rc and:(77)A=-β2R0c2n+4βR0cn-1,B=3(1-β2R0c2n),C=-2(1-βR0cn)2.
Equating (76) to the four quantities (ϕ0,ϕ1,ϕ2,ϕ3) defined as above, we could, in principle, solve this system of four equations in four unknowns to get (α,β,Rc,n) in terms of (ϕ0,ϕ1,ϕ2,ϕ3) and then, using (41)–(51) as functions of the cosmographic parameters. However, setting ϕ1=1 as required by (42) gives the only trivial solution αnRc=0 so that the HS model reduces to the Einstein-Hilbert Lagrangian f(R)=R. In order to escape this problem, we can relax the condition f'(R0)=1 to f'(R0)=(1+ɛ)-1. As we have discussed in Section 4, this is the same as assuming that the present day effective gravitational constant Geff,0=GN/f'(R0) only slightly differs from the usual Newtonian one which seems to be a quite reasonable assumption. Under this hypothesis, we can analytically solve for (α,β,Rc,n) in terms of (ϕ0,ɛ,ϕ2,ϕ3). The actual values of (ϕ0,ϕ2,ϕ3) will be no more given by (41)–(44), but we have checked that they deviate from those expressions (note that the correct expressions for (ϕ0,ϕ2,ϕ3) may still formally be written as (41)–(44), but the polynomials entering them are now different and also depend on powers of ɛ) much less than 10% for ɛ up to 10% well below any realistic expectation.With this caveat in mind, we first solve(78)f(R0)=ϕ0,f′′(R0)=(1+ɛ)-1
to get (79)α=n(1+ɛ)ɛ(R0Rc)1-n(1-ϕ0R0)2,β=n(1+ɛ)ɛ(R0Rc)-n[1-ϕ0R0-ɛn(1+ɛ)].
Inserting these expressions in (76), it is easy to check that Rc cancels out so that we can no more determine its value. Such a result is, however, not unexpected. Indeed, (75) can trivially be rewritten as(80)f(R)=R-α̃Rn1+β̃Rn
with α̃=αRc1-n and β̃=βRc-n which are indeed the quantities that are determined by the above expressions for (α,β). Reversing the discussion, the present day values of f(i)(R) depend on (α,β,Rc) only through the two parameters (α̃,β̃). As such, the use of cosmographic parameters is unable to break this degeneracy. However, since Rc only plays the role of a scaling parameter, we can arbitrarily set its value without loss of generality.On the other hand, this degeneracy allows us to get a consistency relation to immediately check whether the HS model is viable or not. Indeed, solving the equationf′′(R0)=ϕ2, we get(81)n=(ϕ0/R0)+[(1+ɛ)/ɛ](1-ϕ2R0)-(1-ɛ)/(1+ɛ)1-ϕ0/R0,
which can then be inserted into the equations f′′′(R0)=ϕ3 to obtain a complicated relation among (ϕ0,ϕ2,ϕ3) which we do not report for sake of shortness. Solving such a relation with respect to ϕ3/ϕ0 and Taylor expanding to first order in ɛ, the constraint we get reads(82)ϕ3ϕ0≃-1+ɛɛϕ2R0[R0(ϕ2ϕ0)+ɛϕ0-11+ɛ(1-2ɛ1-ϕ0/R0)].
If the cosmographic parameters (q0,j0,s0,l0) are known with sufficient accuracy, one could compute the values of (R0,ϕ0,ϕ2,ϕ3) for a given ɛ (eventually using the expressions obtained for ɛ=0) and then check if they satisfied this relation. If this is not the case, one can immediately give off the HS model also without the need of solving the field equations and fitting the data. Actually, given the still large errors on the cosmographic parameters, such a test only remains in the realm of (quite distant) future applications. However, the HS model works for other tests as shown in [56] and so a consistent cosmography analysis has to be combined with them.


## 5.1. Double Power Law Lagrangian

As a first interesting example, we set(66)f(R)=R(1+αRn+βR-m)
with n and m two positive real numbers (see, e.g., [98] for some physical motivations). The following expressions are immediately obtained:(67)f(R0)=R0(1+αR0n+βR0-m),f'(R0)=1+α(n+1)R0n-β(m-1)R0-m,f′′(R0)=αn(n+1)R0n-1+βm(m-1)R0-(1+m),f′′′(R0)=αn(n+1)(n-1)R0n-2,-βm(m+1)(m-1)R0-(2+m).
Denoting by ϕi (with i=0,…,3) the values of f(i)(R0) determined through (41)–(51), we can solve(68)f(R0)=ϕ0,f′(R0)=ϕ1,f′′(R0)=ϕ2,f′′′(R0)=ϕ3,
which is a system of four equations in the four unknowns (α,β,n,m) that can be analytically solved proceeding as follows. First, we solve the first and second equation with respect to (α,β) obtaining(69)α=1-mn+m(1-ϕ0R0)R0-n,β=-1+nn+m(1-ϕ0R0)R0m,
while, solving the third and fourth equations, we get(70)α=ϕ2R01-n[1+m+(ϕ3/ϕ2)R0]n(n+1)(n+m),β=ϕ2R01+n[1-n+(ϕ3/ϕ2)R0]m(1-m)(n+m).
Equating the two solutions, we get a systems of two equations in the two unknowns (n,m), namely,(71)n(n+1)(1-m)(1-ϕ0/R0)ϕ2R0[1+m+(ϕ3/ϕ2)R0]=1,m(n+1)(m-1)(1-ϕ0/R0)ϕ2R0[1-n+(ϕ3/ϕ2)R0]=1.
Solving with respect to m, we get two solutions, the first one being m=-n which has to be discarded since makes (α,β) goes to infinity. The only acceptable solution is(72)m=-[1-n+(ϕ3ϕ2)R0]
which, inserted back into the above system, leads to a second order polynomial equation for n with solutions(73)n=12[1+ϕ3ϕ2R0±𝒩(ϕ0,ϕ2,ϕ3)ϕ2R0(1+ϕ0/R0)],
where we have defined(74)𝒩(ϕ0,ϕ2,ϕ3)=(R02ϕ02-2R03ϕ0+R04)ϕ32+6(R0ϕ02-2R02ϕ0+R03)ϕ2ϕ3+9(ϕ02-2R0ϕ0+R02)ϕ22+4(R02ϕ0-R03)ϕ23.
Depending on the values of (q0,j0,s0,l0), (73) may lead to one, two or any acceptable solution, that is, real positive values of n. This solution has then to be inserted back into (72) to determine m and then into (69) or (70) to estimate (α,β). If the final values of (α,β,n,m) are physically viable, we can conclude that the model in (66) is in agreement with the data giving the same cosmographic parameters inferred from the data themselves. Exploring analytically what is the region of the (q0,j0,s0,l0) parameter space which leads to acceptable (α,β,n,m) solutions is a daunting task far outside the aim of the present work.


## 5.2. The Hu and Sawicki Model

One of the most pressing problems off(R) theories is the need to escape the severe constraints imposed by the Solar System tests. A successful model has been recently proposed by Hu and Sawicki [56] (HS) setting (note that such a model does not pass the matter instability test so that some viable generalizations [99–101] have been proposed):(75)f(R)=R-Rcα(R/Rc)n1+β(R/Rc)n.
As for the double power law model discussed above, there are four parameters which we can be expressed in terms of the cosmographic parameters (q0,j0,s0,l0).As a first step, it is trivial to get:(76)f(R0)=R0-RcαR0cn1+βR0cn,f′(R0)=1-αnRcR0cnR0(1+βR0cn)2,f′′(R0)=αnRcR0cn[(1-n)+β(1+n)R0cn]R02(1+βR0cn)3,f′′′(R0)=αnRcR0cn(An2+Bn+C)R03(1+βR0cn)4.
with R0c=R0/Rc and:(77)A=-β2R0c2n+4βR0cn-1,B=3(1-β2R0c2n),C=-2(1-βR0cn)2.
Equating (76) to the four quantities (ϕ0,ϕ1,ϕ2,ϕ3) defined as above, we could, in principle, solve this system of four equations in four unknowns to get (α,β,Rc,n) in terms of (ϕ0,ϕ1,ϕ2,ϕ3) and then, using (41)–(51) as functions of the cosmographic parameters. However, setting ϕ1=1 as required by (42) gives the only trivial solution αnRc=0 so that the HS model reduces to the Einstein-Hilbert Lagrangian f(R)=R. In order to escape this problem, we can relax the condition f'(R0)=1 to f'(R0)=(1+ɛ)-1. As we have discussed in Section 4, this is the same as assuming that the present day effective gravitational constant Geff,0=GN/f'(R0) only slightly differs from the usual Newtonian one which seems to be a quite reasonable assumption. Under this hypothesis, we can analytically solve for (α,β,Rc,n) in terms of (ϕ0,ɛ,ϕ2,ϕ3). The actual values of (ϕ0,ϕ2,ϕ3) will be no more given by (41)–(44), but we have checked that they deviate from those expressions (note that the correct expressions for (ϕ0,ϕ2,ϕ3) may still formally be written as (41)–(44), but the polynomials entering them are now different and also depend on powers of ɛ) much less than 10% for ɛ up to 10% well below any realistic expectation.With this caveat in mind, we first solve(78)f(R0)=ϕ0,f′′(R0)=(1+ɛ)-1
to get (79)α=n(1+ɛ)ɛ(R0Rc)1-n(1-ϕ0R0)2,β=n(1+ɛ)ɛ(R0Rc)-n[1-ϕ0R0-ɛn(1+ɛ)].
Inserting these expressions in (76), it is easy to check that Rc cancels out so that we can no more determine its value. Such a result is, however, not unexpected. Indeed, (75) can trivially be rewritten as(80)f(R)=R-α̃Rn1+β̃Rn
with α̃=αRc1-n and β̃=βRc-n which are indeed the quantities that are determined by the above expressions for (α,β). Reversing the discussion, the present day values of f(i)(R) depend on (α,β,Rc) only through the two parameters (α̃,β̃). As such, the use of cosmographic parameters is unable to break this degeneracy. However, since Rc only plays the role of a scaling parameter, we can arbitrarily set its value without loss of generality.On the other hand, this degeneracy allows us to get a consistency relation to immediately check whether the HS model is viable or not. Indeed, solving the equationf′′(R0)=ϕ2, we get(81)n=(ϕ0/R0)+[(1+ɛ)/ɛ](1-ϕ2R0)-(1-ɛ)/(1+ɛ)1-ϕ0/R0,
which can then be inserted into the equations f′′′(R0)=ϕ3 to obtain a complicated relation among (ϕ0,ϕ2,ϕ3) which we do not report for sake of shortness. Solving such a relation with respect to ϕ3/ϕ0 and Taylor expanding to first order in ɛ, the constraint we get reads(82)ϕ3ϕ0≃-1+ɛɛϕ2R0[R0(ϕ2ϕ0)+ɛϕ0-11+ɛ(1-2ɛ1-ϕ0/R0)].
If the cosmographic parameters (q0,j0,s0,l0) are known with sufficient accuracy, one could compute the values of (R0,ϕ0,ϕ2,ϕ3) for a given ɛ (eventually using the expressions obtained for ɛ=0) and then check if they satisfied this relation. If this is not the case, one can immediately give off the HS model also without the need of solving the field equations and fitting the data. Actually, given the still large errors on the cosmographic parameters, such a test only remains in the realm of (quite distant) future applications. However, the HS model works for other tests as shown in [56] and so a consistent cosmography analysis has to be combined with them.


## 6. Constraints onf(R)-Derivatives from  the Data

Equations (41)–(51) relate the present day values of f(R) and its first three derivatives to the cosmographic parameters (q0,j0,s0,l0) and the matter density ΩM. In principle, therefore, a measurement of these latter quantities makes it possible to put constraints on f(i)(R0), with i={0,…,3}, and hence on the parameters of a given fourth order theory through the method shown in the previous section. Actually, the cosmographic parameters are affected by errors which obviously propagate onto the f(R) quantities. Actually, the covariance matrix for the cosmographic parameters is not diagonal so that one has also take care of this to estimate the final errors on f(i)(R0). A similar discussion also holds for the errors on the dimensionless ratios η20 and η30 introduced above. As a general rule, indicating with g(ΩM,p) a generic f(R) related quantity depending on ΩM and the set of cosmographic parameters p, its uncertainty reads:(83)σg2=|∂g∂ΩM|2σM2+∑i=1i=4|∂g∂pi|2σpi2+∑i≠j2∂g∂pi∂g∂pjCij,
where Cij are the elements of the covariance matrix (being Cii=σpi2), we have set (p1,p2,p3,p4)=(q0,j0,s0,l0). and assumed that the error σM on ΩM is uncorrelated with those on p. Note that this latter assumption strictly holds if the matter density parameter is estimated from an astrophysical method (such as estimating the total matter in the Universe from the estimated halo mass function). Alternatively, we will assume that ΩM is constrained by the CMBR related experiments. Since these latter mainly probes the very high redshift Universe (z≃zlss≃1089), while the cosmographic parameters are concerned with the present day cosmo, one can argue that the determination of ΩM is not affected by the details of the model adopted for describing the late Universe. Indeed, we can reasonably assume that, whatever is the dark energy candidate or f(R) theory, the CMBR era is well approximated by the standard GR with a model comprising only dust matter. As such, we will make the simplifying (but well motivated) assumption that σM may be reduced to very small values and is uncorrelated with the cosmographic parameters.Under this assumption, the problem of estimating the errors ong(ΩM,p) reduces to estimating the covariance matrix for the cosmographic parameters given the details of the data set used as observational constraints. We address this issue by computing the Fisher information matrix (see, e.g., [102] and references therein) defined as (84)Fij=〈∂2L∂θi∂θj〉
with L=-2lnℒ(θ1,…,θn), ℒ(θ1,…,θn) the likelihood of the experiment, (θ1,…,θn) the set of parameters to be constrained, and 〈…〉 denotes the expectation value. Actually, the expectation value is computed by evaluating the Fisher matrix elements for fiducial values of the model parameters (θ1,…,θn), while the covariance matrix C is finally obtained as the inverse of F.A key ingredient in the computation ofF is the definition of the likelihood which depends, of course, of what experimental constraint one is using. To this aim, it is worth remembering that our analysis is based on fifth order Taylor expansion of the scale factor a(t) so that we can only rely on observational tests probing quantities that are well described by this truncated series. Moreover, since we do not assume any particular model, we can only characterize the background evolution of the Universe, but not its dynamics which, being related to the evolution of perturbations, unavoidably need the specification of a physical model. As a result, the SNeIa Hubble diagram is the ideal test (see the conclusions for further discussion on this issue) to constrain the cosmographic parameters. We therefore defined the likelihood as(85)ℒ(H0,p)∝exp-χ2(H0,p)2,χ2(H0,p)=∑n=1𝒩SNeIa[μobs(zi)-μth(zn,H0,p)σi(zi)]2,
where the distance modulus to redshift z reads:(86)μth(z,H0,p)=25+5log(cH0)+5logdL(z,p),
and dL(z) is the Hubble free luminosity distance:(87)dL(z)=(1+z)∫0zdzH(z)/H0.
Using the fifth order Taylor expansion of the scale factor, we get for dL(z,p) an analytical expression (reported in Appendix A) so that the computation of Fij does not need any numerical  integration (which makes the estimate faster). As a last ingredient, we need to specify the details of the SNeIa survey giving the redshift distribution of the sample and the error on each measurement. Following [103], we adopt (note that, in [103], the authors assume the data are separated in redshift bins so that the error becomes σ2=σsys2/𝒩bin+𝒩bin(z/zmax)2σm2 with 𝒩bin the number of SNeIa in a bin. However, we prefer to not bin the data so that 𝒩bin=1):(88)σ(z)=σsys2+(zzmax)2σm2
with zmax the maximum redshift of the survey, σsys an irreducible scatter in the SNeIa distance modulus and σm to be assigned depending on the photometric accuracy.In order to run the Fisher matrix calculation, we have to set a fiducial model which we set according to theΛCDM predictions for the cosmographic parameters. For ΩM=0.3 and h=0.72 (with h the Hubble constant in units of 100km/s/Mpc), we get(89)(q0,j0,s0,l0)=(-0.55,1.0,-0.35,3.11).
As a first consistency check, we compute the Fisher matrix for a survey mimicking the recent database in [8] thus setting (𝒩SNeIa,σm)=(192,0.33). After marginalizing over h (which, as well known, is fully degenerate with the SNeIa absolute magnitude ℳ), we get for the uncertainties:(90)(σ1,σ2,σ3,σ4)=(0.38,5.4,28.1,74.0),
where we are still using the indexing introduced above for the cosmographic parameters. These values compare reasonably well with those obtained from a cosmographic fitting of the Gold SNeIa dataset (actually, such estimates have been obtained computing the mean and the standard deviation from the marginalized likelihoods of the cosmographic parameters. As such, the central values do not represent exactly the best fit model, while the standard deviations do not give a rigorous description of the error because the marginalized likelihoods are manifestly non Gaussian. Nevertheless, we are mainly interested in an order of magnitude estimate so that we do not care about such statistical details) [104, 105]:(91)q0=-0.90±0.65,j0=2.7±6.7,s0=36.5±52.9,l0=142.7±320.
Because of the Gaussian assumptions it relies on, the Fisher matrix forecasts are known to be lower limits to the accuracy a given experiment can attain on the determination of a set of parameters. This is indeed the case with the comparison suggesting that our predictions are quite optimistic. It is worth stressing, however, that the analysis in [104, 105] used the Gold SNeIa dataset which is poorer in high redshift SNeIa than the [8] one we are mimicking so that larger errors on the higher order parameters (s0,l0) are expected.Rather than computing the errors onf(R0) and its first three derivatives, it is more interesting to look at the precision attainable on the dimensionless ratios (η20,η30 introduced above since they quantify how much deviations from the linear order are present. For the fiducial model we are considering, both η20 and η30 vanish, while, using the covariance matrix for a present day survey and setting σM/ΩM≃10%, their uncertainties read:(92)(σ20,σ30)=(0.04,0.04).
As an application, we can look at Figures 1 and 2 showing how (η20,η30) depend on the present day EoS w0 for f(R) models sharing the same cosmographic parameters of a dark energy model with constant EoS. As it is clear, also considering only the 1σ range, the full region plotted is allowed by such large constraints on (η20,η30) thus meaning that the full class of corresponding f(R) theories is viable. As a consequence, we may conclude that the present day SNeIa data are unable to discriminate between a Λ dominated Universe and this class of fourth order gravity theories.As a next step, we consider an SNAP-like survey [106] thus setting (𝒩SNeIa,σm)=(2000,0.02). We use the same redshift distribution in [103, Table 1] and add 300 nearby SNeIa in the redshift range (0.03,0.08). The Fisher matrix calculation gives for the uncertainties on the cosmographic parameters:(93)(σ1,σ2,σ3,σ4)=(0.08,1.0,4.8,13.7).
The significant improvement of the accuracy in the determination of (q0,j0,s0,l0) translates in a reduction of the errors on (η20,η30) which now read(94)(σ20,σ30)=(0.007,0.008)
having assumed that, when SNAP data will be available, the matter density parameter ΩM has been determined with a precision σM/ΩM~1%. Looking again at Figures 1 and 2, it is clear that the situation is improved. Indeed, the constraints on η20 makes it possible to narrow the range of allowed models with low matter content (the dashed line), while models with typical values of ΩM are still viable for w0 covering almost the full horizontal axis. On the other hand, the constraint on η30 is still too weak so that almost the full region plotted is allowed.Finally, we consider an hypothetical future SNeIa survey working at the same photometric accuracy as SNAP and with the same redshift distribution, but increasing the number of SNeIa up to𝒩SNeIa=6×104 as expected from, for example, DES [107], PanSTARRS [108], SKYMAPPER [109], while still larger numbers may potentially be achieved by ALPACA [110] and LSST [111]. Such a survey can achieve(95)(σ1,σ2,σ3,σ4)=(0.02,0.2,0.9,2.7)
so that, with σM/ΩM~0.1%, we get (96)(σ20,σ30)=(0.0015,0.0016).
Figure 1 shows that, with such a precision on η20, the region of w0 values allowed essentially reduces to the ΛCDM value, while, from Figure 2, it is clear that the constraint on η30 definitively excludes models with low matter content further reducing the range of w0 values to quite small deviations from the w0=-1. We can therefore conclude that such a survey will be able to discriminate between the concordance ΛCDM model and all the f(R) theories giving the same cosmographic parameters as quiessence models other than the ΛCDM itself.A similar discussion may be repeated forf(R) models sharing the same (q0,j0,s0,l0) values as the CPL model even if it is less intuitive to grasp the efficacy of the survey being the parameter space multivalued. For the same reason, we have not explored what is the accuracy on the double power-law or HS models, even if this is technically possible. Actually, one should first estimate the errors on the present day value of f(R) and its three time derivatives and then propagate them on the model parameters using the expressions obtained in Section 6. The multiparameter space to be explored makes this exercise quite cumbersome so that we leave it for a furthcoming work where we will explore in detail how these models compare to the present and future data.


## 7. What We Have Learnt from Cosmography

The recent amount of good quality data have given a new input to the observational cosmology. As often in science, new and better data lead to unexpected discoveries as in the case of the nowadays accepted evidence for cosmic acceleration. However, a fierce and strong debate is still open on what this cosmic speed up implies for theoretical cosmology. The equally impressive amount of different (more or less) viable candidates have also generated a great confusion so that model independent analyses are welcome. A possible solution could come from cosmography rather than assumingad hoc solutions of the cosmological Friedmann equations. Present day and future SNeIa surveys have renewed the interest in the determination of the cosmographic parameters so that it is worth investigating how these quantities can constrain cosmological models.Motivated by this consideration, in the framework of metric formulation off(R) gravity, we have here derived the expressions of the present day values of f(R) and its first three derivatives as function of the matter density parameter ΩM, the Hubble constant H0 and the cosmographic parameters (q0,j0,s0,l0). Although based on a third order Taylor expansion of f(R), we have shown that such relations hold for a quite large class of models so that they are valid tools to look for viable f(R) models without the need of solving the mathematically difficult nonlinear fourth order differential field equations.Notwithstanding the common claim that we live in the era ofprecision cosmology, the constraints on (q0,j0,s0,l0) are still too weak to efficiently apply the program we have outlined above. As such, we have shown how it is possible to establish a link between the popular CPL parameterization of the dark energy equation of state and the derivatives of f(R), imposing that they share the same values of the cosmographic parameters. This analysis has lead to the quite interesting conclusion that the only f(R) function able to give the same values of (q0,j0,s0,l0) as the ΛCDM model is indeed f(R)=R+2Λ. If future observations will tell us that the cosmographic parameters are those of the ΛCDM model, we can therefore rule out all f(R) theories satisfying the hypotheses underlying our derivation of (41)-(44). Actually, such a result should not be considered as a no way out for higher order gravity. Indeed, one could still work out a model with null values of f′′(R0) and f′′′(R0) as required by the above constraints, but nonvanishing higher order derivatives. One could well argue that such a contrived model could be rejected on the basis of the Occam razor, but nothing prevents from still taking it into account if it turns out to be both in agreement with the data and theoretically well founded.If new SNeIa surveys will determine the cosmographic parameters with good accuracy, acceptable constraints on the two dimensionless ratiosη20∝f′′(R0)/f(R0) and η30∝f′′′(R0)/f(R0) could be obtained thus allowing to discriminate among rival f(R) theories. To investigate whether such a program is feasible, we have pursued a Fisher matrix based forecasts of the accuracy future SNeIa surveys can achieve on the cosmographic parameters and hence on (η20,η30). It turns out that a SNAP-like survey can start giving interesting (yet still weak) constraints allowing to reject f(R) models with low matter content, while a definitive improvement is achievable with future SNeIa survey observing ~104 objects thus making it possible to discriminate between ΛCDM and a large class of fourth order theories. It is worth stressing, however, that the measurement of ΩM should come out as the result of a model independent probe such as the gas mass fraction in galaxy clusters which, at present, is still far from the 1% requested precision. On the other hand, one can also rely on the ΩM estimate from the CMBR anisotropy and polarization spectra even if this comes to the price of assuming that the physics at recombination is strictly described by GR so that one has to limit its attention to f(R) models reducing to f(R)∝R during that epoch. However, such an assumption is quite common in many f(R) models available in literature so that it is not a too restrictive limitation.A further remark is in order concerning what kind of data can be used to constrain the cosmographic parameters. The use of the fifth order Taylor expansion of the scale factor makes it possible to not specify any underlying physical model thus relying on the minimalist assumption that the Universe is described by the flat Robertson-Walker metric. While useful from a theoretical perspective, such a generality puts severe limitations to the dataset one can use. Actually, we can only resort to observational tests depending only on the background evolution so that the range of astrophysical probes reduces to standard candles (such as SNeIa and possibly Gamma Ray Bursts [112]) and standard rods (such as the angular size-redshift relation for compact radio sources). Moreover, pushing the Hubble diagram to z~2 may rise the question of the impact of gravitational lensing amplification on the apparent magnitude of the adopted standard candle. The magnification probability distribution function depends on the growth of perturbations [113–117] so that one should worry about the underlying physical model in order to estimate whether this effect biases the estimate of the cosmographic parameters. However, it has been shown [11, 118–121] that the gravitational lensing amplification does not alter significantly the measured distance modulus for z~1 SNeIa. Although such an analysis has been done for GR based models, we can argue that, whatever is the f(R) model, the growth of perturbations finally leads to a distribution of structures along the line of sight that is as similar as possible to the observed one so that the lensing amplification is approximately the same. We can therefore argue that the systematic error made by neglecting lensing magnification is lower than the statistical ones expected by the future SNeIa surveys. On the other hand, one can also try further reducing this possible bias using the method of flux averaging [122] even if, in such a case, our Fisher matrix calculation should be repeated accordingly. It is also worth noting that the constraints on the cosmographic parameters may be tightened by imposing some physically motivated priors in the parameter space. For instance, we can impose that the Hubble parameter H(z) stays always positive over the full range probed by the data or that the transition from past deceleration to present acceleration takes place over the range probed by the data (so that we can detect it). Such priors should be included in the likelihood definition so that the Fisher matrix should be recomputed which is left for a forthcoming work.Although the present day data are still too limited to efficiently discriminate among rivalf(R) models, we are confident that an aggressive strategy aiming at a very precise determination of the cosmographic parameters could offer stringent constraints on higher order gravity without the need of solving the field equations or addressing the complicated problems related to the growth of perturbations. Almost 80 years after the pioneering distance-redshift diagram by Hubble, the old cosmographic approach appears nowadays as a precious observational tool to investigate the new developments of cosmology.


## 8. The Weak-Field Limit off(R)-Gravity

Before facing the problem of galaxy clusters byf(R)-gravity, a discussion is due on the weak-field limit of such a theory which, being of fourth order in metric formalism, could lead to results radically different with respect to the case f(R)=R, the standard second order General Relativity.Let us consider the general action:(97)𝒜=∫d4x-g[f(R)+𝒳ℒm],
where f(R) is an analytic function of the Ricci scalar R, g is the determinant of the metric gμν, 𝒳=16πG/c4 is the coupling constant and ℒm is the standard perfect-fluid matter Lagrangian. Such an action is the straightforward generalization of the Hilbert-Einstein action of GR obtained for f(R)=R. Since we are considering the metric approach, field equations are obtained by varying (97) with respect to the metric(98)f′Rμν-12fgμν-f;μν′+gμν□f'=𝒳2Tμν.
where Tμν=(-2/-g)(δ(-gℒm)/δgμν) is the energy momentum tensor of matter, the prime indicates the derivative with respect to R and □=;σ;σ. We adopt the signature (+,-,-,-).As discussed in details in [123], we deal with the Newtonian and the post-Newtonian limit of f(R)-gravity on a spherically symmetric background. Solutions for the field equations can be obtained by imposing the spherical symmetry [124]:(99)ds2=g00(x0,r)dx02+grr(x0,r)dr2-r2dΩ,
where x0=ct and dΩ is the angular element.To develop the post-Newtonian limit of the theory, one can consider a perturbed metric with respect to a Minkowski backgroundgμν=ημν+hμν. The metric coefficients can be developed as(100)gtt(t,r)≃1+gtt(2)(t,r)+gtt(4)(t,r),grr(t,r)≃-1+grr(2)(t,r),gθθ(t,r)=-r2,gϕϕ(t,r)=-r2sin2θ,
where we put, for the sake of simplicity, c=1,x0=ct→t. We want to obtain the most general result without imposing particular forms for the f(R)-Lagrangian. We only consider analytic Taylor expandable functions(101)f(R)≃f0+f1R+f2R2+f3R3+⋯.
To obtain the post-Newtonian approximation of f(R)-gravity, one has to plug the expansions (100) and (101) into the field equations (98) and then expand the system up to the orders O(0),O(2) and O(4). This approach provides general results and specific (analytic) Lagrangians are selected by the coefficients fi in (101) [123].If we now consider theO(2)-order of approximation, the field equations (98), in the vacuum case, results to be(102)f1rR(2)-2f1gtt,r(2)+8f2R,r(2)-f1rgtt,rr(2)+4f2rR(2)=0,f1rR(2)-2f1grr,r(2)+8f2R,r(2)-f1rgtt,rr(2)=0,2f1grr(2)-r[f1rR(2)-f1gtt,r(2)-f1grr,r(2)+4f2R,r(2)+4f2rR,rr(2)]=0,f1rR(2)+6f2[2R,r(2)+rR,rr(2)]=0,2grr(2)+r[2gtt,r(2)-rR(2)+2grr,r(2)+rgtt,rr(2)]=0.
It is evident that the trace equation (the fourth in the system (102)), provides a differential equation with respect to the Ricci scalar which allows to solve the system at O(2)-order. One obtains the general solution:(103)gtt(2)=δ0-2GMf1r-δ1(t)e-r-ξ3ξr+δ2(t)er-ξ6(-ξ)3/2r,grr(2)=-2GMf1r+δ1(t)[r-ξ+1]e-r-ξ3ξr,-δ2(t)[ξr+-ξ]er-ξ6ξ2r,R(2)=δ1(t)e-r-ξr-δ2(t)-ξer-ξ2ξr,
where ξ≐f1/6f2, f1 and f2 are the expansion coefficients obtained by the f(R)-Taylor series. In the limit f→R, for a point-like source of mass M we recover the standard Schwarzschild solution. Let us notice that the integration constant δ0 is dimensionless, while the two arbitrary time-functions δ1(t) and δ2(t) have, respectively, the dimensions of lenght-1 and lenght-2; ξ has the dimension lenght-2. As extensively discussed in [123], the functions δi(t) (i=1,2) are completely arbitrary since the differential equation system (102) depends only on spatial derivatives. Besides, the integration constant δ0 can be set to zero, as in the standard theory of potential, since it represents an unessential additive quantity. In order to obtain the physical prescription of the asymptotic flatness at infinity, we can discard the Yukawa growing mode in (103) and then the metric is:(104)ds2=[1-2GMf1r-δ1(t)e-r-ξ3ξr]dt2-[1+2GMf1r-δ1(t)(r-ξ+1)e-r-ξ3ξr]dr2-r2dΩ.
The Ricci scalar curvature is(105)R=δ1(t)e-r-ξr.The solution can be given also in terms of gravitational potential. In particular, we have an explicit Newtonian-like term into the definition. The first of (103) provides the second order solution in term of the metric expansion (see the definition (100)). In particular, it is gtt=1+2ϕgrav=1+gtt(2) and then the gravitational potential of an analytic f(R)-theory is(106)ϕgrav=-GMf1r-δ1(t)e-r-ξ6ξr.
Among the possible analytic f(R)-models, let us consider the Taylor expansion where the cosmological term (the above f0) and terms higher than second have been discarded. For the sake of simplicity, we rewrite the Lagrangian (101) as(107)f(R)~a1R+a2R2+⋯
and specify the above gravitational potential (106), generated by a point-like matter distribution, as:(108)ϕ(r)=-3GM4a1r(1+13e-r/L),
where(109)L≡L(a1,a2)=(-6a2a1)1/2.L can be defined as the interaction length of the problem due to the correction to the Newtonian potential. (Such a length is function of the series coefficients, a1 and a2, and it is not a free independent parameter in the following fit procedure.) We have changed the notation to remark that we are doing only a specific choice in the wide class of potentials (106), but the following considerations are completely general.


## 9. Extended Systems

The gravitational potential (108) is a point-like one. Now we have to generalize this solution for extended systems. Let us describe galaxy clusters as spherically symmetric systems and then we have to extend the above considerations to this geometrical configuration. We simply consider the system composed by many infinitesimal mass elements dm each one contributing with a point-like gravitational potential. Then, summing up all terms, namely integrating them on a spherical volume, we obtain a suitable potential. Specifically, we have to solve the integral:(110)Φ(r)=∫0∞r'2dr'∫0πsinθ'dθ'∫02πdω'ϕ(r').
The point-like potential (108) can be split in two terms. The Newtonian component is(111)ϕN(r)=-3GM4a1r.
The extended integral of such a part is the well-known (apart from the numerical constant 3/4a1) expression. It is(112)ΦN(r)=-34a1GM(<r)r,
where M(<r) is the mass enclosed in a sphere with radius r. The correction term(113)ϕC(r)=-GM4a1e-r/Lr
considering some analytical steps in the integration of the angular part  gives the expression(114)ΦC(r)=-2πG4·L∫0∞dr′r′ρ(r′)·e-|r-r′|/L-e-|r+r′|/Lr.
The radial integral is numerically estimated once the mass density is given. We underline a fundamental difference between such a term and the Newtonian one: while in the latter, the matter outside the spherical shell of radius r does not contribute to the potential, in the former external matter takes part to the integration procedure. For this reason we split the corrective potential in two terms:  (i)
ifr'<r:(115)ΦC,int(r)=-2πG4·L∫0rdr'r'ρ(r')·e-|r-r′|/L-e-|r+r′|/Lr=-2πG4·L∫0rdr'r'ρ(r')·e-(r+r′)/L(-1+e2r′/Lr)(ii)
ifr'>r:(116)ΦC,ext(r)=-2πG4·L∫r∞dr'r'ρ(r')·e-|r-r′|/L-e-|r+r′|/Lr=-2πG4·L∫r∞dr′r′ρ(r′)·e-(r+r′)/L(-1+e2r/Lr).The total potential of the spherical mass distribution will be(117)Φ(r)=ΦN(r)+ΦC,int(r)+ΦC,ext(r).
As we will show below, for our purpose, we need the gravitational potential derivative with respect to the variable r; the two derivatives may not be evaluated analytically so we estimate them numerically, once we have given an expression for the total mass density ρ(r). While the Newtonian term gives the simple expression:(118)-dΦNdr(r)=-34a1GM(<r)r2.
The internal and external derivatives of the corrective potential terms are much longer. We do not give them explicitly for sake of brevity, but they are integral functions of the form(119)ℱ(r,r')=∫α(r)β(r)dr'f(r,r')
from which one has: (120)dℱ(r,r')dr=∫α(r)β(r)dr'df(r,r')dr-f(r,α(r))dαdr(r)+f(r,β(r))dβdr(r).
Such an expression is numerically derived once the integration extremes are given. A general consideration is in order at this point. Clearly, the Gauss theorem holds only for the Newtonian part since, for this term, the force law scales as 1/r2. For the total potential (108), it does not hold anymore due to the correction. From a physical point of view, this is not a problem because the full conservation laws are determined, for f(R)-gravity, by the contracted Bianchi identities which assure the self-consistency. For a detailed discussion, see [48, 78, 125].


## 10. The Cluster Mass Profiles

Clusters of galaxies are generally considered self-bound gravitational systems with spherical symmetry and in hydrostatic equilibrium if virialized. The last two hypothesis are still widely used, despite of the fact that it has been widely proved that most clusters show more complex morphologies and/or signs of strong interactions or dynamical activity, especially in their innermost regions [126, 127].Under the hypothesis of spherical symmetry in hydrostatic equilibrium, the structure equation can be derived from the collisionless Boltzmann equation:(121)ddr(ρgas(r)σr2)+2ρgas(r)r(σr2-σθ,ω2)=-ρgas(r)·dΦ(r)dr,
where Φ is the gravitational potential of the cluster, σr and σθ,ω are the mass-weighted velocity dispersions in the radial and tangential directions, respectively, and ρ is gas mass-density. For an isotropic system, it is(122)σr=σθ,ω.
The pressure profile can be related to these quantities by(123)P(r)=σr2ρgas(r).
Substituting (122) and (123) into (121), we have, for an isotropic sphere,(124)dP(r)dr=-ρgas(r)dΦ(r)dr.
For a gas sphere with temperature profile T(r), the velocity dispersion becomes(125)σr2=kT(r)μmp,
where k is the Boltzmann constant, μ≈0.609 is the mean mass particle and mp is the proton mass. Substituting (123) and (125) into (124), we obtain(126)ddr(kT(r)μmpρgas(r))=-ρgas(r)dΦdr
or, equivalently,(127)-dΦdr=kT(r)μmpr[dlnρgas(r)dlnr+dlnT(r)dlnr].
Now the total gravitational potential of the cluster is(128)Φ(r)=ΦN(r)+ΦC(r)
with(129)ΦC(r)=ΦC,int(r)+ΦC,ext(r).
It is worth underlining that if we consider only the standard Newtonian potential, the total cluster mass Mcl,N(r) is composed by gas mass + mass of galaxies + cD-galaxy mass + dark matter and it is given by the expression: (130)Mcl,N(r)=Mgas(r)+Mgal(r)+MCDgal(r)+MDM(r)=-kT(r)μmpGr[dlnρgas(r)dlnr+dlnT(r)dlnr].Mcl,N means the standard estimated Newtonian mass. Generally the galaxy part contribution is considered negligible with respect to the other two components so we have(131)Mcl,N(r)≈Mgas(r)+MDM(r)≈-kT(r)μmpr[dlnρgas(r)dlnr+dlnT(r)dlnr].
Since the gas-mass estimates are provided by X-ray observations, the equilibrium equation can be used to derive the amount of dark matter present in a cluster of galaxies and its spatial distribution.Inserting the previously definedextended-corrected potential of  (128) into  (127), we obtain(132)-dΦNdr-dΦCdr=kT(r)μmpr[dlnρgas(r)dlnr+dlnT(r)dlnr]
from which the extended-corrected mass estimate follows: (133)Mcl,EC(r)+4a13Gr2dΦCdr(r)=4a13[-kT(r)μmpGr(dlnρgas(r)dlnr+dlnT(r)dlnr)].
Since the use of a corrected potential avoids, in principle, the additional requirement of dark matter, the total cluster mass, in this case, is given by(134)Mcl,EC(r)=Mgas(r)+Mgal(r)+MCDgal(r)
and the mass density in the ΦC term is(135)ρcl,EC(r)=ρgas(r)+ρgal(r)+ρCDgal(r)
with the density components derived from observations.In this work, we will use  (133) to compare the baryonic mass profile Mcl,EC(r), estimated from observations, with the theoretical deviation from the Newtonian gravitational potential, given by the expression -(4a1/3G)r2(dΦC/dr)(r). Our goal is to reproduce the observed mass profiles for a sample of galaxy clusters.


## 11. The Galaxy Cluster Sample

The formalism described in Section10 can be applied to a sample of 12 galaxy clusters. We will use the cluster sample studied in [128, 129] which consists of 13 low-redshift clusters spanning a temperature range 0.7÷9.0keV derived from high quality Chandra archival data. In all these clusters, the surface brightness and the gas temperature profiles are measured out to large radii, so that mass estimates can be extended up to r500 or beyond.


### 11.1. The Gas Density Model

The gas density distribution of the clusters in the sample is described by the analytic model proposed in [129]. Such a model modifies the classical β-model to represent the characteristic properties of the observed X-ray surface brightness profiles, that is, the power-law-type cusps of gas density in the cluster center, instead of a flat core and the steepening of the brightness profiles at large radii. Eventually, a second β-model, with a small core radius, is added to improve the model close to the cluster cores. The analytical form for the particle emission is given by(136)npne=n02·(r/rc)-α(1+r2/rc2)3β-α/2·1(1+rγ/rsγ)ϵ/γ+n022(1+r2/rc22)3β2
which can be easily converted to a mass density using the relation(137)ρgas=nT·μmp=1.41.2nemp,
where nT is the total number density of particles in the gas. The resulting model has a large number of parameters, some of which do not have a direct physical interpretation. While this can often be inappropriate and computationally inconvenient, it suits well our case, where the main requirement is a detailed qualitative description of the cluster profiles.In [129],  (136) is applied to a restricted range of distances from the cluster center, that is, between an inner cutoff rmin, chosen to exclude the central temperature bin (≈10 ÷ 20 kpc) where the ICM is likely to be multi-phase, and rdet, where the X-ray surface brightness is at least 3σ significant. We have extrapolated the above function to values outside this restricted range using the following criteria:  (i)
forr<rmin, we have performed a linear extrapolation of the first three terms out to r=0 kpc;(ii)
forr>rdet, we have performed a linear extrapolation of the last three terms out to a distance r̅ for which ρgas(r̅)=ρc, ρc being the critical density of the Universe at the cluster redshift: ρc=ρc,0·(1+z)3. For radii larger than r̅, the gas density is assumed constant at ρgas(r̅).We point out that, in Table1, the radius limit rmin is almost the same as given in the previous definition. When the value given by [129] is less than the cD-galaxy radius, which is defined in the next section, we choose this last one as the lower limit. On the contrary, rmax is quite different from rdet: it is fixed by considering the higher value of temperature profile and not by imaging methods.Table 1
Column 1: Cluster name. Column 2: Richness. Column 3: cluster total mass. Column 4: gas mass. Column 5: galaxy mass. Column 6: cD-galaxy mass. All mass values are estimated atr=rmax. Column 7: ratio of total galaxy mass to gas mass. Column 8: minimum radius. Column 9: maximum radius.
NameRMcl,N (M⊙)Mgas (M⊙)Mgal (M⊙)McDgal(M⊙)galgasrmin (kpc)rmax (kpc)A13304.35874·10142.73866·10135.20269·10121.10568·10120.23861060A26204.45081·10132.76659·10121.71305·10115.16382·10120.2561316A38322.79785·10142.82467·10135.88048·10121.09217·10120.2552751A47828.51832·10141.05583·10142.15567·10131.67513·10120.22591580A90714.87657·10146.38070·10131.34129·10131.66533·10120.245631226A141331.09598·10159.32466·10132.30728·10131.67345·10120.26571506A179525.44761·10145.56245·10134.23211·10121.93957·10120.11791151A199111.24313·10141.00530·10131.24608·10121.08241·10120.2355618A202928.92392·10141.24129·10143.21543·10131.11921·10120.27621771A239012.09710·10152.15726·10144.91580·10131.12141·10120.23831984MKW4—4.69503·10132.83207·10121.71153·10115.29855·10110.2560434RXJ1159—8.97997·10134.33256·10127.34414·10115.38799·10110.2964568We then compute the gas massMgas(r) and the total mass Mcl,N(r), respectively, for all clusters in our sample, substituting  (136) into (137) and (130), respectively; the gas temperature profile has been described in details in § XI B. The resulting mass values, estimated at r=rmax, are listed in Table 1.


### 11.2. The Temperature Profiles

As stressed in Section11.1, for the purpose of this work, we need an accurate qualitative description of the radial behavior of the gas properties. Standard isothermal or polytropic models, or even the more complex one proposed in [129], do not provide a good description of the data at all radii and for all clusters in the present sample. We hence describe the gas temperature profiles using the straightforward X-ray spectral analysis results, without the introduction of any analytic model.X-ray spectral values have been provided by Vikhlinin (private communication). A detailed description of the relative spectral analysis can be found in [128].


### 11.3. The Galaxy Distribution Model

The galaxy density can be modelled as proposed by [86]. Even if the galaxy distribution is a point-distribution instead of a continuous function, assuming that galaxies are in equilibrium with gas, we can use a β-model, ∝r-3, for r<Rc from the cluster center, and a steeper one, ∝r-2.6, for r>Rc, where Rc is the cluster core radius (its value is taken from Vikhlinin 2006). Its final expression is:(138)ρgal(r)={ρgal,1·[1+(rRc)2]-3/2r<Rc,ρgal,2·[1+(rRc)2]-2.6/2r>Rc,
where the constants ρgal,1 and ρgal,2 are chosen in the following way:  (i)
[86] provides the central number density of galaxies in rich compact clusters for galaxies located within a 1.5 h-1 Mpc radius from the cluster center and brighter than m3+2m (where m3 is the magnitude of the third brightest galaxy): ngal,0~103h3 galaxies Mpc-3. Then we fix ρgal,1 in the range ~1034÷ 1036 kg/kpc3. For any cluster obeying the condition chosen for the mass ratio gal-to-gas, we assume a typical elliptical and cD galaxy mass in the range 1012÷1013M⊙.(ii)
the constantρgal,2 has been fixed with the only requirement that the galaxy density function has to be continuous at Rc.We have tested the effect of varying galaxy density in the above range ~1034÷ 1036 kg/kpc3 on the cluster with the lowest mass, namely A262. In this case, we would expect great variations with respect to other clusters; the result is that the contribution due to galaxies and cD-galaxy gives a variation ≤1% to the final estimate of fit parameters. The cD galaxy density has been modelled as described in [130]; they use a Jaffe model of the form(139)ρCDgal=ρ0,J(r/rc)2(1+r/rc)2,
where rc is the core radius while the central density is obtained from MJ=(4/3)πRc3ρ0,J. The mass of the cD galaxy has been fixed at 1.14×1012M⊙, with rc=Re/0.76, with Re=25 kpc being the effective radius of the galaxy. The central galaxy for each cluster in the sample is assumed to have approximately this stellar mass.We have assumed that the total galaxy-component mass (galaxies plus cD galaxy masses) is≈20 ÷ 25% of the gas mass: in [131], the mean fraction of gas versus the total mass (with dark matter) for a cluster is estimated to be 15 ÷ 20%, while the same quantity for galaxies is 3÷5%. This means that the relative mean mass ratio gal-to-gas in a cluster is ≈20÷25%. We have varied the parameters ρgal,1, ρgal,2 and MJ in their previous defined ranges to obtain a mass ratio between total galaxy mass and total gas mass which lies in this range. Resulting galaxy mass values and ratios gal/gas, estimated at r=rmax, are listed in Table 1.In Figure1, we show how each component is spatially distributed. The CD-galaxy is dominant with respect to the other galaxies only in the inner region (below 100 kpc). As already stated in Section 11.1, cluster innermost regions have been excluded from our analysis and so the contribution due to the cD-galaxy is practically negligible in our analysis. The gas is, as a consequence, clearly the dominant visible component, starting from innermost regions out to large radii, being galaxy mass only 20 ÷ 25% of gas mass. A similar behavior is shown by all the clusters considered in our sample.


### 11.4. Uncertainties on Mass Profiles

Uncertainties on the cluster total mass profiles have been estimated performing Monte-Carlo simulations [132]. We proceed to simulate temperature profiles and choose random radius-temperature values couples for each bin which we have in our temperature data given by [128]. Random temperature values have been extracted from a Gaussian distribution centered on the spectral values, and with a dispersion fixed to its 68% confidence level. For the radius, we choose a random value inside each bin. We have performed 2000 simulations for each cluster and perform two cuts on the simulated profile. First, we exclude those profiles that give an unphysical negative estimate of the mass: this is possible when our simulated couples of quantities give rise to too high temperature-gradient. After this cut, we have ≈1500 simulations for any cluster. Then we have ordered the resulting mass values for increasing radius values. Extreme mass estimates (outside the 10 ÷ 90% range) are excluded from the obtained distribution, in order to avoid other high mass gradients which give rise to masses too different from real data. The resulting limits provide the errors on the total mass. Uncertainties on the electron-density profiles has not been included in the simulations, being them negligible with respect to those of the gas-temperature profiles.


### 11.5. Fitting the Mass Profiles

In the above sections, we have shown that, with the aid of X-ray observations, modelling theoretically the galaxy distribution and using  (133), we obtain an estimate of the baryonic content of clusters.We have hence performed a best-fit analysis of the theoretical  (133) (140)Mbar,th(r)=4a13[-kT(r)μmpGr(dlnρgas(r)dlnr+dlnT(r)dlnr)]-4a13Gr2dΦCdr(r)
versus the observed mass contributions(141)Mbar,obs(r)=Mgas(r)+Mgal(r)+MCDgal(r)
Since not all the data involved in the above estimate have measurable errors, we cannot perform an exactχ-square minimization: Actually, we can minimize the quantity:(142)χ2=1N-np-1·∑i=1N(Mbar,obs-Mbar,theo)2Mbar,theo,
where N is the number of data and np=2 the free parameters of the model. We minimize the χ-square using the Markov Chain Monte Carlo Method (MCMC). For each cluster, we have run various chains to set the best parameters of the used algorithm, the Metropolis-Hastings one: starting from an initial parameter vector p (in our case p=(a1,a2)), we generate a new trial point p' from a tested proposal density q(p',p), which represents the conditional probability to get p', given p. This new point is accepted with probability(143)α(p,p′)=min{1,L(d∣p′)P(p′)q(p′,p)L(d∣p)P(p)q(p,p′)},
where d are the data, L(d∣p′)∝exp(-χ2/2) is the likelihood function, and P(p) is the prior on the parameters. In our case, the prior on the fit parameters is related to (109): being L a length, we need to force the ratio a1/a2 to be positive. The proposal density is Gaussian symmetric with respect of the two vectors p and p', namely q(p,p')∝exp(-Δp2/2σ2), with Δp=p-p'; we decide to fix the dispersion σ of any trial distribution of parameters equal to 20% of trial a1 and a2 at any step. This means that the parameter α reduces to the ratio between the likelihood functions.We have run one chain of105 points for every cluster; the convergence of the chains has been tested using the power spectrum analysis from [133]. The key idea of this method is, at the same time, simple and powerful: if we take the power spectra of the MCMC samples, we will have a great correlation on small scales but, when the chain reaches convergence, the spectrum becomes flat (like a white noise spectrum); so that, by checking the spectrum of just one chain (instead of many parallel chains as in Gelmann-Rubin test) will be sufficient to assess the reached convergence. Remanding to [133] for a detailed discussion of all the mathematical steps. Here we calculate the discrete power spectrum of the chains:(144)Pj=|aNj|2
with(145)aNj=1N∑n=0N-1xnexp[i2πjNn],
where N and xn are the length and the element of the sample from the MCMC, respectively, j=1,…,N/2-1. The wavenumber kj of the spectrum is related to the index j by the relation kj=2πj/N. Then we fit it with the analytical template:(146)P(k)=P0(k*/k)α1+(k*/k)α
or in the equivalent logarithmic form(147)lnPj=lnP0+ln[(k*/kj)α1+(k*/kj)α]-γ+rj,
where γ=0.57216 is the Euler-Mascheroni number and rj are random measurement errors with 〈rj〉=0 and 〈rirj〉=δijπ2/6. From the fit, we estimate the two fundamental parameters, P0 and j* (the index corresponding to k*). The first one is the value of the power spectrum extrapolated for k→0 and, from it, we can derive the convergence ratio from r≈(P0/N); if r<0.01, we can assume that the convergence is reached. The second parameter is related to the turning point from a power-law to a flat spectrum. It has to be >20 in order to be sure that the number of points in the sample, coming from the convergence region, are more than the noise points. If these two conditions are verified for all the parameters, then the chain has reached the convergence and the statistics derived from MCMC well describes the underlying probability distribution (typical results are shown in Figures 2 and 3). Following [133] prescriptions, we perform the fit over the range 1≤j≤jmax, with jmax~10j*, where a first estimation of j* can be obtained from a fit with jmax=1000, and then performing a second iteration in order to have a better estimation of it. Even if the convergence is achieved after few thousand steps of the chain, we have decided to run longer chains of 105 points to reduce the noise from the histograms and avoid under- or over estimations of errors on the parameters. The i-σ confidence levels are easily estimated deriving them from the final sample the 15.87-th and 84.13-th quantiles (which define the 68% confidence interval) for i=1, the 2.28-th and 97.72-th quantiles (which define the 95% confidence interval) for i=2 and the 0.13-th and 99.87-th quantiles (which define the 99% confidence interval) for i=3.After the description of the method, let us now comment on the achieved results.


## 11.1. The Gas Density Model

The gas density distribution of the clusters in the sample is described by the analytic model proposed in [129]. Such a model modifies the classical β-model to represent the characteristic properties of the observed X-ray surface brightness profiles, that is, the power-law-type cusps of gas density in the cluster center, instead of a flat core and the steepening of the brightness profiles at large radii. Eventually, a second β-model, with a small core radius, is added to improve the model close to the cluster cores. The analytical form for the particle emission is given by(136)npne=n02·(r/rc)-α(1+r2/rc2)3β-α/2·1(1+rγ/rsγ)ϵ/γ+n022(1+r2/rc22)3β2
which can be easily converted to a mass density using the relation(137)ρgas=nT·μmp=1.41.2nemp,
where nT is the total number density of particles in the gas. The resulting model has a large number of parameters, some of which do not have a direct physical interpretation. While this can often be inappropriate and computationally inconvenient, it suits well our case, where the main requirement is a detailed qualitative description of the cluster profiles.In [129],  (136) is applied to a restricted range of distances from the cluster center, that is, between an inner cutoff rmin, chosen to exclude the central temperature bin (≈10 ÷ 20 kpc) where the ICM is likely to be multi-phase, and rdet, where the X-ray surface brightness is at least 3σ significant. We have extrapolated the above function to values outside this restricted range using the following criteria:  (i)
forr<rmin, we have performed a linear extrapolation of the first three terms out to r=0 kpc;(ii)
forr>rdet, we have performed a linear extrapolation of the last three terms out to a distance r̅ for which ρgas(r̅)=ρc, ρc being the critical density of the Universe at the cluster redshift: ρc=ρc,0·(1+z)3. For radii larger than r̅, the gas density is assumed constant at ρgas(r̅).We point out that, in Table1, the radius limit rmin is almost the same as given in the previous definition. When the value given by [129] is less than the cD-galaxy radius, which is defined in the next section, we choose this last one as the lower limit. On the contrary, rmax is quite different from rdet: it is fixed by considering the higher value of temperature profile and not by imaging methods.Table 1
Column 1: Cluster name. Column 2: Richness. Column 3: cluster total mass. Column 4: gas mass. Column 5: galaxy mass. Column 6: cD-galaxy mass. All mass values are estimated atr=rmax. Column 7: ratio of total galaxy mass to gas mass. Column 8: minimum radius. Column 9: maximum radius.
NameRMcl,N (M⊙)Mgas (M⊙)Mgal (M⊙)McDgal(M⊙)galgasrmin (kpc)rmax (kpc)A13304.35874·10142.73866·10135.20269·10121.10568·10120.23861060A26204.45081·10132.76659·10121.71305·10115.16382·10120.2561316A38322.79785·10142.82467·10135.88048·10121.09217·10120.2552751A47828.51832·10141.05583·10142.15567·10131.67513·10120.22591580A90714.87657·10146.38070·10131.34129·10131.66533·10120.245631226A141331.09598·10159.32466·10132.30728·10131.67345·10120.26571506A179525.44761·10145.56245·10134.23211·10121.93957·10120.11791151A199111.24313·10141.00530·10131.24608·10121.08241·10120.2355618A202928.92392·10141.24129·10143.21543·10131.11921·10120.27621771A239012.09710·10152.15726·10144.91580·10131.12141·10120.23831984MKW4—4.69503·10132.83207·10121.71153·10115.29855·10110.2560434RXJ1159—8.97997·10134.33256·10127.34414·10115.38799·10110.2964568We then compute the gas massMgas(r) and the total mass Mcl,N(r), respectively, for all clusters in our sample, substituting  (136) into (137) and (130), respectively; the gas temperature profile has been described in details in § XI B. The resulting mass values, estimated at r=rmax, are listed in Table 1.


## 11.2. The Temperature Profiles

As stressed in Section11.1, for the purpose of this work, we need an accurate qualitative description of the radial behavior of the gas properties. Standard isothermal or polytropic models, or even the more complex one proposed in [129], do not provide a good description of the data at all radii and for all clusters in the present sample. We hence describe the gas temperature profiles using the straightforward X-ray spectral analysis results, without the introduction of any analytic model.X-ray spectral values have been provided by Vikhlinin (private communication). A detailed description of the relative spectral analysis can be found in [128].


## 11.3. The Galaxy Distribution Model

The galaxy density can be modelled as proposed by [86]. Even if the galaxy distribution is a point-distribution instead of a continuous function, assuming that galaxies are in equilibrium with gas, we can use a β-model, ∝r-3, for r<Rc from the cluster center, and a steeper one, ∝r-2.6, for r>Rc, where Rc is the cluster core radius (its value is taken from Vikhlinin 2006). Its final expression is:(138)ρgal(r)={ρgal,1·[1+(rRc)2]-3/2r<Rc,ρgal,2·[1+(rRc)2]-2.6/2r>Rc,
where the constants ρgal,1 and ρgal,2 are chosen in the following way:  (i)
[86] provides the central number density of galaxies in rich compact clusters for galaxies located within a 1.5 h-1 Mpc radius from the cluster center and brighter than m3+2m (where m3 is the magnitude of the third brightest galaxy): ngal,0~103h3 galaxies Mpc-3. Then we fix ρgal,1 in the range ~1034÷ 1036 kg/kpc3. For any cluster obeying the condition chosen for the mass ratio gal-to-gas, we assume a typical elliptical and cD galaxy mass in the range 1012÷1013M⊙.(ii)
the constantρgal,2 has been fixed with the only requirement that the galaxy density function has to be continuous at Rc.We have tested the effect of varying galaxy density in the above range ~1034÷ 1036 kg/kpc3 on the cluster with the lowest mass, namely A262. In this case, we would expect great variations with respect to other clusters; the result is that the contribution due to galaxies and cD-galaxy gives a variation ≤1% to the final estimate of fit parameters. The cD galaxy density has been modelled as described in [130]; they use a Jaffe model of the form(139)ρCDgal=ρ0,J(r/rc)2(1+r/rc)2,
where rc is the core radius while the central density is obtained from MJ=(4/3)πRc3ρ0,J. The mass of the cD galaxy has been fixed at 1.14×1012M⊙, with rc=Re/0.76, with Re=25 kpc being the effective radius of the galaxy. The central galaxy for each cluster in the sample is assumed to have approximately this stellar mass.We have assumed that the total galaxy-component mass (galaxies plus cD galaxy masses) is≈20 ÷ 25% of the gas mass: in [131], the mean fraction of gas versus the total mass (with dark matter) for a cluster is estimated to be 15 ÷ 20%, while the same quantity for galaxies is 3÷5%. This means that the relative mean mass ratio gal-to-gas in a cluster is ≈20÷25%. We have varied the parameters ρgal,1, ρgal,2 and MJ in their previous defined ranges to obtain a mass ratio between total galaxy mass and total gas mass which lies in this range. Resulting galaxy mass values and ratios gal/gas, estimated at r=rmax, are listed in Table 1.In Figure1, we show how each component is spatially distributed. The CD-galaxy is dominant with respect to the other galaxies only in the inner region (below 100 kpc). As already stated in Section 11.1, cluster innermost regions have been excluded from our analysis and so the contribution due to the cD-galaxy is practically negligible in our analysis. The gas is, as a consequence, clearly the dominant visible component, starting from innermost regions out to large radii, being galaxy mass only 20 ÷ 25% of gas mass. A similar behavior is shown by all the clusters considered in our sample.


## 11.4. Uncertainties on Mass Profiles

Uncertainties on the cluster total mass profiles have been estimated performing Monte-Carlo simulations [132]. We proceed to simulate temperature profiles and choose random radius-temperature values couples for each bin which we have in our temperature data given by [128]. Random temperature values have been extracted from a Gaussian distribution centered on the spectral values, and with a dispersion fixed to its 68% confidence level. For the radius, we choose a random value inside each bin. We have performed 2000 simulations for each cluster and perform two cuts on the simulated profile. First, we exclude those profiles that give an unphysical negative estimate of the mass: this is possible when our simulated couples of quantities give rise to too high temperature-gradient. After this cut, we have ≈1500 simulations for any cluster. Then we have ordered the resulting mass values for increasing radius values. Extreme mass estimates (outside the 10 ÷ 90% range) are excluded from the obtained distribution, in order to avoid other high mass gradients which give rise to masses too different from real data. The resulting limits provide the errors on the total mass. Uncertainties on the electron-density profiles has not been included in the simulations, being them negligible with respect to those of the gas-temperature profiles.


## 11.5. Fitting the Mass Profiles

In the above sections, we have shown that, with the aid of X-ray observations, modelling theoretically the galaxy distribution and using  (133), we obtain an estimate of the baryonic content of clusters.We have hence performed a best-fit analysis of the theoretical  (133) (140)Mbar,th(r)=4a13[-kT(r)μmpGr(dlnρgas(r)dlnr+dlnT(r)dlnr)]-4a13Gr2dΦCdr(r)
versus the observed mass contributions(141)Mbar,obs(r)=Mgas(r)+Mgal(r)+MCDgal(r)
Since not all the data involved in the above estimate have measurable errors, we cannot perform an exactχ-square minimization: Actually, we can minimize the quantity:(142)χ2=1N-np-1·∑i=1N(Mbar,obs-Mbar,theo)2Mbar,theo,
where N is the number of data and np=2 the free parameters of the model. We minimize the χ-square using the Markov Chain Monte Carlo Method (MCMC). For each cluster, we have run various chains to set the best parameters of the used algorithm, the Metropolis-Hastings one: starting from an initial parameter vector p (in our case p=(a1,a2)), we generate a new trial point p' from a tested proposal density q(p',p), which represents the conditional probability to get p', given p. This new point is accepted with probability(143)α(p,p′)=min{1,L(d∣p′)P(p′)q(p′,p)L(d∣p)P(p)q(p,p′)},
where d are the data, L(d∣p′)∝exp(-χ2/2) is the likelihood function, and P(p) is the prior on the parameters. In our case, the prior on the fit parameters is related to (109): being L a length, we need to force the ratio a1/a2 to be positive. The proposal density is Gaussian symmetric with respect of the two vectors p and p', namely q(p,p')∝exp(-Δp2/2σ2), with Δp=p-p'; we decide to fix the dispersion σ of any trial distribution of parameters equal to 20% of trial a1 and a2 at any step. This means that the parameter α reduces to the ratio between the likelihood functions.We have run one chain of105 points for every cluster; the convergence of the chains has been tested using the power spectrum analysis from [133]. The key idea of this method is, at the same time, simple and powerful: if we take the power spectra of the MCMC samples, we will have a great correlation on small scales but, when the chain reaches convergence, the spectrum becomes flat (like a white noise spectrum); so that, by checking the spectrum of just one chain (instead of many parallel chains as in Gelmann-Rubin test) will be sufficient to assess the reached convergence. Remanding to [133] for a detailed discussion of all the mathematical steps. Here we calculate the discrete power spectrum of the chains:(144)Pj=|aNj|2
with(145)aNj=1N∑n=0N-1xnexp[i2πjNn],
where N and xn are the length and the element of the sample from the MCMC, respectively, j=1,…,N/2-1. The wavenumber kj of the spectrum is related to the index j by the relation kj=2πj/N. Then we fit it with the analytical template:(146)P(k)=P0(k*/k)α1+(k*/k)α
or in the equivalent logarithmic form(147)lnPj=lnP0+ln[(k*/kj)α1+(k*/kj)α]-γ+rj,
where γ=0.57216 is the Euler-Mascheroni number and rj are random measurement errors with 〈rj〉=0 and 〈rirj〉=δijπ2/6. From the fit, we estimate the two fundamental parameters, P0 and j* (the index corresponding to k*). The first one is the value of the power spectrum extrapolated for k→0 and, from it, we can derive the convergence ratio from r≈(P0/N); if r<0.01, we can assume that the convergence is reached. The second parameter is related to the turning point from a power-law to a flat spectrum. It has to be >20 in order to be sure that the number of points in the sample, coming from the convergence region, are more than the noise points. If these two conditions are verified for all the parameters, then the chain has reached the convergence and the statistics derived from MCMC well describes the underlying probability distribution (typical results are shown in Figures 2 and 3). Following [133] prescriptions, we perform the fit over the range 1≤j≤jmax, with jmax~10j*, where a first estimation of j* can be obtained from a fit with jmax=1000, and then performing a second iteration in order to have a better estimation of it. Even if the convergence is achieved after few thousand steps of the chain, we have decided to run longer chains of 105 points to reduce the noise from the histograms and avoid under- or over estimations of errors on the parameters. The i-σ confidence levels are easily estimated deriving them from the final sample the 15.87-th and 84.13-th quantiles (which define the 68% confidence interval) for i=1, the 2.28-th and 97.72-th quantiles (which define the 95% confidence interval) for i=2 and the 0.13-th and 99.87-th quantiles (which define the 99% confidence interval) for i=3.After the description of the method, let us now comment on the achieved results.


## 12. Results

The numerical results of our fitting analysis are summarized in Table2; we give the best fit values of the independent fitting parameters a1 and a2, and of the gravitational length L, considered as a function of the previous two quantities. In Figures 3–5, we give the typical results of fitting, with histograms and power spectrum of samples derived by the MCMC, to assess the reached convergence (flat spectrum at large scales).Table 2
Column 1: Cluster name. Column 2: first derivative coefficient,a1, of f(R) series. Column 3: 1σ confidence interval for a1. Column 4: second derivative coefficient, a2, of f(R) series. Column 5: 1σ confidence interval for a2. Column 6: characteristic length, L, of the modified gravitational potential, derived from a1 and a2. Column 7: 1σ confidence interval for L.
Namea1[a1-1σ, a1+1σ]a2[a2-1σ, a2+1σ]L[L-1σ, L+1σ](kpc2)(kpc2)(kpc)(kpc)A1330.085[0.078, 0.091]-4.98·103[-2.38·104, -1.38·103]591.78[323.34, 1259.50]A2620.065[0.061, 0.071]-10.63[-57.65, -3.17]31.40[17.28, 71.10]A3830.099[0.093, 0.108]-9.01·102[-4.10·103, -3.14·102]234.13[142.10, 478.06]A4780.117[0.114, 0.122]-4.61·103[-1.01·104, -2.51·103]484.83[363.29, 707.73]A9070.129[0.125, 0.136]-5.77·103[-1.54·104, -2.83·103]517.30[368.84, 825.00]A14130.115[0.110, 0.119]-9.45·104[-4.26·105, -3.46·104]2224.57[1365.40, 4681.21]A17950.093[0.084, 0.103]-1.54·103[-1.01·104, -2.49·102]315.44[133.31, 769.17]A19910.074[0.072, 0.081]-50.69[-3.42·102, -13]64.00[32.63, 159.40]A20290.129[0.123, 0.134]-2.10·104[-7.95·104, -8.44·103]988.85[637.71, 1890.07]A23900.149[0.146, 0.152]-1.40·106[-5.71·106, -4.46·105]7490.80[4245.74, 15715.60]MKW40.054[0.049, 0.060]-23.63[-1.15·102, -8.13]51.31[30.44, 110.68]RXJ11590.048[0.047, 0.052]-18.33[-1.35·102, -4.18]47.72[22.86, 125.96](a) histogram of the sample points for parametera1 in Abell 383 coming out the MCMC implementation used to estimate best fit values and errors for our fitting procedure as described in  Section 11.5. Binning (horizontal axis) and relative frequencies (vertical axis) are given by automatic procedure from Mathematica6.0. (b) power spectrum test on sample chain for parameter a1 using the method described in Section 11.5. Black line is the logarithm of the analytical template  (146) for power spectrum; gray line is the discrete power spectrum obtained using  (144)–(145).
(a)(b)The main property of our results is  the presence of atypical scale for each cluster above which our model works really good (typical relative differences are less than 5%), while for lower scale there is a great difference. It is possible to see, by a rapid inspection, that this turning-point is located at a radius ≈150 kpc. Except for very large clusters, it is clear that this value is independent of the cluster, being approximately the same for any member of the considered sample.There are two main independent explanations that could justify this trend: limits due to a break in the state of hydrostatic equilibrium or limits in the series expansion of thef(R)-models.If the hypothesis of hydrostatic equilibrium is not correct, then we are in a regime where the fundamental relations (121)–(127), are not working. As discussed in [128], the central (70 kpc) region of every cluster is strongly affected by radiative cooling and thus it cannot directly be related to the depth of the cluster potential well. This means that, in this region, the gas is not in hydrostatic equilibrium but in a multi-phase, turbulent state, mainly driven by some astrophysical, nongravitational interaction. In this case, the gas cannot be used as a good standard tracer.We have also to consider another limit of our modelling: the requirement that thef(R)-function is Taylor expandable. The corrected gravitational potential which we have considered is derived in the weak field limit, which means(148)R-R0≪a1a2,
where R0 is the background value of the curvature. If this condition is not satisfied, the approach does not work (see [123] for a detailed discussion of this point). Considering that a1/a2 has the dimension of length-2 this condition defines the length scale where our series approximation can work. In other words, this indicates the limit in which the model can be compared with data.For the considered sample, the fit of the parametersa1 and a2, spans the length range {19;200} kpc (except for the biggest cluster). It is evident that every galaxy cluster has a proper gravitational length scale. It is worth noticing that a similar situation, but at completely different scales, has been found out for low surface brightness galaxies modelled by f(R)-gravity [78].Considering the data at our disposal and the analysis which we have performed, it is not possible to quantify exactly the quantitative amount of these two different phenomena (i.e., the radiative cooling and the validity of the weak field limit). However, they are not mutually exclusive but should be considered in details in view of a more refined modelling. (Other secondary phenomena as cooling flows, merger and asymmetric shapes have to be considered in view of a detailed modelling of clusters. However, in this work, we are only interested to show that extended gravity could be a valid alternative to dark matter in order to explain the cluster dynamics.)Similar issues are present also in [134]: they use the the Metric-Skew-Tensor-Gravity (MSTG) as a generalization of the Einstein General Relativity and derive the gas mass profile of a sample of clusters with gas being the only baryonic component of the clusters. They consider some clusters included in our sample (in particular, A133, A262, A478, A1413, A1795, A2029, MKW4) and they find the same different trend for r≤200 kpc, even if with a different behavior with respect of us: our model gives lower values than X-ray gas mass data while their model gives higher values with respect to X-ray gas mass data. This stresses the need for a more accurate modelling of the gravitational potential.However, our goal is to show that potential (108) is suitable to fit the mass profile of galaxy clusters and that it comes from a self-consistent theory.In general, it can be shown that the weak field limit of extended theories of gravity has Yukawa-like corrections [77, 135]. Specifically, given theory of gravity of order (2n+2), the Yukawa corrections to the Newtonian potential are n [136]. This means that if the effective Lagrangian of the theory is(149)ℒ=f(R,□R,…,□kR,…,□nR)-g
we have(150)ϕ(r)=-GMr[1+∑k=1nαke-r/Lk].
Standard General Relativity, where Yukawa corrections are not present, is recovered for n=0 (second order theory) while the f(R)-gravity is obtained for n=1 (fourth-order theory). Any □ operator introduces two further derivation orders in the field equations. This kind of Lagrangian comes out when quantum field theory is formulated on curved spacetime [137]. In the series (150), G is the value of the gravitational constant considered at infinity, Lk is the interaction length of the k-th component of the nonNewtonian corrections. The amplitude αk of each component is normalized to the standard Newtonian term; the sign of αk tells us if the corrections are attractive or repulsive (see [138] for details). Moreover, the variation of the gravitational coupling is involved. In our case, we are taking into account only the first term of the series. It is the the leading term. Let us rewrite (108) as(151)ϕ(r)=-GMr[1+α1e-r/L1].
The effect of nonNewtonian term can be parameterized by {α1,L1} which could be a useful parameterisation which respect to our previous {a1,a2} or {Geff,L} with Geff=3G/(4a1). For large distances, where r≫L1, the exponential term vanishes and the gravitational coupling is G. If r≪L1, the exponential becomes 1 and, by differentiating (151) and comparing with the gravitational force measured in laboratory, we get(152)Glab=G[1+α1(1+rL1)e-r/L1]≃G(1+α1),
where Glab=6.67×10-8 g-1cm3s-2 is the usual Newton constant measured by Cavendish-like experiments. Of course, G and Glab coincide in the standard Newtonian gravity. It is worth noticing that, asymptotically, the inverse square law holds but the measured coupling constant differs by a factor (1+α1). In general, any correction introduces a characteristic length that acts at a certain scale for the self-gravitating systems as in the case of galaxy cluster which we are examining here. The range of Lk of the kth-component of nonNewtonian force can be identified with the mass mk of a pseudo-particle whose effective Compton's length can be defined as(153)Lk=ℏmkc.
The interpretation of this fact is that, in the weak energy limit, fundamental theories which attempt to unify gravity with the other forces introduce, in addition to the massless graviton, particles with mass which also carry the gravitational interaction [139]. See, in particular, [140] for f(R)-gravity. These masses are related to effective length scales which can be parameterized as(154)Lk=2×10-5(1eVmk)cm.
There have been several attempts to experimentally constrain Lk and αk (and then mk) by experiments on scales in the range 1cm<r<1000km, using different techniques [141–143]. In this case, the expected masses of particles which should carry the additional gravitational force are in the range 10-13eV<mk<10-5eV. The general outcome of these experiments, even retaining only the term k=1, is that geophysical window between the laboratory and the astronomical scales has to be taken into account. In fact, the range(155)|α1|~10-2,L1~102÷103m
is not excluded at all in this window. An interesting suggestion has been given by Fujii [144], which proposed that the exponential deviation from the Newtonian standard potential could arise from the microscopic interaction which couples the nuclear isospin and the baryon number.The astrophysical counterparts of these nonNewtonian corrections seemed ruled out till some years ago due to the fact that experimental tests of General Relativity seemed to predict the Newtonian potential in the weak energy limit, “inside” the Solar System. However, as it has been shown, several alternative theories seem to evade the Solar System constraints (see [140, 145] and the references therein for recent results) and, furthermore, indications of an anomalous, long-range acceleration revealed from the data analysis of Pioneer 10/11, Galileo, and Ulysses spacecrafts (which are now almost outside the Solar System) makes these Yukawa-like corrections come again into play [146]. Besides, it is possible to reproduce phenomenologically the flat rotation curves of spiral galaxies considering the values(156)α1=-0.92,L1~40kpc.
The main hypothesis of this approach is that the additional gravitational interaction is carried by some ultra-soft boson whose range of mass is m1~10-27÷10-28 eV. The action of this boson becomes efficient at galactic scales without the request of enormous amounts of dark matter to stabilize the systems [147].Furthermore, it is possible to use a combination of two exponential correction terms and give a detailed explanation of the kinematics of galaxies and galaxy clusters, again without dark matter model [143].It is worthwhile to note that both the spacecrafts measurements and galactic rotation curves indications come from “outside” the usual Solar System boundaries used up to now to test General Relativity. However, the above resultsdo not come from any fundamental theory to explain the outcome of Yukawa corrections. In their contexts, these terms are phenomenological.Another important remark in this direction deserves the fact that some authors [148] interpret also the experiments on cosmic microwave background like the experiment BOOMERANG and WMAP [20, 23] in the framework of modified Newtonian dynamics again without invoking any dark matter model.All these facts point towards the line of thinking that also corrections to the standard gravity have to be seriously taken into account beside dark matter searches.In our case, the parametersa1,2, which determine the gravitational correction and the gravitational coupling, come out “directly” from a field theory with the only requirement that the effective action of gravity could be more general than the Hilbert-Einstein theory f(R)=R. This main hypothesis comes from fundamental physics motivations due to the fact that any unification scheme or quantum field theory on curved space have to take into account higher order terms in curvature invariants [137]. Besides, several recent results point out that such corrections have a main role also at astrophysical and cosmological scales. For a detailed discussion, see [48, 125, 149].With this philosophy in mind, we have plotted the trend ofa1 as a function of the density in Figure 6. As one can see, its values are strongly constrained in a narrow region of the parameter space, so that a1 can be considered a “tracer” for the size of gravitational structures. The value of a1 range between {0.8÷0.12} for larger clusters and {0.4÷0.6} for poorer structures (i.e., galaxy groups like MKW4 and RXJ1159). We expect a particular trend when applying the model to different gravitational structures. In Figure 6, we give characteristic values of density which range from the biggest structure, the observed Universe (bold dashed vertical line), to the smallest one, the Sun (vertical dotted line), through intermediate steps like clusters (vertical dashed line) and galaxies (vertical dot-dashed line). The bold black horizontal line represents the Newtonian limit a1=3/4 and the boxes indicate the possible values of a1 that we obtain by applying our theoretical model to different structures.Figure 6
Density versusa1: predictions on the behavior of a1. The horizontal black bold line indicates the Newtonian-limit, a1→3/4 which we expect to be realized on scales comparable with Solar System. Vertical lines indicate typical approximated values of matter density (without dark matter) for different gravitational structures: Universe (bold dashed) with critical density ρcrit≈10-26kg/m3; galaxy clusters (dashed line) with ρcl≈10-23kg/m3; galaxies (dot-dashed) with ρgal≈10-11kg/m3; sun (dotted) with ρsun≈103kg/m3. Arrows and boxes show the predicted trend for a1.Similar considerations hold also for the characteristic gravitational lengthL directly related to both a1 and a2. The parameter a2 shows a very large range of variation {-106÷-10} with respect to the density (and the mass) of the clusters. The value of L changes with the sizes of gravitational structure (see Figure 7), so it can be considered, beside the Schwarzschild radius, a sort of additional gravitational radius. Particular care must be taken when considering Abell 2390, which shows large cavities in the X-ray surface brightness distribution, and whose central region, highly asymmetric, is not expected to be in hydrostatic equilibrium. All results at small and medium radii for this cluster could hence be strongly biased by these effects [129]; the same will hold for the resulting exceptionally high value of L. Figure 7 shows how observational properties of the cluster, which well characterize its gravitational potential (such as the average temperature and the total cluster mass within r500, plotted in the left and right panel, resp.), well correlate with the characteristic gravitational length L.Single temperature fit to the total cluster spectrum (a) and total cluster mass withinr500 (given as a function of M⊙) (b) are plotted as a function of the characteristic gravitational length L. Temperature and mass values are from [129].
(a)(b)For clusters, we can define a gas-density-weighted and a gas-mass-weighted mean, both depending on the series parametersa1,2. We have(157)〈L〉ρ=318kpc〈a2〉ρ=-3.40·104,〈L〉M=2738kpc〈a2〉M=-4.15·105.
It is straightforward to note the correlation with the sizes of the cluster cD-dominated-central region and the "gravitational" interaction length of the whole cluster. In other words, the parameters a1,2, directly related to the first and second derivative of a given analytic f(R)-model determine the characteristic sizes of the self gravitating structures. We have plotted the baryonic mass vs radii for the clusters of the sample in Figures 8, 9, 10, 11, and 12.Figure 8
As an example of the above results, we have plotted the baryonic mass vs radii for Abell A133. Dashed line is the experimental-observed estimation  (141) of baryonic matter component (i.e., gas, galaxies and cD-galaxy); solid line is the theoretical estimation  (140) for baryonic matter component. Dotted lines are the 1-σ confidence levels given by errors on fitting parameters plus statistical errors on mass profiles as discussed in Section 11.4 in the right panel.Figure 9
As the above case, for cluster Abell 383.Figure 10
As the above cases, for cluster Abell 478.Figure 11
As the above cases, for cluster Abell 1413.Figure 12
As the above cases, for cluster Abell 2029.


## 13. What We Have Learnt from Clusters

We have investigated the possibility that the high observational mass-to-light ratio of galaxy clusters could be addressed byf(R)-gravity without assuming huge amounts of dark matter. We point out that this proposal comes out from the fact that, up to now, no definitive candidate for dark-matter has been observed at fundamental level and then alternative solutions to the problem should be viable. Furthermore, several results in f(R)-gravity seem to confirm that valid alternatives to ΛCDM can be achieved in cosmology. Besides, as discussed in the introduction, the rotation curves of spiral galaxies can be explained in the weak field limit of f(R)-gravity. Results of our analysis go in this direction.We have chosen a sample of relaxed galaxy clusters for which accurate spectroscopic temperature measurements and gas mass profiles are available. For the sake of simplicity, and considered the  sample at our disposal, every cluster has been modelled as a self-bound gravitational system with spherical symmetry and in hydrostatic equilibrium. The mass distribution has been described by a corrected gravitational potential obtained from a generic analyticf(R)-theory. In fact, as soon as f(R)≠R, Yukawa-like exponential corrections emerge in the weak field limit while the standard Newtonian potential is recovered only for f(R)=R, the Hilbert-Einstein theory.Our goal has been to analyze if the dark-matter content of clusters can be addressed by these correction potential terms. As discussed in detail in the previous sections and how it is possible to see by a rapid inspection of figures, the clusters of the sample are consistent with the proposed model at1σ confidence level. This shows, at least qualitatively, that the high mass-to-light ratio of clusters can be explained by using a modified gravitational potential. The good agreement is achieved on distance scales starting from 150 kpc up to 1000 kpc. The differences observed at smaller scales can be ascribed to nongravitational phenomena, such as cooling flows, or to the fact that the gas mass is not a good tracer at this scales. The remarkable result is that we have obtained a consistent agreement with data only using the corrected gravitational potential in a large range of radii. In order to put in evidence this trend, we have plotted the baryonic mass versus radii considering, for each cluster, the scale where the trend is clearly evident.In our knowledge, the fact thatf(R)-gravity could work at these scales has been only supposed but never achieved by a direct fitting with data (see [82, 85] for a review). Starting from the series coefficients a1 and a2, it is possible to state that, at cluster scales, two characteristic sizes emerge from the weak field limit of the theory. However, at smaller scales, for example, Solar System scales, standard Newtonian gravity has to be dominant in agreement with observations and experiments.In summary, if our considerations are right, gravitational interaction depends on the scale and theinfrared limit is led by the series coefficient of the considered effective gravitational Lagrangian. Roughly speaking, we expect that starting from cluster scale to galaxy scale, and then down to smaller scales as Solar System or Earth, the terms of the series lead the clustering of self-gravitating systems beside other nongravitational phenomena. In our case, the Newtonian limit is recovered for a1→3/4 and L(a1,a2)≫r at small scales and for L(a1,a2)≪r at large scales. In the first case, the gravitational coupling has to be redefined, in the second G∞≃G. In these limits, the linear Ricci term is dominant in the gravitational Lagrangian and the Newtonian gravity is restored [136]. Reversing the argument, this could be the starting point to achieve a theory capable of explaining the strong segregation in masses and sizes of gravitationally-bound systems.


## 14. Conclusions

The present status of art of cosmology shows that the Standard Cosmological Model, based on General Relativity, nucleosynthesis, cosmic abundances and large scale structure, has some evident difficulties. These ones, first of all, rely on some lack of a self-consistent formulation of missing matter and cosmic acceleration issues; such shortcomings give rise to further difficulties in interpreting observational data. With an aphorism, one can say thatwe have a book, but not the alphabet to read it.Nowadays there two main philosophical approaches aimed to solve this problem. From one side, there are researchers which try to solve shortcomings of Standard Cosmological Model assuming that General Relativity is right but we need some exotic, invisible kinds of energy and matter to explain cosmic dynamics and large scale structure. On the other side, there are people which believe that General Relativity is not the definitive and comprehensive theory of gravity, and that it should be revised at ultraviolet scales (quantum gravity) and infrared scales (extragalactic and cosmic scales). In the latter case, dark energy and dark matter could be nothing else but the signals that we need a more general theory at large scales, also if General Relativity works very well up to Solar System scales. To some extent, this could be seen as a sort of philosophical debate without solution, but there are possibilities to move the question toward a physical viewpoint.Thef(R)-gravity is strictly related to the second point of view. It is a fruitful approach to generalize General Relativity towards the solution also if, most of the models in literature are nothing else but phenomenological models. It is interesting to note that as soon as Einstein formulated his General Relativity, many authors (and Einstein himself) started to explore other possibilities (see [48] for a review). At the beginning, these researches were mainly devoted to check the mathematical consistency of General Relativity but the issues to achieve the unification of gravity with the other interactions (e.g., electromagnetism) pushed several authors to develop alternative gravity theories. Today, one of the goals of alternative gravity is to understand the effective content and dynamics of the Universe. This question is recently become dramatic since assuming that more than 95% of cosmic matter-energy is unknown at fundamental level is highly disturbing. Alternative gravity could be a way out to this situation. The present status of observations, also if we are living in the era of Precision Cosmology, does not allow in discriminating between alternative gravity, from one side, and the presence of dark energy and dark matter, from the other side (the forthcoming LHC experiments should aid in this sense if new fundamental particles will be detected).However, as discussed in this review, cosmography may be a useful tool to discriminate among different cosmological models being, by definition, a model-independent approach: any cosmographic parameter can be estimated without assigning ana priori cosmological model. So cosmography can be used in two ways  (i)
One can use it to discriminate between General Relativity and alternative theories. This issue strictly depends on the possibility to have good quality data at disposal. We need some minimum sensibility and error requirements on data surveys to solve this question. At the moment, we have not them and we are not able to do this since standard candles are not available at very high red shifts [112].(ii)
We can use the cosmographic parameters to constraint cosmological models as we have done in this paper forf(R)-gravity. Being these parameters model-independent, they results natural “priors” to any theory. As above, the accuracy in estimating them is a crucial issue.We have used, essentially, SNeIa but other classes of objects have to be considered in order to improve such an accuracy (e.g., CMBR, bright galaxies, GRBs, BAOs, weak lensing and so on). Forthcoming space missions will be extremely useful in this sense.Beside cosmography, we have discussed also iff(R)-gravity could be useful to address the problem of mass profile and dynamics of galaxy clusters. This issue is crucial in view of achieving any correct model for large scale structure.Taking into account the weak-field limit of a generic analyticf(R)-function, it is possible to obtain a scale-dependent gravity, where scales of self-gravitating systems could naturally emerge. In this way, one could successfully explain dark matter profiles ranging from galaxies to clusters of galaxies. The results are preliminary but seems to indicate a way in which the dark matter puzzle could be completely solved.In conclusion, the main lesson of this work is that since it is very difficult to discriminate among the huge amount of cosmological models which try to explain the data (deductive approach), it could be greatly fruitful to “reconstruct” the final cosmological model by an inductive approach, that is without imposing it a priori but adopting the philosophy to use the minimum number of parameters. (Following the Occam razor prescriptions:“Entia non sunt multiplicanda praeter necessitatem.”) This “inverse scattering approach” could be not fully satisfactory but could lead to self-consistent results.



---

*Source: 217420-2010-03-25.xml*