f(R) theories

Over the past decade, f(R) theories have been extensively studied as one of the simplest modifications to General Relativity. In this article we review various applications of f(R) theories to cosmology and gravity - such as inflation, dark energy, local gravity constraints, cosmological perturbations, and spherically symmetric solutions in weak and strong gravitational backgrounds. We present a number of ways to distinguish those theories from General Relativity observationally and experimentally. We also discuss the extension to other modified gravity theories such as Brans-Dicke theory and Gauss-Bonnet gravity, and address models that can satisfy both cosmological and local gravity constraints.


Introduction
General Relativity (GR) [225,226] is widely accepted as a fundamental theory to describe the geometric properties of spacetime. In a homogeneous and isotropic spacetime the Einstein field equations give rise to the Friedmann equations that describe the evolution of the universe. In fact, the standard big-bang cosmology based on radiation and matter dominated epochs can be well described within the framework of General Relativity.
However, the rapid development of observational cosmology which started from 1990s shows that the universe has undergone two phases of cosmic acceleration. The first one is called inflation [564,339,291,524], which is believed to have occurred prior to the radiation domination (see [402,391,71] for reviews). This phase is required not only to solve the flatness and horizon problems plagued in big-bang cosmology, but also to explain a nearly flat spectrum of temperature anisotropies observed in Cosmic Microwave Background (CMB) [541]. The second accelerating phase has started after the matter domination. The unknown component giving rise to this latetime cosmic acceleration is called dark energy [310] (see [517,141,480,485,171,32] for reviews). The existence of dark energy has been confirmed by a number of observations -such as supernovae Ia (SN Ia) [490,506,507], large-scale structure (LSS) [577,578], baryon acoustic oscillations (BAO) [227,487], and CMB [560,561,367].
These two phases of cosmic acceleration cannot be explained by the presence of standard matter whose equation of state = / satisfies the condition ≥ 0 (here and are the pressure and the energy density of matter, respectively). In fact, we further require some component of negative pressure, with < −1/3, to realize the acceleration of the universe. The cosmological constant Λ is the simplest candidate of dark energy, which corresponds to = −1. However, if the cosmological constant originates from a vacuum energy of particle physics, its energy scale is too large to be compatible with the dark energy density [614]. Hence we need to find some mechanism to obtain a small value of Λ consistent with observations. Since the accelerated expansion in the very early universe needs to end to connect to the radiation-dominated universe, the pure cosmological constant is not responsible for inflation. A scalar field with a slowly varying potential can be a candidate for inflation as well as for dark energy.
Although many scalar-field potentials for inflation have been constructed in the framework of string theory and supergravity, the CMB observations still do not show particular evidence to favor one of such models. This situation is also similar in the context of dark energy -there is a degeneracy as for the potential of the scalar field ("quintessence" [111,634,267,263,615,503,257,155]) due to the observational degeneracy to the dark energy equation of state around = −1. Moreover it is generally difficult to construct viable quintessence potentials motivated from particle physics because the field mass responsible for cosmic acceleration today is very small ( ≃ 10 −33 eV) [140,365]. While scalar-field models of inflation and dark energy correspond to a modification of the energy-momentum tensor in Einstein equations, there is another approach to explain the acceleration of the universe. This corresponds to the modified gravity in which the gravitational theory is modified compared to GR. The Lagrangian density for GR is given by ( ) = − 2Λ, where is the Ricci scalar and Λ is the cosmological constant (corresponding to the equation of state = −1). The presence of Λ gives rise to an exponential expansion of the universe, but we cannot use it for inflation because the inflationary period needs to connect to the radiation era. It is possible to use the cosmological constant for dark energy since the acceleration today does not need to end. However, if the cosmological constant originates from a vacuum energy of particle physics, its energy density would be enormously larger than the today's dark energy density. While the Λ-Cold Dark Matter (ΛCDM) model ( ( ) = − 2Λ) fits a number of observational data well [367,368], there is also a possibility for the time-varying equation of state of dark energy [10,11,450,451,630].
One of the simplest modifications to GR is the f (R) gravity in which the Lagrangian density is an arbitrary function of [77,512,102,106]. There are two formalisms in deriving field equations from the action in f (R) gravity. The first is the standard metric formalism in which the field equations are derived by the variation of the action with respect to the metric tensor . In this formalism the affine connection Γ depends on . Note that we will consider here and in the remaining sections only torsion-free theories. The second is the Palatini formalism [481] in which and Γ are treated as independent variables when we vary the action. These two approaches give rise to different field equations for a non-linear Lagrangian density in , while for the GR action they are identical with each other. In this article we mainly review the former approach unless otherwise stated. In Section 9 we discuss the Palatini formalism in detail.
The model with ( ) = + The spherically symmetric solutions mentioned above have been derived under the weak gravity backgrounds where the background metric is described by a Minkowski space-time. In strong gravitational backgrounds such as neutron stars and white dwarfs, we need to take into account the backreaction of gravitational potentials to the field equation. The structure of relativistic stars in f (R) gravity has been studied by a number of authors [349,350,594,43,600,466,42,167]. Originally the difficulty of obtaining relativistic stars was pointed out in [349] in connection to the singularity problem of f (R) dark energy models in the high-curvature regime [266]. For constant density stars, however, a thin-shell field profile has been analytically derived in [594] for chameleon models in the Einstein frame. The existence of relativistic stars in f (R) gravity has been also confirmed numerically for the stars with constant [43,600] and varying [42] densities. In this review we shall also discuss this issue.
It is possible to extend f (R) gravity to generalized BD theory with a field potential and an arbitrary BD parameter BD . If we make a conformal transformation to the Einstein frame [213,609,408,611,249,268], we can show that BD theory with a field potential corresponds to the coupled quintessence scenario [23] with a coupling between the field and non-relativistic matter. This coupling is related to the BD parameter via the relation 1/(2 2 ) = 3 + 2 BD [343,596]. One can recover GR by taking the limit → 0, i.e., BD → ∞. The f (R) gravity in the metric formalism corresponds to = −1/ √ 6 [28], i.e., BD = 0. For large coupling models with | | = (1) it is possible to design scalar-field potentials such that the chameleon mechanism works to reduce the effective matter coupling, while at the same time the field is sufficiently light to be responsible for the late-time cosmic acceleration. This generalized BD theory also leaves a number of interesting observational and experimental signatures [596].
In addition to the Ricci scalar , one can construct other scalar quantities such as and from the Ricci tensor and Riemann tensor [142]. For the Gauss-Bonnet (GB) curvature invariant defined by ≡ 2 − 4 + , it is known that one can avoid the appearance of spurious spin-2 ghosts [572,67,302] (see also [98,465,153,447,110,181,109]). In order to give rise to some contribution of the GB term to the Friedmann equation, we require that (i) the GB term couples to a scalar field , i.e., ( ) or (ii) the Lagrangian density is a function of , i.e., ( ). The GB coupling in the case (i) appears in lowenergy string effective action [275] and cosmological solutions in such a theory have been studied extensively (see [34,273,105,147,588,409,468] for the construction of nonsingular cosmological solutions and [463,360,361,593,523,452,453,381,25] for the application to dark energy). In the case (ii) it is possible to construct viable models that are consistent with both the background cosmological evolution and local gravity constraints [458,188,189] (see also [165,180,178,383,633,599]). However density perturbations in perfect fluids exhibit negative instabilities during both the radiation and the matter domination, irrespective of the form of ( ) [383,182]. This growth of perturbations gets stronger on smaller scales, which is difficult to be compatible with the observed galaxy spectrum unless the deviation from GR is very small. We shall review such theories as well as other modified gravity theories.
This review is organized as follows. In Section 2 we present the field equations of f (R) gravity in the metric formalism. In Section 3 we apply f (R) theories to the inflationary universe. Section 4 is devoted to the construction of cosmologically viable f (R) dark energy models. In Section 5 local gravity constraints on viable f (R) dark energy models will be discussed. In Section 6 we provide the equations of linear cosmological perturbations for general modified gravity theories including metric f (R) gravity as a special case. In Section 7 we study the spectra of scalar and tensor metric perturbations generated during inflation based on f (R) theories. In Section 8 we discuss the evolution of matter density perturbations in f (R) dark energy models and place constraints on model parameters from the observations of large-scale structure and CMB. Section 9 is devoted to the viability of the Palatini variational approach in f (R) gravity. In Section 10 we construct viable dark energy models based on BD theory with a potential as an extension of f (R) theories.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 In Section 11 the structure of relativistic stars in f (R) theories will be discussed in detail. In Section 12 we provide a brief review of Gauss-Bonnet gravity and resulting observational and experimental consequences. In Section 13 we discuss a number of other aspects of f (R) gravity and modified gravity. Section 14 is devoted to conclusions.
There are other review articles on f (R) gravity [556,555,618] and modified gravity [171,459,126,397,217]. Compared to those articles, we put more weights on observational and experimental aspects of f (R) theories. This is particularly useful to place constraints on inflation and dark energy models based on f (R) theories. The readers who are interested in the more detailed history of f (R) theories and fourth-order gravity may have a look at the review articles by Schmidt [531] and Sotiriou and Faraoni [556].
In this review we use units such that = = = 1, where is the speed of light, is reduced Planck's constant, and is Boltzmann's constant. We define 2 = 8 = 8 / 2 pl = 1/ 2 pl , where is the gravitational constant, pl = 1.22 × 10 19 GeV is the Planck mass with a reduced value pl = pl / √ 8 = 2.44 × 10 18 GeV. Throughout this review, we use a dot for the derivative with respect to cosmic time and " , " for the partial derivative with respect to the variable , e.g.,

Field Equations in the Metric Formalism
We start with the 4-dimensional action in f (R) gravity: where 2 = 8 , is the determinant of the metric , and ℒ is a matter Lagrangian 1 that depends on and matter fields Ψ . The Ricci scalar is defined by = , where the Ricci tensor is In the case of the torsion-less metric formalism, the connections Γ are the usual metric connections defined in terms of the metric tensor , as This follows from the metricity relation, ∇ = / − Γ − Γ = 0.

Equations of motion
The field equation can be derived by varying the action (2.1) with respect to : ). Einstein gravity, without the cosmological constant, corresponds to ( ) = and ( ) = 1, so that the term ( ) in Eq. (2.7) vanishes. In this case we have = − 2 and hence the Ricci scalar is directly determined by the matter (the trace ). In modified gravity the term ( ) does not vanish in Eq. (2.7), which means that there is a propagating scalar degree of freedom, ≡ ( ). The trace equation (2.7) determines the dynamics of the scalar field (dubbed "scalaron" [564]).
The field equation (2.4) can be written in the following form [568] = 2 (︁ ( ) + ( ) )︁ , (2.8) where ≡ − (1/2) and (2.9) Since ∇ = 0 and ∇ ( ) = 0, it follows that Hence the continuity equation holds, not only for Σ , but also for the effective energy-momentum tensor ( ) defined in Eq. (2.9). This is sometimes convenient when we study the dark energy equation of state [306,568] as well as the equilibrium description of thermodynamics for the horizon entropy [53]. There exists a de Sitter point that corresponds to a vacuum solution ( = 0) at which the Ricci scalar is constant. Since ( ) = 0 at this point, we obtain ( ) − 2 ( ) = 0 . (2.11) The model ( ) = 2 satisfies this condition, so that it gives rise to the exact de Sitter solution [564]. In the model ( ) = + 2 , because of the linear term in , the inflationary expansion ends when the term 2 becomes smaller than the linear term (as we will see in Section 3). This is followed by a reheating stage in which the oscillation of leads to the gravitational particle production. It is also possible to use the de Sitter point given by Eq. (2.11) for dark energy.
We consider the spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime with a time-dependent scale factor ( ) and a metric where is cosmic time. For this metric the Ricci scalar is given by = 6(2 2 +˙) , (2.13) where ≡˙/ is the Hubble parameter and a dot stands for a derivative with respect to . The present value of is given by where ℎ = 0.72 ± 0.08 describes the uncertainty of 0 [264]. The energy-momentum tensor of matter is given by ( ) = diag (− , , , ), where is the energy density and is the pressure. The field equations (2.4) in the flat FLRW background give 16) where the perfect fluid satisfies the continuity equatioṅ
Note that there are some works about the Einstein static universes in f (R) gravity [91,532]. Although Einstein static solutions exist for a wide variety of f (R) models in the presence of a barotropic perfect fluid, these solutions have been shown to be unstable against either homogeneous or inhomogeneous perturbations [532].

Equivalence with Brans-Dicke theory
The f (R) theory in the metric formalism can be cast in the form of Brans-Dicke (BD) theory [100] with a potential for the effective scalar-field degree of freedom (scalaron). Let us consider the following action with a new field , Meanwhile the action in BD theory [100] with a potential ( ) is given by where BD is the BD parameter and (∇ ) 2 ≡ . Comparing Eq. (2.21) with Eq. (2.23), it follows that f (R) theory in the metric formalism is equivalent to BD theory with the parameter BD = 0 [467,579,152] (in the unit 2 = 1). In Palatini f (R) theory where the metric and the connection Γ are treated as independent variables, the Ricci scalar is different from that in metric f (R) theory. As we will see in Sections 9.1 and 10.1, f (R) theory in the Palatini formalism is equivalent to BD theory with the parameter BD = −3/2.

Conformal transformation
The action (2.1) in f (R) gravity corresponds to a non-linear function in terms of . It is possible to derive an action in the Einstein frame under the conformal transformation [213,609,408,611,249,268,410]:˜= Ω 2 , (2.24) where Ω 2 is the conformal factor and a tilde represents quantities in the Einstein frame. The Ricci scalars and˜in the two frames have the following relation  (2.33) Hence the Lagrangian density of the field is given by ℒ = − 1 2˜− ( ) with the energy-momentum tensor (2.34) The conformal factor Ω 2 = = exp( √︀ 2/3 ) is field-dependent. From the matter action (2.32) the scalar field is directly coupled to matter in the Einstein frame. In order to see this more explicitly, we take the variation of the action (2.32) with respect to the field :

35)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Using Eq. (2.24) and the relations √ −˜= 2 √ − and˜= −1 , the energy-momentum tensor of matter is transformed as˜( (2.37) The energy-momentum tensor of perfect fluids in the Einstein frame is given bỹ The derivative of the Lagrangian density ℒ = ℒ ( ) = ℒ ( −1 ( )˜) with respect to is (2.39) The strength of the coupling between the field and matter can be quantified by the following quantity 40) which is constant in f (R) gravity [28]. It then follows that This shows that the field is directly coupled to matter apart from radiation (˜= 0). Let us consider the flat FLRW spacetime with the metric (2.12) in the Jordan frame. The metric in the Einstein frame is given by Equations (2.48) and (2.49) show that the field and matter interacts with each other, while the total energy density˜=˜+˜and the pressure˜=˜+˜satisfy the continuity equation d˜/d˜+ 3˜(˜+˜) = 0. More generally, Eqs. (2.48) and (2.49) can be expressed in terms of the energy-momentum tensors defined in Eqs. (2.34) and (2.37): which correspond to the same equations in coupled quintessence studied in [23] (see also [22]).
In the absence of a field potential ( ) (i.e., massless field) the field mediates a long-range fifth force with a large coupling (| | ≃ 0.4), which contradicts with experimental tests in the solar system. In f (R) gravity a field potential with gravitational origin is present, which allows the possibility of compatibility with local gravity tests through the chameleon mechanism [344,343].
In f (R) gravity the field is coupled to non-relativistic matter (dark matter, baryons) with a universal coupling = −1/ √ 6. We consider the frame in which the baryons obey the standard continuity equation ∝ −3 , i.e., the Jordan frame, as the "physical" frame in which physical quantities are compared with observations and experiments. It is sometimes convenient to refer the Einstein frame in which a canonical scalar field is coupled to non-relativistic matter. In both frames we are treating the same physics, but using the different time and length scales gives rise to the apparent difference between the observables in two frames. Our attitude throughout the review is to discuss observables in the Jordan frame. When we transform to the Einstein frame for some convenience, we go back to the Jordan frame to discuss physical quantities.

Inflation in f (R) Theories
Most models of inflation in the early universe are based on scalar fields appearing in superstring and supergravity theories. Meanwhile, the first inflation model proposed by Starobinsky [564] is related to the conformal anomaly in quantum gravity 3 . Unlike the models such as "old inflation" [339,291,524] this scenario is not plagued by the graceful exit problem -the period of cosmic acceleration is followed by the radiation-dominated epoch with a transient matter-dominated phase [565,606,426]. Moreover it predicts nearly scale-invariant spectra of gravitational waves and temperature anisotropies consistent with CMB observations [563,436,566,355,315]. In this section we review the dynamics of inflation and reheating. In Section 7 we will discuss the power spectra of scalar and tensor perturbations generated in f (R) inflation models.

Inflationary dynamics
We consider the models of the form which include the Starobinsky's model [564] as a specific case ( = 2). In the absence of the matter fluid ( = 0), Eq. (2.15) gives (3. 2) The cosmic acceleration can be realized in the regime = 1 + −1 ≫ 1. Under the approximation ≃ −1 , we divide Eq. During inflation the Hubble parameter evolves slowly so that one can use the approximation |˙/ 2 | ≪ 1 and |¨/(˙)| ≪ 1. Then Eq. (3.3) reduces tȯ . (3.4) Integrating this equation for 1 > 0, we obtain the solution The cosmic acceleration occurs for 1 < 1, i.e., > (1+ √ 3)/2. When = 2 one has 1 = 0, so that is constant in the regime ≫ 1. The models with > 2 lead to super inflation characterized by˙> 0 and ∝ | 0 − | −1/| 1 | ( 0 is a constant). Hence the standard inflation with decreasing occurs for (1 + √ 3)/2 < < 2. In the following let us focus on the Starobinsky's model given by ( ) = + 2 /(6 2 ) , (3.6) where the constant has a dimension of mass. The presence of the linear term in eventually causes inflation to end. Without neglecting this linear term, the combination of Eqs. (2.15) and (2.16) gives¨−˙2 During inflation the first two terms in Eq. (3.7) can be neglected relative to others, which giveṡ ≃ − 2 /6. We then obtain the solution where and are the Hubble parameter and the scale factor at the onset of inflation ( = ), respectively. This inflationary solution is a transient attractor of the dynamical system [407]. The accelerated expansion continues as long as the slow-roll parameter is smaller than the order of unity, i.e., 2 2 . One can also check that the approximate relation 3˙+ 2 ≃ 0 holds in Eq. (3.8) by using ≃ 12 2 . The end of inflation (at time = ) is characterized by the condition ≃ 1, i.e., ≃ / √ 6. From Eq. (3.11) this corresponds to the epoch at which the Ricci scalar decreases to ≃ 2 . As we will see later, the WMAP normalization of the CMB temperature anisotropies constrains the mass scale to be ≃ 10 13 GeV. Note that the phase space analysis for the model (3.6) was carried out in [407,24,131].
We define the number of e-foldings from = to = : Since inflation ends at ≃ + 6 / 2 , it follows that where we used Eq. (3.12) in the last approximate equality. In order to solve horizon and flatness problems of the big bang cosmology we require that 70 [391], i.e., 1 ( ) 7×10 −3 . The CMB temperature anisotropies correspond to the perturbations whose wavelengths crossed the Hubble radius around = 55 -60 before the end of inflation.

Dynamics in the Einstein frame
Let us consider inflationary dynamics in the Einstein frame for the model (3.6) in the absence of matter fluids (ℒ = 0). The action in the Einstein frame corresponds to (2.32) with a field defined by Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Using this relation, the field potential (2.33) reads [408,61,63] (3.16) In Figure 1 we illustrate the potential (3.16) as a function of . In the regime ≫ 1 the potential is nearly constant ( ( ) ≃ 3 2 /(4 2 )), which leads to slow-roll inflation. The potential in the regime ≪ 1 is given by ( ) ≃ (1/2) 2 2 , so that the field oscillates around = 0 with a Hubble damping. The second derivative of with respect to is which changes from negative to positive at = 1 ≡ √︀ 3/2(ln 2)/ ≃ 0.169 pl . Since ≃ 4 2 / 2 during inflation, the transformation (2.44) gives a relation between the cosmic time˜in the Einstein frame and that in the Jordan frame:
The field equations for the action (2.32) are given by Using the slow-roll approximations (d /d˜) 2 ≪ ( ) and |d 2 /d˜2| ≪ |˜d /d˜| during inflation, one has 3˜2 ≃ 2 ( ) and 3˜(d /d˜) + , ≃ 0. We define the slow-roll parameters For the potential (3.16) it follows that where, in the expression of˜2, we have dropped the terms of the order of 1/˜2. The results (3.27) will be used to estimate the spectra of density perturbations in Section 7.

Reheating after inflation
We discuss the dynamics of reheating and the resulting particle production in the Jordan frame for the model (3.6). The inflationary period is followed by a reheating phase in which the second derivative¨can no longer be neglected in Eq. (3.8). Introducing^= 3/2 , we havë Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Since 2 ≫ { 2 , |˙|} during reheating, the solution to Eq. (3.28) is given by that of the harmonic oscillator with a frequency . Hence the Ricci scalar exhibits a damped oscillation around = 0: Let us estimate the evolution of the Hubble parameter and the scale factor during reheating in more detail. If we neglect the r.h.s. of Eq. (3.7), we get the solution ( ) = const × cos 2 ( /2). Setting ( ) = ( ) cos 2 ( /2) to derive the solution of Eq. (3.7), we obtain [426] , (3.30) where os is the time at the onset of reheating. The constant is determined by matching Eq. (3.30) with the slow-roll inflationary solution˙= − 2 /6 at = os . Then we get = 3/ and Taking the time average of oscillations in the regime ( − os ) ≫ 1, it follows that ⟨ ⟩ ≃ (2/3)( − os ) −1 . This corresponds to the cosmic evolution during the matter-dominated epoch, i.e., ⟨ ⟩ ∝ ( − os ) 2/3 . The gravitational effect of coherent oscillations of scalarons with mass is similar to that of a pressureless perfect fluid. During reheating the Ricci scalar is approximately given by ≃ 6˙, i.e.
In the regime ( − os ) ≫ 1 this behaves as In order to study particle production during reheating, we consider a scalar field with mass . We also introduce a nonminimal coupling (1/2) 2 between the field and the Ricci scalar [88]. Then the action is given by where ( ) = + 2 /(6 2 ). Taking the variation of this action with respect to gives We decompose the quantum field in terms of the Heisenberg representation: where^and^ † are annihilation and creation operators, respectively. The field can be quantized in curved spacetime by generalizing the basic formalism of quantum field theory in the flat spacetime. See the book [88] for the detail of quantum field theory in curved spacetime. Then each Fourier mode ( ) obeys the following equation of motion + 3˙+ where the conformal coupling correspond to = 1/6. This result states that, even though = 0 (that is, the field is minimally coupled to gravity), still gives a contribution to the effective mass of . In the following we first review the reheating scenario based on a minimally coupled massless field ( = 0 and = 0). This corresponds to the gravitational particle production in the perturbative regime [565,606,426]. We then study the case in which the nonminimal coupling | | is larger than the order of 1. In this case the non-adiabatic particle production preheating [584,353,538,354] can occur via parametric resonance. In this case there is no explicit coupling among the fields and . Hence the particles are produced only gravitationally. In fact, Eq. (3.38) reduces to where = 2 /6. Since is of the order of ( ) 2 , one has 2 ≫ for the mode deep inside the Hubble radius. Initially we choose the field in the vacuum state with the positive-frequency solution [88]: The presence of the time-dependent term ( ) leads to the creation of the particle . We can write the solution of Eq. (3.39) iteratively, as [626] After the universe enters the radiation-dominated epoch, the term becomes small so that the flat-space solution is recovered. The choice of decomposition of into^and^ † is not unique. In curved spacetime it is possible to choose another decomposition in term of new ladder operatorŝ and^ † , which can be written in terms of^and^ † , such as^=^+ *^ † − . Provided that * ̸ = 0, even though^|0⟩ = 0, we have^|0⟩ ̸ = 0. Hence the vacuum in one basis is not the vacuum in the new basis, and according to the new basis, the particles are created. The Bogoliubov coefficient describing the particle production is The typical wavenumber in the -coordinate is given by , whereas in the -coordinate it is / . Then the energy density per unit comoving volume in the -coordinate is [426] where in the last equality we have used the fact that the term approaches 0 in the early and late times.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 During the oscillating phase of the Ricci scalar the time-dependence of is given by = ( ) sin( ∫︀ 0 d¯), where ( ) = ( ) 1/2 and = ( is a constant). When we evaluate the term d /d in Eq. (3.42), the time-dependence of ( ) can be neglected. Differentiating Eq. (3.42) in terms of and taking the limit ∫︀ 0 d¯≫ 1, it follows that . In [426] it was found that the transition from the oscillating phase to the radiation-dominated epoch occurs slower compared to the estimation given above. Since the epoch of the transient matter-dominated era is about one order of magnitude longer than the analytic estimation [426], we take the value ≃ os + 400 2 pl /( * 3 ) to estimate the reheating temperature . Since the particle energy density ( ) is converted to the radiation energy density = * 2 4 /30, the reheating temperature can be estimated as 4 In [426] the reheating temperature is estimated by taking the maximum value of reached around the ten oscillations of . Meanwhile we estimate at the epoch where becomes a dominant contribution to the total energy density (as in [364]).
As we will see in Section 7, the WMAP normalization of the CMB temperature anisotropies determines the mass scale to be ≃ 3 × 10 −6 pl . Taking the value * = 100, we have 5 × 10 9 GeV. For > the universe enters the radiation-dominated epoch characterized by ∝ 1/2 , = 0, and ∝ −2 .

