Differential operators preserving the space of harmonic functions (aka higher symmetries of the Laplacian) The article http://arxiv.org/abs/hep-th/0206233 (published in Ann. of Math. (2)  161  (2005),  no. 3) deals with linear differential operators $D$ for which there exists another linear differential operator $\delta$ such that $\Delta D = \delta \Delta$. Obviously these operators preserve the kernel of $\Delta$, i.e. the space of harmonic functions. The mentioned article finds essentially all such operators $D$. The result is that up to trivial operators $D = P\Delta$ all the operators $D$ have polynomial coefficients and are generated by sums of compositions of first order operators of this kind.
First question: Let $D$ be any differential operator preserving the space of harmonic functions. It is easy to see that the operator $\delta = \Delta D (\Delta)^{-1}$ is well defined and satisfies $\Delta D = \delta \Delta$. Is $\delta$ also a differential operator?
Second question: Is it true that all differential operators, which preserve the space of harmonic functions, are generated by first order ones with this property?
One can also ask these questions only for linear differential operators or for operators from the Weyl algebra (i.e. linear differential operators with polynomial coefficients). For example, by a theorem of Peetre, the answer to the first question is affirmative if the operator $\delta = \Delta D (\Delta)^{-1}$ is local (i.e. the support of $\delta u$ is contained in the support of $u$).
Third question: What makes the linked article so interesting that it was published in Annals?
 A: The answer to your second question (unless I somehow misread it) is yes precisely because of the result of the paper you refer to (you may also wish to look at this paper and the preprint math-ph/0506002 which address the same subject). This is the case because if $D$ is a differential operator that preserves the space of harmonic functions then there indeed exists a differential operator $\delta$ such that $\Delta D = \delta \Delta$. The latter holds (see e.g. the discussion at p.290 near Eq.(5.5) of the book Applications of Lie groups to Differential Equations by P.J. Olver) because the equation $\Delta f=0$ is totally nondegenerate in the sense of Definition 2.83 of the same book. In spite of the rather technical language the idea behind all this is very simple: if you have a submanifold $N$ of an manifold $M$ defined by the equations $F_1=0, \dots, F_k=0$ with smooth $F$'s and $k<\mathrm{dim}\ M$, then a smooth function $h$ vanishes on $N$ iff there exist smooth functions hj on $M$ such that 
$$h=h_{1} F_1+\cdots+h_k F_k$$ provided $dF_1\wedge \dots \wedge dF_k\neq 0$ on $N$ (see Proposition 2.10 of the same book). In a sense, this is a smooth counterpart of the famous Hilbert's Nullstellensatz in the form stated e.g. here. This result is then applied to the case when $M$ is a jet bundle and $N$ is a submanifold thereof defined by a system of differential equations and all its differential consequences (more precisely, one should rather consider the consequences only up to a certain order, to avoid dealing with infinitely many equations), et voila.
