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.
Best Answer
The answer is positive. Consider a small disk $D(a,r)$ in the domain of uniform convergence (the center is at $a$, radius $r$), and expand $f_n(a+te^{i\theta})$ into a harmonic Fourier series: $$f_k(a+te^{i\theta})=\sum_{n=-\infty}^\infty c_{k,n}t^ne^{in\theta},\quad t<r.$$ Then we have $c_{k,n}\to c_n$, as $k\to\infty$ for each $n$, in view of uniform convergence. Therefore, the analytic parts $$g_k(a+te^{i\theta})=\sum_{n=0}^\infty c_{k,n}te^{in\theta}$$ converge uniformly on each smaller disk. Then $h_k=\overline{f_k-g_k}$ also converge.
Remark. The Hilbert space $L^2$ of functions on the circle $\partial D(a,r)$ is a direct sum $L^2=H^2_++H^2_-$, where $H^2_+$ and $H^2_-$ are "analytic" and "anti-analytic" parts. So the map $f\mapsto g$ is just a projection onto $H^2_+$. This gives an explicit estimate of norms of $h_n$ in terms of norms of $f$.