How to look at the Lie derivative as a partial differential operator? In Gromov's book Partial Differential Relations, a partial differential operator (PDO) is defined as follows. Given a fibration $X\to V$ and a vector bundle $G\to V$ over a manifold $V$, a PDO of order $k$ is a map $D:\Gamma(X)\to \Gamma(G)$, satisfying $D(x)=\Delta\circ j^k_x$ for every $x\in \Gamma(X)$, where $\Delta:X^{(k)}\to G$ is a bundle map. Here $X^{(k)}$ is the $k^{\text{th}}$ jet bundle of sections of $X\to V$ and $j^k_x$ is the $k^{\text{th}}$ jet of $x$.
Now, say $G$ is some tensor bundle on $V$ and $X=TV$. Then for a fixed section $\omega\in\Gamma(G)$, we have an operator,
\begin{align*}
D_\omega:\Gamma(TV)&\to \Gamma(G)\\
x&\mapsto \mathcal{L}_x(\omega)
\end{align*}
where $\mathcal{L}_x(\omega)$ is the Lie derivative of the tensor $\omega$ along the vector field $x$.
How do I see that this $D_\omega$ is indeed a partial differential operator of order 1, as defined by Gromov? In particular, I would like to find out the bundle map $\Delta=\Delta_w$ such that $D_\omega(x)=\Delta\circ j^1_x$.
For the case $G=\Lambda^k T^*V$, I can use the Cartan's formula and get, $$D_\omega(x)=\mathcal{L}_x(\omega)=d\iota_x(\omega)+\iota_xd\omega$$ I do know that the exterior derivative map $d$ can be seen as a PDO. But I still can't see how to get to the map $\Delta$. And I'm in complete loss when $G$ is an arbitrary tensor bundle, since in that case the Lie derivative is defined using flows.
Any help regarding this appreciated.
 A: 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.
A: 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.$
A: Here, $X$ is $\Lambda ^kT^*V$, $G=\Lambda ^kT^*V$, and in local coordinates, one see that $L_X$ is a bundle map, once $X$ is fixed, namely the coordiantes of $L_X \omega$ linearly depend on the one jet of the coordinates of $\omega$. For instance if $\omega= \sum \omega _idx^i$ then$L_X \omega= \sum _{i\not = j} X_j \partial _j \omega _i dx^i\wedge dx ^j$
