Laplace-Beltrami and the isometry group H$\vphantom{a}$i. Consider the Laplacian on $\mathbb R^n$,
$$
\Delta=\partial_i^2
$$
It is easy to prove that the most general differential operator that commutes with rotations and translations is of the form
$$
\sum_k a_k\Delta^k
$$
for some constants $a_k$.
Is there a similar result for non-trivial manifolds?
In other words, consider the Laplacian on a manifold $\mathcal M$,
$$
\Delta=\nabla_i\nabla^i
$$
Under what conditions is $\Delta$ the most general differential operator that commutes with all Killing vectors of $\mathcal M$? It is clear that $[\Delta,\mathcal L_\xi]=0$ for all $\xi\in\mathrm{Iso}(\mathcal M)$, but can we rule out the existence of any other quadratic differential operator that commutes with $\xi$?
My gut tells me that this should be true at least for maximally symmetric spaces. Is this correct? A good reference will be very appreciated.
 A: For the maximally symmetric case with the second order operator, the proof is very simple. 
A second order differential operator acting on scalars can be written in the form 
$$ b^{ij} \nabla^2_{ij} + c^i \nabla_i + e $$
Consider first the action of all isometries that fixes a point $p\in M$. You must have that $b^{ij}(p)$ is a symmetric two tensor that is fixed by the action of the entire special orthogonal group, and so must be proportional to the inverse metric. Now as the principal part is the Laplacian and is invariant under the isometries, you can then conclude similarly that $c^i$ is a vector that is invariant under the action of the entire special orthogonal group, and hence vanishes. 
This implies that your differential operator is 
$$ \lambda \triangle_g + e $$
for scalar functions $\lambda$ and $e$. But then translation symmetry tells you immediately that both $\lambda$ and $e$ are constant. 

The maximal symmetry is used crucially in the argument: for example, consider the standard torus $\mathbb{T}^2$; it is homogeneous but does not admit any rotation symmetry. In standard coordinates any second order differential operator of the form $a_1 \partial^2_{11} + b \partial^2_{12} + a_2 \partial^2_{22}$ commute with all the Killing vector fields (being spanned by $\partial_1$ and $\partial_2$). 

Notice that in addition to the "natural" differential operators mentioned in the  link in Phillip Andreae's comments, there are also other "less natural" differential operators that can commute with the Killing vector fields. 
For example, let $(M,g)$ be any Riemannian manifold with only one Killing vector field $X$ (isometry group being $\mathbb{R}$ or $U(1)$; think a generic surface of revolution). Then quite clearly the second order operator $\mathcal{L}_X^2: f \mapsto X(X(f))$ commutes with the (one and only) Killing vector field, yet is not the Laplacian. 
