Is every harmonic polynomial a linear combination of these?

This is a theorem of Maxwell and Sylvester, see for example https://arxiv.org/abs/math-ph/0408046 , https://arxiv.org/abs/0805.1904 and https://arxiv.org/abs/astro-ph/0412231 but, it is in 3d only.

In higher dimension, this is the object of a conjecture of Shubin, see page 10 of https://arxiv.org/abs/0704.1174

Great question BTW, related to some interesting and beautiful mathematics.

I think there is an error in my calculations in my comments. The following is not a complete answer but it gives an hint on how to proceed.

Define $P(N,K)$ to be the space of homogeneous polynomials of degree $K$ in $N$ variables. It has dimension $\binom{N+K-1}{N-1}=\binom{N+K-1}{K}$ the way of distributing the degree $K$ to the $N$ variables. The space $H(N,K) \subset P(N,K)$ is the subspace of harmonic polynomials. One can show that $\dim H(N, K) = \dim P(N, K) - \dim P(N, K-2) $, because the direct sum of the small spaces is isomorphic to the big space.

Note that the expression you wrote is a map from homogeneous polynomials of degree $K$ to harmonic polynomials of degree $K$ by linear extension. For example, in the case $N=5, K=2$ we have

$$ xy +yz \mapsto r^7\nabla_x \nabla_y \frac{1}{r^3} + r^7 \nabla_y \nabla_z \frac{1}{r^3} $$

Let us call this map $\rho: P(N, K) \to H(N, K) $. Now the relationship you wrote in the comments $\sum_a \nabla_a \nabla_a \frac{1}{r^{\bullet}} = 0$ rewrites as $ \rho(r^2) =0$, where $ r^2 = x_1^2 +\ldots + x_N^2$ is the "norm" polynomial.

We consider the map $r^2: P(N, K-2) \to P(N, K) $ that multiplies a polynomial by $r^2$. By a similar argument, we get that the composition of the two maps

$$ P(N, K-2) \underset{r^2}{\to} P(N, K) \underset{\rho}{\to} H(N, K) $$

is zero. Also, the dimension of the first and the third space sum to the dimension of the middle space. This means that the second map is surjective if and only if the kernel is the image of the first map. Here I can't go on. In other words, we have to show that the only relations that can occur are the ones that we mentioned: contracting two indices and permuting derivatives. I think it is remarkable to mention that $P(N, K) $ writes as the direct sum of $H(N, K) $ and the image of the first map, $r^2 P(N, K-2) $; the thesis is thus equivalent to the fact that every harmonic polynomial $f$ is uniquely expressible as $\rho(f_v) $ for some harmonic polynomial $f_v$. Is this $f$ itself, or a scalar multiple? Which harmonic polynomial we could expect to be canonically associated to any $f$?