{'fraction_non_alphanumeric': 0.05099111244258039, 'fraction_numerical': 0.023479628520071898, 'mean_word_length': 3.8778007793472966, 'pattern_counts': {'":': 0, '<': 63, '<?xml version=': 0, '>': 40, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Density perturbations in f(R) gravity theories in metric and Palatini formalisms

We make a detailed study of matter density perturbations in both metric and Palatini formalisms in theories whose Lagrangian density is a general function, f(R), of the Ricci scalar. We derive these equations in a number of gauges. We show that for viable models that satisfy cosmological and local gravity constraints (LGC), matter perturbation equations derived under a sub-horizon approximation are valid even for super-Hubble scales provided the oscillating mode (scalaron) does not dominate over the matter-induced mode. Such approximate equations are especially reliable in the Palatini formalism because of the absence of scalarons. Using these equations we make a comparative study of the behaviour of density perturbations as well as gravitational potentials for a number of classes of theories. In the metric formalism the parameter m=Rf_{,RR}/f_{,R} characterising the deviation from the Lambda CDM model is constrained to be very small during the matter era in order to ensure compatibility with LGC, but the models in which m grows to the order of 10^{-1} around the present epoch can be allowed. These models also suffer from an additional fine tuning due to the presence of scalaron modes which are absent in the Palatini case. In Palatini formalism LGC and background cosmological constraints provide only weak bounds on |m| by constraining it to be smaller than ~ 0.1. This is in contrast to matter density perturbations which, on galactic scales, place far more stringent constraints on the present deviation parameter m of the order of |m|<10^{-5} - 10^{-4}. This is due to the peculiar evolution of matter perturbations in the Palatini case which exhibits a rapid growth or a damped oscillation depending on the sign of m.


I. INTRODUCTION
Recent high-precision observations by the Wilkinson Microwave Anisotropy Probe (WMAP) [1] together with high redshift supernovae surveys [2], observations of large scale structure [3] and baryon acoustic oscillations [4] have provided strong evidence that the Universe is at present undergoing a phase of accelerated expansion [5,6]. Such an accelerating phase poses a serious problem for cosmology since it is difficult to explain within the usual general relativistic framework.
Phenomenologically, the simplest way to generate such an accelerating phase is through the addition of a cosmological constant to the Einstein's field equations. At a more fundamental level, however, the problem is how to account for such a constant within a candidate theory of quantum gravity. This has motivated a large number of alternative models (see Ref. [6] for a recent review). These models can mostly be divided into two broad groups: those involving an exotic matter source and those modifying the gravitational sector of the theory. An important subset of the latter are the so called generalised gravity theories, based on non-linear lagrangians of the form f (R), where f is a general differentiable function of the Ricci scalar R. Such modifications to the (linear) Einstein-Hilbert action typically arise in effective actions derived from string/M-theory [7,8,9,10].
There are two approaches that may be taken in studying generalised f (R) theories of gravity, depending upon the choice of the independent fields with respect to which the action is varied. In the first ('metric' approach) only variations with respect to the metric are considered, whereas in the second ('Palatini' approach) the action is varied with respect to both the metric and connection. Both approaches result in both metric and Palatini formalisms provided by background cosmological evolution and local gravity constraints. In Section III we give a brief account of metric perturbations and matter density perturbations which are applicable to both metric and Palatini formalisms. In Sections IV and V we derive the evolution equations for matter density perturbations in metric and Palatini formalisms respectively and in each case discuss the constraints they provide by considering a number of theories that are viable according to background cosmological and local gravity constraints. Finally we conclude in Section VI.
In what follows we shall use units such that 8πG = 1, but we restore bare gravitational constant G when it is needed.

II. LOCAL GRAVITY CONSTRAINTS AND VIABLE BACKGROUND COSMOLOGICAL EVOLUTION
The classes of f (R) general gravity theories we consider have actions of the form: where f is a general differentiable function of the Ricci scalar R and S m corresponds to a matter action, which depends upon the metric g µν and matter fields ψ m . The Ricci scalar R is defined by R = g µν R µν , where the Ricci tensor R µν is In the case of the metric formalism, the connections Γ α βγ are the usual metric connections defined in terms of the metric tensor g µν . In this case the field equations are obtained by varying the action (1) with respect to g µν to give where F = ∂f /∂R. Note that Eq. (4) corresponds to a trace of Eq. (3). The third and fourth terms on the left hand side of Eq. (3) come from the variation of R µν with respect to g µν . In the case of the Palatini formalism Γ α βγ and g µν are treated as independent variables. Varying the action (1) with respect to g µν gives where R µν (Γ) is the Ricci tensor corresponding to the connections Γ α βγ , which is in general different from the Ricci tensor corresponding to the metric connections R µν (g). Taking the trace of this equation, we obtain where R(T ) = g µν R µν (Γ) is directly related to T and is different from the Ricci scalar R(g) = g µν R µν (g) in the metric case. Taking the variation with respect to the connection, and using Eq. (5), we find It is useful to express the generalised theories based on action (1) in terms of generalised scalar-tensor theories of Brans-Dicke type with a scalar field φ and the corresponding potential V (φ) thus [38]: where ω BD is the Brans-Dicke parameter. Note that the original Brans-Dicke theory [39] corresponds to V (φ) = 0. To see the correspondence of the action (8) with (1), we recall that the variation of (8) with respect to g µν and φ leads to the following equations Now consider the following correspondence where R = R(g) in the metric case and R = R(T ) in the Palatini case. Comparing Eqs. (9)-(10) with Eqs. (3)- (4), it can readily be seen that the f (R) theories in the metric formalism correspond to the above generalised Brans-Dicke-type theories with ω BD = 0. Similarly the f (R) theories in the Palatini formalism correspond to the Brans-Dicke-type theories with ω BD = −3/2. The Newtonian effective gravitational constant in the Brans-Dicke theory (8) can be derived under a weak-field approximation by considering a spherically symmetric body with a mass M ⊙ , constant density ρ and a radius r ⊙ and a vanishing density (ρ = 0) outside the body. Decomposing the field φ into background and perturbation parts (φ = φ 0 + δφ) and using a linear perturbation theory in the Minkowski background with a perturbation h µν , the effective gravitational constant is given by [13]: Here ℓ is a distance from the center of the body and the scalar field mass squared is defined by whereφ is a local field in Minkowski spacetime. We should emphasise here that the expression (12) is only valid subject to M r ⊙ ≪ 1 [14,15]. If the scalar field mass M is large, it can happen that the condition for the linear perturbation theory (δφ ≪ φ 0 ) becomes invalid. Moreover, this validity depends on the distribution of scalar field mass inside and outside the body. When the mass in the region ℓ < r ⊙ is much larger than the corresponding mass in the region ℓ > r ⊙ , a thin-shell can be formed inside the body to satisfy local gravity constraints through a chameleon mechanism [40]. The formation of the thin-shell occurs in a non-linear region in which the above linear result (12) ceases to be valid [14]. An important point to note here is that the Palatini case, corresponding to ω BD = −3/2, is rather special in a number of fundamental ways. For example, the φ field kinetic term in Eq. (10) vanishes in this case, whereas it is non-zero in the metric case with ω BD = 0. As we shall see below this has the important consequence that the oscillatory scalaron mode is absent in f (R) theories based on the Palatini formalism, whereas it is present in all other models with ω BD = −3/2, including f (R) theories based on the metric formalism.
Also the scalar field mass M blows up for finite potential-dependent terms in the parenthesis of Eq. (13) as ω BD approaches −3/2. For the theories with ω BD = −3/2, the scalar field mediates a fifth force with an interaction range 1/M because the field has a kinetic term in Eq. (10). In the Palatini formalism the field lacks a kinetic term in Eq. (10), which implies that the notion of the usual interaction range determined by mass M does not hold. Thus the Palatini case should be treated separately compared to other theories.
In the usual Brans-Dicke theory with V (φ) = 0 and ω BD = −3/2, the mass M vanishes because the field φ propagates freely. Then the Yukawa-correction term e −Mℓ in Eq. (12) becomes one, in which case the Brans-Dicke parameter ω BD is constrained by local gravity experiments to be larger than 40000 [41]. This constraint does not, however, necessarily apply to f (R) gravity theories, because the presence of the potential V (φ) can make such theories compatible with local gravity constraints under certain conditions (which we shall see later).
Before proceeding to consider matter density perturbations and the constraints arising from them, we shall first briefly review the cosmological and local gravity constraints for f (R) theories in both formalisms. The constraints from background dynamics have been studied in both metric [11,12,20] and Palatini [23,34] formalisms. To discuss these constraints, it is useful to introduce the following dimensionless parameters [11] In the metric case insight into the cosmological dynamics can be obtained by considering trajectories in the (r, m) plane. In the Palatini case it is more convenient to use other dimensionless variables for the background dynamics [23], but the deviation parameter m is still important in order to characterise the deviation from the ΛCDM model.

A. Metric formalism
We first review the cosmological viability in the metric case and then proceed to consider the LGC. In the (r, m) plane the matter point corresponds to P M : (r, m) ≈ (−1, 0). In order to have a saddle matter era followed by a late-time acceleration, we require the following conditions [11] De-Sitter fixed points P A lie on the line r = −2. They are stable provided that As long as conditions (15) and (16) are satisfied, we can realize a successful cosmological sequence (radiation, matter, and de-Sitter epochs).
There are a number of models in the recent literature that satisfy the above cosmological constraints. [32], and (ii) f (R) = R − αR β with α > 0 and 0 < β < 1 [12]. For these models, the parameters m and r satisfy the relation m = C(−r − 1), where C is a positive constant in the neighbourhood of r = −1. Using observational constraints on the background dynamics from SN Ia and the sound horizon of CMB, the parameter m has been shown to be constrained to be m < O(0.1) [12].
If the information from LGC is also included the constraints on the model parameters become very strong. The usual procedure to determine the local gravity constraints for f (R) theories is to consider their Brans-Dicke representations (11) and expand the equations of motion around a background Minkowski metric [13]. Since the metric formalism corresponds to ω BD = 0, the scalar field mass defined in Eq. (13) is given by where we have used the relation (11). If M 2 < 0 the Yukawa correction e −Mℓ is replaced by an oscillating function cos(|M |ℓ), but this case is excluded by the experimental requirement that γ ≈ 1. Hence the mass squared M 2 is required to be positive. Clearly we require the mass M to be heavy for consistency with local gravity experiments. In that case, however, the effective gravitational constant (12) obtained under the linear approximation ceases to be valid. As was already mentioned above, a thin-shell begins to form through a chameleon effect in this non-linear regime. To consider this chameleon effect in f (R) gravity, it is convenient to introduce a new conformally related metric and a scalar field [42]:g µν = φg µν , ϕ = 3/2 ln φ .
Then the action (1) in the Einstein frame takes the form: where the coupling β in f (R) models and the potential U are given by In a spherically symmetric setting with an energy density ρ (= −T ), the field ϕ satisfies the following equation [43] wherel is the distance from the center of symmetry in the Einstein frame and Here we have introduced an energy density ρ * ≡ e 3βϕ ρ, which is conserved in the Einstein frame [40].
As an example, let us consider a spherically symmetric body that has an energy density ρ * = ρ * A inside the body (l <r ⊙ ≡ e −βϕ r ⊙ ). Let the energy density outside the body (l >r ⊙ ) be given by ρ * = ρ * B , which is much smaller than ρ * A . Then the effective potential (22) has two minima at ϕ = ϕ A and ϕ = ϕ B satisfying the relations The effective mass at the potential minima are given by m 2 , which comes from the condition ρ * A ≫ ρ * B . Imposing appropriate boundary conditions atl = 0 andl =r ⊙ , the solution to Eq. (21) can be approximated by [14,40,43] where and Φ ⊙ = GM ⊙ /r ⊙ . A thin-shell is developed under the condition ∆r ⊙ /r ⊙ ≪ 1. In this case the effective coupling |β eff | becomes much smaller than |β| so that the models can be consistent with local gravity constraints [40]. Models that can satisfy these conditions have recently been proposed by (i) Hu & Sawicki [27] and (ii) Starobinsky [28]: where n, λ and R c are positive constants. In both models the cosmological constant disappears in a flat spacetime, i.e., f (R = 0) = 0. Note that R c is roughly of the order of the present cosmological Ricci scalar R 0 for n = O(1) and λ = O(1). In high curvature regimes R ≫ R c these models behave as with where C is a positive constant, and m and r are defined in Eq. (14). Thus they are very close to the ΛCDM model with suppressed values of m during matter and radiation eras (r ≃ −1).
In the regimes R ≫ R c one can show that the term |ϕ B − ϕ A | in Eq. (26) is of the order of m(R B ) for n = O(1), where R B is the Ricci scalar in the neighbourhood of ϕ B (which is generally much larger than R c in an environment where local tests of gravity are carried out). Hence the thin-shell is developed under the condition This can be regarded as a criterion for the compatibility with local gravity constraints. In the case of the earth, the condition (31) corresponds to m(R B ) ≪ Φ ⊙ ∼ 10 −9 . Since Φ ⊙ ≪ 1 in most of local gravity experiments, the parameter m is constrained to be much smaller than the order of unity in the region where the Ricci scalar R B is much larger than the present cosmological one (R 0 ∼ R c ).
Cosmologically the condition (31) implies that viable models need to be very close to the ΛCDM model in radiation and matter dominated epochs (R ≫ R 0 ). However, deviations from the ΛCDM model are allowed around the present acceleration epoch (R ∼ R 0 ). Thus in viable models the parameter m is negligibly small during the radiation and matter eras, but continues to grow by the present epoch.
For the theories of the type (29) the corresponding Brans-Dicke field φ, the potential V (φ) and the mass squared M 2 are given by In these regimes the field is stuck around φ = 1 because of the presence of a ρ-dependent term. When R decreases to the order of R c , the field begins to evolve along the potential V (φ) with a lighter mass M which is not very much different from R c . Thus in the Brans-Dicke description, the departure from the point φ = 1 amounts to deviation from the ΛCDM model. The models (27) and (28) are constructed to satisfy the stability conditions where R 1 is a Ricci scalar at the late-time de-Sitter point. The first condition is required to avoid repulsive gravity, whereas the second ensures the absence of tachyons or ghosts. The second condition is also required for the consistency with LGC (as was shown above) as well as to ensure the stability of density perturbations [22,31] (as we shall see below). We also note that the requirements (31) and (35) are perfectly consistent with the condition 0 < m(R) ≪ 1 derived in Ref. [11] for the existence of the matter era.
To summarise, the conditions (31) and (35), together with the existence of the de-Sitter point (16), are required for the viability of f (R) models in the metric formalism. The condition for the existence of the saddle matter era given in Eq. (15) is automatically satisfied under the requirements (31) and (35).

B. Palatini formalism
Let us next consider cosmological and local gravity constraints for f (R) theories in the Palatini formalism. We first discuss the conditions for the cosmological viability of f (R) models at the background level. In Ref. [23] it was shown that radiation (P r ), matter (P m ) and de-Sitter (P A ) points exist irrespective of the forms of f (R) provided that the function is well-behaved (i.e., it does not show discontinuous or divergent behaviour). Note that effective equations of state corresponding to points P r , P m and P A are given by w eff = 1/3, 0, −1, respectively. The de-Sitter point P A corresponds to F R − 2f = 0, i.e., and C(R) = 0. This de-Sitter solution exists on the same line as in the metric case. Since the eigenvalues of the Jacobian matrix for perturbations around the point P d are (λ 1 , λ 2 ) = (−3 − C(R), −4 − C(R)), the de-Sitter point on the line r = −2 is always a stable attractor. This situation is different from the metric case in which the stability of the de-Sitter point requires the additional condition 0 < m(r = −2) ≤ 1. The stability of radiation and matter points, on the other hand, depends upon the particular f (R) models chosen. The eigenvalues for perturbations are given by P r : (λ 1 , λ 2 ) = (4 + C(R), 1) and P m : (λ 1 , λ 2 ) = (3 + C(R), −1) [23]. Thus the models with C(R) > −3 give rise to the unstable node for P r and the saddle point for P m . Thus models satisfying the condition C(R) > −3 lead to a sequence of radiation, matter and de-Sitter epochs. As an example, let us consider the following model [17,18] where µ and n are constants. In this case one has C(R) = 3n in the regime R n+1 ≫ µ 2(n+1) , which means that a successful background trajectory is realised for n > −1. Note that a stable de-Sitter solution exists with R 1+n 1 = (2 + n)µ 2(n+1) and C(R) = 0. Obviously the constraints for a successful trajectory, at least at the background level, are not so severe compared to the metric formalism. One does not even require the condition m > 0 for the existence of a viable matter era.
Let us next discuss LGC in the Palatini case. Considering the Brans-Dicke picture of f (R) theories, the theories in the Palatini formalism correspond to ω BD = −3/2. Now since the usual description of the interaction range determined by the inverse of mass M can not be applied in this case, one needs to proceed in a different way [13]. From Eq. (6) the field φ = F (R) depends upon the value of the trace T , i.e., φ = φ(T ). We expand the field around the vacuum: ]. Note that we use the non-relativistic approximation under which the velocity v of matter is much smaller than the speed of light c. Carrying out a post-Newtonian expansion around the Minkowski vacuum (g µν = η µν + h µν ) in the solar system, the solutions for the second-order perturbation equations are given by [13] h (2) where V 0 = V (φ 0 ). Here the effective gravitational constant and the post-Newtonian parameter are where M ⊙ and M V are given by and ρ is the energy density of the sun. To ensure LGC, three conditions need to be satisfied [13]: • (iii) The contribution of the term log (φ/φ 0 ) is negligible.
The first condition comes from the experimental requirement γ ≈ 1. Since it is not easy to interpret this requirement directly, we shall elucidate this by considering a specific f (R) model later. Concerning condition (ii), setting T = 0 in Eq. (6) and using (11) When the deviation from the ΛCDM model is not large, the term f (R 0 )/f ,R0 (R 0 ) is of the order of the present cosmological Ricci scalar R 0 ∼ H 2 0 . Hence on the scales of the solar system this condition is well satisfied.
Regarding condition (iii), the presence of the term log (φ/φ 0 ) in Eqs. (39) and (40) leads to an additional acceleration of particles that should be small in order to be consistent with experiments. From the validity of classical Euler equations, the condition (iii) translates to [13] where ρ s is an energy density of the local structure. This means that the field φ(T ) should not have a strong dependence on T . Using the relations T = 2V − φV ,φ and φV ,φφ − V ,φ = f ,R /f ,RR − R, this condition translates to It is sometimes useful to rewrite this condition in terms of the variable m thus: When |m(R s )| ≪ 1, this is well satisfied since both f ,Rs and ρ s /R s are of the order of unity. Note that this constraint is not so restrictive compared to the metric formalism. This can be understood by recalling that in the Palatini case the field is non-dynamical without an interaction range. In the metric formalism one needs a large scalaron mass M to satisfy the thin-shell condition, which leads to a very small value of m(R s ) satisfying Eq. (31). We also note that in the Palatini case the condition f ,RR > 0 is not required in order to satisfy LGC. As a concrete example, let us apply the above constraints to the theories given by Eq. (38) with n > −1. In order to give rise to a late-time acceleration, µ needs to be of the order of the present Hubble radius H 0 . The field φ and the potential V (φ) defined in Eq. (11) are in this case given by Now in the de-Sitter case [vacuum (T = 0)] the solution R 0 satisfies which for the model (38) gives and φ 0 = 2(n + 1) n + 2 , V 0 = n + 1 (n + 2) n/(n+1) µ 2 .
In settings where local gravity experiments are carried out, the parameter ǫ ≡ µ 2 /R s ∼ ρ 0 /ρ s is much smaller than unity. For example if we take the mean density ρ s = 10 −11 g/cm 3 and use the typical values µ 2 ∼ H 2 0 ∼ ρ 0 = 10 −29 g/cm 3 and R s ∼ ρ s , then ǫ is of the order of 10 −18 . When n > 0, then in the limit ǫ → 0, we have φ → 1 and V (φ) → 0. Thus in the expression of M ⊙ given in Eq. (42) the term V 0 /φ 0 dominates over the term V /φ, thus giving Now since µ 2 ∼ ρ 0 ≪ ρ s , then the condition (i) is well satisfied. When −1 < n < 0, then as ǫ → 0 one has φ → 1 and the potential V becomes of the order V ∼ where M ⊙ is the same as that in Eq. (52). The ratio of the integrands in the expressions for M V and M ⊙ can be estimated to be (ρ 0 /ρ s ) n+1 ≪ 1, which means that the condition |M V | ≪ M ⊙ is again satisfied. The parameter m in this case is given by Now since ǫ is much smaller than 1 we obtain |m(R s )| ≪ 1. Hence theories of type (38) with n > −1 can satisfy local gravity constraints. The above discussion shows that it is easier to satisfy the local gravity constraints in the Palatini case than in the metric case. In the latter case we also require the condition f ,RR > 0 to ensure that the scalaron mass squared M 2 is positive. Moreover the requirement of the heavy mass M gives very small values for m(R s ), which imposes the condition that viable f (R) models need to be very close to the ΛCDM model during matter and radiation epochs. We also note that even though the condition |m(R s )| ≪ 1 is also required in the Palatini case, the absolute values of |m(R s )| do not need to be vanishingly small. In fact even models (38) with n > 0 can satisfy the correct Newtonian limit, while they are excluded in the metric formalism because f ,RR is negative in those cases. Thus in the Palatini formalism the models of the type f (R) = R − g(R) can be consistent with local gravity tests provided that the contribution of the term g(R) is not significant relative to the linear term.
In subsequent sections we discuss the evolution of density perturbations for f (R) theories in both metric and Palatini formalisms. We shall carry out a detailed analysis for a number of f (R) models that can satisfy both the cosmological and local gravity constraints and use the evolution of density perturbations to place constraints on the model parameters as well as their deviation from the ΛCDM model.

III. MATTER PERTURBATIONS AND GAUGE ISSUES
In this section we present the equations for matter perturbations that are applicable to f (R) theories in both metric and Palatini formalisms. As our background spacetime we shall consider a flat Friedmann-Lemaitre-Robertson-Walker (FLRW). The perturbed FLRW metric which includes linear scalar metric perturbations α, b, ϕ and E can be written in the form [44] where a(t) is a scale factor. We shall consider a pressure-less matter source. Now since the flow is irrotational for scalar perturbations, we can introduce a velocity potential V in terms of which the components of the energy momentum tensor of the pressure-less matter can be decomposed as where v m is related to the velocity potential V through [45]: Note that the definition of b in Ref. [45] has an opposite sign to that used here. We recall that in both metric and Palatini formalisms the matter energy density ρ m satisfies the standard continuity equatioṅ where H =ȧ/a is the Hubble parameter and a dot represents a derivative with respect to t. The matter perturbation can then be shown to satisfy the following equations of motion in the Fourier space [46,47]: where k is a comoving wavenumber and Defining the following variables where v is a covariant velocity perturbation [48]. Eqs. (59) and (60) can be written as Now choosing a comoving hypersurface, the density perturbation can be expressed in a gauge-invariant way as [45]: We shall define the density contrast on comoving orthogonal hypersurfaces as Now since the right hand side of Eq. (66) is gauge-invariant, δ m can be evaluated in any gauge. The evolution equation for δ m is given byδ where B = Hv − ϕ.
In the following we shall consider three different gauges to fix the gauge degree of freedom: • Comoving gauge: in which the 3-velocity and the scalar shift function vanish (i.e., v = 0). This implies that along with the 3-velocity the momentum vanishes as well. Thus the gauge-invariant δ m in this gauge becomes • Longitudinal gauge: in which the shift vector b and the anisotropic potential E both vanish, resulting in χ = 0. The gauge-invariant δ m in this gauge becomes • Uniform density gauge: in which we have constant density hypersurfaces, i.e. δρ m = 0. The gauge-invariant δ m in this gauge becomes In what follows we discuss the evolution of matter perturbations in the above three gauges, while noting that the physics does not depend upon the choice of gauges. Also while the above discussions hold in both metric and Palatini formalisms, the equation of matter perturbations is different in each case as we shall see in the following sections.

IV. DENSITY PERTURBATIONS IN THE METRIC FORMALISM
In the metric formalism the background equations are given by where R = 6(2H 2 +Ḣ). Note that we only take into account a non-relativistic matter. In Fourier space the scalar metric perturbations satisfy the following equations of motion, in the so-called gauge-ready form [46] − In the following we derive the perturbation equations in the above gauges, in both exact forms as well as using a sub-horizon approximation.

A. Comoving gauge (v = 0)
We first derive the equation of matter perturbations in the comoving gauge (v = 0). When v = 0 we have (63) and (64). Hence from Eq. (76) we find whereas from Eq. (77), the perturbation δF satisfies The evolution of the matter perturbations δ (v) m can then be obtained by solving Eqs. (78) and (79) numerically. For models that satisfy local gravity constraints the mass squared term defined in Eq. (17) can be well approximated by M 2 ≃ f,R 3f,RR , a term which appears on the left hand side of Eq. (79). We are mainly interested in the evolution of perturbations on sub-horizon scales, i.e., We also recall that for the models that satisfy LGC the mass of the scalar field squared M 2 is much larger than R ∼ H 2 ∼ |Ḣ|. Hence either k 2 /a 2 or M 2 is dominant in the parenthesis on the left hand side of Eq. (79). Let us first consider the case in which the time-derivative terms in δF are neglected, i.e., The condition (81) amounts to neglecting the term δF that leads to the oscillation of δF . This is the approximation used in scalar tensor models in Refs. [6,49,50]. Later we explore the validity of such an approximation paying particular attention to the conditions that should be satisfied. Under the conditions (80) and (81), Eq. (79) gives Using the approximation (81) in Eq. (78), we obtain where the "cosmological" effective gravitational constant is given by Note that we have restored the bare gravitational constant G.
Introducing a physical wavelength ℓ = a/k, the parameter ξ defined in Eq. (83) can be written as where in the last approximate equality we have used the approximate relation M 2 ≃ f,R 3f,RR . In the regimes ξ ≪ 1, i.e., (M ℓ) 2 ≫ 1, Eq. (85) gives G cosmo eff ≃ G/F . In this case m ≪ 1 for sub-horizon modes (k ≫ aH). Thus the deviation from the ΛCDM model is small, i.e., |Ḟ /HF | ≪ 1 in Eq. (84). Hence the evolution of matter perturbations is similar to the one in the standard GR case. We recall again that this General Relativistic behaviour can be realized even for ω BD = 0 because of the presence of a potential with a heavy scalar-field mass (M 2 ≫ k 2 /a 2 ).
In the regimes ξ ≫ 1, i.e., (M ℓ) 2 ≪ 1, Eq. (85) gives G cosmo eff ≃ 4G/3F . Thus in this case the evolution of matter perturbations is different from the one in the GR case because of the appearance of the 4/3 factor. If the mass of the Brans-Dicke scalar field is light (M 2 ≪ k 2 /a 2 ), the cosmological effective gravitational constant in Brans-Dicke theory is given by G cosmo eff ≃ G φ 4+2ωBD 3+2ωBD [50]. Thus in the regime ξ ≫ 1, the f (R) theories in the metric formalism behave as the Brans-Dicke theory (with ω BD = 0), with a light scalar-field mass (M 2 ≪ k 2 /a 2 ).

B. Longitudinal gauge (χ = 0)
We shall also derive the approximate equations in the longitudinal gauge (χ = 0) for sub-horizon modes satisfying Eq. (80). We also use the notation α = Φ and ϕ = −Ψ, which then gives the relation Ψ = Φ+δF/F from Eq. (75). In addition to Eq. (81), we impose the following conditions and If the deviation from the ΛCDM model is not significant, the condition (87) is well satisfied. The condition (88) is also satisfied for sub-horizon modes given in Eq. (80) provided that Φ, Ψ and B are of the same order. Under these approximations we obtain, from Eqs. (67), (73), (76) and (77), the following relations From Eq. (74) the term v is of the order of HΦ/ρ m provided that the deviation from the ΛCDM model is not significant. Using Eq. (90) we find that the ratio 3Hv/(δρ m /ρ m ) is of the order of (aH) 2 /k 2 , which is much smaller than unity for sub-horizon modes. This gives δ Compared to the comoving gauge the difference appears only in the friction term. Since viable f (R) models satisfy the condition |Ḟ /HF | ≪ 1, Eq. (84) reduces to Eq. (91). We have also checked that in uniform density gauge (δρ m = 0) the perturbation δ (δ) m satisfies the same approximate equation as Eq. (91). Before ending this Subsection, we shall introduce a number of parameters which can be useful below. A useful parameter is the effective gravitational potential which characterises the deviation of light rays. This is directly linked with the Integrated Sachs-Wolfe (ISW) effect in the CMB [22,31] and weak lensing of distant galaxies [50,51]. From Eq. (90) we can approximate this parameter by We introduce an anisotropic parameter which behaves as η → 1 for ξ ≫ 1 and η → 2ξ for ξ ≪ 1. We also define another variable m . Using the above expressions Σ can be approximated by Note that Σ is directly linked with Φ eff . The parameters (Σ, η) can be especially important in future survey of weak lensing [50,51].

C. The appearance of scalarons
Among the approximations we have used in the previous two subsections, the conditions (35) and (81) can be violated if an oscillating mode (scalaron) dominates over the matter induced mode discussed above. Let us clarify when the oscillating mode becomes important for viable f (R) models satisfying the conditions m ≪ 1 and |Ḟ /HF | ≪ 1. For the sub-horizon modes, Eq. (79) is approximately given by The solution of this equation is the sum of the matter induced mode δF ind and the oscillatory scalaron mode δF osc satisfying Under the condition {M 2 , k 2 /a 2 } ≫ H 2 this equation reduces to the form (a 3/2 δF osc )¨+ ω 2 (a 3/2 δF osc ) ≃ 0, where ω = k 2 /a 2 + M 2 . In the adiabatic regime characterised by |ω/ω 2 | ≪ 1 we obtain the following WKB solution where c is a constant. Hence the solution of the perturbation δR is expressed by For viable f (R) models, the scale factor a and the background Ricci scalar R (0) evolve as a ∝ t 2/3 and R (0) ≃ 4/(3t 2 ) during the matter era. Then the amplitude of δR osc relative to R (0) has a time-dependence Let us consider the models m(r) = C(−r − 1) p (p > 0) for which the mass M evolves as M ∝ t −(p+1) during the matter-dominated epoch. When ξ ≪ 1 and ξ ≫ 1 we have |δR osc |/R (0) ∝ t −(3p+1)/2 and |δR osc |/R (0) ∝ t −2(p+1/3) , respectively. Hence the amplitude of the oscillating mode decreases faster than the background Ricci scalar. Thus if the scalaron is over-produced in the early Universe such that |δR| > R (0) , the stability condition (35) can be violated. This property persists in the radiation-dominated epoch as well [28,29]. Thus in order to ensure the viability of the f (R) theories of gravity in metric formalism, we need to esnure that |δR| is smaller than R (0) at the beginning of the radiation era. This can be achieved by choosing the constant c in Eq. (99) to be sufficiently small which amounts to a fine tuning for these theories. We note that this fine tuning concerns the stability of these theories and is an additional constraint to those usually imposed on the parameters of these theories by observations. Under the condition that the scalaron mode δR osc is negligible relative to the mater-induced mode δR ind , one can derive the evolution for the matter perturbation δ m as well as the effective gravitational potential Φ eff . When ξ ≪ 1 the evolutions of δ m and Φ eff during the matter era are given by Note that the ratio of the matter induced mode relative to the background Ricci scalar evolves as For the models that satisfy cosmological and local gravity constraints, the Universe typically starts from the regime ξ ≪ 1 and evolved into the regime with ξ ≫ 1 during the matterdominated epoch [28,29]. When ξ ≫ 1, δ m and Φ eff evolve as For the models m(r) = C(−r − 1) p , we have the time-dependence |δR ind |/R (0) ∝ t −2p+( √ 33−5)/6 in the regime ξ ≫ 1. This decreases more slowly relative to the ratio |δR osc |/R (0) ∝ t −2(p+1/3) , so the scalaron mode tends to be unimportant with time.
In what follows we shall numerically solve the exact perturbation equations in order to check the validity of approximations used to reach Eqs. (84), (91) and (93). We choose initial conditions such that scalaron mode is suppressed relative to the matter induced mode, i.e. |δR i osc | < |δR i ind |. We refer the reader to Ref. [29] for a comprehensive and detailed study of the scalaron mode. This study also gives the conditions under which the scalaron mode dominates over the matter induced mode at the initial stages.

D. Numerical study of the validity of approximations
In order to study the dynamics of matter perturbations in the metric formalism we shall introduce the following dimensionless variables [11] In terms of these variables, the energy fraction Ω m of the pressureless matter and the effective equation of state w eff are given by The evolution equations for the background dynamics can then be expressed as [11] x where a prime denotes a derivative with respect to the number of e-folding N = log (a). For a later use we also introduce the variable x 4 ≡ aH, which satisfies The matter epoch corresponds to the critical point Since m needs to be much smaller than unity during the matter era, we have P M : (x 1 , x 2 , x 3 ) ≃ (0, −1/2, 1/2) for viable f (R) models. We shall consider a case in which the evolution proceeds from the matter point P M to the de-Sitter point given by: which can easily be shown to lie on the line r = −2 and is stable for 0 < m ≤ 1.

Comoving gauge
In the comoving gauge, the perturbation Eqs. (78) and (79) can be rewritten in terms of the above variables thus: where δF ≡ δF/F . The exact evolution of the matter perturbation can be obtained by solving these equations together with the background equations (106)-(109) for x 1 , x 2 , x 3 and x 4 . Meanwhile, the approximate equation (84) can be expressed in terms of these variables as where Let us consider the case in which the condition M 2 ≫ k 2 /a 2 (i.e., ξ ≪ 1) is satisfied. Since M needs to be large during the matter-dominated epoch to satisfy LGC, this condition holds in viable f (R) models at the beginning of the matter era for the modes relevant to large scale structure [28,29]. Then the term 2x 3 /m dominates over the term k 2 /x 2 4 in Eq. (113), which gives δF ∼ mδ  The perturbations can enter the regime M 2 ≪ k 2 /a 2 (i.e., ξ ≫ 1) before reaching the present epoch, depending on the mode k and on the evolution of M [28,29]. In models m(r) = (−r − 1) 3 this occurs for the modes k/a 0 H 0 > 3.5, where the subscript 0 represent present values. In the case k/a 0 H 0 = 300, the redshift at k/a = M corresponds to z k = 4.83. Since M 2 is always larger than H 2 in the past because of the requirement m ≪ 1, the modes are inside the Hubble radius (k 2 /a 2 H 2 > 1) after the perturbations enter the regime M 2 < k 2 /a 2 . Hence the approximation we used to reach Eq. (114) is valid in this regime. In the regime M 2 < k 2 /a 2 the term (k 2 /x 2 4 )δF in Eq. (113), balances the term ( m , which gives rise to an additional contribution on the right hand side of Eq. (112). This then leads to the approximate equation (114) with ξ ≫ 1, which has a growing-mode solution δ m ∝ t ( √ 33−1)/6 . In Fig. 1 we plot the evolution of δ (v) m for the model m(r) = (−r − 1) 3 with the mode k/a 0 H 0 = 10. Initial conditions are chosen so that the scalaron mode does not dominate over the matter-induced mode. In this case the transition from the regime M 2 > k 2 /a 2 to the region M 2 < k 2 /a 2 occurs at the redshift z k = 1.62. We find that the approximate equation (114) shows an excellent agreement with the results obtained by numerically solving the exact equations (112) and (113). The argument also holds for modes that are initially outside the Hubble radius. Thus the approximate equation (84) is reliable to estimate the growth of matter perturbations and the resulting matter power spectrum, provided that the scalaron does not dominate in the early Universe.

Longitudinal gauge
In the longitudinal gauge the combination of Eqs. (73)-(77) leads to the following perturbation equations δF ′′ + (x 3 + 2)δF ′ + 4 3 where we have used Ψ = Φ + δF . The effective potential defined in Eq. (92) is given by In order to understand the evolution of perturbations at the initial stages of the matter era, let us consider the regime ξ ≪ 1 without assuming the sub-horizon condition k/(aH) ≫ 1. We have in mind viable f (R) models with vanishingly small values of m deep inside the matter epoch. Equation (117) then becomes Note that under sub-horizon approximation we have δF ≃ −2mk 2 Φ/3(aH) 2 , which agrees with Eq. (90). Using Eq. (121) we find that the right hand side of Eq. (116) can be neglected relative to the left hand side, thus giving the solution Φ = constant (together with a decaying mode proportional to t −5/3 ). From Eqs. (92) and (118) we obtain Φ eff ≃ Φ and Equation (122) agrees with the expression (93) obtained under the sub-horizon approximation (k/a ≫ H).
Since Φ eff is a constant, the matter perturbation can be seen from Eq. (122) to evolve as δ (χ) m ∝ a. This is consistent with the approximate equation (91), i.e., which has the growing mode solution δ m ∝ a in the regime ξ ≪ 1. We may ask why the above method reproduces the result derived under the sub-horizon approximation, without employing the approximation k/a ≫ H. In the regime ξ ≪ 1 the perturbation δF is suppressed relative to Φ as given in Eq. (121). This allows us to neglect the right hand side of Eq. (116), giving a constant Φ. This mimics the situation in General Relativity where δF = 0 and Φ = constant together with Eq. (122), resulting in δ (χ) m ∝ a. Moreover, from Eq. (119), the quantity B = Hv + Ψ is well approximated by B ≃ 5Φ/3 = constant. Hence the right hand side of Eq. (67) can be neglected even without assuming the sub-horizon approximation. Thus using the relation (122) we can obtain Eq. (91) in the regime ξ ≪ 1 without assuming k/a ≫ H. The above approximation corresponds to the limit of large M (M 2 ≫ k 2 /a 2 ), which gives rise to the evolution of perturbations close to the case of General Relativity. In General Relativity (δF = 0 andḞ = 0), one has the exact equation (122) from Eqs. (73) and (74). Thus the perturbations in the large M case (ξ ≪ 1) mimic those in General Relativity, apart from the fact that the scalaron is present in the former but not in the latter.
When ξ ≫ 1 one has k 2 /a 2 ≫ M 2 ≫ H 2 , which means that the sub-horizon type approximation we used in the subsection B holds well in this regime. This situation is similar to the case of the comoving gauge. For the modes that start from the regime M 2 ≫ k 2 /a 2 and enter the regime M 2 ≪ k 2 /a 2 before the end of the matter era, the evolution of perturbations changes from the standard general relativistic form (102) to the non-standard form (103).
In Fig. 2 we plot the evolution of δ (χ) m and Φ eff in the model m(r) = (−r − 1) 3 for the mode k = a 0 H 0 that lies outside the Hubble radius at the start of integration (z = 28.9). Together with numerically integrating Eqs. (116)-(118), we also solve the approximate equation (124) with Φ eff derived by (93). From Fig. 2 we find that the approximate equations agree well with the exact numerical results, even if the mode is initially slightly outside the Hubble radius. We caution, however, that for large-scale modes far outside the Hubble radius the scalaron can be important. In fact we have numerically checked that the oscillating mode appears for such super-Hubble modes unless the coefficient of the scalaron in Eq. (100) is fine-tuned to be small. In Fig. 2 the growth of the gravitational potential is not seen in the region z < z k , since the transition redshift is small (z k = 0.36). It can, however, be observed if we consider modes on smaller scales.  (123). Initial conditions for the background quantities were chosen to be the same as in Fig. 1. We also plot δ (χ) m and Φ eff obtained by solving the approximate equations (124) and (93). The approximation is valid even when the mode is initially outside the Hubble radius.
In summary, for viable f (R) models that satisfy the cosmological and local gravity constraints, the approximate Eqs. (124) and (93) are trustable even for the modes outside the Hubble radius, as long as the scalaron is suppressed relative to the matter-induced mode.

E. Constraints on the model m(r) = C(−r − 1) p
Finally we consider the current and future constraints on models m(r) = C(−r − 1) p with 0 < C ≤ 1. At the background level, compatibility with the SNIa observations could result in the divergence of the equation of state of dark energy [12,29]. Interestingly the redshift at which such a divergence may occur could be of order unity. However the current SNIa observations are not yet sufficiently accurate to rule out such cases. Some constraints on the model parameters can be obtained from the present equation of state of dark energy, but even models with p = 1.5 and C = 1 are allowed [29]. Thus the background does not provide strong constraints on the model parameters. However this situation can change in the future when higher-redshift data will become available from the observations of SNIa and Gamma Ray bursts.
There are a number of additional observational constraints on the growth rate of matter perturbations. At the redshift z ∼ 3, McDonald et al. [52] obtained the constraint δ ′ m /δ m = 1.46 ± 0.49 from the measurement of the matter power spectrum from the Lyman-α forests. Taking into account the more recent data reported by Viel and Haehnelt [53] in the redshift range 2 < z < 4, the maximum value of the growth rate allowed by the current observations is given by [54] The current data still have large error bars and some data even allow the parameter range −1 < s < 2 [53].
However it is expected that in future observations the growth rate will be constrained more severely. From Eqs. (102) and (103) we have s = 1 for M 2 ≫ k 2 /a 2 and s = ( √ 33 − 1)/4 = 1.186 for M 2 ≪ k 2 /a 2 . In Fig. 3 we plot the evolution of the growth rate for models m(r) = (−r − 1) 3 for a number of different values of k. The increase of s from unity implies that the perturbations enter the regime M 2 ≪ k 2 /a 2 . For smaller scale modes this transition occurs earlier, which leads to the larger maximum value of s. The growth rate begins to decrease once the Universe enters the late-time accelerated epoch. As estimated analytically, the growth rate is bounded by s < 1.186. Hence the current observational constraint (125) is still too weak to place constraints on m(r) = C(−r − 1) p models.  Fig. 2. The transition redshift z k defined as the redshift where k/a = M , becomes larger for smaller scales. After the matter perturbation enters the region z < z k the growth rate begins to increase toward the value s = 1.186, but it starts to decrease once the Universe enters the stage of accelerated expansion.
However, these models exhibit peculiar features in the matter power spectrum. This is a consequence of the fact that there is a transition redshift z k at which the growth rate begins to change from s = 1 to s = 1.186. For the modes relevant to galaxy clusters (k/a 0 H 0 = O(10 2 )), this transition typically occurs during the matter-dominated epoch (see Fig. 3). Since the time t k at z = z k depends upon the modes k (t k ∝ k −3/(3p+1) ), this leads to the change in the slope of the matter power spectrum. The difference between the slopes of the matter power spectrum determined from galaxy surveys and the CMB spectrum, on the scales k/a 0 H 0 = O(10 2 ), is given by [28,29] ∆n ≃ √ 33 − 5 3p + 1 .
This analytic result agrees well with numerical results except for models with p ≫ 1 [29]. Observationally no significant differences have so far been found between the two power spectra. If we take the bound ∆n < 0.05 we obtain the constraint p ≥ 5. To place further constraints on models, a likelihood analysis is required which employs the data from both the galaxy power spectrum and the CMB. Such observational constraints, including the data from the Lyman-α forest and Gamma Ray bursts, are under consideration [55]. Numerically, we find that models m(r) = (−r − 1) 5 have constraints on the present value of the deviation parameter given by m(z = 0) < ∼ 10 −1 . Thus even though m is constrained to be very small during the matter era, a notable deviation from the ΛCDM model can occur around the present epoch.
Finally the ISW effect in the CMB power spectrum is important on large scale modes with k/a 0 H 0 = O(1). As can be seen from Fig. 2, even models with p = 3 and C = 1 do not give rise to a significant amplification of the gravitational potential. The models with p ≥ 2 are consistent with the low multipoles in the CMB data [29]. Thus this effect does not generally provide stronger additional constraints.

V. DENSITY PERTURBATIONS IN THE PALATINI FORMALISM
In this section we discuss the evolution of density perturbations in the Palatini formalism and the resulting observational consequences. The background equations are given by where the Ricci scalar R satisfies the following relations The matter perturbations satisfy Eqs. (63) and (64) as in the case of the metric formalism. Other perturbation equations are [47] − Note that we corrected several typos found in Ref. [47]. Given the non-dynamical nature of Eq. (135), it is clear that the scalaron mode does not exist in the Palatini case. This is associated with the fact that the Palatini formalism corresponds to generalised Brans-Dicke theory (Eq. (10)) with ω BD = −3/2. The perturbation δF is directly determined by the matter perturbation δρ m , as where m is defined in Eq. (14). Below we shall consider these perturbations in various gauges.

A. Comoving gauge
In the comoving gauge (v = 0) one has α = 0 and κ =δ Combining this with Eq. (136) we obtain an exact second-order equation for δ (v) m , which can be numerically solved (See Appendix A for details). On the other hand, since we are mostly interested in the evolution of modes on sub-horizon scales, it makes sense to consider the approximate equations similar to those considered in the metric case.
Using a sub-horizon type approximation, such that only the terms containing k 2 /a 2 and δρ m are considered on the right hand side of Eq. (137), together with Eq. (136), we obtain the following approximate perturbation equationδ where ξ is defined in Eq. (83). Alternatively we may study the case in which the deviation from the ΛCDM model is small, i.e., as required from LGC (46). The derivative terms such as |m ′ | and |m ′′ | are also assumed to be much smaller than unity. Using the fact that from Eqs. (128) and (136) (138) is valid even for the modes outside the Hubble radius. This situation is similar to the case in the metric formalism. In fact we have confirmed this property by numerically solving the exact equation (137) and comparing it with the solutions of the approximate equation (138).
One can estimate the order of the term RδF on the r.h.s. of Eq. (134) by using Eq. (136), i.e., RδF = mδρ m /(1 − m). This gives rise to the contribution of the order of (ρ m /2F )mδ (v) m in the third term of Eq. (138), which is negligible under the condition (139). As long as we neglect this contribution, we should approximate ξ/(1 − m) ≃ ξ in the third term of Eq. (138). In the following, we implicitly assume this when we write the term (1 − m) in the denominator.
In the limit ξ = k 2 a 2 R m ≪ 1, Eq. (138) agrees with Eq. (84) in the metric formalism. However a significant difference appears in the regime ξ ≫ 1. In that case there is a strong amplification of the matter perturbation in the Palatini case because of the growth of the term ξ in Eq. (138). We shall estimate this growth rate for a number of concrete models in Subsection C below.

B. Longitudinal gauge
We next consider the Longitudinal gauge (χ = 0), and as in the metric case we use the notation α = Φ and ϕ = −Ψ. Under the sub-horizon type approximation used in the comoving case above, the evolution equation reduces to Eq. (89) obtained in the metric case. Using Eqs. (131) and (133) together with Eq. (136), these approximations also give Hence the matter perturbations satisfy the following approximate equation The effective gravitational potential Φ eff defined in Eq. (92) satisfies which is the same as in the metric case. Similarly the parameters η and Σ defined in Eqs. (94) and (95) in this case become We note that while the expression for η is different from that in the metric case, Σ remains the same. The above approximate equations (141) and (142) are valid under the conditions (139) and ξ ≪ 1 even without the sub-horizon approximation. In fact the argument is similar to the metric case in which Eqs. (91) and (93) reduce to the corresponding GR equations for ξ ≪ 1.
We also note that in regimes ξ ≫ 1 the perturbation modes are inside the Hubble radius, which shows that the sub-horizon approximation above is still valid. Thus, as long as the condition (139) is satisfied, we can safely use Eqs. (141) and (142) even for super-Hubble modes. Also since in the Palatini case the perturbation δR is sourced by the matter induced mode only, we do not need to worry about the dominance of the scalaron oscillations for super-Hubble modes, unlike in the metric case.
In Ref. [47] the equation for matter perturbations was derived in the uniform density gauge (δρ m = 0). This is an exact equation as in the comoving gauge (see also Ref. [56]). We shall show in Appendix B that similar to other gauges Eq. (139) is sufficient to reduce the exact equation to the approximate one. We shall also show that the approximate equation is the same as the one in the longitudinal gauge (141).
Since the evolution of matter perturbations hardly depends on the gauge chosen, we shall in what follows denote the matter perturbations simply by δ m .

C. Analytic estimate for the growth of perturbations
As we mentioned above the evolution of perturbations in the regime ξ ≪ 1 is similar to the standard GR case where δ m ∝ t 2/3 , s = δ ′ m /δ m = 1 and Φ ∝ constant. In this subsection we shall estimate the growth rate of perturbations after the system enters the regime ξ > 1. We shall consider models with |m| ≪ 1, for consistency with LGC (46).
During the matter era in which the Ricci scalar evolves as R ∝ t −2 the parameter ξ is given by ξ = ±ma/m k a k , where the subscript "k" denotes the values when the system crosses ξ = 1. Here we note that the plus sign corresponds to a positive m, whereas the negative sign to a negative m. As we showed in Sec. II, the latter case can be allowed unlike the metric case. Under the condition |m| ≪ 1, the matter perturbation (141) satisfies the following equation Lets us consider the case in which the evolution of the parameter m is given by where p is a constant. The values of p in several different f (R) models are given by With the above choice of m, Eq. (144) reduces to For the positive sign in Eq. (146), i.e., for m > 0, the solution of Eq. (146) can be written in terms of a linear combination of Bessel functions J ν and Y ν : where c 1 , c 2 are constants and For the negative sign in Eq. (146), i.e., for m < 0, the solution of Eq. (146) is given by where x and ν are as given in (148). In the following we shall discuss the positive and negative m cases in turn.

m > 0
As an example, we consider the constant m models (p = 0). In this case, the Bessel function J ν (ix) has a growing mode solution Then, in the regimes ξ ≫ 1, the evolution of the matter perturbations and its growth rate are given by where we have used √ 6e (N −N k )/2 ≫ (N − N k )/2. Thus the growth rate s of the matter perturbations increases very rapidly. Also from Eq. (142), in the regimes ξ ≫ 1, the effective gravitational potential grows (double) exponentially as which leads to a strong and observable ISW effect. Similarly, in models with p = 0, one can estimate the evolution of perturbations in the regime ξ ≫ 1: This shows that for models with p > 0 the growth rate increases faster than in the constant m models. When p < −1/3 the above instability can be avoided, but in that case m increases towards the past. Thus unless the present value of m is negligibly small, the condition |m| ≪ 1 required for LGC can be violated during the matter era. Thus these models are hardly distinguishable from the ΛCDM model in the present Universe. In what follows we shall concentrate on the positive p case. We note that for these positive m models the violent growth of δ m and Φ eff comes to an end after the Universe enters the accelerated stage, since ξ and Ω m begin to decrease.
If the Universe crosses the critical point |ξ| = 1 around the end of the matter era, it does not necessarily reach the regime |ξ| ≫ 1. In such cases one can not fully use the above approximate solutions. We shall show later that, in some cases, the Universe can enter the accelerated stage without oscillations in δ m up to the present epoch. The oscillations in δ m can be seen as we choose larger values of |m| and k. The frequency of oscillations tends to grow for larger values of p. The models that enter the regimes |ξ| ≫ 1 generally have difficulty in being consistent with observations, since they lead to largely negative values of s as given by Eq. (155).

Constraints on |m| from the requirement |ξ| < 1
The f (R) models can be consistent with observations if the Universe does not enter the regime |ξ| > 1 until the end of the matter-dominated epoch. One can estimate the ratio of the comoving Hubble radius (aH) −1 during the matter era to its present value thus: where c = 1 in the absence of the dark energy dominated epoch. The presence of a dark energy era leads to a change in the value of c. Numerically this factor is around c = 1.7-1.9. Using the relation R ≃ 3H 2 that holds during the matter era for |m| ≪ 1, we find that |ξ| crosses 1 at a critical redshift If z c is smaller than the order of unity, the Universe does not enter the regime |ξ| > 1 during the matter dominated epoch. This gives the following constraint: The matter power spectrum, in the linear regime, has been observed for the scales 0.02h Mpc −1 < ∼ k < ∼ 0.2h Mpc −1 . Non-linear effects can be important for smaller scales with k > 0.2h Mpc −1 . Taking the value k = 0.2h Mpc −1 ≃ 600a 0 H 0 , below which linear perturbation theory is valid, we obtain the constraint |m(z)| < ∼ 3 × 10 −6 during the matter era. Of course this is a rough estimate and the actual constraints on m(z) depend upon the particular models considered. For example, even if |ξ| crosses 1 during the matter era, the models can be consistent with observations provided that |ξ| does not grow rapidly after the crossing. Whether or not ξ reaches the regime |ξ| ≫ 1 depends on the particular models chosen. Hence to place constraints on m, we need a detailed analysis for each particular model. In the next subsection we shall provide a numerical investigation of a number of f (R) models presented above and place constraints on present values of m as well as model parameters.

D. Constraints on model parameters
In this subsection we shall employ the information provided by the growth of the matter density perturbations to place constraints on the parameters of the f (R) models presented in subsection C above. We do this by numerically solving the exact perturbation (165) together with the background equations (127) and (128). See Appendix C for equations written in a form convenient for numerical integration.
1. Constant m models: f (R) = αR 1+m − Λ Compared to other models considered here, the growth of |ξ| is rather mild in the constant m models, being of the form |ξ| ∝ a = e N . Thus, in order for |ξ| to grow from 0.1 to 10, one would require an increase in the number of e-foldings by 4.6.
We shall first consider the positive m case. In the left panel of Fig. 4 we plot the evolution of growth rate s = δ ′ m /δ m for the mode k = 600a 0 H 0 for several values of m. For m = 3 × 10 −5 we numerically obtain z c ∼ 11, denoted by a black dot in Fig. 4. This almost agrees with the analytical estimate (157) which gives z c ≈ 10. In the regimes ξ ≪ 1 the evolution of matter perturbations is given by δ m = δ ′ m ∝ a, which results in s ≃ 1. The growth rate s begins to move away from unity as ξ becomes of order 0.1, and then continues to grow before the Universe enters the stage of accelerated expansion. For this model we find s max ∼ 2.06 and ξ max ∼ 3.13, which shows that the model does not enter the regime ξ ≫ 1 where the evolution of perturbations is described by Eqs. (150) and (151).
For a model with m = 1.5×10 −5 , the critical redshift occurs at around z c ∼ 5 with s ∼ 1.4. The maximum value of the growth rate is s max ∼ 1.57, which corresponds to the marginal case satisfying the criterion (125). For a model with m = 2.0 × 10 −6 , the evolution of perturbations is not much different from the general relativistic case. To show the variation of the growth rate as a function of scale, we depict in the right panel of Fig. 4 the evolution of s for the model m = 1.5 × 10 −5 for three different values of k. As can be seen, the maximum value of the growth rate s decreases as k is decreased (i.e. the scales become larger). In particular, for the mode k = 100a 0 H 0 (corresponding to k = 0.33h Mpc −1 ), the evolution of perturbations exhibits no difference compared to the corresponding evolution in the general relativistic case. Hence the matter power spectrum is enhanced on small scales (k = 0.1h-0.2h Mpc −1 ), while the spectrum remains similar to the standard general relativistic case on larger scales (k = 0.02h-0.04h Mpc −1 ). This results in different spectral indices on different scales. Placing more precise constraints on m, would require performing a likelihood analysis using the data from the matter power spectrum. However, in order to obtain an order of magnitude for the maximum value of m, it is sufficient to use the criterion (125) for the mode k = 600a 0 H 0 . For the constant m models we find the constraint to be m < ∼ 10 −5 . When m is negative, the growth rate s decreases unlike the positive m case. In the left panel of Fig. 5 we plot the evolution of s for three different negative values of m for the mode k/a 0 H 0 = 600. As can be seen s tends to decrease more rapidly with increasing |m|. If m = −2.0 × 10 −5 the present value of s becomes very small (s < −1). As we see in the right panel of Fig. 5, when m = −2.0 × 10 −5 , there is a significant fall in the values of s between k/a 0 H 0 = 300 and k/a 0 H 0 = 600. This can lead to large differences in the spectral indices of the matter power spectrum for small and large scale modes. From the above argument |m| should be smaller than the order of 10 −5 , which has an upper bound similar to the positive m case.
When m = −2.0 × 10 −5 the Universe crosses the point |ξ| = 1 at the redshift z c ∼ 7.4, but the increase of |ξ| for z < z c is mild. Moreover the quantity |ξ| begins to decrease after the Universe enters the accelerated stage. Numerically we obtain the value ξ ∼ −0.77 at present (z = 0). Thus the system does not reach the regime |ξ| ≫ 1, and hence not a single period of oscillation occurs by the present epoch. However, for larger values of |m|, we have numerically checked that the oscillations of δ m indeed occur. We also recall that the growth of effective gravitational potential Φ eff leads to an ISW effect in the CMB spectrum on large scales (k/a 0 H 0 ∼a few). However, when |m| ∼ 10 −5 , Φ eff does not grow for these modes. As a result the ISW effect does not provide stronger constraints m than those provided by the matter power spectrum.
The observational constraints on the parameters of this model were studied in Ref. [47]. [Note that R c is not very much different from the present value of the cosmological Ricci scalar R 0 .] Here, we shall obtain constraints on the parameter m which for this model is given by and make a comparison between our results. The late-time de-Sitter point (R = R 1 ) is obtained from the constraint equation F R − 2f = 0, to give (R 1 /R c ) 1−β = λ(2 − β). Thus at this de-Sitter point the variable m satisfies m(R 1 ) = β/2 .
For β < 1, the parameter m in the regime R ≫ R c is given by which decreases towards the past. If β (< 1) is of the order of unity, the quantity m is too large to satisfy the requirement (158) for the mode k = 600a 0 H 0 during the matter era (recall that from Eq. (160) the present value of m is of the order of β). This is basically associated with the fact that, in the regimes R ≫ R c , the model gives a linear relation between m and r [m = C(−r − 1)]. Thus we need the condition |β| ≪ 1 in order to be compatible with the criterion (158).
To determine the changes in the behaviour of this model as a function of β, we considered three distinct values of β and calculated the corresponding growth rate s and the parameter m in each case. Our results are summarised in Fig. 6. The left hand panel shows the evolution of the growth rate s for λ = 1 and k = 600a 0 H 0 with the three different values of β. For β = 1.5 × 10 −4 the present value of the parameter m is around m(z 0 ) ∼ 6.7 × 10 −5 , which is close to the value of m at the de-Sitter point (m(R 1 ) = 7.5 × 10 −5 ). We also find that the parameter ξ crosses 1 at a redshift z c ∼ 3 with m(z c ) ∼ 1.2 × 10 −5 . Furthermore, we find that the growth rate s is larger for these models than in the case of constant m models. This is due to the fact that ξ in this case evolves faster, as ξ ∝ t 2 (4/3−β) . The maximum growth rate reached for β = 1.5 × 10 −4 corresponds to s max ∼ 1.88 with ξ ∼ 4. As expected, models with smaller values of β possess growth rates which are more compatible with observational constraints. Employing the criterion (125) for the mode k = 600a 0 H 0 , we find the constraint β < 8.2 × 10 −5 . This is slightly larger than the constraint β < 3.0 × 10 −5 obtained in [35] from the likelihood analysis of the SDSS data [57]. In the left panel of Fig. 6 we also consider this case in order to find the corresponding evolution of s. The maximum value of the growth rate in this case is found to be s max ∼ 1.095, which indicates that the constraint (125) is rather weak. Nevertheless, the criterion (125) is certainly sufficient in order to extract the order of magnitude of the bound on β.
In the right panel of Fig. 6 we plot the evolution of the parameter m for the case with λ = 1 and k = 600a 0 H 0 , for three values of β. As can be seen m increases from the past to the present. Using the criterion (125) we obtain the bound m(z = 0) < 3.5 × 10 −5 . If we adopt the severer criterion s < 1.1, the constraint becomes m(z = 0) < 1.3 × 10 −5 . Thus the deviation from the ΛCDM model is constrained to be small (m(z = 0) < ∼ 10 −5 ). We also studied the effects of changing the parameter λ in the action on the bounds on β. We considered two cases with λ = 10 and λ = 100. We found that these changes in λ have negligible effects on the constraints imposed on β and m(z = 0), compared to that obtained from the case with λ = 1. The reason for this lack of sensitivity is that a change in parameter λ is compensated by corresponding changes to the values of R i , a 0 and H 0 .
When β < 0 the parameter m is negative from Eq. (161). In the left panel of Fig. 7 we plot the evolution of s for three different values of β with k/a 0 H 0 = 600. We find that the present values of s become smaller than −1 for |β| > ∼ 1.2 × 10 −4 , in which case |m(z = 0)| is smaller than the order of 5.3 × 10 −5 (see the right panel of Fig. 7). Thus if we use the criterion s(z = 0) > ∼ −1 for the validity of the models, the upper bounds of |β| and |m(z = 0)| are similar to those in the positive β case.
Finally we consider the above model (where n > 0) recently discussed by Starobinksy [28]. The parameter m for this model is given by and the de-Sitter point at R = R 1 corresponds to Once we fix the value of λ, x 1 is known accordingly. In the regime R ≫ R c the parameter m behaves as Because of the presence of a larger power of (R c /R) in the expression for m in this case, m decreases more rapidly towards the past compared to the model f (R) = R − λR c (R/R c ) β discussed above. For the mode k = 600a 0 H 0 , the bound (158) implies that m has to be smaller than the order of 10 −6 -10 −5 by the end of the matter-dominated epoch, in order to ensure that the model does not enter the regime with ξ > 1.
In Fig. 8 we plot, for mode k = 600a 0 H 0 , the evolution of s and m for λ = 2.5 with three different values of n. When n = 3.07 the critical redshift is given by z c ∼ 1.05 with m ∼ 1 × 10 −5 . The rapid increase of s occurs in the regime ξ > 1, after which the growth rate reaches a maximum value s max ∼ 2. The present value of m is found to be m = 4.5 × 10 −4 , which is an order of magnitude larger than its corresponding value at ξ = 1. Using the criterion (125), we obtain the constraints n > 3.23 and m(z = 0) < 2.9 × 10 −4 for λ = 2.5. The present value of m in this model is one order of magnitude larger than the corresponding values in the constant m models as well as f (R) = R − λR c (R/R c ) β model. We also find that in contrast to the model f (R) = R−λR c (R/R c ) β the constraints on n for the Starobinsky model are sensitive to the values of the parameter λ. For larger values of λ the constraints on n is weaker. For example, for λ = 10 and λ = 50 we find the corresponding constraints on n imposed by (125) to be n > 1.74 (m(z = 0) ∼ 1.4 × 10 −4 ) and n > 1.09 (m(z = 0) ∼ 1.1 × 10 −4 ) respectively. This can be understood in the following way. When λ is increased, we obtain a larger ratio R 1 /R c from Eq. (163), which also leads to a larger ratio R/R c in the past. Then from Eq. (164) a smaller value of n is sufficient to realize the condition |m| ≪ 1. It can also be seen from the form of the action that the values of R c can also affect the constraints on n. We find that for small λ values, R c has a small effect on the constraint, whereas for large values of λ the affect of changing R c is negligible.
From Eq. (164) we find that m can be negative for −1/2 < n < 0 (and λ > 0) in the regime R ≫ R c . When n is close to 0, the models are close to the model f (R) = R − λR c (R/R c ) β discussed above. We find that s(z = 0) is larger than −1 for |n| < 9.3 × 10 −5 , in which case we have |m(z = 0)| < 4.5 × 10 −5 . When n is close to −1/2, Eq. (164) seems to suggest that the models should be close to the constant m models. However, care needs to be taken in this case since m changes sign from negative to positive at (R/R c ) 2 = 1/(2n + 1) in the deep matter dominated epoch. As a result, for n close to −1/2, we numerically find that the growth rate s shows a rapid growth for (R/R c ) 2 < 1/(2n + 1). Thus, in the limit n → −1/2, the models do not behave as constant m models and they are excluded observationally.
In summary, the present values of m are constrained to be m(z = 0) < ∼ 10 −4 from the bound (125) in both Starobinsky and Hu & Sawicki models.

VI. CONCLUSIONS
We have made a detailed study of the evolution of density perturbations in f (R) gravity theories in both metric and Palatini formalisms and employed them to study the viability of models in each case. To study the viability of concrete models we considered three sets of constraints, provided by the background cosmological evolution, local gravity experiments and the evolution of matter density perturbations respectively.
We began by considering the cosmological and local gravity constraints. For models satisfying these constraints, we proceeded to study the additional constraints provided by the evolution of density perturbations to further constrain the model parameters as well as their deviation from the ΛCDM model.
The f (R) theories in the metric formalism are equivalent to generalised Brans-Dicke theories with a scalar-field potential V (φ) and Brans-Dicke parameter ω BD = 0. The presence of the field potential, allows the construction of f (R) models that satisfy the local gravity constraints under the use of a chameleon mechanism. We find that for typical models of the forms (27) and (28) to satisfy the cosmological and local gravity constraints, the parameter m is required to be much smaller than unity during the radiation and matter eras but can grow to values of order of 0.1 in the accelerated epoch. Models in the metric formalism also suffer from an additional fine tuning due to the presence of scalaron oscillating modes (which is absent in the Palatini case). Finally, to be stable these theories require f ,RR to be positive.
On the other hand, the f (R) theories in the Palatini formalism correspond to generalised Brans-Dicke theories with a scalar-field potential V (φ) and Brans-Dicke parameter ω BD = −3/2. This makes these theories special in the sense that the oscillating scalar degree of freedom (scalaron) is absent in these theories and therefore the corresponding fine tuning to the metric case does not exist. Also unlike the case of the metric formalism, there is in this case no notion of field mass M that determines an interaction length mediated by a fifth force. Thus the LGC for these theories need to be analyzed separately in contrast to theories with ω BD = −3/2. The main condition required in this case in order to satisfy the LGC is that |m| is smaller than the order of unity. Moreover the requirement for the cosmological viability in the Palatini formalism is not severe compared to the metric case. Thus in contrast to the metric case, to satisfy the cosmological and local gravity constraints, we do not require vanishingly small values of m during radiation and matter dominated epochs and furthermore f ,RR can be negative in this case. As a result, even models of the type f (R) = R − µ 2(n+1) /R n with n > 0 can be allowed at the background level unlike the metric case.
We then studied the constraints provided by the evolution of density perturbations in each case. In the case of the metric formalism we derived the equations for matter perturbations under sub-horizon approximations in several different gauges. In regimes M 2 ≫ k 2 /a 2 (i.e., ξ ≪ 1), we found the approximate perturbation equations to be valid even without using sub-horizon approximations, provided that the scalaron mode is not dominant relative to the matter-induced mode. This is a consequence of the fact that, when M 2 ≫ k 2 /a 2 , the evolution of perturbations mimics that in General Relativity. After the Universe enters the regime M 2 ≪ k 2 /a 2 , the modes are inside the Hubble radius due to the fact that the condition M 2 ≫ R is required for the compatibility with LGC. Thus, for the models that satisfy LGC, as long as the scalarons do not dominate over the matter-induced mode, approximate perturbation equations are valid even for the modes that initially lie outside the Hubble radius. In the Palatini case the approximate equations are even more reliable because of the absence of scalarons.
In the metric formalism, most viable f (R) models take the form m(r) = C(−r − 1) p (p > 1) in the regimes where the Ricci scalar is larger than the order of the present cosmological value. In these models, the modes relevant to the observed matter power spectrum correspond to the regimes M 2 ≫ k 2 /a 2 with the growth rate s = δ ′ m /δ m = 1 at the beginning of the matter era. These models typically enter the regime M 2 ≪ k 2 /a 2 during the matter era in which the growth rate of matter perturbations is given by s = 1.186. If we use the present observational bound s < ∼ 1.5, we do not obtain strong constraints on these models. However, since the transition time at k/a = M depends upon the mode k, there is a difference in the spectral indices between the matter power spectrum and the CMB spectrum [see Eq. (126)]. If we take the bound ∆n < 0.05, the models with p ≥ 5 are allowed. The present value of the parameter m is constrained to be m(z = 0) < ∼ 10 −1 . Thus, while m needs to be negligibly small during the radiation and matter eras, one can have appreciable deviation from the ΛCDM model around the present epoch.
In the Palatini formalism the approximate matter perturbation equations are valid even for super-Hubble modes, for models satisfying LGC (|m| ≪ 1). If m is positive, there is a strong amplification of δ m in the regime ξ ≫ 1, whereas if m < 0 the matter perturbation exhibits a damped oscillation for |ξ| ≫ 1. When the quantity m evolves as m ∝ t 2p during the matter era, we have analytically estimated the growth rate s in both positive and negative m cases [see Eqs. (152) and (155), respectively]. From the requirement that the Universe does not enter the regimes |ξ| > 1 during the matter era, we obtain the constraint |m(z)| < ∼ (a 0 H 0 /k) 2 for z > O (1). While this is a good criterion to avoid non-standard evolution of matter perturbations, one needs to carry out a more detailed analysis to place constraints on the quantity m for each f (R) model. When m is positive, we have obtained the constraint m < ∼ 10 −5 by considering the modes k relevant to the matter power spectrum. We also studied the evolution of perturbations for the models f (R) = R− λR c (R/R c ) β and f (R) = R− λR c [1 − (1 + R 2 /R 2 c ) −n ]. For these models we found the constraints m(z = 0) < ∼ 10 −5 and m(z = 0) < ∼ 10 −4 , respectively, from the requirement s < ∼ 1.5. Thus, unlike the metric case, the deviation from the ΛCDM model at the present epoch is small even when m grows from the matter era to the accelerated epoch. This situation does not change for negative values of m.
In summary, for viable models in the metric formalism, the quantity m is constrained to be very much smaller than the order of unity during the matter era from LGC, but it can grow to the order of 0.1 around the present epoch. In the Palatini formalism, LGC and background cosmological constraints do not place strong bounds on m (only requiring |m| < ∼ 10 −1 ), but the density perturbations can provide stringent constraints: |m| < ∼ 10 −5 -10 −4 . Thus in the Palatini case the f (R) theories are hardly distinguishable from the ΛCDM model even at the present epoch. This follows from a peculiar evolution of the matter perturbations in the Palatini case, in the regime |ξ| > 1, that exhibits rapid growth (when m > 0) or damped oscillations (when m < 0).
While the constraints obtained here are sufficient to give the orders of magnitude of the allowed model parameters, it will be of interest to obtain more precise constraints by using recent and upcoming observational data including large scale structure, CMB, Supernova Ia, gamma ray bursts and weak lensing.
In this appendix we present the exact matter perturbation equation in the comoving gauge in the Palatini formalism. As in the metric case this equation needs to be solved simultaneously with the background equations (127)-(128). Unlike the metric case, however, it is not easy in this case to find dimensionless variables in terms of which both sets of equations close. As a result we proceed to integrate the equations directly. Using the background equations and ignoring the radiation, the perturbation equation (137) can be written as where the coefficients P 1 , P 2 , P 3 are given by where J is defined as All the terms in the coefficients P 1 , P 2 , P 3 can be expressed in terms of the scale factor a (or equivalently N ), which thus allows Eq. (165) to close and be readily integrated numerically. In the uniform density gauge (δρ m = 0) we have δF = 0 = δR from Eq. (136) andv (δ) = α, κ = 3Hv (δ) + k 2 a 2 v (δ) from Eqs. (63) and (64). Substituting these relations for Eq. (134) we obtain where we have used the relation Then the matter perturbation, δ Introducing a dimensionless quantity we get the differential equation for y [23]: where C(R) is defined in Eq. (36).
The following relations also hold Specifying the value of y, the initial Ricci scalar R is determined by Eq. (182). Solving Eq. (181), we obtain y, R, H and Ω m from Eqs. (182), (180) and (183). The effective equation of state of dark energy is given by As long as the deviation from the ΛCDM model is small (|m| ≪ 1), we have w eff ≃ −y 1 . The perturbation equations (141) and (142) are given by Although we solve exact perturbation equations, the above perturbation equations are found to be very accurate.