It is a matter of calculation. Take successive derivatives of the conformal Killing equation, until you notice that you can solve for all derivatives of $X$ at a certain order (3rd order here) in terms of lower order ones. I don't know if there is a canonical reference, but a calculation equivalent to such an equation on $X$ is summarized in Eqs. (69.2--5) of

Eisenhart, L. P., Riemannian geometry, Princeton: Princeton University Press. vii, 306 p. (1949). ZBL0041.29403. Probably (69.2--5) was already there in the original (1926) edition.

More explicitly, from Eisenhart: \begin{gather} \tag{69.1}\label{69.1} \nabla_j \xi_i + \nabla_i \xi_j = \psi g_{ij} \\ \tag{69.2}\label{69.2} \nabla_k \nabla_j \xi_i = -\xi_m R^m{}_{kij} + \frac{1}{2} (g_{ij} \nabla_k \psi + g_{ik} \nabla_j \psi - g_{jk} \nabla_i \psi) \\ \tag{69.4}\label{69.4} g^{il} \xi_m \nabla_l R^m{}_{kij} - \xi_m \nabla_k R^m{}_j - \nabla_k \xi_m R^m{}_j - \nabla_j \xi_m R^m{}_k + \frac{(n-2)}{2} \nabla_k \nabla_j \psi + \frac{1}{2} g_{jk} \nabla^l \nabla_l \psi = 0 \\ \tag{69.5}\label{69.5} \nabla^l \nabla_l \psi = \frac{2}{n-1} (\xi_m \nabla_i R^{mi} + \nabla_i \xi_m R^{mi}) \end{gather} For completeness, note that basically $\psi = \frac{2}{n} \nabla^i \xi_i$. To get a 3rd order equation for $\xi_i$, start by using \eqref{69.5} to eliminate $\nabla^l \nabla_l \psi$ from \eqref{69.4}, and use the result to solve for $\nabla_j \nabla_i \psi$. Then, take an extra derivative of \eqref{69.2} and eliminate $\nabla_j \nabla_i \psi$.

Feel free to complete the calculation and add the resulting equation to this answer.

Probably, this result has been rederived independently multiple times in the literature, or referenced in passing as a "well-known" result without citing any convenient source, as in the article referenced by the OP.

More generally, any PDE on an unknown $X$ which has a differential consequence of the form $\nabla^k X = F(X, \ldots, \nabla^{k-1} X)$ is known as a PDE of finite type. While many geometric equations (like variants of the Killing and conformal Killing equations) are known by folklore to be of finite type, it can be surprisingly non-trivial to locate an original or convenient reference for such facts.

It is a matter of calculation. Take successive derivatives of the conformal Killing equation, until you notice that you can solve for all derivatives of $X$ at a certain order (3rd order here) in terms of lower order ones. I don't know if there is a canonical reference, but a calculation equivalent to such an equation on $X$ is summarized in Eqs. (69.2--5) of

Eisenhart, L. P., Riemannian geometry, Princeton: Princeton University Press. vii, 306 p. (1949). ZBL0041.29403. Probably (69.2--5) was already there in the original (1926) edition.

More explicitly, from Eisenhart: \begin{gather} \tag{69.1}\label{69.1} \nabla_j \xi_i + \nabla_i \xi_j = \psi g_{ij} \\ \tag{69.2}\label{69.2} \nabla_k \nabla_j \xi_i = -\xi_m R^m{}_{kij} + \frac{1}{2} (g_{ij} \nabla_k \psi + g_{ik} \nabla_j \psi - g_{jk} \nabla_i \psi) \\ \tag{69.4}\label{69.4} g^{il} \xi_m \nabla_l R^m{}_{kij} - \xi_m \nabla_k R^m{}_j - \nabla_k \xi_m R^m{}_j - \nabla_j \xi_m R^m{}_k + \frac{(n-2)}{2} \nabla_k \nabla_j \psi + \frac{1}{2} g_{jk} \nabla^l \nabla_l \psi = 0 \\ \tag{69.5}\label{69.5} \nabla^l \nabla_l \psi = \frac{2}{n-1} (\xi_m \nabla_i R^{mi} + \nabla_i \xi_m R^{mi}) \end{gather} For completeness, note that basically $\psi = \frac{2}{n} \nabla^i \xi_i$. To get a 3rd order equation for $\xi_i$, start by using \eqref{69.5} to eliminate $\nabla^l \nabla_l \psi$ from \eqref{69.4}, and use the result to solve for $\nabla_j \nabla_i \psi$. Then, take an extra derivative of \eqref{69.2} and eliminate $\nabla_j \nabla_i \psi$.

Feel free to complete the calculation and add the resulting equation to this answer.

A quick sketch of the required steps to prove (69.2-5) are as follows:

