How to look at the Lie derivative as a partial differential operator?

There is another way to show something is a linear partial differential operator. One can characterize partial differential operators in a coordinate-free manner using commutators:

Let $E\to M$ and $F\to M$ be two vector bundles over $M$. Consider the set, $\mathcal{L}(\Gamma(E),\Gamma(F))$, of $\mathbf{R}$-linear operators between sections of $E$ and sections of $F$. Now define $\mathbf{PDO}_0(E,F):=\mathrm{hom}(E,F)$; note that this coincides with the local definition of an order $0$ differential operator. The key observation is that by considering a differential operator $D$ of order $k$ and looking at the commutator $[D,f](\sigma):=D(f\sigma)-fD(\sigma)$ that the degree of the differential operator decreases by $1$. In some sense, it cancels out the Leibniz rule.

For example, if $D\in\mathbf{PDO}_0(E,F)$ we have for any $f\in\mathscr{C}^\infty(M)$ and $\sigma\in\Gamma(E)$ that $[D,f](\sigma)=D(f\sigma)-fD(\sigma)=0$.

Now we can equivalently give an inductive definition of order $k$ partial differential operators as $$\mathbf{PDO}_k(E,F):=\{D\in\mathcal{L}(\Gamma(E),\Gamma(F))~:~[D,f]\in\mathbf{PDO}_{k-1}(E,F)~\text{ for all }f\in\mathscr{C}^\infty(M)\}.$$ It shouldn't be too difficult to show the equivalence of this definition with that given by factorization through the $k^\mathrm{th}$ Jet space of $E$. The advantage of this definition is that it makes it very easy to show that a operator is a order $k$ partial differential operator $-$ we simply have to show that $k+1$ successive commutators of the operator is zero.

Now given this characterization of partial differential operators, we see immediately see that the Lie derivative is a differential operator of order $1$: Let $E=F=T^{s}_{r}M$ denote the $r$-times contravariant and $s$-times covariant tensor product. Fix a vector field $X\in\Gamma(T^0_1M)$, and consider any smooth functions $f,g\in\mathscr{C}^\infty(M)$ and a tensor field $T\in\Gamma(E).$ We have by using the Leibniz rule of the Lie derivative, \begin{align*} [[\mathscr{L}_X,f],g](T)& =[\mathscr{L}_X,f](gT)-g[\mathscr{L}_X,f](T)\\& =\mathscr{L}_X(fgT)-f\mathscr{L}_X(gT)-g\mathscr{L}_X(fT)+fg\mathscr{L}_X(T)\\ & = (\mathscr{L}_X(fg))T+(fg)\mathscr{L}_X(T)-f\mathscr{L}_X(g)T-(fg)\mathscr{L}_X(T)\\ &\qquad -g\mathscr{L}_X(f)T-(fg)\mathscr{L}_X(T)+(fg)\mathscr{L}_X(T)\\ & = \mathscr{L}_X(fg)T-f\mathscr{L}_X(g)T-g\mathscr{L}_X(f)T\\ & = 0. \end{align*} EDIT: I made the mistake of considering the differential operator $\mathscr{L}_X:\Gamma(E)\to\Gamma(E)$ for a fixed vector field $X\in\Gamma(TM)$, instead of what was asked in the question: Fix a tensor field $\omega\in T^s_rM$ for some $r,s\in\mathbf{N}$, and consider the operator $D_\omega:\Gamma(TM)\to\Gamma(T^s_rM)$ given by $X\mapsto\mathscr{L}_X\omega$. Note that in the case $\omega\in\Omega_k(M)$ is a differential $k$-form, it is very straightforward to see that $D_\omega\in\mathbf{PDO}_1(TM,\bigwedge^kTM^*)$ $-$ recall the identity $$\mathscr{L}_{fX}\omega=f\mathscr{L}_X\omega + df\wedge\iota_X\omega,\qquad f\in\mathscr{C}^\infty(M).$$ Now fix $X\in\Gamma(TM)$ and $f,g\in\mathscr{C}^\infty(M).$ A direct computation shows that \begin{align*} [[D_\omega,f],g](X) & = [D_\omega,f](gX) - g[D_\omega,f](X) \\ & = D_\omega(fg X) - fD_\omega(gX) - gD_\omega(fX) + fgD_\omega(X) \\ & = \mathscr{L}_{fgX}\omega - f\mathscr{L}_{gX}\omega - g\mathscr{L}_{fX}\omega + fg\mathscr{L}_{X}\omega \\ & = fg\mathscr{L}_{X}\omega +d(fg)\wedge\iota_X\omega - fg\mathscr{L}_{X}\omega - fdg\wedge\iota_X\omega \\ & \qquad -fg\mathscr{L}_{X}\omega - gdf\wedge\iota_X\omega +fg\mathscr{L}_{X}\omega \\ & = d(fg)\wedge\iota_X\omega - fdg\wedge\iota_X\omega - gdf\wedge\iota_X\omega \\ & = 0. \end{align*} So we have demonstrated that $D_\omega\in\mathbf{PDO}_1(TM,\bigwedge^kTM^*).$

