Why is the Laplacian important in Riemannian geometry? As I've learned more Riemannian geometry, many of my teachers have said that studying the Laplacian (and its eigenvalues) is very important.  But I must admit, I've never fully understood why.

Fundamentally, I would like to know why the Laplacian is important among all differential operators on a Riemannian manifold.  I would also like to know what geometric information the Laplacian is supposed to encode.

That being said, I have spent a little time thinking about all this, and my current understanding is as follows:


*

*I've heard that the Laplacian is the "simplest" isometrically-invariant "scalar differential operator" on a Riemannian manifold.  If true, this statement would convince me of its importance. However, I don't know to what extent this is true.

*An isometric immersion $f \colon S \to M$ is harmonic iff it is a minimal submanifold of $M$.  In particular, an isometrically immersed submanifold of $\mathbb{R}^n$ is minimal iff its coordinate functions are harmonic.

*The Euler-Lagrange equation for the Dirichlet energy is $\Delta f = 0$.  (But why we care about minimzing energy is also somewhat mysterious to me.)

*Weitzenböck formulas comparing two elliptic second-order differential operators (and especially Laplacians) give Bochner-type vanishing theorems.
I should point out that I'm aware that harmonic functions satisfy many of the nice properties that complex-analytic functions do (by virtue of elliptic regularity and maximum principle magic).  Still, this doesn't quite tell me why I should care about the Laplace operator itself.
Note: I'm aware of this related question on the eigenvalues of the Laplacian.  But again, my interest is in Riemannian geometry; matters of applied mathematics (while interesting) are not my focus right now.
 A: I think the principal reason is Laplacian is the simplest second order elliptic operator available. So once one prove something non-trivial for the Laplacian, it is useful in other settings as well by consider a generalized Laplacian over the manifold. It is natural to consider a differential operator of any order in general(especially if one works with pseudo-differential operators), but higher order ones are not easy to study (one need techniques like Moser iteration or energy methods in general). Another reason is its association with the heat kernel, which I assume you must already know judged by the comments. 
