Quadratic from Riemann tensor with zero divergence I am looking for the following tensor field $T^i_j$ on closed riemann manifold $M$ which is dependent only from metric tensor $g_{ij}$ and its Riemann cirvature tensor $R^i_{jkl}$ and moreover quadratic by $R^i_{jkl}$ which have zero divergence:
$$
\nabla_i T^i_j = 0,
$$
For example, the Einstein tensor have zero divergence but it is linear by Ricci tensor and hence linear by Riemann curvature tensor. On the other hand the following tensor:
$$
T^i_j = R_j^i R,
$$
is quadratic by $R^i_{jkl}$ but in general have non zero divergence in the sense given above.
 A: There might be two ways to construct the expected tensor.
First, perhaps you may want to have a look at $f(S)$ gravity, where $S$ means the scalar curvature. This method considers the following functional of the metric $g$
$$
I[g]=\int_{M}f(S)\,{\rm d}V,
$$
where $M$ stands for the closed Riemannian manifold, ${\rm d}V$ is its invariant volume form, and $f$ denotes some smooth function. In the special case $f(S)=S$, this functional reduces to the Einstein-Hilbert action. Now, apply calculus of variation to $I[g]$ with respect to $g$, and Noether's theorem guarantees its functional derivative
$$
\frac{\delta I[g]}{\delta g^{\mu\nu}}=f'(S)\,{\rm Ric}_{\mu\nu}-\frac{1}{2}\,g_{\mu\nu}\,f(S)+\left(g_{\mu\nu}\,g^{\sigma\tau}\nabla_{\sigma}\nabla_{\tau}-\nabla_{\mu}\nabla_{\nu}\right)f'(S)
$$
to be divergence free, where ${\rm Ric}$ is the Ricci tensor, and $\nabla$ means the Levi-Civita connection. Finally, by choosing, e.g., $f(S)=S^2$, you may obtain a divergence-free tensor with its main part (first two terms in the last formulation) quadratic in the Riemann curvature tensor. Unfortunately, as you will see, this result includes additional derivative terms, hence is not a fully quadratic expression.
Second, perhaps you may want to explore a most general form of the expected quadratic tensor
$$
T=\left<C,R\otimes R\right>,
$$
where $C$ is a $g$-isotropic tensor, $R$ means the Riemann curvature tensor, and $\left<\cdot,\cdot\right>$ denotes the tensor contraction operator. Its entrywise form reads
$$
T_{\xi\eta}=C_{\tau\rho\xi\eta}^{\mu\nu\sigma\alpha\beta\gamma}\,R_{\mu\nu\sigma}^{\tau}\,R_{\alpha\beta\gamma}^{\rho}.
$$
Here $g$-isotropic means that $C$ remains identical under all isometries that preserves $g$. Therefore, if $C$ is rank-two, it must be
$$
C_{\mu\nu}=K\,g_{\mu\nu}
$$
with $K$ being a constant; if $C$ is rank-four, it must be
$$
C_{\mu\nu\alpha\beta}=K_1\,g_{\mu\nu}\,g_{\alpha\beta}+K_2\,g_{\mu\alpha}\,g_{\nu\beta}+K_3\,g_{\mu\beta}\,g_{\nu\alpha}
$$
with $K_1$, $K_2$, $K_3$ being constants. Now that $C$ is rank-ten from above, and you will have 945 $K_j$'s (of course, provided that $R$ is partly anti-symmetric, many of these constants vanish). Nevertheless, since $C$ is made up of $g$, it appears to be a "constant" tensor as you apply $\nabla$ to it, i.e., $\nabla_{\mu}C_{\cdots}=0$. Hence your goal reduces to explore the symmetric property of this rank-ten $C$ upon the requirement
$$
\nabla^{\xi}T_{\xi\eta}=0.
$$
Unfortunately, I am afraid this is so huge a task that even a special trial takes much more than expected.
