Let $\mathbb{S}=\{x\in\mathbb{R}^n|x_1^2+\ldots +x_n^2=1\}$ be the unit sphere in $\mathbb{R}^n$, $\mathbb{C}[x]=\mathbb{C}[x_1,\ldots ,x_n]$ the complex-valued polynomial functions on $\mathbb{R}^n$, and $\Delta=\partial_1^2+\ldots +\partial_n^2$ the Laplacian. Let $H_n=\ker(\Delta)$ be the space of harmonic polynomials. Consider the restriction to the sphere
$$\rho:\mathbb{C}[x_1,\ldots ,x_n]\to\mathbb{C}[x_1,\ldots ,x_n]/(x_1^2+\ldots +x_n^2-1)$$
and the image $\bar{H}_n=\rho(H_n)\subseteq\mathbb{C}[x_1,\ldots ,x_n]/(x_1^2+\ldots +x_n^2-1)$.
Now, there is the curious fact that $\bar{H}_n$ is in fact the whole of $\mathbb{C}[x_1,\ldots ,x_n]/(x_1^2+\ldots +x_n^2-1)$. In other words, for every $f\in\mathbb{C}[x_1,\ldots ,x_n]$, the polynomial restriction $f|_{\mathbb{S}}$ on the sphere has a harmonic representative $h\in H_n$ such that $f|_{\mathbb{S}}=h|_{\mathbb{S}}$.
This follows from a $\sum(\mathrm{radial})\cdot (\mathrm{harmonic})$ decomposition
Theorem. Let $f$ be a polynomial. Then $f(x)=\sum_i f_i(x_1^2+\ldots +x_n^2)\cdot h_i(x)$ where $h_i\in H_n$ and $f_i\in\mathbb{C}[t]$.
Indeed: $f|_\mathbb S =(\sum_i f_i(x_1^2+\ldots+x_n^2) \cdot h_i)|_\mathbb S=\sum_i f_i(1)\cdot h_i|_\mathbb S = (\sum_i f_i(1)h_i)|_\mathbb S$ and now $h:=\sum_i f_i(1)h_i$ is harmonic.
In particular, harmonic restrictions on the sphere form an algebra (equal to the algebra of all polynomial restrictions). From this, I think, one can prove the spherical harmonics are a Hilbert basis on the sphere by using Stone-Weierstrass (and density of continuous functions in $L^2$) without proving that they form a basis of eigenfunctions of the Laplace-Beltrami operator $\Delta_\mathbb{S}$ on $\mathbb S$ (which by the way they do, and it's kinda interesting in itself).
Q. Is this stuff part of some theory, or just a coincidence for the usual Laplacian?
For example, what if $D$ is a differential operator with symbol $\sigma_D(x)$ and $X_D\subseteq\mathbb{R}^n$ the variety $\{x|\sigma_D(x)=1\}$. Is there sometimes a polynomial decomposition
$$f=\sum_i f_i(\sigma_D(x))\cdot h_i(x)$$
with $Dh_i=0$? Or an "intrinsic" operator $P$ on $X_D$, related to $D$ in a similar way as $\Delta_\mathbb{S}$ is related to $\Delta$, such that the homogeneous $h\in\ker D$ are eigenfunctions of $P$?