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>$$
Observes that \begin{eqnarray*} d^\nabla A(X, Y, Z)&:=& (\nabla_X A)(Y,Z) - (\nabla_Y A)(X, Z)\\ &=& X\left(A(Y, Z) \right)-A\left(\nabla_X Y, Z\right)-A\left(Y,\nabla_X Z \right)\\ &-& Y\left(A(X, Z) \right)+A\left(\nabla_Y X, Z\right)+A\left(X, \nabla_Y Z\right) \end{eqnarray*}
The Hessian of $\varphi$ is a (0,2)-tensor and is symmetric, that is, $\nabla^{2}\varphi(X,Y)=\nabla^{2}\varphi (Y,X)$ (see the Kuhnel book, p. 246).
Then,
\begin{eqnarray*} d^\nabla\nabla^{2}\varphi(X,Y,Z)&=& X\left( \nabla^{2}\varphi(Y,Z)\right)-\nabla^{2}\varphi\left(\nabla_X Y,Z \right)-\nabla^{2}\varphi\left(Y,\nabla_X Z \right)\\ &-& Y\left( \nabla^{2}\varphi(X,Z)\right)+\nabla^{2}\varphi\left(\nabla_Y X,Z \right)+\nabla^{2}\varphi\left(X, \nabla_Y Z \right) \\ \end{eqnarray*}
Now, remember that $\nabla^{2}\varphi(X,Y)=\left<\nabla_X\operatorname{grad}\varphi, Y\right>$, besides its symmetry:
\begin{eqnarray*} d^\nabla\nabla^{2}\varphi(X,Y,Z)&=& X\left(\left< \nabla_Y\operatorname{grad}\varphi, Z\right> \right)-\left<\nabla_Z\operatorname{grad}\varphi, \nabla_X Y\right>-\left<\nabla_Y\operatorname{grad}\varphi, \nabla_X Z\right>\\ &-&Y\left(\left< \nabla_X\operatorname{grad}\varphi, Z\right> \right)+ \left<\nabla_Z\operatorname{grad}\varphi, \nabla_Y X\right>+ \left<\nabla_X\operatorname{grad}\varphi, \nabla_Y Z\right>\\ \end{eqnarray*}
Using the derivative of the metric tensor: $X\left(\left< U,V\right>\right)=\left<\nabla_X U,V\right>+\left<U, \nabla_X V\right>$, you get
\begin{eqnarray*} d^\nabla\nabla^{2}\varphi(X,Y,Z)&=& \left<\nabla_X\nabla_Y\operatorname{grad}\varphi, Z\right>+\left<\nabla_Y\operatorname{grad}\varphi,\nabla_XZ\right>- \left<\nabla_Z\operatorname{grad}\varphi, \nabla_XY\right>-\left<\nabla_Y\operatorname{grad}\varphi, \nabla_X Z\right>\\ &-& \left<\nabla_Y\nabla_X\operatorname{grad}\varphi, Z\right>-\left<\nabla_X\operatorname{grad}\varphi, \nabla_YZ\right>+\left<\nabla_Z\operatorname{grad}\varphi, \nabla_YX\right>+\left<\nabla_X\operatorname{grad}\varphi, \nabla_YZ\right> \end{eqnarray*}
Check that some terms cancel. Use Hessian's symmetry again to get $\left<\nabla_Z\operatorname{grad}\varphi, \left[X,Y \right]\right> =\left<\nabla_{\left[X,Y \right]}\operatorname{grad}\varphi,Z \right>$. Using the tensor definition it has to be $\left<R(X,Y)W,Z\right>=\left< \nabla_X\nabla_YW-\nabla_Y\nabla_XW-\nabla_{\left[X,Y\right]}W, Z\right>$. Now, just take $W=\operatorname{grad}\varphi$, then
$$ d^\nabla \nabla^{2} \varphi (X, Y, Z) =\left<R(X, Y)\text{grad}\varphi, Z \right> $$
as desired. The other identities follow through similar calculations.