The material mentioned by Yuri in the comment can also be found in "Topics in differential geometry", available here
http://www.mat.univie.ac.at/~michor/dgbook.pdf
by the same author, section 8.15. It may be a bit simpler than "Natural operations in differential geometry".
I was looking at the material right now. From what I have understood you need to be able to lift the flow of the vector field on the base manifold to the total space of the vector bundle in such a way that
Given all this the Lie derivative can be defined in the usual why as the $\mathrm{d}/\mathrm{d}t \,|_{t=0}\, \phi^*_t$ where $\phi$ is the lifted flow.
Of course for this to be useful you need to have some preferred/unique way of lifting the flow on the base to the total space.
It is completely correct.
The given formula for the Lie derivative of a one-form follows from Cartan's identity:
$$\mathcal{L}_X\alpha = i_X(d\alpha) + d(i_X\alpha).$$
Let's compute $\mathcal{L}_X\alpha$ using this identity and check we get the same result.
First, we can compute the exterior derivative of $\alpha$:
$$d\alpha = 2ydy\wedge dx + 2xdx\wedge dy = (2x-2y)dx\wedge dy.$$
Using the formula $i_X(\beta\wedge\gamma) = (i_X\beta)\wedge\gamma + (-1)^{|\beta|}\beta\wedge(i_X\gamma)$, we can compute the interior product of $d\alpha$:
\begin{align*}
i_X(d\alpha) &= (2x - 2y)i_X(dx)\wedge dy - (2x - 2y)dx\wedge i_X(dy)\\
&= (2x - 2y)dx(X)dy - (2x - 2y)dy(X)dx\\
&= (2x - 2y)dx\left(\frac{\partial}{\partial x} + xy\frac{\partial}{\partial y}\right)dy - (2x-2y)dy\left(\frac{\partial}{\partial x} + xy\frac{\partial}{\partial y}\right)dx\\
&= (2x-2y)dx\left(\frac{\partial}{\partial x}\right)dy - (2x-2y)(xy)dy\left(\frac{\partial}{\partial y}\right)dx\\
&= (2xy^2 - 2x^2y)dx + (2x-2y)dy.
\end{align*}
For the second term, the interior product of $\alpha$ is
\begin{align*}
i_X\alpha &= y^2i_X(dx) + x^2i_X(dy)\\
&= y^2dx(X) + x^2dy(X)\\
&= y^2dx\left(\frac{\partial}{\partial x} + xy\frac{\partial}{\partial y}\right) + x^2dy\left(\frac{\partial}{\partial x} + xy\frac{\partial}{\partial y}\right)\\
&= y^2dx\left(\frac{\partial}{\partial x}\right) + x^2(xy)dy\left(\frac{\partial}{\partial y}\right)\\
&= y^2 + x^3y
\end{align*}
and taking the exterior derivative, we have
$$d(i_X\alpha) = 2ydy + 3x^2ydx + x^3dy = 3x^2ydx + (x^3 + 2y)dy.$$
Combining, we see that
\begin{align*}
\mathcal{L}_X\alpha &= i_X(d\alpha) + d(i_X\alpha)\\
&= (2xy^2 - 2x^2y)dx + (2x - 2y)dy + 3x^2ydx + (x^3 + 2y)dy\\
&= (2xy^2 + x^2y)dx + (2x+x^3)dy
\end{align*}
which is the same result you obtained.
Best Answer
First of all, it is important to notice that $\theta$ and $\varphi$ stand for the same thing here. In her paper, the authors refers to $\theta$ when $\Bbb S^1$ is thought of as the coordinate on the group $\Bbb S^1$, acting on a manifold $M$. She refers to $\varphi$ when it is thought of as the manifold $\Bbb S^1$. For the action by translation of the circle on itself, they coincide and hence, $J = -\frac{\partial}{\partial \theta} = - \frac{\partial}{\partial \varphi}$.
$\Omega^1(\Bbb S^1)$ is a $\mathcal{C}^{\infty}(\Bbb S^1)$-module of rank $1$. Since $d\varphi$ is nonzero everywhere, $\{d\varphi\}$ is a basis of $\Omega^1(\Bbb S^1)$, that is, any differential form of degree $1$ on the circle is of the form $gd\varphi$.
Since $\mathcal{L}(J)$ is a derivation, $\mathcal{L}(J)[g d\varphi] = (\mathcal{L}(J)g)d\varphi + g \mathcal{L}(J) d\varphi$. You can check that $\mathcal{L}(J) d\varphi = 0$ since the action of $\mathbb{S}^1$ onto itself is by translation and hence the coordinate $\varphi$ is invariant under this action. Therefore, $\mathcal{L}(J)[gd\varphi] = 0 \iff \mathcal{L}(J)g = 0$. The problem reduces to solving the equation $\mathcal{L}(J)g = 0$, which is in fact just the equation $dg = 0$, and since $\Bbb S^1$ is connected, its solutions are constant maps.
Comment The author isn't a 'he', but a 'she', since Michèle Vergne is a woman although she's a mathematician ;)