Existence of a special vector field on Riemannian manifolds? In a Riemannian Manifold $(M,g)$ a vector field $X$ is said to be Killing vector field if $L_X g$=0 and is said to be conformal if $L_X g= fg$ for some smooth real function $f$ on $M$.
Also, the notion of bi-Killing vector field is defined (in which we have $L_X (L_X g)=0$. So it is natural to define bi-conformal vector field as $L_X (L_X g)=fg$. The second Lie derivative of such vector field is proportional to the metric tensor.
But, I need to find a vector field $X\in \mathcal{X}M $ in a Riemannian manifold $(M,g)$ such that the second Lie derivative be proportional to $X$, in fact

Is there exist a vector field $X$ such that $L_X (L_X g)=\lambda X^\flat \otimes X^\flat$?

Where, $\lambda$ is a real non-zero constant.
Any suggestion is highly appreciated.
 A: In dimensions greater than $1$, this is an overdetermined equation for $X$, if $g$ is assumed given, so one doesn't expect most metrics $g$ to admit such a vector field (other than the zero vector field).
In dimension $1$, one can assume that $g = {\mathrm{d}x}^2$.  If $X^\flat = f(x)\,\mathrm{d}x$, then $L_X(L_Xg) = \lambda\,X^\flat\otimes X^\flat$ becomes the equation
$$
2f(x)\,f''(x) + 4\,f'(x)^2 - \lambda\,f(x)^2=0.
$$
Multiplying by $\tfrac32 f(x)$, this equation becomes
$$
\bigl(f(x)^3\bigr)'' - \tfrac32\lambda\,f(x)^3 = 0, 
$$
so, if $\lambda = \tfrac23a^2>0$, then $f(x)^3$ is a linear combination of $\cosh ax$ and $\sinh ax$ while, if $\lambda = -\tfrac23a^2<0$, $f(x)^3$ is a linear combination of $\cos ax$ and $\sin ax$.
Meanwhile, if one wants to know the local generality (up to diffeomorphism) of the set of pairs $(g,X)$ satisfying $L_X(L_Xg) = \lambda\,X^\flat\otimes X^\flat$, where $X$ is nonvanishing, one can choose 'flow box' coordinates, in which $X = \partial/\partial x^1$.  The equation for the unknown $g = g_{ij}(x)\,\mathrm{d}x^i\mathrm{d}x^j$ then becomes a second-order, determined system that has unique local solutions when one specifies
$$
a_{ij}(x^2,\ldots,x^n)=g_{ij}(0,x^2,\ldots,x^n)\quad\text{and}\quad
b_{ij}(x^2,\ldots,x^n)=\frac{\partial g_{ij}}{\partial x^1}(0,x^2,\ldots,x^n).
$$
Thus, the local solutions appear to depend on $n(n{+}1)$ arbitrary functions of $n{-}1$ variables. However, one must bear in mind that flow-box coordinates for $X$ depend on $n$ functions of $n{-}1$ variables. Taking this 'flexibility' in mind, one finds that the local solution pairs $(g,X)$ up to diffeomorphism depend on $n^2$ functions of $n{-}1$ variables.  Since the general metric in $n$ dimensions up to diffeomorphism depends on $\tfrac12n(n{-}1)$ functions of $n$ variables, it follows that, when $n>1$, for the generic metric $g$, the only $X$ for which $(g,X)$ satisfies the given equation is $X\equiv0$.
For the case $n=2$, it might be interesting to determine the conditions that the geometric invariants of a surface metric must satisfy in order to support a nonzero $X$ that satisfies $L_X(L_Xg) = \lambda\,X^\flat\otimes X^\flat$.
