I am working on the following theorem
Theorem (J. A. Schouten) For $n\geq 4$ the metric $g$ is conformally flat if and only if $W=0$. For $n=3$, the metric $g$ is conformally flat if and only if the relation $$ (\nabla_X C) (Y, Z) = (\nabla_Y C) (X, Z)$$ holds for all $X, Y, Z$. Here $C$ and $W$ denote the Schouten tensor and the Weyl tensor with $R=C\bullet g+W$.
I am following Theorem 8.31, pg. 352 in Differential Geometry: Curves - Surfaces - Manifolds by Wolfgang Kühnel. For the demonstration he uses the following derivation:
$$d^\nabla A(X, Y, Z):= (\nabla_X A)(Y,Z) - (\nabla_Y A)(X, Z)$$
for an arbitrary symmetric (0,2)-tensor. So he concludes, for a scalar function $\varphi$:
\begin{eqnarray*} d^\nabla \nabla^{2} \varphi (X, Y, Z) &=& \left<R(X, Y)\text{grad}\varphi, Z \right>\\ d^\nabla\left( \nabla\varphi\cdot\nabla \varphi \right)(X, Y, Z) &=& (Y\varphi)\nabla^{2}\varphi(X, Z) - (X\varphi)\nabla^{2}\varphi(Y, Z)\\ d^\nabla\left(\frac{1}{2}\vert\vert \text{grad}\varphi\vert\vert\cdot g\right)(X,Y,Z) &=& \nabla^{2}\varphi (X,\text{grad}\varphi)\left<Y,Z\right>-\nabla^{2}\varphi(Y,\text{grad}\varphi)\left<X,Z\right> \end{eqnarray*}
where $R$ is the curvature tensor. Similarly, for a one-form $\alpha=d\varphi$ he concludes:
\begin{eqnarray*} d^\nabla\alpha(X, Y, Z) &=& - \alpha(R(X, Y)Z)\\ d^\nabla(\alpha\cdot\alpha)(X,Y,Z) &=& d\alpha(X,Y)\alpha(Z)+\alpha(Y)\nabla\alpha(X,Z)-\alpha(X)\nabla\alpha(Y,Z)\\ d^\nabla\left(\frac{1}{2}\vert\vert \alpha\vert\vert^{2}\cdot g\right)(X,Y,Z) &=& \left<\nabla_X\alpha, \alpha\right>\left<Y,Z\right> - \left<\nabla_Y\alpha, \alpha\right>\left<X,Z\right> \end{eqnarray*}
I am having difficulties in obtaining these six equalities above. I tried several times. For each of the identities above, the desired term appears but still contains some other terms that have not been canceled. I looked for other references that may present the proof of this theorem but I haven't found it.
For the first equality, using the definition of $d^\nabla$ in $d^\nabla\nabla^{2}\varphi(X,Y,Z)$, I got $$d^\nabla\nabla^{2}\varphi(X,Y,Z) = \left<R(X,Y)\text{grad}\varphi, Z\right>+\left<\nabla_{[X,Y]}\text{grad}\varphi, Z\right>+\left<\nabla_{Y}\text{grad}\varphi, \nabla_X Z\right>-\left<\nabla_{X}\text{grad}\varphi,\nabla_Y Z\right> $$
And for the second:
$$ d^\nabla\left( \nabla\varphi\cdot\nabla \varphi \right)(X, Y, Z) = (Y\varphi)\nabla^{2}\varphi(X, Z) - (X\varphi)\nabla^{2}\varphi(Y, Z)+Z(\varphi)\left<\text{grad}\varphi, [X,Y] \right> + Y(\varphi)\left<\text{grad}\varphi, \nabla_X Z\right>-X(\varphi)\left<\text{grad}\varphi, \nabla_Y Z \right>$$