Case: | | 1
If | | is larger than the order of unity, one can expect the explosive particle production called preheating prior to the perturbative regime discussed above. Originally the dynamics of such gravitational preheating was studied in [70,592] for a massive chaotic inflation model in Einstein gravity. Later this was extended to the f (R) model (3.6) [591]. Introducing a new field = 3/2 , Eq. (3.37) reads As long as | | is larger than the order of unity, the last two terms in the bracket of Eq. (3.51) can be neglected relative to . Since the Ricci scalar is given by Eq. (3.33) in the regime ( − os ) ≫ 1, it follows that¨+ The oscillating term gives rise to parametric amplification of the particle . In order to see this we introduce the variable defined by ( − os ) = 2 ± /2, where the plus and minus signs correspond to the cases > 0 and < 0 respectively. Then Eq. (3.52) reduces to the Mathieu equation . (3.54) The strength of parametric resonance depends on the parameters and . This can be described by a stability-instability chart of the Mathieu equation [419,353,591]. In the Minkowski spacetime the parameters and are constant. If and are in an instability band, then the perturbation grows exponentially with a growth index , i.e., ∝ . In the regime ≪ 1 the resonance occurs only in narrow bands around = ℓ 2 , where ℓ = 1, 2, ..., with the maximum growth index = /2 [353]. Meanwhile, for large (≫ 1), a broad resonance can occur for a wide range of parameter space and momentum modes [354].
In the expanding cosmological background both and vary in time. Initially the field is in the broad resonance regime ( ≫ 1) for | | ≫ 1, but it gradually enters the narrow resonance regime ( 1). Since the field passes many instability and stability bands, the growth index stochastically changes with the cosmic expansion. The non-adiabaticity of the change of the frequency 2 = 2 / 2 + 2 − 4 sin{ ( − os )}/( − os ) can be estimated by the quantity where the non-adiabatic regime corresponds to na 1. For small and we have na ≫ 1 around ( − os ) = , where are positive integers. This corresponds to the time at which the Ricci scalar vanishes. Hence, each time crosses 0 during its oscillation, the non-adiabatic particle production occurs most efficiently. The presence of the mass term tends to suppress Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 the non-adiabaticity parameter na , but still it is possible to satisfy the condition na 1 around = 0. For the model (3.6) it was shown in [591] that massless particles are resonantly amplified for | | 3. Massive particles with of the order of can be created for | | 10. Note that in the preheating scenario based on the model ( , ) = (1/2) 2 2 + (1/2) 2 2 2 the parameter decreases more rapidly ( ∝ 1/ 2 ) than that in the model (3.6) [354]. Hence, in our geometric preheating scenario, we do not require very large initial values of [such as > (10 3 )] to lead to the efficient parametric resonance.
While the above discussion is based on the linear analysis, non-linear effects (such as the mode-mode coupling of perturbations) can be important at the late stage of preheating (see, e.g., [354,342]). Also the energy density of created particles affects the background cosmological dynamics, which works as a backreaction to the Ricci scalar. The process of the subsequent perturbative reheating stage can be affected by the explosive particle production during preheating. It will be of interest to take into account all these effects and study how the thermalization is reached at the end of reheating. This certainly requires the detailed numerical investigation of lattice simulations, as developed in [255,254].
At the end of this section we should mention a number of interesting works about gravitational baryogenesis based on the interaction (1/ 2 * ) ∫︀ d 4 √ − between the baryon number current and the Ricci scalar ( * is the cut-off scale characterizing the effective theory) [179,376,514]. This interaction can give rise to an equilibrium baryon asymmetry which is observationally acceptable, even for the gravitational Lagrangian ( ) = with close to 1. It will be of interest to extend the analysis to more general f (R) gravity models.

Dark Energy in f (R) Theories
In this section we apply f (R) theories to dark energy. Our interest is to construct viable f (R) models that can realize the sequence of radiation, matter, and accelerated epochs. In this section we do not attempt to find unified models of inflation and dark energy based on f (R) theories.
Originally the model ( ) = − / ( > 0, > 0) was proposed to explain the late-time cosmic acceleration [113,120,114,143] (see also [456,559,17,223,212,16,137,62] for related works). However, this model suffers from a number of problems such as matter instability [215,244], the instability of cosmological perturbations [146,74,544,526,251], the absence of the matter era [28,29,239], and the inability to satisfy local gravity constraints [469,470,245,233,154,448,134]. The main reason why this model does not work is that the quantity , ≡ 2 / 2 is negative. As we will see later, the violation of the condition , > 0 gives rise to the negative mass squared 2 for the scalaron field. Hence we require that , > 0 to avoid a tachyonic instability. The condition , ≡ / > 0 is also required to avoid the appearance of ghosts (see Section 7.4). Thus viable f (R) dark energy models need to satisfy [568] , > 0 , where 0 is the Ricci scalar today.
In the following we shall derive other conditions for the cosmological viability of f (R) models. This is based on the analysis of [26]. For the matter Lagrangian ℒ in Eq. (2.1) we take into account non-relativistic matter and radiation, whose energy densities and satisfẏ (4.60)

Dynamical equations
We introduce the following variables together with the density parameters It is straightforward to derive the following equations  [160,382,488,252,31,198,280,72,41,159,235,1,279,483,321,432].
The effective equation of state of the system is defined by which is equivalent to eff = −(2 3 − 1)/3. In the absence of radiation ( 4 = 0) the fixed points for the above dynamical system are )︂ , The points 5 and 6 are on the line ( ) = − − 1 in the ( , ) plane. The matter-dominated epoch (Ω ≃ 1 and eff ≃ 0) can be realized only by the point 5 for close to 0. In the ( , ) plane this point exists around ( , ) = (−1, 0). Either the point 1 or 6 can be responsible for the late-time cosmic acceleration. The former is a de Sitter point ( eff = −1) with = −2, in which case the condition (2.11) is satisfied. The point 6 can give rise to the accelerated expansion In order to analyze the stability of the above fixed points it is sufficient to consider only timedependent linear perturbations ( ) ( = 1, 2, 3) around them (see [170,171] for the detail of such analysis). For the point 5 the eigenvalues for the 3 × 3 Jacobian matrix of perturbations are where 5 ≡ ( 5 ) and ′ 5 ≡ d d ( 5 ) with 5 ≈ −1. In the limit that | 5 | ≪ 1 the latter two eigenvalues reduce to −3/4 ± √︀ −1/ 5 . For the models with 5 < 0, the solutions cannot remain for a long time around the point 5 because of the divergent behavior of the eigenvalues Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 as 5 → −0. The model ( ) = − / ( > 0, > 0) falls into this category. On the other hand, if 0 < 5 < 0.327, the latter two eigenvalues in Eq. (4.77) are complex with negative real parts. Then, provided that ′ 5 > −1, the point 5 corresponds to a saddle point with a damped oscillation. Hence the solutions can stay around this point for some time and finally leave for the late-time acceleration. Then the condition for the existence of the saddle matter era is The first condition implies that viable f (R) models need to be close to the ΛCDM model during the matter domination. This is also required for consistency with local gravity constraints, as we will see in Section 5.
The eigenvalues for the Jacobian matrix of perturbations about the point 1 are . This shows that the condition for the stability of the de Sitter point 1 is [440,243,250,26] The trajectories that start from the saddle matter point 5 satisfying the condition (4.78) and then approach the stable de Sitter point 1 satisfying the condition (4.80) are, in general, cosmologically viable.
One can also show that 6 is stable and accelerated for (a) ′ 6 < −1, ( Since both 5 and 6 are on the line = − − 1, only the trajectories from ′ 5 > −1 to ′ 6 < −1 are allowed (see Figure 2). This means that only the case (a) is viable as a stable and accelerated fixed point 6 . In this case the effective equation of state satisfies the condition eff > −1.
From the above discussion the following two classes of models are cosmologically viable.

Viable f (R) dark energy models
We present a number of viable f (R) models in the ( , ) plane. First we note that the ΛCDM model corresponds to = 0, in which case the trajectory is the straight line (i) in Figure 2. The trajectory (ii) in Figure 2 represents the model ( ) = ( − Λ) [31], which corresponds to the straight line ( ) = [(1 − )/ ] + − 1 in the ( , ) plane. The existence of a saddle matter epoch demands the condition ≥ 1 and ≃ 1. The trajectory (iii) represents the model [26,382] which corresponds to the curve = (1 + )/ . The trajectory (iv) represents the model ( ) = − ( + 1)( 2 + + ), in which case the late-time accelerated attractor is the point 6 with ( In [26] it was shown that needs to be close to 0 during the radiation domination as well as the matter domination. with > 0, 0 < < 1, and (iv) ( ) = − ( + 1)( 2 + + ). From [31].
as it does not oscillate around = 0. The model ( ) = − / ( > 0, > 0) is not viable because the condition , > 0 is violated. As we will see in Section 5, the local gravity constraints provide tight bounds on the deviation parameter in the region of high density ( ≫ 0 ), e.g., ( ) 10 −15 for = 10 5 0 [134,596]. In order to realize a large deviation from the ΛCDM model such as ( ) > (0.1) today ( = 0 ) we require that the variable changes rapidly from the past to the present. The f (R) model given in Eq. (4.81), for example, does not allow such a rapid variation, because evolves as ≃ (− −1) in the region ≫ 0 . Instead, if the deviation parameter has the dependence it is possible to lead to the rapid decrease of as we go back to the past. The models that behave as Eq. (4.82) in the regime ≫ 0 are The models (A) and (B) have been proposed by Hu and Sawicki [306] and Starobinsky [568], respectively. Note that roughly corresponds to the order of 0 for = (1). This means that = 2 + 1 for ≫ 0 . In the next section we will show that both the models (A) and (B) are consistent with local gravity constraints for 1. In the model (A) the following relation holds at the de Sitter point:
Similarly the model (B) satisfies [568] (1 + 2 ) +2 ≥ 1 + ( + 2) 2 + ( + 1)(2 + 1) 4 , (4.87) . Another model that leads to an even faster evolution of is given by [587] (C) ( A similar model was proposed by Appleby and Battye [35]. In the region ≫ the model (4.89) behaves as ( ) ≃ − [1 − exp(−2 / )], which may be regarded as a special case of (4.82) in the limit that ≫ 1 5 . The Ricci scalar at the de Sitter point is determined by , as The models (A), (B) and (C) are close to the ΛCDM model for ≫ , but the deviation from it appears when decreases to the order of . This leaves a number of observational signatures such as the phantom-like equation of state of dark energy and the modified evolution of matter density perturbations. In the following we discuss the dark energy equation of state in f (R) models. In Section 8 we study the evolution of density perturbations and resulting observational consequences in detail.

Equation of state of dark energy
In order to confront viable f (R) models with SN Ia observations, we rewrite Eqs. (4.59) and (4.60) as follows: where is some constant and  Defining DE and DE in the above way, we find that these satisfy the usual continuity equatioṅ where the last approximate equality is valid in the regime where the radiation density is negligible relative to the matter density . The viable f (R) models approach the ΛCDM model in the past, i.e., → 1 as → ∞. In order to reproduce the standard matter era (3 2 ≃ 2 ) for ≫ 1, we can choose = 1 in Eqs. (4.92) and (4.93). Another possible choice is = 0 , where 0 is the present value of . This choice may be suitable if the deviation of 0 from 1 is small (as in scalar-tensor theory with a nearly massless scalar field [583,93]). In both cases the equation of state DE can be smaller than −1 before reaching the de Sitter attractor [306,31,587,435], while the effective equation of state eff is larger than −1. This comes from the fact that the denominator in Eq. (4.97) becomes smaller than 1 in the presence of the matter fluid. Thus f (R) gravity models give rise to the phantom equation of state of dark energy without violating any stability conditions of the system. See [210,417,136,13] for observational constraints on the models (4.83) and (4.84) by using the background expansion history of the universe. Note that as long as the late-time attractor is the de Sitter point the cosmological constant boundary crossing of eff reported in [52,50] does not typically occur, apart from small oscillations of eff around the de Sitter point.
There are some works that try to reconstruct the forms of f (R) by using some desired form for the evolution of the scale factor ( ) or the observational data of SN Ia [117,130,442,191,621,252]. We need to caution that the procedure of reconstruction does not in general guarantee the stability of solutions. In scalar-tensor dark energy models, for example, it is known that a singular behavior sometimes arises at low-redshifts in such a procedure [234,271]. In addition to the fact that the reconstruction method does not uniquely determine the forms of f (R), the observational data of the background expansion history alone is not yet sufficient to reconstruct f (R) models in high precision.
Finally we mention a number of works [115,118,119,265,319,515,542,90] about the use of metric f (R) gravity as dark matter instead of dark energy. In most of past works the power-law f (R) model = has been used to obtain spherically symmetric solutions for galaxy clustering. In [118] it was shown that the theoretical rotation curves of spiral galaxies show good agreement with observational data for = 1.7, while for broader samples the best-fit value of the power was found to be = 2.2 [265]. However, these values are not compatible with the bound | − 1| < 7.2 × 10 −19 derived in [62,160] from a number of other observational constraints. Hence, it is unlikely that f (R) gravity works as the main source for dark matter.

Local Gravity Constraints
In this section we discuss the compatibility of f (R) models with local gravity constraints (see [469,470,245,233,154,448,251] for early works, and [31,306,134] for experimental constraints on viable f (R) dark energy models, and [101,210,330,332,471,628,149,625,329,45,511,277,534,133,445,309,89] for other related works). In an environment of high density such as Earth or Sun, the Ricci scalar is much larger than the background cosmological value 0 . If the outside of a spherically symmetric body is a vacuum, the metric can be described by a Schwarzschild exterior solution with = 0. In the presence of non-relativistic matter with an energy density , this gives rise to a contribution to the Ricci scalar of the order 2 .
If we consider local perturbations on a background characterized by the curvature 0 , the validity of the linear approximation demands the condition ≪ 0 . We first derive the solutions of linear perturbations under the approximation that the background metric (0) is described by the Minkowski metric . In the case of Earth and Sun the perturbation is much larger than 0 , which means that the linear theory is no longer valid. In such a non-linear regime the effect of the chameleon mechanism [344,343] becomes important, so that f (R) models can be consistent with local gravity tests.

Linear expansions of perturbations in the spherically symmetric background
First we decompose the quantities , ( ), and into the background part and the perturbed part: = 0 + , = 0 (1 + ), and = (0) + about the approximate Minkowski background ( (0) ≈ ). In other words, although we consider close to a mean-field value 0 , the metric is still very close to the Minkowski case. The linear expansion of Eq. (2.7) in a time-independent background gives [470,250,154,448] The variable is defined in Eq. (4.67). Since 0 < ( 0 ) < 1 for viable f (R) models, it follows that 2 > 0 (recall that 0 > 0).
We consider a spherically symmetric body with mass , constant density (= − ), radius , and vanishing density outside the body. Since is a function of the distance from the center of the body, Eq. (5.1) reduces to the following form inside the body ( < ): whereas the r.h.s. vanishes outside the body ( > ). The solution of the perturbation for positive 2 is given by

5)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where ( = 1, 2, 3, 4) are integration constants. The requirement that ( ) > → 0 as → ∞ gives 4 = 0. The regularity condition at = 0 requires that 2 = − 1 . We match two solutions (5.4) and (5.5) at = by demanding the regular behavior of ( ) and ′ ( ). Since ∝ , this implies that is also continuous. If the mass satisfies the condition ≪ 1, we obtain the following solutions As we have seen in Section 2.3, the action (2.1) in f (R) gravity can be transformed to the Einstein frame action by a transformation of the metric. The Einstein frame action is given by a linear action in˜, where˜is a Ricci scalar in the new frame. The first-order solution for the perturbation ℎ of the metric˜= 0 ( + ℎ ) follows from the first-order linearized Einstein equations in the Einstein frame. This leads to the solutions ℎ 00 = 2 /( 0 ) and ℎ = 2 /( 0 ) . Including the perturbation to the quantity , the actual metric is given by [448] =˜≃ Using the solution (5.7) outside the body, the (00) and ( ) components of the metric are where ( ) eff and are the effective gravitational coupling and the post-Newtonian parameter, respectively, defined by For the f (R) models whose deviation from the ΛCDM model is small ( ≪ 1), we have 2 ≃ 0 /[3 ( 0 )] and ≃ 8 . This gives the following estimate where Φ = /( 0 ) = 4 2 /(3 0 ) is the gravitational potential at the surface of the body. The approximation ≪ 1 used to derive Eqs. (5.6) and (5.7) corresponds to the condition The validity of the linear expansion requires that ≪ 0 , which translates into ≪ ( 0 ). Since ≃ 2 /(3 0 ) = 2Φ /3 at = , one has ≪ ( 0 ) ≪ 1 under the condition (5.12). Hence the linear analysis given above is valid for ( 0 ) ≫ Φ .
For the distance close to the post Newtonian parameter in Eq. (5.10) is given by ≃ 1/2 (i.e., because ≪ 1). The tightest experimental bound on is given by [616,83,617]:

14)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 which comes from the time-delay effect of the Cassini tracking for Sun. This means that f (R) gravity models with the light scalaron mass ( ≪ 1) do not satisfy local gravity constraints [469,470,245,233,154,448,330,332]. The mean density of Earth or Sun is of the order of ≃ 1 -10 g/cm 3 , which is much larger than the present cosmological density (0) ≃ 10 −29 g/cm 3 . In such an environment the condition ≪ 0 is violated and the field mass becomes large such that ≫ 1. The effect of the chameleon mechanism [344,343] becomes important in this nonlinear regime ( ≫ 0 ) [251,306,134,101]. In Section 5.2 we will show that the f (R) models can be consistent with local gravity constraints provided that the chameleon mechanism is at work.

Chameleon mechanism in f (R) gravity
Let us discuss the chameleon mechanism [344,343] in metric f (R) gravity. Unlike the linear expansion approach given in Section 5.1, this corresponds to a non-linear effect arising from a large departure of the Ricci scalar from its background value 0 . The mass of an effective scalar field degree of freedom depends on the density of its environment. If the matter density is sufficiently high, the field acquires a heavy mass about the potential minimum. Meanwhile the field has a lighter mass in a low-density cosmological environment relevant to dark energy so that it can propagate freely. As long as the spherically symmetric body has a thin-shell around its surface, the effective coupling between the field and matter becomes much smaller than the bare coupling | |. In the following we shall review the chameleon mechanism for general couplings and then proceed to constrain f (R) dark energy models from local gravity tests.

Field profile of the chameleon field
The action (2.1) in f (R) gravity can be transformed to the Einstein frame action (2.32) with the coupling = −1/ √ 6 between the scalaron field = √︀ 3/(2 2 ) ln and non-relativistic matter. Let us consider a spherically symmetric body with radius˜in the Einstein frame. We approximate that the background geometry is described by the Minkowski space-time. Varying the action (2.32) with respect to the field , we obtain where˜is a distance from the center of symmetry that is related to the distance in the Jordan frame via˜= √ = − . The effective potential eff is defined by where * is a conserved quantity in the Einstein frame [343]. Recall that the field potential ( ) is given in Eq. (2.33). The energy density˜in the Einstein frame is related with the energy density in the Jordan frame via the relation˜= / 2 = 4 . Since the conformal transformation gives rise to a coupling between matter and the field,˜is not a conserved quantity. Instead the quantity * = 3 = −˜c orresponds to a conserved quantity, which satisfies˜3 * = 3 . Note that Eq. (5.15) is consistent with Eq. (2.42).
In the following we assume that a spherically symmetric body has a constant density * = inside the body (˜<˜) and that the energy density outside the body (˜>˜) is * = (≪ ). The mass of the body and the gravitational potential Φ at the radius˜are given by = (4 /3)˜3 and Φ = /˜, respectively. The effective potential has minima at the field values and : The former corresponds to the region of high density with a heavy mass squared 2 ≡ eff, ( ), whereas the latter to a lower density region with a lighter mass squared 2 ≡ eff, ( ). In the case of Sun, for example, the field value is determined by the homogeneous dark matter/baryon density in our galaxy, i.e., ≃ 10 −24 g/cm 3 .  . The field value, at which the inverted effective potential has a maximum, is different depending on the density * , see Eq. (5.22). In the upper panel "de Sitter" corresponds to the minimum of the potential, whereas "singular" means that the curvature diverges at = 0.
When > 0 the effective potential has a minimum for the models with , < 0, which occurs, e.g., for the inverse power-law potential ( ) = 4+ − . The f (R) gravity corresponds to a negative coupling ( = −1/ √ 6), in which case the effective potential has a minimum for , > 0. As an example, let us consider the shape of the effective potential for the models (4.83) and (4.84). In the region ≫ both models behave as Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 For this functional form it follows that The r.h.s. of Eq. (5.20) is smaller than 1, so that < 0. The limit → ∞ corresponds to → −0.
In the limit → −0 one has → /(2 2 ) and , → ∞. This property can be seen in the upper panel of Figure 3, which shows the potential ( ) for the model (4.84) with parameters = 1 and = 2. Because of the existence of the coupling term − / √ 6 * , the effective potential eff ( ) has a minimum at Since ∼ 2 * ≫ in the region of high density, the condition | | ≪ 1 is in fact justified (for and of the order of unity). The field mass about the minimum of the effective potential is given by . (5.23) This shows that, in the regime ∼ 2 * ≫ , is much larger than the present Hubble parameter 0 (∼ √ ). Cosmologically the field evolves along the instantaneous minima characterized by Eq. (5.22) and then it approaches a de Sitter point which appears as a minimum of the potential in the upper panel of Figure 3.
In order to solve the "dynamics" of the field in Eq. (5.15), we need to consider the inverted effective potential (− eff ). See the lower panel of Figure 3 for illustration [which corresponds to the model (4.84)]. We impose the following boundary conditions: The boundary condition (5.25) can be also understood as lim˜→ ∞ d /d˜= 0. The field is at rest at˜= 0 and starts to roll down the potential when the matter-coupling term in Eq. (5.15) becomes important at a radius˜1. If the field value at˜= 0 is close to , the field stays around in the region 0 <˜<˜1. The body has a thin-shell if˜1 is close to the radiusõ f the body.
In the region 0 <˜<˜1 one can approximate the r.h.s. of Eq.
In fact, this satisfies the boundary condition (5.24).
In the region˜1 <˜<˜the field | (˜)| evolves toward larger values with the increase of˜. In the lower panel of Figure 3 the field stays around the potential maximum for 0 <˜<˜1, but in the regime˜1 <˜<˜it moves toward the left (largely negative region). Since | , | ≪ | | in this regime we have that d eff /d ≃ in Eq. (5.15), where we used the condition ≪ 1. Hence we obtain the following solution where and are constants. Since the field acquires a sufficient kinetic energy in the region˜1 <˜<˜, the field climbs up the potential hill toward the largely negative region outside the body (˜>˜). The shape of the effective potential changes relative to that inside the body because the density drops from to . The kinetic energy of the field dominates over the potential energy, which means that the term d eff /d in Eq.  If the mass outside the body is small to satisfy the condition˜≪ 1 and ≫ , we can neglect the contribution of the -dependent terms in Eqs. (5.29) -(5.32). Then the field profile is given by [575] , (˜>˜) . Originally a similar field profile was derived in [344,343] by assuming that the field is frozen at = in the region 0 <˜<˜1. The radius 1 is determined by the following condition where Φ = 2 /(8˜) = 2˜2 /6 is the gravitational potential at the surface of the body. Using this relation, the field profile (5.37) outside the body reduces to If the field value at˜= 0 is away from , the field rolls down the potential for˜> 0. This corresponds to taking the limit˜1 → 0 in Eq. (5.40), in which case the field profile outside the body is given by This shows that the effective coupling is of the order of and hence for | | = (1) local gravity constraints are not satisfied.

Thin-shell solutions
Let us consider the case in which˜1 is close to˜, i.e.
Δ˜≡˜−˜1 ≪˜. (5.42) This corresponds to the thin-shell regime in which the field is stuck inside the star except around its surface. If the field is sufficiently massive inside the star to satisfy the condition˜≫ 1, Eq. (5.39) gives the following relation where th is called the thin-shell parameter [344,343]. Neglecting second-order terms with respect to Δ˜/˜and 1/(˜) in Eq. (5.40), it follows that where eff is the effective coupling given by Since th ≪ 1 under the conditions Δ˜/˜≪ 1 and 1/(˜) ≪ 1, the amplitude of the effective coupling eff becomes much smaller than 1. In the original papers of Khoury and Weltman [344,343] the thin-shell solution was derived by assuming that the field is frozen with the value = in the region 0 <˜<˜1. In this case the thin-shell parameter is given by th ≃ Δ˜/˜, which is different from Eq. (5.43). However, this difference is not important because the condition Δ˜/˜≫ 1/(˜) is satisfied for most of viable models [575].

Post Newtonian parameter
We derive the bound on the thin-shell parameter from experimental tests of the post Newtonian parameter in the solar system. The spherically symmetric metric in the Einstein frame is described by [251] d˜2 Under the condition | | ≪ 1 we obtain the following relations In the following we use the approximation ≃˜, which is valid for | | ≪ 1. Using the thin-shell solution (5.44), it follows that where we have used the approximation | | ≫ | | and hence ≃ 6 Φ th / . The term in Eq. (5.48) is smaller than ( ) = / under the condition / < (6 2 th ) −1 . Provided that the field reaches the value with the distance satisfying the condition / < (6 2 th ) −1 , the metric ( ) does not change its sign for < . The post-Newtonian parameter is given by The experimental bound (5.14) can be satisfied as long as the thin-shell parameter th is much smaller than 1. If we take the distance = , the constraint (5.14) translates into where th,⊙ is the thin-shell parameter for Sun. In f (R) gravity ( = −1/ √ 6) this corresponds to th,⊙ < 2.3 × 10 −5 .

Experimental bounds from the violation of equivalence principle
Let us next discuss constraints on the thin-shell parameter from the possible violation of equivalence principle (EP). The tightest bound comes from the solar system tests of weak EP using the free-fall acceleration of Earth ( ⊕ ) and Moon ( Moon ) toward Sun [343]. The experimental bound on the difference of two accelerations is given by [616,83,617] Provided that Earth, Sun, and Moon have thin-shells, the field profiles outside the bodies are given by Eq. (5.44) with the replacement of corresponding quantities. The presence of the field Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 ( ) with an effective coupling eff gives rise to an extra acceleration, fifth = | eff ∇ ( )|. Then the accelerations ⊕ and Moon toward Sun (mass ⊙ ) are [343] where th,⊕ is the thin-shell parameter of Earth, and Φ ⊙ ≃ 2.1 × 10 −6 , Φ ⊕ ≃ 7.0 × 10 −10 , Φ Moon ≃ 3.1 × 10 −11 are the gravitational potentials of Sun, Earth and Moon, respectively. Hence the condition (5.52) translates into [134,596] th,⊕ < 8.
which corresponds to th,⊕ < 2.2 × 10 −6 in f (R) gravity. This bound provides a tighter bound on model parameters compared to (5.51).

Constraints on model parameters in f (R) gravity
We place constraints on the f (R) models given in Eqs. (4.83) and (4.84) by using the experimental bounds discussed above. In the region of high density where is much larger than , one can use the asymptotic form (5.19) to discuss local gravity constraints. Inside and outside the spherically symmetric body the effective potential eff for the model (5.19) has two minima at The bound (5.56) translates into where ≡ 1 / and 1 is the Ricci scalar at the late-time de Sitter point. In the following we consider the case in which the Lagrangian density is given by (5.19) for ≥ 1 . If we use the models (4.83) and (4.84), then there are some modifications for the estimation of 1 . However this change should be insignificant when we place constraints on model parameters.
At the de Sitter point the model (5.19) satisfies the condition = 2 +1 /[2( 2 − − 1)]. Substituting this relation into Eq. (5.58), we find For the stability of the de Sitter point we require that ( 1 ) < 1, which translates into the condition 2 > 2 2 + 3 + 1. Hence the term /[2( 2 − − 1)] in Eq. (5.59) is smaller than 0.25 for > 0. We use the approximation that 1 and are of the orders of the present cosmological density 10 −29 g/cm 3 and the baryonic/dark matter density 10 −24 g/cm 3 in our galaxy, respectively. From Eq. (5.59) we obtain the bound [134] > 0.9 . Under this condition one can see an appreciable deviation from the ΛCDM model cosmologically as decreases to the order of .
If we consider the model (4.81), it was shown in [134] that the bound (5.56) gives the constraint < 3 × 10 −10 . This means that the deviation from the ΛCDM model is very small. Meanwhile, for the models (4.83) and (4.84), the deviation from the ΛCDM model can be large even for = (1), while satisfying local gravity constraints. We note that the model (4.89) is also consistent with local gravity constraints.

Cosmological Perturbations
The f (R) theories have one extra scalar degree of freedom compared to the ΛCDM model. This feature results in more freedom for the background. As we have seen previously, a viable cosmological sequence of radiation, matter, and accelerated epochs is possible provided some conditions hold for f (R). In principle, however, one can specify any given = ( ) and solve Eqs. (2.15) and (2.16) for those ( ( )) compatible with the given ( ).
Therefore the background cosmological evolution is not in general enough to distinguish f (R) theories from other theories. Even worse, for the same ( ), there may be some different forms of f (R) which fulfill the Friedmann equations. Hence other observables are needed in order to distinguish between different theories. In order to achieve this goal, perturbation theory turns out to be of fundamental importance. More than this, perturbations theory in cosmology has become as important as in particle physics, since it gives deep insight into these theories by providing information regarding the number of independent degrees of freedom, their speed of propagation, their time-evolution: all observables to be confronted with different data sets.
The main result of the perturbation analysis in f (R) gravity can be understood in the following way. Since it is possible to express this theory into a form of scalar-tensor theory, this should correspond to having a scalar-field degree of freedom which propagates with the speed of light. Therefore no extra vector or tensor modes come from the f (R) gravitational sector. Introducing matter fields will in general increase the number of degrees of freedom, e.g., a perfect fluid will only add another propagating scalar mode and a vector mode as well. In this section we shall provide perturbation equations for the general Lagrangian density ( , ) including metric f (R) gravity as a special case.

Perturbation equations
We start with a general perturbed metric about the flat FLRW background [57,352,231,232,437] where , , , are scalar perturbations, , are vector perturbations, and ℎ is the tensor perturbations, respectively. In this review we focus on scalar and tensor perturbations, because vector perturbations are generally unimportant in cosmology [71].
For generality we consider the following action where ( , ) is a function of the Ricci scalar and the scalar field , ( ) and ( ) are functions of , and is a matter action. We do not take into account an explicit coupling between the field and matter. The action (6.2) covers not only f (R) gravity but also other modified gravity theories such as Brans-Dicke theory, scalar-tensor theories, and dilaton gravity. We define the quantity ( , ) ≡ / . Varying the action (6.2) with respect to and , we obtain the following field equations
We decompose and into homogeneous and perturbed parts, =¯+ and =¯+ , respectively. In the following we omit the bar for simplicity. The energy-momentum tensor of an ideal fluid with perturbations is where characterizes the velocity potential of the fluid. The conservation of the energy-momentum tensor (∇ = 0) holds for the theories with the action (6.2) [357].
For the action (6.2) the background equations (without metric perturbations) are given by where is given in Eq. (2.13).
For later convenience, we define the following perturbed quantities Perturbing Einstein equations at linear order, we obtain the following equations [316,317] (see also [436,566,355,438,312,313,314,492,138,33,441,328]) ]︂ , (6.14) Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where is given by We shall solve the above equations in two different contexts: (i) inflation (Section 7), and (ii) the matter dominated epoch followed by the late-time cosmic acceleration (Section 8).

