The bracket of two Killing fields is Killing

By the smooth bracket construction, $B=[X,Y]$ is a smooth vector field. The tensor Lie derivative applies to $g$ and to the smooth covariant two-tensors ${\mathcal{L}}_{X}g$ and ${\mathcal{L}}_{Y}g$; its flow definition agrees with the metric Lie derivative in the Killing condition. Thus both iterated derivatives and ${\mathcal{L}}_{B}g$ are defined.

Substitute the formula $({\mathcal{L}}_{V}S{)}_{ij}=V^a{\mathrm{\partial }}_a{S}_{ij}+S_{aj}{\mathrm{\partial }}_i{V}^a+S_{ia}{\mathrm{\partial }}_j{V}^a$ twice into ${\mathcal{L}}_X{\mathcal{L}}_{Y}g-{\mathcal{L}}_Y{\mathcal{L}}_{X}g$. The smooth metric and vector-field coordinate components are $C^2$ on each chart, so their mixed second partial derivatives commute. The mixed second metric derivatives cancel. The coefficients of ${\mathrm{\partial }}_a{g}_{ij}$ become $B^a=X^b{\mathrm{\partial }}_b{Y}^a-Y^b{\mathrm{\partial }}_b{X}^a$, and the coefficients multiplying $g_{aj}$ and $g_{ia}$ are ${\mathrm{\partial }}_i{B}^a$ and ${\mathrm{\partial }}_j{B}^a$, respectively. The remaining terms cancel in pairs. By the same formula applied to $B$, this is exactly ${\mathcal{L}}_{B}g$.

Both iterated Lie derivatives vanish because each inner derivative of $g$ vanishes and the coordinate formula sends the zero tensor to zero. The identity in step 2.1 therefore gives ${\mathcal{L}}_{[X,Y]}g=0$, so the smooth bracket field is Killing.