Let $\mathcal{M}$ be a compact Riemannian manifold and let $\Delta$ be the (scalar) Laplace-Beltrami operator on $\mathcal{M}$. Then $\Delta$ has a discrete spectrum and if we order its distinct eigenvalues $\lambda_i$ by magnitude then some very simple examples suggest that the magnitude of $\lambda_i$ might be roughly quadratic in $i$. For instance, on the circle $S^1$ eigenfunctions have the form $\cos(nx)$ or $\sin(nx)$ for $n \in \mathbb{N}_0$; hitting these functions with $-\frac{\partial}{\partial x^2}$ yields $n^2\cos(nx)$ and $n^2\sin(nx)$, respectively. Similar analysis can be done for the geometrically flat torus $T^2$. On the 2-sphere, we have $\lambda_i=i(i+1)$. This rough idea of "differentiating twice leads to a square" makes me suspect that a similar relationship might hold for other domains -- what can be said in general? I'm particularly interested in smooth surfaces embedded in $\mathbb{R}^3$.
Update: Weyl's formula provides some valuable information about the Laplace spectrum, but does not determine the asymptotic growth of $\lambda_i$. For instance, suppose we have a manifold such that $N(R) \approx R$, i.e., the number of eigenvalues with value no greater than $R$ is roughly equal to $R$ itself. Letting $n_i$ be the multiplicity of $\lambda_i$, this relationship holds for, say, $\lambda_i = i(i+1)$ and $n_i = 2i+1$ (which is the situation on the sphere), since
$$ N(\lambda_i) = \sum_{j=0}^i n_j = \sum_{j=0}^i 2i+1 = i^2 + 2i + 1 \approx i(i+1) = \lambda_i. $$
But it also holds for $\lambda_i = i$ and $n_i = 1$ since then
$$ N(\lambda_i) = \sum_{j=0}^i n_j = \sum_{j=0}^i 1 = i + 1 \approx i = \lambda_i. $$
