You’re staring at the “variation-of-Ricci” mess that always shows up when deriving the field eqs from an R2R^2R2–type action. The cure is three identities; once you use them in the right order, the whole block collapses to
gμν□R+R Rμν−∇μ∇νR.g_{\mu\nu}\Box R + R\,R_{\mu\nu} - \nabla_\mu\nabla_\nu R .gμν□R+RRμν−∇μ∇νR.
Here’s the minimal path.
1) Palatini + variation of the connection
Palatini identity:
δRαβ=∇ρδΓραβ−∇βδΓραρ.\delta R_{\alpha\beta}=\nabla_\rho \delta\Gamma^{\rho}{}_{\alpha\beta}-\nabla_\beta \delta\Gamma^{\rho}{}_{\alpha\rho}.δRαβ=∇ρδΓραβ−∇βδΓραρ.
Metric variation of the Levi–Civita connection:
δΓραβ=12gρλ (∇αδgλβ+∇βδgλα−∇λδgαβ).\delta\Gamma^{\rho}{}_{\alpha\beta}=\tfrac12 g^{\rho\lambda}\!\left(\nabla_\alpha \delta g_{\lambda\beta}+\nabla_\beta \delta g_{\lambda\alpha}-\nabla_\lambda \delta g_{\alpha\beta}\right).δΓραβ=21gρλ(∇αδgλβ+∇βδgλα−∇λδgαβ).
Insert this in Palatini and simplify (use ∇γgρλ=0\nabla_\gamma g^{\rho\lambda}=0∇γgρλ=0). After a standard symmetrization you get the textbook result:
δRαβ=12(∇ρ∇αδgρβ+∇ρ∇βδgρα−□ δgαβ−∇α∇βh) (h≡gρσδgρσ).\boxed{\,\delta R_{\alpha\beta}=\tfrac12\Big(\nabla_\rho\nabla_\alpha \delta g^\rho{}_\beta+\nabla_\rho\nabla_\beta \delta g^\rho{}_\alpha-\Box\,\delta g_{\alpha\beta}-\nabla_\alpha\nabla_\beta h\Big)\,}\qquad (h\equiv g^{\rho\sigma}\delta g_{\rho\sigma}).δRαβ=21(∇ρ∇αδgρβ+∇ρ∇βδgρα−□δgαβ−∇α∇βh)(h≡gρσδgρσ).
Equivalently, the “functional derivative” object you wrote (those ∇δΓ\nabla \delta\Gamma∇δΓ terms divided by δgμν\delta g^{\mu\nu}δgμν) is nothing but the kernel that maps δgμν\delta g^{\mu\nu}δgμν to δRαβ\delta R_{\alpha\beta}δRαβ above.
2) Vary
The piece you posted is precisely the part linear in δRαβ\delta R_{\alpha\beta}δRαβ:
2 Rαβ gαγgβδ δRγδ.2\,R_{\alpha\beta}\,g^{\alpha\gamma}g^{\beta\delta}\,\delta R_{\gamma\delta}.2RαβgαγgβδδRγδ.
Now insert the boxed formula for δRγδ\delta R_{\gamma\delta}δRγδ, integrate by parts twice (or, since you’re effectively taking a functional derivative, move derivatives off δg\delta gδg onto RαβR_{\alpha\beta}Rαβ; throw away total divergences). You’ll produce derivatives of RαβR_{\alpha\beta}Rαβ.
Use the contracted Bianchi identity:
∇αRαβ=12 ∇βR,∇α∇βRαβ=12 □R.\nabla^\alpha R_{\alpha\beta}=\tfrac12\,\nabla_\beta R,\qquad\nabla^\alpha\nabla^\beta R_{\alpha\beta}=\tfrac12\,\Box R.∇αRαβ=21∇βR,∇α∇βRαβ=21□R.
After the dust settles, the “∇∇ δg\nabla\nabla\,\delta g∇∇δg” block yields
gμν□R−∇μ∇νR .\boxed{\ g_{\mu\nu}\Box R-\nabla_\mu\nabla_\nu R\ }. gμν□R−∇μ∇νR .
3) The algebraic (no-derivative) part
The two “Ricci×Ricci” contractions at the front,
RαμRγνgαγ+RμβRνδgβδ,R_{\alpha\mu}R_{\gamma\nu}g^{\alpha\gamma}+R_{\mu\beta}R_{\nu\delta}g^{\beta\delta},RαμRγνgαγ+RμβRνδgβδ,
are the same term written twice:
RαμRαν+RμβRβν=2 RμαRαν.R_{\alpha\mu}R^{\alpha}{}_{\nu}+R_{\mu\beta}R^{\beta}{}_{\nu}=2\,R_{\mu}{}^{\alpha}R_{\alpha\nu}.RαμRαν+RμβRβν=2RμαRαν.
In the full variation of ∫ −g R2\int \! \sqrt{-g}\,R^2∫−gR2 this combines with the piece from δ(−g)\delta(\sqrt{-g})δ(−g) and δ(gαγgβδ)\delta(g^{\alpha\gamma}g^{\beta\delta})δ(gαγgβδ) to trade RμαRανR_\mu{}^\alpha R_{\alpha\nu}RμαRαν for R RμνR\,R_{\mu\nu}RRμν and −12gμνR2-\tfrac12 g_{\mu\nu}R^2−21gμνR2. Since your target line shows only the part that comes with δR\delta RδR, you keep just
R Rμν .\boxed{\ R\,R_{\mu\nu}\ }. RRμν .
Result (the bit you wanted)
Collecting the derivative and algebraic pieces that come from the δRαβ\delta R_{\alpha\beta}δRαβ-block gives exactly
gμν □R + R Rμν − ∇μ∇νR .\boxed{\ g_{\mu\nu}\,\Box R \;+\; R\,R_{\mu\nu} \;-\; \nabla_\mu\nabla_\nu R\ }. gμν□R+RRμν−∇μ∇νR .
sanity check: for the R2R^2R2 Lagrangian one indeed gets
Hμν(R2)=2RRμν−12gμνR2−2∇μ∇νR+2gμν□RH_{\mu\nu}^{(R^2)}=2R R_{\mu\nu}-\tfrac12 g_{\mu\nu}R^2-2\nabla_\mu\nabla_\nu R+2 g_{\mu\nu}\Box RHμν(R2)=2RRμν−21gμνR2−2∇μ∇νR+2gμν□R;
your line is the “half” that comes from the δRαβ\delta R_{\alpha\beta}δRαβ part, as advertised.
If you want, I can write this out explicitly index-by-index with every integration-by-parts step, but the three boxes above are the whole trick.
