Why is the Laplacian ubiquitous? The title says it all. 
I'm wondering why the Laplacian appears everywhere, e.g. number theory, Riemannian geometry, quantum mechanics, and representation theory. And people seems to care about their eigenvalues. Is there any deep reason behind this?
 A: To build on what Andrey has mensioned about Brownian motion. The Laplacian is one of the points of connection between stochastic processes and analysis. The Laplacian appears as the infinitesimal generator of Brownian motion and conversely a self-adjoint operator that has some of the properties of the Laplacian can be used to define a `Brownian motion' on spaces other than $\mathbb{R}^{d}$. For example one of the early proofs of  the existence of a Brownian motion on the Sierpinski carpet centers on first creating a Laplacian on the Spierpinski carpet. 
In this context the eigenvalues of the Laplacian can be used to determine the properties of the heat kernel associated to the Brownian motion, like continuity or Gaussian/sub-Gaussian/or other bounds. 
I hope this gives you some part of what you asked for. 
A: The Laplacian of a function $u$ at a point $x$ measures the average extent to which the value of $u$ at $x$ deviates from the value of $u$ at nearby points to $x$ (cf. the mean value theorem for harmonic functions).  As such, it naturally occurs in any system in which some quantity at a point is influenced by the value of the same quantity at nearby points.  (This also explains the link between the Laplacian and Brownian motion (or random walks), in which one repeatedly travels from a point $x$ to a randomly selected nearby point to $x$.)
The notion of "nearby points" requires only that one have a Riemannian metric structure on the underlying space, and so the Laplacian is a natural Riemannian invariant (and also a conformal invariant in two dimensions).  This makes it a useful operator-theoretic proxy for Riemannian or conformal structure, in particular allowing one to use spectral theory to start controlling the Riemannian or conformal geometry of a domain.  Being invariant, it also has a good chance of commuting (or nearly commuting) with other interesting invariant operators (e.g. partial derivatives, Hecke operators, heat or wave operators, etc.), and so the spectral theory of the Laplacian is often useful in illuminating the spectral theory and dynamics of many other invariant objects.
A: The importance of the Laplacian is a reflection of the importance of Riemannian geometry, both for its own sake and in these other fields.  (Obviously, absent a generalization, you can't have a Laplacian without a Riemannian metric.)  In Riemannian geometry, the Laplacian is the first scalar linear differential operator available which is "covariant", i.e., that depends only on the Riemannian structure and not on extra choices such as coordinates.
It's natural for the first non-trivial possible structure of a given type to be fundamental.  One reason is that it can be an approximation to something else with higher-order terms.  For example, in the wave equation, if it is meant as a realistic model of sound waves, actually there are all kinds of higher order, non-linear effects; but one begins with the Laplace operator as the correct approximation for small waves.  It has to be correct because it's the only one available.  The same thing happens with heat and the heat equation.
Actually, something interesting happens if you relax the use of the word "scalar".  If you have a spin manifold, then there is a Dirac operator, which is a Laplacian-like operator that is just as important.  But even apart from that construction, there are various other operators, such as the Hodge Laplacian, that act on vector and tensor fields rather than on scalar functions.  These are not exactly the same as the original Laplacian, but they tend to get the same name.
A: The Laplacian is a second order differential operator which is both linear
and invariant under any rigid body motion (rotations and translations). 
There is no simpler non-trivial operator with these properties. In modeling many
physical phenomena, one expects these two properties to arise, so sometimes
the simplest approximation based on physical principles leads to the Laplacian.
A: In Physics, the essential reason boils down to symmetry. One expects the fundamental laws of physics to be independent of where you are or how you are oriented in space, and if they are described by scalar fields then there are good reasons to expect these fields to be extrema of a first order Lagrangian  that also has these symmetries, i.e., that are invariant under Euclidean motions. This leads naturally (see Deane Yang's answer) to field equations that involve the Laplacian. This is in particular how the field equations for the electric and (Newtonian) gravitational potential arise. (Of course this is a VERY condensed answer---a full version would take many pages.)
A: The Laplace operator is to Analysis and PDEs is (almost) what a sum of squares is to Linear Algebra and Statistics.


*

*The Laplacian is the simplest differential quadratic form corresponding via the Fourier transform to the square of the Euclidean distance. This may explain in part its fundamental role in the Harmonic Analysis on Euclidean spaces.

*The Laplacian (or, more precisely, $\frac{1}{2}\Delta$) is the infinitesimal generator of a Brownian motion on $\mathbb R^n$, which is the simplest and most ubiquitous of the continuous-time stochastic  processes. 

*The Laplace-Beltrami operator on a Riemannian manifold is conformal invariant. 
The connection between harmonic functions, complex analysis and probability is the most tight in dimension 2 (Levy's theorem, Schramm-Loewner evolution, ...).

*The Laplace operator is the trace of the coefficient of the quadratic term in a local Taylor expansion of a function (the Hessian matrix). This implies that it will pop up (together with the determinant of the Hessian matrix) in many problems related to optimization. 
A: My take on this question: it's not the operator, it's the solutions. More explicitly, the solutions of $\Delta u=0$ enjoy the fundamental property that $u$ at a point is the average of the values of $u$ at surrounding points. You will agree that this property is very physical, and sounds reasonable for many quantities which tend to diffuse: heat, tension, energy, etc. And not only physics, also pure math: minimal surfaces, discrete variants which have to do with averaging processes, and so on. Basically, every time there is a quantity which diffuses at a fast time scale, so you may speak of an 'equilibrium', you see that this averaging property pops out (or something similar).
Now, the average condition is nonlocal, and the way to express it as a local (differential) equation is precisely by saying that $\Delta u=0$. Add this to all the already mentioned properties (symmetry, positivity, being the nicest possible Fourier multiplier, etc etc) and you will be slightly less surprised.
A: Being the question very ample, its ampleness reflects itself in the possible answers.
So, even if trivial, I give a basic observation just in order to compose the complessive frame:
A linear differential operator L on the Schwartz space in $\mathbb{R}^n$ is invariant under the euclidean group, i.e. $L\circ {R^{\ast}}={R^{\ast}}\circ L$, for any $R\in E(n)$, if and only if there is a polynomial in one indeterminate $P(T)$ such that $L=P(\Delta)$.
For a reference see Theorem 8.51 in Folland G., Real Analysis, 2nd Ed..
Excuse me if this observation is too elementary.
A: It is the Euler-Lagrange operator for the $L^2$ norm of the gradient of a function, curve, vector field, whatever, which turns out to be a natural and important functional in physics, geometry, and many other settings.
