Let $(M^n,g)$ be a Riemannian manifold, and $T$ a symmetric $(1,1)$-tensor field, i.e., $\langle T(X),Y\rangle = \langle X,T(Y)\rangle $. For convenience, denote $$\Delta\_Tu=\sum\_i\langle \nabla\_{e\_i}\nabla u, Te\_i\rangle $$
and
$$\mathrm{Ric}\_T(X,Y)=\sum\_i\langle R(X,e\_i)(Te\_i), Y\rangle , $$
where $u$ is a smooth function on $M$ and $\{e\_i\}$ is a local ON frame field.


Now assume $T$ is a Codazzi operator, i.e., for any $X,Y\in \Gamma(TM)$, $(\nabla\_XT)Y=(\nabla\_YT)X$. We choose $\{e\_i\}\_{i=1}^n$ be a local orthonormal frame field of $M$ such that $\nabla\_{\star }e\_i=0$ at the considered point. For the distance function r(x) from a fixed point $x\_0$, by the definition, we have ($\nabla\_XT$ is symmetric since $T$ is symmetric)
 \begin{equation\*}
 \begin{split}
 \Delta\_{\nabla\_{\partial\_r}T}r=&\sum\_{i=1}^n\langle \nabla\_{e\_i}\partial\_r,(\nabla\_{\partial\_r}T)e\_i\rangle=\sum\_{i=1}^n\langle \nabla\_{e\_i}\partial\_r,(\nabla\_{e\_i}T)\partial\_r\rangle \\
 =&\sum\_{i=1}^ne\_i\langle \partial\_r,(\nabla\_{e\_i}T)\partial\_r\rangle -\sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}\nabla\_{e\_i}T)\partial\_r\rangle -\sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}T)(\nabla\_{e\_i}\partial\_r)\rangle .
 \end{split}
 \end{equation\*}
 However, 
 \begin{equation\*}
 \begin{split}
 \sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}T)(\nabla\_{e\_i}\partial\_r)\rangle =&\sum\_{i=1}^n\langle (\nabla\_{e\_i}T)\partial\_r,\nabla\_{e\_i}\partial\_r\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{\partial\_r}T)e\_i,\nabla\_{e\_i}\partial\_r\rangle =\Delta\_{\nabla\_{\partial\_r}T}r.
 \end{split}
 \end{equation\*}
 Hence, we obtain 
 \begin{equation}
 \begin{split}
 \Delta\_{\nabla\_{\partial\_r}T}r=\frac{1}{2}\sum\_{i=1}^ne\_i\langle \partial\_r,(\nabla\_{e\_i}T)\partial\_r\rangle -\frac{1}{2}\sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}\nabla\_{e\_i}T)\partial\_r\rangle 
 \end{split}
 \end{equation}
 We now compute the two terms of the R.H.S. of the above equality. Firstly, notice that $\nabla\_{\partial\_r}\partial\_r=0$, we have
 \begin{equation}
 \begin{split}
 \sum\_{i=1}^ne\_i\langle \partial\_r,(\nabla\_{e\_i}T)\partial\_r\rangle =&\sum\_{i=1}^ne\_i\langle \partial\_r,(\nabla\_{\partial\_r}T)e\_i\rangle =\sum\_{i=1}^ne\_i\langle (\nabla\_{\partial\_r}T)\partial\_r,e\_i\rangle \\
 =&\sum\_{i=1}^ne\_i (\partial\_r\langle T\partial\_r, e\_i\rangle )-\sum\_{i=1}^ne\_i\langle T\partial\_r,\nabla\_{\partial\_r}e\_i\rangle \\
 =&\sum\_{i=1}^n\partial\_r(e\_i\langle T\partial\_r,e\_i\rangle )-\sum\_{i=1}^n\langle T\partial\_r, \nabla\_{e\_i}\nabla\_{\partial\_r}e\_i\rangle \\ 
 =&\sum\_{i=1}^n\partial\_r\langle (\nabla\_{e\_i}T)\partial\_r,e\_i\rangle +\sum\_{i=1}^n\partial\_r\langle T\nabla\_{e\_i}\partial\_r, e\_i\rangle \\
 &+\sum\_{i=1}^n\partial\_r\langle T\partial\_r,\nabla\_{e\_i}e\_i\rangle -\sum\_{i=1}^n\langle T\partial\_r, \nabla\_{e\_i}\nabla\_{\partial\_r}e\_i\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{\partial\_r}\nabla\_{e\_i}T)\partial\_r,e\_i\rangle +\partial\_r(\Delta\_Tr)\\
 &+\sum\_{i=1}^n\langle T\partial\_r,\nabla\_{\partial\_r}\nabla\_{e\_i}e\_i\rangle -\sum\_{i=1}^n\langle T\partial\_r, \nabla\_{e\_i}\nabla\_{\partial\_r}e\_i\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{\partial\_r}\nabla\_{e\_i}T)\partial\_r,e\_i\rangle +\partial\_r(\Delta\_Tr)+\mathrm{Ric}(\partial\_r, T\partial\_r).
 \end{split}
 \end{equation}
 Secondly, 
 \begin{equation}
 \begin{split}
 \sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}\nabla\_{e\_i}T)\partial\_r\rangle =&\sum\_{i=1}^n\langle \partial\_r,\nabla\_{e\_i}((\nabla\_{e\_i}T)\partial\_r)\rangle -\sum\_{i=1}^n\langle \partial\_r,(\nabla\_{e\_i}T)\nabla\_{e\_i}\partial\_r\rangle \\
 =& \sum\_{i=1}^n\langle \partial\_r,\nabla\_{e\_i}((\nabla\_{\partial\_r}T)e\_i)\rangle -\sum\_{i=1}^n\langle \partial\_r,(\nabla\_{\nabla\_{e\_i}\partial\_r}T)e\_i\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{e\_i}\nabla\_{\partial\_r}T)e\_i-(\nabla\_{\nabla\_{e\_i}\partial\_r}T)e\_i,\partial\_r\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{\partial\_r}\nabla\_{e\_i}T)e\_i,\partial\_r\rangle -\sum\_{i=1}^n\langle (R(\partial\_r,e\_i)T)e\_i,\partial\_r\rangle \\
 =&\sum\_{i=1}^n\langle (\nabla\_{\partial\_r}\nabla\_{e\_i}T)e\_i,\partial\_r\rangle +\mathrm{Ric}(\partial\_r,T\partial\_r)-\mathrm{Ric}\_T(\partial\_r,\partial\_r).
 \end{split}
 \end{equation}
 From the above three equalities we obtain
 \begin{equation\*}
 \begin{split}
 \Delta\_{\nabla\_{\partial\_r}T}r
 =\frac{1}{2}\partial\_r(\Delta\_Tr)
 +\frac{1}{2}\mathrm{Ric}\_T(\partial\_r,\partial\_r).
 \end{split}
 \end{equation\*} 


