Generating all rotationally invariant linear differential operators It is relatively easy to show that the Laplacian
$$
\Delta = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}
$$
Is the unique second order linear differential operator that is invariant under rotations in the sense that
$$
\Delta (f(R\mathbf{x})) = (\Delta f)(R\mathbf{x}).
$$
The way I remember proving this was to write down a general 
$$
D = \sum_i a_i \frac{\partial}{\partial x_i} + \sum_{ij} b_{ij} \frac{\partial^2}{\partial x_i \partial x_j},
$$
and then demand the invariance property. A multiple of the Laplacian will then fall out.
I am wondering if there is a way to general all such operators (up to a certain degree). Are they all powers of the Laplacian? What happens for the vector Laplacian?
 A: They are not all powers of the Laplacian.  For example, if a function $g$ is invariant under rotation, then the 0th order differential operator $D_g(f) = g f$ is invariant under rotation.  There is also the (first-order) radial derivative operator
$$
R(f) = xf_x + yf_y\,.
$$
Any (non-commutative) polynomial in $D_g$, $R$ and $\Delta$ will be a rotationally invariant differential operator, but these are not independent.  For example, we have the identity
$$
\left[\Delta, D_{x^2+y^2}\right] = D_4 + 4R,
$$
so $R$ is already in the ring generated by $\Delta$ and the $D_g$ where $g$ varies through rotationally invariant functions.  
Added comment:  Of course, I should have mentioned the other first order linear rotationally invariant operator, the angular derivative operator:
$$
A(f) = x\,f_y - y\,f_x\,,
$$
with the identity
$$
R^2 + A^2 = D_{x^2+y^2}\,\Delta
$$
A: Polynomials in the Laplacian are invariant under all isometries, not just rotations. This characterizes them by e.g. Helgason (1959, Thm 11) $=$ (1962, Prop. X.2.10) $=$ (1984, Prop. II.4.11), or Folland (1995, Thm 2.1).
A: More a comment than an answer - since $\nabla $ transforms as a vector, you could construct scalars by contracting an arbitrary number (say $n$) of $\nabla $s with an arbitrary rank-$n$ tensor. So it seems that to generate all operators in question, you'd at least have to generate all tensors on your underlying space (which implicitly seems to be 2-dimensional?). You'd discard some, such as antisymmetric of rank greater than 2, since the contraction would vanish, but it still seems there is a large number of possibilities. Perhaps some irreducibility specification is needed to make the question sharper?
