general relativity - Deriving Gauss-Bonnet Gravity (Or just higher order corrections)

I have been working for some time now on deriving the equations of motion (EOM) for the Gauss-Bonnet Gravity, which is given by the action:

$$\int d^D x \sqrt{|g|} (R^2-4R_{ab}R^{ab}+R_{abcd}R^{abcd}).$$

I've tried for some time to derive the first-order variation with respect to $\delta g_{ab}$. I have always been stuck with the derivation for the $\delta (R_{ab}R^{ab})$ component. Suppose I have only a Riemann tensor that contains only constants,

$$\begin{align*} \delta(R_{ab}R^{ab}) & = (\delta R_{ab})R^{ab} + (\delta R^{ab})R_{ab} \\ & = (\delta R_{ab})R^{ab}+ \delta(R_{ecfd}g^{ac}g^{bd})g^{ef}R_{ab} \\ & = (\delta R_{ab})R^{ab}+ (\delta R_{ecfd})g^{ac}g^{bd}g^{ef}R_{ab} + (\delta g^{ac}) R_{cd}g^{bd}R_{ab}+(\delta g^{bd})R^{a}{}_{d}R_{ab}\\&\quad+ (\delta g^{ef}) R_{e}{}^{a}{}_{f}{}^{b}R_{ab}. \end{align*}$$

This however, gives us:

$$\begin{align*} & = - (\delta g^{hi}) R_{ichd}R^{cd} +( \delta g^{hi}) R_{hcid}R^{cd} + 2(\delta g^{ac}) R_{cd}R_{a}{}^{d} \\ & = 2(\delta g^{ac}) R_{cd}R_{a}{}^{d}. \end{align*}$$

The correct answer as we know it, is $2R^{ab}R_{acbd}\delta g^{cd}$.

Answer

The way to resolving this is as follows. First, don't assume that the Riemann tensor is a constant:

$$\begin{align*} \delta(R_{\mu\nu}R^{\mu\nu}) & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha} + \square R_{\mu\nu} -\nabla_{\alpha}\nabla_{\mu}R_{\nu}{}^\alpha - \nabla_{\alpha}\nabla_{\nu}R_{\mu}{}^\alpha + g_{\mu\nu} \nabla_{\alpha}\nabla_{\beta}R^{\alpha\beta}]\delta g^{\mu\nu} \\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}-((\nabla_{\mu}\nabla_{\alpha}+ \square R_{\mu\nu} -[\nabla_\mu,\nabla_\alpha])R_{\nu}{}^\alpha \\&\quad+(\nabla_{\nu}\nabla_{\alpha}-[\nabla_\nu,\nabla_\alpha])R_{\mu}{}^\alpha) + g_{\mu\nu} \nabla_{\alpha}\nabla_{\beta}R^{\alpha\beta}]\delta g^{\mu\nu}\\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}+ \square R_{\mu\nu} \\&\quad-\nabla_\mu\nabla_\nu R + [\nabla_\mu,\nabla_\alpha]R_{\nu}{}^\alpha + [\nabla_\nu,\nabla_\alpha]R_{\mu}{}^\alpha + \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu}\\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} + 2 R_{\mu}{}^{\alpha}R_{\nu\alpha}+ \square R_{\mu\nu} -\nabla_\mu\nabla_\nu R +2(-R_{\alpha \nu}R^{\alpha}{}_\mu \\& \quad+ R^{\alpha\beta}R_{\mu\alpha\nu\beta})+ \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu} \\ & = [-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta} -\nabla_\mu\nabla_\nu R + \square R_{\mu\nu} +2R^{\alpha\beta}R_{\mu\alpha\nu\beta}+ \frac{1}{2}g_{\mu\nu} \square R]\delta g^{\mu\nu}. \end{align*}$$

Dropping the constant terms, we arrive at $-\frac{1}{2} g_{\mu\nu}R_{\alpha\beta}R^{\alpha\beta}+2R^{\alpha\beta}R_{\mu\alpha\nu\beta}$, including the term from varying $\sqrt{|g|}$. The problem with the previous solution is in the transition from the second to the third line, because one still has to vary the metric, and one should vary everything first without assuming constants, since the commutator of the covariant derivatives do not commute. I think the problem arises because even if the Riemann tensor is a constant, the metric that comes up with it does not necessarily (e.g. $S^2$) that it is constant, so we can't just say $\delta R_{abcd} = 0$, and we have to check the parts of it.