• (69.1) is the conformal Killing equation
• To obtain (69.2), observe that the first Bianchi identity implies $$\nabla_{[i}\nabla_j \xi_{k]} = 0 $$ for any vector field $\xi$; here brackets denote full antisymmetrization. Write $\pi_{ij} = \nabla_i \xi_j + \nabla_j \xi_i$ for the symmetric part, you can write it as $$ \nabla_i (\nabla_j \xi_k - \frac12 \pi_{jk}) + \nabla_j(\nabla_k \xi_i - \frac12 \pi_{ki}) - \nabla_k(\nabla_j \xi_i - \frac12 \pi_{ji}) = 0 $$ move all the $\pi$ terms to the RHS and the second and third $\nabla^2\xi$ terms provide the Riemann tensor.
• (69.4) is trickier. Start by taking the divergence $\nabla^i$ of (69.2). The LHS looks like $\nabla^i\nabla_k\nabla_j \xi_i = \nabla_k \nabla_j \nabla^i \xi_i + [\nabla^i, \nabla_k\nabla_j] \xi_i$. The commutator generates terms of the form $\nabla \mathrm{Ric} \cdot \xi$ and $\mathrm{Ric} \cdot \nabla X$. Apply also $\nabla^i \xi_i = \frac{n}{2} \psi$.
• (69.5) is simply the trace of (69.4).

Probably, this result has been rederived independently multiple times in the literature, or referenced in passing as a "well-known" result without citing any convenient source, as in the article referenced by the OP.

More generally, any PDE on an unknown $X$ which has a differential consequence of the form $\nabla^k X = F(X, \ldots, \nabla^{k-1} X)$ is known as a PDE of finite type. While many geometric equations (like variants of the Killing and conformal Killing equations) are known by folklore to be of finite type, it can be surprisingly non-trivial to locate an original or convenient reference for such facts.

It is a matter of calculation. Take successive derivatives of the conformal Killing equation, until you notice that you can solve for all derivatives of $X$ at a certain order (3rd order here) in terms of lower order ones. I don't know if there is a canonical reference, but a calculation equivalent to such an equation on $X$ is summarized in Eqs. (69.2--5) of

Eisenhart, L. P., Riemannian geometry, Princeton: Princeton University Press. vii, 306 p. (1949). ZBL0041.29403. Probably (69.2--5) was already there in the original (1926) edition.

Probably, this result has been rederived independently multiple times in the literature, or referenced in passing as a "well-known" result without citing any convenient source, as in the article referenced by the OP.

More generally, any PDE on an unknown $X$ which has a differential consequence of the form $\nabla^k X = F(X, \ldots, \nabla^{k-1} X)$ is known as a PDE of finite type. While many geometric equations (like variants of the Killing and conformal Killing equations) are known by folklore to be of finite type, it can be surprisingly non-trivial to locate an original or convenient reference for such facts.

It is a matter of calculation. Take successive derivatives of the conformal Killing equation, until you notice that you can solve for all derivatives of $X$ at a certain order (3rd order here) in terms of lower order ones. I don't know if there is a canonical reference, but a calculation equivalent to such an equation on $X$ is summarized in Eqs. (69.2--5) of

Eisenhart, L. P., Riemannian geometry, Princeton: Princeton University Press. vii, 306 p. (1949). ZBL0041.29403. Probably (69.2--5) was already there in the original (1926) edition.

More explicitly, from Eisenhart: \begin{gather} \tag{69.1}\label{69.1} \nabla_j \xi_i + \nabla_i \xi_j = \psi g_{ij} \\ \tag{69.2}\label{69.2} \nabla_k \nabla_j \xi_i = -\xi_m R^m{}_{kij} + \frac{1}{2} (g_{ij} \nabla_k \psi + g_{ik} \nabla_j \psi - g_{jk} \nabla_i \psi) \\ \tag{69.4}\label{69.4} g^{il} \xi_m \nabla_l R^m{}_{kij} - \xi_m \nabla_k R^m{}_j - \nabla_k \xi_m R^m{}_j - \nabla_j \xi_m R^m{}_k + \frac{(n-2)}{2} \nabla_k \nabla_j \psi + \frac{1}{2} g_{jk} \nabla^l \nabla_l \psi = 0 \\ \tag{69.5}\label{69.5} \nabla^l \nabla_l \psi = \frac{2}{n-1} (\xi_m \nabla_i R^{mi} + \nabla_i \xi_m R^{mi}) \end{gather} For completeness, note that basically $\psi = \frac{2}{n} \nabla^i \xi_i$. To get a 3rd order equation for $\xi_i$, start by using \eqref{69.5} to eliminate $\nabla^l \nabla_l \psi$ from \eqref{69.4}, and use the result to solve for $\nabla_j \nabla_i \psi$. Then, take an extra derivative of \eqref{69.2} and eliminate $\nabla_j \nabla_i \psi$.

Feel free to complete the calculation and add the resulting equation to this answer.

Probably, this result has been rederived independently multiple times in the literature, or referenced in passing as a "well-known" result without citing any convenient source, as in the article referenced by the OP.

More generally, any PDE on an unknown $X$ which has a differential consequence of the form $\nabla^k X = F(X, \ldots, \nabla^{k-1} X)$ is known as a PDE of finite type. While many geometric equations (like variants of the Killing and conformal Killing equations) are known by folklore to be of finite type, it can be surprisingly non-trivial to locate an original or convenient reference for such facts.