Wikipedia mentions (here and here) that the Lie derivative has the following appealing commutator: $$[\mathcal{L}_X,\iota_Y]=\iota_{[X,Y]}$$ The only way I know to demonstrate this identity relies on coordinate-based manipulations, which are summarized (with a mistake) in this problem. I would like a coordinate-free (purely algebraic, one might say) proof of the same result.
Here's one possible route to a proof, but I can't seem to stick the landing. Let $\phi_t$ be the diffeomorphisms infinitesimally generated by $X$ and recall that $\mathcal{L}_X\omega=\left.\partial_t(\phi_t^*\omega)\right|_{t=0}$. Similarly, $[X,Y]=\mathcal{L}_XY=\left.\partial_t(\phi_{-t}^*Y)\right|_{t=0}$.
Now $\iota_{(\cdot_1)}(\cdot_2)$ is bilinear and smooth, so it commutes with the (Gateaux) time derivative when the equation is interpreted sufficiently weakly. So it suffices to show that $$\left.\partial_t\left([\phi_t^*,\iota_Y]-\iota_{\phi_{-t}^*Y}\right)\right|_{t=0}=0$$ But, applying LHS to a test form, I don't see how to deduce the result.
This answer suggests that it should be obvious from the Leibniz rule, but I don't see why the contraction operator $C$ in that answer should commute with $\mathcal{L}_X$ (as it does, moving from the first to second line).
How can I prove the commutator claim in a coordinate-independent manner?