Gauge-invariant quantities
Before discussing the detail for the evolution of cosmological perturbations, we construct a number of gauge-invariant quantities. This is required to avoid the appearance of unphysical modes. Let us consider the gauge transformation = + ,^= + , (6.20) where and characterize the time slicing and the spatial threading, respectively. Then the scalar metric perturbations , , and transform as [57,71,412] = −˙, (6.21) Matter perturbations such as and obey the following transformation rulê Note that the quantity is also subject to the same transformation:^= −˙. We express the scalar part of the 3-momentum energy-momentum tensor 0 as For the scalar field and the perfect fluid one has = −˙and = −( + ) , respectively. This quantity transforms as^= + ( + ) . (6.28) Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 One can construct a number of gauge-invariant quantities unchanged under the transformation (6.20): In f (R) gravity one can introduce a scalar field as in Eq. (2.31), so that ℛ = ℛ . From the gauge-invariant quantity (6.31) it is also possible to construct another gauge-invariant quantity for the matter perturbation of perfect fluids: We note that the tensor perturbation ℎ is invariant under the gauge transformation [412]. We can choose specific gauge conditions to fix the gauge degree of freedom. After fixing a gauge, two scalar variables and are determined accordingly. The Longitudinal gauge corresponds to the gauge choice^= 0 and^= 0, under which = ( +˙) and = . In this gauge one haŝ Φ =^andΨ = −^, so that the line element (without vector and tensor perturbations) is given by where we omitted the hat for perturbed quantities. The uniform-field gauge corresponds to^= 0 which fixes = /˙. The spatial threading is fixed by choosing either^= 0 or^= 0 (up to an integration constant in the former case). For this gauge choice one hasR =^. Since the spatial curvature (3) ℛ on the constant-time hypersurface is related to via the relation (3) ℛ = −4∇ 2 / 2 , the quantity ℛ is often called the curvature perturbation on the uniform-field hypersurface. We can also choose the gauge condition = 0 or^= 0.

Perturbations Generated During Inflation
Let us consider scalar and tensor perturbations generated during inflation for the theories (6.2) without taking into account the perfect fluid ( = 0). In f (R) gravity the contribution of the field such as is absent in the perturbation equations (6.11) -(6.16). One can choose the gauge condition = 0, so that ℛ = . In scalar-tensor theory in which is the function of alone (i.e., the coupling of the form ( ) without a non-linear term in ), the gauge choice = 0 leads to ℛ = . Since = , = 0 in this case, we have ℛ = ℛ = . We focus on the effective single-field theory such as f (R) gravity and scalar-tensor theory with the coupling ( ) , by choosing the gauge condition = 0 and = 0. We caution that this analysis does not cover the theory such as ℒ = ( ) + 2 [500], because the quantity depends on both and (in other words, = , + , ). In the following we write the curvature perturbations ℛ and ℛ as ℛ.

Curvature perturbations
Plugging Eq. (7.34) into Eq. (6.11), we have where we have used the background equation (6.7). Plugging Eqs. (7.34) and (7.35) into Eq. (7.36), we find that the curvature perturbation satisfies the following simple equation in Fourier spacë where is a comoving wavenumber and Introducing the variables = √ and = ℛ, Eq. (7.37) reduces to where a prime represents a derivative with respect to the conformal time = ∫︀ −1 d . In General Relativity with a canonical scalar field one has = 1 and = 1, which corresponds to =˙2/ 2 . Then the perturbation corresponds to = [− + (˙/ ) ]. In the spatially flat gauge ( = 0) this reduces to = − , which implies that the perturbation corresponds to a canonical scalar field = .
In modified gravity theories it is not clear at this stage that Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 the perturbation = √ ℛ corresponds a canonical field that should be quantized, because Eq. (7.37) is unchanged by multiplying a constant term to the quantity defined in Eq. (7.38). As we will see in Section 7.4, this problem is overcome by considering a second-order perturbed action for the theory (6.2) from the beginning.
In order to derive the spectrum of curvature perturbations generated during inflation, we introduce the following variables [315] Then the quantity can be expressed as Then the solution to Eq. (7.39) can be expressed as a linear combination of Hankel functions, (1) where 1 and 2 are integration constants. During inflation one has | | ≪ 1, so that ′′ / ≈ ( ) 2 . For the modes deep inside the Hubble radius ( ≫ , i.e., | | ≫ 1) the perturbation satisfies the standard equation of a canonical field in the Minkowski spacetime: ′′ + 2 ≃ 0. After the Hubble radius crossing ( = ) during inflation, the effect of the gravitational term ′′ / becomes important. In the super-Hubble limit ( ≪ , i.e., | | ≪ 1) the last term on the l.h.s. of Eq. (7.37) can be neglected, giving the following solution where 1 and 2 are integration constants. The second term can be identified as a decaying mode, which rapidly decays during inflation (unless the field potential has abrupt features). Hence the curvature perturbation approaches a constant value 1 after the Hubble radius crossing ( < ). In the asymptotic past ( → −∞) the solution to Eq. (7.39) is determined by a vacuum state in quantum field theory [88], as → − / √ 2 . This fixes the coefficients to be 1 = 1 and 2 = 0, giving the following solution We define the power spectrum of curvature perturbations, Using the solution (7.45), we obtain the power spectrum [317]

47)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where we have used the relations (1) Since the curvature perturbation is frozen after the Hubble radius crossing, the spectrum (7.47) should be evaluated at = . The spectral index of ℛ, which is defined by where ℛ is given in Eq. (7.42). As long as | | ( = 1, 2, 3, 4) are much smaller than 1 during inflation, the spectral index reduces to where we have ignored those terms higher than the order of 's. Provided that | | ≪ 1 the spectrum is close to scale-invariant ( ℛ ≃ 1). From Eq. (7.47) the power spectrum of curvature perturbations can be estimated as A minimally coupled scalar field in Einstein gravity corresponds to 3 = 0, 4 = 0 and = 2 / 2 , in which case we obtain the standard results ℛ − 1 ≃ −4 1 − 2 2 and ℛ ≃ 4 /(4 2˙2 ) in slow-roll inflation [573,390].

Tensor perturbations
Tensor perturbations ℎ have two polarization states, which are generally written as = +, × [391]. In terms of polarization tensors + and × they are given by If the direction of a momentum is along the -axis, the non-zero components of polarization tensors are given by + = − + = 1 and × = × = 1. For the action (6.2) the Fourier components ℎ ( = +, ×) obey the following equation [314] ℎ + This is similar to Eq. (7.37) of curvature perturbations, apart from the difference of the factor instead of . Defining new variables = √ and = ℎ / √ 16 , it follows that We have introduced the factor 16 to relate a dimensionless massless field ℎ with a massless scalar field having a unit of mass. If˙= 0, we obtain We follow the similar procedure to the one given in Section 7.1. Taking into account polarization states, the spectrum of tensor perturbations after the Hubble radius crossing is given by

55)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 which should be evaluated at the Hubble radius crossing ( = ). The spectral index of is where is given in Eq. (7.54). If | | ≪ 1, this reduces to Then the amplitude of tensor perturbations is given by We define the tensor-to-scalar ratio For a minimally coupled scalar field in Einstein gravity, it follows that ≃ −2 1 , ≃ 16 2 /( 2 pl ), and ≃ 16 1 .

The spectra of perturbations in inflation based on f (R) gravity
Let us study the spectra of scalar and tensor perturbations generated during inflation in metric f (R) gravity. Introducing the quantity = 3˙2/(2 2 ), we have 4 =¨/(˙) and Since the field kinetic term˙2 is absent, one has 2 = 0 in Eqs. (7.42) and (7.49). Under the conditions | | ≪ 1 ( = 1, 3, 4), the spectral index of curvature perturbations is given by In the absence of the matter fluid, Eq. (2.16) translates into From Eqs. (7.50) and (7.60), the amplitude of ℛ is estimated as (7.63) Using the relation 1 ≃ − 3 , the spectral index (7.57) of tensor perturbations is given by Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Let us consider the inflation model: ( ) = ( > 0). From the discussion given in Section 3.1 the slow-roll parameters ( = 1, 3, 4) are constants: In this case one can use the exact results (7.48) and (7.56) with ℛ and given in Eqs. (7.42) and (7.54) (with 2 = 0). Then the spectral indices are we obtain the scale-invariant spectra with ℛ = 1 and = 0. Even the slight deviation from = 2 leads to a rather large deviation from the scale-invariance. If = 1.7, for example, one has ℛ − 1 = = −0.13, which does not match with the WMAP 5-year constraint: For the model ( ) = + 2 /(6 2 ), the spectrum of the curvature perturbation ℛ shows some deviation from the scale-invariance. Since inflation occurs in the regime ≫ 2 and |˙| ≪ 2 , one can approximate ≃ /(3 2 ) ≃ 4 2 / 2 . Then the power spectra (7.63) and (7.58) yield where we have employed the relation 3 ≃ − 1 .
Recall that the evolution of the Hubble parameter during inflation is given by Eq. (3.9). As long as the time at the Hubble radius crossing ( = ) satisfies the condition ( 2 /6)( − ) ≪ , one can approximate ( ) ≃ . Using Eq. (3.9), the number of e-foldings from = to the end of inflation can be estimated as Then the amplitude of the curvature perturbation is given by The WMAP 5-year normalization corresponds to ℛ = (2.445 ± 0.096) × 10 −9 at the scale = 0.002 Mpc −1 [367]. Taking the typical value = 55, the mass is constrained to be Using the relation ≃ 4 2 / 2 , it follows that 4 ≃ − 1 . Hence the spectral index (7.62) reduces to For = 55 we have ℛ ≃ 0.964, which is in the allowed region of the WMAP 5-year constraint ( ℛ = 0.960 ± 0.013 at the 68% confidence level [367]). The tensor-to-scalar ratio (7.65) can be estimated as

73)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 which satisfies the current observational bound < 0.22 [367]. We note that a minimally coupled field with the potential ( ) = 2 2 /2 in Einstein gravity (chaotic inflation model [393]) gives rise to a larger tensor-to-scalar ratio of the order of 0.1. Since future observations such as the Planck satellite are expected to reach the level of = (10 −2 ), they will be able to discriminate between the chaotic inflation model and the Starobinsky's f (R) model.

The power spectra in the Einstein frame
Let us consider the power spectra in the Einstein frame. Under the conformal transformatioñ = , the perturbed metric (6.1) is transformed as We decompose the conformal factor into the background and perturbed parts, as In what follows we omit a bar from . We recall that the background quantities are transformed as Eqs.
Since the tensor perturbation is also invariant, the tensor-to-scalar ratio˜in the Einstein frame is identical to that in the Jordan frame. For example, let us consider the model ( ) = + 2 /(6 2 ). Since the action in the Einstein frame is given by Eq. (2.32), the slow-roll parameters˜3 and˜4 vanish in this frame. Using Eqs. (7.49) and (3.27), the spectral index of curvature perturbations is given by˜ℛ where we have ignored the term of the order of 1/˜2. Since˜≃ in the slow-roll limit (|˙/(2 )| ≪ 1), Eq. (7.78) agrees with the result (7.72) in the Jordan frame. Since = (d /d˜) 2 /˜2 in the Einstein frame, Eq. (7.59) gives the tensor-to-scalar ratiõ where the background equations (3.21) and (3.22) are used with slow-roll approximations. Equation (7.79) is consistent with the result (7.73) in the Jordan frame. The equivalence of the curvature perturbation between the Jordan and Einstein frames also holds for scalar-tensor theory with the Lagrangian ℒ = ( ) /(2 2 ) − (1/2) ( ) − ( ) [411,240]. For the non-minimally coupled scalar field with ( ) = 1 − 2 2 [269,241] the spectral indices of scalar and tensor perturbations have been derived by using such equivalence [366,590].

