Harmonic polynomials on the sphere 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$?
 A: There is a similar separation of the variables for ellipsoids
$x^2/a^2+y^2/b^2+z^2/c^2=1$ related to ellipsoidal harmonic functions,
see Whittaker Watson, Chap. 23. There is a generalization to
arbitrary dimension, see, for example
https://arxiv.org/pdf/math/0610718.pdf
A: I view this as a concatenation of two facts:

Fact 1: Let $k$ be a field, let $I$ be an ideal of $k[x_1, \ldots, x_n]$ and let $J$ be associated graded ideal of $I$, meaning that a degree $d$ homogenous polynomial $g(x)$ is in $J$ if and only if there is an element of $I$ of the form $g(x) + (\text{terms of degree $< d$})$. Let $\{ p_a \}$ be a set of homogenous polynomials which maps to a basis of $k[x]/J$. Then $\{ p_a \}$ also maps to a basis of $k[x]/I$.
The proof of Fact 1 is just an upper triangularity argument.
Apply Fact 1 with $I$ the ideal $\langle \sum x_i^2 -1 \rangle$ and $J$ the ideal $\langle \sum x_i^2 \rangle$. So this reduces us to the question of why harmonic polynomials are a basis for $\mathbb{R}[x]/\langle \sum x_i^2 \rangle$.

Turn $\mathbb{R}[x_1, x_2, \ldots, x_n]$ into a module over itself by 
$$g(x_1, \ldots, x_n) \ast f(x_1, \ldots, x_n) := g(\tfrac{\partial}{\partial x_1}, \ldots, \tfrac{\partial}{\partial x_n}){\big(} f(x_1, \ldots, x_n) {\big)}$$
To be clear, this is a differential operator, encoded by $g$, acting on the polynomial $f$.
This action induces an inner product on the degree $d$ homogenous polynomials by $\langle g,f \rangle := g \ast f$, since $g \ast f$ is a degree $0$ polynomial and can be thought of an element of $\mathbb{R}$.
Let $J \subseteq \mathbb{R}[x_1, \ldots, x_n]$ be any graded ideal of $\mathbb{R}[x_1, \ldots, x_n]$. We write $J = \bigoplus J_d$ for its composition into graded pieces. Let $J^{\perp} = \bigoplus J_d^{\perp}$, where the orthogonal complement is taken with the above inner product on $\mathbb{R}[x_1, \ldots, x_n]_d$. By generalities of inner products, the map $J^{\perp} \hookrightarrow R[x_1, \ldots, x_n] \twoheadrightarrow R[x_1,\ldots, x_n]/J$ is an isomorphism.
Fact 2 We have $$J^{\perp} = \{ f : g \ast f=0 \ \forall g \in J \}.$$ Moreover, if $g_1$, $g_2$, ..., $g_N$ is a list of generators for $J$, it is sufficient to impose that $g_i \ast f=0$ for each generator $g_i$.
Thus, the set of polynomials with $\sum (\tfrac{\partial}{\partial x_i})^2 f =0$ is a basis for $\mathbb{R}[x]/(\sum x_i^2)$.

I learned this by reading Mark Haiman's papers on "diagonal harmonics", and trying to figure out why he was calling them harmonic functions. I'll try to find some better references.
