Consider Einstein-Hilbert action for gravitational field,
\[ S = \int \mathrm{d}^4 x \sqrt{-g} \left\{ \frac{1}{2\kappa}R - \frac{1}{\kappa}\Lambda + \mathcal{L}_M \right\}\] where $R = R^{\mu\nu} g_{\mu\nu}$ is the curvature scalar, $R^{\mu\nu} = R^\lambda_{\mu\lambda\nu}$ is the Racci tensor, $\kappa = 8 \pi G/c^2$.

Einstein-Hilbert action is the most simple Lorentzian invariant action that encodes space-time curvature and gives correct gravity (Newtonian gravity) in weak field limit. A straightforward generalization of Einstein-Hilbert action is the $f(R)$ gravity: $R \to f(R) = a_0 + a_1 R+ a_2 R^2 + \cdots$. From the aesthetic point of view, one may argue the higher order terms shouldn't be there. But there terms may arise for a good reason. We know that at high energy $\sim T_\text{Pl}$, general relativity (GR) will be superseded by quantum gravity. This in fact may even happen before Planck scale. General relativity is merely a low energy effective theory. Effective theories are not always neat. In fact, it is well known that quantized general relativity is non-renormalizable. Therefore, quantum fluctuation would bring infinite terms into the effective Lagrangian.

A natural question is, if these term has any observable effect. Let's restrict ourselves to classical theory only. The question is, how well has general relativity been tested [1]? $f(R)$ gravity and many other competing theories of gravity shares the same foundations with general relativity, except for the amount of gravity produced by the same energy. The relevant test includes deflection of light ray and time delay (Fig. 1), Mercury perihelion shift and spin precession (Table 4), change of Newton's constant (Table 5), mass of graviton (or propagation distance of gravity) etc.

Note that most of these tests are conducted with in solar system. They cannot exclude gravity theory with only large scale effects. It has been known [2-3] $R^2$ and/or $R^{\mu\nu}R_{\mu\nu}$ term can produce a massive scalar field in addition to graviton. This field couples with matter through Yukawa potential. Its mass can be chosen such than too heavy to propagate a observable distance, but too light $\ll M_\text{Pl}$ to modified gravity. In other words, it is a massive, very weakly interaction particle. Naturally, it is a dark matter candidate, since it only manifests itself in large scales.

Now let's see how it happens. Consider $R^n$ term: \[ S_n = \frac{1}{\kappa}\int \mathrm{d}^4 x \sqrt{-g} R^n \] Now let's take functional derivative with respect to $g^{\mu\nu}$. \[ \frac{\delta \sqrt{-g} R^n}{\delta g^{\mu\nu}} = \sqrt{-g} R^{n-1} \left( n \frac{\delta R}{\delta g^{\mu\nu}} + \frac{R}{\sqrt{-g}}\frac{\delta\sqrt{-g}}{\delta g^{\mu\nu}} \right)

\] Now, according to Jacobi's formular, $\delta g = g g_{\mu\nu} \delta g^{\nu\mu}$. So \[ \frac{1}{\sqrt{-g}}\frac{\delta\sqrt{-g}}{\delta g^{\mu\nu}} = -\frac{1}{2}g_{\mu\nu}.

\] \[ \frac{\delta R}{\delta g^{\mu\nu} } = R_{\mu\nu} + g^{\alpha\beta} \left( \frac{\delta \Gamma_{\alpha\beta;\lambda}^\lambda}{\delta g^{\mu\nu}} - \frac{\delta\Gamma_{\alpha\lambda;\beta}^\lambda}{\delta g^{\mu\nu}} \right) \] Note that $ \delta R^\lambda_{\;\;\alpha\lambda\beta} = \delta\Gamma^\lambda_{\alpha\beta;\lambda} - \delta\Gamma^\lambda_{\alpha\lambda;\beta}$. $\delta \Gamma$'s are in fact tensors: \[ \delta\Gamma^\lambda_{\alpha\beta} = \frac{1}{2} g^{\lambda\rho}\left\{ \delta g_{\rho\alpha;\beta} + \delta g_{\rho\beta;\alpha} - \delta g_{\alpha\beta;\rho} \right\}. \]

Hence \[ \frac{\delta R}{\delta g^{\mu\nu} } = R_{\mu\nu} + \frac{\delta g^{\alpha\beta}_{ \;\;\;\;;\alpha;\beta}}{\delta g^{\mu\nu}} - \frac{\delta g^{\alpha \;\;;\beta}_{\;\; \alpha \;\;;\beta}}{\delta g^{\mu\nu}} \]

Therefore, the action \[ \delta S_n = \frac{1}{\kappa} \int \mathrm{d}^4 x \sqrt{-g} \delta g^{\mu\nu} R^{n-1} \left\{ n R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \frac{n-1}{R} \left(R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) + \frac{(n-1)(n-2)}{R^2} \left( R_{;\mu}R_{;\nu} - g_{\mu\nu} R^{;\alpha} R_{;\alpha} \right) \right\} \] yields the equation of motion: \[ \sum_n a_n R^{n-1} \left\{ n R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \frac{n-1}{R} \left(R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) + \frac{(n-1)(n-2)}{R^2} \left( R_{;\mu}R_{;\nu} - g_{\mu\nu} R^{;\alpha} R_{;\alpha} \right) \right\} = \frac{\kappa}{2} T_{\mu\nu} \]

