show/hide this revision's text 3 added 15 characters in body

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 $h_i$ 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.
show/hide this revision's text 2 added 31 characters in body

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 a an manifold $M$ defined by the equations $F_1=0,\dots, F_1=0, \dots, F_k=0$ , $k<m$, with smooth $F$'s, F$'s and $k<\mathrm{dim}\ M$, then a smooth function $h$ vanishes on $N$ iff there exist smooth functions $h_i$ on $M$ such that $f=\sum_{i=1}^k h_i F_i$ $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.
show/hide this revision's text 1

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 a manifold $M$ defined by the equations $F_1=0,\dots, F_k=0$, $k<m$, with smooth $F$'s, then a function $h$ vanishes on $N$ iff there exist smooth functions $h_i$ on $M$ such that $f=\sum_{i=1}^k h_i F_i$ 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.