Invariance of the space of harmonic functions under derivation associated to a non-vanishing vector field Let  $X$  be  a  non-vanishing real analytic vector  field  on an open  manifold $M$. What  kind  of  obstructions would  appear when we search for  a Riemannian metric  on $M$ such that the space  of  harmonic  functions would  be  invariant under the derivation operator $f \mapsto X.f$? 
A harmonic function is a function $f$ which satisfy $\Delta_g (f)=0$ where $\Delta_g$ is the Laplace operator associated to the metric $g$.
 A: (A partial answer; the full picture is likely to be very complicated, with special behaviors in dimensions 1 and 2 compared to higher ones.)
Dimension 1
Given a Riemannian metric on a one dimensional manifold there is, up to a sign, a unique "unit-length" vector field $Y$ tangent to $M$, and in this case you must have $\Delta_g = Y^2$. It is easy to see that $X$ preserves harmonic functions if and only if $[X,Y] = 0$, which implies $X = cY$ for some constant $c$. 
Therefore the only obstruction to to the existence of a suitable metric is whether $X$ vanishes: if $X$ is a non-vanishing vector field, then you can define a co-metric by $X\otimes X$. If $X$ vanishes at some point, then such a metric doesn't exist. 
General Dimensions
Given a metric $g$, it is well-known that 
$$ [X, \Delta_g]f = - {}^{(0)}\pi^{ab} \nabla^2_{ab} f - g^{ab} \cdot {}^{(1)}\pi_{ab}{}^c \nabla_c f \tag{*}$$
where
$$ -{}^{(0)}\pi^{ab} = \mathcal{L}_X g^{ab} = - \nabla^a X^b - \nabla^b X^a \tag{1}$$
and
$${}^{(1)}\pi_{ab}{}^c = \frac12 g^{cd} (\nabla_a {}^{(0)}\pi_{bd} + \nabla_b {}^{(0)}\pi_{ad} - \nabla_d {}^{(0)} \pi_{ab} ) \tag{2}
$$ 
For $X$ to preserve the space of harmonic functions, a sufficient condition (necessary and sufficient if you work only locally) can be found by examining the equations above. 
First note that it is not necessary that $[X,\Delta_g] = 0$. If we are only interested in preserving the harmonic functions, it suffices that $[X, \Delta_g] \propto \Delta_g$. Examining (*) this means that:


*

*We want ${}^{(0)}\pi^{ab} = k g^{ab}$ for some function $k$. This means that $X$ is a conformal isometry of $g$.

*We want ${}^{(1)}\pi_{ab}{}^c g^{ab}$ to vanish. Together with the previous statement we want
$$ 0 = (2-n) \nabla^c k,$$
where $n$ is the number of spatial dimensions. So (as is well known), we conclude that in 2 dimensions a sufficient condition is that $X$ is a conformal isometry of $g$, and in higher dimensions we require that $X$ is a homothety. 


Now, homotheties also obey the Jacobi equation, in the form 
$$ \nabla_a \nabla_b X_c+ R_{bcad} X^d = 0 $$
(for one choice of sign of the Riemann curvature tensor). A particular consequence is that if $X$ and $\nabla X$ both vanish at a point $p$, then $X$ vanish everywhere. This means that we have an obstruction:
Obstruction: Let $n \geq 3$. If there exists a point $p$ such that the vector field $X$ vanishes to second order (in the sense that for any function $f$ and any other vector field $Y$, both $Xf|_p = 0 = YXf|_p$), then there does not exist any metric $g$ on $M$ such that $X$ preserves the harmonic functions. 
Remark: In $n = 2$, this is not an obstruction. Consider the vector field $(y^2 - x^2) \partial_x - 2 xy \partial_y$ on $\mathbb{R}^2$. It is simple to check that this vector field preserves the harmonic functions on $\mathbb{R}^2$, but this vector field vanishes to second order at the origin. 
