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.