"Noncommutative heat equation" -- a strange generalization of Killing vectors for a flat metric Let $(M,g)$ be a smooth (pseudo)Riemannian manifold with a flat metric $g$, and $X$, $Y$ be vector fields on $M$ such that
$$
L_X^2 (g)=L_Y(g). \hspace{70mm} \mbox{(1)}
$$ 
where $L_Z$ is the Lie derivative along $Z$ and $L_X^2(g)\equiv L_X(L_X(g))$.
If $X=0$ then $Y$ is just a Killing vector for $g$ but 

is there any (geometric?) interpretation for $X$ and $Y$ in the general case?

In fact I wonder whether Eq.(1), be it for fixed $g$ and unknown $X,Y$ or conversely for fixed $X,Y$ and unknown $g$, was studied systematically at all: I failed to find any relevant references.
Motivation: Equation (1) follows from the last formula in Proposition 5 from a paper on classification of compatible Hamiltonian structures that I have recently come across. It looks a bit like the heat equation for the metric, hence the title.
 A: To understand the local geometry of this equation, I think one should apply the Calabi resolution of the Killing equation.  (See E. Calabi, On compact, Riemannian manifolds with constant curvature. I, in Differential Geometry, C. B. Allendoerfer, ed., vol. 3 of Proceedings of Symposia in Pure Mathematics, pp. 155–180. AMS, Providence, RI, 1961.  For a useful modern exposition, see http://arxiv.org/pdf/1409.7212.pdf.)
The basic point is that, when $g$ is flat (or, more generally, has constant sectional curvature) the Killing operator $X\mapsto D_1(X) =  L_X(g)$ as a map from $C^\infty(TM)$ to $C^\infty(S^2(T^*M))$ can be embedded into a locally exact sequence of sheaves
$$
0 \longrightarrow \mathscr{K}_g\longrightarrow C^\infty(TM)\longrightarrow C^\infty(S^2(T^*M))\longrightarrow C^\infty\bigl(C^{(2,2)}(TM)\bigr)\longrightarrow\cdots\tag1
$$
that resolves the sheaf $\mathscr{K}_g$ of $g$-Killing vector fields on $M$ and where the sheaf mappings are given by known linear differential operators.  An element $h$ of $C^\infty(S^2(T^*M))$ can be written (locally) in the form $h = L_Y(g)$ for some vector field $Y$ if and only if $D_2(h) = 0$.  
This shows that if $L_X(L_X(g)) = L_Y(g)$, then $X$ must satisfy the system of differential equations
$$
D_2\bigl(L_X(L_X(g))\bigr) = 0.\tag2
$$
Conversely, if a vector field $X$ satisfies this equation, then, at least locally (and globally if $M$ is simply-connected), there will exist a vector field $Y$ such that $L_X(L_X(g)) = L_Y(g)$, and this $Y$ will be unique up to the addition of a Killing vector field.
The differential equation (2) for $X$ is generally overdetermined (and nonlinear) when the dimension $n$ of $M$ is greater than $2$, and one expects the space of its local solutions to be finite dimensional.  Thus, one expects that, even locally, the space of solutions of the original system for $(X,Y)$ will be finite dimensional when $n>2$.  To check this, one should compute the characteristic variety of (2).  If there is interest, I can do this and report the result.
A: Here is an observation which does not answer the question, but does at least tell where not to look for examples. It is shown that for compact Riemannian manifolds with non-positive Ricci curvature there are no interesting solutions to the stated equations, so to find interesting solutions one must look: a. in positive curvature, b. at indefinite metrics, c. on noncompact manifolds, or d. in infinite dimensions.
Suppose $M$ is compact and $g$ is Riemannian. It is claimed that if $g$ has non-positive Ricci curvature, then either the Ricci curvature of $g$ is somewhere negative and both $X$ and $Y$ are $0$, or $g$ is Ricci flat and both $X$ and $Y$ are parallel. To begin with, do not suppose anything about the Ricci curvature of $g$. Let $D$ be its Levi-Civita connection and raise and lower indices using $g_{ij}$ and the inverse symmetric bivector $g^{ij}$. In what follows I use the abstract index notations, so square brackets (resp. parentheses) denote anti-symmetrization (resp. symmetrization) over the enclosed indices. 
For any metric $(L_{X}g)_{ij} = 2D_{(i}X_{j)}$ and 
$$\begin{equation}
(L_{X}^{2}g)_{ij} =  X^{p}D_{p}(L_{X}g)_{ij} + (D_{i}X^{p})(L_{X}g)_{pj} + (D_{j}X^{p})(L_{X}g)_{ip}. 
\end{equation}$$
Trace this to obtain
$$\begin{equation}
g^{ij}(L_{X}^{2}g)_{ij} = 2X^{p}D_{p}D^{q}X_{q} + 4D^{(p}X^{q)}D_{(p}X_{q)}. 
\end{equation}$$
By assumption this equals $g^{ij}(L_{Y}g)_{ij} = 2D^{p}Y_{p}$. Since by assumption $M$ is compact and without boundary, integration by parts yields
$$\begin{align}
0 = 4\int_{M}D^{(p}X^{q)}D_{(p}X_{q)} - 2\int_{M}(D_{p}X^{p})^{2}.
\end{align}$$
In general, this yields no obvious conclusions. Go back to the equation preceeding the integration and commute derivatives to obtain 
$$\begin{align}
g^{ij}(L_{X}^{2}g)_{ij} = 2X^{p}D^{q}D_{p}X_{q} + 4D^{(p}X^{q)}D_{(p}X_{q)} -2R_{pq}X^{p}X^{q},
\end{align}$$
in which $R_{ij}$ is the Ricci curvature of $g_{ij}$. Now integrating by parts yields
$$\begin{align}
0 &= \int_{M}\left(4D^{(p}X^{q)}D_{(p}X_{q)} - 2D^{q}X^{p}D_{p}X_{q} - 2R_{pq}X^{p}X^{q}\right)\\
 & = 2\int_{M}\left(D^{(p}X^{q)}D_{(p}X_{q)} + D^{[p}X^{q]}D_{[p}X_{q]} -R_{pq}X^{p}X^{q}\right).
\end{align}$$
Because $g$ is Riemannian, if the Ricci curvature is non-positive this implies that $D_{(i}X_{j)} = 0$ and $D_{[i}X_{j]} = 0$ so that $D_{i}X_{j} = 0$ and $X$ is parallel. In the original equation this implies $Y$ is Killing. By the original Bochner argument, if the Ricci curvature is non-positive and somewhere negative then there is no non-zero Killing field, so the only possibility for non-trivial solutions is that $g$ be Ricci flat, in which case the Bochner argument forces $Y$ to be parallel. 