A: Two points which I haven't seen in the other answers:
For a high-level explanation of the Laplacian's ubiquity, see Peter Stredder "Natural differential operators on Riemannian manifolds and representations of the orthogonal and special orthogonal groups" in Journal of Differential Geometry, 1975. Particularly the end of section 4.
The quasi-historical (I'm not a historian, I may be wrong) explanation is that the Laplacian on $\mathbb{R}^n$ is quite important. One high-level reason is that every constant-coefficient linear partial differential operator $L:C^\infty(\mathbb{R}^n)\to C^\infty(\mathbb{R}^n)$ which is rotationally invariant is a polynomial in the laplacian. The more honest reason (which some people might argue is a consequence of the high-level reason) is that it just came up in a number of contexts. One specific context, which you might find compelling if you think complex analysis is automatically important, is Riemann and Dirichlet's construction of harmonic functions with prescribed boundary values as the fundamental part of the proof of the Riemann mapping theorem. If you prefer physics, then it should be compelling that the Newtonian potential in gravitation is related to density by the Poisson equation $\Delta u=\rho$. Or, if you accept that the matrix of second derivatives of a function $u$ is linear-algebraically natural, then $\Delta u$ warrants its recognition as the first symmetric polynomial - along with the other symmetric polynomials (just like in Riemannian geometry).
Anyhow, if you accept that $\Delta:C^\infty(\mathbb{R}^n)\to C^\infty(\mathbb{R}^n)$ is natural, then you should be compelled that $\Delta^g$ is automatically given as its geometric version, a fact which is totally obscured in most textbooks since they start with objects like "derivations" and "parallel transport". What this means is the following:
Let $y$ be a non-standard coordinate system on $\mathbb{R}^n.$ Let $x$ denote the standard coordinates. Then two applications of the chain rule shows
\begin{align*}
\frac{\partial f}{\partial x_i}&=\frac{\partial f}{\partial y_p}\frac{\partial y_p}{\partial x_i}\\
\frac{\partial^2f}{\partial x_i^2}&=\frac{\partial^2f}{\partial y_p\partial y_q}\frac{\partial y_p}{\partial x_i}\frac{\partial y_q}{\partial x_i}+\frac{\partial f}{\partial y_p}\frac{\partial^2 y_p}{\partial x_i^2}\\
\Delta f&=\sum_{p=1}^n\sum_{q=1}^n\frac{\partial^2f}{\partial y_p\partial y_q}\sum_{i=1}^n\frac{\partial y_p}{\partial x_i}\frac{\partial y_q}{\partial x_i}+\sum_{p=1}^n\frac{\partial f}{\partial y_p}\Delta y_p.\end{align*}
Define $g_{pq}=\sum_{i=1}^n\frac{\partial x^i}{\partial y_p}\frac{\partial x^i}{\partial y_q}$ and $g^{pq}=\sum_{i=1}^n\frac{\partial y_p}{\partial x_i}\frac{\partial y_q}{\partial x_i}$, with one actually defined via the other by the observation that as matrix-valued functions, they are inverses of one another. Now a small calculation shows, with an implicit summation over all non-$r$ indices,
\begin{align*}
\frac{1}{2}g^{ij}g^{rs}\Big(\frac{\partial g_{js}}{\partial y^i}+\frac{\partial g_{is}}{\partial y^j}-\frac{\partial g_{ij}}{\partial y^s}\Big)&=\frac{\partial y_i}{\partial x_\beta}\frac{\partial y_j}{\partial x_\beta}\frac{\partial y_r}{\partial x_\gamma}\frac{\partial y_s}{\partial x_\gamma}\frac{\partial^2x_\alpha}{\partial y_i\partial y_j}\frac{\partial x_\alpha}{\partial y_s}\\
&=\frac{\partial y_i}{\partial x_\beta}\frac{\partial y_j}{\partial x_\beta}\frac{\partial y_r}{\partial x_\gamma}\frac{\partial^2x_\gamma}{\partial y_i\partial y_j}
\end{align*}
This is identical with $-\Delta y^r$, as can be seen by differentiating $\sum_{a=1}^n\frac{\partial x^i}{\partial y^a}\frac{\partial y^a}{\partial x^j}=\delta_j^i$ with respect to $x_j$ and summing over $j$. The point is that, just by considering the standard laplacian $\sum\frac{\partial^2}{\partial x_i^2}$ and changing coordinates arbitrarily, we have arrived at
$$\Delta f=g^{pq}\frac{\partial^2f}{\partial y_p\partial y_q}-\frac{1}{2}g^{ij}g^{rs}\Big(\frac{\partial g_{js}}{\partial y^i}+\frac{\partial g_{is}}{\partial y^j}-\frac{\partial g_{ij}}{\partial y^s}\Big)\frac{\partial f}{\partial y_r}$$
with $g_{ij}$ and $g^{ij}$ as defined above. If you accept that the definition of $g_{ij}$ coincides with what we call "the Euclidean metric relative to the coordinate system $(y_1,\ldots,y_n)$" then I think you must accept that if we're given a general Riemannian metric and a function $f$, then the function
$$g^{pq}\frac{\partial^2f}{\partial y_p\partial y_q}-\frac{1}{2}g^{ij}g^{rs}\Big(\frac{\partial g_{js}}{\partial y^i}+\frac{\partial g_{is}}{\partial y^j}-\frac{\partial g_{ij}}{\partial y^s}\Big)\frac{\partial f}{\partial y_r}$$
is just as fundamental as the ordinary laplacian is on $\mathbb{R}^n.$ Actually, in a sense, the above calculation shows that this is true even if you're not given a general Riemannian metric, and are content to study $\mathbb{R}^n.$ In this sense the "Riemannian laplacian" is even more fundamental than Riemannian geometry.
A: There is one important piece of evidence that has not been mentioned in the "comments" above, namely, de Rham's theorem for compact manifolds, and the idea that harmonic forms uniquely represent (real) cohomology classes. That certainly testifies to the significance of the Laplacian (the kernel of which consists of harmonic forms).
A: I've learned a few things over the last couple of years and would like to share.  Much of this has been said above, but I'd like to collect things in one place.  Here are a few big reasons the Laplacian is important to Riemannian geometry:
$(1)$ Hodge Theory.
The Hodge Theorem states that on a compact oriented Riemannian manifold $(M,g)$, one has:
$$\{\alpha \in \Omega^k(M) \colon \Delta \alpha = 0\} \cong H^k_{\text{dR}}(M).$$
This is remarkable: the left-hand side is the space of harmonic $k$-forms, which a priori depends on the metric.  But the right-hand side depends only on the homotopy type of the manifold!  This provides a link between analysis and geometry on the left, and topology on the right.
This theorem has been generalized in several directions.  For example: If $(M,g)$ is a compact Hermitian manifold, then
$$\{\alpha \in \Omega^{p,q}(M) \colon \Delta_{\overline{\partial}}\alpha = 0\} \cong H^{p,q}(M).$$
Again, the left side a priori depends on both the complex structure and the metric, whereas the right side depends only the complex structure.
As mentioned in original post, the Hodge Theorem can be coupled with Weitzenbock formulas to prove vanishing theorems.  The classic example: If $\omega$ is a harmonic $1$-form and $X = \omega^\#$ is its dual vector field, then
$$\textstyle \Delta \frac{1}{2}|X|^2 = |\nabla X|^2 + \text{Ric}(X,X).$$
The Bochner method then shows that if $(M^n,g)$ is compact, oriented, and has $\text{Ric} \geq 0$, then the harmonic $1$-forms coincide with the parallel $1$-forms, and hence (by the Hodge Theorem) $b^1(M) \leq n$.  If we had the stronger $\text{Ric} > 0$, then $b^1(M) = 0$.  This is one example of many.
$(2)$ Spectral Geometry.
Given a linear operator, we should examine its eigenvalues (and eigenspaces).
Example: On an isometrically immersed surface in $\mathbb{R}^3$, the eigenvalues of the shape operator are the principal curvatures of the surface.
Example: On an oriented Riemannian $4$-manifold, the eigenspaces of the Hodge $\ast$ operator give the decomposition of $2$-forms $\Lambda^2 = \Lambda^2_+ \oplus \Lambda^2_-$ into self-dual and anti-self-dual parts.
The eigenvalues of the Laplacian provide invariants of the Riemannian manifold, and so encode geometric information.  Now, I still don't have an intuitive grasp of this information (at the time of writing), but I've learned a few cool things.
First, on a compact Riemannian manifold $(M^n,g)$, the lowest eigenvalue $\lambda_1$ of $\Delta$ is related to the Ricci curvature.  That is, if $\text{Ric} \geq k > 0$, then (Lichnerowicz)
$$\lambda_1 \geq \frac{n}{n-1}k$$
and this estimate is sharp (Obata): equality holds iff $(M^n,g)$ is the round sphere.
Second, again on a compact Riemannian manifold $(M^n,g)$, the sequence of eigenvalues grow according to an asymptotic formula (Weyl's Law) that depends only on the dimension $n$ and the volume of $(M,g)$.
In general, the eigenvalues cannot be computed explicitly (except in very special cases), and it is estimates on the eigenvalues that provide links to geometry.  That said, the sequence of eigenvalues alone is not enough to completely determine the geometry (Milnor tori).
$(3)$ Euler-Lagrange Equations.
A significant part of Riemannian geometry is tied up with the calculus of variations -- that is, with various functionals and their critical points and extrema.  The basic examples are the volume functional and energy functionals -- forming the theory of geodesics, minimal submanifolds, and harmonic maps -- but there are many others beyond these.
The critical points of a functional are described by "Euler-Lagrange equations," and these equations frequently involve the Laplacian.  For example, if $F \colon (M,g) \to (N,h)$ is an isometric immersion, then the condition that $F(M) \subset N$ have zero mean curvature is the PDE system given in local coordinates by
$$\Delta_M F^j - \sum_{i,k} \sum_{\alpha, \beta} g^{\alpha \beta}(x) \Gamma^j_{ik}(F(x)) \frac{\partial F^i}{\partial x^\alpha} \frac{\partial F^k}{\partial x^\beta} = 0,$$
for $1 \leq j \leq \dim(N)$.

Remarks
(a) The Laplacian can be placed in the more general context of (squares of) Dirac operators.  I'd like to say more about this some other time.
(b) The condition $\Delta \alpha = 0$ is equivalent to having both $d\alpha = 0$ and $d^*\alpha = 0$.  This recasts the 2nd-order equation $\Delta \alpha = 0$ as a 1st-order system.
This "closed and co-closed" pairing seems to come up a lot.  Examples include the Yang-Mills equation (for connections) and the integrability condition for $\text{G}_2$-structures.
(c) Even in the setting of flat $\mathbb{R}^n$, harmonic functions (solutions to $\Delta u = 0$) are pretty interesting.  For one thing, solving $\Delta u = 0$ implies that the function $u$ enjoys many properties similar to holomorphic functions:


*

*Mean Value Property

*An integral representation formula

*Liouville-type theorems

*Maximum Principle

*Regularity results


This suggests that many of the nice properties of holomorphic functions in $\mathbb{C}$ are really artifacts of harmonicity.
For another thing, solutions of $\Delta u = 0$ are exactly the critical points of the Dirichlet energy functional.  This links harmonic functions to the calculus of variations.  More geometrically, it suggests a link between harmonic functions and minimal submanifolds.  There are several results which realize this link, but I'll leave it at that.
