What is the significance of differential operators over other operators in group theory?
Good question.
The Laplace-Beltrami on $\mathbb S^n$ is not "any" operator, of course. It is related to the group $SO(n)$, like the Euclidean Laplacian is related to the group of isometries of $\mathbb R^n$. The precise mathematical notion that captures this should be that of Casimir operator, but that is something that I do not understand, I must admit.
So I will give you my down-to-earth point of view. The starting point is Schur's lemma of representation theory, according to which, if $\phi\colon V\to V$ is an intertwining map of an irreducible representation of a group, then $\phi$ is a scalar multiple of the identity operator.
So what does it mean? It means that a linear operator $\phi$, that commutes with the action of a group, will be diagonalized by the decomposition of the Hilbert space in irreducible representations (irreps).
And that's precisely what the various Laplacians do: the Laplace-Beltrami commutes with $SO(n)$, therefore the spherical harmonics must be its eigenfunctions, because they are the basis functions of the irreps of $SO(n)$. Similarly, the Euclidean Laplacian commutes with translations, so the plane waves $e_\mathbf{p}(\mathbf{x}):=\exp(i\, \mathbf{x}\cdot\mathbf{ p})$ must be its eigenfunctions.
EDIT (after comments). I have deliberately been sloppy here, because there's a number of technical points that I do not want to address (mostly because I can't, to be honest).
- I speak of "action of the group", while the precise term is "representation", as you point out. Which representation, and on which Hilbert space? For both the Laplace-Beltrami and the Euclidean Laplacian, the Hilbert space is $L^2(X)$, with $X=\mathbb S^n$ and $X=\mathbb R^n$, respectively. The group, which we denote by $G$, is $G=SO(n)$ and $G=(\mathbb R^n ,+)$, respectively. And the representation is always given by $$ \rho(g) f (x)=f(g^{-1}\cdot x),$$ where $g^{-1}\cdot x$ reads $x-g$ in the case of $G=(\mathbb R^n, +)$. To say that the Laplacian is an intertwining map means that $$\Delta\rho(g)=\rho(g)\Delta.$$
- Schur's lemma says that if $\mathcal H\subset L^2(X)$ is an irreducible subspace for $\rho$ then there exists a scalar $\lambda$ such that $\left.\Delta\right|_{\mathcal H} =\lambda I$. In other words, $$ \Delta f = \lambda f,\qquad \forall f\in\mathcal H.$$ Notice that the eigenvalue $\lambda$ is dependent on $\mathcal H$. Actually, each irreducible subspace is labeled by these eigenvalues. We have therefore realized the irreducible subspaces, that are a rather slippery concept (see definition below), as eigenspaces of the Laplacian, which is much more concrete. Recall that irreducible subspace means that $\rho(g)\mathcal H \subset \mathcal H$ for all $g\in G$ and that if $\mathcal H'$ is a subspace of $\mathcal H$ with the same property then $\mathcal H'$ is either $\{0\}$ or coincides with $\mathcal H$.
I am sweeping something under the rug. Indeed, I am acting as if $\Delta$ were defined on $L^2(X)$, which is not the case: there exist functions whose Laplacian is not defined, or it is not an element of $L^2(X)$. This is especially bad in the case $X=\mathbb R^n$, in which case even the "eigenfunctions" that we discover are not eigenfunctions at all. I am talking about the plane waves $e_\mathbf{p}$, whose $L^2$ norm is $\infty$. There is a lot of theory to study to fully understand how to go around these technical difficulties.
But, whatever, I think that what we did in this post is enough as an illuminating guiding principle.