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.