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 $\leq d$})$$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.