On equation $\Delta \circ \partial/\partial X=\partial/\partial X \circ \Delta$ on a Riemannian manifold Assume that $M$ is  a compact Riemannian manifold whose  Laplacian is  denoted  by $\Delta$.  Assume  that the  Euler  characteristic of  $M$ is  zero. Does  $M$ admit a  non vanishing  vector  field  $X$ which satisfy $$(*) \qquad \Delta \circ \partial/\partial X=\partial/\partial X \circ \Delta$$
What  can be said  about the  structure of  the  Lie  algebra of  all  vector  fields $X$ with the property $(*)$?
As  a  second  question: Every  vector  field  $X$ on $M$ defines a second order differential operator  on $C^{\infty}(M)$  with $$D(f)=\Delta(X.f)-X.\Delta(f)$$
This  is a  second order operator  since the third order terms  cancel  each others.
What is  the principal symbol of this operator , precisely? Can this PDE  be an elliptic operator when $M$ is a compact manifold?(I mean:is there  an  example  of this  situation in compact case?) Does every  compact  manifold admit a vector  field $X$ for which this  PDE would be  an elliptic operator?What would  be a  dynamical interpretation for the  index of this PDE. This  is  a  dynamical motivation for the later question..
 A: This is an expansion on Igor's comment and fixes my previous mistake (see also Willie's answer).  A direct computation gives
$$ D(f) = 2f^{ab}\nabla_{(a}X_{b)} + X_a(\nabla^b\nabla_b\nabla^a-\nabla^a\nabla_b\nabla^b)f + f^a\nabla^b\nabla_bX_a . $$
The middle summand is the same as $R_{ab}f^a$, while the formula $\nabla^b\nabla_aX_b=\nabla_a\nabla^bX_b+R_a{}^bX_b$ yields
$$ D(f) = 2f^{ab}\nabla_{(a}X_{b)} + (\nabla^b\nabla_bX_a + \nabla^b\nabla_aX_b - \nabla_a\nabla^bX_b)f^a . $$
This can be rewritten
$$ D(f) = \langle L_Xg, \nabla^2f \rangle + \langle \nabla f, \delta L_Xg - d\delta X\rangle , $$
which easily gives the formula for $D$ when $X$ is a conformal Killing field.  The interpretations for your questions are as follows:
Case 1: $\nabla_{(a}X_{b)}$ is not identically zero.  This is when $D$ is a second-order operator.  Note that the trace of $\nabla_{(a}X_{b)}$ is the divergence. By the divergence theorem, the integral of the trace is zero, so the bilinear form $\nabla_{(a}X_{b)}$ cannot be positive definite.  Thus there is no example of the type requested in your second question.  (Note that on noncompact manifolds there are examples: e.g. Euclidean space with the Euler vector field $X=\sum x^i\partial_{x^i}$, so $D$ is proportional to the Laplacian.)
Case 2: $\nabla_{(a}X_{b)}\equiv0$.  (That is, $X$ is Killing.)  Since the trace of $L_Xg$ is $2\delta X$, we see that $D\equiv0$.  That is, the Lie algebra of vector fields satisfying ($\ast$) is the Lie algebra of Killing vector fields.
A: The following formula is known among the experts but hard to find in the literature, so I figure I will document it here. Throughout $(M,g)$ denote an arbitrary pseudo-Riemannian manifold, and $\nabla$ its Levi-Civita connection. 

Definition Given a vector field $X$, its corresponding 0th order deformation tensor is defined to be ${}^{(X,0)}\pi := \mathcal{L}_X g$, where $\mathcal{L}_X$ is Lie differentiation with respect to $X$.
  The corresponding 1st order deformation tensor is defined using a formula similar to that of Christoffel symbols:
  $$ {}^{(X,1)}\pi_{ab}{}^c := \frac12 g^{cd} \left[ \nabla_a ( {}^{(X,0)}\pi_{bd}) + \nabla_b ({}^{(X,0)}\pi_{ad}) - \nabla_d ({}^{(X,0)}\pi_{ab}) \right] $$
Lemma Let $\Xi$ be an arbitrary $k$-covariant tensor field. And let $X$ be a vector field. The following formula holds for the commutation:
  $$ [ \nabla_a, \mathcal{L}_X ] \Xi_{b_1\cdots b_k} = \sum_{j = 1}^k {}^{(X,1)}\pi_{a b_j}{}^c \Xi_{b_1 \cdots b_{j-1} c b_{j+1} \cdots b_k} $$

With the aid of these formulas, we have immediately that, writing $\triangle_g$ for the Laplace-Beltrami operator, first
$$ [ \nabla_X, \triangle_g] f = [\mathcal{L}_X, \triangle_g ] f $$
because Lie derivation and covariant differentiation act identically on scalars, and then
$$ [\mathcal{L}_X, g^{ab}\nabla_a\nabla_b] f = \mathcal{L}_X (g^{ab}) \nabla^a \nabla_b f + g^{ab} [\mathcal{L}_X, \nabla_a] \nabla_b f + g^{ab} \nabla_a [\mathcal{L}_X, \nabla_b ]f $$
The first factor we can compute to get
$$ \mathcal{L}_X(g^{ab}) = - {}^{(X,0)}\pi^{ab} $$
using that $g^{ab} g_{bc} = \delta^a_c$. The third factor vanishes because Lie differentiation commutes with exterior differentiation. And we use our Lemma for the second term. We get, finally
$$ [\nabla_X, \triangle_g] f = - {}^{(X,0)}\pi^{ab} \nabla_a\nabla_b f - g^{ab} ~{}^{(X,1)}\pi_{ab}{}^c \nabla_c f. $$
Remarks: 


*

*Notice that the first order deformation tensor is defined in terms of the 0th order one. So that when $X$ is Killing, automatically both the ${}^{(X,0)}\pi$ and ${}^{(X,1)}\pi$ vanish, and differentiation with $X$ commutes with the Laplacian. 

*In the case ${}^{(X,0)}\pi = \phi g$ for some scalar function $\phi$ (so $X$ is conformally Killing), one can check that the formula reduces to the one I gave in a comment above. 

*When the function $\phi$ in the previous item is a non-zero constant (which some people refer to as $X$ being a homothetic vector field) one gets the special case
$$ [ \nabla_X, \triangle_g] f = \phi \triangle_g f $$
