I'm trying to understand the equivalence Killing vector fields $\iff$ parallel sections on $\Lambda^2(M)\oplus TM$, for a Riemannian manifold $M$. I suppose that there exists a morphism $\phi:TM\rightarrow \Lambda^2(M)\oplus TM$ and a connection $D$ on $\Lambda^2(M)\oplus TM$ such that $\phi$ yields an isomorphism between the space of Killing vector fields and the space of parallel sections w.r.t $D$ (is it true? it can be that $\phi$ can only be defined for the Killing vector fields, not for all the vector fields...).
How do I construct $\phi$ and $D$?
I might have to use the Levi-Civita connection $\nabla$ on $M$. Applying $\nabla$ on $X\in TM$ repeatedly I obtain $\nabla X\in T^*M\otimes TM$ and $\nabla^2 X\in \Lambda^2(M)\otimes TM$. How can I use them to construct $\phi(X)$? Maybe I need to pass to local coordinates?
What you are looking for is referred to as the prolongation of the Killing equation or as rewriting it in first order closed form. Indeed, you can construct a first order differential operator $\phi$ and a connection $D$ in such a way that $X$ is a Killing field if and only if $D\phi(X)=0$. You can view the following as local coordinate computations but they actually make sense globally if you interpret them as abstract index notation.
The key is to understand consequences of the Killing equation. Converting a Vector field $X^i$ to a one form $X_i$, the Killing equation is that $\nabla_iX_j$ is skew symmetric in $i$ and $j$. But then it coincides with the exterior derivative of the one-form $X_j$, so can exploit the fact that $d^2=0$. Using once more that $d$ can be obtained as the alternation of the covariant derivative, you see that for a Killing field, you get $$ 0=\nabla_i\nabla_jX_k+\nabla_k\nabla_iX_j+\nabla_j\nabla_kX_i. $$ Using the Killing equation once more, the last summand can be rewritten as $-\nabla_j\nabla_iX_k$ and then by definition of the curvature, the sum with the first term becomes $R_{ij}{}^\ell{}_kX_\ell$. So what you conclude is that for a Killing field you get $$\nabla_k\nabla_iX_j=-R_{ij}{}^\ell{}_kX_\ell. \qquad (*) $$ Now using a vector notation for the sum $TM\oplus\Lambda^2(M)$, you can define $$ \phi(X)=\binom{X^i}{\nabla_{[j}X_{k]}} \qquad D_a\binom{X^b}{\mu_{cd}}=\binom{\nabla_aX^b-\mu_{a}{}^b}{\nabla_a\mu_{cd}+R_{cd}{}^e{}_aX_e}. $$ One immediately verifies that $D$ indeed is a connection.
Then from $(*)$, we immediately conclude that for a Killing field $X$, $D\phi(X)=0$. Conversely $D_a\binom{X^b}{\mu_{cd}}=0$ first implies that $\mu_{ab}=\nabla_aX_b$ and since $\mu$ is skew symmetric, $X$ is a Killing field. Inserting this into the second row, you see that this vanishes by $(*)$.