When computing the "formal" derivative of the action functional $$\mathcal{A}\colon C^{\infty}(S^1,M) \to \mathbb{R}$$ $$x \mapsto \int_{S^1}x^* \alpha $$ on a contact manifold $M$ with contact form $\alpha$, we pick a variation $x_s$ of $x$ with variational vector field $Y$ and aim to compute $$\frac{d}{ds}\mathcal{A}(x_s).$$ In the computations I see, the first steps are \begin{align*} \left(\frac{d}{ds}\mathcal{A}(x_s)\right)_{|_{s=0}}&=\left(\frac{d}{ds}\int_{S^1}x_s^* \alpha\right)_{|_{s=0}} \\ &= \int_{S^1} \left(\frac{\partial}{\partial s}x_s^* \alpha\right)_{|_{s=0}} \\ &\stackrel{!}{=} \int_{S^1} x^* \mathcal{L}_Y \alpha \\ &=\cdots \end{align*} and from then on I can follow. However, I don't see why the Lie derivative appears here. The $x_s^*$ are pull-backs by maps of the circle. I don't see the derivative of the pull-back of a flow on $M$, only the derivative of the pull-back of the maps $S^1 \to M$. In fact, I don't see the flow of $Y$ appearing.
My question is: How do I see why the equation $(!)$ holds?
It is not immediately clear to me how to interpret the Lie derivative $\mathcal{L}_{Y} \alpha$, since $Y$ is only a vector field along $x$. One way to deal with this is to pull everything back to the domain of the variation. We see the variation $x_s$ as a map $\xi: S^1 \times (-\epsilon, \epsilon) \to M$, with $\xi( \cdot, s) = x_s$. On the cylinder $S^1 \times (-\epsilon, \epsilon)$ there is the vector field $\frac{\partial}{\partial s}$ which points along the $(-\epsilon, \epsilon)$ direction. Denote by $\Phi_t$ its flow. Let $i_s: S^1 \to S^1\times (-\epsilon, \epsilon)$ be the inclusion at level $s$. Then of course $i_s = \Phi_s \circ i_0$, and $x_s = \xi \circ i_s = \xi \circ \Phi_s \circ i_0$.
Now we can compute: \begin{align} \frac{d}{ds}\bigg|_{s=0} \, x_s^* \alpha &= \frac{d}{ds}\bigg|_{s=0} \, i_0^* \; \Phi_s^* \;\xi^* \alpha \\ &= i_0^* \;\frac{d}{ds}\bigg|_{s=0} \, \Phi_s^* \;\xi^* \alpha \\ &= i_0^* \;\mathcal{L}_{\frac{\partial}{\partial s}} \xi^*\alpha. \end{align}
So here instead of the (a priori) ill-defined $\mathcal{L}_Y \alpha$ we have the derivative of $\xi^* \alpha$ in the direction of $\frac{\partial}{\partial s}$, which pretty much means "differentiate $\alpha$ in the direction transverse to the variation", as we would want intuitively. From there, the next step is to use Cartan's formula, but I think you had this part figured out.