The Lagrangian for cosmological perturbations
In Section 7.1 we used the fact that the field which should be quantized corresponds to = √ ℛ. This can be justified by writing down the action (6.1) expanded at second-order in the perturbations [437]. We recall again that we are considering an effective single-field theory such as f (R) gravity and scalar-tensor theory with the coupling ( ) . Carrying out the expansion of the action (6.2) in second order, we find that the action for the curvature perturbation ℛ (either ℛ or ℛ ) is given by [311] ( where is given in Eq. (7.38). In fact, the variation of this action in terms of the field ℛ gives rise to Eq. (7.37) in Fourier space. We note that there is another approach called the Hamiltonian formalism which is also useful for the quantization of cosmological perturbations. See [237,209,208,127] for this approach in the context of f (R) gravity and modified gravitational theories.
Introducing the quantities = ℛ and = √ , the action (7.80) can be written as where a prime represents a derivative with respect to the conformal time = ∫︀ −1 d . The action (7.81) leads to Eq. (7.39) in Fourier space. The transformation of the action (7.80) to (7.81) gives rise to the effective mass 6 We have seen in Eq. (7.42) that during inflation the quantity ′′ / can be estimated as ′′ / ≃ 2( ) 2 in the slow-roll limit, so that 2 ≃ −2 2 . For the modes deep inside the Hubble radius ( ≫ ) the action (7.81) reduces to the one for a canonical scalar field in the flat spacetime. Hence the quantization should be done for the field = √ ℛ, as we have done in Section 7.1. From the action (7.81) we understand a number of physical properties in f (R) theories and scalar-tensor theories with the coupling ( ) listed below.
1. Having a standard d'Alambertian operator, the mode has speed of propagation equal to the speed of light. This leads to a standard dispersion relation = / for the high-modes in Fourier space.
2. The sign of corresponds to the sign of the kinetic energy of ℛ. The negative sign corresponds to a ghost (phantom) scalar field. In f (R) gravity (with˙= 0) the ghost appears for < 0. In Brans-Dicke theory with ( ) = 2 and ( ) = BD / [100] (where > 0) the condition for the appearance of the ghost (˙2 + 3˙2/(2 2 ) < 0) translates into BD < −3/2. In these cases one would encounter serious problems related to vacuum instability [145,161]. 6 If we define = √ ℛ and plugging it into Eq. (7.80), we obtain the perturbed action for the field after the partial integration: where √︀ − (0) = 3 and is defined in Eq. (7.82). Then, for the field , we obtain the Klein-Gordon equation = 2 in the large-scale limit ( → 0), which defines the mass in an invariant way in the FLRW background.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 3. The field has the effective mass squared given in Eq. (7.82). In f (R) gravity it can be written as where we used the background equation (2.16) to write˙in terms of and 2 . In Fourier space the perturbation obeys the equation of motion For 2 / 2 ≫ 2 , the field propagates with speed of light. For small satisfying 2 / 2 ≪ 2 , we require a positive 2 to avoid the tachyonic instability of perturbations. Recall that the viable dark energy models based on f (R) theories need to satisfy , ≪ (i.e., = , / , ≪ 1) at early times, in order to have successful cosmological evolution from radiation domination till matter domination. At these epochs the mass squared is approximately given by

Observational Signatures of Dark Energy Models in f (R) Theories
In this section we discuss a number of observational signatures of dark energy models based on metric f (R) gravity. Our main interest is to distinguish these models from the ΛCDM model observationally. In particular we study the evolution of matter density perturbations as well as the gravitational potential to confront f (R) models with the observations of large-scale structure (LSS) and Cosmic Microwave Background (CMB). The effect on weak lensing will be discussed in Section 13.1 in more general modified gravity theories including f (R) gravity.

Matter density perturbations
Let us consider the perturbations of non-relativistic matter with the background energy density and the negligible pressure ( = 0). In Fourier space Eqs. (6.17) and (6.18) givė where ≡ − and we used the relation = 3( −˙) + ( 2 / 2 ) . In the following we consider the evolution of perturbations in f (R) gravity in the Longitudinal gauge (6.33). Since = 0, = Φ, = −Ψ, and = 3( Φ +Ψ) in this case, Eqs. (6.11), (6.13), (6.15), and (8.89) give a minimally coupled scalar field [567], which was numerically confirmed in [403]. This was further extended to scalar-tensor theories [93,171,586] and f (R) gravity [586,597]. Precisely speaking, in f (R) gravity, this approximation corresponds to We also define the effective gravitational potential This quantity characterizes the deviation of light rays, which is linked with the Integrated Sachs-Wolfe (ISW) effect in CMB [544] and weak lensing observations [27]. From Eq. (8.97) we have From Eq. (6.12) the term is of the order of 2 Φ/( 2 ) provided that the deviation from the ΛCDM model is not significant. Using Eq. (8.97) we find that the ratio 3 /( / ) is of the order of ( / ) 2 , which is much smaller than unity for sub-horizon modes. Then the gaugeinvariant perturbation given in Eq. (8.88) can be approximated as ≃ / . Neglecting the r.h.s. of Eq. (8.93) relative to the l.h.s. and using Eq. (8.97) with ≃ , we get the equation for matter perturbations: where eff is the effective (cosmological) gravitational coupling defined by [586,597] eff ≡ We recall that viable f (R) dark energy models are constructed to have a large mass in the region of high density ( ≫ 0 ). During the radiation and deep matter eras the deviation parameter = , / , is much smaller than 1, so that the mass squared satisfies If grows to the order of 0.1 by the present epoch, then the mass today can be of the order of 0 . In the regimes 2 ≫ 2 / 2 and 2 ≪ 2 / 2 the effective gravitational coupling has the asymptotic forms eff ≃ / and eff ≃ 4 /(3 ), respectively. The former corresponds to the "General Relativistic (GR) regime" in which the evolution of mimics that in GR, whereas the Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 latter corresponds to the "scalar-tensor regime" in which the evolution of is non-standard. For the f (R) models (4.83) and (4.84) the transition from the former regime to the latter regime, which is characterized by the condition 2 = 2 / 2 , can occur during the matter domination for the wavenumbers relevant to the matter power spectrum [306,568,587,270,589].
In order to derive Eq. (8.100) we used the approximation that the time-derivative terms of on the l.h.s. of Eq. (8.92) is neglected. In the regime 2 ≫ 2 / 2 , however, the large mass can induce rapid oscillations of . In the following we shall study the evolution of the oscillating mode [568]. For sub-horizon perturbations Eq. (8.92) is approximately given bÿ + 3˙+ The solution of this equation is the sum of the matter induce mode ind ≃ ( 2 /3) /( 2 / 2 + 2 ) and the oscillating mode osc satisfying As long as the frequency = √︀ 2 / 2 + 2 satisfies the adiabatic condition |˙| ≪ 2 , we obtain the solution of Eq. (8.104) under the WKB approximation: where is a constant. Hence the solution of the perturbation is expressed by [568,587] ≃ 1 3 , For viable f (R) models, the scale factor and the background Ricci scalar (0) evolve as ∝ 2/3 and (0) ≃ 4/(3 2 ) during the matter era. Then the amplitude of osc relative to (0) has the time-dependence The f (R) models (4.83) and (4.84) behave as ( ) = (− − 1) with = 2 + 1 in the regime ≫ . During the matter-dominated epoch the mass evolves as ∝ −( +1) . In the regime 2 ≫ 2 / 2 one has | osc |/ (0) ∝ −(3 +1)/2 and hence the amplitude of the oscillating mode decreases faster than (0) . However the contribution of the oscillating mode tends to be more important as we go back to the past. In fact, this behavior was confirmed in the numerical simulations of [587,36]. This property persists in the radiation-dominated epoch as well. If the condition | | < (0) is violated, then can be negative such that the condition , > 0 or , > 0 is violated for the models (4.83) and (4.84). Thus we require that | | is smaller than (0) at the beginning of the radiation era. This can be achieved by choosing the constant in Eq. (8.106) to be sufficiently small, which amounts to a fine tuning for these models.
For the models (4.83) and (4.84) one has = 1−2 ( / ) −2 −1 in the regime ≫ . Then the field defined in Eq. (2.31) rapidly approaches 0 as we go back to the past. Recall that in the Einstein frame the effective potential of the field has a potential minimum around = 0 because of the presence of the matter coupling. Unless the oscillating mode of the field perturbation is strongly suppressed relative to the background field (0) , the system can access the curvature singularity at = 0 [266]. This is associated with the condition | | < (0) discussed above. This Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 curvature singularity appears in the past, which is not related to the future singularities studied in [461,54]. The past singularity can be cured by taking into account the 2 term [37], as we will see in Section 13.3. We note that the f (R) models proposed in [427] [e.g., ( ) = − ln(1+ / )] to cure the singularity problem satisfy neither the local gravity constraints [580] nor observational constraints of large-scale structure [194].
As long as the oscillating mode osc is negligible relative to the matter-induced mode ind , we can estimate the evolution of matter perturbations as well as the effective gravitational potential Φ eff . Note that in [192,434] the perturbation equations have been derived without neglecting the oscillating mode. As long as the condition | osc | < | ind | is satisfied initially, the approximate equation (8.100) is accurate to reproduce the numerical solutions [192,589]. Equation (8.100) can be written as . The matter-dominated epoch corresponds to eff = 0 and Ω = 1. In the regime 2 ≫ 2 / 2 the evolution of and Φ eff during the matter dominance is given by where we used Eq. (8.99). The matter-induced mode ind relative to the background Ricci scalar (0) evolves as | ind |/ (0) ∝ 2/3 ∝ . At late times the perturbations can enter the regime 2 ≪ 2 / 2 , depending on the wavenumber and the mass . When 2 ≪ 2 / 2 , the evolution of and Φ eff during the matter era is [568] ∝ ( For the model ( ) = (− − 1) , the evolution of the matter-induced mode in the region 2 ≪ 2 / 2 is given by | ind |/ (0) ∝ −2 +( √ 33−5)/6 . This decreases more slowly relative to the ratio | osc |/ (0) [587], so the oscillating mode tends to be unimportant with time.

The impact on large-scale structure
We have shown that the evolution of matter perturbations during the matter dominance is given by ∝ 2/3 for 2 ≫ 2 / 2 (GR regime) and ∝ ( √ 33−1)/6 for 2 ≪ 2 / 2 (scalar-tensor regime), respectively. The existence of the latter phase gives rise to the modification to the matter power spectrum [146,74,544,526,251] (see also [597,493,494,94,446,278,435] for related works). The transition from the GR regime to the scalar-tensor regime occurs at 2 = 2 / 2 . If it occurs during the matter dominance ( ≃ 3 2 ), the condition 2 = 2 / 2 translates into [589] where we have used the relation 2 ≃ /(3 ) (valid for ≪ 1). We are interested in the wavenumbers relevant to the linear regime of the galaxy power spectrum [577,578]: If the transition characterized by the condition (8.111) occurs during the deep matter era ( ≫ 1), we can estimate the critical redshift at the transition point. In the following let us consider the models (4.83) and (4.84). In addition to the approximations 2 ≃ 2 0 Ω (0) (1 + ) 3 and ≃ 3 2 during the matter dominance, we use the the asymptotic forms ≃ (− − 1) 2 +1 and ≃ −1 − / with = 2 (2 + 1)/ 2 . Since the dark energy density today can be approximated as We caution that, when is close to Λ (the redshift at = Λ ), the estimation (8.115) begins to lose its accuracy. The ratio of the two power spectra today, i.e., ( 0 )/ ΛCDM ( 0 ) is in general different from Eq. (8.115). However, numerical simulations in [587] show that the difference is small for of the order of unity.
The modified evolution (8.110) of the effective gravitational potential for < leads to the integrated Sachs-Wolfe (ISW) effect in CMB anisotropies [544,382,545]. However this is limited to very large scales (low multipoles) in the CMB spectrum. Meanwhile the galaxy power spectrum is directly affected by the non-standard evolution of matter perturbations. From Eq. (8.115) there should be a difference between the spectral indices of the CMB spectrum and the galaxy power spectrum on the scale (8.112) [568]: Observationally we do not find any strong signature for the difference of slopes of the two spectra. If we take the mild bound Δ < 0.05, we obtain the constraint > 2. Note that in this case the local gravity constraint (5.60) is also satisfied. In order to estimate the growth rate of matter perturbations, we introduce the growth index defined by [484] ≡˙= (Ω ) , (8.117) Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 whereΩ = 2 /(3 2 ) = Ω . This choice ofΩ comes from writing Eq. (4.59) in the form we have ignored the contribution of radiation. Since the viable f (R) models are close to the ΛCDM model in the region of high density, the quantity approaches 1 in the asymptotic past. Defining DE andΩ in the above way, the Friedmann equation can be cast in the usual GR form with non-relativistic matter and dark energy [568,270,589].
The growth index in the ΛCDM model corresponds to ≃ 0.55 [612,395], which is nearly constant for 0 < < 1. In f (R) gravity, if the perturbations are in the GR regime ( 2 ≫ 2 / 2 ) today, is close to the GR value. Meanwhile, if the transition to the scalar-tensor regime occurred at the redshift larger than 1, the growth index becomes smaller than 0.55 [270]. Since 0 <Ω < 1, the smaller implies a larger growth rate. In Figure 4 we plot the evolution of the growth index in the model (4.83) with = 1 and = 1.55 for a number of different wavenumbers. In this case the present value of is degenerate around 0 ≃ 0.41 independent of the scales of our interest. For the wavenumbers = 0.1 ℎ Mpc −1 and = 0.01 ℎ Mpc −1 the transition redshifts correspond to = 5.2 and = 2.7, respectively. Hence these modes have already entered the scalar-tensor regime by today.
From Eq. (8.114) we find that gets smaller for larger and . If the mode = 0.2 ℎ Mpc −1 crossed the transition point at > (1) and the mode = 0.01 ℎ Mpc −1 has marginally entered (or has not entered) the scalar-tensor regime by today, then the growth indices should be strongly dispersed. For sufficiently large values of and one can expect that the transition to the regime 2 ≪ 2 / 2 has not occurred by today. The following three cases appear depending on the values of and [589]:    Figure 4), in which case the mode = 0.01 ℎ Mpc −1 entered the scalar-tensor regime for > (1). The regions (i), (ii), (iii) can be found numerically by solving the perturbation equations. In Figure 5 we plot those regions for the model (4.84) together with the bounds coming from the local gravity constraints as well as the stability of the late-time de Sitter point. Note that the result in the model (4.83) is also similar to that in the model (4.84). The parameter space for 3 and = (1) is dominated by either the region (ii) or the region (iii). While the present observational constraint on is quite weak, the unusual converged or dispersed spectra found above can be useful to distinguish metric f (R) gravity from the ΛCDM model in future observations. We also note that for other viable f (R) models such as (4.89) the growth index today can be as small as 0 ≃ 0.4 [589]. If future observations detect such unusually small values of 0 , this can be a smoking gun for f (R) models.

Non-linear matter perturbations
So far we have discussed the evolution of linear perturbations relevant to the matter spectrum for the scale 0.01 -0.2 ℎ Mpc −1 . For smaller scale perturbations the effect of non-linearity becomes important. In GR there are some mapping formulas from the linear power spectrum to Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 the non-linear power spectrum such as the halo fitting by Smith et al. [540]. In the halo model the non-linear power spectrum ( ) is defined by the sum of two pieces [169]: where ( ) is a linear power spectrum and In modified gravity theories, Hu and Sawicki (HS) [307] provided a fitting formula to describe a non-linear power spectrum based on the halo model. The mass function d /d ln and the halo profile depend on the root-mean-square ( ) of a linear density field. The Sheth-Tormen mass function [535] and the Navarro-Frenk-White halo profile [449] are usually employed in GR. Replacing for GR obtained in the GR dark energy model that follows the same expansion history as the modified gravity model, we obtain a non-linear power spectrum ( ) according to Eq. (8.118). In [307] this non-linear spectrum is called ∞ ( ). It is also possible to obtain a nonlinear spectrum 0 ( ) by applying a usual (halo) mapping formula in GR to modified gravity. This approach is based on the assumption that the growth rate in the linear regime determines the nonlinear spectrum. Hu and Sawicki proposed a parametrized non-linear spectrum that interpolates between two spectra ∞ ( ) and 0 ( ) [307]: where nl is a parameter which controls whether ( ) is close to 0 ( ) or ∞ ( ). In [307] they have taken the form Σ 2 ( ) = 3 ( )/(2 2 ). The validity of the HS fitting formula (8.120) should be checked with -body simulations in modified gravity models. In [478,479,529] -body simulations were carried out for the f (R) model (4.83) with = 1/2 (see also [562,379] for -body simulations in other modified gravity models). The chameleon mechanism should be at work on small scales (solar-system scales) for the consistency with local gravity constraints. In [479] it was found that the chameleon mechanism tends to suppress the enhancement of the power spectrum in the non-linear regime that corresponds to the recovery of GR. On the other hand, in the post Newtonian intermediate regime, the power spectrum is enhanced compared to the GR case at the measurable level.
Koyama et al. [371] studied the validity of the HS fitting formula by comparing it with the results of -body simulations. Note that in this paper the parametrization (8.120) was used as a fitting formula without employing the halo model explicitly. In their notation 0 corresponds to " non−GR " derived without non-linear interactions responsible for the recovery of GR (i.e., gravity is modified down to small scales in the same manner as in the linear regime), whereas ∞ corresponds to " GR " obtained in the GR dark energy model following the same expansion history as that in the modified gravity model. Note that nl characterizes how the theory approaches GR by the chameleon mechanism. Choosing Σ as where is the linear power spectrum in the modified gravity model, they showed that, in the f (R) model (4.83) with = 1/2, the formula (8.120) can fit the solutions in perturbation theory Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Figure 6: Comparison between -body simulations and the two fitting formulas in the f (R) model (4.83) with = 1/2. The circles and triangles show the results of -body simulations with and without the chameleon mechanism, respectively. The arrow represents the maximum value of (= 0.08 ℎ Mpc −1 ) by which the perturbation theory is valid. (Left) The fitting formula by Smith et al. [540] is used to predict non−GR and GR. The solid and dashed lines correspond to the power spectra with and without the chameleon mechanism, respectively. For the chameleon case nl ( ) is determined by the perturbation theory with nl ( = 0) = 0.085. (Right) The -body results in [479] are interpolated to derive non−GR without the chameleon mechanism. The obtained non−GR is used for the HS fitting formula to derive the power spectrum in the chameleon case. From [371].
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 very well by allowing the time-dependence of the parameter nl in terms of the redshift . In the regime 0 < < 1 the parameter nl is approximately given by nl ( = 0) = 0.085.
In the left panel of Figure 6 the relative difference of the non-linear power spectrum ( ) from the GR spectrum GR ( ) is plotted as a dashed curve ("no chameleon" case with nl = 0) and as a solid curve ("chameleon" case with non-zero nl derived in the perturbative regime). Note that in this simulation the fitting formula by Smith et al. [540] is used to obtain the non-linear power spectrum from the linear one. The agreement with -body simulations is not very good in the non-linear regime ( > 0.1 ℎ Mpc −1 ). In [371] the power spectrum non−GR in the no chameleon case (i.e., nl = 0) was derived by interpolating the -body results in [479]. This is plotted as the dashed line in the right panel of Figure 6. Using this spectrum non−GR for nl ̸ = 0, the power spectrum in -body simulations in the chameleon case can be well reproduced by the fitting formula (8.120) for the scale < 0.5 ℎ Mpc −1 (see the solid line in Figure 6). Although there is some deviation in the regime > 0.5 ℎ Mpc −1 , we caution that -body simulations have large errors in this regime. See [530] for clustered abundance constraints on the f (R) model (4.83) derived by the calibration of -body simulations.
In the quasi non-linear regime a normalized skewness, 3 = ⟨ 3 ⟩/⟨ 2 ⟩ 2 , of matter perturbations can provide a good test for the picture of gravitational instability from Gaussian initial conditions [79]. If large-scale structure grows via gravitational instability from Gaussian initial perturbations, the skewness in a universe dominated by pressureless matter is known to be 3 = 34/7 in GR [484]. In the ΛCDM model the skewness depends weakly on the expansion history of the universe (less than a few percent) [335]. In f (R) dark energy models the difference of the skewness from the ΛCDM model is only less than a few percent [576], even if the growth rate of matter perturbations is significantly different. This is related to the fact that in the Einstein frame dark energy has a universal coupling = −1/ √ 6 with all non-relativistic matter, unlike the coupled quintessence scenario with different couplings between dark energy and matter species (dark matter, baryons) [30].

Cosmic Microwave Background
The effective gravitational potential (8.98) is directly related to the ISW effect in CMB anisotropies. This contributes to the temperature anisotropies today as an integral [308,214] where is the optical depth, = ∫︀ −1 d is the conformal time with the present value 0 , and ℓ [ ( 0 − )] is the spherical Bessel function for CMB multipoles ℓ and the wavenumber . In the limit ℓ ≫ 1 (i.e., small-scale limit) the spherical Bessel function has a dependence ℓ ( ) ≃ (1/ℓ)( /ℓ) ℓ−1/2 , which is suppressed for large ℓ. Hence the dominant contribution to the ISW effect comes from the low ℓ modes (ℓ = (1)). In the ΛCDM model the effective gravitational potential is constant during the matter dominance, but it begins to decay after the Universe enters the epoch of cosmic acceleration (see the left panel of Figure 7). This late-time variation of Φ eff leads to the contribution to Θ ISW , which works as the ISW effect.
For viable f (R) dark energy models the evolution of Φ eff during the early stage of the matter era is constant as in the ΛCDM model. After the transition to the scalar-tensor regime, the effective gravitational potential evolves as Φ eff ∝ (  The CMB power spectrum ℓ(ℓ + 1) ℓ /(2 ) for the ΛCDM model and f (R) models with 0 = 0.5, 1.5, 3.0, 5.0. As 0 increases, the ISW contributions to low multipoles decrease, reach the minimum around 0 = 1.5, and then increase. The black points correspond to the WMAP 3-year data [561]. From [545].
model. In order to quantify the difference from the ΛCDM model at the level of perturbations, [628,544,545] defined the following quantity  Figure 7 shows that, as we increase the values of today (= 0 ), the evolution of Φ eff at late times tends to be significantly different from that in the ΛCDM model. This comes from the fact that, for increasing , the transition to the scalar-tensor regime occurs earlier.
From the right panel of Figure 7 we find that, as 0 increases, the CMB spectrum for low multipoles first decreases and then reaches the minimum around 0 = 1.5. This comes from the reduction in the decay rate of Φ eff relative to the ΛCDM model, see the left panel of Figure 7. Around 0 = 1.5 the effective gravitational potential is nearly constant, so that the ISW effect is almost absent (i.e., Θ ISW ≈ 0). For 0 1.5 the evolution of Φ eff turns into growth. This leads to the increase of the large-scale CMB spectrum, as 0 increases. The spectrum in the case 0 = 3.0 is similar to that in the ΛCDM model. The WMAP 3-year data rule out 0 > 4.3 at the 95% confidence level because of the excessive ISW effect [545].
There is another observational constraint coming from the angular correlation between the CMB temperature field and the galaxy number density field induced by the ISW effect [544]. The f (R) models predict that, for 0 1, the galaxies are anticorrelated with the CMB because of the sign change of the ISW effect. Since the anticorrelation has not been observed in the observational data of CMB and LSS, this places an upper bound of 0 1 [545]. This is tighter than the bound Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 0 < 4.3 coming from the CMB angular spectrum discussed above.
Finally we briefly mention stochastic gravitational waves produced in the early universe [421,172,122,123,174,173,196,20]. For the inflation model ( ) = + 2 /(6 2 ) the primordial gravitational waves are generated with the tensor-to-scalar ratio of the order of 10 −3 , see Eq. (7.73). It is also possible to generate stochastic gravitational waves after inflation under the modification of gravity. Capozziello et al. [122,123] studied the evolution of tensor perturbations for a toy model = 1+ in the FLRW universe with the power-law evolution of the scale factor. Since the parameter is constrained to be very small (| | < 7.2 × 10 −19 ) [62,160], it is very difficult to detect the signature of f (R) gravity in the stochastic gravitational wave background. This property should hold for viable f (R) dark energy models in general, because the deviation from GR during the radiation and the deep matter era is very small.

Palatini Formalism
In this section we discuss f (R) theory in the Palatini formalism [481]. In this approach the action (2.1) is varied with respect to both the metric and the connection Γ . Unlike the metric approach, and Γ are treated as independent variables. Variations using the Palatini approach [256,607,608,261,262,260] lead to second-order field equations which are free from the instability associated with negative signs of , [422,423]. We note that even in the 1930s Lanczos [378] proposed a specific combination of curvature-squared terms that lead to a secondorder and divergence-free modified Einstein equation.
The background cosmological dynamics of Palatini f (R) gravity has been investigated in [550,553,21,253,495], which shows that the sequence of radiation, matter, and accelerated epochs can be realized even for the model ( ) = − / with > 0 (see also [424,457,495]). The equations for matter density perturbations were derived in [359]. Because of a large coupling between dark energy and non-relativistic matter dark energy models based on Palatini f (R) gravity are not compatible with the observations of large-scale structure, unless the deviation from the ΛCDM model is very small [356,386,385,597]. Such a large coupling also gives rise to nonperturbative corrections to the matter action, which leads to a conflict with the Standard Model of particle physics [261,262,260] (see also [318,472,473,475,55]).
There are also a number of works [470,471,216,552] about the Newtonian limit in the Palatini formalism (see also [18,19,107,331,511,510]). In particular it was shown in [55,56] that the nondynamical nature of the scalar-field degree of freedom can lead to a divergence of non-vacuum static spherically symmetric solutions at the surface of a compact object for commonly-used polytropic equations of state. Hence Palatini f (R) theory is difficult to be compatible with a number of observations and experiments, as long as the models are constructed to explain the late-time cosmic acceleration. Moreover it is also known that in Palatini gravity the Cauchy problem [609] is not well-formulated due to the presence of higher derivatives of matter fields in field equations [377] (see also [520,135] for related works). We also note that the matter Lagrangian (such as the Lagrangian of Dirac particles) cannot be simply assumed to be independent of connections. Even in the presence of above mentioned problems it will be useful to review this theory because we can learn the way of modifications of gravity from GR to be consistent with observations and experiments.

3)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where a prime represents a derivative in terms of ( ). The variation of the action (2.1) with respect to the connection leads to the following equation In Einstein gravity ( ( ) = − 2Λ and ( ) = 1) the field equations (9.2) and (9.4) are identical to the equations (2.7) and (2.4), respectively. However, the difference appears for the f (R) models which include non-linear terms in . While the kinetic term is present in Eq. (2.7), such a term is absent in Palatini f (R) gravity. This has the important consequence that the oscillatory mode, which appears in the metric formalism, does not exist in the Palatini formalism. As we will see later on, Palatini f (R) theory corresponds to Brans-Dicke (BD) theory [100] with a parameter BD = −3/2 in the presence of a field potential. Such a theory should be treated separately, compared to BD theory with BD ̸ = −3/2 in which the field kinetic term is present.
As we have derived the action (2.21) from (2.18), the action in Palatini f (R) gravity is equivalent to Since the derivative of in terms of is , = /(2 2 ), we obtain the following relation from Eq. (9.2): 4 − 2 , = . (9.7) Using the relation (9.3), the action (9.5) can be written as Comparing this with Eq. (2.23) in the unit 2 = 1, we find that Palatini f (R) gravity is equivalent to BD theory with the parameter BD = −3/2 [262,470,551]. As we will see in Section 10.1, this equivalence can be also seen by comparing Eqs. (9.1) and (9.4) with those obtained by varying the action (2.23) in BD theory. In the above discussion we have implicitly assumed that ℒ does not explicitly depend on the Christoffel connections Γ . This is true for a scalar field or a perfect fluid, but it is not necessarily so for other matter Lagrangians such as those describing vector fields. There is another way for taking the variation of the action, known as the metric-affine formalism [299,558,557,121]. In this formalism the matter action depends not only on the metric but also on the connection Γ . Since the connection is independent of the metric in this approach, one can define the quantity called hypermomentum [299], as Δ ≡ (−2/ √ − ) ℒ / Γ . The usual assumption that the connection is symmetric is also dropped, so that the antisymmetric quantity called the Cartan torsion tensor, ≡ Γ [ ] , is defined. The non-vanishing property of allows the presence of torsion in this theory. If the condition Δ [ ] = 0 holds, it follows that the Cartan torsion tensor vanishes ( = 0) [558]. Hence the torsion is induced by matter fields with the anti-symmetric hypermomentum. The f (R) Palatini gravity belongs to f (R) theories in the metric-affine formalism with Δ = 0. In the following we do not discuss further f (R) theory in the metric-affine formalism. Readers who are interested in those theories may refer to the papers [557,556].

Background cosmological dynamics
We discuss the background cosmological evolution of dark energy models based on Palatini f (R) gravity. We shall carry out general analysis without specifying the forms of f (R). We take into account non-relativistic matter and radiation whose energy densities are and , respectively. In the flat FLRW background (2.12) we obtain the following equations 6 (︃ +2 together with the continuity equations,˙+ 3 = 0 and˙+ 4 = 0. Combing Eqs. (9.9) and (9.10) together with continuity equations, it follows thaṫ In order to discuss cosmological dynamics it is convenient to introduce the dimensionless variables: by which Eq. (9.12) can be written as Differentiating 1 and 2 with respect to = ln , we obtain [253] where The following constraint equation also holds Hence the Ricci scalar can be expressed in terms of 1 and 2 . Differentiating Eq. (9.11) with respect to , it follows thaṫ
In Figure 8 we plot the evolution of eff as well as 1 and 2 for the model ( ) = − / with = 0.02. This shows that the sequence of ( ) radiation domination ( eff = 1/3), ( ) matter domination ( eff = 0), and de Sitter acceleration ( eff = −1) is realized. Recall that in Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 does not matter because there is no propagating degree of freedom with a mass associated with the second derivative , [554].
In [21,253] the dark energy model ( ) = − / was constrained by the combined analysis of independent observational data. From the joint analysis of Super-Nova Legacy Survey [39], BAO [227] and the CMB shift parameter [561], the constraints on two parameters and are ∈ [−0.23, 0.42] and ∈ [2.73, 10.6] at the 95% confidence level (in the unit of 0 = 1) [253]. Since the allowed values of are close to 0, the above model is not particularly favored over the ΛCDM model. See also [116,148,522,46,47] for observational constraints on f (R) dark energy models based on the Palatini formalism.

Matter perturbations
We have shown that f (R) theory in the Palatini formalism can give rise to the late-time cosmic acceleration preceded by radiation and matter eras. In this section we study the evolution of matter density perturbations to confront Palatini f (R) gravity with the observations of large-scale structure [359,356,357,598,380,597]. Let us consider the perturbation of non-relativistic matter with a homogeneous energy density . Koivisto and Kurki-Suonio [359] derived perturbation equations in Palatini f (R) gravity. Using the perturbed metric (6.1) with the same variables Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 as those introduced in Section 6, the perturbation equations are given by   In the above approximation we do not need to worry about the dominance of the oscillating mode of perturbations in the past. Note also that the same approximate equation of as Eq. (9.35) can be derived for different gauge choices [597].
The parameter is a crucial quantity to characterize the evolution of perturbations. This quantity can be estimated as ≈ ( / ) 2 , which is much larger than for sub-horizon modes ( ≫ ). In the regime ≪ 1 the matter perturbation evolves as ∝ 2/3 . Meanwhile the evolution of in the regime ≫ 1 is completely different from that in GR. If the transition characterized by = 1 occurs before today, this gives rise to the modification to the matter spectrum compared to the GR case.
In the regime ≫ 1, let us study the evolution of matter perturbations during the matter dominance. We shall consider the case in which the parameter (with | | ≪ 1) evolves as ∝ , where is a constant. For the model ( ) = − ( / ) ( < 1) the power corresponds to = 1 + , whereas for the models (4.83) and (4.84) with > 0 one has = 1 + 2 . During the matter dominance the parameter evolves as = ±( / ) 2 +2/3 , where the subscript " " denotes the value at which the perturbation crosses = ±1. Here + and − signs correspond to the cases > 0 and < 0, respectively. Then the matter perturbation equation (9.35) reduces to When > 0, the growing mode solution to Eq. (9.38) is given by This shows that the perturbations exhibit violent growth for > −1/3, which is not compatible with observations of large-scale structure. In metric f (R) gravity the growth of matter perturbations is much milder. When < 0, the perturbations show a damped oscillation: where = √ 6 (3 +1)( − )/2 /(3 + 1), and is a constant. The averaged value of the growth rate is given by¯= −(3 + 2)/4, but it shows a divergence every time changes by . These negative values of are also difficult to be compatible with observations. The f (R) models can be consistent with observations of large-scale structure if the universe does not enter the regime | | > 1 by today. This translates into the condition [597] | ( = 0)| ( 0 0 / ) 2 . (9.41) Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Let us consider the wavenumbers 0.01 ℎ Mpc −1 0.2 ℎ Mpc −1 that corresponds to the linear regime of the matter power spectrum. Since the wavenumber = 0.2 ℎ Mpc −1 corresponds to ≈ 600 0 0 (where "0" represents present quantities), the condition (9.41) gives the bound | ( = 0)| 3 × 10 −6 .
If we use the observational constraint of the growth rate, 1.5 [418,605,211], then the deviation parameter today is constrained to be | ( = 0)| 10 −5 -10 −4 for the model ( ) = − ( / ) ( < 1) as well as for the models (4.83) and (4.84) [597]. Recall that, in metric f (R) gravity, the deviation parameter can grow to the order of 0.1 by today. Meanwhile f (R) dark energy models based on the Palatini formalism are hardly distinguishable from the ΛCDM model [356,386,385,597]. Note that the bound on ( = 0) becomes even severer by considering perturbations in non-linear regime. The above peculiar evolution of matter perturbations is associated with the fact that the coupling between non-relativistic matter and a scalar-field degree of freedom is very strong (as we will see in Section 10.1).
The above results are based on the fact that dark matter is described by a cold and perfect fluid with no pressure. In [358] it was suggested that the tight bound on the parameter can be relaxed by considering imperfect dark matter with a shear stress. Although the approach taken in [358] did not aim to explain the origin of a dark matter stress Π that cancels the -dependent term in Eq. (9.35), it will be of interest to further study whether some theoretically motivated choice of Π really allows the possibility that Palatini f (R) dark energy models can be distinguished from the ΛCDM model.

Shortcomings of Palatini f (R) gravity
Since ≈ 2 2 it follows that ≈ 2 2 /( 2 2 pl ). Perturbing the Jordan-frame action (9.8) [which is equivalent to the action in Palatini f (R) gravity] to second-order and using the solution Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 ≈ 1 + 2 2 /( 2 2 pl ), we find that the effective action of the Higgs field for an energy scale much lower than (= 100 -1000 GeV) is given by [55] ≃ ∫︁ Since ≈ for ≪ , the correction term can be estimated as In order to give rise to the late-time acceleration we require that ≈ 0 ≈ 10 −42 GeV. For the Higgs mass = 100 GeV it follows that ≈ 10 56 ≫ 1. This correction is too large to be compatible with the Standard Model of particle physics.
The above result is based on the models ( ) = − 2( +1) / with = (1). Having a look at Eq. (9.44), the only way to make the perturbation small is to choose very close to 0. This means that the deviation from the ΛCDM model is extremely small (see [388] for a related work). In fact, this property was already found by the analysis of matter density perturbations in Section 9.3. While the above analysis is based on the calculation in the Jordan frame in which test particles follow geodesics [55], the same result was also obtained by the analysis in the Einstein frame [261,262,260,318].
Another unusual property of Palatini f (R) gravity is that a singularity with the divergent Ricci scalar can appear at the surface of a static spherically symmetric star with a polytropic equation of state = Γ 0 with 3/2 < Γ < 2 (where is the pressure and 0 is the rest-mass density) [56,55] (see also [107,331]). Again this problem is intimately related with the particular algebraic dependence (9.2) in Palatini f (R) gravity. In [56] it was claimed that the appearance of the singularity does not very much depend on the functional forms of f (R) and that the result is not specific to the choice of the polytropic equation of state.
The Palatini gravity has a close relation with an effective action which reproduces the dynamics of loop quantum cosmology [477]. [474] showed that the model ( ) = + 2 /(6 2 ), where is of the order of the Planck mass, is not plagued by a singularity problem mentioned above, while the singularity typically arises for the f (R) models constructed to explain the late-time cosmic acceleration (see also [504] for a related work). Since Planck-scale corrected Palatini f (R) models may cure the singularity problem, it will be of interest to understand the connection with quantum gravity around the cosmological singularity (or the black hole singularity). In fact, it was shown in [60] that non-singular bouncing solutions can be obtained for power-law f (R) Lagrangians with a finite number of terms.
Finally we note that the extension of Palatini f (R) gravity to more general theories including Ricci and Riemann tensors was carried out in [384,387,95,236,388,509,476]. While such theories are more involved than Palatini f (R) gravity, it may be possible to construct viable modified gravity models of inflation or dark energy.

Extension to Brans-Dicke Theory
So far we have discussed f (R) gravity theories in the metric and Palatini formalisms. In this section we will see that these theories are equivalent to Brans-Dicke (BD) theory [100] in the presence of a scalar-field potential, by comparing field equations in f (R) theories with those in BD theory. It is possible to construct viable dark energy models based on BD theory with a constant parameter BD . We will discuss cosmological dynamics, local gravity constraints, and observational signatures of such generalized theory.

Brans-Dicke theory and the equivalence with f (R) theories
Let us start with the following 4-dimensional action in BD theory where BD is the BD parameter, ( ) is a potential of the scalar field , and is a matter action that depends on the metric and matter fields Ψ . In this section we use the unit 2 = 8 = 1/ 2 pl = 1, but we recover the gravitational constant and the reduced Planck mass pl when the discussion becomes transparent. The original BD theory [100] does not possess the field potential ( ).
Taking the variation of the action (10.1) with respect to and , we obtain the following field equations where ( ) is the Ricci scalar in metric f (R) gravity, and is the energy-momentum tensor of matter. In order to find the relation with f (R) theories in the metric and Palatini formalisms, we consider the following correspondence Recall that this potential (which is the gravitational origin) already appeared in Eq. (2.28). We then find that Eqs.  [262,470,551]. Recall that we also showed this by rewriting the action (2.1) in the form (9.8).
One can consider more general theories called scalar-tensor theories [268] in which the Ricci scalar is coupled to a scalar field . The general 4-dimensional action for scalar-tensor theories can be written as where = / 2 . We have introduced a new scalar field to make the kinetic term canonical: We define a quantity that characterizes the coupling between the field and non-relativistic matter in the Einstein frame: Recall that, in metric f (R) gravity, we introduced the same quantity in Eq. (2.40), which is constant ( = −1/ √ 6). For theories with =constant, we obtain the following relations from Eqs. (10.7) and (10.8): In this case the action (10.5) in the Jordan frame reduces to [596] = ∫︁ (10.10) In the limit that → 0 we have ( ) → 1, so that Eq. (10.10) recovers the action of a minimally coupled scalar field in GR.
Let us compare the action (10.10) with the action (10.1) in BD theory. Setting = = −2 , the former is equivalent to the latter if the parameter BD is related to via the relation [343,596] 3 + 2 BD = 1 2 2 .
(10.11) This shows that the GR limit ( BD → ∞) corresponds to the vanishing coupling ( → 0). Since = −1/ √ 6 in metric f (R) gravity one has BD = 0, as expected. The Palatini f (R) gravity corresponds to BD = −3/2, which corresponds to the infinite coupling ( 2 → ∞). In fact, Palatini gravity can be regarded as an isolated "fixed point" of a transformation involving a special conformal rescaling of the metric [247]. In the Einstein frame of the Palatini formalism, the scalar field does not have a kinetic term and it can be integrated out. In general, this leads to a matter action which is non-linear, depending on the potential ( ). This large coupling poses a number of problems such as the strong amplification of matter density perturbations and the conflict with the Standard Model of particle physics, as we have discussed in Section 9.
Note that BD theory is one of the examples in scalar-tensor theories and there are some theories that give rise to non-constant values of . For example, the action of a nonminimally coupled scalar field with a coupling corresponds to ( ) = 1− 2 and ( ) = 1, which gives the field-dependent coupling ( ) = /[1 − 2 (1 − 6 )] 1/2 . In fact the dynamics of dark energy in such a theory has been studied by a number of authors [22,601,151,68,491,44,505]. In the following we shall focus on the constant coupling models with the action (10.10). We stress that this is equivalent to the action (10.1) in BD theory.

Cosmological dynamics of dark energy models based on Brans-Dicke theory
The first attempt to apply BD theory to cosmic acceleration is the extended inflation scenario in which the BD field is identified as an inflaton field [374,571]. The first version of the inflation model, which considered a first-order phase transition in BD theory, resulted in failure due to the graceful exit problem [375,613,65]. This triggered further study of the possibility of realizing inflation in the presence of another scalar field [394,78]. In general the dynamics of such a multifield system is more involved than that in the single-field case [71]. The resulting power spectrum of density perturbations generated during multi-field inflation in BD theory was studied by a number of authors [570,272,156,569].
In the context of dark energy it is possible to construct viable single-field models based on BD theory. In what follows we discuss cosmological dynamics of dark energy models based on the action (10.10) in the flat FLRW background given by (2.12) (see, e.g., [596,22,85,289,5,327,139,168] for dynamical analysis in scalar-tensor theories). Our interest is to find conditions under which a sequence of radiation, matter, and accelerated epochs can be realized. This depends upon the form of the field potential ( ). We first carry out general analysis without specifying the forms of the potential. We take into account non-relativistic matter with energy density and radiation with energy density . The Jordan frame is regarded as a physical frame due to the usual conservation of non-relativistic matter ( ∝ −3 ). Varying the action (10.10) with respect to and , we obtain the following equations where = −2 . We introduce the following dimensionless variables and also the density parameters Taking the derivatives of 1 , 2 and 3 with respect to = ln , we find

Name
If is a constant, i.e., for the exponential potential = 0 − , one can derive fixed points for Eqs.  Table 1 we list the fixed points of the system in the absence of radiation ( 3 = 0). Note that the radiation point corresponds to ( 1 , 2 , 3 ) = (0, 0, 1). The point (a) is the -matter-dominated epoch ( MDE) during which the density of non-relativistic matter is a non-zero constant. Provided that 2 ≪ 1 this can be used for the matter-dominated epoch. The kinetic points (b1) and (b2) are responsible neither for the matter era nor for the accelerated epoch (for | | 1). The point (c) is the scalar-field dominated solution, which can be used for the late-time acceleration for eff < −1/3. When 2 ≪ 1 this point yields the cosmic acceleration for − √ 2 + 4 < < √ 2 + 4 . The scaling solution (d) can be responsible for the matter era for | | ≪ | |, but in this case the condition eff < −1/3 for the point (c) leads to 2 2. Then the energy fraction of the pressureless matter for the point (d) does not satisfy the condition Ω ≃ 1. The point (e) gives rise to the de Sitter expansion, which exists for the special case with = 4 [which can be also regarded as the special case of the point (c)]. From the above discussion the viable cosmological trajectory for constant is the sequence from the point (a) to the scalar-field dominated point (c) under the conditions 2 ≪ 1 and − √ 2 + 4 < < √ 2 + 4 . The analysis based on the Einstein frame action (10.6) also gives rise to the MDE followed by the scalar-field dominated solution [23,22].
Let us consider the case of non-constant . The fixed points derived above may be regarded as "instantaneous" points 7 [195,454] varying with the time-scale smaller than −1 . As in metric f (R) gravity ( = −1/ √ 6) we are interested in large coupling models with | | of the order of unity. In order for the potential ( ) to satisfy local gravity constraints, the field needs to be heavy in the region ≫ 0 ∼ 2 0 such that | | ≫ 1. Then it is possible to realize the matter era by the point (d) with | | ≪ | |. Moreover the solutions can finally approach the de Sitter solution (e) with = 4 or the field-dominated solution (c). The stability of the point (e) was analyzed in [596,250,242] by considering linear perturbations 1 , 2 and . One can easily show that the 7 In strict mathematical sense the instantaneous fixed point is not formally defined because it varies with time.
During the radiation and deep matter eras one has = 6(2 2 +˙) ≃ / from Eqs. (10.12) -(10.13) by noting that 0 is negligibly small relative to the background fluid density. From Eq. (10.14) the field is nearly frozen at a value satisfying the condition , + ≃ 0. Then the field evolves along the instantaneous minima given by As long as ≫ 2 0 we have that | | ≪ 1. In this regime the slope in Eq. (10.24) is much larger than 1. The field value | | increases for decreasing and hence the slope decreases with time.
Since ≫ 1 around = 0, the instantaneous fixed point (d) can be responsible for the matterdominated epoch provided that | | ≪ . The variable = −2 decreases in time irrespective of the sign of the coupling and hence 0 < < 1. The de Sitter point is characterized by = 4 , i.e., The de Sitter solution is present as long as the solution of this equation exists in the region 0 < 1 < 1. From Eq. (10.24) the derivative of in terms of is given by The above discussion shows that for the model (10.23) the matter point (d) can be followed by the stable de Sitter solution (e) for 0 < < 1. In fact numerical simulations in [596] show that the sequence of radiation, matter and de Sitter epochs can be in fact realized. Introducing the energy density DE and the pressure DE of dark energy as we have done for metric f (R) gravity, the dark energy equation of state DE = DE / DE is given by the same form as Eq. (4.97). Since for the model (10.23) increases toward the past, the phantom equation of state ( DE < −1) as well as the cosmological constant boundary crossing ( DE = −1) occurs as in the case of metric f (R) gravity [596].
As we will see in Section 10.3, for a light scalar field, it is possible to satisfy local gravity constraints for | | 10 −3 . In those cases the potential does not need to be steep such that ≫ 1 in the region ≫ 0 . The cosmological dynamics for such nearly flat potentials have been discussed by a number of authors in several classes of scalar-tensor theories [489,451,416,271]. It is also possible to realize the condition DE < −1, while avoiding the appearance of a ghost [416,271].

Local gravity constraints
We study local gravity constraints (LGC) for BD theory given by the action (10.10). In the absence of the potential ( ) the BD parameter BD is constrained to be BD > 4 × 10 4 from solar-system experiments [616,83,617]. This bound also applies to the case of a nearly massless field with the potential ( ) in which the Yukawa correction − is close to unity (where is a scalar-field mass and is an interaction length). Using the bound BD > 4 × 10 4 in Eq. (10.11), we find that This is a strong constraint under which the cosmological evolution for such theories is difficult to be distinguished from the ΛCDM model ( = 0). If the field potential is present, the models with large couplings (| | = (1)) can be consistent with local gravity constraints as long as the mass of the field is sufficiently large in the region of high density. For example, the potential (10.23) is designed to have a large mass in the highdensity region so that it can be compatible with experimental tests for the violation of equivalence principle through the chameleon mechanism [596]. In the following we study conditions under which local gravity constraints can be satisfied for the model (10.23).
As in the case of metric f (R) gravity, let us consider a configuration in which a spherically symmetric body has a constant density inside the body with a constant density = (≪ ) outside the body. For the potential = / 2 in the Einstein frame one has , ≃ −2 0 (2 ) −1 under the condition | | ≪ 1. Then the field values at the potential minima inside and outside the body are The field mass squared 2 ≡ , at = ( = , ) is approximately given by Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Recall that 0 is roughly the same order as the present cosmological density 0 ≃ 10 −29 g/cm 3 .
The baryonic/dark matter density in our galaxy corresponds to ≃ 10 −24 g/cm 3 . The mean density of Sun or Earth is about = (1) g/cm 3 . Hence and are in general much larger than 0 for local gravity experiments in our environment. For˜≫ 1 the chameleon mechanism we discussed in Section 5.2 can be directly applied to BD theory whose Einstein frame action is given by Eq. (10.6) with = −2 .
The bound (5.56) coming from the possible violation of equivalence principle in the solar system translates into Let us consider the case in which the solutions finally approach the de Sitter point (e) in Table 1.

Evolution of matter density perturbations
Let us next study the evolution of perturbations in non-relativistic matter for the action (10.10) with the potential ( ) and the coupling ( ) = −2 . As in metric f (R) gravity, the matter perturbation satisfies Eq. (8.93) in the Longitudinal gauge. We define the field mass squared as 2 ≡ , . For the potential consistent with local gravity constraints [such as (10.23)], the mass is much larger than the present Hubble parameter 0 during the radiation and deep matter eras. Note that the condition 2 ≫ is satisfied in most of the cosmological epoch as in the case of metric f (R) gravity.
The perturbation equations for the action (10.10) are given in Eqs. (6.11) -(6.18) with = ( ) , = (1 − 6 2 ) , and = . We use the unit 2 = 1, but we restore 2 when it is necessary. In the Longitudinal gauge one has = 0, = Φ, = −Ψ, and = 3( Φ +Ψ) in these equations. Since we are interested in sub-horizon modes, we use the approximation that the terms containing 2 / 2 , , , and 2 are the dominant contributions in Eqs. (6.11) -(6.19). We shall neglect the contribution of the time-derivative terms of in Eq. (6.16). As we have discussed for metric f (R) gravity in Section 8.1, this amounts to neglecting the oscillating mode of perturbations. The initial conditions of the field perturbation in the radiation era need to be chosen so that the oscillating mode osc is smaller than the matter-induced mode ind . In Fourier space Eq. (6. 16) gives

36)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Using this relation together with Eqs. (6.13) and (6.18), it follows that Combing this equation with Eqs. (6.11) and (6.13), we obtain [596,547] (see also [84,632,631]) where the effective gravitational coupling is In the regime 2 / ≫ 2 / 2 ("GR regime") this reduces to eff = / , so that the evolution of and Φ eff during the matter domination (Ω = /(3 2 ) ≃ 1) is standard: In the regime 2 / ≪ 2 / 2 ("scalar-tensor regime") we have where we used the relation (10.11) between the coupling and the BD parameter BD . Since BD = 0 in f (R) gravity, it follows that eff = 4 /(3 ). Note that the result (10.42) agrees with the effective Newtonian gravitational coupling between two test masses [93,175]. The evolution of and Φ eff during the matter dominance in the regime 2 / ≪ 2 / 2 is Hence the growth rate of for ̸ = 0 is larger than that for = 0. As an example, let us consider the potential (10.23). During the matter era the field mass squared around the potential minimum (induced by the matter coupling) is approximately given by which decreases with time. The perturbations cross the point 2 / = 2 / 2 at time = , which depends on the wavenumber . Since the evolution of the mass during the matter domination is given by ∝ − 2− 1− , the time has a scale-dependence: . More precisely the critical redshift at time can be estimated as [596] ≃ where the subscript "0" represents present quantities. For the scales 30 0 0 600 0 0 , which correspond to the linear regime of the matter power spectrum, the critical redshift can be in the region > 1. Note that, for larger , decreases.
When < and > the matter perturbation evolves as ∝ 2/3 and ∝ ( √ 25+48 2 −1)/6 , respectively (apart from the epoch of the late-time cosmic acceleration). The matter power spectrum at time = Λ (at which¨= 0) shows a difference compared to the ΛCDM model, which is given by The CMB power spectrum is also modified by the non-standard evolution of the effective gravitational potential Φ eff for > . This mainly affects the low multipoles of CMB anisotropies through of the ISW effect. Hence there is a difference between the spectral indices of the matter power spectrum and of the CMB spectrum on the scales (0.01 ℎ Mpc −1 0.2 ℎ Mpc −1 ) [596]: Note that this covers the result (8.116) in f (R) gravity ( = −1/ √ 6 and = 2 /(2 + 1)) as a special case. Under the criterion Δ ( Λ ) < 0.05 we obtain the bounds > 0.957 for = 1 and > 0.855 for = 0.5. As long as is close to 1, the model can be consistent with both cosmological and local gravity constraints. The allowed region coming from the bounds Δ ( Λ ) < 0.05 and (10.35) are illustrated in Figure 9. The growth rate of for > is given by =˙/( ) = ( √︀ 25 + 48 2 − 1)/4. As we mentioned in Section 8, the observational bound on is still weak in current observations. If we use the criterion < 2 for the analytic estimation = ( √︀ 25 + 48 2 − 1)/4, we obtain the bound < 1.08 (see Figure 9). The current observational data on the growth rate as well as Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 its growth index is not enough to place tight bounds on and , but this will be improved in future observations. Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3

Relativistic Stars in f (R) Gravity and Chameleon Theories
In Section 5 we discussed the existence of thin-shell solutions in metric f (R) gravity in the Minkowski background, i.e., without the backreaction of metric perturbations. For the f (R) dark energy models (4.83) and (4.84), Frolov [266] anticipated that the curvature singularity at = 0 (shown in Figure 3) can be accessed in a strong gravitational background such as neutron stars. Kobayashi and Maeda [349,350] studied spherically symmetric solutions for a constant density star with a vacuum exterior and claimed the difficulty of obtaining thin-shell solutions in the presence of the backreaction of metric perturbations. In [594] thin-shell solutions were derived analytically in the Einstein frame of BD theory (including f (R) gravity) under the linear expansion of the gravitational potential Φ at the surface of the body (valid for Φ < 0.3). In fact, the existence of such solutions was numerically confirmed for the inverse power-law potential ( ) = 4+ − [594]. For the f (R) models (4.83) and (4.84), it was numerically shown that thin-shell solutions exist for Φ 0.3 by the analysis in the Jordan frame [43,600,42] (see also [167]). In particular Babichev and Langlois [43,42] constructed static relativistic stars both for constant energy density configurations and for a polytropic equation of state, provided that the pressure does not exceed one third of the energy density. Since the relativistic pressure tends to be stronger around the center of the spherically symmetric body for larger Φ , the boundary conditions at the center of the body need to be carefully chosen to obtain thin-shell solutions numerically. In this sense the analytic estimation of thin-shell solutions carried out in [594] can be useful to show the existence of static star configurations, although such analytic solutions have been so far derived only for a constant density star.
In the following we shall discuss spherically symmetric solutions in a strong gravitational background with Φ 0.3 for BD theory with the action (10.10). This analysis covers metric f (R) gravity as a special case (the scalar-field degree of freedom defined in Eq. (2.31) with = −1/ √ 6). While field equations will be derived in the Einstein frame, we can transform back to the Jordan frame to find the corresponding equations (as in the analysis of Babichev and Langlois [42]). In addition to the papers mentioned above, there are also a number of works about spherically symmetric solutions for some equation state of matter [330,332,443,444,300,533].

Field equations
We already showed that under the conformal transformation˜= −2 the action (10.10) is transformed to the Einstein frame action: Recall that in the Einstein frame this gives rise to a constant coupling between non-relativistic matter and the field . We use the unit 2 = 8 = 1, but we restore the gravitational constant when it is required. Let us consider a spherically symmetric static metric in the Einstein frame: where Ψ(˜) and Φ(˜) are functions of the distance˜from the center of symmetry. For the action (11.1) the energy-momentum tensors for the scalar field and the matter are given, respec-Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 tively, by˜( For the metric (11.2) the (00) and (11) components for the energy-momentum tensor of the field are˜0 where a prime represents a derivative with respect to˜. The energy-momentum tensor of matter in the Einstein frame is given by˜= diag (−˜,˜,˜,˜), which is related to ( ) in the Jordan frame via˜( ) = 4 ( ) . Hence it follows that˜= 4 and˜= 4 . Variation of the action (11.1) with respect to gives where ℒ = −(∇ ) 2 /2 − ( ) is the field Lagrangian density. Since the derivative of ℒ in terms of is given by Eq. (2.41), i.e., ℒ / = √ −˜(−˜+ 3˜), we obtain the equation of the field [594,42]: where a tilde represents a derivative with respect to˜. From the Einstein equations it follows that If the field potential ( ) is responsible for dark energy, we can neglect both ( ) and ′2 relative to˜in the local region whose density is much larger than the cosmological density ( 0 ∼ 10 −29 g/cm 3 ). In this case Eq. (11.8) is integrated to give 4¯2˜d¯. (11.12) Substituting Eqs. (11.8) and (11.9) into Eq. (11.7), it follows that The gravitational potential Φ around the surface of a compact object can be estimated as Φ ≈˜˜2, where˜is the mean density of the star and˜is its radius. Provided that Φ ≪ 1, Eq. (11.13) reduces to Eq. (5.15) in the Minkowski background (note that the pressure˜is also much smaller than the density˜for non-relativistic matter).

Constant density star
Let us consider a constant density star with˜=˜. We also assume that the density outside the star is constant,˜=˜. We caution that the conserved density˜( ) in the Einstein frame is given by˜( ) = −˜ [ 343]. However, since the condition ≪ 1 holds in most cases of our interest, we do not distinguish between˜( ) and˜in the following discussion.
Inside the spherically symmetric body (0 <˜<˜), Eq. (11.12) gives (11.14) Neglecting the field contributions in Eqs. (11.8) -(11.11), the gravitational background for 0 < <˜is characterized by the Schwarzschild interior solution. Then the pressure˜(˜) inside the body relative to the density˜can be analytically expressed as (11.15) where Φ is the gravitational potential at the surface of the body: (11.16) Here = 4˜3˜/3 is the mass of the spherically symmetric body. The density˜is much smaller than˜, so that the metric outside the body can be approximated by the Schwarzschild exterior solution Φ(˜) ≃˜= Φ˜,˜(˜) ≃ 0 (˜>˜) . (11.17) In the following we shall derive the analytic field profile by using the linear expansion in terms of the gravitational potential Φ . This approximation is expected to be reliable for Φ < (0.1). From Eqs. (11.14) -(11.16) it follows that Substituting these relations into Eq. (11.13), the field equation inside the body is approximately given by If is close to at˜= 0, the field stays around in the region 0 <˜<˜1. The body has a thin-shell if˜1 is close to the radius˜of the body.
In the region 0 <˜<˜1 the field derivative of the effective potential around = can be approximated by d eff /d = , +˜≃ 2 ( − ). The solution to Eq. (11.19) can be obtained by writing the field as = 0 + , where 0 is the solution in the Minkowski background and is the perturbation induced by Φ . At linear order in and Φ we obtain In the region˜1 <˜<˜the field | (˜)| evolves towards larger values with increasing˜. Since the matter coupling term˜dominates over , in this regime, it follows that d eff /d ≃˜.
In the region outside the body (˜>˜) the field climbs up the potential hill after it acquires sufficient kinetic energy in the regime˜1 <˜<˜. Provided that the field kinetic energy dominates over its potential energy, the r.h.s. of Eq. (11.13) can be neglected relative to its l.h.s. of it. Moreover the terms that include˜and˜in the square bracket on the l.h.s. of Eq. (11.13) is much smaller than the term (1 + 2Φ )/˜. Using Eq. (11.17), it follows that 24) whose solution satisfying the boundary condition (˜→ ∞) = is where is a constant.
The coefficients , , , are known by matching the solutions (11.21), (11.23), (11.25) and their derivatives at˜=˜1 and˜=˜. If the body has a thin-shell, then the condition Δ˜= −˜1 ≪˜is satisfied. Under the linear expansion in terms of the three parameters Δ˜/˜, Φ , Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 and 1/(˜) we obtain the following field profile [594]: (11.28) where th = ( − )/(6 Φ ) is the thin-shell parameter, and . (11.31) As long as˜1Φ ≫ 1, the parameter is much smaller than 1. In order to derive the above field profile we have used the fact that the radius˜1 is determined by the condition 2 [ (˜1) − ] =˜, and hence (11.32) where is defined by (11.33) which is much smaller than 1. Using Eq. (11.32) we obtain the thin-shell parameter In terms of a linear expansion of , , Δ˜/˜, Φ , the field profile (11.28) outside the body is

35)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where the effective coupling is To leading-order this gives eff = 3 [Δ˜/˜+ 1/(˜)] = 3 th , which agrees with the result (5.45) in the Minkowski background. As long as Δ˜/˜≪ 1 and 1/(˜) ≪ 1, the effective coupling | eff | can be much smaller than the bare coupling | |, even in a strong gravitational background. From Eq. (11.26) the field value and its derivative around the center of the body with radius ≪ 1/ are given by . (11.38) In the Minkowski background (Φ = 0), Eq. (11.38) gives ′ (˜) > 0 for > 0 (or ′ (˜) < 0 for < 0). In the strong gravitational background (Φ ̸ = 0) the second term in the square bracket of Eq. (11.38) can lead to negative ′ (˜) for > 0 (or positive ′ (˜) for < 0), which leads to the evolution of | (˜)| toward 0. This effects comes from the presence of the strong relativistic pressure around the center of the body. Unless the boundary conditions at˜= 0 are appropriately chosen the field tends to evolve toward | (˜)| = 0, as seen in numerical simulations in [349,350] for the f (R) model (4.84). However there exists a thin-shell field profile even for ′ (˜) > 0 (and = −1/ √ 6) around the center of the body. In fact, the derivative ′ (˜) can change its sign in the regime 1/ <˜<˜1 for thin-shell solutions, so that the field does not reach the curvature singularity at = 0 [594].
For the inverse power-law potential ( ) = 4+ − , the existence of thin-shell solutions was numerically confirmed in [594] for Φ < 0.3. Note that the analytic field profile (11.26) was used to set boundary conditions around the center of the body. In Figure 10 we show the normalized field = / versus˜/˜for the model ( ) = 6 −2 with Φ = 0.2, Δ˜/˜= 0.1,˜= 20, and = 1. While we have neglected the term , relative to˜to estimate the solution in the region˜1 <˜<˜analytically, we find that this leads to some overestimation for the field value outside the body (˜>˜). In order to obtain a numerical field profile similar to the analytic one in the region˜>˜, we need to choose a field value slightly larger than the analytic value around the center of the body. The numerical simulation in Figure 10 corresponds to the choice of such a boundary condition, which explicitly shows the presence of thin-shell solutions even for a strong gravitational background.

