Why is a diffeomorphism an isometry if and only if it commutes with the Laplacian? I came across the following statement in a book on automorphic forms:

In general, on a Riemannian manifold, the Laplace-Beltrami operator $\Delta$ is characterized by the property that a diffeomorphism is an isometry if and only if it leaves $\Delta$ invariant.

For this question, I am interested in the last part of the statement, that "a diffeomorphism is an isometry if and only if it leaves $\Delta$ invariant." What is a precise statement of this result? We're viewing $\Delta$ as an operator on $C^\infty(M)$, right? Is there a way to prove it just from the description of $\Delta$ in coordinates:
$$
\frac{-1}{\sqrt{\det{g}}}\sum_{i,j}\partial_j(g^{ij}\sqrt{\det{g}}\partial_i f)?
$$
I'd also be interested in a way to prove it from the definition of $\Delta$ as $-\text{div}\circ\nabla$.
 A: Let $\phi: M \to N$ and suppose $(M,g)$ and $(N,h)$ are Riemannian manifolds.  If $\phi^* \circ \Delta_h = \Delta_g \circ \phi^*$, then, as Jack Lee indicated, $\phi$ is an isometry. Here is an elementary and complete proof using the $-div \circ \nabla$ point of view:
By hypothesis, for each pair of smooth functions
$u,v: N \to {\mathbb R}$ we have
\begin{eqnarray*}
%
\int_M (\Delta_g (u \circ \phi)) \cdot (v \circ \phi) \, dV_g
&=&
\int_M \left((\Delta_h u) \circ \phi \right) \cdot (v \circ \phi)  \, dV_g \\
&=&
\int_N (\Delta_h u)  \cdot v \cdot  (\phi^{-1})^{*}(dV_g) \\
&=&
\int_N (\Delta_h u)  \cdot  v \cdot  \rho\, dV_h \\
%
\end{eqnarray*}
where $(\phi^{-1})^{*}(dV_g)  = \rho \cdot dV_h$.
If we assume that $u$ is compactly supported, then by integrating
each side by parts we obtain
$$  
\int_M \bar{g} (d( u \circ \phi), d(v \circ \phi) ) \, dV_g~
=~
\int_N \bar{h} (d u, d(v \cdot \rho) ) \, dV_h.
$$
where $\bar{g}$ and $\bar{h}$ are the co-metrics associated to $g$ and $h$ respectively.
By changing variables in the integral on the left hand side, we obtain
$$  
\int_N (\phi^{-1})^*(\bar{g}) (du, dv )\cdot \rho \cdot dV_h~
=~
\int_N \bar{h} (d u, d(v \cdot \rho) ) \, dV_h.
$$
If we choose $v=1$, then the left hand side equals zero and so
$$ 
0~
=~
\int_N \bar{h} (d u, d \rho ) \, dV_h.
$$
Since $u$ is an arbitrary smooth function, this implies
that $d \rho=0$ and so $\rho$ is constant.
Therefore we may divide both sides of the previous identity by
$\rho$ and obtain
$$ 
\int_N (\phi^{-1})^*(\bar{g}) (du, dv )\cdot dV_h~
=~
\int_N \bar{h} (d u, dv )\cdot dV_h.
$$
It follows that $(\phi^{-1})^*(\bar{g})= \bar{h}$ and so $\phi^*(h)=g$.
A: Here's one way to make it precise.
Theorem. Suppose $(M_1,g_1)$ and $(M_2,g_2)$ are Riemannian manifolds, $\Delta_1,\Delta_2$ are their respective Laplace operators, and $\phi\colon M_1\to M_2$ is a diffeomorphism. Then $\phi$ is an isometry if and only if for every $f\in C^\infty(M_2)$ we have $\Delta_1(f\circ\phi) = (\Delta_2 f)\circ \phi$.
(Of course, if you're only interested in diffeomorphisms from a Riemannian manifold to itself, you can take $M_1=M_2=M$ and $g_1=g_2=g$.)
I don't have time to write down a complete proof, but the idea is that the principal symbol of a differential operator is a diffeomorphism-invariant quadratic function on the cotangent bundle, and the principal symbol of the Laplace operator is given by the squared norm function on covectors.  By the polarization identity, the squared norm function determines the metric on covectors (the "dual metric"), and by the musical isomorphism between $TM$ and $T^*M$, this in turn determines the metric.