I found it to be a little bit messier to consider an arbitrary $(r,s)$-type tensor. I don't know if there are any nifty identities for $\mathscr{L}_{fX}T$ for smooth functions $f$, vector fields $X$, and tensor fields $T$. Nevertheless there isn't really anything deep going on in the following computations. Recall that we can consider $(r,s)$-type tensors as multilinear maps from $(\Omega_1(M))^r\times(\Gamma(TM))^s\to\mathbf{R}$. Under this natural isomorphism it isn't hard to show \begin{align*} (\mathscr{L}_XT)(\alpha^1,\dots,\alpha^r,Y_1,\dots,Y_s)&=X(T(\alpha^1,\dots,\alpha^r,Y_1,\dots,Y_s))-\sum_{j=1}^{r}T(\alpha^1,\dots,\mathscr{L}_X\alpha^j,\dots,\alpha^r,Y_1,\dots,Y_s) \\&\quad - \sum_{k=1}^{s}T(\alpha^1,\dots,\alpha^r,Y_1,\dots,\mathscr{L}_XY_k,\dots,Y_s), \end{align*} where $X\in\Gamma(TM)$, $T\in\Gamma(T^s_rM)$, $\alpha_j\in\Omega_1(M)$, and $Y_k\in\Gamma(TM).$ For simplicity, I will write $\vec{\alpha}=(\alpha^1,\dots,\alpha^r)$ and $\vec{Y}=(Y_1,\dots,Y_s).$ In particular, we see that for some $f\in\mathscr{C}^\infty(M)$ that \begin{align*} (\mathscr{L}_{fX}T)(\vec{\alpha},\vec{Y}) & = f(\mathscr{L}_XT(\vec{\alpha},\vec{Y})) - \sum_{j=1}^{r}T(\alpha^1,\dots,df\wedge\iota_X\alpha^j,\dots,\alpha^r,\vec{Y}) + \sum_{k=1}^{s}df(Y_k)T(\vec{\alpha},Y_1,\dots,X,\dots,Y_s). \end{align*} Now we compute for arbitrary $f,g\in\mathscr{C}^\infty(M)$, $T\in\Gamma(T^s_rM)$, $\vec{\alpha}\in(\Omega_1(M))^r$, and $\vec{Y}\in(\Gamma(TM))^s$: \begin{align*} & ([[D_T,f],g](T))(\vec{\alpha},\vec{Y}) = D_T(fg X) - fD_T(gX) - gD_T(fX) + fgD_T(X) \\ & = fg(\mathscr{L}_X(T)(\vec{\alpha},\vec{Y})) - \sum_{j=1}^{r}T(\alpha^1,\dots,d(fg)\wedge\iota_X\alpha^j,\dots,\alpha^r,\vec{Y}) + \sum_{k=1}^{s}d(fg)(Y_k)T(\vec{\alpha},Y_1,\dots,X,\dots,Y_s) \\ & \quad - fg(\mathscr{L}_X(T)(\vec{\alpha},\vec{Y})) + f\sum_{j=1}^{r}T(\alpha^1,\dots,dg\wedge\iota_X\alpha^j,\dots,\alpha^r,\vec{Y}) - f\sum_{k=1}^{s}dg(Y_k)T(\vec{\alpha},Y_1,\dots,X,\dots,Y_s) \\ & \quad - fg(\mathscr{L}_X(T)(\vec{\alpha},\vec{Y})) + g\sum_{j=1}^{r}T(\alpha^1,\dots,df\wedge\iota_X\alpha^j,\dots,\alpha^r,\vec{Y}) - g\sum_{k=1}^{s}df(Y_k)T(\vec{\alpha},Y_1,\dots,X,\dots,Y_s) \\ & \quad + fg(\mathscr{L}_X(T)(\vec{\alpha},\vec{Y})) \\ & = 0. \end{align*} To get that it's equal to zero I simply used the multilinearity of all of the mappings involved and that $d(fg)=fdg + gdf$. I tried to make it easier to read by having each line include only one expansion of the terms above In particular, since $\vec{\alpha}$ and $\vec{Y}$ were arbitrary we deduce that $[[D_T,f],g]=0,$ and so $D_T\in\mathbf{PDO}_1(TM,T^s_rM)$ as desired.

In light of all of the above examples, we see that this characterization does indeed make it very straightforward to show that an operator is indeed an order $k$ partial differential operators. Another example: it's easy to use the above to see that the exterior derivative is a PDO of order at most $1.$

The statement can be proven by induction.

Let $\omega_1$ and $\omega_2$ be tensors which satisfy the property. That is, the operators $X\mapsto\mathcal{L}_X\omega_1$ and $X\mapsto\mathcal{L}_X\omega_2$ are differential operators of order $1$. By the Leibniz rule,$$\mathcal{L}_X(\omega_1\otimes\omega_2)=(\mathcal{L}_X\omega_1)\otimes\omega_2+\omega_1\otimes(\mathcal{L}_X\omega_2),$$and it follows immediately that $X\mapsto\mathcal{L}_X(\omega_1\otimes\omega_2)$ is also a differential operator of order $1$.

You already know the statement is true for $1$-forms. It is also true for vector fields, as $\mathcal{L}_XY=[X,Y]$. Since the tensor algebra is generated by vector fields and $1$-forms, we are done.