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


[Vector Fields in a Riemannian Manifold](https://mathoverflow.net/questions/215288/vector-fields-in-a-riemannian-manifold/215300#comment671908_215300)



Is the space of harmonic functions invariant under the derivational operator associated with a geodesible flow?
Take a [warped product](https://arxiv.org/pdf/1307.0236.pdf) $M=\mathbb S^1\times\_f\mathbb S^1$ for a nonconstant smooth function $f$.


The horizontal vector field $V$ satisfies your condition,
but for harmonic function is not invariant $f$, the function $Vf$ is not harmonic.


To see this, consider the $V$-flow $\Phi^t\colon M\to M$.
If $Vf$ is always harmonic then so is $f\circ\Phi^t$ for any $t$.
The latter fails evidently.


