I'm dealing with some extension of GR, with action:
$S=\int d^4x\Big[\sqrt{-g} f(R,R_{\mu\nu}R^{\mu\nu})$
Varying this action gives:
$\delta S=\int d^4x\Big[\delta\sqrt{-g} f(R,R_{\mu\nu}R^{\mu\nu})+\sqrt{-g} (\delta R f_R+ \delta(R_{\mu\nu}R^{\mu\nu})f_Y)$ where $Y=R_{\mu\nu}R^{\mu\nu}$ and $f_Y=\frac{\partial f}{\partial{Y}}$ etc
$f(R)$ part is relatively simple; while variation of the second term is:
$\delta(R_{\mu\nu}R^{\mu\nu})f_Y=2(R^{\mu\nu}\delta R_{\mu\nu}+g^{\rho\sigma}R_{\mu\rho}R_{\nu\sigma}\delta g^{\mu\nu})f_Y$
Using relation:
$ \delta R_{\mu\nu}=\frac{1}{2}g^{\alpha\beta}[\nabla_{\alpha}\nabla_{\nu}\delta g_{\mu\beta}+\nabla_{\alpha}\nabla_{\mu}\delta g_{\nu\beta}-\nabla_{\mu}\nabla_{\nu}\delta g_{\alpha\beta}-\nabla_{\alpha} \nabla_{\beta}\delta g_{\mu\nu}] $
Now, integration by parts 2 times, gives:
$\int...g^{\alpha\beta}\nabla_{\alpha}\nabla_{\nu}(\delta g_{\mu\beta})R^{\mu\nu}f_Y= \int ...g^{\alpha\beta}\nabla_{\alpha}\nabla_{\nu}(R^{\mu\nu}f_Y)\delta g_{\mu\beta} $
and
$\int ...g^{\alpha\beta}\nabla_{\alpha}\nabla_{\mu}(\delta g_{\nu\beta})R^{\mu\nu}f_Y= \int ...g^{\alpha\beta}\nabla_{\alpha}\nabla_{\mu}(R^{\mu\nu}f_Y)\delta g_{\nu\beta}$
$\int...-g^{\alpha\beta}\nabla_{\mu}\nabla_{\nu}(\delta g_{\alpha\beta})R^{\mu\nu}f_Y=\int ...-g^{\alpha\beta}\nabla_{\mu}\nabla_{\nu}(R^{\mu\nu}f_Y) \delta g_{\alpha\beta}$
$\int ... -g^{\alpha\beta}\nabla_{\alpha} \nabla_{\beta}(\delta g_{\mu\nu})R^{\mu\nu}f_Y=\int ...-g^{\alpha\beta}\nabla_{\alpha} \nabla_{\beta}(R^{\mu\nu}f_Y)\delta g_{\mu\nu}$
Where boundary term is assumed to vanish. Apllying $\delta g_{\mu\nu}=-g_{\mu\rho} g_{\nu\sigma}\delta g^{\rho\sigma}$
leads to
$g^{\alpha\beta}\nabla_{\alpha}\nabla_{\nu}(R^{\mu\nu}f_Y)\delta g_{\mu\beta}=-g^{\alpha\beta}\delta g^{\rho\sigma}g_{\mu\rho} g_{\beta\sigma}\nabla_{\alpha}\nabla_{\nu}(R^{\mu\nu}f_Y)=-\delta g^{\rho\sigma}\nabla_{\sigma}\nabla_{\nu}(R_{\rho}^{\;\;\nu}f_Y)$
$ g^{\alpha\beta}\nabla_{\alpha}\nabla_{\mu}(R^{\mu\nu}f_Y)\delta g_{\nu\beta}=-g^{\alpha\beta}\delta g^{\rho\sigma}g_{\nu\rho}g_{\beta\sigma}\nabla_{\alpha}\nabla_{\mu}(R^{\mu\nu}f_Y)=-\delta g^{\rho\sigma}\nabla_{\sigma}\nabla_{\mu}(R^{\mu}_{\;\;\rho}f_Y ) $
$ -g^{\alpha\beta}\nabla_{\mu}\nabla_{\nu}(R^{\mu\nu}f_Y) \delta g_{\alpha\beta}=g^{\alpha\beta}\delta g^{\rho\sigma}g_{\alpha\rho}g_{\beta\sigma} \nabla_{\mu}\nabla_{\nu}(R^{\mu\nu}f_Y) =\delta g^{\rho\sigma} g_{\rho\sigma}\nabla_{\mu}\nabla_{\nu}(R^{\mu\nu}f_Y)$
$-g^{\alpha\beta}\nabla_{\alpha} \nabla_{\beta}(R^{\mu\nu}f_Y)\delta g_{\mu\nu}=\delta g^{\rho\sigma}\Box(R_{\sigma\rho}f_Y)$
while
$g^{\rho\sigma}R_{\mu\rho}R_{\nu\sigma}\delta g^{\mu\nu}f_Y=\delta g^{\rho\sigma}R^{\nu}_{\rho}R_{\nu\sigma}f_Y$
Then, variation will be:
$ \delta(R_{\mu\nu}R^{\mu\nu})f_Y=2(R^{\mu\nu}\delta R_{\mu\nu}+g^{\rho\sigma}R_{\mu\rho}R_{\nu\sigma}\delta g^{\mu\nu})f_Y=\delta g^{\rho\sigma}[\Box(R_{\sigma\rho}f_Y)+g_{\rho\sigma}\nabla_{\mu}\nabla_{\nu}(R^{\mu\nu}f_Y)+2R^{\nu}_{\rho}R_{\nu\sigma}f_Y-2\frac{1}{2}\nabla_{\sigma}\nabla_{\mu}(R^{\mu}_{\;\;\rho}f_Y )-2\frac{1}{2}\nabla_{\sigma}\nabla_{\nu}(R_{\rho}^{\;\;\nu}f_Y)]$.
Terms in bracket will appear in final field quations when other $...$ terms are present.
Where i have some problems with antysymetrization $C_{(ab)}=\frac{1}{2}(C_{ab}+C_{ba})$ - i know that since i am dealing with symmetric tensors i can switch indices:
$-\big(\frac{1}{2}\nabla_{\sigma}\nabla_{\mu}(R^{\mu}_{\;\;\rho}f_Y )+\frac{1}{2}\nabla_{\sigma}\nabla_{\nu}(R_{\rho}^{\;\;\nu}f_Y)\big)=-\frac{1}{2}(\nabla_{\sigma}\nabla_{\mu}(R^{\mu}_{\;\;\rho}f_Y )+\nabla_{\rho}\nabla_{\nu}(R_{\sigma}^{\;\;\nu}f_Y)\big))$
Hovewer, shouldn't antisymmetrization of the indices come from actual variational principle? I mean in paper https://arxiv.org/abs/1503.08751 they've could just use $-2\nabla_{\mu}\nabla_{\alpha}(R^{\alpha}_{\nu}f_Y)$ - am i missing something in my derivation of this part of the field equations?