I have been trying to understand more about the norm squared of the Hessian on a Riemannian manifold, that is $$ |\nabla^2f|^2 $$ This quantity shows up in the Bochner formula, for instance. On $\mathbb{R}^n$, the induced fiber metric on the 2-tensor bundle is just a dot product of matrices (viewed as $\mathbb{R}^{n\times n}$ in this case), so viewing the contribution from the diagonal terms, we get $$ \sum_{i=1}^n\frac{\partial^2 f}{\partial x_i^2} $$ These can be used to estimate the square of the Laplacian via the Cauchy-Schwartz inequality.
I'm wondering how to show that we have an analogous property in the more general case, that is, showing that the norm of the Hessian squared is equal to a sum of terms that resemble the Laplacian in coordinates squared plus some other non-negative terms.
So far, I have tried the following: take a symmetric two-tensor $T$ on a manifold $(M^n,g)$; it induces a self-adjoint operator $A$ satisfying $T(v,w)=\langle A(v),w\rangle$. Fix a point $p\in M$ and choose a local frame about $p$, call it $E_j$, which diagonalizes $A$. Using the canonical association between endomorphisms and $(1,1)$ tensors, we can write $$ A=\lambda_i \epsilon^i\otimes E_i $$ where there is a sum in $i$ and $\epsilon^i$ is the coframe dual to $E_i$. Then $$ \langle A,A\rangle =\sum_{i,j}\lambda_i\lambda_j\langle\epsilon^i,\epsilon^j\rangle\langle E_i,E_j\rangle $$ Here, I want to claim that the last quantity is a sum of squares of eigenvalues plus a non-negative error term, but I think I am missing some relevant identity from linear algebra. The sum of squares of eigenvalues would then specialize in our Hessian case to be related to the Laplacian terms (I think). I'm interested in advice either for how to make this approach work, or for a different approach to thinking about this that is simpler.