Solution (source code)

= Solution

Write $\mathcal L_XY=[X,Y]$. Expanding the definition on $fZ$ gives terms proportional to derivatives of $f$ with coefficient
$$
X(Yf)-[X,Y]f-Y(Xf)=0.
$$
All remaining terms carry an overall factor $f$, so
$$
(\mathcal L_X\nabla)_Y(fZ)=f(\mathcal L_X\nabla)_YZ.
$$
The same argument gives linearity in $Y$, confirming that the <Lie derivative of an affine connection> is tensorial.

Using the torsion-free identities $[X,Y]=\nabla_XY-\nabla_YX$ and $[X,Z]=\nabla_XZ-\nabla_ZX$, expand
$$
\begin{aligned}
(\mathcal L_X\nabla)_YZ
={}&[X,\nabla_YZ]-\nabla_{[X,Y]}Z-\nabla_Y[X,Z]\\
={}&R(X,Y)Z+\nabla_Y(\nabla_ZX)-\nabla_{\nabla_YZ}X.
\end{aligned}
$$
This is the required formula.

Solved by gpt-5.6-sol high.