The action in $R^2$ gravity is given by \begin{equation} S[g^*_\mathrm{\mu\nu}] = \frac{1}{16\pi}\int d^4x \sqrt{-g^*}(R^*+aR{^*}^2) \end{equation} Where $g^* = \det(g^*_\mathrm{\mu\nu})$ and $a$ is a constant. For my problem I neglect the matter term. I've allready transformed this action in the Einstein frame by using a conformal transformation of the metric $g^*_\mathrm{\mu\nu}=A^2g_\mathrm{\mu\nu}$ \begin{equation} S[g_\mathrm{\mu\nu},\phi] = \frac{1}{16\pi}\int d^4x \sqrt{-g}(R-2\partial_\mu\phi\partial^\mu\phi-V(\phi)) \end{equation} Here is $\phi$ a scalar field and $V(\phi)$ a potential. Now I want to derive the filed equation from this action in the Einstein frame. For this I define $f(R) = R-2\partial_\mu\phi\partial^\mu\phi-V(\phi)$. So the variation respect to the metric is \begin{align} \delta S &= \int d^4x \left(f(R)\delta\sqrt{-g}+\sqrt{-g}\delta f(R)\right)\\ &= \int d^4x \left(f(R)\delta\sqrt{-g}+\sqrt{-g}\frac{df}{dR}\delta R\right) = 0 \end{align} After some calculations there are those two terms in the action \begin{equation} \delta S \propto \int d^4x\sqrt{-g}\frac{df}{dR}\left(g_\mathrm{\mu\nu}g^\mathrm{\gamma\alpha}\nabla_\alpha\nabla\gamma\delta g^\mathrm{\mu\nu}-\nabla_\mu\nabla_\nu\delta g^\mathrm{\mu\nu}\right) \end{equation} The problem is that the derivatives relate to the variation in the metric. To work around this problem I want to integrate by parts twice. So that the derivatives are on $\frac{df}{dR}$. But what is this term $\frac{df}{dR}$? I think it is equal to one. But the derivatives on this constant vanishes, right? Nevertheless, In the final field equations there are these terms \begin{equation} R_\mathrm{\mu\nu}-\frac{1}{2}g_\mathrm{\mu\nu}R = 2\partial_\mu\phi\partial_\nu\phi-g_\mathrm{\mu\nu}\partial^\sigma\phi\partial_\sigma\phi-\frac{1}{2}V(\phi)g_\mathrm{\mu\nu} \end{equation} $\\$
Edit: The variation of the action respect to the metric is \begin{align} \delta S &= \int d^4x \sqrt{-g}\left[(R_\mathrm{\mu\nu}-\frac{1}{2}g_\mathrm{\mu\nu}R)\delta g^\mathrm{\mu\nu}+\frac{1}{2}g_\mathrm{\mu\nu}V(\phi)\delta g^\mathrm{\mu\nu}+g_\mathrm{\mu\nu}\partial_\gamma\phi\partial^\gamma\phi\delta g^\mathrm{\mu\nu}\right]\\ &+\int d^4x \sqrt{-g}\left(g_\mathrm{\mu\nu}\Box\delta g^\mathrm{\mu\nu}-\nabla_\mu\nabla_\nu\delta g^\mathrm{\mu\nu}\right) \end{align} where I set $\frac{df}{dR}=1$. I don't know how to calculate the second integral.