$f(R)$ gravity is an extension of Einstein's General Relativity derived from relaxing the hypothesis that the Hilbert-Einstein action for the gravitational field is strictly linear in the Ricci curvature scalar $R$, i.e. $f(R)= R$. In this sense, $f(R)$ gravity represents a class of theories defined as arbitrary functions of $R$. It can be considered as the simplest example of Extended Theory of Gravity (Capozziello and De Laurentis, 2011).
In $f(R)$ gravity, the Hilbert-Einstein action: \begin{eqnarray}\tag{1}{\cal A}[g]=\int R \sqrt{-g}\,\mathrm{d}^4x\,, \end{eqnarray} is generalized as: \begin{eqnarray}\tag{2}{\cal A}[g]=\int f(R) \sqrt{-g} \, \mathrm{d}^4x\,, \end{eqnarray} where we are using physical units $8\pi G c^{-4}=1$, the determinant of the metric tensor is $g=|g_{\mu\nu}|$ and f(R) is a function of the Ricci curvature scalar $R$. This approach is aimed to address problems and shortcomings coming from the Standard Cosmological Model and from the lack of a final theory of Quantum Gravity. Due to these issues, alternative theories of gravity have been considered in order to attempt, at least, a semi-classical scheme where General Relativity and its positive results has to be recovered. One of the most fruitful approaches is that of Extended Theories of Gravity (see (Capozziello and De Laurentis, 2011) and references therein) which have become a sort of paradigm in the study of gravitational interaction. They are based on corrections and enlargements of the Einstein theory. The paradigm consists, essentially, in adding higher-order curvature invariants and minimally or non-minimally coupled scalar fields into dynamics related to the effective actions of some fundamental theory (Nojiri and Odintsov, 2011) (Nojiri and Odintsov, 2007) (Capozziello and Francaviglia, 2007) (De Felice and Tsujikawa, 2010) (Sotiriou and Faraoni, 2010). This approach is coherent to the fact that these generalized theories emerge, like Einstein's gravity, from the Gauge Theory, and can be framed into a general bundle structure (Capozziello and De Laurentis, 2010). Other motivations to modify General Relativity come from the Mach Principle which leads to assume a varying gravitational coupling. This principle states that the local inertial frame is determined by some average of the motion of distant astronomical objects (Bondi, 1952). This fact implies that the gravitational coupling can be scale-dependent and related to some scalar field. As a consequence, the concept of inertia and the Equivalence Principle have to be revised. Besides, any unification scheme as Superstrings, Supergravity or Grand Unified Theories, takes into account effective actions where non-minimal couplings to the geometry or higher-order terms of the curvature invariants are present. Such contributions are due to one-loop or higher-loop corrections in the strong gravity regimes near the full (not yet available) Quantum Gravity (Nojiri and Odintsov, 2007). Specifically, this scheme was adopted in order to deal with the quantization on curved space-times and the result was that the interactions among quantum scalar fields and background geometry or the gravitational self-interactions yield corrective terms in the Hilbert-Einstein Lagrangian (Birrell and Davies, 1982). Moreover, it has been realized that such corrective terms are inescapable in order to obtain the effective action of Quantum Gravity at scales closed to the Planck one (Vilkovisky, 1992). All these approaches are not the full Quantum Gravity but are needed as working schemes towards it. In summary, higher-order terms in curvature invariants (such as $R^{2}$, $R_{\mu\nu} R^{\mu\nu}$, $R^{\mu\nu\alpha\beta}R_{\mu\nu\alpha\beta}$, $R \,\Box R$, or $R \,\Box^{k}R$) or non-minimally coupled terms between scalar fields and geometry (such as $\phi^{2}R$) have to be added to the effective Lagrangian of gravitational field if quantum corrections are considered. For instance, one can notice that such terms occur in the effective Lagrangian of strings or in Kaluza-Klein theories, if the mechanism of dimensional reduction is adopted (Gasperini and Veneziano, 1992).
On the other hand, from a conceptual point of view, there are no a priori reason to restrict the gravitational Lagrangian to a linear function of the Ricci scalar $R$ minimally coupled with matter (Capozziello and Francaviglia, 2007). Furthermore, the idea that there are no exact laws of physics could be taken into serious account: in such a case, the effective Lagrangians of physical interactions are stochastic functions. This feature means that the local gauge invariances (i.e. conservation laws) are well approximated only in the low energy limit and the fundamental physical constants can vary (Barrow and Ottewil, 1983).
Beside fundamental physics motivations, all these theories acquired a huge interest in cosmology due to the fact that they naturally exhibit inflationary behaviours able to overcome the shortcomings of Standard Cosmological Model (based on General Relativity). The related cosmological models seem realistic and capable of matching with the Cosmic Microwave Background Radiation observations (Starobinsky, 1980) (Duruisseau and Kerner, 1983) (La and Steinhardt, 1989). Furthermore, it is possible to show that, via conformal transformations, the higher-order and non-minimally coupled terms correspond to the Einstein gravity plus one or more than one minimally coupled scalar fields (Teyssandier and Tourrenc, 1983) (Maeda, 1989) (Wands, 1994) (Capozziello, de Ritis and Marino, 1998) (Adams, Freese and Guth, 1991).
More precisely, higher-order terms appear as contributions of order two in the derivatives of field equations. For example, a term like $R^{2}$ gives fourth order equations (Ruzmaikina and Ruzmaikin, 1970), $R \ \Box R$ gives sixth order equations (Amendola et al., 1993) (Adams, Freese and Guth, 1991) (Buchdahl, 1951), $R\,\Box^{2}R$ gives eighth order equations (Battaglia-Mayer and Schmidt, 1993) and so on. By a conformal transformation, any 2nd-order derivative term corresponds to a scalar field. The dynamics of any scalar field is given by the corresponding Klein-Gordon equation, which is second order. For example, fourth-order gravity gives Einstein's equations plus one scalar field, sixth-order gravity gives Einstein's equations plus two scalar fields and so on (Schmidt, 1990).
Considering the mathematical point of view, the problem of reducing more general theories to the Einstein standard form has been extensively treated; one can see that, through a Legendre transformation, higher-order theories, under suitable regularity conditions on the Lagrangian, take the form of the Einstein one in which a scalar field (or more than one) is the source of the gravitational field (see for example (Capozziello and Francaviglia, 2007) (Sokolowski, 1989) (Ferraris, Francaviglia and Magnano, 1989) (Magnano and Sokolowski, 1994) (Jakubiec and Kijowski, 1988)).
In any case, the debate on the physical meaning of conformal transformations is far to be solved (see (Capozziello and Faraoni, 2010) (Capozziello and De Laurentis, 2010) (Capozziello, De Laurentis and Faraoni, 2009) (Faraoni, 2008), and references therein for a comprehensive review). Several authors claim for a true physical difference between Jordan frame (higher-order theories and/or variable gravitational coupling) since there are experimental and observational evidences which point out that the Jordan frame could be suitable to better match solutions with data. Others state that the true physical frame is the Einstein one according to the energy theorems (Magnano and Sokolowski, 1994). However, the discussion is open and no definitive statement has been formulated up to now.
The problem can be faced from a more general point of view and the Palatini approach to gravity, (see below) could be useful in this context (Einstein, 1925) (Buchdahl, 1979). In (Olmo, 2011), this approach is widely discussed for Extended Theories of Gravity and several important applications are reported. The fundamental idea of the Palatini formalism is to consider the connection $\Gamma$, entering the definition of the Ricci tensor, to be independent of the metric $g$ defined on the space-time ${\cal M}$. The Palatini formulation for the Hilbert-Einstein theory results to be equivalent to the purely metric theory: this follows from the fact that the field equations for the connection $\Gamma$, firstly considered to be independent of the metric, give the Levi-Civita connection of the metric $g$. As a consequence, there is no reason to impose the Palatini variational principle in the Hilbert-Einstein theory instead of the metric variational principle. However, the situation completely changes if we consider theories of gravity depending on functions of curvature invariants, as $f(R)$, or non-minimally coupled to some scalar field. In these cases, the Palatini and the metric variational principle provide different field equations and the theories thus derived differ (Magnano and Sokolowski, 1994) (Ferraris, Francaviglia and Volovich, 1994). The relevance of the Palatini approach, in this framework, has been recently proven in relation to cosmological applications (Vollick, 2003) (Li and Chu, 2006) (Li, Chan and Chu, 2007) (Capozziello, De Laurentis, Francaviglia and Mercadante, 2009) (Olmo, 2011).
From a physical point of view, considering the metric $g$ and the connection $\Gamma$ as independent fields means to decouple the metric structure of space-time and its geodesic structure (being, in general, the connection $\Gamma$ not the Levi-Civita connection of $g$). The chronological structure of space-time is governed by $g$ while the trajectories of particles, moving in the space-time, are governed by $\Gamma$. This decoupling enriches the geometric structure of space-time and generalizes the purely metric formalism. The metric-affine structure of space-time is naturally translated, by the Palatini field equations, into a bi-metric structure. Beside the physical metric $g$, another metric $\hat{g}$ is involved. This new metric is related, in the case of $f(R)$ gravity, to the connection. As a matter of fact, the connection $\Gamma$ results to be the Levi-Civita connection of $\hat{g}$ and thus provides the geodesic structure (Allemandi, Capone, Capozziello and Francaviglia, 2006).
If we consider the case of non-minimally coupled interaction in the gravitational Lagrangian (scalar-tensor theories), the new metric $\hat{g}$ is related to the non-minimal coupling. The new metric $\hat{g}$ can be thus related to different geometric and physical aspects of the gravitational theory. Thanks to the Palatini formalism, the non-minimal coupling and the scalar field, entering the evolution of the gravitational fields, are separated from the metric structure of space-time. The situation mixes when we consider the case of higher-order-scalar-tensor theories. Due to these features, the Palatini approach could greatly contribute to clarify the physical meaning of conformal transformation (Allemandi, Capone, Capozziello and Francaviglia, 2006).
In metric $f(R)$ gravity, the field equations are obtained by varying with respect to the metric and not treating the connection independently. The main steps are the same as in the case of the variation of the Hilbert-Einstein action but there are also some important differences.
The variation of the determinant is \begin{eqnarray}\tag{3}\delta \sqrt{-g}= -\frac{1}{2} \sqrt{-g} g_{\mu\nu} \delta g^{\mu\nu}\,.\end{eqnarray} The Ricci scalar is defined as \begin{eqnarray}\tag{4} R = g^{\mu\nu} R_{\mu\nu}.\end{eqnarray} Therefore the variation with respect to the inverse metric gμν is given by \[ \begin{align}\tag{5} \delta R &= R_{\mu\nu} \delta g^{\mu\nu} + g^{\mu\nu} \delta R_{\mu\nu} = R_{\mu\nu} \delta g^{\mu\nu} + g^{\mu\nu}(\nabla_\rho \delta \Gamma^\rho_{\nu\mu} - \nabla_\nu \delta \Gamma^\rho_{\rho\mu})\,. \end{align} \] Since δΓλμν is the difference of two connections, it should transform as a tensor. Therefore, it can be written as \[\tag{6}\delta \Gamma^\lambda_{\mu\nu}=\frac{1}{2}g^{\lambda \alpha}\left(\nabla_\mu\delta g_{\alpha\nu}+\nabla_\nu\delta g_{\alpha\mu}-\nabla_\alpha\delta g_{\mu\nu} \right).\] Substituting into the above equation: \[\tag{7}\delta R= R_{\mu\nu} \delta g^{\mu\nu}+g_{\mu\nu}\Box \delta g^{\mu\nu}-\nabla_\mu \nabla_\nu \delta g^{\mu\nu}\,,\] where ∇μ is the covariant derivative and $\Box= g^{\mu\nu} \nabla_{\mu}\nabla_{\nu}$ is the D'Alembert operator. We adopt the metric signature $(+,-,-,-)$. Let us consider now a generic analytic function $f(R)$ obeying the variational principle $ \delta \int d^{4}x \sqrt{-g} \, f(R) =0$. We have \begin{eqnarray} && \delta \int d^{4}x \sqrt{-g} \, f(R) = \int d^{4}x \left[ \delta \left( \sqrt{-g} \, f(R) \right) + \sqrt{-g} \, \delta\left(f(R) \right)\right] \nonumber\\ &&\nonumber\\ && = \int d^{4}x \sqrt{-g} \, {\textstyle\left[f'(R)R_{\mu\nu} -\frac{1}{2}g_{\mu\nu}f(R)\right]}\delta g^{\mu\nu} + \int d^{4}x \sqrt{-g} \, f'(R)g^{\mu\nu}\delta R_{\mu\nu} \, ,\nonumber\\ && \tag{8} \end{eqnarray} where the prime denotes differentiation with respect to $R$. These integrals can be evaluated in the local inertial frame. By using \begin{eqnarray}\tag{9} g^{\mu\nu}\delta R_{\mu\nu} =g^{\mu\nu}\partial_{\sigma} \left(\delta G_{\mu\nu}^{\sigma}\right) -g^{\mu\sigma}\partial_{\sigma}\left(\delta G_{\mu\nu}^{\nu}\right) \equiv\partial_{\sigma}W^{\sigma}\,, \end{eqnarray} where \begin{eqnarray}\tag{10} W^{\sigma}\equiv g^{\mu\nu} \delta G_{\mu\nu}^{\sigma} -g^{\mu\sigma}\delta G_{\mu\nu}^{\nu} \,, \end{eqnarray} the second integral in Eq. (8) can be written as \begin{eqnarray}\tag{11} \int d^{4}x \sqrt{-g} \, f'(R)g^{\mu\nu}\delta R_{\mu\nu} =\int d^{4}x \sqrt{-g} \, f'(R) \partial_{\sigma}W^{\sigma} \, . \end{eqnarray} Integration by parts yields \begin{eqnarray}\tag{12} \int d^{4}x \sqrt{-g} \, f'(R)g^{\mu\nu}\delta R_{\mu\nu} &= & \int d^{4}x \frac{\partial}{\partial x^{\sigma} } \left[ \sqrt{-g} \, f'(R)W^{\sigma}\right] - \int d^{4}x \partial_{\sigma}\left[\sqrt{-g}f'(R)\right]W^{\sigma}\, . \end{eqnarray} The first integrand is a total divergence and can be discarded assuming that fields vanish at infinity. We have \begin{eqnarray}\tag{13} \int d^{4}x \sqrt{-g} \, f'(R) g^{\mu\nu}\delta R_{\mu\nu} = -\int d^{4}x \partial_{\sigma} \left[\sqrt{-g} \, f'(R)\right]W^{\sigma} \,. \end{eqnarray}
Let us evaluate $W^{\sigma}$ appearing in Eq. (13). We have \begin{eqnarray}\tag{14} \delta G_{\mu\nu}^{\sigma} = \delta\left[\frac{1}{2}g^{\sigma\alpha}\left(\partial_{\mu}g_{\alpha\nu} +\partial_{\nu}g_{\mu\alpha}-\partial_{\alpha}g_{\mu\nu}\right)\right] = \frac{1}{2}g^{\sigma\alpha}\left[\partial_{\mu}\left(\delta g_{\alpha\nu}\right)+\partial_{\nu}\left(\delta g_{\mu\alpha}\right)-\partial_{\alpha}\left(\delta g_{\mu\nu}\right)\right] \,. \end{eqnarray}
In the local inertial frame, it is \begin{eqnarray}\tag{15} \partial_{\alpha}g_{\mu\nu}=\nabla_{\alpha}g_{\mu\nu}= 0\,, \end{eqnarray} and \begin{eqnarray}\tag{16} \delta G_{\mu\nu}^{\nu}=\frac{1}{2} \, g^{\nu\alpha}\partial_{\mu}\left(\delta g_{\nu\alpha}\right) \,. \end{eqnarray} By combining Eqs. (15) and (16), one obtains \begin{eqnarray} g^{\mu\nu}\delta G_{\mu\nu}^{\sigma} & = & \frac{1}{2} \, g^{\mu\nu}\left[ -\partial_{\mu}\left(g_{\alpha\nu}\delta g^{\alpha\sigma}\right)-\partial_{\nu}\left(g_{\mu\alpha}\delta g^{\sigma\alpha}\right)-g^{\sigma\alpha}\partial_{\alpha}\left(\delta g_{\mu\nu}\right)\right] = \frac{1}{2} \, \partial^{\sigma}\left(g_{\mu\nu}\delta g^{\mu\nu} \right)-\partial^{\mu}\left(g_{\alpha\mu}\delta g^{\nu\alpha}\right) \,, \tag{17} \\ &&\nonumber\\ g^{\mu\sigma}\delta G_{\mu\nu}^{\nu} & = & -\frac{1}{2} \, \partial^{\sigma}\left(g_{\nu\alpha}\delta g^{\nu\alpha}\right) \,,\tag{18} \end{eqnarray} from which it follows that \begin{eqnarray}\tag{19} W^{\sigma}=\partial^{\sigma}\left(g_{\mu\nu}\delta g^{\mu\nu}\right)-\partial^{\mu}\left(g_{\mu\nu}\delta g^{\sigma\nu}\right) \,. \end{eqnarray} Using this equation, one can write \begin{eqnarray} \int d^{4}x \sqrt{-g} \, f'(R)g^{\mu\nu}\delta R_{\mu\nu} = \int d^{4}x \partial_{\sigma}\left[\sqrt{-g} \, f'(R)\right] \left[\partial^{\mu}\left(g_{\mu\nu}\delta g^{\sigma\nu}\right)-\partial^{\sigma}\left(g_{\mu\nu}\delta g^{\mu\nu}\right)\right] \,. \tag{20} \end{eqnarray} Integrating by parts and discarding total divergences, one obtains \begin{eqnarray} \int d^{4}x \, \sqrt{-g} \, f'(R)g^{\mu\nu}\delta R_{\mu\nu} = \int d^{4}x \, g_{\mu\nu} \partial^{\sigma}\partial_{\sigma}\left[\sqrt{-g}f'(R)\right]\delta g^{\mu\nu} - \int d^{4}x \, g_{\mu\nu}\partial^{\mu}\partial_{\sigma}\left[\sqrt{-g} \, f'(R)\right] \delta g^{\sigma\nu} \,. \tag{21} \end{eqnarray} The variation of the action is then \begin{eqnarray}\tag{22} \delta\int d^{4}x \sqrt{-g} \, f(R) =\int d^{4}x \sqrt{-g}\,\left[ f'(R)R_{\mu\nu}-\frac{1}{2} f(R) g_{\mu\nu}\right]\delta g^{\mu\nu}+\int d^{4}x \left[ g_{\mu\nu}\partial^{\sigma}\partial_{\sigma} \left(\sqrt{-g}\, f'(R)\right)-g_{\sigma\nu}\partial^{\mu} \partial_{\sigma}\left(\sqrt{-g}f'(R)\right) \right]\delta g^{\mu\nu} \,. \end{eqnarray} The vanishing of the variation implies the fourth order vacuum field equations \begin{eqnarray}\tag{23} f'(R)R_{\mu\nu}-\frac{f(R)}{2} \, g_{\mu\nu} =\nabla_{\mu} \nabla_{\nu}f'(R)-g_{\mu\nu}\Box f'(R)\,. \end{eqnarray} These equations can be rearranged in the Einstein-like form \begin{eqnarray}\tag{24} f'(R) R_{\mu\nu}-\frac{ f'(R)}{2} \, g_{\mu\nu} R +\frac{f'(R)}{2} \, g_{\mu\nu}R-\frac{f(R)}{2} \, g_{\mu\nu}= \nabla_{\mu} \nabla_{\nu} f'(R)- g_{\mu\nu}\Box f'(R)\,, \end{eqnarray} and then \begin{eqnarray}\tag{25} G_{\mu\nu}=\frac{1}{f'(R)} \left\{ \nabla_{\mu}\nabla_{\nu} f'(R) - g_{\mu\nu}\Box f'(R) + g_{\mu\nu} \frac{ \left[ f(R)-f'(R) R \right]}{2}\right\}\,. \end{eqnarray} The right hand side of Eq. (25) is then regarded as an effective stress-energy tensor $T_{\mu\nu}^{(curv)}$, which we call curvature energy-momentum tensor, sourcing the effective Einstein equations. Although this interpretation is questionable, in principle, because the field equations describe a theory different from General Relativity, and one is forcing upon them the interpretation as effective Einstein equations, this approach is quite useful in practice.
Quadratic corrections in the Ricci scalar, motivated by the attempts to renormalize General Relativity, constitute a straightforward extension of General Relativity and are particularly relevant in cosmology since they allow the construction of a self-consistent inflationary model (Starobinsky, 1980). Let us derive the field equations for the Lagrangian density \begin{eqnarray} \tag{26} {\cal L}=R+\alpha R^{2}+{\cal L}^{(m)}\,, \end{eqnarray} from the variational principle $\delta\int d^{4}x \sqrt{-g} \, {\cal L}=0 $. We consider first the vacuum case. The variation gives \begin{eqnarray} \tag{27} \int d^{4}x \sqrt{-g} \, G_{\alpha\beta}\delta g^{\alpha\beta} + \alpha \, \delta\int d^{4}x \sqrt{-g} \, R^{2} =0\,, \end{eqnarray} where the variation of $R\sqrt{-g} $ produces the Einstein tensor. The second term on the right hand side of Eq. (27) gives \begin{eqnarray} \tag{28} \delta \int d^{4}x \sqrt{-g} \, R^{2} = -\frac{1}{2}\int d^{4}x \sqrt{-g} \,\, g_{\alpha\beta}\delta g^{\alpha\beta}R^{2} + 2\int d^{4}x \sqrt{-g} \, R\delta R \end{eqnarray} and \begin{eqnarray} \tag{29} \int d^{4}x \sqrt{-g} \, R\delta R = \int d^{4}x \sqrt{-g} \, R\left(\delta g^{\alpha\beta} R_{\alpha\beta} + g^{\alpha\beta}\delta R_{\alpha\beta} \right) \,. \end{eqnarray} By using \begin{eqnarray}\tag{30} g^{\alpha\beta}\delta R_{\alpha\beta} = \nabla_{\alpha}\nabla_{\beta} h^{\alpha\beta} -\Box h \,, \end{eqnarray} where \begin{eqnarray} \tag{31} h^{\alpha\beta} \equiv -\delta g^{\alpha\beta},\ \;\;\;\;\; h \equiv -g_{\alpha\beta}\delta g^{\alpha\beta} \,, \end{eqnarray} one has \begin{eqnarray}\tag{32} \int d^{4}x \sqrt{-g} \, R \, g^{\alpha\beta}\delta R_{\alpha\beta} = \int d^{4}x \sqrt{-g} \, R \left( \nabla_{\alpha}\nabla_{\beta} h^{\alpha\beta} -\Box h\right) \,. \end{eqnarray} Integrating by parts twice, the operators $\nabla_{\alpha}\nabla_{\beta}$ and $\Box$ acting on $ h^{\alpha\beta}$ and $h$, respectively, transfer their action onto $R$ and \begin{eqnarray}\tag{33} \int d^{4}x \sqrt{-g} \, R \, g^{\alpha\beta}\delta R_{\alpha\beta} = \int d^{4}x \sqrt{-g} \, \left(h ^{\alpha\beta} \nabla_{\alpha}\nabla_{\beta} R -h\Box R\right) \,. \end{eqnarray} Using Eq.(31), Eq. (33) becomes \begin{eqnarray}\tag{34} \int d^{4}x \sqrt{-g} \, R \, g^{\alpha\beta}\delta R_{\alpha\beta} =\int d^{4}x \sqrt{-g} \left(-\delta g^{\alpha\beta} \nabla_{\alpha}\nabla_{\beta} R +g_{\alpha\beta}\Box R\delta g^{\alpha\beta}\right) \,. \end{eqnarray} Upon substitution of Eq.(34) into Eq.(29), one obtains \begin{eqnarray} \tag{35} \int d^{4}x \sqrt{-g} \, R\delta R = \int d^4x \, \sqrt{-g} \left(R\delta g^{\alpha\beta}R_{\alpha\beta}-\delta g^{\alpha\beta} \nabla_{\alpha}\nabla_{\beta} R +g_{\alpha\beta}\Box R\delta g^{\alpha\beta}\right)\,, \end{eqnarray} and Eq.(28) takes the form \begin{eqnarray} \delta \int d^{4}x \sqrt{-g} \, R^{2} & = & -\frac{1}{2}\int d^{4}x \sqrt{-g} \, g_{\alpha\beta}\delta g^{\alpha\beta} R^{2}+ 2\int d^{4}x \sqrt{-g} \, \left(R\delta g^{\alpha\beta}R_{\alpha\beta}-\delta g^{\alpha\beta} \nabla_{\alpha}\nabla_{\beta} R +g_{\alpha\beta} \Box R\delta g^{\alpha\beta}\right) \nonumber\\ &&\nonumber\\ & =& \int d^{4}x \sqrt{-g} \, {\textstyle\left(2RR_{\alpha\beta} -\frac{1}{2} \, g_{\alpha\beta}R^{2}\right)}\delta g ^{\alpha\beta} + 2 \int d^{4}x \sqrt{-g} \, \left(g_{\alpha\beta}\Box R- \nabla_{\alpha}\nabla_{\beta} R \right)\delta g^{\alpha\beta}\,. \tag{36} \end{eqnarray} Substituting this equation into Eq.(27) and including the matter part of the Lagrangian ${\cal L}^{(m)}$ which produces the energy-momentum tensor $T_{\mu\nu}^{(m)}$, the field equations \begin{eqnarray}\tag{37} G_{\alpha\beta} + \alpha\left[ 2R\left( R_{\alpha\beta} -\frac{1}{4} g_{\alpha\beta}R\right) + 2\left( g_{\alpha\beta}\Box R-\nabla_{\alpha}\nabla_{\beta} R\right)\right] = T_{\alpha\beta}^{(m)}\,, \end{eqnarray} are obtained.
The trace of Eq. (37) is \begin{eqnarray}\tag{38} \Box R -\frac{1}{6\alpha}\left(R+T^{(m)} \right)=0\,. \end{eqnarray} One can also define an angular frequency $\omega$ (equivalent to a mass $m$) so that \begin{eqnarray}\tag{39} -\frac{1}{6\alpha}=\omega^{2}=m^{2} \,. \end{eqnarray} Following this definition, Eq. (38) becomes \begin{eqnarray}\tag{40} \Box R + m^{2}\left(R+T^{(m)} \right)=0 \,. \end{eqnarray} Eq. (40) can be seen as an effective Klein-Gordon equation for the effective scalar field degree of freedom $R$ (sometimes called scalaron).
We can simplify the analysis of f(R) theories by introducing an auxiliary field Φ. Assuming $f''(R)\neq 0$ for all R, let V(Φ) be the Legendre transform of f(R) so that Φ = f′(R) and R = V′(Φ). One obtains the action (O'Hanlon, 1972)
\[\tag{41}{\cal A} = \int d^4x \sqrt{-g} \left[ \Phi R - V(\Phi) + \mathcal{L}_{\text{m}}\right].\]
The Euler–Lagrange equations are \[\tag{42}V'(\Phi)=R\,,\] \[\tag{43}\Phi \left( R_{\mu\nu} - \frac{1}{2}g_{\mu\nu} R \right) + \left(g_{\mu\nu}\Box -\nabla_\mu \nabla_\nu \right) \Phi + \frac{1}{2} g_{\mu\nu}V(\Phi) = T_{\mu\nu}\,.\] Substituting Φ with f′(R) , we obtain exactly the same equations as before. However, the equations are now second order in the derivatives, instead of fourth order. By performing a conformal rescaling \[\tag{44}\tilde{g}_{\mu\nu}=\Phi g_{\mu\nu},\] we transform to the Einstein frame: \[\tag{45}R=\Phi^{-1} \left[ \tilde{R} + \frac{3\tilde{\Box} \Phi}{\Phi} +\frac{9}{2}\left(\frac{\tilde{\nabla} \Phi}{\Phi}\right)^2 \right]\,,\] and then \[\tag{46}{\cal A} = \int d^4x \sqrt{-\tilde{g}}\left[ -\frac{1}{2}\tilde{R} + \frac{3}{2}\left( \frac{\tilde{\nabla}\Phi}{\Phi} \right)^2 - \frac{V(\Phi)}{\Phi^2} \right]\,,\] after integrating by parts.
We define \[\tag{47}\tilde{\Phi} = \sqrt{3} \ln{\Phi}\,,\] and substitute \[\tag{48}{\cal A} = \int \mathrm{d}^4x \sqrt{-\tilde{g}}\left[ -\frac{1}{2}\tilde{R} + \frac{1}{2}\left(\tilde{\nabla}\tilde{\Phi}\right)^2 - \tilde{V}(\tilde{\Phi}) \right]\,,\] where \[\tag{49}\tilde{V}(\tilde{\Phi}) = e^{-2/\sqrt{3}\;\tilde{\Phi}}V(e^{\tilde{\Phi}/\sqrt{3}}).\] This is General Relativity minimally coupled to a scalar field.
The fundamental idea of the Palatini formalism is to consider the connection $\Gamma^{\mu}_{\nu\alpha}$, entering the definition of the Ricci tensor, as a variable independent of the space-time metric $ g_{\mu\nu}$. The Palatini formulation of General Relativity is equivalent to the metric version of this theory as a consequence of the fact that the field equations for the connection $\Gamma^{\alpha}_{\mu\nu}$ give exactly the Levi-Civita connection of the metric $g_{\mu\nu}$ (Wald, 1984). As a consequence, there is no particular reason to impose the Palatini variational principle in General Relativity instead of the metric variational principle.
The situation is different in Extended Theories of Gravity depending on functions of curvature invariants or non-minimally coupled scalar fields. In these cases, the Palatini and the metric variational principle yield different field equations (Ferraris, Francaviglia and Volovich, 1994) (Magnano and Sokolowski, 1994).
The Palatini approach, in the context of Extended Theories of Gravity has been the subject of much interest in cosmological applications (Olmo, 2011). As discussed above, considering the metric $g_{\mu\nu}$ and the connection $\Gamma^{\alpha}_{\mu\nu}$ as independent fields amounts to decoupling the metric structure of space-time and its geodesic structure with the connection $\Gamma^{\alpha}_{\mu\nu}$ being distinct from the Levi-Civita connection of $g_{\mu\nu}$. In principle, this decoupling enriches the geometric structure of space-time and generalizes the purely metric formalism. By means of the Palatini field equations, this dual structure of space-time is naturally translated into a bimetric structure of the theory: instead of a metric and an independent connection, the Palatini formalism can be seen as containing two independent metrics $g_{\mu\nu}$ and ${\hat g}_{\mu\nu}=f'(R) \, g_{\mu\nu}$. In Palatini $f(R)$ gravity, the new metric ${\hat g}_{\mu\nu}$ determining the geodesics is related to the connection $\Gamma^{\alpha}_{\mu\nu}$ by the fact that the latter turns out to be the Levi-Civita connection of ${\hat g}_{\mu\nu}$. In scalar-tensor gravity, the second metric ${\hat g}_{\mu\nu}$ is related to the non-minimal coupling of the Brans-Dicke-like scalar. In the Palatini formalism the non-minimal coupling and the scalar field are separated from the metric structure of space-time. Physical consequences of this fact are discussed in (Capozziello, De Laurentis, Francaviglia and Mercadante, 2009) (Amendola, Enqvist and Koivisto, 2011). However, also other geometrical invariants, besides $R$, can be considered in the Palatini formalism. In (Li, Mota and Shaw, 2008), microscopic and macroscopic behaviors of Palatini modified gravity theories are discussed, in particular a detailed study of $f(R,R^{\mu\nu}R_{\mu\nu})$ models is reported. In (Olmo, Sanchis-Alepuz and Tripathi, 2009) dynamical aspects of Palatini $f(R)$, $f(R^{\mu\nu}R_{\mu\nu})$, and $f(R,R^{\mu\nu}R_{\mu\nu})$ theories are studied. Isotropic and anisotropic bouncing cosmologies in Palatini $f(R)$ and $f(R^{\mu\nu}R_{\mu\nu})$ theories are discussed in (Barragan and Olmo , 2010). A Lagrangian of the type $\mathcal{L}=R+\alpha R_{\mu\nu}R^{\mu\nu} $ is also studied in Palatini formalism in the paper (Buchdahl, 1979).
In the Palatini formalism, the Ricci scalar in $f({\cal R})$ is ${\cal R} \equiv {\cal R}( g,\Gamma) \equiv g^{\alpha\beta} {\cal R}_{\alpha \beta}(\Gamma )$ and is a generalized Ricci scalar, whereas $ {\cal R}_{\mu \nu}(\Gamma )$ is the Ricci tensor of a torsion-free connection $\Gamma^{\alpha}_{\mu\nu}$ which, a priori, has no relations with the space-time metric $g_{\mu\nu}$ which gives the Ricci scalar $R$. The gravitational sector of the theory is described by the function $f({\cal R})$, while $\sqrt{-g}$ denotes the usual scalar density. The field equations derived with the Palatini variational principle are \begin{eqnarray} && f^{\prime }({\cal R}) {\cal R}_{(\mu\nu)}(\Gamma)-\frac{ f({\cal R})}{2}\, g_{\mu \nu }=T^{(m)}_{\mu\nu}\,, \tag{50} \\ &&\nonumber\\ && \nabla _{\alpha }^{\Gamma } \left[ \sqrt{-g} \, f' ({\cal R} )g^{\mu \nu } \right] =0\,, \tag{51} \end{eqnarray} where $\nabla^{\Gamma}_{\mu} $ is the covariant derivative of the non-metric connection $\Gamma^{\alpha}_{\mu\nu}$.
It is important to stress that Eq.(51) is obtained under the assumption that the matter sector described by ${\cal L}^{(m)}$ is functionally independent of the (non-metric) connection $\Gamma^{\alpha}_{\mu\nu}$; however it may contain metric covariant derivatives $\stackrel{g}{\nabla}$ of the matter fields. This means that the matter stress-energy tensor $ T^{(m)}_{\mu\nu} \left[ g,\Psi \right]$ depends on the metric $g_{\mu\nu}$ and on the matter fields collectively denoted by $\Psi$, together with their covariant derivatives with respect to the Levi-Civita connection of $g_{\mu\nu}$. It is easy to see, from Eq.(51), that $\sqrt{-g} \, f' ({\cal R} ) g^{\mu \nu }$ is a symmetric tensor density of weight $1$, which naturally leads to the introduction of a new metric ${\hat g}_{\mu \nu}$ conformally related to $g_{\mu\nu}$ by (Allemandi, Capone, Capozziello and Francaviglia, 2006) (Ferraris, Francaviglia and Volovich, 1994)
\begin{equation}\tag{52} \sqrt{-g} \, f' ({\cal R}) \, g^{\mu \nu} = \sqrt{-{\hat g}} \,\, {\hat g}^{\mu \nu }\,. \end{equation} With this definition $\Gamma^{\alpha}_{\mu\nu} $ is the Levi-Civita connection of the metric ${\hat g}_{\mu\nu}$, with the only restriction that the conformal factor $\sqrt{-g} \, f' ({\cal R}) g^{\mu \nu}$ relating $g_{\mu\nu}$ and ${\hat g}_{\mu\nu}$ be non-degenerate. In the case of the Hilbert-Einstein Lagrangian, it is $f' ({\cal R})=1$ and the statement is trivial.
The conformal transformation \begin{equation}\tag{53} g_{\mu\nu} \longrightarrow {\hat g}_{\mu \nu }=f' ({\cal R}) \, g_{\mu \nu } \end{equation} implies that ${\cal R}_{(\mu\nu)}(\Gamma)={\cal R}_{\mu \nu}({\hat g})$. It is useful to consider the trace of the field equations (50) \begin{equation} \tag{54} f' ({\cal R}) {\cal R} -2f( {\cal R})=g^{\alpha\beta} T^{(m)}_{\alpha\beta}\equiv T^{(m)} \,, \end{equation} which controls the solutions of Eq.(51). We refer to this scalar equation as the structural or master equation of space-time. In vacuo and in the presence of conformally invariant matter with $T^{(m)}=0$, this scalar equation admits constant solutions. In these cases, Palatini $f({\cal R})$ gravity reduces to General Relativity with a cosmological constant (Ferraris, Francaviglia and Volovich, 1994).
In the case of interaction with matter fields, the structural equation (53), if explicitly solvable, provides, in principle an expression ${\cal R}=f( T^{(m)})$ and, as a result, both $f({\cal R})$ and $f'({\cal R})$ can be expressed in terms of $T^{(m)}$. This fact allows one to express, at least formally, ${\cal R}$ in terms of $T^{(m)}$, which has deep consequences for the description of physical systems. Matter rules the bi-metric structure of space-time and, consequently, both the geodesic and metric structures which are intrinsically different.
$f(R)$ gravity can be achieved also in Hamiltonian formulation in some different approaches (Capozziello and Garattini, 2007). For example, the Ostrogradsky treatment of higher derivative models can be considered or the conformal equivalence between $f(R)$ gravity and Einstein’s General Relativity coupled to a scalar field, as we will see below. However, all these formulations are related by canonical transformations as shown in (Deruelle, Sendouda, and Youssef 2009).
In particular, in the framework of the Arnowitt-Deser-Misner ($\mathcal{ADM}$) formalism (Arnowitt, Deser and Misner, 1963), one can investigate the possibility to find out cosmological terms as eigenvalues of $f(R)$-Hamiltonians in a Sturm-Liouville-like problem (Garattini, 2006). This issue is particularly relevant from several viewpoints. First of all, the vacuum energy of gravitational field is not a particular feature of General Relativity where the cosmological constant has to by added by hand into dynamics. At a classical level, it is well known that $f(R)$ gravity, for the Ricci scalar $R$ equal to a constant, exhibit several de Sitter solutions (Barrow and Ottewill,1983 ) but a definite discussion, at a fundamental level, considering cosmological terms as eigenvalues of such theories is necessary. Besides, the computation of the Casimir energy, the seeking for zero point energy in different backgrounds (for different $f(R)$, we expect different zero point energies and, obviously, different vacuum states), give a track to achieve one-loop energy regularization and renormalization for this kind of theories (Elizalde, Odintsov, Romeo, Bytsenko, Zerbini 1994). On the other hand, these issues can be considered in a multigravity approach to spacetime foam if the $\mathcal{N}$ spacetimes constituting the foam are supposed to evolve, in general, with different curvature laws and ground states (cosmological constants) (Garattini, 2002).
Let us briefly report how to compute the Hamiltonian constraint for General Relativity considering the standard Hilbert-Einstein theory $f\left(R\right)=R$ and the Arnowitt-Deser-Misner ($\mathcal{ADM}$) $3+1$ decomposition (Arnowitt, Deser and Misner, 1963). In terms of these variables, the line element is \[ ds^{2}=g_{\mu\nu}\left( x\right) dx^{\mu}dx^{\nu}=\left( -N^{2}+N_{i} N^{i}\right) dt^{2}+2N_{j}dtdx^{j}+g_{ij}dx^{i}dx^{j}.\tag{55} \] $N$ is the lapse function, while $N_{i}$ the shift functions. In this section we adopt the signature $(-,+,+,+)$ and $\kappa =8\pi c^{-4}$ to point out the role of coupling in the various contributions. In terms of these variables, the gravitational Lagrangian, with the boundary terms neglected, can be written as \begin{equation} \mathcal{L}\left[ N,N_{i},g_{ij}\right] =\sqrt{-g}R=\frac{N{}\,\sqrt{^{3}g} }{2\kappa}\text{ }\left[ K_{ij}K^{ij}-K^{2}+\,\left( ^{3}R-2\Lambda _{c}\right) \right],\tag{56} \end{equation} where $K_{ij}$ is the second fundamental form, $K=$ $g^{ij}K_{ij}$ is the trace, $^{3}R$ is the three dimensional scalar curvature and $\sqrt{^{3}g}$ is the three dimensional determinant of the metric. The conjugate momentum is simply \begin{equation} \pi^{ij}=\frac{\delta\mathcal{L}}{\delta\left( \partial_{t}g_{ij}\right) }=\left( ^{3}g^{ij}K-K^{ij}\text{ }\right) \frac{\sqrt{^{3}g}}{2\kappa }.\tag{57} \end{equation} By a Legendre transformation, we calculate the Hamiltonian \begin{equation} H={\displaystyle\int} d^{3}x\left[ N\mathcal{H+}N_{i}\mathcal{H}^{i}\right], \tag{58} \end{equation} where \begin{equation} \mathcal{H}=\left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}-\frac{\sqrt {^{3}g}}{2\kappa}\left( ^{3}R-2\Lambda_{c}\right) \tag{59} \end{equation} and \begin{equation} \mathcal{H}^{i}=-2\nabla_{j}\pi^{ji}.\tag{60} \end{equation} where $\Lambda_{c}$ is the bare cosmological constant. The equations of motion lead to two classical constraints \begin{equation}\tag{61} \left\{ \begin{array} [c]{c} \mathcal{H}=0\\ \mathcal{H}^{i}=0 \end{array} \right. \end{equation} representing invariance under time re-parameterization and invariance under diffeomorphism, respectively. $G_{ijkl}$ is the supermetric defined as \begin{equation} G_{ijkl}=\frac{1}{2\sqrt{-g}}(g_{ik}g_{jl}+g_{il}g_{jk}-g_{ij}g_{kl}).\tag{62} \end{equation} When $\mathcal{H}$ and $\mathcal{H}^{i}$ are considered as operators acting on some wave function, we have \begin{equation} \mathcal{H}\Psi\left[ g_{ij}\right] =0\tag{63} \end{equation} and \begin{equation} \mathcal{H}_{i}\Psi\left[ g_{ij}\right] =0.\tag{64} \end{equation} Eq. (63) is the Wheeler-DeWitt equation (WDW) (De Witt, 1967). Eqs. (63) and (64) describe the wave function of the universe $\Psi\left[ g_{ij}\right] $. The WDW equation represents invariance under time reparameterization in an operatorial form. This standard lore can be applied to a generic $f(R)$ theory of gravity with the aim to achieve a cosmological term as an eigenvalue of the WDW equation.
Let us consider now the Lagrangian density describing a generic $f(R)$ model, namely \begin{eqnarray}\tag{65} \mathcal{L}=\sqrt{-g}\left( f\left( R\right) -2\Lambda_{c}\right) ,\qquad\text{with}\;f^{\prime\prime}\neq 0. \end{eqnarray} A cosmological term is added also in this case for the sake of generality. Obviously $f^{\prime\prime}=0$ corresponds to General Relativity. The generalized Hamiltonian density for the $f\left( R\right)$ theory assumes the form \begin{equation}\tag{66} \mathcal{H}=\frac{1}{2\kappa}\left[\frac{\mathcal{P}}{6}\left( {}^{\left( 3\right) }R-2\Lambda_{c}-3K_{ij}K^{ij}+K^{2}\right) +V(\mathcal{P})-\frac {1}{3}g^{ij}\mathcal{P}_{\mid ij}-2p^{ij}K_{ij}\right], \end{equation} where \begin{equation} V(\mathcal{P})=\sqrt{-g}\left[ Rf^{\prime}\left( R\right) -f\left( R\right) \right].\tag{67} \end{equation} Henceforth, the superscript $3$ indicating the spatial part of the metric will be omitted on the metric itself. When $f\left( R\right) =R$, $V(\mathcal{P} )=0$ as it should be. Since \begin{equation}\tag{68} \mathcal{P}^{ij}=-2\sqrt{-g}g^{ij}f^{\prime}\left( R\right) \qquad \Longrightarrow\qquad\mathcal{P=}-6\sqrt{-g}f^{\prime}\left( R\right), \end{equation} we have \begin{equation} \mathcal{H}=\frac{1}{2\kappa}\left[ -\sqrt{-g}f^{\prime}\left( R\right) \left( {}^{\left( 3\right) }R-2\Lambda_{c}-3K_{ij}K^{ij}+K^{2}\right) +V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}-2p^{ij}K_{ij}\right].\tag{69} \end{equation} With the help of Eq. (57), Eq.(83) becomes \begin{equation}\tag{70} \mathcal{H}=f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +\frac{1}{2\kappa}\left[ \sqrt {-g}f^{\prime}\left( R\right) \left( 2K_{ij}K^{ij}\right) +V(\mathcal{P} )+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right)_{\mid ij}-2p^{ij}K_{ij}\right]. \end{equation} However \begin{equation}\tag{71} p^{ij}=\sqrt{-g}K^{ij}, \end{equation} then we obtain \begin{equation}\tag{72} \mathcal{H}=f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +\frac{1}{2\kappa}\left[ 2\sqrt {-g}K_{ij}K^{ij}\left( f^{\prime}\left( R\right) -1\right) +V(\mathcal{P} )+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right)\right)_{\mid ij}\right] \end{equation} and transforming into canonical momenta, one gets \begin{equation} \mathcal{H}=f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +2\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4}^{2}\right] \left( f^{\prime}\left( R\right) -1\right) +\frac{1}{2\kappa}\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right)_{\mid ij}\right].\tag{73} \end{equation} By imposing the Hamiltonian constraint, we obtain \begin{eqnarray}\tag{74} f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right)}R\right]+2\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \left( f^{\prime}\left( R\right) -1\right) +\frac{1}{2\kappa }\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}\right] =-f^{\prime}\left( R\right) \sqrt{-g} \frac{\Lambda_{c}}{\kappa} \end{eqnarray} If we assume that $f^{\prime}\left( R\right) \neq0$ the previous expression becomes \begin{equation}\tag{75} \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g} }{2\kappa}{}^{\left( 3\right) }R\right] +\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4}^{2}\right] \frac{2\left( f^{\prime }\left( R\right) -1\right) }{f^{\prime}\left( R\right) }+\frac{1}{2\kappa f^{\prime}\left( R\right) }\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt {-g}f^{\prime}\left( R\right) \right) _{\mid ij}\right]=-\sqrt{-g} \frac{\Lambda_{c}}{\kappa}{}. \end{equation} Now, we integrate over the hypersurface $\Sigma$ to obtain \begin{eqnarray}\tag{76} \int_{\Sigma}d^{3}x\left\{\left[\left(2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right)}R\right] +\left( 2\kappa\right)\left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \frac{2\left(f^{\prime}\left( R\right)-1\right)} {f^{\prime}\left( R\right)}\right\} +\int_{\Sigma}d^{3}x\frac{1}{2\kappa f^{\prime}\left( R\right) }\left[ V(\mathcal{P})+2g^{ij}\left(\sqrt{-g}f^{\prime}\left( R\right)\right) _{\mid ij}\right] =-\frac{\Lambda_{c}}{\kappa}\int_{\Sigma}d^{3}x\sqrt{-g}. \end{eqnarray} The term \begin{equation}\tag{77} \frac{1}{\kappa}\int_{\Sigma}d^{3}x\frac{1}{f^{\prime}\left( R\right) }g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right)\right)_{\mid ij} \end{equation} appears to be a three-divergence and therefore will not contribute to the computation. The remaining equation simplifies into \begin{equation} \int_{\Sigma}d^{3}x\left\{ \left[ \left( 2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right) }R\right] +\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \frac{2\left( f^{\prime}\left( R\right) -1\right) } {f^{\prime}\left( R\right) }+\frac{V(\mathcal{P})}{2\kappa f^{\prime}\left( R\right) }\right\} =-\frac{\Lambda_{c}}{\kappa}\int_{\Sigma}d^{3}x\sqrt {g}.\tag{78} \end{equation} By a canonical procedure of quantization, one can obtain the vacuum state of a generic $f(R)$ theory.
Let us consider now the Lagrangian density describing a generic $f(R)$ theory of gravity, namely \begin{equation} \mathcal{L}=\sqrt{-g}\left( f\left( R\right) -2\Lambda_{c}\right) ,\qquad\text{with}\;f^{\prime\prime}\neq0.\tag{79} \end{equation} A cosmological term is added also in this case for the sake of generality. The generalized Hamiltonian density for the $f\left(R\right)$ theory assumes the form \begin{equation}\tag{80} \mathcal{H}=\frac{1}{2\kappa}\left[ \frac{\mathcal{P}}{6}\left( {}^{\left( 3\right) }R-2\Lambda_{c}-3K_{ij}K^{ij}+K^{2}\right) +V(\mathcal{P})-\frac {1}{3}g^{ij}\mathcal{P}_{\mid ij}-2p^{ij}K_{ij}\right], \end{equation} where \begin{equation} V(\mathcal{P})=\sqrt{-g}\left[ Rf^{\prime}\left( R\right) -f\left( R\right) \right].\tag{81} \end{equation} Henceforth, the superscript $3$ indicating the spatial part of the metric will be omitted on the metric itself. When $f\left( R\right) =R$, $V(\mathcal{P} )=0$ as it should be. Since \begin{equation}\tag{82} \mathcal{P}^{ij}=-2\sqrt{-g}g^{ij}f^{\prime}\left( R\right) \qquad \Longrightarrow\qquad\mathcal{P=}-6\sqrt{-g}f^{\prime}\left( R\right), \end{equation} we have \begin{equation} \mathcal{H}=\frac{1}{2\kappa}\left[ -\sqrt{-g}f^{\prime}\left( R\right) \left( {}^{\left( 3\right) }R-2\Lambda_{c}-3K_{ij}K^{ij}+K^{2}\right) +V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}-2p^{ij}K_{ij}\right].\tag{83} \end{equation} With the help of Eq. (57), Eq. (83) becomes \begin{equation}\tag{84} \mathcal{H}=f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +\frac{1}{2\kappa}\left[ \sqrt {g}f^{\prime}\left( R\right) \left( 2K_{ij}K^{ij}\right) +V(\mathcal{P} )+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}-2p^{ij}K_{ij}\right]. \end{equation} However \begin{equation}\tag{85} p^{ij}=\sqrt{-g}K^{ij}, \end{equation} then we obtain \begin{equation}\tag{86} \mathcal{H}=f^{\prime}\left( R\right) \left[ \left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +\frac{1}{2\kappa}\left[ 2\sqrt {g}K_{ij}K^{ij}\left( f^{\prime}\left( R\right) -1\right) +V(\mathcal{P} )+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right)\right) _{\mid ij}\right] \end{equation} and transforming into canonical momenta, one gets \begin{eqnarray} \mathcal{H}=f^{\prime}\left( R\right) \left[\left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}\left( ^{\left( 3\right) }R-2\Lambda_{c}\right) \right] +2\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4}^{2}\right] \left( f^{\prime}\left( R\right) -1\right) +\frac{1}{2\kappa}\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}\right].\tag{87} \end{eqnarray} By imposing the Hamiltonian constraint, we obtain \begin{eqnarray}\tag{88} f^{\prime}\left(R\right)\left[\left(2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right) }R\right]+2\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \left( f^{\prime}\left( R\right)-1\right)+\frac{1}{2\kappa }\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left(R\right) \right) _{\mid ij}\right] =-f^{\prime}\left( R\right)\sqrt{-g} \frac{\Lambda_{c}}{\kappa} \end{eqnarray} If we assume that $f^{\prime}\left( R\right)\neq0$ the previous expression becomes \begin{equation}\tag{89} \left[\left( 2\kappa\right) G_{ijkl}\pi^{ij}\pi^{kl}{}-\frac{\sqrt{-g} }{2\kappa}{}^{\left( 3\right) }R\right] +\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4}^{2}\right] \frac{2\left( f^{\prime }\left( R\right) -1\right) }{f^{\prime}\left( R\right) }+\frac{1}{2\kappa f^{\prime}\left( R\right) }\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt {g}f^{\prime}\left( R\right)\right)_{\mid ij}\right] =-\sqrt{-g} \frac{\Lambda_{c}}{\kappa}{}. \end{equation} Now, we integrate over the hypersurface $\Sigma$ to obtain \begin{eqnarray}\tag{90} \int_{\Sigma}d^{3}x\left\{ \left[ \left( 2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right) }R\right] +\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \frac{2\left( f^{\prime}\left( R\right) -1\right)} {f^{\prime}\left( R\right)}\right\}+\int_{\Sigma}d^{3}x\frac{1}{2\kappa f^{\prime}\left( R\right) }\left[ V(\mathcal{P})+2g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij}\right]=-\frac{\Lambda_{c}}{\kappa}\int_{\Sigma}d^{3}x\sqrt{-g}. \end{eqnarray} The term \begin{equation}\tag{91} \frac{1}{\kappa}\int_{\Sigma}d^{3}x\frac{1}{f^{\prime}\left( R\right) }g^{ij}\left( \sqrt{-g}f^{\prime}\left( R\right) \right) _{\mid ij} \end{equation} appears to be a three-divergence and therefore will not contribute to the computation. The remaining equation simplifies into \begin{equation} \int_{\Sigma}d^{3}x\left\{\left[\left( 2\kappa\right) G_{ijkl}\pi ^{ij}\pi^{kl}{}-\frac{\sqrt{-g}}{2\kappa}{}^{\left( 3\right) }R\right] +\left( 2\kappa\right) \left[ G_{ijkl}\pi^{ij}\pi^{kl}+\frac{\pi}{4} ^{2}\right] \frac{2\left( f^{\prime}\left( R\right) -1\right)} {f^{\prime}\left( R\right) }+\frac{V(\mathcal{P})}{2\kappa f^{\prime}\left( R\right) }\right\} =-\frac{\Lambda_{c}}{\kappa}\int_{\Sigma}d^{3}x\sqrt {-g}.\tag{92} \end{equation} Starting from these considerations, by a canonical procedure of quantization, one obtains the vacuum state of a generic $f(R)$ theory (Capozziello and Garattini, 2007). In this context, the cosmological constant (vacuum state) emerges as a Wheeler-De Witt eigenvalue. The related wave functional can be split by an orthogonal decomposition and then, constructing the transverse-traceless propagator, it is possible to obtain, after a variational minimization, the total one-loop energy density for the TT tensors. Such a quantity explicitly depends on the form of $f(R)$. The one-loop energy regularization and renormalization are achieved by the zeta function regularization method. The resulting renormalized $\Lambda_{0}$ is a running constant which can be set to zero depending on the value of an arbitrary mass scale parameter $\mu$.
Conformal transformations are a useful tool to deal with Extended Theories of Gravity. Let the pair $\{{\cal M}, g_{\mu\nu}\}$ be a space-time, with ${\cal M}$ a smooth manifold of dimension $n \geq 2$ and $g_{\mu\nu} $ a Lorentzian metric on ${\cal M}$. The point-dependent rescaling of the metric tensor is \begin{eqnarray} \tag{93} g_{\mu\nu} \longrightarrow \tilde{g}_{\mu\nu}=\Omega^2 g_{\mu\nu} \,, \end{eqnarray} where the conformal factor $\Omega(x) $ is a nowhere vanishing, regular function, is called a Weyl or conformal transformation (see (Bronnikov, 2001) (Bronnikov, 2002) (Bronnikov and Shikin, 2002)). Due to this metric rescaling, the lengths of spacelike and timelike intervals and the norms of spacelike and timelike vectors are changed, while null vectors and null intervals of the metric $g_{\mu\nu}$ remain null in the rescaled metric $\tilde{g}_{\mu\nu}$. The light cones are left unchanged by the transformation (93) and the space-times $\{{\cal M}, g_{\mu\nu}\}$ and $\{{\cal M}, \tilde{g}_{\mu\nu}\} $ exhibit the same causal structure (Wald, 1984). A vector that is timelike, spacelike, or null with respect to the metric $g_{\mu\nu}$ has the same character with respect to $\tilde{g}_{\mu\nu}$, and vice-versa.
The transformation properties of geometrical quantities are (Synge, 1955) (Wald, 1984).
\begin{eqnarray} \tag{94} \tilde{g}^{\mu\nu}=\Omega^{-2}\, g^{\mu\nu} \,, \;\;\;\;\;\;\;\;\;\; \tilde{g}=\Omega^{2n} \, g \,, \end{eqnarray} for the inverse metric and the metric determinant, \begin{eqnarray} \tag{95} \tilde{\Gamma}^{\alpha}_{\beta \gamma}= \Gamma^{\alpha}_{\beta\gamma}+\Omega^{-1}\left( \delta^{\alpha}_{\beta} \nabla_{\gamma}\Omega + \delta^{\alpha}_{\gamma} \nabla_{\beta} \Omega -g_{\beta\gamma}\nabla^{\alpha} \Omega \right) \,, \end{eqnarray} for the Christoffel symbols, \begin{eqnarray} \tilde{R}_{\alpha\beta\gamma}^{\delta} &=& {R_{\alpha \beta \gamma}}^{\delta}+2 \, \delta^{\delta}_{[\alpha} \nabla_{\beta]}\nabla_{\gamma} ( \ln \Omega ) -2g^{\delta \sigma} g_{\gamma [ \alpha}\nabla_{\beta ]}\nabla_{\sigma} ( \ln \Omega )+ 2 \nabla_{[ \alpha} ( \ln \Omega ) \, \delta^{\delta}_{\beta ]} \nabla_{\gamma}( \ln \Omega )+\nonumber\\&& -2\nabla_{[ \alpha}( \ln \Omega ) \, g_{\beta ] \gamma} \, g^{\delta \sigma} \, \nabla_{\sigma}( \ln \Omega ) - 2g_{\gamma [ \alpha} \delta^{\delta}_{\beta ]} \, g^{\sigma \rho } \nabla_{\sigma} ( \ln \Omega ) \, \nabla_{\rho} ( \ln \Omega ) \,, \tag{96} \end{eqnarray} for the Riemann tensor, \begin{eqnarray} {\tilde R}_{\alpha\beta } & =& R_{\alpha\beta } -(n-2) \nabla_{\alpha}\nabla_{\beta } ( \ln \Omega ) -g_{\alpha \beta } g^{\rho\sigma } \, \nabla_{\sigma} \nabla_{\rho} ( \ln \Omega )+\left( n-2 \right) \nabla_{\alpha} ( \ln \Omega ) \nabla_{\beta}( \ln \Omega ) -\left( n-2 \right) g_{\alpha\beta }\, g^{\rho\sigma} \, \nabla_{\rho}( \ln \Omega ) \nabla_{\sigma}( \ln \Omega ) \,, \tag{97} \end{eqnarray} for the Ricci tensor, and \begin{eqnarray} \tilde{R} & \equiv & \tilde{g}^{\alpha\beta} \tilde{R}_{\alpha\beta }= \Omega^{-2} \left[ R-2 \left( n-1 \right) \Box \left( \ln \Omega \right) - \left( n-1 \right) \left( n-2 \right) \, \frac{g^{\alpha\beta} \nabla_{\alpha} \Omega \, \nabla_{\beta} \Omega}{\Omega^2} \right] \,, \tag{98} \end{eqnarray} for the Ricci scalar. In the case of $n=4$ space-time dimensions, the transformation property of the Ricci scalar can be written as \begin{eqnarray} \tilde{R} &=& \Omega^{-2} \left[ R-\frac{6 \Box \Omega}{\Omega} \right] = \Omega^{-2} \left[ R-\frac{12 \Box ( \sqrt{\Omega})}{\sqrt{\Omega}} + \frac{3g^{\alpha\beta} \nabla_{\alpha} \Omega \nabla_{\beta} \Omega} {\Omega^2}\right]\,. \tag{99} \end{eqnarray} The Weyl tensor ${C_{\alpha\beta\gamma}}^{\delta}$, with the last index controvariant, is conformally invariant \begin{eqnarray} \tag{100} {\tilde C}_{\alpha\beta\gamma}^{\delta}= {C_{\alpha\beta\gamma}}^{\delta}\,, \end{eqnarray} but the same tensor with indices raised or lowered with respect to ${C_{\alpha\beta\gamma}}^{\delta}$ is not. This property explains the name conformal tensor used for ${C_{\alpha\beta\gamma}}^{\delta}$ (Lorentz, 1937). If the original metric $g_{\alpha\beta} $ is Ricci-flat i.e., $R_{\alpha\beta}=0$), the conformally transformed metric $\tilde{g}_{\alpha\beta}$ is not given by Eq.(97). In the conformally transformed world, the conformal factor $\Omega$ plays the role of an effective matter and this fact has consequences for the physical interpretation of the theory. A vacuum metric in the Jordan frame is not such in the Einstein frame, and the interpretation of what is matter and what is gravity becomes frame-dependent (Sotiriou, Liberati, and Faraoni, 2008). However, if the Weyl tensor vanishes in one frame, it also vanishes in the conformally related frame. Conformally flat metrics are mapped into conformally flat metrics. Since, in general, tensorial quantities are not invariant under conformal transformations, neither are the tensorial equations describing geometry and physics. An equation involving a tensor field $\psi$ is said to be conformally invariant if there exists a number $w$ (the conformal weight of $\psi$) such that, if $\psi$ is a solution of a tensor equation with the metric $g_{\mu\nu}$ and the associated geometrical quantities, $\tilde{\psi} \equiv \Omega^w \psi$ is a solution of the corresponding equation with the metric $\tilde{g}_{\mu\nu}$ and the associated geometry. In addition to geometric quantities, one needs to consider the behavior of standard forms of matter under conformal transformations. Most forms of matter and fields are not conformally invariant: invariance under conformal transformations is a very special property. In general, the covariant conservation equation for a (symmetric) stress-energy tensor $T_{\alpha\beta}^{(m)} $ representing ordinary matter, \begin{eqnarray} \tag{101} \nabla^{\beta} \, T_{\alpha\beta}^{(m)} =0\,, \end{eqnarray} is not conformally invariant (Wald, 1984). The conformally transformed $\tilde{T}_{\alpha\beta}^{(m)} $ satisfies the equation \begin{eqnarray} \tag{102} \tilde{\nabla}^{\beta} \, \tilde{T}_{\alpha\beta}^{(m)} = - \tilde{T}^{(m)} \,\tilde{\nabla}_{\alpha} \left( \ln \Omega \right)\,. \end{eqnarray} Clearly, the conservation equation (101) is conformally invariant only for a matter component that has vanishing trace $T^{(m)}$ of the energy-momentum tensor. This feature is associated with light-like behavior; examples are the electromagnetic field and a radiative fluid with equation of state $p^{(m)}=\rho^{(m)}/3$. Unless $T^{(m)}=0$, Eq.(102) describes an exchange of energy and momentum between matter and the scalar field $\Omega$, reflecting the fact that matter and the geometric factor $\Omega$ are directly coupled in the Einstein frame description. Since the geodesic equation ruling the motion of free particles in General Relativity can be derived from the conservation equation (101) (geodesic hypothesis), it follows that timelike geodesics of the original metric $g_{\alpha\beta}$ are not geodesics of the rescaled metric $\tilde{g}_{\alpha\beta}$ and vice-versa. Particles in free fall in the world $\left\{{\cal M}, g_{\alpha\beta} \right\}$ are subject to a force proportional to the gradient $\tilde{\nabla}^{\alpha} \Omega$ in the rescaled world $\left\{ {\cal M}, \tilde{g}_{\alpha\beta} \right\}$ (this is often identified as a fifth force acting on all massive particles and, therefore, it can be said that no massive test particles exist in the Einstein frame). The stress-energy tensor definition in terms of the matter action ${\cal A}^{(m)}=\int d^4x\, \sqrt{-g} \, {\cal L}^{(m)} $ is \begin{eqnarray} \tag{103} \tilde{T}_{\alpha\beta}^{(m)}=\frac{-2}{\sqrt{ -\tilde{g} } } \, \frac{ \delta \left( \sqrt{-\tilde{g}} \, \, {\cal L}^{(m)} \right) }{\delta \tilde{g}^{\alpha\beta} } \,, \end{eqnarray} together with the rescaling (93) of the metric, yields \begin{equation}\tag{104} \tilde{T}_{\alpha\beta}^{(m)}= \Omega^{-2}\, T_{\alpha\beta}^{(m)}\,,\quad {\tilde{T}_{\alpha}^{\beta}}^{(m)} =\Omega^{-4}\,{T_{\alpha}^{\beta}}^{(m)}\,,\quad {\tilde{T}^{\alpha\beta}}^{(m)}=\Omega^{-6}\, {T^{\alpha\beta}}^{(m)}\,,\quad \tilde{T}^{(m)}=\Omega^{-4}\, T^{(m)}\,. \end{equation} The last relation makes it clear that the trace vanishes in the Einstein frame if and only if it vanishes in the Jordan frame. Conformal invariance corresponds to the absence of a characteristic length (or mass) scale in the physics. In general, the effective potential $V(\phi)$, coming from conformal transformations, contains dimensional parameters (such as a mass $m$, that is a further characteristic gravitational length). This means that the further degrees of freedom coming from Extended Theories of Gravity give rise to features that could play a fundamental role in the dynamics of astrophysical structures addressing, in principle, the dark matter problem.
Cosmology is one of the straightforward applications that motivates the introduction of $f(R)$ gravity. As we will discuss below, $f(R)$ gravity could constitute a geometrical explanation for the dark side of the universe (Capozziello, 2002). The main ingredient of such an interpretation lies on the fact that the further degrees of freedom of $f(R)$ gravity can be assembled into an effective curvature stress-energy tensor (Capozziello, De Laurentis and Lambiase, 2012) which could give rise to dark energy effects.
Let us start from the action (2) where we add also matter. We have \begin{equation}\tag{105} {\cal A}=\int d^4x \sqrt{-g} \left[f(R)+{\cal L}_{(m)} \right]\,{.} \end{equation}
The field are
\begin{equation}\tag{106} f'(R)R_{\alpha\beta}-\frac{1}{2}f(R)g_{\alpha\beta}= f'(R)^{;\alpha\beta}(g_{\alpha\mu}g_{\beta\nu}-g_{\alpha\beta}g_{\mu\nu})+ \tilde{T}^{(m)}_{\alpha\beta}\,, \end{equation} which can be recast in the more expressive form \begin{equation}\tag{107} G_{\alpha\beta}=R_{\alpha\beta}-\frac{1}{2}g_{\alpha\beta}R=T^{(curv)}_{\alpha\beta}+T^{(m)}_{\alpha\beta}\,, \end{equation} where \begin{equation} \tag{108} T^{(curv)}_{\alpha\beta}=\frac{1}{f'(R)}\left\{\frac{1}{2}g_{\alpha\beta}\left[f(R)-Rf'(R)\right]+ f'(R)^{;\alpha\beta}(g_{\alpha\mu}g_{\beta\nu}-g_{\alpha\beta}g_{\mu\nu}) \right\} \end{equation} and \begin{equation}\tag{109} T^{(m)}_{\alpha\beta}=\frac{1}{f'(R)}\tilde{T}^{(m)}_{\alpha\beta}\,, \end{equation} is the stress-energy tensor of matter where we have taken into account the nontrivial coupling to geometry.
Reducing the action to a point-like, Friedmann-Robertson-Walker one, we can write for the geometric part \begin{equation}\tag{110} {\cal A}_{(curv)}=\int dt {\cal L}(a, \dot{a}, R, \dot{R})\,{,} \end{equation} where dot means derivative with respect to the cosmic time. The scale factor $a$ and the Ricci scalar $R$ can be considered as canonical variables. This position could seem arbitrary since $R$ depends on $a, \dot{a}, \ddot{a}$, but it is generally used in canonical quantization. The definition of $R$ in terms of $a, \dot{a}, \ddot{a}$ introduces a constraint which eliminates second and higher order derivatives in action (110), and gives a system of second order differential equations in $\{a, R\}$. Action (110) can be written as \begin{equation}\tag{111} {\cal A}_{(curv)}=2\pi^2\int dt \left\{ a^3f(R)-\lambda\left [ R+6\left ( \frac{\ddot{a}}{a}+\frac{\dot{a}^2}{a^2}+\frac{k}{a^2}\right)\right]\right\}\,{,} \end{equation} where the Lagrange multiplier $\lambda$ is derived by varying with respect to $R$. It is \begin{equation}\tag{112} \lambda=a^3f'(R)\,{.} \end{equation} The point-like Lagrangian is then \begin{equation} \tag{113} {\cal L}={\cal L}_{(curv)}+{\cal L}_{(m)}=a^3\left[f(R)-R f'(R)\right]+6a\dot{a}^2f'(R)+6a^2\dot{a}\dot{R}f''(R)-6ka f'(R)+a^3p_{(m)}\,, \end{equation} where we have taken into account also the fluid matter contribution which is a pressure term.
The Euler-Lagrange equations are \begin{equation} \tag{114} 2\left(\frac{\ddot{a}}{a}\right)+\left(\frac{\dot{a}}{a}\right)^2+ \frac{k}{a^2}=-p_{(tot)}\,, \end{equation} and \begin{equation} \tag{115} f''(R)\left[R+6\left(\frac{\ddot{a}}{a}+\frac{\dot{a}^2}{a}^2+\frac{k}{a^2}\right)\right]=0\,. \end{equation} The dynamical system is completed by the energy condition \begin{equation} \tag{116} \left(\frac{\dot{a}}{a}\right)^2+\frac{k}{a^2}=\frac{1}{3}\rho_{(tot)}\,. \end{equation}
Combining Eq.(114) and Eq.(116), we obtain the Friedmann equation \begin{equation} \tag{117} \left(\frac{\ddot{a}}{a}\right)=-\frac{1}{6}\left[\rho_{(tot)}+3p_{(tot)} \right]\,, \end{equation} where it is clear that the accelerated or decelerated behaviour depends on the r.h.s.. It is \begin{equation} \tag{118} p_{(tot)}=p_{(curv)}+p_{(m)}\;\;\;\;\;\rho_{(tot)}=\rho_{(curv)}+\rho_{(m)}\,, \end{equation} where we have put in evidence the curvature and matter contributions.
From the curvature-stress-energy tensor, we can define a curvature pressure \begin{equation} \tag{119} p_{(curv)}=\frac{1}{f'(R)}\left\{2\left(\frac{\dot{a}}{a}\right)\dot{R}f''(R)+\ddot{R}f''(R)+\dot{R}^2f'''(R) -\frac{1}{2}\left[f(R)-Rf'(R)\right] \right\}\,, \end{equation} and a curvature density \begin{equation} \tag{120} \rho_{(curv)}=\frac{1}{f'(R)}\left\{\frac{1}{2}\left[f(R)-Rf'(R)\right] -3\left(\frac{\dot{a}}{a}\right)\dot{R}f''(R) \right\}\,. \end{equation}
From Eq. (117), the accelerated behaviour is achieved if \begin{equation} \tag{121} \rho_{(tot)}+ 3p_{(tot)}< 0\,, \end{equation} which means \begin{equation} \tag{122} \rho_{(curv)}> \frac{1}{3}\rho_{(tot)}\,, \end{equation} assuming that all matter components have non-negative pressure. In other words, conditions to obtain acceleration depends on the relation \begin{equation} \tag{123} \rho_{(curv)}+3p_{(curv)}=\frac{3}{f'(R)}\left\{\dot{R}^2f'''(R)+\left(\frac{\dot{a}}{a}\right)\dot{R}f''(R) +\ddot{R}f''(R)-\frac{1}{3}\left[f(R)-Rf'(R)\right]\right\}\,, \end{equation} which has to be compared with matter contribution. However, it has to be \begin{equation} \tag{124} \frac{p_{(curv)}}{\rho_{(curv)}}=w_{(curv)}\,, \qquad -1\leq w_{(curv)}<0\,. \end{equation} The form of $f(R)$ is the main ingredient to obtain curvature quintessence.
As simple choice in order to fit the above prescriptions, we ask for solutions of the form \begin{equation} \tag{125} f(R)=f_0 R^n\,,\qquad a(t)=a_0\left(\frac{t}{t_0}\right)^{\beta}\,. \end{equation} The interesting cases are for $n\neq 1$ and $\beta\geq 1$ (accelerated behaviour). Inserting Eqs. (125) into the above dynamical system, we obtain the exact solutions \begin{equation} \tag{126} \beta=2\,;\qquad n=-1,\frac{3}{2}\,;\qquad k=0\,. \end{equation} In both cases, the deceleration parameter is \begin{equation} \tag{127} q_0=-\frac{1}{2}\,, \end{equation} in agreement with the cosmic accelerated behavior.
The case $n=3/2$ deserves a further discussion. Considering conformal transformation from Jordan frame to Einstein frame, it is possible to give explicit form for the scalar field potential. It is \begin{equation}\tag{128} \tilde{g}_{\alpha\beta}\equiv f'(R)g_{\alpha\beta}\,{,}\qquad \varphi=\sqrt{\frac{3}{2}}\ln f'(R)\,{.} \end{equation} The conformal equivalence of the Lagrangians gives \begin{equation}\tag{129} {\cal L}=\sqrt{-g}\,f_0R^{3/2}\longleftrightarrow \tilde{\cal L}=\sqrt{-\tilde{g}}\left[-\frac{\tilde{R}}{2}+ \frac{1}{2}\nabla_{\mu}\varphi\nabla^{\mu}\varphi-V_0\exp\left( \sqrt{\frac{2}{3}}\varphi\right)\right]\,{,} \end{equation} in our physical units. This is the so-called Liouville field theory and it is one of the few cases where a fourth-order Lagrangian can be expressed, in the Einstein frame, in terms of elementary functions under a conformal transformation. It is possible to obtain the general cosmological solution
\begin{equation}\tag{130} a(t)=a_0\sqrt{c_4t^4+c_3t^3+c_2t^2+c_1t+c_0}\,. \end{equation}
The constants $c_i$ are combinations of the initial conditions. Their values determine the cosmological evolution. For example, $c_4\neq 0$ gives a power law inflation while, if the regime is dominated by the linear term in $c_1$, we get a radiation-dominated stage. The above one is just a simple model showing how the dark energy issue can be recovered into the framework of $f(R)$ gravity. However more realistic models can be worked out as reported in literature (Capozziello and Francaviglia,2008) (Nojiri and Odintsov,2007)(Starobinsky,1979). In particular, in (Nojiri and Odintsov, 2003), it has been proposed the possibility of unifying the early-time inflation with late-time acceleration including quintessence. A detailed review in this topic is given in (Nojiri and Odintsov, 2011). It is also worth noticing that, as for other dark energy models, $f(R)$ gravity may show up the finite-time future singular behaviour of one of the known four types as classified in (Nojiri, Odintsov and Tsujikawa, 2005). Several works discussed the finite-time future singularities in F(R) gravity (for review, see (Nojiri and Odintsov, 2011)). Furthermore (Bamba, Nojiri and Odintsov, 2008) have demonstrated that the removal of future singularities in $f(R)$ gravity is easily achieved by the addition of $R^2$ term. Simultaneously with such removal, this term straightforwardly leads to the unification of inflation with dark energy. In other words, a universe model which consistently contains inflation and dark energy, as predicted by $f(R)$ gravity, does not present future singularities.
The phase space portrait of cosmological models, deduced from $f(R)$ gravity, can be discussed by analytical and numerical methods (see (Capozziello, Occhionero and Amendola, 1993) (Carloni, Dunsby, Capozziello and Troisi,2005 ) (Carloni and Dunsby, 2009) (Abdelwahab, Goswami and Dunsby,2012 ). According to this approach, peculiar structures, such as attractors and singular points, emerge. Trajectories of interest are those undergoing an inflationary/dark energy expansion and then reaching Friedmannian asymptotic stable phases. These features can be moreover discussed through effective potentials deduced starting from the form of $f(R)$ gravity (see in particular (Capozziello, Occhionero and Amendola, 1993)). Different kinds of potential regions of phase space can be recognized distinguishing various physical behaviors. There are, in general, allowed regions where trajectories reach the Friedman phases after undergoing an inflationary period, disconnected regions where trajectories, although physical, never reach the Friedman stages and forbidden regions where there are no physical solutions. A survey of the global phase space can be given via Poincaré projections where suitable variables are adopted. This picture allows to represent singularities and attractors at finite and at infinite. The phase-space view can be considered one of the most powerful approaches to deal with dynamical systems in cosmology in order to extract physical information from solutions (see e.g. (Wainwright and Ellis, 1997)).
Cosmological perturbations can be considered also in the framework of $f(R)$ gravity. The evolution of the scalar perturbations since the radiation epoch until the present time can be compared with the evolution of perturbations in the concordance model, i.e., the $\Lambda$CDM model. This method allows to compare the effects due to $f(R)$ modifications in the matter power spectrum as measured today. Specifically, while cosmography is an approach to test the background universe, the matter power spectrum allows to test the clumpy universe. The matter power spectrum provides information on clustering processes coming from the early radiation up to the current epoch. Furthermore, it is a tool to discriminate between General Relativity and modified gravity at the background level. Here, we will briefly sketch the scalar perturbations for metric $f(R)$ gravity giving the main quantities that have to be considered in the perturbation analysis.
The scalar part of cosmological perturbations can be written in the so called Newtonian gauge assuming the conformal invariance, that is (Mukhanov, Feldeman and Branderberger, 1992) (Dodelson, 2003) \begin{eqnarray}\tag{131} ds^2 = a^2\left[\left(1+2\Phi\right)d\eta^2 - \left(1-2\Psi\right)\delta_{ij}dx^idx^j\right], \end{eqnarray} where $\Phi\left(\eta,x^i\right)$ and $\Psi\left(\eta,x^i\right)$ are the gauge invariant Bardeen potentials (Bardeen, 1980) and $\eta$ the conformal time. At first order, the energy-momentum tensor consists of the quantities (Mukhanov, Feldeman and Branderberger, 1992) (Dodelson, 2003) \begin{eqnarray}\tag{132} T^0_0 &= (\rho+\delta\rho), \qquad T^0_i &= -(\rho+p)v_i, \qquad T^i_j &=-(p+\delta p)\delta^i_j, \end{eqnarray} where $\delta\rho$ is the energy density perturbation, $\delta p$ is the pressure perturbation and $v$ is the peculiar velocity. In the simplest approach, we discard the effects of anisotropic stress terms. In the metric $f(R)$ gravity, the $(i,j)$ components relate the potentials $\Phi$ and $\Psi$ with the perturbation of $f'(R)$ (Pogosian and Silvestri, 2007) \begin{align}\tag{133} \left(\Phi-\Psi\right)f'(R) = \delta f'(R) = f''(R)\delta R. \end{align} It is important to stress that the equality between the two potentials, valid in General Relativity, is no longer true in $f(R)$ gravity; the derivative $ f'(R)$ can be interpreted as a new degree of freedom in the background evolution. This means that $\delta f'(R)$ is a new degree of freedom for perturbations. Let us define the new variables $\Psi^+$ and $\Xi$, \begin{eqnarray}\tag{134} \Psi^+ &\equiv \frac{\Phi+\Psi}{2}, \qquad \Xi &\equiv \frac{\delta f'(R)}{ f'(R)} = \Phi-\Psi. \end{eqnarray} $\Psi^+$ corresponds to $\Phi^+$ in Ref. (Pogosian and Silvestri, 2007), while $\Xi$ corresponds to the variable $\chi$ of the same reference divided by $ f'(R)$. We choose the variable $\Xi$ instead of $\chi$, because it represents the difference between the two potentials and allows us to control the integration when the theory deviates from General Relativity ($ f'(R)\gg 1$). Using these variables, we can write the $(0,0)$ and $(i,0)$ components of the perturbed Einstein equations as (de la Cruz-Dombriz, Dobado and Maroto, 2008) \begin{eqnarray} \tag{135} &&\left[1+\frac{1}{2}\frac{\left(f'(R)\right)_N}{ f'(R)}\right]\left(\Psi^+\right)_N +\left[1+\frac{\left( f'(R)\right)_N}{ f'(R)} + \frac{k^2}{3\mathcal{H}^2}\right]\Psi^+ +\frac{1}{4}\frac{\left( f'(R)\right)_N}{ f'(R)}\left(\Xi\right)_N -\frac{1}{2}\left[1 - \mathcal{H}_N +2\frac{\left( f'(R)\right)_N}{ f'(R)}\right]\Xi=-\frac{a^2\kappa^2}{6\mathcal{H}^2 f'(R)}\delta\rho, \end{eqnarray} \begin{eqnarray}\tag{136} \left(\Psi^+\right)_N + \left[1+\frac{1}{2}\frac{\left( f'(R)\right)_N}{ f'(R)}\right]\Psi^+ -\frac{3}{4}\frac{\left( f'(R)\right)_N}{ f'(R)}\Xi &= - \frac{a^2\kappa^2\left(\rho+p\right)}{2 f'(R)}\frac{v}{\mathcal{H}}. \end{eqnarray} The $N$-subscript indicates the derivative with respect to $N\equiv\log(a)$. ${\mathcal{H}}$ is the Hubble parameter written in the conformal time $\eta$. The above equations can be combined to obtain the evolution equations for $\Psi^+ $ and $\Xi$, that is \begin{align} \tag{137} \left(\Psi^+ \right)_N &= -\Psi^+ - \frac{1}{4}\frac{\left( f'(R)\right)_N}{ f'(R)}\left(2\Psi^+ - 3\Xi\right) - \frac{a^2\kappa^2(\rho+p)}{2 f'(R)}\frac{v}{\mathcal H}, \\ \tag{138} \left(\Xi\right) _N &= \Xi + \frac{a^2\kappa^2(\rho+p)}{ f'(R)}\frac{v}{\mathcal H} + \frac{1}{2}\frac{\left( f'(R)\right)_N}{ f'(R)}\left(2\Psi^+ - 3\Xi\right) - \frac{2}{3\mathcal H^2}\frac{f_R}{\left( f'(R)\right)_N}\left[2k^2\Psi^+ + \frac{a^2\kappa^2}{ f'(R)}\left(\delta\rho-3\mathcal{H}(\rho+p)v\right)\right] +\frac{2}{\mathcal{H}}\frac{ f'(R)}{\left( f'(R)\right)_N}\left(\mathcal H - \mathcal H_N\right)\Xi. \end{align} Eqs. (137) and (138) are equivalent to those obtained in Ref. (Pogosian and Silvestri, 2007). The evolution equation for $\Psi^+$ (137) can be used to eliminate the dependence of $(\Xi)_N$ on $(\Psi^+)_N$ in Eq.~(138). This is, finally, the evolution equation for the metric perturbations.
We need now the evolution equations for perturbed matter quantities that are $\delta\rho$ and $v$. For $L$ perfect fluids, each with an energy-momentum tensor as (132), we can define the total, $\delta$, and individual, $\delta^{(i)}$, relative density perturbation, respectively, as \begin{eqnarray}\tag{139} \delta = \sum_{i=1}^L\frac{\rho^{(i)}}{\rho}\delta^{(i)}, \qquad \delta^{(i)} \equiv \frac{\delta\rho^{(i)}}{\rho^{(i)}}, \end{eqnarray} while the total and individual peculiar velocities are related by \begin{align}\tag{140} v &= \sum_{i=1}^L \frac{\rho^{(i)}+p^{(i)}}{\rho+p}v^{(i)}. \end{align} By perturbing the conservation equations of the energy momentum tensor, we find that, for adiabatic and non-interacting fluids, each pair $\delta^{(i)}$ and $v^{(i)}$ satisfies the equations (Mukhanov, Feldeman and Branderberger, 1992) \begin{align} \tag{141} &\left(\delta^{(i)}\right)_N + 3\left(c_s^{(i)2}-w^{(i)}\right) - (1+w^{(i)})k^2\frac{v^{(i)}}{\mathcal H}= 3(1+w^{(i)})\left(\Psi^+ + \frac{1}{2}\Xi\right)_N, \\ \tag{142} &\left(v^{(i)}\right)_N + \left(1-3c_s^{(i)2} \right)v^{(i)} + \frac{c_s^{(i)2}}{1+w^{(i)}}\frac{\delta^{(i)}}{\mathcal H} = - \frac{1}{\mathcal H}\left(\Psi^+ -\frac{1}{2} \Xi\right). \end{align} Here, $w^{(i)}\equiv p^{(i)}/\rho^{(i)}$ and $c_s^{(i)2}\equiv p^{(i)}_N/\rho^{(i)}_N$ are, respectively, the state parameter and the squared speed of sound of the $(i)$-fluid.
Eqs. (137), (138), (141), and (142) form a closed set of equations that allows the evolution of perturbations since the early radiation epoch until the present time. Once the present day value of the perturbation variables is computed, they can be related with observable quantities.
However also other approaches are possible for perturbation theory in $f(R)$ gravity. For example, the evolution of density perturbations can be dealt with the Covariant and Gauge Invariant theory. In (Carloni, Dunsby and Troisi, 2008), $f(R)\sim R^n$ power-law models are are dealt with general perturbations equations obtaining exact solutions for scales much bigger than the Hubble radius. These solutions have a some interesting features. In particular, all values of $n$ show growing modes for the density contrast, even if the universe undergoes an accelerated expansion. Such a behavior does not occur in General Relativity, where as soon as dark energy dominates, the density contrast experiences a decay. This peculiarity is a particular feature of fourth order gravity models.
The evolution of scalar perturbations in$f(R)$ gravity can be tracked since the radiation epoch until the present time, and compared with the evolution of perturbations in the $\Lambda$CDM model. In particular, the modifications in the matter power spectrum due to $f(R)$ gravity can be compared, in principle, with today measured cosmological parameters.
While cosmography is one of the several approaches to test the background universe, the matter power spectrum is a tool to test the clumpy universe. Such a spectrum provides information about the process of clustering of matter in the universe at all scales, i.e., from those corresponding to the early radiation dominated epoch up to those which have recently reentered the horizon. In addition, it is a formidable tool to break the possible degeneracy between General Relativity and modified gravity at the background level. In order to obtain the matter power spectrum, we need the evolution of the scalar perturbations as deduced above. In other words, the perturbative analysis carried above allows to see how the $f(R)$ corrections affect the matter power spectrum $\mathcal{P}_{\delta_{m}}$. Notice however that the correct definition of $\mathcal{P}_{\delta_{m}}$ uses the energy density perturbation $\delta_m$ in the comoving gauge (for a discussion on this point see e.g. Refs. (Wands and Slozar,2009) (Bruni, Crittenden, Koyama, Maartens, Pitrou and Wands,2011), (Clifton, Ferreira,Padilla and Skordis2012)), while the variables used here correspond to the Newton gauge. In order to account for this gauge difference, the matter power spectrum can be calculated as \begin{equation}\tag{143} \mathcal{P}_{\delta_{m}} \equiv \langle|\delta_{m}^{(com)}|^2\rangle = \langle|\delta_{m} - 3\mathcal{H}v_{m}|^2\rangle. \end{equation}
As an example, let as consider a cosmological model filled with dust and an effective energy density given by $f(R)$ gravity. In principle, the $f(R)$ model can account for the radiation content of the universe and drives its late time acceleration. The initial conditions can be chosen as \begin{align} \tag{144} \Psi^+(N_{ini})=1, ~~~~~~ \left(\Psi^+\right)_N(N_{ini})=0, ~~~~~~ \Xi(N_{ini})=0, ~~~~~~ \left(\Xi\right)_N(N_{ini})=0, \end{align} to match the initial conditions in the $\Lambda$CDM model with a radiation component and a single inflationary field, i.e., the field $\Psi^+$ is initially constant and the modes are well outside of the horizon. During the radiation dominated epoch, the $f(R)$ corrections, in general, increase the value of the metric perturbations by several orders of magnitude, in contrast with what happens in General Relativity. Furthermore, from the above analysis, the initial value of the density perturbation $\delta$ are related to $\Psi^+(N_{ini})$ as \begin{equation} \tag{145} \delta_m(N_{ini}) = - \frac{6\mathcal{H}^2f_R}{a^2\kappa^2\rho}\left[1+\frac{\left(f_R\right)_N}{f_R} + \frac{k^2}{3\mathcal{H}^2}\right]\Psi^+ (N_{ini})\,. \end{equation} Initially one can have $(f_R)_N/f_R\ll1$, furthermore, given that all the relevant modes are outside of the horizon, we have that $k^2/(3\mathcal{H}^2)\ll1$. Therefore, the relation between $\delta$ and $\Psi^+$ is ruled by the factor $6\mathcal{H}^2f_R/a^2\kappa^2\rho$. In General Relativity, the Friedmann equation gives $3\mathcal{H}^2=a^2\kappa^2\rho$ ($f_R=1)$ and the matter and metric perturbations have the same order of magnitude. However, in $f(R)$ gravity, if the matter energy density is initially sub-dominant with respect to the $f(R)$ effective energy density, the factor on the right hand side of Eq.~(145) can become large, and so does $\delta$. Considering the General Relativity limit on Eq.~(145) we obtain $\delta_m=-2\Psi^+$. This relation refers to the total matter energy density perturbation $\delta$. For a General Relativity model, beside matter, we have also radiation, therefore we need two initial conditions for $\delta_m$ and $\delta_r$. Those are obtained by imposing the adiabaticity conditions, see for example (Dodelson, 2003), ( Amendola and Tsujikawa, 2010) \begin{equation}\tag{146} \frac{\delta_m}{1+w_m} = \frac{\delta_r}{1+w_r} = \frac{\delta}{1+w}, \end{equation} where $w_m=0$, $w_r=1/3$, and initially $w\approx1/3$. Therefore, in a General Relativity model with matter, radiation and initial adiabatic perturbations, we have the initial relation between the matter and the metric perturbations \begin{align}\tag{147} \delta_m = \frac{3}{4}\delta_r=\frac{3}{4}\delta = -\frac{3}{2}\Psi^+. \end{align} The presence of radiation and the adiabaticity of the initial perturbations thus corrects the General Relativity limit of Eq.~(145) by a factor of $3/4$.
In $f(R)$ gravity, the initial value of $\delta_m$ can be orders of magnitude different than the one in General Relativity due to the effective stress-energy tensor. As a result, the matter power spectrum can be very different from the one of the $\Lambda$CDM model. Besides the overall increase/decrease in the amplitude of the spectrum, the $f(R)$ effects change the shape of the spectrum, in particular on the high $k$ regime. This is consistent with the scale-dependent effects of $f(R)$ on the matter perturbation (Pogosian and Silvestri, 2007) (de la Cruz-Dombriz, Dobado and Maroto, 2008).
In summary, one can study the evolution of the first order metric and matter perturbations since the radiation dominated epoch until the present time and obtain the theoretical power spectrum of the matter perturbation. The numerical methods to achieve results can be based, for example, on the so called "quasi static approximation" (de la Cruz-Dombriz, Dobado and Maroto, 2008). In general, the effects of the $f(R)$ corrections on the matter perturbation are stronger on the modes with higher wave-numbers, as is expected from the corrections to the closed evolution equation of $\delta_m$ (de la Cruz-Dombriz, Dobado and Maroto, 2008). More precisely, this can be seen directly on the matter power spectra obtained: for low $k$ end of the spectrum, the shape of the spectra obtained is consistent with the one of the $\Lambda$CDM model (up to a multiplicative factor); for modes with $k> k_{eq}\approx0.01$ Mpc$^{-1}$ the shape of the spectra obtained starts to diverge from the one in the $\Lambda$CDM model. Notice however that the matter power spectrum can be obtained assuming that the initial amplitude profile of the initial conditions is that of $\Phi^+$ for single field inflation, that is \begin{equation}\tag{148} \mathcal P _{\Psi^+} = \langle |\Psi^+|^2\rangle = \frac{8\pi^2}{9k^3}A_s\left(\frac{k}{k_0}\right)^{n_s-1}, \end{equation} which could be modified due to $f(R)$ effects. Here, $A_s$ is the curvature $\mathcal R$ power spectrum at some pivot scale $k_0$ ; $n_s$ is the scalar spectrum index that can be obtained by observations (e.g. the PLANCK mission (Ade et al.,2014). The relation $\mathcal{R}=(3/2)\Psi^+$ can be used during the radiation epoch and for modes outside of the horizon. As a final remark, Large Scale Structure and perturbation analysis vs observations can put, in principle, strong constraints on the form and the evolution of any $f(R)$ gravity model.
The above considerations can be worked out in view of addressing the problem of the dark side of the universe. Dark matter issues come from dynamical mass estimates of self-gravitating systems. In several astrophysical observations, there is more matter dynamically inferred than that can be accounted for from luminous components. This mass discrepancy is usually attributed to additional (missing or dark) matter, assuming the validity of Newton law of gravity at astrophysical scales. Oort was the first that posed the missing matter problem (Oort, 1932) (Oort, 1960). By observing the Doppler red-shift values of stars moving near the plane of our Galaxy, he asserted that he could calculate how fast the stars were moving. He found that there had to be enough matter inside the Galaxy such that the central gravitational force was strong enough to keep the stars from escaping, much as Sun's gravitational pull keeps a planet in its orbit. But when the calculation was made, it turned out that there was not enough mass in the Galaxy. The discrepancy was not small: the Galaxy had to be at least twice as massive as the sum of the mass of all its visible components combined. In addition, in the 1960's the radial profile of the tangential velocity of stars in their orbits around the Galactic Center, as a function of their distance from that center, was measured. It was found that typically, once we get away from the Galactic Center, all the stars travel with the same velocity independent of their distance out from the Galactic Center.
There were problems, too, at a larger scale. In 1933, Zwicky announced that when he measured the individual velocities of a large group of galaxies, known as the Coma Cluster, he found that all of the galaxies that he considered were moving so rapidly relative to one another that the cluster should have come apart long ago. The visible mass of the galaxies making up the cluster was far too little to produce enough gravitational force to hold the cluster together. So not only our own Galaxy was lacking mass, but also the whole Coma Cluster of galaxies was suffering the same problem at a different scale (Zwicky, 1933). Initially, the problem was only approached by leaving Newton's law inviolated and postulating the existence of some invisible dark entities to make up the missing mass. At the beginning, it has never came to mind anyone to go back and examine the basic assumption that only gravity was at work in these cases. It was easier to patch up the theory introducing invisible entities. Many names have been coined to define this invisible entity, a bit as in the days of ether.
There are the Massive Compact Halo Objects (MACHOs), objects like black holes, and neutron stars that purportedly populate the outer reaches of galaxies like the Milky Way. Then there are the Weakly Interacting Massive Particles (WIMPs), which possess mass, yet do not interact with ordinary matter (baryons such as protons and neutrons) because they are composed by something unknown. Dark (missing) matter (DM) even comes in two flavors, hot (HDM) and cold (CDM). The CDM is supposedly to be in dead stars, planets, brown dwarfs ("failed stars") etc., while HDM is postulated to be fast moving in particles floating throughout the universe. It should be constituted by neutrinos, tachyons etc. But where is all of this missing matter? The truth is that after many years of looking for it, there is still no definitive proof that WIMPs exist, or that MACHOs will ever make up more than five percent of the total reserve of missing dark stuff.
Besides, by adding a further ingredient, the cosmological constant $\Lambda$, such a model (now $\Lambda$CDM) has become the new cosmological paradigm usually called the concordance model. In fact, high quality data coming from the measurements of cluster properties as the mass, the correlation function and the evolution with redshift of their abundance, the Hubble diagram of Type Ia Supernovae, the optical surveys of large scale structure, the anisotropies of the cosmic microwave background, the cosmic shear measured from weak lensing surveys and the Lyman-$\alpha$ forest absorption are evidences toward a spatially flat universe with a subcritical matter content and undergoing a phase of accelerated expansion. Interpreting all this information in a self-consistent model is the main task of modern cosmology and $\Lambda$CDM model provides a good fit to the most part of the data giving a reliable picture of the today observed universe.
Nevertheless, it is affected by serious theoretical shortcomings that have motivated the search for alternative candidates generically referred to as dark energy or quintessence. Such models range from scalar fields rolling down self interaction potentials to phantom fields, from phenomenological unified models of dark energy and dark matter to alternative gravity theories.
Essentially, dark energy (or any alternative component) has to act as a negative pressure fluid which gives rise to an overall acceleration of the Hubble fluid. Despite of the clear mechanisms generating the observed cosmological dynamics, the nature and the fundamental properties of dark energy remain essentially unknown notwithstanding the great theoretical efforts made up to now.
The situation for dark matter is similar: its clustering and distribution properties are fairly well known at every scale but its nature is unknown, up to now, at fundamental level.
On the other hand, the need of unknown components as dark energy and dark matter could be considered nothing else but as a signal of the breakdown of Einstein General Relativity at astrophysical (galactic and extragalactic) and cosmological scales. In this context, Extended Theories of Gravity could be, in principle, an interesting alternative to explain cosmic acceleration and large scale structure without any missing components. In their simplest version, the Ricci curvature scalar $R$, linear in the Hilbert-Einstein action, could be replaced by a generic function $f(R)$ whose true form could be "reconstructed" by the data. In fact, there is no a priori reason to consider the gravitational Lagrangian linear in the Ricci scalar while observations and experiments could contribute to define and constrain the "true" theory of gravity (Capozziello and Faraoni, 2010) (Capozziello and De Laurentis, 2010).
Coming to the weak-field limit, which essentially means considering Solar System scales, any alternative relativistic theory of gravity is expected to reproduce General Relativity which, in any case, is firmly tested only in this limit and at these scales (Will,1993). Even this limit is a matter of debate since several relativistic theories do not reproduce exactly the Einstein results in their Newtonian limit but, in some sense, generalize them.
In general, all these efforts can be included in the so called Modified Newtonian Dynamics (MOND), first proposed by Milgrom (Milgrom, 1983) (Milgrom and Sanders, 1983), in the attempt to explain missing matter without dark matter but assuming a change into dynamics at scales larger than Solar System's ones. In general, any relativistic theory of gravitation yields corrections to the weak-field gravitational potentials which, at the post-Newtonian level and in the Parametrized Post-Newtonian formalism, could constitute a test of these theories (Will,1993). This point deserves a deep discussion. Beside the fundamental physics motivations coming from Quantum Gravity and unification theories (Capozziello and De Laurentis, 2011), Extended Theories of Gravity pose the problem that there are further gravitational degrees of freedom (related to higher order terms, non-minimal couplings and scalar fields in the field equations) and gravitational interaction is not invariant at any scale. This means that, besides the Schwarzschild radius, other characteristic gravitational scales come out from dynamics. Such scales, in the weak field approximation, should be responsible of characteristic lengths of astrophysical structures that should result confined in this way (Capozziello and De Laurentis,2012).
We can deal with the Newtonian and the post-Newtonian limit of $f(R)$ gravity adopting the spherical symmetry. The solution of field equations can be obtained considering the general spherically symmetric metric: \begin{eqnarray}\tag{149} ds^2\,=\,g_{\sigma\tau}dx^\sigma dx^\tau=g_{00}(x^0,r)d{x^0}^2-g_{rr}(x^0,r)dr^2-r^2d\Omega \end{eqnarray} where $x^0\,=\,ct$ and $d\Omega$ is the solid angle. In order to develop the Newtonian limit, let us consider the perturbed metric with respect to a Minkowskian background $g_{\mu\nu}\,=\,\eta_{\mu\nu}+h_{\mu\nu}$. The metric entries can be developed as: \begin{eqnarray}\tag{150} \left\{\begin{array}{ll} g_{tt}(t, r)\simeq1+g^{(2)}_{tt}(t,r)+g^{(4)}_{tt}(t,r)\,, \\\\ g_{rr}(t,r)\simeq-1+g^{(2)}_{rr}(t,r)\,,\\\\ g_{\theta\theta}(t,r)=-r^2\,,\\\\ g_{\phi\phi}(t,r)=-r^2\sin^2\theta\,, \end{array}\right.\, \end{eqnarray} where we are assuming $c\,=\,1$,$x^0=ct\rightarrow t$ and expansion is at order $c^{-2}$ and $c^{-4}$ . Since we want to obtain the most general result, we does not provide any specific form for $f(R)$ Lagrangian. We assume, however, analytic Taylor expandable $f(R)$ functions with respect to a certain value $R\,=\,R_0$:
\begin{eqnarray}\tag{151} f(R)=\sum_{n}\frac{f^n(R_0)}{n!}(R-R_0)^n\simeq f_0+f_1R+f_2R^2+f_3R^3+...\,, \end{eqnarray} where we have assumed $R_0=0$. In order to obtain the weak field approximation, one has to insert expansions (150) and (151) into field Eqs. (23) and expand the system up to the orders ${\mathcal O}(0)$, ${\mathcal O}(2)$ e ${\mathcal O}(4)$ (that is, as stated above at order $c^{-2}$ and $c^{-4}$). This approach provides general results and specific (analytic) theories are selected by the coefficients $f_i$ in Eq.(151). It is worth noticing that, at the order ${\mathcal O}(0)$, the field equations give the condition $f_0 =0$ and then the solutions at further orders do not depend on this parameter. If we consider the ${\mathcal O}(2)$-order approximation, the field equations in vacuum, results to be \begin{eqnarray}\tag{152} \left\{\begin{array}{ll} f_1rR^{(2)}-2f_1g^{(2)}_{tt,r}+8f_2R^{(2)}_{,r}-f_1rg^{(2)}_{tt,rr}+4f_2rR^{(2)}=0\,, \\\\ f_1rR^{(2)}-2f_1g^{(2)}_{rr,r}+8f_2R^{(2)}_{,r}-f_1rg^{(2)}_{tt,rr}=0\,, \\\\ 2f_1g^{(2)}_{rr}-r\left[f_1rR^{(2)}-f_1g^{(2)}_{tt,r}-f_1g^{(2)}_{rr,r}+4f_2R^{(2)}_{,r}+4f_2rR^{(2)}_{,rr}\right]=0\,, \\\\ f_1rR^{(2)}+6f_2\left[2R^{(2)}_{,r}+rR^{(2)}_{,rr}\right]=0\,, \\\\ 2g^{(2)}_{rr}+r\left[2g^{(2)}_{tt,r}-rR^{(2)}+2g^{(2)}_{rr,r}+rg^{(2)}_{tt,rr}\right]=0\,. \end{array} \right.\end{eqnarray} It is evident that the trace equation (the fourth in the system (152)), provides a differential equation with respect to the Ricci scalar which allows to solve exactly the system (152) at ${\mathcal O}(2)$-order. Finally, one gets the general solution:
\begin{eqnarray}\tag{153} \left\{\begin{array}{ll} g^{(2)}_{tt}=\delta_0-\frac{Y}{f_1r}-\frac{\delta_1(t)e^{-r\sqrt{-\xi}}}{3\xi r}+\frac{\delta_2(t)e^{r\sqrt{-\xi}}}{6({-\xi)}^{3/2}r} \\\\ g^{(2)}_{rr}=-\frac{Y}{f_1r}+\frac{\delta_1(t)[r\sqrt{-\xi}+1]e^{-r\sqrt{-\xi}}}{3\xi r}-\frac{\delta_2(t)[\xi r+\sqrt{-\xi}]e^{r\sqrt{-\xi}}}{6\xi^2r} \\\\ R^{(2)}=\frac{\delta_1(t)e^{-r\sqrt{-\xi}}}{r}-\frac{\delta_2(t)\sqrt{-\xi}e^{r\sqrt{-\xi}}}{2\xi r}\end{array} \right. \end{eqnarray} where $\xi\doteq\displaystyle\frac{f_1}{6f_2}$, and $Y$ is an arbitrary integration constant. When we consider the limit $f(R)\rightarrow R$, in the case of a point-like source of mass $M$, we recover the standard Schwarzschild solution. Let us notice that the integration constant $\delta_0$ is dimensionless, while the two arbitrary functions of time $\delta_1(t)$ and $\delta_2(t)$ have respectively the dimensions of $length^{-1}$ and $length^{-2}$. The functions of time $\delta_i(t)$ ($i=1,2$) are completely arbitrary since the differential equation system (152) contains only spatial derivatives and can be fixed to constant values. Besides, the integration constant $\delta_0$ can be set to zero since it represents an unessential additive quantity for the potential. It is possible to write the general solution of the problem considering the previous expression (149). In order to match at infinity the Minkowskian prescription for the metric, one can discard the Yukawa growing mode and then we have: \begin{eqnarray}\tag{154} \left\{\begin{array}{ll}ds^2\,=\,\biggl[1-\frac{2GM}{f_1r}-\frac{\delta_1(t)e^{-r\sqrt{-\xi}}}{3\xi r}\biggr]dt^2- \biggl[1+\frac{2GM}{f_1r}-\frac{\delta_1(t)(r\sqrt{-\xi}+1)e^{-r\sqrt{-\xi}}}{3\xi r}\biggr]dr^2-r^2d\Omega\\\\R\,=\,\frac{\delta_1(t)e^{-r\sqrt{-\xi}}}{r}\end{array}\right. \end{eqnarray} At this point, one can provide the gravitational potential. The first of (153) gives the second order solution in term of the metric expansion (see the definition (150)). This term coincides with the gravitational potential at the Newton order. In particular, since $g_{tt}\,=\,1+2\Phi_{grav}\,=\,1+g_{tt}^{(2)}$, the gravitational potential of $f(R)$-gravity, analytic in the Ricci scalar $R$, is \begin{eqnarray}\tag{155} \Phi_{grav}\,=\,-\left[\frac{GM}{f_1r}+\frac{\delta_1(t)e^{-r\sqrt{-\xi}}}{6\xi r}\right]\,. \end{eqnarray} This general result means that the standard Newton potential is achieved only in the particular case $f(R)=R$ while it is not so for any analytic $f(R)$ models. Eq.(155) deserves some comments. The parameters $f_{1,2}$ and the function $\delta_1$ represent the deviations with respect the standard Newton potential. To test these theories of gravity inside the Solar System, we need to compare such quantities with respect to the current experiments, or, in other words, Solar System constraints should be evaded fixing such parameters. On the other hand, these parameters could acquire non-trivial values (e.g. $f_1\neq 1,\,\delta_1(t)\neq 0,\,\xi\neq 1$) at scales different from the Solar System ones. We note that the $\xi$ parameter can be related to an effective mass being $ m^2= (3\xi)^{-1} $ and can be interpreted also as an effective length $L$. Eq. (155) can be recast as \begin{equation} \tag{156} \Phi(r) = -\frac{G M}{ (1+\delta) r}\left(1+\delta e^{-\frac{r}{L}}\right)\,, \end{equation} where the first term is the Newtonian-like part of the potential for a point-like mass $\displaystyle{\frac{M}{1+\delta}}$ and the second term is a modification of gravity including a scale length $L$ associated to the above coefficients of the Taylor expansion. If $\delta=0$ the Newtonian potential and the standard gravitational coupling are recovered. Comparing Eqs.(155) and (156), we obtain that $\displaystyle{1+\delta=f_1}$, and $\delta$ is related to $\delta_1(t)$ through \begin{equation}\tag{157} \delta_1=-\frac{6GM}{L^2}\left(\frac{\delta}{1+\delta}\right) \end{equation} where $\displaystyle{\frac{6GM}{L^2}}$ and $\delta_1$ can be assumed quasi-constant. Under these assumptions, the scale length $L$ could naturally reproduce several phenomena that range from Solar System to cosmological scales. Understanding on which scales the modifications to General Relativity are working or what is the weight of corrections to the Newtonian gravitational potential is a crucial point that could confirm or rule out these alternative approaches (Capozziello and De Laurentis, 2012).
Solar System tests for relativistic theories of gravity include gravitational redshift, deflection of light by the Sun, planetary orbit precession at perihelion. General Relativity is consistent with these experimental tests. Considering them gives a useful framework where predictions of different relativistic theories can be parameterized in a systematic way (Capozziello and De Laurentis, 2011). In particular, the Post Newtonian limit framework (Will,1993) has become a basic tool to connect gravitational theories with Solar System experiments. Furthermore, Post Newtonian formalism can be applied to all metric theories of gravitation where self-gravitating bodies satisfy the Einstein Equivalence Principle. In general, it is possible to find constraints on $f(R)$ models (and any Extended Theories of Gravity) from Solar System experiments combined with experiments on the violation of Equivalence Principle (Capozziello and Tsujikawa, 2008).
In an environment of high density such as Earth or Sun, the Ricci scalar $R$ is larger than the cosmological background. If the outside of a spherically symmetric body is the vacuum, the metric can be described by a Schwarzschild-like exterior solution with $R=0$. In the presence of non-relativistic matter with a matter-energy density $\rho_m$, this gives rise to a contribution to the Ricci scalar $R$ of the order $\rho_m$. The task is now to test if the above result involving the Yukawa correction to the potential, is compatible or not with Solar System tests and Equivalence Principle. We will deal with the further degrees of freedom, coming from $f(R)$, that can be dealt under the standard of an additional scalar field coming into dynamics. This approach is allowed by conformal transformations.
When the mass of an equivalent scalar field degree of freedom is heavy in a region with high density, a spherically symmetric body can be represented as a thin-shell, so that an effective coupling due to the fifth force (the Yukawa correction) is suppressed through a Chameleon Mechanism ( Khoury and Weltman, 2004). In this mechanism, the mass of the effective scalar field depends on the density environment. If the matter density is sufficiently high, the field acquires a heavy mass about the potential minimum. On the other hand, the field has lighter mass in a low-density cosmological environment relevant for dark energy. In this case, it can propagate freely. As long as a spherically symmetric body has thin-shell around its surface, the effective coupling between the field and matter becomes much smaller than the coupling (De Felice and Tsujikawa, 2010).
Starting from the action for $f(R)$, introducing the new metric $\tilde{g}_{\mu \nu}$ and a scalar field $\phi$, conformal transformations allow to write
\begin{eqnarray}\tag{158} \tilde{g}_{\mu \nu}=\psi g_{\mu \nu}\,, \quad \phi=\sqrt{\frac{3}{2}}\,{\rm ln}\,\psi\,, \end{eqnarray}
where, clearly, the conformal transformation is specified by $\psi=\partial f/\partial R$. We can rewrite the $f(R)$ action in the Einstein frame as (Maeda, 1989) (Capozziello and Faraoni, 2010) (Capozziello and De Laurentis, 2010), in physical units \begin{eqnarray} \tag{159} {\cal S} &=&\int{\rm d}^{4}x\sqrt{-\tilde{g}}\left[ -\frac{1}{2}\tilde{R}+\frac{1}{2} g^{\mu\nu} \phi_{;\mu}\phi_{;\nu}-V(\phi) \right] +{\cal L}^{(m)} (\tilde{g}_{\mu \nu} e^{2\beta \phi}, \Psi_m)\,, \end{eqnarray}
where
\begin{eqnarray}\tag{160} \beta=-\frac{1}{\sqrt{6}}\,,\quad V(\phi(\psi))=\frac{R(\psi)\psi-f}{2\psi^2}\,. \end{eqnarray}
The field $\phi$ is directly coupled to the non-relativistic standard matter by a coupling $\beta$ that we will specify below. In a spherically symmetric space-time, the variation of the action (159) with respect to the scalar field $\phi$ gives
\begin{eqnarray} \tag{161} \frac{d^2 \phi}{d {\tilde{r}}^2}+ \frac{2}{{\tilde{r}}} \frac{ d\phi}{d {\tilde{r}}}= \frac{d V_{eff}}{d\phi}\,, \end{eqnarray}
where ${\tilde{r}}$ is the distance from the center of symmetry and
\begin{eqnarray} \tag{162} V_{eff}(\phi)=V(\phi)+ e^{\beta \phi}\rho^*\,. \end{eqnarray}
Here $\rho^*$ is a conserved quantity in the Einstein frame. It is related to the matter-energy density $\rho$, in the Jordan frame, via the relation $\rho^*=e^{3\beta \phi}\rho$.
Assuming that a spherically symmetric body has a constant density $\rho^*=\rho_A^*$ inside the body (${\tilde r}<{\tilde {r}}_c$) and that the energy density outside the body (${\tilde r}>{\tilde {r}}_c$) is $\rho^*=\rho_B^*$, the mass $M_c$ of the body and the gravitational potential $\Phi_c$ at the radius ${\tilde r}_c$ are given by $M_c=(4\pi/3){\tilde r}_c^3 \rho_A^*$ and $\Phi_c=M_c/8\pi {\tilde{r}}_c$, respectively. The effective potential $V_{eff} (\phi)$ has two minima at the field values $\phi_A$ and $\phi_B$ satisfying $V_{eff}' (\phi_A)=0$ and $V_{eff}' (\phi_B)=0$, respectively. The former corresponds to the region with a high density that gives rise to a heavy mass squared $m_A^2 \equiv V_{eff}''(\phi_A)$, whereas the latter to the lower density region with a lighter mass squared $m_B^2 \equiv V_{eff}''(\phi_B)$. In the high-density regime (with a heavy field mass), the spherically symmetric body acquires a thin-shell under the Chameleon Mechanism. When the thin-shell develops inside the body, the following thin-shell parameter is much smaller than the order of unity (see (Capozziello and Tsujikawa, 2008) and references therein):
\begin{eqnarray} \tag{163} \frac{\Delta {\tilde{r}}_c}{{\tilde{r}}_c} =\frac{\phi_B-\phi_A}{6\beta \Phi_c}\,. \end{eqnarray}
Solving Eq.(161) with appropriate boundary conditions, the field profile outside the body (${\tilde{r}}>{\tilde{r}}_c$) is given by
\begin{eqnarray} \tag{164} \phi({\tilde{r}}) \simeq -\frac{\beta_{eff}}{4\pi} \frac{M_c e^{-m_B({\tilde{r}}-{\tilde{r}}_c)}}{{\tilde{r}}}+\phi_B\,, \end{eqnarray}
where the magnitude of the effective coupling, $\beta_{eff}= (3\beta) (\Delta {\tilde{r}}_c/{\tilde{r}}_c)$, is much smaller than unity when the thin-shell is formed. The bound on the thin-shell parameter from experimental tests of the Post Newtonian parameter in Solar Systems can be derived. The spherically symmetric metric in Einstein frame is given by
\begin{eqnarray} \tag{165} d{\tilde s}^2=\psi d s^2={\tilde{g}}_{\mu\nu} d{\tilde{x}}^\mu d{\tilde{x}}^\nu=\left[1-2{\cal {\tilde{A}}}({\tilde{r}})\right]dt^2-\left[1+2{\cal {\tilde{B}}}({\tilde r})\right]d{\tilde{r}}^2-{\tilde{r}}^2d\Omega^2\,, \end{eqnarray} where ${\tilde{A}}({\tilde{r}})$ and ${\tilde B}({\tilde r})$ are functions of ${\tilde r}$ and, as above, $d\Omega^2=d\theta^2+(\sin^2 \theta)d\phi^2$. In the weak field limit, it is ${\tilde{A}}({\tilde{r}})\ll 1$ and ${\tilde B}({\tilde{r}})\ll 1$. The metric outside the spherically symmetric body with mass $M_c$ is given by ${\tilde A}({\tilde {r}})\simeq{\tilde B}({\tilde{r}})\simeq \frac{GM_c}{{\tilde {r}}}$. Let us transform back the metric (165) into the Jordan frame. It is \begin{eqnarray}\tag{166} {\rm d}s^2=[1-2{\cal{A}}(r)]{\rm d}t^2 -[1+2{\cal{B}}(r)]{\rm d}r^2-r^2 {\rm d} \Omega^2\,, \end{eqnarray} Comparing the two metrics under the condition $|\phi| \ll 1$, we find that the following relations hold (Faulkner et al, 2007)
\begin{eqnarray}\tag{167} r = \psi^{-1/2}\tilde{r}\,, \quad {\cal{A}}(r) \simeq \tilde{{\cal{A}}}(\tilde{r})+ \frac{\phi(\tilde{r})}{\sqrt{6}}\,, \quad {\cal{B}}(r) \simeq \tilde{{\cal{B}}}(\tilde{r})+ \frac{\tilde{r}}{\sqrt{6}} \frac{{\rm d}\phi(\tilde{r})}{{\rm d} \tilde{r}}\,. \end{eqnarray}
Under the weak-field approximation in the Einstein frame, the metric components $\tilde{{\cal{A}}}(\tilde{r})$ and $\tilde{{\cal{B}}}(\tilde{r})$ outside a spherically symmetric body with mass $M_c$ are given by $\tilde{{\cal{A}}}(\tilde{r}) \simeq \Phi_c=M_c/(8\pi \tilde{r})$. Provided that $|\phi| \ll 1$, one has $\psi \simeq 1$ and hence $\tilde{r} \simeq r$. In the following, we omit a tilde for the quantities in the Einstein frame because the condition $|\phi| \ll 1$ always holds for the analytic $f(R)$ models we are studying. Using the thin-shell solution (164), we find that the metrics ${\cal{A}}(r)$ and ${\cal{B}}(r)$ in the Jordan frame are
\begin{eqnarray}\tag{168} {\cal{A}}(r)=\frac{G_{\rm eff}M_c}{r}\,,\quad {\cal{B}}(r)=\gamma \frac{G_{\rm eff}M_c}{r}\,, \end{eqnarray}
where the effective gravitational coupling $G_{\rm eff}$ and the post-Newtonian parameter $\gamma$ are given by
\begin{eqnarray} \tag{169} G_{\rm eff} &\simeq& G \left[ 1-\frac{\sqrt{6}}{3}\beta_{\rm eff} e^{-m_B (r-r_c)} \right]\,,\\ \\ \tag{170} \gamma &\simeq& \frac{1+(\sqrt{6}\beta_{\rm eff}/3)(1+m_{B} r) e^{-m_B (r-r_c)}}{1-(\sqrt{6}\beta_{\rm eff}/3) e^{-m_B (r-r_c)}}\,. \end{eqnarray}
If the condition $m_Br \ll 1$ is satisfied in an environment where local gravity experiments are carried out, we obtain
\begin{eqnarray}\tag{171} G_{\rm eff} &\simeq& G \left( 1+\frac{\Delta r_c}{r_c} \right)\,, \\ \\ \tag{172} \gamma &\simeq& \frac{1-\Delta r_c/r_c} {1+\Delta r_c/r_c}\,. \end{eqnarray}
As long as the body has a thin-shell ($|\Delta r_c/r_c| \ll 1$), one can satisfy the experimental bounds $|\gamma-1| \ll 1$ (Will,1993). In the thick-shell regime, we just need to change $\beta_{\rm eff}$ in Eqs.(169)-(170) to $\beta=-1/\sqrt{6}$. Under the condition $m_Br \ll 1$, we obtain $G_{\rm eff}=4G/3$ and $\gamma=1/2$, which contradicts with the experimental bound mentioned above. If the mass $m_B$ is heavy so that $m_Br \gg 1$ for $r \gtrsim r_c$, the mass $m_A$ inside the body is typically much heavier to satisfy the relation $m_A r_c \gg 1$. In such cases the body has a thin-shell, which means that the post-Newtonian parameter is given by Eq.(172).
In Ref. (Faulkner et al, 2007), it was shown that under the Chameleon Mechanism the post-Newtonian parameter, $\gamma={\cal B}(r)/{\cal A}(r)$, is approximately given by (172) provided that the condition $m_B r \ll 1$ holds at Solar System scales. The present tightest constraint on $\gamma$ is $|\gamma-1| <2.3 \times 10^{-5}$ (Will,1993), which translates into
\begin{eqnarray} \tag{173} \frac{\Delta r_c}{r_c}<1.15 \times 10^{-5}\,. \end{eqnarray} Let us now place experimental bounds from a possible violation of Equivalence Principle. The thin-shell condition around the Earth under the Chameleon Mechanism has to be considered ( Khoury and Weltman, 2004). The Earth has a radius $r_{\oplus}=6 \times 10^{3}$ km with a mean density $\rho_{\oplus} \simeq 5.5$ g/cm$^3$. The atmosphere exists in the region $r_{\oplus}<r<r_{atm}$ with a homogeneous density $\rho_{atm} \simeq 10^{-3}$ g/cm$^3$. The region outside the atmosphere ($r>r_{atm}$) has a homogenous density $\rho_G \simeq 10^{-24}$ g/cm$^3$. The atmosphere exists in the region $r_{\oplus}<r<r_{atm}$ with a homogeneous density $\rho_{atm} \simeq 10^{-3}$ g/cm$^3$. The region outside the atmosphere ($r>r_{atm}$) has a homogenous density $\rho_G \simeq 10^{-24}$ g/cm$^3$. Defining the gravitational potentials as $\Phi_{\oplus}=\rho_{\oplus}r_{\oplus}^2/6$ and $\Phi_{\rm atm}=\rho_{atm}r_{atm}^2/6$, we have that $\Phi_{\oplus} \simeq 5.5 \times 10^3 \Phi_{\rm atm}$ because $\rho_{\oplus} \simeq 5.5 \times 10^3 \rho_{atm}$ and $r_{\oplus} \simeq r_{\rm atm}$. Recalling the relation $\Delta r_{atm}/r_{atm}= (\phi_G-\phi_{atm})/(6\beta \Phi_{atm})$, where $\phi_G$ and $\phi_{atm}$ correspond to the field values at the local minima of the effective potential (162) in the regions $r>r_{atm}$ and $r_{\oplus}<r<r_{atm}$ respectively, we find $\Delta r_{\oplus}/r_{\oplus} \equiv -(\phi_G-\phi_{atm}) /\sqrt{6} \Phi_{\oplus} \simeq 2.0 \times 10^{-4} (\Delta r_{atm}/r_{atm})$. When the atmosphere has a thin-shell, then the thickness of the shell ($\Delta r_{atm}$) is smaller than that of the atmosphere: $r_s=10$-10$^2$km. Taking the value $r_s=10^2$ km and $r_{atm} =6.5 \times 10^3$ km, we obtain $\Delta r_{\rm atm}/r_{atm}<1.6 \times 10^{-2}$. Hence the condition for the atmosphere to have a thin-shell is estimated as
\begin{eqnarray} \tag{174} \frac{\Delta r_{\oplus}}{r_{\oplus}} \lesssim 10^{-6}\,. \end{eqnarray}
Solar System tests of Equivalence Principle consider, essentially, free-fall acceleration of the Moon and the Earth toward the Sun (see (Capozziello and De Laurentis, 2011) for a discussion on the Equivalence Principle). The constraint on the difference of the two accelerations is given by
\begin{eqnarray} \tag{175} \eta \equiv 2\frac{|a_{\rm Moon}-a_{\oplus}|} {a_{\rm Moon}+a_{\oplus}}<10^{-13}\,. \end{eqnarray}
The Sun and the Moon have the thin-shells like the Earth (Capozziello and Tsujikawa, 2008), where the field profiles outside the spheres are given as in Eq.(164) with the replacement of corresponding quantities. We note that the acceleration induced by a fifth force with the field profile $\phi(r)$ and the effective coupling $\beta_{eff}$ is $a^{fifth}= |\beta_{eff}\phi(r)|$. Then the accelerations $a_{\oplus}$ and $a_{Moon}$ are (De Felice and Tsujikawa, 2010) (Capozziello and Tsujikawa, 2008)
\begin{eqnarray}\tag{176} a_{\oplus} &\simeq& \frac{GM_{\odot}}{r^2} \left[ 1+3 \left(\frac{\Delta r_{\oplus}} {r_{\oplus}}\right)^2 \frac{\Phi_{\oplus}}{\Phi_{\odot}} \right]\,,\\ \\\tag{177} a_{\rm Moon} &\simeq& \frac{GM_{\odot}}{r^2} \left[ 1+3 \left(\frac{\Delta r_{\oplus}} {r_{\oplus}}\right)^2 \frac{\Phi_{\oplus}^2} {\Phi_{\odot}\Phi_{\rm Moon}}\right]\,, \end{eqnarray}
where $\Phi_{\odot} \simeq 2.1 \times 10^{-6}$, $\Phi_{\oplus} \simeq 7.0 \times 10^{-10}$ and $\Phi_{Moon} \simeq 3.1 \times 10^{-11}$ are the gravitational potentials of Sun, Earth and Moon, respectively. Hence the condition (175) translates into
\begin{eqnarray} \tag{178} \frac{\Delta r_{\oplus}}{r_{\oplus}}<2 \times 10^{-6}\,, \end{eqnarray}
which gives the same order of the upper bound as in the thin-shell condition (174) for the atmosphere. The constraint coming from the violation of strong equivalence principle (Will,1993) provides a bound $\Delta r_{\oplus}/r_{\oplus}<10^{-4}$ (Capozziello and Tsujikawa, 2008), which is weaker than (178).
In summary, any correction that one can take into account for the Newtonian potential has to satisfy the above constraints on Equivalence Principle and Solar System. This statement has to hold also for $f(R)$ gravity.
$f(R)$ gravity can be linearized producing a third massive mode of gravitational radiation. Taking the trace of (25) one gets
\begin{equation} 3\square f'(R)+Rf'(R)-2f(R)=0,\tag{179}\end{equation}
and, with the identifications (Capozziello, Corda and De Laurentis, 2008)
\begin{equation} \begin{array}{ccccc} \Phi\rightarrow f'(R) & & \textrm{and } & & \frac{dV}{d\Phi}\rightarrow\frac{2f(R)-Rf'(R)}{3}\end{array}\tag{180}\end{equation}
a Klein-Gordon equation for the effective $\Phi$ scalar field is obtained
\begin{equation} \square\Phi=\frac{dV}{d\Phi}.\tag{181}\end{equation} We assume $\Phi_{0}$ to be a minimum for $V$
\begin{equation} V\simeq\frac{1}{2}\alpha\delta\Phi^{2}\Rightarrow\frac{dV}{d\Phi}\simeq m^{2}\delta\Phi,\tag{182}\end{equation}
and the constant $m$ has the mass dimension.
Assuming
\begin{equation} \begin{array}{c} g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu}\,, \qquad \Phi=\Phi_{0}+\delta\Phi.\end{array}\tag{183}\end{equation}
to first order in $h_{\mu\nu}$ and $\delta\Phi$. Considering $\widetilde{R}_{\mu\nu\rho\sigma}$, $\widetilde{R}_{\mu\nu}$ and $\widetilde{R}$ as the linearized quantities that correspond to $R_{\mu\nu\rho\sigma}$, $R_{\mu\nu}$ and $R$, the linearized field equations are
\begin{equation}\tag{184} \begin{array}{c} \widetilde{R}_{\mu\nu}-\frac{\widetilde{R}}{2}\eta_{\mu\nu}=(\partial_{\mu}\partial_{\nu}h_{f}-\eta_{\mu\nu}\square h_{f})\\ \\{}\square h_{f}=m^{2}h_{f},\end{array}\end{equation}
where
\begin{equation} h_{f}\equiv\frac{\delta\Phi}{\Phi_{0}}.\tag{185}\end{equation}
$\widetilde{R}_{\mu\nu\rho\sigma}$ and Eqs. (184) are invariants for gauge transformations
\begin{equation}\tag{186} \begin{array}{c} h_{\mu\nu}\rightarrow h'_{\mu\nu}=h_{\mu\nu}-\partial_{(\mu}\epsilon_{\nu)}\,, \quad\delta\Phi\rightarrow\delta\Phi'=\delta\Phi;\end{array}\end{equation}
then
\begin{equation} \bar{h}_{\mu\nu}\equiv h_{\mu\nu}-\frac{h}{2}\eta_{\mu\nu}+\eta_{\mu\nu}h_{f}\tag{187}\end{equation}
can be defined. Considering the transformation for the parameter $\epsilon^{\mu}$
\begin{equation} \square\epsilon_{\nu}=\partial^{\mu}\bar{h}_{\mu\nu},\tag{188}\end{equation} a gauge can be chosen, that is
\begin{equation} \partial^{\mu}\bar{h}_{\mu\nu}=0.\tag{189}\end{equation}
In this way, field equations read like
\begin{equation} \square\bar{h}_{\mu\nu}=0\tag{190}\end{equation}
\begin{equation} \square h_{f}=m^{2}h_{f}\tag{191}\end{equation}
Solutions of Eqs. (190) and (191) are plane waves
\begin{equation} \bar{h}_{\mu\nu}=A_{\mu\nu}(\overrightarrow{p})\exp(ip^{\alpha}x_{\alpha})+c.c.\tag{192}\end{equation}
\begin{equation} h_{f}=a(\overrightarrow{p})\exp(iq^{\alpha}x_{\alpha})+c.c.\tag{193}\end{equation}
where
\begin{equation} \begin{array}{ccc} k^{\alpha}\equiv(\omega,\overrightarrow{p}) & & \omega=p\equiv|\overrightarrow{p}|\\ \\q^{\alpha}\equiv(\omega_{m},\overrightarrow{p}) & & \omega_{m}=\sqrt{m^{2}+p^{2}}.\end{array}\tag{194}\end{equation}
In Eqs. (190) and (192) the equation and the solution for the standard waves of General Relativity (Landau and Lifsits, 1999) (Misner, Thorne and Wheeler, 1973) have been obtained, while Eqs. (191) and (193) are respectively the equation and the solution for the massive mode related to $f(R)$ gravity (see also (Bogdanos, Capozziello, De Laurentis and Nesseris,2010)).
The fact that the dispersion law for the modes of the massive field $h_{f}$ is not linear has to be emphatized. The velocity of any massles mode (arising from General Relativity) $\bar{h}_{\mu\nu}$ is the light speed $c$. On the other hand, the dispersion law (the second of Eq. (194)) for the modes of $h_{f}$ is that of a massive wave-packet. The group-velocity is
\begin{equation} \overrightarrow{v_{G}}=\frac{\overrightarrow{p}}{\omega},\tag{195}\end{equation}
that is the velocity of a massive particle with mass $m$ and momentum $\overrightarrow{p}$.
From Eqs. (194) and Eq. (195), we obtain
\begin{equation} v_{G}=\frac{\sqrt{\omega^{2}-m^{2}}}{\omega}.\tag{196}\end{equation}
The mass is
\begin{equation} m=\sqrt{(1-v_{G}^{2})}\omega.\tag{197}\end{equation}
From Eq. (187) it is
\begin{equation} h_{\mu\nu}=\bar{h}_{\mu\nu}-\frac{\bar{h}}{2}\eta_{\mu\nu}+\eta_{\mu\nu}h_{f}.\tag{198}\end{equation}
Considering the massles case one obtains
\begin{equation} \begin{array}{c} \square\epsilon^{\mu}=0\\ \\\partial_{\mu}\epsilon^{\mu}=-\frac{\bar{h}}{2}+h_{f},\end{array}\tag{199}\end{equation}
which gives the total transversality of the field. In the massive case
\begin{equation} \square\epsilon^{\mu}=m^{2}h_{f}.\tag{200}\end{equation}
In the same way, it is possible to show that it does not exist any linear relation between the tensorial field $\bar{h}_{\mu\nu}$ and the massive field $h_{f}$. After some gauge considerations it is possible to show that the gravitational waves can be written as
\begin{equation} h_{\mu\nu}(t,z)=A^{+}(t-z)e_{\mu\nu}^{(+)}+A^{\times}(t-z)e_{\mu\nu}^{(\times)}+h_{f}(t-v_{G}z)\eta_{\mu\nu}.\tag{201}\end{equation}
The term $A^{+}(t-z)e_{\mu\nu}^{(+)}+A^{\times}(t-z)e_{\mu\nu}^{(\times)}$ describes the two standard polarizations of gravitational waves which arise from General Relativity, while the term $h_{f}(t-v_{G}z)\eta_{\mu\nu}$ is the massive field arising from $f(R$) gravity. In other words, the derivative function $f'(R$) generates a third massive mode for gravitational waves which is not present in standard General Relativity.
$f(R)$ gravity is a class of effective theories representing a new approach to the gravitational interaction. The paradigm is that Einstein's General Relativity has to be extended in order to address several shortcomings emerging at ultraviolet and infrared scales. These are essentially due to the lack of a final, self-consistent theory of quantum gravity. From the astrophysical and cosmological viewpoints, the goal is to encompass phenomena like dark energy and dark matter under a geometric standard related to the possibility that gravitational interaction depends on the scales.
This geometric view, in principle, does not need the introduction of further particle ingredients and preserves all the well-posed results of General Relativity, being based on the same fundamental principles (Equivalence Principle, diffeomorphism invariance, gauge invariance, etc.).
The main criticism to this approach is that, until now, no $f(R)$ model, or any Extended Theory of Gravity, succeeds in addressing the whole phenomenology ranging from quantum to cosmological scales. Besides, the $f(R)$ description of dark side of the universe is substantially equivalent to that related to the hypothesis of dark material constituents. The need of one or more than one experimentum crucis, capable of retaining or ruling out one of the two concurring pictures, is pressing to solve the debate.