Let $M$ be a Riemannian manifold. When does the Laplace-Beltrami operator differ from the trace of the Hessian? Let $(M^n,g)$ be some Riemannian manifold, and let $(x^1,\ldots,x^n)$ be local (oriented) coordinates. Then we may define the Laplace-Beltrami operator $\Delta$ on a function $f\in C^\infty(M)$ by setting $$\Delta f = d^* df = -\operatorname\ast d\operatorname\ast df = -\frac{1}{\sqrt{\det(g_{ij})}}\sum_{j=1}^n\frac{\partial}{\partial x^j}\left(\mathrm{grad}^j\, f\sqrt{\det(g_{ij})}\right)$$ while we can also consider minus the trace of the Hessian,
$$ -\operatorname{tr}\operatorname{Hess}(f) = -\operatorname{tr}\nabla^2 f = -g^{ij}\left(\frac{\partial^2 f}{\partial x^i\partial x^j} - \Gamma_{ij}^k\frac{\partial f}{\partial x^k}\right).
$$
Now, on flat space, these are exactly identical. But are they the same in general? They both seem to be fairly natural generalizations of the notion of the Laplacian on $\Bbb R^n$.
 A: They are identical. Consider the "divergence formula" used often in Ricci-calculus, given by $$ \nabla_\mu X^\mu=\frac{1}{\sqrt{\det g}}\partial_\mu\left(X^\mu\sqrt{\det(g)}\right). $$
Now let $X^\mu=\nabla^\mu f$. Then we have $$ \nabla_\mu\nabla^\mu f= \frac{1}{\sqrt{\det g}}\partial_\mu\left(\nabla^\mu f\sqrt{\det(g)}\right),$$ which is the trace of the covariant Hessian, and is the same as minus your expression for the Laplace-Beltrami operator.

Edit: About the proof of the divergence formula:
The connection coefficients for the Levi-Civita connection are given by $$ \Gamma^\sigma_{\mu\nu}=\frac{1}{2}g^{\sigma\lambda}(\partial_\mu g_{\nu\lambda}+\partial_\nu g_{\mu\lambda}-\partial_\lambda g_{\mu\nu}). $$
The divergence (with respect to the connection) is given by $$ \nabla_\mu X^\mu=\partial_\mu X^\mu+\Gamma^{\mu}_{\mu\nu}X^\nu, $$ so we need to caclulate the contraction of the connection coefficients:
$$ \Gamma^\mu_{\mu\nu}=\frac{1}{2}g^{\mu\lambda}(\partial_{\mu}g_{\nu\lambda}+\partial_\nu g_{\mu\lambda}-\partial_\lambda g_{\mu\nu}). $$ Due to symmetry, the first and third terms in the brackets cancel, so we're left with $$ \Gamma^\mu_{\mu\nu}=\frac{1}{2}g^{\mu\lambda}\partial_\nu g_{\mu\lambda}. $$
Recall Jacobi's formula (proof is on the wiki page): For a differentiable function $A:\mathbb{R}\rightarrow GL(n,\mathbb{R})$
, we have $$ \frac{d}{dt}\det(A)=\det(A)\text{Tr}\left(A^{-1}\frac{dA}{dt}\right). $$
We apply this formula to the metric tensor with $g=(g_{\mu\nu})$
, so we have $$ \partial_\nu\det(g)=\det(g)\text{Tr}(g^{-1}\partial_\nu g)=\det(g)g^{\mu\lambda}\partial_\nu g_{\mu\lambda}=2\det(g)\Gamma^\mu_{\mu\nu}. $$
We divide by $2\det(g)$: $$ \frac{1}{2}\partial_\nu\ln\det(g)=\Gamma^\mu_{\mu\nu}=\partial_\nu\ln(\sqrt{\det(g)})=\frac{1}{\sqrt{\det(g)}}\partial_\nu\sqrt{\det(g)}. $$
We now insert this expression for the contracted connection coefficients into the expression for the divergence:
$$ \nabla_\mu X^\mu=\partial_\mu X^\mu+\frac{1}{\sqrt{\det(g)}}\partial_\mu\sqrt{\det(g)}X^\mu=\frac{1}{\sqrt{\det(g)}}\left(\sqrt{\det(g)}\partial_\mu X^\mu+\partial_\mu\sqrt{\det(g)}X^\mu\right) \\ =\frac{1}{\sqrt{\det g}}\partial_\mu(\sqrt{\det g}X^\mu).$$
A: This matter is naturally resolved by using an orthonormal frame, rather than a coordinate frame. To be more precise, let $p\in M$, and let $e_1,\ldots,e_n\in T_pM$ be an orthonormal basis. Extend this basis to a local frame by parallel transport along geodesic rays from $p$. This way, we obtain a local orthonormal frame which satisfies $\nabla e_i|_p=0,\quad i=1,\ldots,n$. Let $\eta^1,\ldots,\eta^n$ denote the dual frame of the cotangent bundle.
Let $f:M\to\mathbb{R}$. Then $df=e_i(f)\eta^i,$ and so, $$*df=(-1)^{i+1}e_i(f)\eta^1\wedge\ldots\widehat{\eta^i}\ldots\wedge\eta^n.$$As each $\eta^i$ satisfies $\nabla\eta^i|_p=0$, we also have $d\eta^i|_p=0$, and so, by the Leibniz rule, we have at $p$ $$d*df=\sum e_i(e_i(f))\eta^1\wedge\ldots\wedge\eta^n=\sum e_i(e_i(f))\mathrm{vol},$$which finally yields$$*d*df|_p=\sum e_i(e_i(f)).$$On the other hand, for the Hessian we have$$\nabla^2f|_p(e_i,e_j)=\nabla_{e_i} df(e_j)=e_i(e_j(f))-df(\nabla_{e_i}e_j)=e_i(e_j(f)),$$or in other words,$$\nabla^2f|_p=e_i(e_j(f))\eta^i\otimes\eta^j.$$The trace is thus equal to $\sum e_i(e_i(f)),$ and we are done.
