Solution (source code)

= Solution

Let $D_{abc}=\nabla_a\nabla_bK_c$. Differentiating the <Killing equation> gives $D_{abc}=-D_{acb}$. For a covector, the <Ricci identity> is
$$
D_{abc}-D_{bac}=-R^d{}_{cab}K_d=R_{cdab}K^d.
$$
Call this difference $C_{abc}$. The differentiated <Killing equation> then gives
$$
C_{abc}=D_{abc}+D_{bca},\quad
C_{bca}=D_{bca}+D_{cab},\quad
C_{cab}=D_{cab}+D_{abc}.
$$
Taking the first minus the second plus the third yields
$$
2D_{abc}=(R_{cdab}-R_{adbc}+R_{bdca})K^d.
$$
The pair symmetries of the <Riemann curvature tensor> and the <first Bianchi identity> reduce the coefficient to $2R_{cbad}$. Therefore the <second covariant derivative of a Killing vector> is
$$
\boxed{\nabla_a\nabla_bK_c=R_{cbad}K^d,\qquad
\nabla_a\nabla_bK^c=R^c{}_{bad}K^d.}
$$

For an <affine parameter> $\lambda$ on a <geodesic>, its tangent obeys $V^b\nabla_bV^a=0$. Thus
$$
\frac{d}{d\lambda}(K_aV^a)
=V^bV^a\nabla_bK_a+K_aV^b\nabla_bV^a=0.
$$
The second term vanishes by the <geodesic equation>; the first contracts a symmetric product $V^aV^b$ with the antisymmetric derivative from the <Killing equation>. This establishes the <geodesic conserved quantity from a Killing vector>:
$$
\boxed{K_aV^a\text{ is constant along every affinely parametrized geodesic}.}
$$