Relativistic stars in metric f (R) gravity
The results presented above are valid for BD theory including metric f (R) gravity with the coupling = −1/ √ 6. While the analysis was carried out in the Einstein frame, thin-shell solutions were numerically found in the Jordan frame of metric f (R) gravity for the models (4.83) and (4.84) [43,600,42]. In these models the field = √︀ 3/2 ln in the region of high density ( ≫ ) is very close to the curvature singularity at = 0. Originally it was claimed in [266,349] that relativistic stars are absent because of the presence of this accessible singularity. However, as we have discussed in Section 11.2, the crucial point for obtaining thin-shell solutions is not the Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 existence of the curvature singularity but the choice of appropriate boundary conditions around the center of the star. For the correct choice of boundary conditions the field does not reach the singularity and thin-shell field profiles can be instead realized. In the Starobinsky's model (4.84), static configurations of a constant density star have been found for the gravitational potential Φ smaller than 0.345 [600].
For the star with an equation of state˜< 3˜, the effective potential of the field (in the presence of a matter coupling) does not have an extremum, see Eq. (11.7). In those cases the analytic results in Section 11.2 are no longer valid. For the equation of state˜< 3˜there is a tachyonic instability that tends to prevent the existence of a static star configuration [42]. For realistic neutron stars, however, the equation of state proposed in the literature satisfies the condition˜> 3˜throughout the star.
Babichev and Langlois [43,42] chose a polytropic equation of state for the energy density and the pressure in the Jordan frame:  [43,42] showed that 3˜can remain smaller thanf or realistic neutron stars. Note that the energy density is a decreasing function with respect to the distance from the center of star. Even for such a varying energy density, static star configurations have been shown to exist [43,42].
The ratio between the central density center and the cosmological density at infinity is parameterized by the quantity 0 = 2 pl / center . Realistic values of 0 are extremely small and it is a challenging to perform precise numerical simulations in such cases. We also note that the field mass Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 in the relativistic star is very much larger than its cosmological mass and hence a very high accuracy is required for solving the field equation numerically [600,581]. The authors in [43,600,42] carried out numerical simulations for the values of 0 of the order of 10 −3 -10 −4 . Figure 11 illustrates an example of the thin-shell field profile for the polytropic equation of state (11.39) in the model (4.84) with = 1 and 0 = 10 −4 [43]. In the regime 0 <˜< 1.5 the field is nearly frozen around the extremum of the effective potential, but it starts to evolve toward its asymptotic value = for˜> 1.5. Although the above analysis is based on the f (R) models (4.83) and (4.84) having a curvature singularity at = 0, such a singularity can be cured by adding the 2 term [350]. The presence of the 2 term has an advantage of realizing inflation in the early universe. However, the f (R) models (4.83) and (4.84) plus the 2 term cannot relate the epoch of two accelerations smoothly [37]. An example of viable models that can allow a smooth transition without a curvature singularity is [37] ( ) = (1 − ) + ln where , (0 < < 1/2), , and are constants. In [42] a static field profile was numerically obtained even for the model (11.40).
Although we have focused on the stellar configuration with Φ 0.3, there are also works of finding static or rotating black hole solutions in f (R) gravity [193,497]. Cruz-Dombriz et al. [193] derived static and spherically symmetric solutions by imposing that the curvature is constant. They also used a perturbative approach around the Einstein-Hilbert action and found that only solutions of the Schwarzschild-Anti de Sitter type are present up to second order in perturbations. The existence of general black hole solutions in f (R) gravity certainly deserves for further detailed study. It will be also of interest to study the transition from neutron stars to a strong-scalar-field Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 state in f (R) gravity [464]. While such an analysis was carried out for a massless field in scalartensor theory, we need to take into account the field mass in the region of high density for realistic models of f (R) gravity.
Pun et al. [498] studied physical properties of matter forming an accretion disk in the spherically symmetric metric in f (R) models and found that specific signatures of modified gravity can appear in the electromagnetic spectrum. In [92] the virial theorem for galaxy clustering in metric f (R) gravity was derived by using the collisionless Boltzmann equation. In [398] the construction of traversable wormhole geometries was discussed in metric f (R) gravity. It was found that the choice of specific shape functions and several equations of state gives rise to some exact solutions for f (R).

