Differential operators are coKleisli morphisms of the jet co-monad The following statement may be "well known but not well known enough", and my question is which reference would state it explicitly:
The construction of Jet bundles is a comonad on suitable bundles over a given base $X$. A differential operator $D \colon \Gamma_X(E_1) \to \Gamma_X(E_2)$ is equivalently a morphism of bundles of the form $\tilde D \colon Jet(E_1) \to E_2$ (the associated differential operator is $\phi \mapsto \tilde D\circ j_\infty(\phi)$). Under this identification, the composition $D_2\circ D_1$of two differential operators is given by the coKleisli composite, namely by the composite
$$
  Jet(E_1) \to Jet(Jet(E_1)) \stackrel{Jet(\tilde D_1)}{\longrightarrow} Jet(E_2) \stackrel{\tilde D_2}{\longrightarrow} E_3
$$
where the first morphism is the coproduct of the Jet comonad.
In the context of differential geometry, the article 


*

*Joseph Krasil'shchik, Alexander Verbovetsky, Homological Methods in Equations of Mathematical Physics (arXiv:math/9808130)


has essentially all the ingredients for this statement (p. 13,14,17), but does not make it explicit in the above form. In the algebraic context there may be more references that state it this way or almost state it this way (the statement that $Jet$ is a comonad for sure, but how about the differential operators being coKleisli maps?). Whatever references you are aware of, please let me know.
 A: It turns out that the statement in question appeared in full beauty in 


*

*Michal Marvan, A note on the category of partial differential equations, in Differential geometry and its applications, Proceedings of the Conference in Brno August 24-30, 1986,  


There in fact is given the more general statement that the full Eilenberg-Moore category of the jet comonad is equivalent to that of partial differential equations in variables in the given base space.
Of course the Kleisli category in question is a full subcategory of that, and the statement that I asked for is the special case of the above made fully explicitly in 


*

*Michal Marvan, section 1.1 of On Zero-Curvature Representations of Partial Differential Equations, 1993


Now since this statement about differential operators is true, of course every other true statement about differential operators that one finds elsewhere in the literature will be compatible with it, but Marvan is the only author I am aware of who explicitly states the neat general abstract statement that I asked for. It's one of those abstract trivialities, if you wish, that are worth making explicit.
On the other hand, Marvan checks that $\mathrm{Jet}$ is a comonad directly, without realizing it as the base change comonad $\mathrm{Jet}_X \simeq i^\ast_X (i_X)_\ast$ along the unit of the de Rham stack monad $i_X \colon X \to X_{dR}$. For that statement of course one needs to be in a model of synthetic differential geometry in order for $X_{dR}$ to exist, such as (but not necessarily) algebraic geometry.
In the algebraic context, the reference stating explicitly $\mathrm{Jet} \simeq i^\ast_X (i_X)_\ast $ that I am aware of is


*

*Jacob Lurie, above prop. 0.9 of Notes on crytsals and algebraic D-modules, 2009


Here it is maybe noteworthy that from this it follows immediately that the Jet-comonad is right adjoint to the "infinitesimal disk bundle" monad $T_{\mathrm{inf}} := i_X^\ast (i_X)_!$. 
This adjunction $(T_{\mathrm{inf}} \dashv \mathrm{Jet})$ itself, without however its origin in the base change adjoint triple of $i_X$, was observed in the context of synthetic differential geometry in 


*

*Anders Kock, above prop. 2.2 of Formal manifolds and synthetic theory of jet bundles, 1980


I'll just add that one immediate neat implication of making the general abstract comonad theory behind this explicit is that it gives in full generality that for any topos (or $\infty$-topos) $\mathbf{H}$ equipped with an "infinitesimal shape modality" $X\mapsto \Im X = X_{dR}$, then since $\mathrm{Jet} := (i_X)^\ast (i_X)_\ast \colon \mathbf{H}_{/X} \to \mathbf{H}_{/X}$ is a right adjoint, a standard fact (here for toposes, here maybe for $\infty$-toposes) gives that its EM-category, hence the category $\mathrm{PDE}(X)$ of PDEs in $X$ is itself a topos, sitting by a geometric morphism
$$\mathrm{PDE}(X) \stackrel{{\longleftarrow}}{\longrightarrow} \mathbf{H}_{/X}$$
over $\mathbf{H}$, whose direct image is (non-fully) that co-Kleisli category of bundles over $X$  with differential operators between them.
Now moreover, using the infinitesimal shape modality $\Im$ it also follows that on $\mathbf{H}_{/X}$ there is an étalification coprojection $Et$ (here). This is such that when $E \in \mathbf{H}$ is a coefficient for a differential cohomology theory (e.g. any stable object in the case that $\mathbf{H}$ is cohesive (here)), then $Et(X^\ast E)$ has the interpretation of being the "bundle of $E$-connections" over $X$. So with the above geometric morphism we may transfer all this to the topos $\mathrm{PDE}(X)$ along the composite
$$
  \iota \colon \mathbf{H} \stackrel{X^\ast}{\longrightarrow} \mathbf{H}_{/X}\stackrel{Et}{\longrightarrow} \mathbf{H}_{/X} \stackrel{}{\longrightarrow} \mathrm{PDE}(X)
  \,.
$$
Unwinding what this all means, in view of Marvan's insight above, it follows, I think, that for any bundle $F \in \mathbf{H}_{/X} \to \mathrm{PDE}(X)$ then maps in $\mathrm{PDE}(X)$ of the form
$$
  F \longrightarrow \iota E
$$
are equivalently horizontal "$E$-differential forms" on the jet bundle of $F$, in the sense of variational bicomplex theory.  
That might count as some non-trivial mileage gained out of making explicit the statement that I was asking for.
A: This is not really an answer since I don't know a good reference, just to elaborate for readers (as the OP mentions) that this description of differential operators is readily equivalent to one that might be more familiar to algebraic geometers. Namely the jet comonad is precisely the descent comonad from X to the de Rham functor of X, namely the quotient of X by the equivalence relation (de Rham groupoid) given by the formal neighborhood of the diagonal --- in other words $J(V)=\pi^*\pi_* V$ where $\pi:X\to X_{dR}$ is the quotient map. Indeed the coalgebra of jets (the dual to the algebra of differential operators) is the groupoid coalgebra of the de Rham groupoid. Comodules for this comonad are sheaves on X with descent data, i.e. (with suitable hypotheses) sheaves on the de Rham space. Thus the definition is expressing the identification of modules for differential operators with sheaves on the de Rham space --- or in Grothendieck terminology, stratifications on X / sheaves on the infinitesimal site (in characteristic 0). This is the modern definition of $D$-modules, eg most modernly in the work of Gaitsgory-Rozenblyum, but going back to Grothendieck's work on crystals. 
The question asks about arbitrary differential operators between vector bundles -- but note (following e.g. M. Saito, Induced D-modules and differential complexes, Bull. Soc. Math. France 117 (1989), 361–387.) that a differential operator between vector bundles $V\to W$ is the same thing as a map of (right) D-modules between the corresponding induced D-modules $V\otimes D\to W\otimes D$ or between the sheaves of jets $J(V)=\pi^*\pi_* V \to J(W)=\pi^*\pi_* W$, whence the coKleisli interpretation.
