Constraining $f(R)$ gravity by the Large Scale Structure

Over the past decades, General Relativity and the concordance $\Lambda$CDM model have been successfully tested using several different astrophysical and cosmological probes based on large datasets ({\it precision cosmology}). Despite their successes, some shortcomings emerge due to the fact that General Relativity should be revised at infrared and ultraviolet limits and to the fact that the fundamental nature of Dark Matter and Dark Energy is still a puzzle to be solved. In this perspective, $f(R)$ gravity have been extensively investigated being the most straightforward way to modify General Relativity and to overcame some of the above shortcomings. In this paper, we review various aspects of $f(R)$ gravity at extragalactic and cosmological levels. In particular, we consider cluster of galaxies, cosmological perturbations, and N-Body simulations, focusing on those models that satisfy both cosmological and local gravity constraints. The perspective is that some classes of $f(R)$ models can be consistently constrained by Large Scale Structure.


Introduction
The concordance ΛCDM cosmological model is based on Einstein's general relativity (GR), the standard model of particles with the inclusion of two new ingredients, which are the cosmological constant Λ and dark matter. It provides a theoretical framework to describe the formation and the evolution of the structures in the Universe. It has been successfully tested using many different datasets, such as the Supernovae Type Ia (SNeIa) [1][2][3][4], the matter power spectrum from the two-degree field (2dF) survey of galaxies [5,6] and from the Sloan Digital Sky Survey (SDSS) [7], the temperature fluctuations of the Cosmic Microwave Background (CMB) radiation [8][9][10][11][12][13][14][15][16][17][18], the baryonic acoustic oscillations (BAO) [19], and so on. Despite its successes, the following questions are still need to be answered: Is GR sufficient to explain all of the gravitational phenomena from collapsed objects to the evolution of the Universe and the formation of the structures? Does the theory need to be changed, or at least modified? There exist two different approaches to solve the debate: preserve the successes of GR by incorporating new particles and/or scalar fields not yet observed and improve the accuracy of the data to further verify the model; or modify the theory of gravity to make it compatible with quantum mechanics and cosmological observations without introducing additional particles and fields. Independently of the preferred approach, having alternatives demands testing the model more deeply. In the concordance ΛCDM model, the most important energy density component is the cosmological constant Λ (or more in general, the dark energy (DE)), required to explain the accelerated expansion of the Universe, with Ω Λ = 0.691 ± 0.006 [13] in units of the critical density. The second in importance is dark matter (DM), with Ω DM = 0.259 ± 0.005 [13], needed to explain the emergence of the large-scale structure (LSS). However, the lack of a full comprehension of the fundamental nature of those two components, whether particles or scalar fields, is completely unknown and has demanded to further test the GR to understand if it is the effective theory of gravity. Additionally, the need for requiring two unknown components to fully explain the astrophysical and cosmological observations has been interpreted as a breakdown of the theory at those scales [20].
As an alternative to introducing extra matter components in the stress-energy tensor, it is possible to change the geometrical description of the gravitational interaction, adding in the Hilbert-Einstein Lagrangian, higher-order curvature invariants, such as R 2 , R µν R µν , R µναβ R µναβ , R R or R k R), and minimally or non-minimally coupled terms between scalar fields and geometry (such as φ 2 R) [20][21][22][23][24][25][26][27][28]. Extended theories of gravity (ETGs) can be classified as: scalar-tensor theories, when the geometry is non-minimally coupled to some scalar field; and higher order theories, when the action contains derivatives of the metric components of order higher than the second. ETGs differ from GR in many concrete aspects and, therefore, are testable with current and forthcoming experiments. The most straightforward way to extend the theory of gravitation is replacing the Einstein-Hilbert Lagrangian, which is linear in the Ricci scalar, R, with a more general function of the curvature f (R). Since the exact functional form of the f (R)-Lagrangian is unknown, theoretical predictions should be matched to all available data (at all scales) to accommodate the observations. f (R) models have been used to test: stellar formation and evolution [29][30][31][32][33][34]; emission of gravitational waves and constraints on the massive graviton modes [35][36][37][38][39][40][41][42][43][44]; the flat rotation curves of spiral galaxies [45]; the velocity dispersion of elliptical galaxies [46]; to describe a cluster of galaxies [47][48][49][50][51][52][53]; the structure formation and the evolution of the Universe [14,[54][55][56][57][58][59][60][61][62][63][64][65]. Although all tests provide clear indications that the f (R)-gravity models could overcame the shortcomings of the ΛCDM model, a final and self-consistent framework has not been reached yet.
In addition to the simple f (R) model, it is possible to construct much more complex Lagrangians. As an example, teleparallel gravity considers Lagrangians that are a function of torsion scalar f (T ) [66]. In such models, the Ricci scalar is assumed to vanish, and its role is played by torsion. Furthermore, the particularity of these f (T ) models is that they include torsion as the source of DE without considering a cosmological constant, and they give rise to second order field equations that are easier to solve than the fourth order ones from f (R) gravity [66][67][68][69][70][71][72][73][74][75][76][77]. These models have been widely investigated from fundamental to cosmological scales [76,[78][79][80][81][82][83][84]. Recently, many efforts been devoted to investigating the f (R, T ) models, where the gravitational Lagrangian consists of an arbitrary function of the Ricci scalar R and the torsion scalar T [85][86][87][88]. Another class of interesting models that have been introduced are, for example, f (G), where G is Gauss-Bonnet topological invariant or combinations of these with the Ricci scalar as the f (R, G) [89][90][91][92][93][94][95][96]. Using these models, it is possible to construct viable cosmological models, preventing ghost contributions and contributing to the regularization of the gravitational action.
In this review, we will focus only on f (R) models. We summarize the main achievements, pointing out the difficulties and the future perspectives in probing the f (R)-gravity using LSS datasets. In Section 2, we report the general formalism and the main features of f (R)-gravity models, and also, we focus on two models/approaches that have led to several results in the last few years. In Section 3, we describe two approaches to probe f (R)-gravity using a cluster of galaxies. In Section 4, we review the most important results obtained using N-body simulations to quantitatively describe the physical effect of f (R)-gravity models on the structure formation. In Section 5, we summarize the main approaches to test alternative models using the expansion history of the Universe. In Section 6, we consider the constraints obtained using the CMB power spectrum. Finally in Section 7, we give the main conclusions and the future perspectives for this field.

f(R) Gravity
In metric formalism, the simplest four-dimensional action in f (R) gravity is: where g is the determinant of the metric g µν , and L matter is the standard perfect fluid matter Lagrangian. Varying the action Equation (1) with respect to the metric tensor, we obtain the field equations and its trace: Here, f ,R = df /dR, = g µν ∇ µ ∇ ν is the d'Alembert operator, and T µν = δg µν is the energy momentum tensor of the matter. We can recast the Equation (2) in Einstein form as follows: where we can reinterpret the first term on the right side of the equation as an extra stress-energy tensor contribution T (curv) µν . f (R) gravity is sufficiently general to include all of the basic features of ETGs, but at the same time, it can be easily connected to the observations. Let us consider now two classes of f (R) models, which are designed to satisfy cosmological and Solar System constraints.

Chameleon Models
One way to modify gravity is to introduce a scalar field (φ), which is coupled to the matter components of the Universe. The range and the length of interactions mediated by the scalar field depend on the density of the environment. In high-density configurations, the effective mass of the scalar field is large enough to suppress it and to accommodate all existing constraints from Solar System tests of gravity. Since the mass depends on the local matter density, the cosmological evolution of the scalar field could be an explanation for the acceleration of the Universe. Many environment-dependent screenings have been proposed, such as the chameleon [97], the dilaton [98], the symmetron [99] or the Vainshtein [100,101] mechanisms. One of the most tested mechanism is the chameleon; it must satisfy the following equation [97]: where V is the potential of the scalar field, β is the coupling to the matter, ρ is the matter density and M Pl = √ 8πG is the Planck mass. The chameleon field leads to the fifth force in the form of: f (R) models show a chameleon mechanism having a fixed value of the coupling β = 1 6 . These models have an additional degree of freedom that mediates the fifth force [27]: where φ ∞ describes the efficiency of the screening mechanism when r → ∞. Thus, the f (R) Lagrangian cannot be considered as a completely free function; it must satisfy many different constraints coming from both theoretical and phenomenological side. For example, in a region with high curvature, the first derivative of the Lagrangian must be f ,R > 0, as well as the second derivative (we define: f ,RR = d 2 f (R)/dR 2 ) f ,RR > 0. These conditions will avoid tachyons and ghost solutions. Furthermore, to satisfy, Solar System constraints must be |f ,R0 | << 1 in the present Universe. Two well-studied choices matching these constraints are the Starobinsky [102] and Hu-Sawicki [103] models. In this review, we will focus only on the latter, which is expressed by the following Lagrangian: from which one can obtain: where the Ricci scalar for the flat ΛCDM model is: Here, the termΩ Λ Ω m matches the standard Ω Λ Ω m when c 1 /c 2 2 → 0. Hu and Sawicki [103] have shown that to recover GR results within the Solar System, one must be |f, R0 | < 10 −6 .

Analytical f(R) Gravity Models and Yukawa-Like Gravitational Potentials.
In this subsection, we will illustrate how an analytic f (R) model expandable in Taylor series gives rise to a Yukawa correction to the Newtonian gravitational potential. Phenomenologically, this kind of correction behaves as the so-called chameleon mechanism, becoming negligible at small scales [104,105]. In the case of the post-Newtonian corrections, such a solution leads to redefining the coupling constants in order to fulfill the experimental observations. In order to derive the correction terms to the Newtonian potential coming from an analytic f (R) model, one must consider both the second and the fourth order in the perturbation expansion of the metric. Let us consider the perturbed metric with respect to a Minkowskian background g µν = η µν + h µν . The metric components can be developed as follows: tt (t, r) + g (4) tt (t, r) g rr (t, r) ≃ 1 + g (2) rr (t, r) g θθ (t, r) = r 2 g φφ (t, r) = r 2 sin 2 θ (11) where c = 1 , x 0 = ct → t, and a general spherically-symmetric metric has been considered: where dΩ is the solid angle. The next step is to introduce an analytic Taylor expandable f (R) functions with respect to a certain value R = R 0 : In order to obtain the weak limit, we must introduce Equations (11) and (13) in the field Equations (2) and (3) and to expand the equations up to the orders O(0), O(2) and O(4). This approach permits us to select the Taylor coefficients in Equation (13). We remember that at zero order O(0), the field equations give the condition f 0 = 0, and thus, the solutions at further orders do not depend on this parameter. The field equations at O(2)-order become: tt,r + 8f ,RR0 R (2) ,r − f ,R0 rg (2) tt,rr + 4f ,RR0 rR (2) rr,r + 8f ,RR0 R (2) ,r − f ,R0 rg (2) tt,rr = 0 2f ,R0 g (2) rr − r f ,R0 rR (2) − f ,R0 g (2) tt,r − f ,R0 g (2) rr,r + 4f ,RR0 R (2) ,r + 4f ,RR0 rR (2) ,rr = 0 ,r + rR (2) ,rr = 0 2g (2) rr + r 2g (2) tt,r − rR (2) + 2g (2) rr,r + rg (2) tt,rr = 0 As we can see, the fourth equation in the above system (the trace) gives us a differential equation with respect to the Ricci scalar, which allows solving the system exactly at O(2)-order. Then, we obtain: where 6f ,RR0 , f ,R0 and f ,RR0 . In the limit f (R) → R, we recover the standard Schwarzschild solution. Let us stress that L has the dimension of length, the integration constant δ 0 is dimensionless, while δ 1 (t) has the dimensions of length −1 and δ 2 (t) the dimensions of length −2 . Since the differential equations in the system Equation (14) contain only spatial derivatives, we can fix the functions of time δ i (t) (i = 1, 2) to arbitrarily-constant values. Furthermore, we can set δ 0 to zero, because it is not an essential additive quantity. If we want a general solution of previous Equation (15), we exclude the Yukawa growing mode, obtaining the following relations: where r g = 2MG. Now, we can look for the solution in terms of the gravitational potential from Equation (15), obtaining the following relation: Let us note that the L parameter is related to the effective mass: and it can be interpreted also as an effective length. Being 1 + δ = f ,R0 and δ 1 = − 6GM L 2 δ 1 + δ quasi-constant, the Equation (17) becomes: Here, the first term is the Newtonian-like part of the gravitational potential to point-like mass, and the second term is the Yukawa correction including a scale length, L. If δ = 0, the Newtonian potential is recovered. With this assumption, the new gravitational scale length L could naturally arise and accommodate several phenomena ranging from Solar System to cosmological scales.

Constraining f(R) Gravity Models Using Clusters of Galaxies
Clusters of galaxies are the largest virialized object in the Universe. They are the intermediate step between the galactic and the cosmological scales; thus, any relativistic theory of gravity must be capable of correctly describing their physics. Clusters typically contain a number of galaxies ranging from a few hundreds to one thousand, grouped in a region of ∼2 Mpc, contributing 3% to the total mass of the cluster. A more important component is represented by the baryons residing in a hot inter cluster (IC) gas. Although IC gas is highly rarefied, the electron number density is n e ∼ 10 −4 − 10 −2 cm −3 ; it makes up 12% of the total mass; and it reaches high temperatures ranging from 10 7 to 10 8 K, becoming a strong X-ray source with a luminosity typically ranging L X ∼ 10 43 − 10 45 erg/s.
One of the most promising tools to study clusters of galaxies is the Sunyaev-Zeldovich (SZ) effect [106,107]. CMB photons cross clusters of galaxies, and they are scattered off by free electrons present in the hot IC gas. This interaction produces secondary temperature fluctuations of the CMB power spectrum due to: (1) the thermal motion of the electrons in the gravitational potential well of the cluster (it is named Thermal Sunyaev-Zeldovich (TSZ) ; [106]); (2) the kinematic motion of the cluster as a whole with respect to the CMB rest frame (it is named Kinetic Sunyaev-Zeldovich KSZ ; [107]). In the direction of a cluster (n), the SZ effect is given by: where dτ = σ T n e dl is the optical depth, σ T is the Thomson cross-section, k B is the Boltzmann constant, m e c 2 is the electron annihilation temperature, c is the speed of light, ν is the frequency of the observation and v cl is the peculiar velocity of the cluster. T 0 = 2.725 ± 0.002 K [108] is the current CMB blackbody temperature, and g(ν) is the frequency dependence of the TSZ effect. In the non-relativistic limit, g(x) = xcoth(x/2) − 4, with x = hν/k B T the reduced frequency. The physical description of the TSZ is commonly given by introducing the Comptonization parameter: where P e (r) is the electron pressure profile that must be specified to predict the TSZ anisotropies. Using X-ray observations and numerical simulations, several cluster profiles (n e T e ) have been proposed.
Isothermal β-model fits the X-ray emitting region of the clusters of galaxies well [109,110], with the electron density given by: where n e,0 is the central electron density and r c is the core radius. Those parameters, together with the electron temperature and the slope β, are determined from observations. Using the X-ray surface brightness of clusters, the slope value ranges in the interval [0.6 − 0.8] [111].
The isothermal β-model over-predicts the TSZ effect in the outskirts of the cluster of galaxies [112]. Recently, a phenomenological parametrization of the electron pressure profile, derived from the numerical simulations, has been proposed [113,114]. The functional form of Universal pressure profile is: where x = r/r 500 , and r 500 is the radius at which the mean overdensity of the cluster is 500-times the critical density of the Universe at the same redshift. Then, c 500 is the concentration parameter at r 500 . The model parameters were constrained using X-ray data [114,115] or CMB data [116]; their best fit values are quoted in Table 1. Table 1. Parameters of universal pressure profile fitted by different groups using X-ray [114,115] and CMB data [116]. Since TSZ does not depend on redshift in the ΛCDM model, clusters are a very important laboratory to test cosmology (see [117][118][119][120] and the references within), even to test the fundamental pillars of the Big Bang scenario [121,122].
To study if these discrepancies are due to the theoretical modeling and, at the same time, to constrain/rule out ETGs at cluster scales, the pressure profile of a cluster of galaxies has been considered. On the one hand, it is interesting to study models that can explain the cluster without requiring any dark component [50]. On the other hand, models constructed to mimic the ΛCDM expansion history by replacing the DE with higher order terms in the gravitational Lagrangian must also describe the emergence of the LSS and the dynamical properties of clusters [51][52][53].

Pressure Profile from Yukawa-Like Gravitational Potential.
One approach to test ETGs using the TSZ pressure profile is focused on the f (R) Lagrangian that is expandable in Taylor's series (Equation (13) in Section 2.2), where the Yukawa-like correction to the Newtonian potential, Equation (19), could be interpreted as the dynamical effect of DM in clusters of galaxies [50]. Making the hypothesis that: the gas is in hydrostatic equilibrium within the modified gravitational potential well (without any DM term): the gas follows a polytropic equation of state: the Equations (19), (24) and (25) can be numerically integrated to compute the pressure profiles by closing the system with the equation for the mass conservation: Thus, the pressure profile will be a function of the two extra gravitational parameters (δ, L) and the polytropic index γ. Let us remark that the density ρ(r) includes only baryonic matter, without resorting to DM particles.

Data and Results.
To constrain the model, the Planck (all Planck data are publicly available and can be downloaded from http://www.cosmos.esa.int/web/planck) 2013 data and an X-ray clusters catalog [136] have been used. In Figure 1, the predicted Comptonization parameter for the universal pressure profile with the parameters given in Table 1  The model is particularized for the Coma cluster, located at redshift z = 0.023, with core radius r c = 0.25 Mpc, X-ray temperature T X = 6.48 keV and electron density n e,0 = 3860 m −3 . The plot shows that the higher order terms in the gravitational f (R)-Lagrangian could potentially explain the cluster without accounting for DM contribution.
The SZ emission was measured on the Spectral Matching Independent Component Analysis (SMICA) map (SMICA is a foreground-cleaned map; it has been constructed using a component separation method by combining the data at all frequencies [137]; the SMICA map has a 5 ′ resolution) at the locations of all X-ray selected clusters, and then, it was compared with the model prediction. The temperature anisotropies were averaged over rings of width θ 500 /2, where θ 500 is the angular scale subtended by the r 500 . The error associated with each data point was computed carrying out 1000 random simulations and evaluating the profile for each simulation. The analysis was performed by computing the likelihood function log L = −χ 2 /2 as: where N is the number of data points, and p = (δ, L, γ). In Equation (27), d(x i ) are the data, and C i,j is the correlation matrix between the average temperature anisotropy on the discs and rings.
Two parameterizations have been tested: (A) L = ζr 500 , the scale length is different for each; clusters scale linearly with r 500 ; (B) L has the same value for all of the clusters. The models were constructing choosing "physical" flat priors on the parameters. Looking at the Yukawa-potential in Equation (19), if δ < −1, it becomes repulsive, and if δ = −1, it diverges; while the polytropic index γ was varied in the range corresponding to an isothermal and an adiabatic state of a monoatomic gas, respectively. The priors are summarized in Table 2.   Figures 2 and 3 show the 2D contours at the 68% and 95% confidence levels of the marginalized likelihoods for both Parameterizations (A) and (B), respectively. Although the contours are opened and only upper limits can be obtained, the value δ = 0 is always excluded at more than the 95% confidence level. Therefore, Newtonian potential without DM cannot fit cluster pressure profile, and thus, either DM or modified gravity must be the right description of the dynamical properties of the cluster of galaxies. The upper limits are summarized in Table 3.    Figure 2. This parameterization also provides opened contours. Therefore, also in this case, one can only give upper limits on the parameters and can not distinguish which parameterization is the best one.

Chameleon Gravity: Hydrostatic and Weak Lensing Mass Profile of Galaxy Cluster.
A second approach in using clusters of galaxies to test ETGs is focused on the chameleon gravity models that introduce a scalar field non-minimally coupled to the matter components. This coupling gives rise to the fifth force, which is tightly constrained at Solar System scales (see Section 2.1). The main difference with the other approach discussed above is that the chameleon gravity only eliminates DE, preserving the DM component as the source of the gravitational field needed to explain the emergence of the LSS.
Because of the screening mechanism, the effect of this force is negligible at small scales, and thus, it does not affect the density distribution of the IC gas in the cores of galaxy clusters, but its effect may not be totally screened in the outskirts. As a consequence, the gas distribution will be more compact in the outskirts under the effect of an additional pressure term that balances the effect of the fifth force. Thus, the hydrostatic mass of a cluster should be different from the mass obtained from weak gravitational lensing that depends only on the distribution of DM along the line of sight and does not assume hydrostatic equilibrium.
The model studied in [51,52] assumes the spherically-symmetric distribution of IC gas and DM. The IC gas is in hydrostatic equilibrium within the potential well generated by the DM distribution.
Departures from hydrostatic equilibrium are parameterized, introducing a non-thermal term in the pressure, while the chameleon field directly contributes to the total amount of mass. The model is compared to lensing observations used to estimate the mass, with X-ray and SZ observations used to reconstruct the gas pressure profile [51][52][53].
Let us start writing down the equation of the hydrostatic equilibrium for the IC gas component: where M(< r) is the total mass enclosed in a sphere of radius r, ρ gas is the gas density distribution within the sphere and P tot is the total gas pressure that can be recast as: including both thermal (P th (r)) and non-thermal (P non−th (r)) pressure, and the term due to the chameleon field (P φ (r)). Thus, the total mass in Equation (28) is: By definition, each term can be expressed as: M φ (r) is the mass associated with the chameleon field, and M non−th represents the departure from the hydrostatic equilibrium mass M th . Introducing the equation of state of the IC gas, P th (r) = kn gas (r)T gas (r) the Equation (31) can be re-written as where the identity ρ gas (r) = µm p n gas (r) has been considered with the mean molecular weight µ and the proton mass m p . According to hydrodynamical simulations of the ΛCDM model [138,139], the non-thermal pressure term can be modeled as: n gas (r)kT gas (r) with: where α nt , β nt , n nt and n M are constants. Therefore, the non-thermal mass term is: This term is modeled without considering the chameleon field, and it has been shown to reproduce accurately the non-thermal contribution, even in f (R)-gravity. Nevertheless, a more accurate procedure would require modeling this term from hydrodynamical simulations under chameleon gravity. Finally, the total mass inferred from hydrostatic equilibrium has to be equal to the the total mass as inferred from weak lensing (WL): The DM density distribution is equally well described in the f (R)-model and in Newtonian gravity using the Navarro-Frenk-White (NFW) profile [140]. This profile has the following functional form [141]: and it has been constructed using N-body simulations of the ΛCDM model. Thus, the DM mass within the radius r is: that is equivalent to the WL mass if the DM component dominates over the baryonic one. The model has been tested using only COMA cluster [52] and using a sample of 58 X-ray selected cluster associated also with weak lensing measurements [53].

Data and Results
The chameleon model is essentially described by the following parameters (β, φ ∞ ) in Equations (5) and (7). In order to normalize them, those parameter have been replaced by: and constrained by carrying out the Monte Carlo Markov chain (MCMC) method using X-ray, SZ and WL observations. Another important parameter for studying clusters of galaxies is the critical radius: where ρ s is the density at this radius. It represents the distance from the DM halo center where the screening mechanism acts and it is not possible to distinguish between chameleon gravity and GR [51].
In [52], the authors have constrained chameleon gravity only using the Coma cluster (z = 0.0231). For this cluster, they implement the MCMC method using the X-ray temperature profile reported by the X-ray Multi-Mirror Mission (XMM) -Newton [142] and Suzaku [143] for the inner and outer region, respectively. Then, the X-ray surface brightness profile was taken from XMM-Newton [144] and the SZ pressure profile measured by Planck [131]. For the WL counterpart, the measurements were reported in [145]. More recently, in [53], the authors have constrained the chameleon model using publicly-available weak lensing data provided by the Canada France Hawaii Lensing Survey (CFHTLenS; [146]), and X-ray data taken from the XMM-Newton archive. Their sample includes 58 clusters spanning the redshift range z = [0.1 − 1.2], with measured X-ray temperatures T x = [0. 2 − 8] keV. For this sample of clusters, they also carried out an MCMC analysis. In Figure 4 are shown the two-dimensional contours for the parameters β 2 and φ ∞,2 as generated by [53]. The red and blue lines are the 95% and 99% confidence limits from [52]. Those results represent the best constraints on chameleon gravity using clusters. On the one hand, at low values of β, the departure from GR is too small to be detected with these observational errors. On the other hand, GR gravity is recovered within the critical radius r crit , and large values of β lead to a lower value for φ ∞ that could keep r crit in the whole cluster region. Although the analysis from [53] using the profiles outside of the critical radius of the cluster could better constrain large values of β, the analysis carried out by [52] also used SZ measurements and could discriminate better between GR and chameleon gravity at lower values of β. In both [52,53], it also constrained the term f ,R0 from Equation (7). The result provides an upper limit today |f ,R0 | < 6 × 10 −5 at a 95% of confidence level. Figure 4. The confidence contours for the renormalized parameters (φ ∞,2 , β 2 ) The 95% and 99% confidence levels are plotted in light gray and medium gray, respectively. The results are from [53]. The 95% and 99% confidence contours from [52] are over-plotted also in red and blue, respectively. The vertical line corresponds to |f ,R0 | < 6 × 10 −5 .

N-Body Hydrodynamical Simulations in f(R) Gravity
All chameleon-like theories of gravity need a specific screening mechanisms (Section 2.1). The effect of all such mechanisms is reflected on the non-linear perturbations of the matter distribution. Therefore, some of the most optimal tools to provide a full description of these mechanisms are modifications of standard N-body algorithms [147][148][149][150][151][152][153].
The aims were related to the study of the non-linear effects on the emergence of the LSS, such as the matter power spectrum [148,150,154], the redshift-space distortions [55] and the physical and statistical properties of the cold DM halos and voids [56,140,[155][156][157][158][159][160]. However, the same scales and structures that could provide a powerful tool for testing the theory of gravitation are also seriously affected by non-linear astrophysical processes [161][162][163][164][165][166]. In the few last years, a very powerful code to investigate the combined effects of astrophysical processes and alternative theories of gravity has been developed [57]. The code is named Modified Gravity (MG)-GADGET. It was used to study the Hu-Sawicki model [103], pointing out the existence of an important degeneracy between modified gravity models and the uncertainties in the baryonic physics particularly, related to AGN feedback [57]. Afterwards, the code was combined with a massive neutrinos algorithm [167], to study the cosmic degeneracy between modified gravity and massive neutrinos. It was found that the deviations from the matter power spectrum in the ΛCDM model are, at most, of the order of 10% at z = 0 taking f (R) models with f ,R0 = −1 × 10 −4 , and the total neutrinos mass Σ i m ν i = 0.4 eV [168]. Using again the Hu-Sawicki parameterization, the effects of the model on the IC gas and its physical properties have been studied. It was found that the DM velocity dispersions in the f (R) model are ∼4/3 larger than in ΛCDM, when low mass halos are considered. Nevertheless, the mass of the objects that are affected by the chameleon field depend on the value of the background field: for |f ,R0 | = 10 −4 the field is never screened, nor for massive clusters; for |f ,R0 | = 10 −5 , the chameleon field is screened for mass above ∼10 14.5 M⊙; and for |f ,R0 | = 10 −6 , the Hu-Sawicki model does not produce any effect in the range of mass explored in [169]. Finally, MD-GADGET and the Hu-Sawicki model have been also used to study the synthetic Lyman-α absorption spectra using simulations that include star formation and cooling effects. The discrepancies with the simulations based on the ΛCDM model were at most 10% with the background field fixed at |f ,R0 | = 10 −4 [58].
All tests carried out using N-body simulations have provided several indications about the presence of a degeneracy between the effects of the complex baryonic physics (highly non-linear) and the effects of the modifications of the theory of gravitation. Even with the high accuracy of future data, without a full understanding of the non-linear astrophysical processes, it would be not possible to break the degeneracy and to constrain or rule out ETGs. Therefore, there is the need to fully explore this degeneracy with N-body simulations to know the level of bias introduced.

Constraining the Expansion History of the Universe in f(R) Gravity
Cosmological models can be tested looking at background evolution and also at the growth of the cosmic structures. Testing the f (R) model at the background level allows finding the so-called "realistic" models that have the radiation, matter and DE eras very close to the ΛCDM ones. Although many realistic f (R) models have been constructed, it should be considered that different models with a comparable late-time accelerated expansion and similar background densities' evolution could differ in the growing of the matter density perturbations [170]. For example, it is well known that the effects at the background level of considering a collisional matter component filling the Universe are negligible. Nevertheless, non-collisional and collisional matter affect in different ways the evolution of matter perturbations, opening the possibility of distinguishing them [171]. Another problem is represented by the need to break the degeneracies between the evolving DE models and modified gravity in order to confirm/rule out models [172]. Thus, the growth factor (GF) data become extremely important to test modified gravity and to find their signatures. However, the available data have less accuracy than the other datasets, such as SNeIa, Baryonic Acoustic Oscillation (BAO) and H(z); thus, it is not possible to fully constrain f (R)-models. This scenario will totally change with the forthcoming Euclid observations. That dataset will have an unprecedented accuracy, and it will allow constraining deviations from the standard model at a few percent level [54].
Let us start by the flat Friedman-Lemaitre-Robertson-Walker metric (FLRW), where a(t) is the rate expansion of the Universe. The Ricci scalar is related to the Hubble function, by the following expression: and using the non-relativistic matter and radiation approximation, the following conservation laws are satisfied:ρ m + 3Hρ m = 0 (47) Using Equations (2), (3), (45) and (46), the modified Friedman equations are obtained [173]: To compute the the evolution of the matter density perturbations, δ m ≡ δρm ρm , one can mainly follow two different approach: • Sub-horizon approximation: The evolution of the matter density perturbations can be obtained solving the following equation [174]: where k is the comoving wave number, and G eff (a, k) is the effective gravitational constant. The explicit forms of G ef f depend on the ETGs model. In the case of f (R) gravity, it is given by [175]: By definition, the growth factor f g is: thus, Equation (51) becomes: Since Equation (51) does not give raise to analytical solutions, it is customary to introduce an efficient parameterization of the matter perturbations. The most common one is the so-called γ-parameterization: is obtained for the value γ ≈ 0.545. Therefore, any departures from this value could indicate the need to use a different model. More in general, the γ function could evolve with the redshift in the case of modified gravity models [176,178,179]: This approach has been widely used to constrain the f (R) gravity, pointing out that: the current data have not enough constraining power to discriminate between different modified gravity models and the ΛCDM model; the so-called growth index is the most efficient way to provide constraints on the model parameters; and the γ-parameterization represents the best way to measure deviation from the standard model [176,177].
• Dynamical variables: Starting from Equations (49) and (50), it is possible to reach a model-independent description of the evolution of the matter density perturbations [54,180,181]. Defining the following dimensionless variables: the density parameters and the effective equation of state of the system can be written as: The background evolution is given by: where the prime indicates the derivative with respect to the scale factor d/d ln a, and m is a parameter related to the theory that indicates how much a specific model is close to the ΛCDM (or how much it deviates from it). This parameter is defined as: and by definition, one also finds the following identity: Equations (65) to (70) represent a closed system that can be integrated numerically, obtaining the evolution of the background densities perturbations. We just need to choose an f (R)-model to compute the two parameters m and r that we need to integrate the system. From the analysis of the critical points, we have several conditions that an f (R) model has to respect to be cosmologically viable [54,180]. In general, only models with m ≥ 0 and very close to the ΛCDM model (f (R) = R − 2Λ) are cosmologically viable. Thus, the equations for the matter density perturbations are [181,182]: where δf ,R ≡ δf ,R /f ,R and x 5 ≡ aH satisfies: The evolution at all scales of the growth of the matter density perturbation, δ m , could be obtained via numerical integration. However, a more straightforward approach is to assume that the growth rate, f g , is a function of time, but not of scale [183][184][185][186][187], and to use the γ-parameterization in Equations (55) and (56).
The detection of possible deviation from ΛCDM model has been studied fitting the γ-parameterization to the forthcoming Euclid dataset [54]. Assuming the ΛCDM model as the reference model, it has been found that: the forthcoming data will be able to constrain the parameter γ and w with an accuracy of 4% and 2% if they do not depend on redshift (standard evolution) and to clearly distinguish the ΛCDM from alternative scenarios at more than 2σ; see Figure 5.

Testing Gravity Using the Cosmic Microwave Background Data
The CMB temperature anisotropies play a key rule to constrain the cosmological model. They are mostly generated at the last-scattering surface, and thus, they provide a way to probe the early Universe [13,17,133]. However, the lack of a full-consistent theoretical framework and the fine tuning problem related to the cosmological constant, Λ, have encouraged studying wide classes of alternative DE or modified gravity models. The most difficult challenge of modern cosmology is to use current and forthcoming datasets to discriminate between those alternative visions of the Universe [54,[188][189][190][191]. Many of those models mimic the ΛCDM at background evolution, but differ in the perturbations, and they could produce several effects on the CMB power spectrum. As an example, different models could have different expansion histories that shift the location of the peaks [192], and could affect the gravitational potential at late time, changing the integrated Sachs-Wolfe (ISW) effect [193,194]. Moreover, modification of the gravitational theory should be reflected in modifications of the lensing potential [195,196], the growth of the structures [197][198][199][200] and also the amplitude of the primordial B-modes [201,202].
Modified gravity models affect the evolution of the Universe at both the background and the perturbation [203][204][205][206]. There are mainly two different approaches used to constrain alternative models of gravity: the parametrized modified gravity approach and the effective field theory.
The parametrized modified gravity approach constructs functions probing the geometry of space-time and the growth of perturbations. Starting from a spatially-flat FLRW metric, the perturbed line element in the Newtonian gauge is given by: where τ is the conformal time. Ψ and Φ are the two scalar gravitational potentials that are a function of the scale Ψ = Ψ(k, a) and Φ = Φ(k, a) and are related to the energy-momentum tensor, T µν , through Einstein's equations. Those potentials can be modeled to fix all degrees of freedom in order to match observational data. Moreover, the difference of those two potentials Φ − Ψ, which is the anisotropic stress, has to be equal to zero in GR. Any departure from this value would indicate a departure from GR [207,208]. One of the most common parameterization is implemented in the Modification of Growth with Code for Anisotropies in the Microwave Background (the code is publicly available at http://www.sfu.ca/ aha25/MGCAMB.html) (MGCAMB) [209,210], and it considers the following: to parameterize the Poisson and the anisotropic stress equations using two scale and time-dependent functions µ(k, a) and γ(k, a). Those two functions are determined by choosing the particular model. For example, some classes of f (R) models lead to [211]: where the parameters β i are dimensionless and represent the couplings, and the parameters λ i are lengths. Those functional forms have been used in early studies to forecast errors of the parameters for the Dark Energy Survey (www.darkenergysurvey.org/) (DES, [212]) and Large Synoptic Survey Telescope (http://www.lsst.org/lsst/) (LSST, [213]) surveys [210]. The main limitation of this approach is the use of the quasi-static regime assuming that the scales at which one is interested are still linear and smaller than the horizon. Furthermore, the time derivatives are neglected. Effective field theory (EFT) [214,215] describes the whole range of the scalar field theories. The Lagrangian preserves the isotropy and homogeneity of the Universe at the background, providing a more general approach without assuming the quasi-static regime, but adding much more free parameters that must be constrained. The f (R)-models are a sub-class of those theories. The action of EFT reads: + m 2 2 (τ ) (g µν + n µ n ν ) ∂ µ a 2 g 00 ∂ ν a 2 g 00 where δR (3) is its spatial perturbation, K µ ν is the extrinsic curvature, m 0 is the bare (reduced) Planck mass and S m is the matter part of the action, including all fluid components, such as baryons, cold DM, radiation, and neutrinos, but it does not include DE. The nine time-dependent functions [216] {Ω, c, Λ,M 3 1 ,M 4 2 ,M 2 3 , M 4 2 ,M 2 , m 2 2 }, have to be fixed to specify the theory. The EFT has been implemented in the publicly-available EFTCAMBcode (http://www.lorentz.leidenuniv.nl/hu/codes/, Version 1.1, October 2014.) [59,217], that numerically solves both the background and perturbation equations. The code has been used to constrain massive neutrinos in f (R) gravity [59].
Planck collaboration has used in both approaches to fit f (R) models to investigate their effect combining early time datasets, such as CMB, and cosmological datasets at much lower redshift [14].
In general, f (R) models require setting the initial condition: and the boundary condition: High values of B 0 , its present value, change the ISW effect and the CMB lensing potential [211,218,219], Marchini2013. Their analysis has shown the presence of a degeneracy between the optical depth τ and B 0 that can be broken adding other probes on the structure formation, such as weak lensing, CMB lensing or redshift-space distortions (see Figure 6). . The 68% and 95% contour plots for the two parameters, {Log 10 (B 0 ), τ }. There is a degeneracy between the two parameters for Planck temperature power spectrum (TT) + BSH (the combination of BAO, SNIa, and H 0 datasets). Adding lensing will break the degeneracy between the two. Here, Planck indicates Planck TT.
The same study carried out using MGCAMB code led to the same results within the uncertainties, improving the previous constraint on B 0 [220].
Nevertheless, both EFTCAMB and MGCAMB assume a specific background cosmology. MGCAMB code assumes for the background evolution the ΛCDM or wCDM, while EFTCAMB assumes a fixed background evolution, and the CMB power spectrum is calculated for an f (R) model re-constructed by the fixed expansion history. A more recent code, named f (R) Code for Anisotropies in the Microwave Background (FRCAMB) [221], allows reconstructing the CMB power spectrum for any f (R) model by specifying the first and second derivatives of the Lagrangian. Although the results obtained are consistent with the one from EFTCAMB and MGCAMB, the code has the advantage of being designed for any f (R) model at both the background and perturbation level.

Discussion and Future Perspectives
In this review, we have outlined the possibilities offered by the current and the forthcoming astrophysical and cosmological dataset to constrain or rule out the f (R)-models. It is mandatory to search for the answer to one of the fundamental questions in modern cosmology: is GR the effective theory of gravitation? This question comes from the evidence that GR is not enough to fully explain the cosmological evolution of the Universe and the emergence of the clustered structures, and it needs the addition of two unknown components. The difficulties in identifying the DM and DE components at a fundamental level have opened the possibility that the answer resides in modifying the theory of gravity. Nevertheless, having an alternative demands testing it study the physical phenomena at much smaller scales than the cosmological ones. To this effect, much analysis has been carried out to ensure that modified gravity could recover the tight constraints at the Solar System scale. Models that survive these probes have been used to test galactic and extragalactic scales. Clusters of galaxies provide an excellent laboratory to probe the different features of the different modified gravity models. The variety of the phenomenology related to the cluster physics allows one to use SZ, X-ray and WL observations in order to constrain f (R)-models. Many efforts have been also devoted to quantitatively describing the physical effects of alternative cosmological scenarios. Those efforts have required developing new algorithms for hydrodynamical N-body simulations, to study the degeneracy between baryonic physics and modifications of gravity and also new Boltzmann codes to predict the CMB power spectrum and to take advantage of the Planck data to constrain modified gravity models.
Nevertheless, a self-consistent picture has not been reached yet. The current datasets do not have enough statistical power to discriminate between modified gravity models and standard cosmology, and when the needed precision is present, several degeneracies have been found. Next generation experiments, such as the LSST [213] and the Wide-Field Infrared Survey Telescope (WFIRST) (http://wfirst.gsfc.nasa.gov/) [222], will undertake large imaging surveys of the sky, while other surveys, such as Euclid [223] and the Big-Baryon Oscillation Spectroscopic Survey (BigBOSS) [224], will make spectroscopic measurements of galaxies. The combination of all forthcoming datasets will hopefully put strong constraints on the evolution of the Universe and the emergence of the LSS, allowing us to discriminate between the concordance ΛCDM model and its alternatives.