Gauss-Bonnet Gravity
So far we have studied modification to the Einstein-Hilbert action via the introduction of a general function of the Ricci scalar. Among the possible modifications of gravity this may be indeed a very special case. Indeed, one could think of a Lagrangian with all the infinite and possible scalars made out of the Riemann tensor and its derivatives. If one considers such a Lagrangian as a fundamental action for gravity, one usually encounters serious problems in the particle representations of such theories. It is well known that such a modification would introduce extra tensor degrees of freedom [635,283,284]. In fact, it is possible to show that these theories in general introduce other particles and that some of them may lead to problems.
For example, besides the graviton, another spin-2 particle typically appears, which however, has a kinetic term opposite in sign with respect to the standard one [572,67,302,465,153,303,99]. The graviton does interact with this new particle, and with all the other standard particles too. The presence of ghosts, implies the existence of particles propagating with negative energy. This, in turn, implies that out of the vacuum a particle (or more than one) and a ghost (or more than one) can appear at the same time without violating energy conservation. This sort of vacuum decay makes each single background unstable, unless one considers some explicit Lorentz-violating cutoff in order to set a typical energy/time scale at which this phenomenon occurs [145,161].
However, one can treat these higher-order gravity Lagrangians only as effective theories, and consider the free propagating mode only coming from the strongest contribution in the action, the Einstein-Hilbert one, for which all the modes are well behaved. The remaining higher-derivative parts of the Lagrangian can be regarded as corrections at energies below a certain fundamental scale. This scale is usually set to be equal to the Planck scale, but it can be lower, for example, in some models of extra dimensions. This scale cannot be nonetheless equal to the dark energy density today, as otherwise, one would need to consider all these corrections for energies above this scale. This means that one needs to re-sum all these contributions at all times before the dark energy dominance. Another possible approach to dealing with the ghost degrees of freedom consists of using the Euclidean-action path formalism, for which, one can indeed introduce a notion of probability amplitude for these spurious degrees of freedom [294,162].
The late-time modifications of gravity considered in this review correspond to those in low energy scales. Therefore we have a correction which begins to be important at very low energy scales compared to the Planck mass. In general this means that somehow these corrections cannot be treated any longer as corrections to the background, but they become the dominant contribution. In this case the theory cannot be treated as an effective one, but we need to assume that the form of the Lagrangian is exact, and the theory becomes a fundamental theory for gravity. In this sense these theories are similar to quintessence, that is, a minimally coupled scalar field with a suitable potential. The potential is usually chosen such that its energy scale matches with the dark energy density today. However, for this theory as well, one needs to consider this potential as fundamental, i.e., it does not get quantum corrections that can spoil the form of the potential itself. Still it may not be renormalizable, but so far we do not know any 4-dimensional renormalizable theory of gravity. In this case then, if we introduce a general modification of gravity responsible for the late-time cosmic acceleration, we should prevent this theory from introducing spurious ghost degrees of freedom.

Lovelock scalar invariants
One may wonder whether it is possible to remove these spin-2 ghosts. To answer this point, one should first introduce the Lovelock scalars [399]. These scalars are particular combinations/contractions of the Riemann tensor which have a fundamental property: if present in the Lagrangian, they only introduce second-order derivative contributions to the equations of motion.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Let us give an example of this property [399]. Soon after Einstein proposed General Relativity [226] and Hilbert found the Lagrangian to describe it [301], Kretschmann [372] pointed out that general covariance alone cannot explain the form of the Lagrangian for gravity. In the action he introduced, instead of the Ricci scalar, the scalar which now has been named after him, the Kretschmann scalar: At first glance this action looks well motivated. The Riemann tensor is a fundamental tensor for gravitation, and the scalar quantity 1 ≡ can be constructed by just squaring it. Furthermore, it is a theory for which Bianchi identities hold, as the equations of motion have both sides covariantly conserved. However, in the equations of motion, there are terms proportional to ∇ ∇ together with its symmetric partner ( ↔ ). This forces us to give in general at a particular slice of spacetime, together with the metric elements , their first, second, and third derivatives. Hence the theory has many more degrees of freedom with respect to GR.
In addition to the Kretschmann scalar there is another scalar 2 ≡ which is quadratic in the Riemann tensor . One can avoid the appearance of terms proportional to ∇ ∇ ( ) for the scalar quantity, which is called the Gauss-Bonnet (GB) term [572,67]. If one uses this invariant in the action of dimensions, as then the equations of motion coming from this Lagrangian include only the terms up to second derivatives of the metric. The difference between this scalar and the Einstein-Hilbert term is that this tensor is not linear in the second derivatives of the metric itself. It seems then an interesting theory to study in detail. Nonetheless, it is a topological property of four-dimensional manifolds that √ − can be expressed in terms of a total derivative [150], as Then the contribution to the equations of motion disappears for any boundaryless manifold in four dimensions.
In order to see the contribution of the GB term to the equations of motion one way is to couple it with a scalar field , i.e., ( ) , where ( ) is a function of . More explicitly the action of such theories is in general given by where ( ), ( ), and ( ) are functions of . The GB coupling of this form appears in the low energy effective action of string-theory [275,273], due to the presence of dilaton-graviton mixing terms.
There is another class of general GB theories with a self-coupling of the form [458],

7)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 where ( ) is a function in terms of the GB term (here we used the unit 2 = 1). The equations of motion, besides the standard GR contribution, will get contributions proportional to ∇ ∇ , [188,189]. This theory possesses more degrees of freedom than GR, but the extra information appears only in a scalar quantity , and its derivative. Hence it has less degrees of freedom compared to Kretschmann gravity, and in particular these extra degrees of freedom are not tensor-like. This property comes from the fact that the GB term is a Lovelock scalar. Theories with the more general Lagrangian density /2 + ( , 1 , 2 ) have been studied by many people in connection to the dark energy problem [142,110,521,420,585,64,166,543,180]. These theories are plagued by the appearance of spurious spin-2 ghosts, unless the Gauss-Bonnet (GB) combination is chosen as in the action (12.7) [465,153,447] (see also [110,181,109]).
Let us go back to discuss the Lovelock scalars. How many are they? The answer is infinite (each of them consists of linear combinations of equal powers of the Riemann tensor). However, because of topological reasons, the only non-zero Lovelock scalars in four dimensions are the Ricci scalar and the GB term . Therefore, for the same reasons as for the GB term, a general function of f (R) will only introduce terms in the equations of motion of the form ∇ ∇ , where ≡ / . Once more, the new extra degrees of freedom introduced into the theory comes from a scalar quantity, .
In summary, the Lovelock scalars in the Lagrangian prevent the equations of motion from getting extra tensor degrees of freedom. A more detailed analysis of perturbations on maximally symmetric spacetimes shows that, if non-Lovelock scalars are used in the action, then new extra tensor-like degrees of freedom begin to propagate [572,67,302,465,153,303,99]. Effectively these theories, such as Kretschmann gravity, introduce two gravitons, which have kinetic operators with opposite sign. Hence one of the two gravitons is a ghost. In order to get rid of this ghost we need to use the Lovelock scalars. Therefore, in four dimensions, one can in principle study the following action This theory will not introduce spin-2 ghosts. Even so, the scalar modes need to be considered more in detail: they may still become ghosts. Let us discuss more in detail what a ghost is and why we need to avoid it in a sensible theory of gravity.

Ghosts
What is a ghost for these theories? A ghost mode is a propagating degree of freedom with a kinetic term in the action with opposite sign. In order to see if a ghost is propagating on a given background, one needs to expand the action at second order around the background in terms of the perturbation fields. After integrating out all auxiliary fields, one is left with a minimal number of gauge-invariant fields ⃗ . These are not unique, as we can always perform a field redefinition (e.g., a field rotation). However, no matter which fields are used, we typically need -for non-singular Lagrangians -to define the kinetic operator, the operator which in the Lagrangian appears as ℒ =⃗⃗ + . . . [186,185]. Then the sign of the eigenvalues of the matrix defines whether a mode is a ghost or not. A negative eigenvalue would correspond to a ghost particle. On a FLRW background the matrix will be in general time-dependent and so does the sign of the eigenvalues. Therefore one should make sure that the extra scalar modes introduced for these theories do not possess wrong signs in the kinetic term at any time during the evolution of the Universe, at least up to today.
An overall sign in the Lagrangian does not affect the classical equations of motion. However, at the quantum level, if we want to preserve causality by keeping the optical theorem to be valid, then the ghost can be interpreted as a particle which propagates with negative energy, as already stated above. In other words, in special relativity, the ghost would have a four-momentum ( , ⃗ ) Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 with < 0. However it would still be a timelike particle as 2 − ⃗ 2 > 0, whether is negative or not. The problem arises when this particle is coupled to some other normal particle, because in this case the process 0 = + 1 + 2 + . . . with < 0 can be allowed. This means in general that for such a theory one would expect the pair creation of ghost and normal particles out of the vacuum. Notice that the energy is still conserved, but the energy is pumped out of the ghost particle.
Since all the particles are coupled at least to gravity, one would think that out of the vacuum particles could be created via the decay of a couple of gravitons emitted in the vacuum into ghosts and non-ghosts particles. This process does lead to an infinite contribution unless one introduces a cutoff for the theory [145,161], for which one can set observational constraints.
We have already seen that, for metric f (R) gravity, the kinetic operator in the FLRW background reduces to given in Eq. (7.60) with the perturbed action (7.80). Since the sign of is determined by , one needs to impose > 0 in order to avoid the propagation of a ghost mode.

( ) gravity
Let us consider the theory (12.7) in the presence of matter, i.e.
where we have recovered 2 . For the matter action we consider perfect fluids with an equation of state . The variation of the action (12.9) leads to the following field equations [178,383] where is the energy-momentum tensor of matter. If ∝ , then it is clear that the theory reduces to GR.

Cosmology at the background level and viable ( ) models
In the flat FLRW background the (00) component of Eq. (12.10) leads to 3 2 = , − − 24 3, + 2 ( + ) , (12.11) where and are the energy densities of non-relativistic matter and radiation, respectively. The cosmological dynamics in ( ) dark energy models have been discussed in [458,165,383,188,633,430]. We can realize the late-time cosmic acceleration by the existence of a de Sitter point satisfying the condition [458] 3 2 where 1 and 1 are the Hubble parameter and the GB term at the de Sitter point, respectively. From the stability of the de Sitter point we require the following condition [188] 0 < 6 1 , The GB term is given by where eff = −1 − 2˙/(3 2 ) is the effective equation of state. We have < 0 and˙> 0 during both radiation and matter domination. The GB term changes its sign from negative to positive during the transition from the matter era ( = −12 4 ) to the de Sitter epoch ( = 24 4 ). Perturbing Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Eq. (12.11) about the background radiation and matter dominated solutions, the perturbations in the Hubble parameter involve the mass squared given by 2 ≡ 1/(96 4 , ) [188]. For the stability of background solutions we require that 2 > 0, i.e., , > 0. Since the term 24 3i n Eq. (12.11) is of the order of 8 , , this is suppressed relative to 3 2 for 6 , ≪ 1 during the radiation and matter dominated epochs. In order to satisfy this condition we require that , approaches 0 in the limit | | → ∞. The deviation from the ΛCDM model can be quantified by the following quantity [182] ≡, =˙, = 72 6 where a prime represents a derivative with respect to = ln . During the radiation and matter eras we have = 192 6 , and = 72 6 , , respectively, whereas at the de Sitter attractor = 0.
The GB term inside and outside a spherically symmetric body (mass ⊙ and radius ⊙ ) with a homogeneous density is given by = −48( ⊙ ) 2 / 6 ⊙ and = 48( ⊙ ) 2 / 6 , respectively ( is a distance from the center of symmetry). In the vicinity of Sun or Earth, | | is much larger than the present cosmological GB term, 0 . As we move from the interior to the exterior of the star, the GB term crosses 0 from negative to positive. This means that ( ) and its derivatives with respect to need to be regular for both negative and positive values of whose amplitudes are much larger than 0 .
The above discussions show that viable ( ) models need to obey the following conditions: 1. ( ) and its derivatives , , , , . . . are regular.
A couple of representative models that can satisfy these conditions are [188] (A) , respectively. They are constructed to give rise to the positive , for all . Of course other models can be introduced by following the same prescription. These models can pass the constraint of successful expansion history that allows the smooth transition from radiation and matter eras to the accelerated epoch [188,633]. Although it is possible to have a viable expansion history at the background level, the study of matter density perturbations places tight constraints on these models. We shall address this issue in Section 12.3.4.