For $n = 1, a_1 = \frac{1}{2} $ it reduces to GR. One of the most distinguish feature of this formulation is for $n=1$ (GR), two extra terms vanish. If we keep terms up to $n=2$, the equation of motion simplifies as:\[

(1 + 2 a_2 R) G_{\mu\nu} + a_2 g_{\mu\nu} R^2 + 2 a_2 \left( R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) = \kappa T_{\mu\nu}

\]

Now, let's identify the particle contents. Define $h_{\mu\nu} = g_{\mu\nu}-\eta_{\mu\nu}$, where $\eta_{\mu\nu}$ is the Minkowski metric. In weak field limit, we can simply take $h_{\mu\nu} \to \delta g_{\mu\nu}$ and $g_{\mu\nu} \to \eta_{\mu\nu}$. But graviton is usually defined as the inverse trace $\bar{h}_{\mu\nu} = h_{\mu\nu} - \frac{1}{2}\eta_{\mu\nu} h$, where $ h = h^\mu_\mu$. Note that, $G_{\mu\nu} = -\frac{1}{2} \Box \bar{h}_{\mu\nu} $ and $ R = \bar{h}^{\mu\nu}_{\;\;\;\; ,\mu\nu} - \frac{1}{2} \Box h $. The lowest order terms are: \[ S = \frac{1}{\kappa}\int \mathrm{d}^4 x \left\{ -\frac{1}{4} \bar{h}^{\mu\nu} \Box \bar{h}_{\mu\nu} + \frac{1}{8} h ( 2 a_2 \Box^2 + \Box ) h - \frac{a_2}{2} \bar{h}^{\mu\nu} \partial_\mu\partial_\nu \Box h \right\} \]

As we can see, there is an additional field besides the conventional graviton field, with mass $m^2_h = -\frac{1}{2 a_2}$. This field interacts with matter via Yukawa coupling: $\sqrt{-g} \simeq 1 + \delta \sqrt{-g} = 1 + \frac{1}{2} h$. Therefore, it propagates $e^{- r /m_h}$. Taking $R \ll m^2_h \ll M_\text{Pl}$, $R^2$ gravity produces a massive weakly interacting particle.

[1]: Clifford M. Will,

[2]: Jose A. R. Cembranos,

[3]: S. Stelle, Gen. Relativ. Gravit. 9, 353 (1978)

.

Einstein-Hilbert action is the most simple Lorentzian invariant action that encodes space-time curvature and gives correct gravity (Newtonian gravity) in weak field limit. A straightforward generalization of Einstein-Hilbert action is the $f(R)$ gravity: $R \to f(R) = a_0 + a_1 R+ a_2 R^2 + \cdots$. From the aesthetic point of view, one may argue the higher order terms shouldn't be there. But there terms may arise for a good reason. We know that at high energy $\sim T_\text{Pl}$, general relativity (GR) will be superseded by quantum gravity. This in fact may even happen before Planck scale. General relativity is merely a low energy effective theory. Effective theories are not always neat. In fact, it is well known that quantized general relativity is non-renormalizable. Therefore, quantum fluctuation would bring infinite terms into the effective Lagrangian.

A natural question is, if these term has any observable effect. Let's restrict ourselves to classical theory only. The question is, how well has general relativity been tested [1]? $f(R)$ gravity and many other competing theories of gravity shares the same foundations with general relativity, except for the amount of gravity produced by the same energy. The relevant test includes deflection of light ray and time delay (Fig. 1), Mercury perihelion shift and spin precession (Table 4), change of Newton's constant (Table 5), mass of graviton (or propagation distance of gravity) etc.

Fig. 1: test of deflection of light, $ \delta \theta = \frac{1+\gamma}{2}\frac{4 G M}{d} $ and test of time delay experiment $ \delta t = 2(1+\gamma ) GM \ln \left( \frac{(r_\oplus + \mathbf{x}_\oplus\cdot \mathbf{n})(r_e - \mathbf{x}_e\cdot \mathbf{n})}{d^2}\right)$. General relativity predicts $\gamma = 1$. |

Table 5: test of change of Newton constant. |

Now let's see how it happens. Consider $R^n$ term: \[ S_n = \frac{1}{\kappa}\int \mathrm{d}^4 x \sqrt{-g} R^n \] Now let's take functional derivative with respect to $g^{\mu\nu}$. \[ \frac{\delta \sqrt{-g} R^n}{\delta g^{\mu\nu}} = \sqrt{-g} R^{n-1} \left( n \frac{\delta R}{\delta g^{\mu\nu}} + \frac{R}{\sqrt{-g}}\frac{\delta\sqrt{-g}}{\delta g^{\mu\nu}} \right)

\] Now, according to Jacobi's formular, $\delta g = g g_{\mu\nu} \delta g^{\nu\mu}$. So \[ \frac{1}{\sqrt{-g}}\frac{\delta\sqrt{-g}}{\delta g^{\mu\nu}} = -\frac{1}{2}g_{\mu\nu}.

