This is a more advanced question than Estimates of $\Delta|\nabla u|$ for harmonic function $u$ . The Bochner formula and refined Kato inequality tells us that on a Riemannian manifold $(M^{n},g)$, if $u$ is a harmonic function in $B_{2}(x)\subset (M^{n},g)$, then $$ \frac{1}{n}|\nabla^{2}u|^{2}\leqslant |\nabla^{2}u|^{2}-|\nabla|\nabla u||^{2}=|\nabla u|\Delta|\nabla u|-Ric(\nabla u,\nabla u)\leqslant |\nabla^{2}u|^{2}. $$ My question is how to obtain an analogous $L^{p}$ estimate for $\Delta|\nabla u|$. Let me fully state my question below.
For any $\delta>0$, does there exist $\epsilon>0$ such that if $(M^{n},g)$ is a manifold with $|Ric|\leqslant\epsilon$ and $u$ is a harmonic function in $B_{2}(x)\subset (M^{n},g)$ with $$ \begin{split} &|\nabla u|\leqslant 1,\\ &\int_{B_{2}}||\nabla u|-1|<\epsilon, \\ &\int_{B_{2}}|\nabla^{2}u|^{2p} < \epsilon \end{split} $$ for some integer $p\geqslant1$, then $$ \int_{B_{1}}(\Delta|\nabla u|)^{p}<\delta \ ? $$ This problem has been solved in the case of Euclidean space. See Estimates of $\Delta|\nabla u|$ for harmonic function $u$ . But on a manifold with $|Ric|\leqslant\epsilon$, that method seems to not work, since we cannot get harmonicity for $\nabla u$ and apply Harnack inequality.
Also, it is easy when $p=1$, since we can multiply the Cheeger-Colding’s cutoff function (which has bounded laplacian provided lower bounds on Ricci, see https://doi.org/10.2307/2118589) and then integrate by parts. But for integer $p\geqslant2$, I can’t see how to do it.
