Variation of products of Riemann tensor $\delta (\sqrt{-g} RR \epsilon \epsilon)$

Question

This time i'd like to vary the following form of action.

I found it from appendix of "One-loop divergencies in the theory of gravitation" by G. 'T Hooft and M. Veltman.

Starting from variation by Lebiniz rule i have \begin{align} &\delta (\sqrt{-g} R_{\mu\nu\rho\sigma} R_{\alpha\beta\gamma\phi} \epsilon^{\mu\nu\alpha\beta} \epsilon^{\rho\sigma\gamma\phi}) \\ & = - \frac{1}{2} \sqrt{-g} g_{ab} \delta g^{ab} R_{\mu\nu\rho\sigma} R_{\alpha\beta\gamma\phi} \epsilon^{\mu\nu\alpha\beta} \epsilon^{\rho\sigma\gamma\phi} + 2\sqrt{-g} R_{\alpha\beta\gamma\phi} \delta R_{\mu\nu\rho\sigma} \epsilon^{\mu\nu\alpha\beta} \epsilon^{\rho\sigma\gamma\phi} + 2\sqrt{-g} R_{\alpha\beta\gamma\phi} R_{\mu\nu\rho\sigma} \epsilon^{\mu\nu\alpha\beta} \delta\epsilon^{\rho\sigma\gamma\phi} \end{align} Here $\epsilon^{abcd} = \frac{1}{\sqrt{-g}} \tilde{\epsilon}^{abcd}$ where $\tilde{\epsilon}$ is just number, $i.e$, $\epsilon^{abcd}$ is tensorial density. what i have in trouble is the last term, i simply noticed that \begin{align} \delta (\epsilon^{abcd}) = \delta \left( \frac{1}{\sqrt{-g}} \tilde{\epsilon}^{abcd}\right) = - \frac{1}{2} g^{\lambda\theta} \delta g_{\lambda\theta} \epsilon^{\alpha\beta\mu\nu} \end{align}

Computing $\delta R_{abcd}$ and from symmetric properties i obtain second term, but having some problem with third term. $i.e$, I obtain the first two terms in paper but having trouble obtaining the third term

The answer in paper says \begin{align} &\delta (\sqrt{g} RR \epsilon\epsilon) \\ &= -\frac{1}{2} g^{\alpha\beta} \delta g_{\alpha\beta} \sqrt{g} (RR \epsilon\epsilon) - 4 \sqrt{g} g_{\lambda\mu} (\nabla_\beta \delta \Gamma_\nu{}^\mu{}_\alpha) R_{\rho\theta \pi\delta} \epsilon^{\gamma \nu\pi \delta} \epsilon^{\rho\theta \alpha \beta} + 2 \sqrt{g} \delta g_{\gamma \mu} g^{\mu\pi} R_{\pi \theta \alpha \beta} R_{\rho \tau \nu\delta} \epsilon^{\gamma \theta \nu \delta} \epsilon^{\rho \tau \alpha \beta} \end{align}

I wonder how the third term came up. Am i missing somewhere?

---

Answer 1

\begin{align} &\delta R_{\mu\nu\rho\sigma} = g_{\mu\xi} ( \nabla_\rho \delta \Gamma_{\nu}{}^\xi{}_\sigma - \nabla_\sigma \Gamma_{\nu}{}^\xi{}_\rho) + \delta g_{\mu\xi} R^{\xi}{}_{\nu\rho\sigma} = g_{\mu\xi} ( \nabla_\rho \delta \Gamma_{\nu}{}^\xi{}_\sigma - \nabla_\sigma \Gamma_{\nu}{}^\xi{}_\rho) + \delta g_{\mu\xi} g^{\xi \eta} R_{\eta\nu\rho\sigma} \\ &\delta \epsilon^{\alpha\beta\mu\nu} = -\frac{1}{2} g^{\lambda\theta} \delta g_{\lambda\theta} \epsilon^{\alpha\beta\mu\nu} \end{align} using this I have same result in paper!

Thanks @AccidentalFourierTransform!