The derivative of the action functional and Lie derivative 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?
 A: After pondering for a while, I can't make sense of the expression $x^*(\mathcal{L}_Y \alpha)$ in general. The way I managed to circumvent this was to look at the cylinder $S^1 \times (-\epsilon, \epsilon)$ itself.
Let $H: S^1 \times (\epsilon, \epsilon) \to M$ be the variation under question (i.e., $H(\cdot,s)=x_s$) and let $j_s:S^1 \to S^1 \times (-\epsilon,\epsilon)$ be the inclusion $j_s(t)=(t,s)$. Then
$$x_s=H\circ j_s =H \circ \Phi_s\circ j_0,$$
where $\Phi_s$ is the flow of the vector field $(0,\partial/\partial s$).
Thus,
\begin{align*}
\left.\frac{\partial }{\partial s}\right|_{s=0}x_s^*\alpha&=\left.\frac{\partial }{\partial s}\right|_{s=0}j_0^*\Phi_s^*H^*\alpha \\
&=j_0^*\mathcal{L}_{\partial/\partial s}H^*\alpha.
\end{align*}
I think the above is the correct computation. The objective of the computation on OP is to get at the fact that
$$\left.\frac{d}{ds}\right|_{s=0}\mathcal{A}(x_s)=\int_0^1d\alpha(Y,\dot{x}). $$
We now see that we arrive at this final result with the above computation. Indeed,
\begin{align*}
\left.\frac{d}{ds}\right|_{s=0}\mathcal{A}(x_s)&=\int_{S^1}j_0^*\mathcal{L}_{\partial/\partial s}H^*\alpha \\
&=\int_{S^1}j_0^*(d\iota_{\partial/\partial_s}H^*\alpha+\iota_{\partial/\partial_s}dH^*\alpha) \\
&=\int_{S^1}j_0^*\iota_{\partial/\partial_s}H^*d\alpha \\
&=\int_{S^1}j_0^*\iota_{\partial/\partial_s}d\alpha(H_*(\cdot),H_*(\cdot)) \\
&=\int_{S^1}j_0^*d\alpha\left(\frac{\partial}{\partial s}H,H_*(\cdot)\right) \\
&=\int_{S^1}d\alpha\left(Y,H_*(\cdot)\right) \\
&=\int_0^1 d\alpha\left(Y,\dot{x}\right).
\end{align*}
A: You need to use the fact that $x_s$ is $\Phi_s \circ x$, where $\Phi_s : M\to M$ is the flow of the vector field $Y$. Then
$$\frac{\partial}{\partial s} x_s^* \alpha \bigg|_{s=0} =  \frac{\partial}{\partial s}(\Phi_s \circ x)^* \alpha \bigg|_{s=0}= \frac{\partial}{\partial s} x^* (\Phi_s^* \alpha) \bigg|_{s=0}= x^* \left(\frac{\partial}{\partial s} \Phi_s^* \alpha\bigg|_{s=0} \right)= x^* (L_Y \alpha). $$
A: 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.
