Variation of the Ricci tensor "squared" and antisymmetrization of the derivatives 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.
This part should be equal to:
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? 
 A: The first three terms in your variation of $\delta(R_{\mu\nu}R^{\mu\nu})f_Y$ can be found in what you say is your result, so let's focus on the last two terms.
First note that, by integrating by parts, you are changing the order of differentiation. What this means is that, if in the original expression you first differentiate with $\nabla_\alpha$ and then with $\nabla_\beta$, integrating by parts results in an expression where you fist differentiate with $\nabla_\beta$ and then with $\nabla_\alpha$. These two orderings are not equivalent, because spacetime is, in general, curved.
Having said that, recall that indices which are summed over are called dummy indices, and can be relabelled to whatever is convenient. The last two terms in your variation of $\delta(R_{\mu\nu}R^{\mu\nu})f_Y$ are therefore equal, as
$$
\nabla_\mu\nabla_\sigma(R^\mu{}_\rho f_Y) = \nabla_\nu\nabla_\sigma(R^\nu{}_\rho f_Y)
=[R^\nu{}_\rho f_Y]_{;\sigma\nu}.
$$
I've taken the liberty of changing the order of differentiation, as I've argued above, and used the semicolon notation for covariant derivatives in order to compare with your reference expression. As you noticed, the only thing that is missing is the symmetrisation over the indices $\rho$ and $\sigma$. As you expected, this symmetrisation is due to the variation: the variations $\delta g_{\mu\nu}$ and $\delta g^{\mu\nu}$ are symmetric. Therefore, the two are related, not just by your expression, but by
$$
\delta g_{\mu\nu} = -g_{\mu(\rho}g_{\sigma)\nu}\delta g^{\rho\sigma}.
$$
Because this expression is contracted with something else such a symmetrisation is normally taken care of automatically. When varying the action, however, you break the contraction. You would then lose this symmetry unless you make it explicit.