\] \[ \frac{\delta R}{\delta g^{\mu\nu} } = R_{\mu\nu} + g^{\alpha\beta} \left( \frac{\delta \Gamma_{\alpha\beta;\lambda}^\lambda}{\delta g^{\mu\nu}} - \frac{\delta\Gamma_{\alpha\lambda;\beta}^\lambda}{\delta g^{\mu\nu}} \right) \] Note that $ \delta R^\lambda_{\;\;\alpha\lambda\beta} = \delta\Gamma^\lambda_{\alpha\beta;\lambda} - \delta\Gamma^\lambda_{\alpha\lambda;\beta}$. $\delta \Gamma$'s are in fact tensors: \[ \delta\Gamma^\lambda_{\alpha\beta} = \frac{1}{2} g^{\lambda\rho}\left\{ \delta g_{\rho\alpha;\beta} + \delta g_{\rho\beta;\alpha} - \delta g_{\alpha\beta;\rho} \right\}. \]

Hence \[ \frac{\delta R}{\delta g^{\mu\nu} } = R_{\mu\nu} + \frac{\delta g^{\alpha\beta}_{ \;\;\;\;;\alpha;\beta}}{\delta g^{\mu\nu}} - \frac{\delta g^{\alpha \;\;;\beta}_{\;\; \alpha \;\;;\beta}}{\delta g^{\mu\nu}} \]

Therefore, the action \[ \delta S_n = \frac{1}{\kappa} \int \mathrm{d}^4 x \sqrt{-g} \delta g^{\mu\nu} R^{n-1} \left\{ n R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \frac{n-1}{R} \left(R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) + \frac{(n-1)(n-2)}{R^2} \left( R_{;\mu}R_{;\nu} - g_{\mu\nu} R^{;\alpha} R_{;\alpha} \right) \right\} \] yields the equation of motion: \[ \sum_n a_n R^{n-1} \left\{ n R_{\mu\nu} - \frac{1}{2} g_{\mu\nu} R + \frac{n-1}{R} \left(R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) + \frac{(n-1)(n-2)}{R^2} \left( R_{;\mu}R_{;\nu} - g_{\mu\nu} R^{;\alpha} R_{;\alpha} \right) \right\} = \frac{\kappa}{2} T_{\mu\nu} \]

For $n = 1, a_1 = \frac{1}{2} $ it reduces to GR. One of the most distinguish feature of this formulation is for $n=1$ (GR), two extra terms vanish. If we keep terms up to $n=2$, the equation of motion simplifies as:\[

(1 + 2 a_2 R) G_{\mu\nu} + a_2 g_{\mu\nu} R^2 + 2 a_2 \left( R_{;\mu;\nu} - g_{\mu\nu} R^{;\alpha}_{\;\;;\alpha} \right) = \kappa T_{\mu\nu}

\]

Now, let's identify the particle contents. Define $h_{\mu\nu} = g_{\mu\nu}-\eta_{\mu\nu}$, where $\eta_{\mu\nu}$ is the Minkowski metric. In weak field limit, we can simply take $h_{\mu\nu} \to \delta g_{\mu\nu}$ and $g_{\mu\nu} \to \eta_{\mu\nu}$. But graviton is usually defined as the inverse trace $\bar{h}_{\mu\nu} = h_{\mu\nu} - \frac{1}{2}\eta_{\mu\nu} h$, where $ h = h^\mu_\mu$. Note that, $G_{\mu\nu} = -\frac{1}{2} \Box \bar{h}_{\mu\nu} $ and $ R = \bar{h}^{\mu\nu}_{\;\;\;\; ,\mu\nu} - \frac{1}{2} \Box h $. The lowest order terms are: \[ S = \frac{1}{\kappa}\int \mathrm{d}^4 x \left\{ -\frac{1}{4} \bar{h}^{\mu\nu} \Box \bar{h}_{\mu\nu} + \frac{1}{8} h ( 2 a_2 \Box^2 + \Box ) h - \frac{a_2}{2} \bar{h}^{\mu\nu} \partial_\mu\partial_\nu \Box h \right\} \]

As we can see, there is an additional field besides the conventional graviton field, with mass $m^2_h = -\frac{1}{2 a_2}$. This field interacts with matter via Yukawa coupling: $\sqrt{-g} \simeq 1 + \delta \sqrt{-g} = 1 + \frac{1}{2} h$. Therefore, it propagates $e^{- r /m_h}$. Taking $R \ll m^2_h \ll M_\text{Pl}$, $R^2$ gravity produces a massive weakly interacting particle.

[1]: Clifford M. Will,

*The Confrontation between General Relativity and Experiment*, Living Rev. Relativity, 9, (2006), 3, http://www.livingreviews.org/lrr-2006-3[2]: Jose A. R. Cembranos,

*Dark Matter from R2 Gravity*, Phys. Rev. Lett. 102, 141301 (2009)[3]: S. Stelle, Gen. Relativ. Gravit. 9, 353 (1978)

.