First of all, I think the formula for the homotopy should be
$$ g^1 \gamma -\gamma = \partial(H\gamma) - H(\partial \gamma).$$
As the left hand side is the top minus the bottom, the right hand side then should be the boundary of $I\times \sigma$ minus the cylindrical part of the boundary ($I\times \partial \gamma$).
Using that, we proceed to the proof. By fundamental theorem of calculus and $g^{s+t} = g^s g^t$, for any $k$-form $\omega$ we have
\begin{align}
(g^1)^*\omega - \omega &= \int_0^1\frac{d}{dt} (g^t)^*\omega|_{t=s} \ ds \\
&=\int_0^1 (g^s)^* \frac{d}{dt} (g^t)^* \omega \bigg|_{t=0} ds \\
&= \int_0^1 (g^s)^* \mathcal L_X \omega \ ds.
\end{align}
If we integrate this on $\gamma$, by Fubini's theorem,
$$\int_\gamma \int_0^1 (g^s)^* \mathcal L_X \omega \ ds = (-1)^k\int_0^1 \int_\gamma H_s^* \mathcal L_X \omega \ ds,$$
where $H_s : \sigma \to M$ is given by $H_s(x) = H(s, x) = g^s \gamma(x)$. Note that $H_0 = \gamma$.
On the other hand,
$$\begin{split}
\int_\gamma (g^1)^* \omega - \omega &= \int_{g^1 \gamma - \gamma} \omega \\
&= \int_{\partial (H\gamma) - H(\partial \gamma)} \omega \\
&= \int_{H\gamma} d\omega - \int_{H(\partial \gamma)} \omega.
\end{split}$$
Now we use that (prove later) for any $r$-form $\alpha$ and for $r-1$ chain $\gamma$ we have
$$\tag{1} \int_{H\gamma}\alpha = (-1)^{r-1}\int_0^1 \int_\gamma H_s^* i_X \alpha\ ds,$$
then
$$\begin{split}
\int_\gamma (g^1)^* \omega - \omega &= (-1)^k\int_0^1 \int_\gamma H_s^* i_X d\omega \ ds - (-1)^{k-1} \int_0^1 \int_{\partial\gamma}H_s^* i_X \omega \ ds\\
&= (-1)^k \int_0^1 \int_\gamma H_s^* (i_X d\omega + di_X \omega) \ ds .
\end{split}$$
Hence
$$\int_0^1 \int_\gamma H_s^* \mathcal L_X \omega \ ds = \int_0^1 \int_\gamma H_s^* (i_X d\omega + di_X \omega) \ ds .$$
Note that if we consider the homotopy $G_t$, which is the restriction of $H\gamma$ to $[0, t]\times \sigma$, then similarly we have
$$\int_0^t \int_\gamma H_s^* \mathcal L_X \omega \ ds = \int_0^t \int_\gamma H_s^* (i_X d\omega + di_X \omega) \ ds $$
for all $t >0$. Thus (using $H_0 = \gamma$)
$$ \int_\gamma \mathcal L_X \omega= \int_\gamma (i_X d\omega + di_X \omega).$$
As the above holds for all $\gamma$, we have
$$ \mathcal L_X \omega = i_X d\omega + di_X \omega.$$
It remains to prove $(1)$. Let $\alpha$ be an $r$ form. By definition, (call $H\gamma = F$)
$$\int_{H\gamma} \alpha = \int_{I \times \sigma} F^*\alpha$$
if we write in local coordinates $(t = x^0, x_1, \cdots, x_{r-1})$ in $I\times \sigma$ and $(y^1, \cdots, y^m)$ of $M$.
$$\begin{split}
F^* \alpha (t, x)&= F^* \left(\sum_{|I| = r} \alpha_{i_1\cdots i_{r}} dy^{i_1} \wedge \cdots \wedge dy^{i_{r}}\right) \\
&= \sum_{|I| =r} \alpha_{i_1\cdots i_{r}} (F(t, x)) dF^{i_1} \wedge \cdots \wedge dF^{i_{r}} \\
\end{split}$$
Note that
$$dF^{i_j} = \frac{\partial F^{i_j}}{\partial t} (t, x) dt + \sum_{p=1}^{r-1} \frac{\partial F^{i_j}}{\partial x^p} (t, x) dx^p = X^{i_j}(F(t,x)) dt + H_t^* dy^{i_j} $$
as $X = \frac{\partial F}{\partial t}$ by definition of $F = H\gamma$. Thus
$$\begin{split}
F^* \alpha&= \sum_{|I| =r} \alpha_{i_1\cdots i_{r}} \left( X^{i_1} dt + H_t^* dy^{i_1}\right) \wedge \cdots \wedge \left( X^{i_{k+1}} dt + H_t^* dy^{i_{r}}\right) \\
&= \sum_{|I| =r} \alpha_{i_1\cdots i_{r}} \ dt \wedge H_t^* \left( \sum_{j=1}^{r}(-1)^{j-1} X^{i_j} dy^{i_1} \wedge \cdots \wedge \hat{dy^{i_j}} \wedge \cdots \wedge dy^{i_{r}}\right) \\
&= \sum_{|I| =k+1} \alpha_{i_1\cdots i_{r}} \ dt \wedge H_t^*\left( i_X (dy^{i_1} \wedge \cdots \wedge dy^{i_{r}})\right) \\
&= dt \wedge H_t^*(i_X \alpha) = (-1)^{r-1} H_t^* (i_X\alpha) \wedge dt.
\end{split}$$
Thus
$$ \int_{H\gamma} \alpha = (-1)^{r-1} \int_{I\times \sigma} H_t^*(i_X \alpha) \wedge dt = (-1)^{r-1} \int_0^1 \int_\gamma H_t^*(i_X \alpha) dt.$$
I'm going to use capital letters for the vector fields and $\mathcal{L}$ for the Lie derivative.
By Cartan's formula, we have
$$(\mathcal{L}_X\alpha)(Y, Z) = (d(i_X\alpha))(Y, Z) + (i_X(d\alpha))(Y, Z) = (d(i_X\alpha))(Y, Z) + (d\alpha)(X, Y, Z).$$
If $\beta$ is a one-form, $(d\beta)(Y, Z) = Y\beta(Z) - Z\beta(Y) - \beta([X, Y])$, so for $\beta = i_X\alpha$ we see that
\begin{align*}
(d(i_X\alpha))(Y, Z) &= Y(i_X\alpha)(Z) - Z(i_X\alpha)(Y) - (i_X\alpha)([Y, Z])\\
&= Y\alpha(X, Z) - Z\alpha(X, Y) - \alpha(X, [Y, Z])\\
&= Y\alpha(X, Z) - Z\alpha(X, Y) + \alpha([Y, Z], X).
\end{align*}
As for the other term, we have
\begin{align*}
&\ (d\alpha)(X, Y, Z)\\
=&\ X\alpha(Y, Z) - Y\alpha(X, Z) + Z\alpha(X, Y) - \alpha([X, Y], Z) + \alpha([X, Z], Y) - \alpha([Y, Z], X).
\end{align*}
Combining the two, we obtain
\begin{align*}
(\mathcal{L}_X\alpha)(Y, Z) &= X\alpha(Y, Z) - \alpha([X, Y], Z) + \alpha([X, Z], Y)\\
&= X\alpha(Y, Z) - \alpha(\mathcal{L}_X Y, Z) + \alpha(\mathcal{L}_X Z, Y)\\
&= X\alpha(Y, Z) - \alpha(\mathcal{L}_X Y, Z) - \alpha(Y, \mathcal{L}_X Z).
\end{align*}
Now note that if $f$ is a smooth function, then $\mathcal{L}_Xf = d(i_Xf) + i_X(df) = d(0) + df(X) = Xf$. In particular,
$$\mathcal{L}_X(\alpha(Y, Z)) = X\alpha(Y, Z).$$
Therefore
$$(\mathcal{L}_X\alpha)(Y, Z) = \mathcal{L}_X(\alpha(Y, Z)) - \alpha(\mathcal{L}_X Y, Z) - \alpha(Y, \mathcal{L}_X Z)$$
and hence
$$\mathcal{L}_X(\alpha(Y, Z)) = (\mathcal{L}_X\alpha)(Y, Z) + \alpha(\mathcal{L}_X Y, Z) + \alpha(Y, \mathcal{L}_X Z).$$
More generally, if $\alpha$ is a contravariant $k$-tensor (for example, a $k$-form),
$$\mathcal{L}_X(\alpha(Y_1, \dots, Y_k)) = (\mathcal{L}_X\alpha)(Y_1, \dots, Y_k) + \sum_{i=1}^k\alpha(Y_1, \dots, Y_{i-1}, \mathcal{L}_XY_i, Y_{i+1}, \dots, Y_k).$$
In fact, one usually defines the Lie derivative to act on tensors by precisely this rule.
An algebraic way of expressing the above equation is: the Lie derivative obeys the Leibniz rule with respect to contractions. In fact, this can be used as one of the axioms that uniquely determines the action of the Lie derivative on tensor fields - see here.
Best Answer
$i_X\circ L_X=i_X\circ (i_X\circ d+d\circ i_X)=i_X\circ d\circ i_X$ since $i_X\circ i_X=0$.
$L_X\circ i_X=(i_X\circ d +d\circ i_X)\circ i_X=i_X\circ d\circ i_X$. You deduce that $i_X\circ L_X=L_X\circ i_X$.