Assume that $V$ is a vector field on a Riemannian manifold $(M,g)$ with natural volume form $\Omega$ arising from $g$. Assume that the solution curves of $V$ are parametrized geodesics of the Riemannian metric $g$.
Is it true to say that the space of harmonic functions is invariant under the derivational operator $D(f)=V.f=df(V)$?
The question is related to the following post