Understanding proof of Cartan's magic formula: $L_X = i_X \circ d+d \circ i_X$

A possible proof of Cartan's magic formula $$L_X = i_X \circ d+d \circ i_X$$ is to follow the steps:

• Show that two derivations on $\Omega^{\bullet}(M)$ commuting with $d$ are equal iff they agree on $\Omega^0(M)$.
• Show that $L_X$ is a derivation on $\Omega^{\bullet}(M)$ commuting with $d$.
• Show that $i_X \circ d + d \circ i_X$ is a derivation on $\Omega^{\bullet}(M)$ commuting with $d$.
• Show that $L_X f = Xf = i_Xdf+ d i_Xf$ for all $f \in C^{\infty}(M)=\Omega^0(M)$.

I followed the sketch of proof, and the only point where I am stuck is when I have to prove that $L_X$ and $d$ commute.

More precisely, if I write $$L_Xd\omega =\frac{d}{dt}_{|t=0} \phi_t^* d\omega = \frac{d}{dt}_{|t=0} d\phi_t^* \omega = \lim\limits_{t \to 0} d \left( \frac{1}{t} \left( \phi_t^* \omega- \omega \right) \right),$$ is there an argument to permute $d$ and $\lim\limits_{t \to 0}$?

Putting in the limit made it more mysterious. You're just using the fact that for smooth functions mixed partials are equal. That is, $$\frac{\partial}{\partial t}\frac{\partial}{\partial x^j}=\frac{\partial}{\partial x^j}\frac{\partial}{\partial t}.$$

EDIT: In particular, if you write it out in local coordinates, it's computed by taking partial derivatives with the respect to the variables $x^1,\dots,x^n$. Writing $\phi_t^*\omega=\sum f_I(x,t)dx^I$ (where $I=(i_1,\dots,i_k)$ ranges over length-$k$ multiindices), then $d\phi_t^*\omega=\sum \frac{\partial f_I(x,t)}{\partial x^j}\,dx^j\wedge dx^I$, etc.

Here is an attempt to answer my question following Ted Shifrin's comment:

First, we may write $\phi_t^* \omega(x)= \sum\limits_{I} f_I(x,t)dx^I$ where $I$ ranges over length-$k$ multiindices. Then

$$\begin{array}{ll} dL_X \omega & = d \frac{d}{d t}_{|t=0} \phi_t^* \omega = d \sum\limits_I \frac{\partial }{\partial t}_{|t=0} f_I dx^I \\ & = \sum\limits_I \sum\limits_{i=1}^n \frac{\partial}{\partial x_i} \frac{\partial}{\partial t}_{|t=0} f_I dx^i \wedge dx^I = \sum\limits_I \sum\limits_{i=1}^n \frac{\partial}{\partial t}_{|t=0} \frac{\partial}{\partial x_i} f_I dx^i \wedge dx^I \\ & = \sum\limits_I \frac{d}{d t}_{|t=0} df_I \wedge dx^I = \frac{d}{dt}_{|t=0} d \phi_t^* \omega = \frac{d}{dt}_{|t=0} \phi_t^* d \omega = L_Xd\omega \end{array}$$

The external differential does not depend on the choice of chart. Therefore $$d\phi_t^*\omega=\phi_t^*d \omega.$$ Taking derivative we get $$dL_X\omega=L_Xd\omega.$$