Now, my question is that when $T=\mathrm{Id}\_{TM}$ the above equation becomes
 \begin{equation\*}
 \begin{split}
 \partial\_r(\Delta\_r)+\mathrm{Ric}(\partial\_r,\partial\_r)=0.
 \end{split}
 \end{equation\*} 
But it is well known that the Bochner formula for the distance function 
\begin{equation\*}
 \begin{split}
 |\mathrm{Hess}r|^2+\partial\_r(\Delta\_r)+\mathrm{Ric}(\partial\_r,\partial\_r)=0.
 \end{split}
 \end{equation\*} 
This obtain a contradiction. 


What is wrong with the above derivation? Thanks in advence.



What is wrong with the derivation?
Let me re-do the computation using index notation. 


Your goal is to compute $\partial\_r (\Delta\_T r)$ which we can rewrite as


$$ \nabla^c( T^{ab} \nabla^2\_{ab} r) \nabla\_c r = \nabla^c T^{ab} \nabla^2\_{ab} r \nabla\_c r + T^{ab} (\nabla^c \nabla\_a \nabla\_b r) \nabla\_c r$$


The first term on the right is $\Delta\_{\nabla\_{\partial\_r} T} r$. So we focus our attention on the second term. 


$$ \nabla^c\nabla\_a \nabla\_b r = \nabla\_a\nabla^c \nabla\_b r + [\nabla^c, \nabla\_a ] \nabla\_b r $$


So (the $\pm$ is just me forgetting which sign is the right one)


$$ T^{ab} (\nabla^c\nabla\_a \nabla\_b r) \nabla\_c r = T^{ab}(\nabla\_a \nabla^c \nabla\_b r) \nabla\_c r \pm \mathrm{Ric}\_T(\partial\_r, \partial\_r) $$


The term in the middle is


$$ \nabla\_a \nabla^c \nabla\_b r \nabla\_c r = \nabla\_a (\underbrace{\nabla^c \nabla\_b r \cdot \nabla\_c r}\_{ = 0}) - \nabla^c \nabla\_b r \nabla\_a \nabla\_c r.$$


The term that $=0$ is so due to the fact that $\nabla r$ is geodesic.
So your identity should be


$$ \partial\_r (\Delta\_T r) + \mathrm{Ric}\_T(\partial\_r, \partial\_r) = \Delta\_{\nabla\_{\partial\_r} T} r - T^{ab} g^{cd} \nabla^2\_{bd} r \nabla^2\_{ac} r $$


which now reduces correctly to Bochner's identity when $T = g$. 


**Remark** Notice that nowhere in this derivation is either of the property you started with, namely that $T$ is symmetric and $T$ is Codazzi, used. 


