The Schouten tensor is given by $$ P_{i j}=\frac{1}{n-2}\left(R_{i j}-\frac{1}{2(n-2)} R g_{i j}\right).$$
I want to compute the divergence of the Bach tensor in dimension $4$. From this post, we know that the Bach tensor is divergence-free in dimension $4.$ I want to show it explicitly.
The divergence of the back tensor
$$ \begin{aligned} \nabla^{j} B_{i j} &=\nabla^{j} ( \nabla^{k} \nabla_{k} P_{i j}-\nabla^{j} \nabla^{k} \nabla_{j} P_{i k})+\nabla^{j}\left(P^{k l} W_{i k j l}\right) \\ &= (\nabla^{j} \nabla^{k} \nabla_{k} P_{i j}-\nabla^{k} \nabla^{j} \nabla_{k} P_{i j})+\nabla^{j}\left(P^{k l} W_{i k j l}\right) \end{aligned} .$$ The second term requires the divergence of the Weyl tensor, and I have a nice formula for the divergence of the Weyl tensor. I'm struggling to compute the following portion:
$$ \begin{aligned} (\nabla^{j} \nabla^{k} \nabla_{k} P_{i j}-\nabla^{k} \nabla^{j} \nabla_{k} P_{i j}) \end{aligned} .$$
On page 126 of Jiaqi Chen's Doctoral thesis, he says the Ricci identity gives us
$$ \begin{aligned} (\nabla^{j} \nabla^{k} \nabla_{k} P_{i j}-\nabla^{k} \nabla^{j} \nabla_{k} P_{i j}) \end{aligned} $$ equals to $$(\nabla^{j} P^{k l}) R_{i k j l}.$$ I could not derive this equality by the Ricci Identity. Could you help me establish the equality that Jiaqi Chen claimed? Thanks so much. Note that here our connection $\nabla$ is torsion-free and metric compatible.
Update: Finally, I did the calculation and it was nasty. I'm still wondering how Dr. Chen did this in just two lines. I'm still looking for a better/short calculation.