Numerical analysis
In order to discuss cosmological solutions in the low-redshift regime numerically for the models (12.16) and (12.17), it is convenient to introduce the following dimensionless quantities where a prime represents a derivative with respect to = ln . The quantities 6 , and ( , − )/ 2 can be expressed by and once the model is specified. Figure 12: The evolution of (multiplied by 10 4 ) and eff versus the redshift = 0/ − 1 for the model (12.16) with parameters = 100 and = 3 × 10 −4 . The initial conditions are chosen to be = −1.499985, = 20, and Ω = 0.99999. We do not take into account radiation in this simulation. From [182]. Figure 12 shows the evolution of and eff without radiation for the model (12.16) with parameters = 100 and = 3 × 10 −4 . The quantity is much smaller than unity in the deep matter era ( eff ≃ 0) and it reaches a maximum value prior to the accelerated epoch. This is followed by the decrease of toward 0, as the solution approaches the de Sitter attractor with eff = −1. While the maximum value of in this case is of the order of 10 −4 , it is also possible to realize larger maximum values of such as max 0.1.
For high redshifts the equations become too stiff to be integrated directly. This comes from the fact that, as we go back to the past, the quantity , (or ) becomes smaller and smaller. In fact, this also occurs for viable f (R) dark energy models in which , decreases rapidly for higher . Here we show an iterative method (known as the "fixed-point" method) [420,188] that can be used in these cases, provided no singularity is present in the high redshift regime [188]. We definē 2 and¯to be¯2 ≡ 2 / 2 0 and¯≡ / 4 0 , where the subscript "0" represents present values. The models (A) and (B) can be written in the form where¯2 Λ = Ω (0) / 3 +Ω (0) / 4 +Λ/3 (which represents the Hubble parameter in the ΛCDM model).
In the following we omit the tilde for simplicity. In Eq. (12.24) there are derivatives of in terms of up to second-order. Then we write Eq. (12.24) in the form At high redshifts ( 0.01) the models (A) and (B) are close to the ΛCDM model, i.e., 2 ≃ 2 Λ . As a starting guess we set the solution to be 2 (0) = 2 Λ . The first iteration is then 2 ′′ )︀ . [︁ = ln for the model (12.16) with = 0, 1, · · · , 6. The model parameters are = 10 and = 0.075. The iterative method provides the solutions with high accuracy in the regime −4. From [188].
If the starting guess is in the basin of a fixed point, 2 ( ) will converge to the solution of the equation after the -th iteration. For the convergence we need the following condition which means that each correction decreases for larger . The following relation is also required to be satisfied: Once the solution begins to converge, one can stop the iteration up to the required/available level of precision. In Figure 13 we plot absolute errors for the model (12.16), which shows that the iterative method can produce solutions accurately in the high-redshift regime. Typically this method stops working when the initial guess is outside the basin of convergence. This happen for low redshifts in which the modifications of gravity come into play. In this regime we just need to integrate Eqs. (12.19) -(12.22) directly.

Solar system constraints
We study local gravity constraints on cosmologically viable ( ) models. First of all there is a big difference between ( ) and f (R) theories. The vacuum GR solution of a spherically symmetric manifold, the Schwarzschild metric, corresponds to a vanishing Ricci scalar ( = 0) outside the star. In the presence of non-relativistic matter, approximately equals to the matter density 2 for viable f (R) models. On the other hand, even for the vacuum exterior of the Schwarzschild metric, the GB term has a non-vanishing value = = 12 2 / 6 [178,185], where = 2 ⊙ / ⊙ is the Schwarzschild radius of the object. In the regime | | ≫ * the models (A) and (B) have a correction term of the order √ * 2 * / 2 plus a cosmological constant term −( + 1) √ * . Since does not vanish even in the vacuum, the correction term 2 * / 2 can be much smaller than 1 even in the absence of non-relativistic matter. If matter is present, this gives rise to the contribution of the order of 2 ≈ ( 2 ) 2 to the GB term. The ratio of the matter contribution to the vacuum value (0) = 12 2 / 6 is estimated as At the surface of Sun (radius ⊙ = 6.96 × 10 10 cm = 3.53 × 10 24 GeV −1 and mass ⊙ = 1.99 × 10 33 g = 1.12 × 10 57 GeV), the density drops down rapidly from the order ≈ 10 −2 g/cm 3 to the order ≈ 10 −16 g/cm 3 . If we take the value = 10 −2 g/cm 3 we have ≈ 4 × 10 −5 (where we have used 1 g/cm 3 = 4.31 × 10 −18 GeV 4 ). Taking the value = 10 −16 g/cm 3 leads to a much smaller ratio: ≈ 4 × 10 −33 . The matter density approaches a constant value ≈ 10 −24 g/cm 3 around the distance = 10 3 ⊙ from the center of Sun. Even at this distance we have ≈ 4 × 10 −31 , which means that the matter contribution to the GB term can be neglected in the solar system we are interested in.
In order to discuss the effect of the correction term 2 * / 2 on the Schwarzschild metric, we introduce a dimensionless parameter = √︀ * / , (12.29) where = 12/ 4 is the scale of the GB term in the solar system. Since √ * is of the order of the Hubble parameter 0 ≈ 70 km sec −1 Mpc −1 , the parameter for the Sun is approximately given by ≈ 10 −46 . We can then decompose the vacuum equations in the form + Σ = 0 , (12.30) where is the Einstein tensor and Here˜is defined by =˜. We introduce the following ansatz for the metric where the functions and are expanded as power series in , as = 0 ( ) + 1 ( ) + ( 2 ) , = 0 ( ) + 1 ( ) + ( 2 ) .
Since 0 and 0 are known, one can solve the differential equations for 1 and 1 . This process can be iterated order by order in . For the model (A) introduced in (12.16), we obtain the following differential equations for 1 and 1 [185]: where ≡ / . The solutions to these equations are Here we have neglected the contribution coming from the homogeneous solution, as this would correspond to an order renormalization contribution to the mass of the system. Although ≪ 1, the term in ln only contributes by a factor of order 10 2 . Since ≫ 1 the largest contributions to 1 and 1 correspond to those proportional to 3 , which are different from the Schwarzschildde Sitter contribution (which grows as 2 ). Hence the model (12.16) gives rise to the corrections larger than that in the cosmological constant case by a factor of . Since is very small, the contributions to the solar-system experiments due to this modification are too weak to be detected. The strongest bound comes from the shift of the perihelion of Mercury, which gives the very mild bound < 2 × 10 5 [185]. For the model (12.17) the constraint is even weaker, (1 + ) < 10 14 . In other words, the corrections look similar to the Schwarzschild-de Sitter metric on which only very weak bounds can be placed.

Ghost conditions in the FLRW background
In the following we shall discuss ghost conditions for the action (12.9). For simplicity let us consider the vacuum case ( = 0) in the FLRW background. The action (12.9) can be expanded at second order in perturbations for the perturbed metric (6.1), as we have done for the action (6.2) in Section 7.4. Before doing so, we introduce the gauge-invariant perturbed quantity This quantity completely describes the dynamics of all the scalar perturbations. Note that for the gauge choice = 0 one has ℛ = . Integrating out all the auxiliary fields, we obtain the second-order perturbed action [186] ( Recall that has been introduced in Eq. (12.15).
In order to avoid that the scalar mode becomes a ghost, one requires that > 0, i.e.
> −1 /4 . (12.44) This relation is dynamical, as one requires to know how and its derivatives change in time. Therefore whatever ( ) is, the propagating scalar mode can still become a ghost. If, > 0 and > 0, then > 0 and hence the ghost does not appear. The quantity characterizes the speed of propagation for the scalar mode, which is again dependent on the dynamics. For any GB theory, one can give initial conditions of and˙such that 2 becomes negative. This instability, if present, governs the high momentum modes in Fourier space, which corresponds to an Ultra-Violet (UV) instability. In order to avoid this UV instability in the vacuum, we require that the effective equation of state satisfies eff < −2/3. At the de Sitter point ( eff = −1) the speed is time-independent and reduces to the speed of light ( = 1). Suppose that the scalar mode does not have a ghost mode, i.e., > 0. Making the field redefinition = ℛ and = √ , the action (12.41) can be written as where a prime represents a derivative with respect to = ∫︀ −1 d and 2 ≡ − ′′ /( 2 ). In order to realize the positive mass squared ( 2 > 0), the condition , > 0 needs to be satisfied in the regime ≪ 1 (analogous to the condition , > 0 in metric f (R) gravity).

Viability of ( ) gravity in the presence of matter
In the presence of matter, other degrees of freedom appear in the action. Let us take into account a perfect fluid with the barotropic equation of state = / . It can be proved that, for small scales (i.e., for large momenta ) in Fourier space, there are two different propagation speeds given by [182]  The first result is expected, as it corresponds to the matter propagation speed. Meanwhile the presence of matter gives rise to a correction term to 2 2 in Eq. (12.43). This latter result is due to the fact that the background equations of motion are different between the two cases. Recall that for viable ( ) models one has | | ≪ 1 at high redshifts. Since the background evolution is approximately given by Hence the UV instability can be avoided for < −1/2. During the radiation era ( = 1/3) and the matter era ( = 0), the large momentum modes are unstable. In particular this leads to the violent growth of matter density perturbations incompatible with the observations of large-scale structure [383,182]. The onset of the negative instability can be characterized by the condition [182] ≈ ( / ) 2 . (12.49) As long as ̸ = 0 we can always find a wavenumber (≫ ) satisfying the condition (12.49). For those scales linear perturbation theory breaks down, and in principle one should look for all Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 higher-order contributions. Hence the background solutions cannot be trusted any longer, at least for small scales, which makes the theory unpredictable. In the same regime, one can easily see that the scalar mode is not a ghost, as Eq. (12.44) is satisfied (see Figure 12). Therefore the instability is purely classical. This kind of UV instability sets serious problems for any theory, including ( ) gravity.

The speed of propagation in more general modifications of gravity
We shall also discuss more general theories given by Eq. (12.8) where we do not take into account the matter term here. It is clear that this function allows more freedom with respect to the background cosmological evolution 8 , as now one needs a two-parameter function to choose. However, once more the behavior of perturbations proves to be a strong tool in order to have a deep insight into the theory. The second-order action for perturbations is given by where we have introduced the gauge-invariant field with ≡ , and ≡ , . The forms of ( ), 1 ( ) and 2 ( ) are given explicitly in [186]. The quantity 2 vanishes either on the de Sitter solution or for those theories satisfying Δ ≡ This result implies that the superluminal propagation is always present in these theories, and the speed is scale-dependent. On the other hand, when Δ = 0, this behavior is not present at all. Therefore, there is a physical property by which different modifications of gravity can be distinguished. The presence of an extra matter scalar field does not change this regime at high [185], because the Laplacian of the gravitational field is not modified by the field coupled to gravity in the form ( , , ).

Gauss-Bonnet gravity coupled to a scalar field
At the end of this section we shall briefly discuss theories with a GB term coupled to a scalar field with the action given in Eq. (12.6). The scalar coupling with the GB term often appears as higher-order corrections to low-energy, tree-level effective string theory based on toroidal compactifications [275,276]. More explicitly the low-energy string effective action in four dimensions is given by where is a dilaton field that controls the string coupling parameter, 2 = . The above action is the string frame action in which the dilaton is directly coupled to a scalar curvature, . The Lagrangian ℒ is that of additional matter fields (fluids, axion, modulus etc.). The Lagrangian ℒ corresponds to higher-order string corrections including the coupling between the GB term and the dilaton. A possible set of corrections include terms of the form [273,105,147] where ′ is a string expansion parameter and ( ) is a general function of . The constant is an additional parameter which depends on the types of string theories: = −1/4, −1/8, and 0 correspond to bosonic, heterotic, and superstrings, respectively. If we require that the full action agrees with the three-graviton scattering amplitude, the coefficients and are fixed to be = −1, = 1, and ( ) = − − [425].
In the Pre-Big-Bang (PBB) scenario [275] the dilaton evolves from a weakly coupled regime ( ≪ 1) toward a strongly coupled region ( 1) during which the Hubble parameter grows in the string frame (superinflation). This superinflation is driven by a kinetic energy of the dilaton field and it is called a PBB branch. There exists another Friedmann branch with a decreasing curvature. If ℒ = 0 these branches are disconnected to each other with the appearance of a curvature singularity. However the presence of the correction ℒ allows the existence of nonsingular solutions that connect two branches [273,105,147].
The corrections ℒ are the sum of the tree-level ′ corrections and the quantum -loop corrections ( = 1, 2, 3, · · · ) with the function ( ) given by ( ) = − ∑︀ =0 ( −1) , where ( ≥ 1) are coefficients of -loop corrections (with 0 = 1). In the context of the PBB cosmology it was shown in [105] there exist regular cosmological solutions in the presence of tree-level and one-loop corrections, but this is not realistic in that the Hubble rate in Einstein frame continues to increase after the bounce. Nonsingular solutions that connect to a Friedmann branch can be obtained by accounting for the corrections up to two-loop with a negative coefficient ( 2 < 0) [105,147]. In the context of Ekpyrotic cosmology where a negative potential ( ) is present in the Einstein frame, it is possible to realize nonsingular solutions by taking into account corrections similar to ℒ given above [588]. For a system in which a modulus field is coupled to the GB term, one can also realize regular solutions even without the higher-derivative term (∇ ) 4 in Eq. (12.57) [34,224,336,337,338,623,12,582]. These results show that the GB term can play a crucial role to eliminate the curvature singularity.
In the context of dark energy there are some works which studied the effect of the GB term on the late-time cosmic acceleration. A simple model that can give rise to cosmic acceleration is provided by the action [463] = where ( ) and ( ) are functions of a scalar field . This can be viewed as the action in the Einstein frame corresponding to the Jordan frame action (12.56). We note that the conformal Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 transformation gives rise to a coupling between the field and non-relativistic matter in the Einstein frame. Such a coupling is assumed to be negligibly small at low energy scales, as in the case of the runaway dilaton scenario [274,176]. For the exponential potential ( ) = 0 − and the coupling ( ) = ( 0 / ) , cosmological dynamics has been extensively studied in [463,360,361,593] (see also [523,452,453,381]). In particular it was found in [360,593] that a scaling matter era can be followed by a late-time de Sitter solution which appears due to the presence of the GB term.
Koivisto and Mota [360] placed observational constraints on the above model using the Gold data set of Supernovae Ia together with the CMB shift parameter data of WMAP. The parameter is constrained to be 3.5 < < 4.5 at the 95% confidence level. In the second paper [361], they included the constraints coming from the BBN, LSS, BAO and solar system data and showed that these data strongly disfavor the GB model discussed above. Moreover, it was shown in [593] that tensor perturbations are subject to negative instabilities in the above model when the GB term dominates the dynamics (see also [290]). Amendola et al. [25] studied local gravity constraints on the model (12.58) and showed that the energy contribution coming from the GB term needs to be strongly suppressed for consistency with solar-system experiments. This is typically of the order of Ω GB 10 −30 and hence the GB term of the coupling ( ) cannot be responsible for the current accelerated expansion of the universe.
In summary the GB gravity with a scalar field coupling allows nonsingular solutions in the high curvature regime, but such a coupling is difficult to be compatible with the cosmic acceleration at low energy scales. Recall that dark energy models based on ( ) gravity also suffers from the UV instability problem. This shows how the presence of the GB term makes it difficult to satisfy all experimental and observational constraints if such a modification is responsible for the latetime acceleration. This property is different from metric f (R) gravity in which viable dark energy models can be constructed.

Other Aspects of f (R) Theories and Modified Gravity
In this section we discuss a number of topics related with f (R) theories and modified gravity. These include weak lensing, thermodynamics and horizon entropy, unified models of inflation and dark energy, f (R) theories in the extra dimensions, Vainshtein mechanism, DGP model, Noether and Galileon symmetries.

Weak lensing
Weak gravitational lensing is sensitive to the growth of large scale structure as well as the relation between matter and gravitational potentials. Since the evolution of matter perturbations and gravitational potentials is different from that of GR, the observations of weak lensing can provide us an important test for probing modified gravity on galactic scales (see [2,527,27,595,528,548] for theoretical aspects and [546,348,322,3,629,373,326,177,73] for observational aspects). In particular a number of wide-field galaxy surveys are planned to measure galaxy counts and weak lensing shear with high accuracy, so these will be useful to distinguish between modified gravity and the ΛCDM model in future observations. Let us consider BD theory with the action (10.10), which includes f (R) gravity as a specific case. Note that the method explained below can be applied to other modified gravity models as well. The equations of matter perturbations and gravitational potentials Φ, Ψ in BD theory have been already derived under the quasi-static approximation on sub-horizon scales ( ≫ ), see Eqs. (10.38), (10.39), and (10.40). In order to discuss weak lensing observables, we define the lensing deflecting potential Φ wl ≡ Φ + Ψ , (13.1) and the effective density field Writing the angular position of a source and the direction of weak lensing observation to be ⃗ and ⃗ , respectively, the deformation of the shape of galaxies can be quantified by the amplification matrix = d ⃗ /d ⃗ . The components of the matrix are given by [66] = − where is the maximum distance to the source and ( ) ≡ ∫︀ The convergence is a function on the 2-sphere and hence it can be expanded in the form where ⃗ ℓ = (ℓ 1 , ℓ 2 ) with ℓ 1 and ℓ 2 integers. We define the power spectrum Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 of the shear to be ⟨^w l ( ⃗ ℓ)^* wl ( ⃗ ℓ ′ )⟩ = (ℓ) (2) ( ⃗ ℓ − ⃗ ℓ ′ ). Then the convergence has the same power spectrum as , which is given by [66,601] (ℓ) = We assume that the sources are located at the distance (corresponding to the redshift ), giving the relations ( ) = ( − ) and ( ) = ( − )/ . From Eq. (13.3) eff can be expressed as where is the matter power spectrum. Hence the convergence spectrum (13.6) reads We recall that, during the matter era, the transition from the GR regime ( ∝ 2/3 and Φ wl = constant) to the scalar-tensor regime ( ∝ ( √ 25+48 2 −1)/6 and Φ wl ∝ ( √ 25+48 2 −5)/6 ) occurs at the redshift characterized by the condition (10.45). Since the early evolution of perturbations is similar to that in the ΛCDM model, the weak lensing potential at late times is given by the formula [214] Φ wl ( , ) = 9 The initial power spectrum generated during inflation is where Φ is the spectral index and 2 is the amplitude of Φ wl [71,214]. Therefore we obtain the power spectrum of matter perturbations, as Plugging Eq. (13.11) into Eq. (13.7), we find that the convergence spectrum is given by  Note that satisfies the differential equation d /d = 1/ ( ). In Figure 14 we plot the convergence spectrum in f (R) gravity with the potential (10.23) for two different values of together with the ΛCDM spectrum. Recall that this model corresponds to the f (R) model ( ) = − [︀ 1 − ( / ) −2 ]︀ with the correspondence = 2 /(2 + 1). Figure 14 shows the convergence spectrum in the linear regime characterized by ℓ 200. The ΛCDM model corresponds to the limit → ∞, i.e., → 1. The deviation from the ΛCDM model becomes more significant for smaller away from 1. Since the evolution of Φ wl changes from Φ wl = constant to Φ wl ∝ ( √ 25+48 2 −5)/6 at the transition time ℓ characterized by the condition 2 / = (ℓ/ ) 2 / 2 , this leads to a difference of the spectral index of the convergence spectrum compared to that of the ΛCDM model [595]: This estimation is reliable for the transition redshift ℓ much larger than 1. In the simulation of Figure 14 the numerical value of Δ for = 0.7 at ℓ = 200 is 0.056 (with ℓ = 3.26), which is slightly smaller than the analytic value Δ = 0.068 estimated by Eq. (13.14). The deviation of the spectral index of from the ΛCDM model will be useful to probe modified gravity in future high-precision observations. Note that the galaxy-shear correlation spectrum will be also useful to constrain modified gravity models [528].
Recent data analysis of the weak lensing shear field from the Hubble Space Telescope's COSMOS survey along with the ISW effect of CMB and the cross-correlation between the ISW and galaxy Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 distributions from 2MASS and SDSS surveys shows that the anisotropic parameter ≡ Ψ/Φ is constrained to be < 1 at the 98% confidence level [73]. For BD theory with the action (10.10) the quasi-static results (10.38) and (10.39) of the gravitational potentials give Since ≃ (1 − 2 2 )/(1 + 2 2 ) in the scalar-tensor regime ( 2 / 2 ≫ 2 / ), one can realize < 1 in BD theory. Of course we need to wait for further observational data to reach the conclusion that modified gravity is favored over the ΛCDM model.
To conclude this session we would like to point out the possibility of using the method of gravitational lensing tomography [574]. This method consists of considering lensing on different redshift data-bins. In order to use this method, one needs to know the evolution of both the linear growth rate and the non-linear one (typically found through a standard linear-to-non-linear mapping). Afterward, from observational data, one can separate different bins in order to make fits to the models. Having good data sets, this procedure is strong enough to further constrain the models, especially together with other probes such as CMB [322,320,632,292].

Thermodynamics and horizon entropy
It is known that in Einstein gravity the gravitational entropy of stationary black holes is proportional to the horizon area , such that = /(4 ), where is gravitational constant [75]. A black hole with mass obeys the first law of thermodynamics, d = d [59], where = /(2 ) is a Hawking temperature determined by the surface gravity [293]. This shows a deep physical connection between gravity and thermodynamics. In fact, Jacobson [324] showed that Einstein equations can be derived by using the Clausius relation d = d on local horizons in the Rindler spacetime together with the relation ∝ , where d and are the energy flux across the horizon and the Unruh temperature seen by an accelerating observer just inside the horizon respectively.
Unlike stationary black holes the expanding universe with a cosmic curvature has a dynamically changing apparent horizon with the radius¯= ( 2 + / 2 ) −1/2 , where is a cosmic curvature [108] (see also [296]). Even in the FLRW spacetime, however, the Friedmann equation can be written in the thermodynamical form d = −d + d , where is the work density present in the dynamical background [8]. For matter contents of the universe with energy density and pressure , the work density is given by = ( − )/2 [297,298]. This method is identical to the one established by Jacobson [324], that is, d = −d + d .
In metric f (R) gravity Eling et al. [228] showed that a non-equilibrium treatment is required such that the Clausius relation is modified to d = d / + d , where = /(4 ) is the Wald horizon entropy [610] and d is the bulk viscosity entropy production term. Note that corresponds to a Noether charge entropy. Motivated by this work, the connections between thermodynamics and modified gravity have been extensively discussed -including metric f (R) gravity [6,7,281,431,619,620,230,103,51,50,157] and scalar-tensor theory [281,619,620,108].
Let us discuss the relation between thermodynamics and modified gravity for the following general action [53] = ∫︁ where ≡ − (1/2) ∇ ∇ is a kinetic term of a scalar field . For the matter Lagrangian ℒ we take into account perfect fluids (radiation and non-relativistic matter) with energy density and pressure . In the FLRW background with the metric d 2 = ℎ d d +¯2dΩ 2 , wherē = ( ) and 0 = , 1 = with the two dimensional metric ℎ = diag(−1, Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 the Friedmann equations are given by where ≡ / and Note that and originate from the energy-momentum tensor ( ) defined by where the Einstein equation is given by Defining the density and the pressure of "dark" components in this way, they obey the following equation˙+ 3 ( + ) = 3 8 ( 2 + / 2 )˙. (13.24) For the theories with˙̸ = 0 (including f (R) gravity and scalar-tensor theory) the standard continuity equation does not hold because of the presence of the last term in Eq. (13.24). In the following we discuss the thermodynamical property of the theories given above. The apparent horizon is determined by the condition ℎ¯¯= 0, which gives¯= (︀ 2 + / 2 )︀ −1/2 in the FLRW spacetime. Taking the differentiation of this relation with respect to and using Eq. ( The apparent horizon has the Hawking temperature = | |/(2 ), where is the surface gravity given by

28)
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Here we have defined ≡ + and ≡ + . For the total equation of state = / less than 1/3, as is the case for standard cosmology, one has ≤ 0 so that the horizon temperature is given by The modified first-law of thermodynamics (13.35) suggests a deep connection between the horizon thermodynamics and Friedmann equations in modified gravity. The term d can be interpreted as a term of entropy production in the non-equilibrium thermodynamics [228]. The theories with = constant lead to d^= 0, which means that the first-law of equilibrium thermodynamics holds. The theories with d ̸ = 0, including f (R) gravity and scalar-tensor theory, give rise to the additional non-equilibrium term (13.36) [6,7,281,619,620,108,50,53].
The main reason why the non-equilibrium term d appears is that the energy density and the pressure defined in Eqs. (13.20) and (13.21) do not satisfy the standard continuity equation for˙̸ = 0. On the other hand, if we define the effective energy-momentum tensor ( ) as Eq. (2.9) in Section 2, it satisfies the continuity equation (2.10). This correspond to rewriting the Einstein equation in the form (2.8) instead of (13.23). Using this property, [53] showed that equilibrium description of thermodynamics can be possible by introducing the Bekenstein-Hawking entropŷ = /(4 ). In this case the horizon entropy^takes into account the contribution of both the Wald entropy in the non-equilibrium thermodynamics and the entropy production term.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 13.3 Curing the curvature singularity in f (R) dark energy models, unified models of inflation and dark energy In Sections 5.2 and 8.1 we showed that there is a curvature singularity for viable f (R) models such as (4.83) and (4.84). More precisely this singularity appears for the models having the asymptotic behavior (5.19) in the region of high density ( ≫ ). As we see in Figure 3, the field potential ( ) = ( − )/(2 2 ) has a finite value /(2 2 ) in the limit = √︀ 3/(16 ) pl ln → 0. In this limit one has , → 0, so that the scalaron mass 1/(3 , ) goes to infinity. This problem of the past singularity can be cured by adding the term 2 /(6 2 ) to the Lagrangian in f (R) dark energy models [37]. Let us then consider the modified version of the model (4.83): For this model one can easily show that the potential ( ) = ( − )/(2 2 ) extends to the region > 0 and that the curvature singularity disappears accordingly. Also the scalaron mass approaches the finite value in the limit → ∞. The perturbation is bounded from above, which can evade the problem of the dominance of the oscillation mode in the past.
Since the presence of the term 2 /(6 2 ) can drive inflation in the early universe, one may anticipate that both inflation and the late-time acceleration can be realized for the model of the type (13.37). This is like a modified gravity version of quintessential inflation based on a single scalar field [486,183,187,392]. However, we have to caution that the transition between two accelerated epochs needs to occur smoothly for successful cosmology. In other words, after inflation, we require a mechanism in which the universe is reheated and then the radiation/matter dominated epochs follow. However, for the model (13.37), the Ricci scalar evolves to the point , = 0 and it enters the region , < 0. Crossing the point , = 0 implies the divergence of the scalaron mass. Moreover, in the region , < 0, the Minkowski space is not a stable vacuum state. This is problematic for the particle creation from the vacuum during reheating. The similar problem arises for the models (4.84) and (4.89) in addition to the model proposed by Appleby and Battye [35]. Thus unified f (R) models of inflation and dark energy cannot be constructed easily in general (unlike a number of related works [456,460,462]). Brookfield et al. [104] studied the viability of the model ( ) = − / + ( , > 0) by using the constraints coming from Big Bang Nucleosynthesis and fifth-force experiments and showed that it is difficult to find a unique parameter range for consistency of this model.
In order to cure the above mentioned problem, Appleby et al. [37] proposed the f (R) model (11.40). Note that the case = 0 corresponds to the Starobinsky inflationary model ( ) = + 2 /(6 2 ) [564] and the case = 1/2 corresponds to the model of Appleby and Battye [35] plus the 2 /(6 2 ) term. Although the above mentioned problem can be evaded in this model, the reheating proceeds in a different way compared to that in the model ( ) = + 2 /(6 2 ) [which we discussed in Section 3.3]. Since the Hubble parameter periodically evolves between = 1/(2 ) and = / , the reheating mechanism does not occur very efficiently [37]. The reheating temperature can be significantly lower than that in the model ( ) = + 2 /(6 2 ). It will be of interest to study observational signatures in such unified models of inflation and dark energy.

f (R) theories in extra dimensions
Although f (R) theories have been introduced mainly in four dimensions, the same models may appear in the context of braneworld [502,501] in which our universe is described by a brane embedded in extra dimensions (see [404] for a review). This scenario implies a careful use of f (R) theories, because a boundary (brane) appears. Before looking at the real working scenario Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 in braneworld, it is necessary to focus on the mathematical description of f (R) models through a sensible definition of boundary conditions for the metric elements on the surface of the brane.
Some works appeared regarding the possibility of introducing f (R) theories in the context of braneworld scenarios [499,40,96,513,97]. In doing so one requires a surface term [222,482,69,48,49,286], which is known as the Hawking-Luttrell term [295] (analogous to the York-Gibbons-Hawking one for General Relativity). The action we consider is given by where ≡ / , is the determinant of the induced metric on the − 1 dimensional boundary, and is the trace of the extrinsic curvature tensor. In this case particular attention should be paid to boundary conditions on the brane, that is, the Israel junction conditions [323]. In order to have a well-defined geometry in five dimensions, we require that the metric is continuous across the brane located at = 0. However its derivatives with respect to can be discontinuous at = 0. The Ricci tensor in Eq. (2.4) is made of the metric up to the second derivatives ′′ with respect to . This means that ′′ have a deltafunction dependence proportional to the energy-momentum tensor at a distributional source (i.e., with a Dirac's delta function centered on the brane) [87,86,536]. In general this also leads to the discontinuity of the Ricci scalar across the brane. Since the discontinuity of can lead to inconsistencies in f (R) gravity, one should add this extra-constraint as a junction condition. In other words, one needs to impose that, although the metric derivative is discontinuous, the Ricci scalar should still remain continuous on the brane. This is tantamount to imposing that the extra scalar degree of freedom introduced is continuous on the brane. We use Gaussian normal coordinates with the metric This has a delta function behavior for the -components, leading to [207] ≡ where is the matter stress-energy tensor on the brane. Hence is continuous, whereas its first derivative is not, in general. This imposes an extra condition on the metric crossing the brane at = 0, compared to General Relativity in which the condition for the continuity of is not present. However, it is not easy to find a solution for which the metric derivative is discontinuous but is not. Therefore some authors considered matter on the brane which is not universally coupled with the induced metric. This approach leads to the relaxation of the condition that is continuous. Such a matter action can be found by analyzing the action in the Einstein frame and introducing a scalar field coupled to the scalaron on the brane as follows [207] The presence of the coupling ( ) with the field modifies the Israel junction conditions. Indeed, if = 0, then must be continuous, but if ̸ = 0, can have a delta function profile. This method may help for finding a solution for the bulk that satisfies boundary conditions on the brane.

Vainshtein mechanism
Modifications of gravity in recent works have been introduced mostly in order to explain the late-time cosmic acceleration. This corresponds to the large-distance modification of gravity, but gravity at small distances is subject to change as well. Modified gravity models of dark energy must pass local gravity tests in the solar system. The f (R) models discussed in Section 4 are designed to satisfy local gravity constraints by having a large scalar-field mass, while at the same time they are responsible for dark energy with a small mass compatible with the Hubble parameter today.
It is interesting to see which modified gravity theories have successful Newton limits. There are two known mechanisms for satisfying local gravity constraints, (i) the Vainshtein mechanism [602], and (ii) the chameleon mechanism [344,343] (already discussed in Section 5.2). Both consist of using non-linearities in order to prevent any other fifth force from propagating freely. The chameleon mechanism uses the non-linearities coming from matter couplings, whereas the Vainshtein mechanism uses the self-coupling of a scalar-field degree of freedom as a source for the non-linear effect.
There are several examples where the Vainshtein mechanism plays an important role. One is the massive gravity in which a consistent free massive graviton is uniquely defined by Pauli-Fierz theory [258,259]. The massive gravity described by the Fierz-Pauli action cannot be studied through the linearization close to a point-like mass source, because of the crossing of the Vainshtein radius, that is the distance under which the linearization fails to study the metric properly [602]. Then the theory is in the strong-coupling regime, and things become obscure as the theory cannot be understood well mathematically. A similar behavior also appears for the Dvali-Gabadadze-Porrati (DGP) model (we will discuss in the next section), in which the Vainshtein mechanism plays a key role for the small-scale behavior of this model.
Besides a standard massive term, other possible operators which could give rise to the Vainshtein mechanism come from non-linear self-interactions in the kinetic term of a matter field . One of such terms is given by ∇ ∇ , (13.43) which respects the Galilean invariance under which 's gradient shifts by a constant [455] (treated in section 13.7.2). This allows a robust implementation of the Vainshtein mechanism in that nonlinear self-interacting term can allow the decoupling of the field from matter in the gravitationally bounded system without introducing ghosts. Another example of the Vainshtein mechanism may be seen in ( ) gravity. Recall that in this theory the contribution to the GB term from matter can be neglected relative to the vacuum value (0) = 12 (2 ) 2 / 6 . In Section 12.3.3 we showed that on the Schwarzschild geometry the modification of gravity is very small for the models (12.16) and (12.17), because the GB term has a value much larger than its cosmological value today. The scalar-field degree of freedom acquires a large mass in the region of high density, so that it does not propagate freely. For the model (12.16) we already showed that at the linear level the coefficients and of the spherically symmetric metric (12.32) are where ≡ /(2 ), 1 ( ) and 1 ( ) are given by Eqs. (12.38) and (12.39), and ≈ 10 −46 for our solar system. Of course this result is trustable only in the region for which 1 ≪ 1/ . Outside this region non-linearities are important and one cannot rely on approximate methods any longer. Therefore, for this model, we can define the Vainshtein radius as For ∼ 1, this value is well outside the region in which solar-system experiments are carried out. This example shows that the Vainshtein radius is generally model-dependent.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 In metric f (R) gravity a non-linear effect coming from the coupling to matter fields (in the Einstein frame) is crucially important, because vanishes in the vacuum Schwarzschild background. The local gravity constraints can be satisfied under the chameleon mechanism rather than the non-linear self coupling of the Vainshtein mechanism.

DGP model
The Dvali-Gabadadze-Porrati (DGP) [220] braneworld model has been considered as a model which could modify gravity because of the existence of the extra-dimensions. In the DGP model the 3-brane is embedded in a Minkowski bulk spacetime with infinitely large 5th extra dimensions. The Newton's law can be recovered by adding a 4-dimensional (4D) Einstein-Hilbert action sourced by the brane curvature to the 5D action [219]. While the DGP model recovers the standard 4D gravity for small distances, the effect from the 5D gravity manifests itself for large distances.
Remarkably it is possible to realize the late-time cosmic acceleration without introducing an exotic matter source [201,203].
The DGP model is given by the action where˜is the metric in the 5D bulk and =˜is the induced metric on the brane with ( ) being the coordinates of an event on the brane labeled by . The first and second terms in Eq. (13.46) correspond to Einstein-Hilbert actions in the 5D bulk and on the brane, respectively. Note that 2 (5) and 2 (4) are 5D and 4D gravitational constants, respectively, which are related with 5D and 4D Planck masses, (5) and (4) , via 2 (5) = 1/ 3 (5) and 2 (4) = 1/ 2 (4) . The Lagrangian ℒ brane describes matter localized on the 3-brane.
The equations of motion read (5)  , (13.48) where is the extrinsic curvature [609] calculated on the brane and is the energy-momentum tensor of localized matter. Since ∇ ( − ) = 0, the continuity equation ∇ = 0 follows from Eq. (13.48). The length scale is defined by .

(13.49)
If we consider the flat FLRW brane ( = 0), we obtain the modified Friedmann equation [201,203]  In the DGP model the modification of gravity comes from a scalar-field degree of freedom, usually called , which is identified as a brane bending mode in the bulk. Then one may wonder if such a field mediates a fifth force incompatible with local gravity constraints. However, this is not the case, as the Vainshtein mechanism is at work in the DGP model for the length scale smaller than the Vainshtein radius * = ( 2 ) 1/3 , where is the Schwarzschild radius of a source. The model can evade solar system constraints at least under some range of conditions on the energymomentum tensor [204,285,496]. The Vainshtein mechanism in the DGP model originates from a non-linear self-interaction of the scalar-field degree of freedom.
Although the DGP model is appealing and elegant, it is also plagued by some problems. The first one is that, although the model does not possess ghosts on asymptotically flat manifolds, at the quantum level, it does have the problem of strong coupling for typical distances smaller than 1000 km, so that the theory is not easily under control [401]. Besides the model typically possesses superluminal modes. This may not directly violate causality, but it implies a non-trivial non-Lorentzian UV completion of the theory [304]. Also, on scales relevant for structure formation (between cluster scales and the Hubble radius), a quasi-static approximation to linear cosmological perturbations shows that the DGP model contains a ghost mode [369]. This linear analysis is valid as long as the Vainshtein radius * is smaller than the cluster scales. The original DGP model has been tested by using a number of observational data at the background level [525,238,405,9,549]. The joint constraints from the data of SN Ia, BAO, and the CMB shift parameter show that the flat DGP model is under strong observational pressure, while the open DGP model gives a slightly better fit [405,549]. Xia [622] showed that the parameter in the modified Friedmann equation 2 − / 2− = 2 (4) /3 [221] is constrained to be = 0.254 ± 0.153 (68% confidence level) by using the data of SN Ia, BAO, CMB, gamma ray bursts, and the linear growth factor of matter perturbations. Hence the flat DGP model ( = 1) is not compatible with current observations.
On the sub-horizon scales larger than the Vainshtein radius, the equation for linear matter perturbations in the DGP model was derived in [400,369] under a quasi-static approximation: where is the non-relativistic matter density on the brane and In the deep matter era one has ≫ 1 and hence ≃ − , so that is largely negative (| | ≫ 1). In this regime the evolution of is similar to that in GR ( ∝ 2/3 ). Since the background solution finally approaches the de Sitter solution characterized by Eq. (13.51), it follows that ≃ 1 − 2 ≃ −1 asymptotically. Since 1 + 1/(3 ) ≃ 2/3, the growth rate in this regime is smaller than that in GR.
The index of the growth rate = (Ω ) is approximated by ≈ 0.68 [395]. This is quite different from the value ≃ 0.55 for the ΛCDM model. If the future imaging survey of galaxies can constrain within 20%, it will be possible to distinguish the ΛCDM model from the DGP Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 model observationally [624]. We recall that in metric f (R) gravity the growth index today can be as small as = 0.4 because of the enhanced growth rate, which is very different from the value in the DGP model.
Comparing Eq. (13.53) with the effective gravitational constant (10.42) in BD theory with a massless limit (or the absence of the field potential), we find that the parameter BD has the following relation with : Since < 0 for the self-accelerating DGP solution, it follows that BD < −3/2. This corresponds to the theory with ghosts, because the kinetic energy of a scalar-field degree of freedom becomes negative in the Einstein frame [175]. There is a claim that the ghost may disappear for the Vainshtein radius * of the order of −1 0 , because the linear perturbation theory is no longer applicable [218]. In fact, a ghost does not appear in a Minkowski brane in the DGP model. In [370] it was shown that the Vainshtein radius in the early universe is much smaller than the one in the Minkowski background, while in the self accelerating universe they agree with each other. Hence the perturbative approach seems to be still possible for the weak gravity regime beyond the Vainshtein radius.
There have been studies regarding a possible regularization in order to avoid the ghost/strong coupling limit. Some of these studies have focused on smoothing out the delta profile of the Ricci scalar on the brane, by coupling the Ricci scalar to some other scalar field with a given profile [363,362]. In [516] the authors included the brane and bulk cosmological constants in addition to the scalar curvature in the action for the brane and showed that the effective equation of state of dark energy can be smaller than −1. A monopole in seven dimensions generated by a SO(3) invariant matter Lagrangian is able to change the gravitational law at its core, leading to a lower dimensional gravitational law. This is a first approach to an explanation of trapping of gravitons, due to topological defects in classical field theory [508,184]. Other studies have focused on re-using the delta function profile but in a higher-dimensional brane [334,333,197]. There is also an interesting work about the possibility of self-acceleration in the normal DGP branch [ = −1 in Eq. (13.50)] by considering an f (R) term on the brane action [97] (see also [4]). All these attempts indeed point to the direction that some mechanism, if not exactly DGP rather similar to it, may avoid a number of problems associated with the original DGP model.

Special symmetries
Since general covariance alone does not restrict the choice of the Lagrangian function, e.g., for f (R) theory, one can try to shrink the set of allowed functions by imposing some extra symmetry. In particular one can assume that the theory possesses a symmetry on some special background. However, allowing some theories to be symmetrical on some backgrounds does not imply these theories are viable by default. Nevertheless, this symmetry helps to give stronger constraints on them, as the allowed parameter space drastically reduces. We will discuss here two of the symmetries studied in the literature: Noether symmetries on a FLRW background and the Galileon symmetry on a Minkowski background.

Noether symmetry on FLRW
The action for metric f (R) gravity can be evaluated on a FLRW background, in terms of the fields ( ) and ( ), see [125,124] (and also [129,128,415,433,429,604,603,199,200,132]). Then the Lagrangian turns out to be non-singular, ℒ( ,˙), where 1 = , and 2 = . Its Euler-Lagrange equation is given by ( ℒ/˙) · − ℒ/ = 0. Contracting these equations with a vector function Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 ( ), (where 1 = and 2 = are two unknown functions of the ), we obtain A symmetry exists if the equation ℒ = 0 has non-trivial solutions. As a byproduct, the form of f (R), not specified in the point-like Lagrangian ℒ, is determined in correspondence to such a symmetry.
It can be proved that such do exist [124], and they correspond to where 1 , 2 , 3 are constants. However, in order that ℒ vanishes, one also needs to set the constraint (provided that 2 ̸ = 0) where 0 has a dimension of mass. For a general it is not possible to solve the Euler-Lagrange equation and the constraint equation (13.59) at the same time. Hence, we have to use the Noether constraint in order to find the subset of those which make this possible. Some partial solutions (only when 0 = 0) were found, but whether this symmetry helps finding viable models of f (R) is still not certain. However, the f (R) theories which possess Noether currents can be more easily constrained, as now the original freedom for the function in the Lagrangian reduced to the choice of the parameters and 0 .

Galileon symmetry
Recently another symmetry, the Galileon symmetry, for a scalar field Lagrangian was imposed on the Minkowski background [455]. This idea is interesting as it tries to decouple light scalar fields from matter making use of non-linearities, but without introducing new ghost degrees of freedom [455]. This symmetry was chosen so that the theory could naturally implement the Vainshtein mechanism. However, the same mechanism, at least in cosmology, seems to appear also in the FLRW background for scalar fields which do not possess such a symmetry (see [539,351,190]).
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 Keeping a universal coupling with matter (achieved through a pure nonminimal coupling with ), Nicolis et al. [455] imposed a symmetry called the Galilean invariance on a scalar field in the Minkowski background. If the equations of motion are invariant under a constant gradient-shift on Minkowski spacetime, that is → + + , (13.61) where both and are constants, we call a Galileon field. This implies that the equations of motion fix the field up to such a transformation. The point is that the Lagrangian must implement the Vainshtein mechanism in order to pass solar-system constraints. This is achieved by introducing self-interacting non-linear terms in the Lagrangian. It should be noted that the Lagrangian is studied only at second order in the fields (having a nonminimal coupling with ) and the metric itself, whereas the non-linearities are fully kept by neglecting their backreaction on the metric (as the biggest contribution should come only from standard matter). The equations of motion respecting the Galileon symmetry contain terms such as a constant, (up to fourth power), and other power contraction of the tensor ∇ ∇ . It is due to these non-linear derivative terms by which the Vainshtein mechanism can be implemented, as it happens in the DGP model [401].
Nicolis et al. [455] found that there are only five terms ℒ with = 1, . . . , 5 which can be inserted into a Lagrangian, such that the equations of motion respect the Galileon symmetry in 4-dimensional Minkowski spacetime. The first three terms are given by All these terms generate second-order derivative terms only in the equations of motion. The approach in the Minkowski spacetime has motivated to try to find a fully covariant framework in the curved spacetime. In particular, Deffayet et al. [205] found that all the previous 5-terms can be written in a fully covariant way. However, if we want to write down ℒ 4 and ℒ 5 covariantly in curved spacetime and keep the equations of motion free from higher-derivative terms, we need to introduce couplings between the field and the Riemann tensor [205]. The following two terms keep the field equations to second-order, where the last terms in Eqs. (13.65) and (13.66) are newly introduced in the curved spacetime. These terms possess the required symmetry in Minkowski spacetime, but mostly, they do not introduce derivatives higher than two into the equations of motion. In this sense, although originated from an implementation of the DGP idea, the covariant Galileon field is closer to the approach of the modifications of gravity in ( , ), that is, a formalism which would introduce only secondorder equations of motion. This result can be extended to arbitrary dimensions [202]. One can find, analogously to the Lovelock action-terms, an infinite tower of terms that can be introduced with the same property of keeping the equations of motion at second order. In particular, let us consider the action Here 1··· is the Levi-Civita tensor. The first product in Eq. (13.69) is defined to be one when = 0 and 0 for < 0, and the second product is one when = 1 + 2 , and 0 when < 2 + 2 . In ℒ ( +1, ) there will be + 1 powers of , and powers of the Riemann tensor. In four dimensions, for example, ℒ (1,0) and ℒ (2,0) are identical to ℒ 1 and ℒ 2 introduced before, respectively. Instead, ℒ (3,0) , ℒ (4,0) −(1/4) ℒ (4,1) , and ℒ (5,0) −(3/4)ℒ (5,1) reduce to ℒ 3 , ℒ 4 , and ℒ 5 , up to total derivatives, respectively.
In general non-linear terms discussed above may introduce the Vainshtein mechanism to decouple the scalar field from matter around a star, so that solar-system constraints can be satisfied. However the modes can have superluminal propagation, which is not surprising as the kinetic terms get heavily modified in the covariant formalism. Some studies have focused especially on the ℒ 3 term only, as this corresponds to the simplest case. For some models the background cosmological evolution is similar to that in the DGP model, although there are ghostlike modes depending on the sign of the time-velocity of the field [158]. There are some works for cosmological dynamics in Brans-Dicke theory in the presence of the non-linear term ℒ 3 [539,351,190] (although the original Galileon symmetry is not preserved in this scenario). Interestingly the ghost can disappear even for the case in which the Brans-Dicke parameter BD is smaller than −2. Moreover this theory leaves a number of distinct observational signatures such as the enhanced growth rate of matter perturbations and the significant ISW effect in CMB anisotropies.
At the end of this section we should mention conformal gravity in which the conformal invariance forces the gravitational action to be uniquely given by a Weyl action [414,340]. Interestingly the conformal symmetry also forces the cosmological constant to be zero at the level of the action [413]. It will be of interest to study the cosmological aspects of such theory, together with the possibility for the avoidance of ghosts and instabilities.

Conclusions
We have reviewed many aspects of f (R) theories studied extensively over the past decade. This burst of activities is strongly motivated by the observational discovery of dark energy. The idea is that the gravitational law may be modified on cosmological scales to give rise to the late-time acceleration, while Newton's gravity needs to be recovered on solar-system scales. In fact, f (R) theories can be regarded as the simplest extension of General Relativity.
The possibility of the late-time cosmic acceleration in metric f (R) gravity was first suggested by Capozziello in 2002 [113]. Even if f (R) gravity looks like a simple theory, successful f (R) dark energy models need to satisfy a number of conditions for consistency with successful cosmological evolution (a late-time accelerated epoch preceded by a matter era) and with local gravity tests on solar-system scales. We summarize the conditions under which metric f (R) dark energy models are viable: 1. , > 0 for ≥ 0 , where 0 is the Ricci scalar today. This is required to avoid a ghost state.

,
> 0 for ≥ 0 . This is required to avoid the negative mass squared of a scalar-field degree of freedom (tachyon).

( ) → − 2Λ for
≥ 0 . This is required for the presence of the matter era and for consistency with local gravity constraints.

0 <
, , ( = −2) < 1 at = − , = −2. This is required for the stability and the presence of a late-time de Sitter solution. Note that there is another fixed point that can be responsible for the cosmic acceleration (with an effective equation of state eff > −1).
We clarified why the above conditions are required by providing detailed explanation about the background cosmological dynamics (Section 4), local gravity constraints (Section 5), and cosmological perturbations (Sections 6 -8).
After the first suggestion of dark energy scenarios based on metric f (R) gravity, it took almost five years to construct viable models that satisfy all the above conditions [26,382,31,306,568,35,587]. In particular, the models (4.83), (4.84), and (4.89) allow appreciable deviation from the ΛCDM model during the late cosmological evolution, while the early cosmological dynamics is similar to that of the ΛCDM. The modification of gravity manifests itself in the evolution of cosmological perturbations through the change of the effective gravitational coupling. As we discussed in Sections 8 and 13, this leaves a number of interesting observational signatures such as the modification to the galaxy and CMB power spectra and the effect on weak lensing. This is very important to distinguish f (R) dark energy models from the ΛCDM model in future high-precision observations.
As we showed in Section 2, the action in metric f (R) gravity can be transformed to that in the Einstein frame. In the Einstein frame, non-relativistic matter couples to a scalar-field degree of freedom (scalaron) with a coupling of the order of unity ( = −1/ √ 6). For the consistency of metric f (R) gravity with local gravity constraints, we require that the chameleon mechanism [344,343] is at work to suppress such a large coupling. This is a non-linear regime in which the linear expansion of the Ricci scalar into the (cosmological) background value 0 and the perturbation is no longer valid, that is, the condition ≫ 0 holds in the region of high density. As long as a spherically symmetric body has a thin-shell, the effective matter coupling eff is suppressed to avoid the propagation of the fifth force. In Section 5 we provided detailed explanation about the chameleon mechanism in f (R) gravity and showed that the models (4.83) and (4.84) are consistent with present experimental bounds of local gravity tests for > 0.9.
Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 The construction of successful f (R) dark energy models triggered the study of spherically symmetric solutions in those models. Originally it was claimed that a curvature singularity present in the models (4.83) and (4.84) may be accessed in the strong gravitational background like neutron stars [266,349]. Meanwhile, for the Schwarzschild interior and exterior background with a constant density star, one can approximately derive analytic thin-shell solutions in metric f (R) and Brans-Dicke theory by taking into account the backreaction of gravitational potentials [594]. In fact, as we discussed in Section 11, a static star configuration in the f (R) model (4.84) was numerically found both for the constant density star and the star with a polytropic equation of state [43,600,42]. Since the relativistic pressure is strong around the center of the star, the choice of correct boundary conditions along the line of [594] is important to obtain static solutions numerically.
The model ( ) = + 2 /(6 2 ) proposed by Starobinsky in 1980 is the first model of inflation in the early universe. Inflation occurs in the regime ≫ 2 , which is followed by the reheating phase with an oscillating Ricci scalar. In Section 3 we studied the dynamics of inflation and (p)reheating (with and without nonminimal couplings between a field and ) in detail. As we showed in Section 7, this model is well consistent with the WMAP 5-year bounds of the spectral index of curvature perturbations and of the tensor-to-scalar ratio . It predicts the values of smaller than the order of 0.01, unlike the chaotic inflation model with = (0.1). It will be of interest to see whether this model continues to be favored in future observations.
Besides metric f (R) gravity, there is another formalism dubbed the Palatini formalism in which the metric and the connection Γ are treated as independent variables when we vary the action (see Section 9). The Palatini f (R) gravity gives rise to the specific trace equation (9.2) that does not have a propagating degree of freedom. Cosmologically we showed that even for the model ( ) = − / ( > 0, > −1) it is possible to realize a sequence of radiation, matter, and de Sitter epochs (unlike the same model in metric f (R) gravity). However the Palatini f (R) gravity is plagued by a number of shortcomings such as the inconsistency with observations of large-scale structure, the conflict with Standard Model of particle physics, and the divergent behavior of the Ricci scalar at the surface of a static spherically symmetric star with a polytropic equation of state = Γ 0 with 3/2 < Γ < 2. The only way to avoid these problems is that the f (R) models need to be extremely close to the ΛCDM model. This property is different from metric f (R) gravity in which the deviation from the ΛCDM model can be significant for of the order of the Ricci scalar today.
In Brans-Dicke (BD) theories with the action (10.1), expressed in the Einstein frame, nonrelativistic matter is coupled to a scalar field with a constant coupling . As we showed in in Section 10.1, this coupling is related to the BD parameter BD with the relation 1/(2 2 ) = 3 + 2 BD . These theories include metric and Palatini f (R) gravity theories as special cases where the coupling is given by = −1/ √ 6 (i.e., BD = 0) and = 0 (i.e., BD = −3/2), respectively. In BD theories with the coupling of the order of unity we constructed a scalar-field potential responsible for the late-time cosmic acceleration, while satisfying local gravity constraints through the chameleon mechanism. This corresponds to the generalization of metric f (R) gravity, which covers the models (4.83) and (4.84) as specific cases. We discussed a number of observational signatures in those models such as the effects on the matter power spectrum and weak lensing.
Besides the Ricci scalar , there are other scalar quantities such as and constructed from the Ricci tensor and the Riemann tensor . For the Gauss-Bonnet (GB) curvature invariant ≡ 2 − 4 + one can avoid the appearance of spurious spin-2 ghosts. There are dark energy models in which the Lagrangian density is given by ℒ = + ( ), where ( ) is an arbitrary function in terms of . In fact, it is possible to explain the late-time cosmic acceleration for the models such as (12.16) and (12.17), while at the same time local gravity constraints are satisfied. However density perturbations in perfect fluids exhibit violent negative instabilities during both the radiation and the matter domination, irrespective of the form of ( ). The growth of perturbations gets stronger on smaller scales, Living Reviews in Relativity http://www.livingreviews.org/lrr-2010-3 which is incompatible with the observed galaxy spectrum unless the deviation from GR is very small. Hence these models are effectively ruled out from this Ultra-Violet instability. This implies that metric f (R) gravity may correspond to the marginal theory that can avoid such instability problems.
In Section 13 we discussed other aspects of f (R) gravity and modified gravity theories -such as weak lensing, thermodynamics and horizon entropy, Noether symmetry in f (R) gravity, unified f (R) models of inflation and dark energy, f (R) theories in extra dimensions, Vainshtein mechanism, DGP model, and Galileon field. Up to early 2010 the number of papers that include the word "f (R)" in the title is over 460, and more than 1050 papers including the words "f (R)" or "modified gravity" or "Gauss-Bonnet" have been written so far. This shows how this field is rich and fruitful in application to many aspects to gravity and cosmology.
Although in this review we have focused on f (R) gravity and some extended theories such as BD theory and Gauss-Bonnet gravity, there are other classes of modified gravity theories, e.g., Einstein-Aether theory [325], tensor-vector-scalar theory of gravity [76], ghost condensation [38], Lorentz violating theories [144,282,389], and Hořava-Lifshitz gravity [305]. There are also attempts to study f (R) gravity in the context of Hořava-Lifshitz gravity [346,347]. We hope that future high-precision observations can distinguish between these modified gravity theories, in connection to solving the fundamental problems for the origin of inflation, dark matter, and dark energy.