When do Harmonic polynomials constitute the kernel of a differential operator? Let $f$ be a real polynomial of two variables. Let $\partial_f=f\left(\frac{\partial}{\partial x},\frac{\partial}{\partial y}\right)$. Let $H$ denote the space of harmonic polynomials, i.e., polynomials $h$ such that $\partial_{x^2+y^2}h=0$.
Is it true that if the kernel of $\partial_f$ contains $H$ then $x^2+y^2|f(x,y)$? 
I would appreciate an elementary proof/disproof.
 A: It seems the answer is yes. We can write $F(x,y)=\sum_{d=0}^Nf_d(x,y)$ where $f_d$ is a homogeneous polynomial of degree $d$. Since for all $n$, $P_n(x,y):=(x+iy)^n+(x-iy)^n$ and $Q_n(x,y):=(x+iy)^n-(x-iy)^n$ are harmonic (because holomorphic), we should have $\partial_FP_n=\partial_FQ_n=0$. If $d>n$, we have $\partial_{f_d}P_n=\partial_{f_d}Q_n=0$. We have to show that for all $2\leq d\leq N$, $(x^2+y^2)\mid f_d(x,y)$ (we have $f_0=f_1=0$, evaluating $\partial_F$ on $1$, $x$ and $y$). We fix $d$ and we assume that $(x^2+y^2)\mid f_k$ for $2\leq k< d$. We have $\partial_{f_d}P_d=\partial_{f_d}Q_d=0$. We can write $f_d=\sum_{k=0}^da_{d,k}x^ky^{d-k}$. So 
\begin{align*}
\partial_{f_d}P_d&=\sum_{k=0}^da_{d,k}\partial_x^k\partial_y^{d-k}\left((x+iy)^d+(x-iy)^d\right)\\
&=\sum_{k=0}^da_{d,k}\frac{d!}{(d-k)!}\partial_y^{d-k}\left((x+iy)^{d-k}+(x-iy)^{d-k}\right)\\
&=\sum_{k=0}^da_{d,k}\frac{d!}{(d-k)!}\left((i^{d-k}+(-i)^{d-k}\right)(d-k)!\\
&=\sum_{k=0}^da_{d,k}d!\left((i^{d-k}+(-i)^{d-k}\right)
\end{align*}
hence 
$$\tag{1}\sum_{k=0}^da_{d,k}i^{d-k}(1+(-1)^{d-k})=0\mbox{ and }\sum_{k=0}^da_{d,k}i^{d-k}(1-(-1)^{d-k})=0.$$
Now, we have to consider the case $d$ odd and $d$ even. If $d=2p$, then we have relationships of the form 
$$a_{2p,2p}-a_{2p,2p-2}+a_{2p,2p-2}-\cdots= 0$$ and 
$$a_{2p,2p-1}-a_{2p,2p-3}+a_{2p,2p-5}-\cdots= 0.$$
Putting that into the expression of $f_{2p}$, we get after computations 
$$f_{2p}(x,y)=\sum_{k=1}^pa_{2p,2k}y^{2d-2k}(x^{2k}-(-1)^ky^{2k}) +\sum_{k=1}^{p-1}a_{2p,2k+1}y^{2d-2k-1}x(x^{2k}-(-1)^ky^{2k}),$$
and since $x^{2k}-(-1)^ky^{2k}=(x^2)^k-(-y^2)^k=(x^2+y^2)\sum_{l=0}^{k-1}x^{2l}(-y^2)^{l-k}$, $(x^2+y^2)\mid f_{2p}$.
The same computations give the result for $d$ odd, since if $d=2p+1$, (1) becomes $$\sum_{k=0}^{2p+1}a_{2p+1,k}i^{-k}(1-(-1)^k)=0\mbox{ and }\sum_{k=0}^{2p+1}a_{2p+1,k}i^{-k}(1+(-1)^k)=0.$$
