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.