Derivation of $f(R)$ field equations, problem with integration by parts I am following the derivation of the field equations on the the Wikipedia page for $f(R)$ gravity.
But I do not understand the following step:
$$
\delta S = \int \frac{1}{2\kappa} \sqrt{-g} \left(\frac{\partial f}{\partial R} (R_{\mu\nu} \delta g^{\mu\nu}+g_{\mu\nu}\Box \delta g^{\mu\nu}-\nabla_\mu \nabla_\nu \delta g^{\mu\nu}) -\frac{1}{2} g_{\mu\nu} \delta g^{\mu\nu} f(R) \right) 
$$
the wiki article says, the next step is to integrate the second and third terms by parts to yield:
$$
\delta S = \int \frac{1}{2\kappa} \sqrt{-g}\delta g^{\mu\nu} \left(\frac{\partial f}{\partial R} R_{\mu\nu}-\frac{1}{2}g_{\mu\nu} f(R)+[g_{\mu\nu}\Box -\nabla_\mu \nabla_\nu] \frac{\partial f}{\partial R} \right)\, \mathrm{d}^4x 
$$
In other words, integrating by parts should yield:
$$
\int \sqrt{-g} \left(\frac{\partial f}{\partial R} (g_{\mu\nu}\Box \delta g^{\mu\nu}-\nabla_\mu \nabla_\nu \delta g^{\mu\nu}) \right)\, d^4x  $$
$$= \int \sqrt{-g}\delta g^{\mu\nu} \left([g_{\mu\nu}\Box -\nabla_\mu \nabla_\nu] \frac{\partial f}{\partial R} \right) \mathrm{d}^4x
$$
From there getting the usual f(R) field equations is trivial. What I'm confused by is how to integrate by parts to get that. 
I have tried many different ways the one I think is most correct is:
assuming $g_{\mu \nu} \Box$ and $\nabla_\mu \nabla_\nu$ are differential operators then $u' = g_{\mu \nu} \Box \delta g^{\mu\nu}$ and $v = f'$, similarly with the $\nabla_\mu \nabla_\nu$ so using the formula for integration by parts:
$$
\int u'v = uv -\int uv'
$$
I get:
$$
\int \sqrt{-g} \left(f' (g_{\mu\nu}\Box \delta g^{\mu\nu}-\nabla_\mu \nabla_\nu \delta g^{\mu\nu}) \right)\, d^4x $$
$$= -\int \sqrt{-g}\delta g^{\mu\nu} \left([g_{\mu\nu}\Box -\nabla_\mu \nabla_\nu] f' \right) \mathrm{d}^4x
$$
because the $uv$ term will disappear.
So can any one explain to me why I have the minus sign and Wikipedia doesn't? 
Is it ok to use $g_{\mu \nu} \Box$ as a differential operator? I have tried other ways such as writing $\Box$ explicitly and using integration by parts twice but I also couldn't get the correct answer as i end up with terms such as $\nabla_\nu \nabla_\mu$ which cant be correct.
There is a similar post on physics forums on this step but it does not answer my question and is now closed. 
 A: You need to use recursive integration by parts to deal with the second derivatives.
See, for example, equation (18) here.
The last term in (18) has a multiplying factor of $(-1)^n$, where $n$ is the order of the derivative. In your case, $n=2$ and the minus sign vanishes.
A: The box operator $\Box = g^{αβ}\nabla_{α}\nabla_{β} = \nabla^{α}\nabla_{α}$ and the two covariant derivatives $\nabla_{μ}\nabla_{ν}$ act on the variation of the inverse metric tensor. You want $δg^{μν}$ to be a multiplying factor so you must integrate by parts twice (because you have two derivatives) to make the derivatives act on $\cfrac{df}{dR}$. The box term will be: 
$f'(R)\nabla^{α}\nabla_{α}(δg^{μν})$
From Leibniz rule we know that: 
$\nabla^{α}(f'(R)\nabla_{α}δg^{μν}) = \nabla^{α}f'(R)\nabla_{α}δg^{μν} + f'(R)\nabla^{α}\nabla_{α}δg^{μν}$
The left hand side term is a total derivative so it is zero which means: 
$f'(R)\nabla^{α}\nabla_{α}δg^{μν} = - \nabla^{α}f'(R)\nabla_{α}δg^{μν}$
Applying Leibniz rule again: 
$\nabla_{α}(-\nabla^{α}f'(R)δg^{μν}) = - \nabla_{α}\nabla^{α}f'(R)δg^{μν} - \nabla_{α}δg^{μν}\nabla^{α}f'(R) $
The left hand side term is a total derivative. 
$- \nabla_{α}δg^{μν}\nabla^{α}f'(R) = \nabla_{α}\nabla^{α}f'(R)δg^{μν} \Rightarrow $ 
$f'(R)\nabla^{α}\nabla_{α}δg^{μν} = \nabla_{α}\nabla^{α}f'(R)δg^{μν} $.
You have to do the same with the $f'(R)\nabla_{μ}\nabla_{ν}δg^{μν}$ term. 
The box term also contains $g_{μν}$ but metric is compatible so: 
$\nabla^{α}\nabla_{α}(g_{μν}f'(R)δg^{μν}) = \nabla^{α}\nabla_{α}(f'(R)δg^{μν})g_{μν}$
EDIT: Another derivation can be found here: https://arxiv.org/abs/1002.0617
