Shuffle (co-)multiplication and generalized Leibniz formula in tensor calculus The headline already says it: Is anybody (except me, UPDATE: plus Gavrilov) aware of this formula for higher total covariant derivatives of tensor products?
It is the simplest application of the commutative shuffle product in differential geometry. (The anticommutative shuffle product is sometimes invoked in a formula for the wedge product of differential forms.)
UPDATE x: In elementary calculus Leibniz' formula is
$(fg)^{(n)} \;=\; \sum_{k=0}^n {n \choose k} f^{(n)}g^{(n-k)}$
for the n-th derivative of a product of functions. There's a generalization to multivariate calculus, with a multi variant of the binomial coefficients. In tensor calculus they are replaced by a more complicated combinatorial thing: shuffle multiplication.
UPDATE: I have asked purely verbally because there are many ways to do "this formula", depending on the tensor calculus formalism and the linear algebra. One could as well use the adjoint of shuffle multiplication, i.e. the Lie-Hopf comultiplication (or, "unshuffle comultiplication").
UPDATE 2: Statement and proof in my formalism is presented in my 2nd answer below.
UPDATE 3: Meanwhile I'm convinced that my formalism of total covariant derivative, presented below, does it the wrong way: Cotangent spaces should be tensored to the left, not to the right. But I'm in good company (from Kobayashi-Nomizu to Jost). There are 3 reasons for doing it to the left:


*

*R.Palais 1965: Seminar on the Atiyah-Singer Index Theorem IV §9: Consistency with jet bundles.

*Pure multilinear algebra consistency considerations of canonical isomorphisms.

*Gavrilov's cocycle identity should come without twist. Here's a well-known special case in wrong notation: $id_E\otimes X\otimes Y\cdot\nabla^{2,E}a=\nabla_Y^E\nabla_X^Ea-\nabla_{\nabla_Y^{T^\ast M}X}^Ea$.


Given this sorry example of the state of the art of tensor calculus, I'm no longer surprised about the scandal that it took until 2009 for something as naturally self-suggesting as Leibniz' formula to surface in the literature (and then "buried first-class" in a Siberian journal.)
-- Having stuck out my neck this far, I will try to copy my two answers into the question, so I can set out a bounty offering most of my reputation points :-)

Answer 1:

Alexey V. Gavrilov: The Leibniz Formula for the Covariant Derivative and Some of Its Applications, Siberian Advances in Mathematics 22 (2), 80-94 (2012) Springer link
The Russian version is on arxiv.org
Here the formula is given in terms of comultiplication and the proof is based on the classical definition of higher covariant derivative (i.e. by iterating a partial covariant derivative formula, cf. e.g. Kobayashi-Nomizu Vol. I §III.2). This makes for a more complicated proof, but exhibits very interesting algebra (related to the algebra of symbols of differential operators).
Gavrilov wrote in an email that he was also wondering if the formula is known: 

I had a hope that someone will write me about it, but you are the first.

Thus I regard my question still unanswered...

Answer 2:

Here comes the statement and proof in my private formalism of "dual total tensor calculus":  Roughly it is the physicist's abstract index notation with abstract indices replaced by uncompromising abstract multilinear algebra. It might look clumsy at first, but it can vastly simplify some extremely complicated classical computations.
In my formalism the Leibniz formula suggests itself from the "combinatorial" definition of the shuffle product. It would also allow for a quick proof based on the Lie-Hopf unshuffle coproduct. 

$\newcommand{\shuffle}{\mathrm{ш}} \newcommand{\id}{\mathrm{id}} \newcommand{\von}{ \mathrm{:}\ }
\newcommand{\nach}{\mathbin{\rightarrow}} \newcommand{\tens}{\mathbin{\otimes}}
\newcommand{\boxtens}{\mathbin{\boxtimes}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\Om}{\Omega}
\newcommand{\C}{\mathcal{C}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\homnabla}{{ }^\mathrm{J}\!\nabla}
\DeclareMathOperator{\Hom}{Hom}
\newcommand{\diff}{\mathrm{d}}$
Algebra Notation
Permutation operators. Let $E_1,\ldots, E_n$ be vector spaces/bundles and $\sigma\von\{1,\ldots,n\}\nach\{1,\ldots,n\}$ a permutation. Denote by
$$\bigl(\begin{smallmatrix}E_1&E_2&\cdots&E_n\\\sigma(1)&\sigma(2)&\cdots&\sigma(n)\end{smallmatrix}\bigr)\von
E_1\tens\cdots\tens E_n\longrightarrow  E_{\sigma^{-1}(1)}\tens\cdots\tens E_{\sigma^{-1}(n)}$$
the isomorphism which puts $E_i$ to the new place $\sigma(i)$.
For the most frequently occuring permutation operator we introduce an extra sign:
Let $E,F,A,B$ be vector spaces/bundles. Set
$$
(E\tens A)\boxtens_{E,F} (F\tens B)  :=  E\tens F\tens A\tens B  = 
(\begin{smallmatrix}E&A&F&B\\1&3&2&4\end{smallmatrix})\cdot E\tens A\tens F\tens B
$$
Tensor algebra. 
We write $
V^\otimes_\bullet := \bigoplus_{p\in\N_0} V^\otimes_p
$
where $V^\otimes_0=\R$ and $V^\otimes_p:=V^{\tens p}$. This helps keeping lines a bit shorter.
The shuffle product
on the tensor algebra $V^\otimes_\bullet$
is the direct sum of the recursively defined maps
$$\shuffle_{p,q}\von V^\otimes_p \tens V^\otimes_q \nach V^\otimes_{p+q}$$
where $\shuffle_{0,0} := 1$, $\shuffle_{p,0} := \id_{V^\otimes_p}$, $\shuffle_{0,q} := \id_{V^\otimes_q}$ and
for $p,q\ge1$,
$$
\shuffle_{p,q}  :=  \shuffle_{p-1,q} \tens \id_V 
         \cdot \left(\begin{smallmatrix}V^\otimes_{p-1}&V&V^\otimes_{q}\\1&3&2\end{smallmatrix}\right)
                          +  \shuffle_{p,q-1} \tens \id_V    %\\
$$
Written multiplicatively, for $\Phi\in V^\otimes_r$, $\Psi\in V^\otimes_s$, $v,w\in V$,
$$
(\Phi\tens v) \shuffle (\Psi\tens w)  :=  (\Phi\shuffle (\Psi\tens w))\tens v  +  ((\Phi\tens v)\shuffle \Psi)\tens w 
$$

Calculus notation
Let $M$ be a differential manifold. Abbreviate $\F:=\C^\infty(M)$ and $\Om:=T^\ast M$ for the cotangent bundle. (It is the primordial bundle in my calculus, like with Zariski's construction for locally ringed space.)
The total covariant derivatives on a tensor bundle $E\nach M$ is a $\R$-linear map of sections
$$\nabla\von \Gamma(E)\nach \Gamma(E\tens\Om)$$ such that
$$\nabla(\varphi e)  =  \varphi\nabla e + e\tens\diff\varphi  \quad \forall  \varphi \in\F, e\in\Gamma(E)$$
For being serious with abstract multilinear algebra a separate notion of total covariant derivative on homomorphism bundles is eventually needed:
$$\homnabla\von \Gamma(\Hom(E,F))\nach \Gamma(\Hom(E,F\tens\Om))$$
(It will be used here only implicitly.)
Total c.d. on bundles $E\nach M$ and $F\nach M$ are extended to (or assumed compatible with) tensor products and homomorphisms by the
total product rules
$$\nabla( e\tens f ) =  \bigl(\nabla e\bigr)\boxtens_{E,F} f +  e\tens\bigl(\nabla  f\bigr)
  =  (\begin{smallmatrix}E&\Om&F\\1&3&2\end{smallmatrix})\cdot\bigl(\nabla e\bigr)\tens f
         +  e\tens\bigl(\nabla f\bigr)$$
$$\nabla(h\cdot e)  =  \bigl(\homnabla h\bigr)\cdot e  +  h\tens\id_\Om\cdot\nabla e$$
where $e\in\Gamma(E)$, $f\in\Gamma(F)$, $h\in\Gamma(\Hom(E,F))$.
Example: From the total product rules it is immediately seen that $\homnabla$ applied to permutation operators gives zero. Thus we get
$$\nabla(a\boxtens_{E,F}b)
 =  \id_{E\tens F}\tens(\begin{smallmatrix}A&\Om&B\\1&3&2\end{smallmatrix})
 \cdot (\nabla a)\boxtens_{E,F}b  +  a\tens\nabla b
$$
This might look a tad ridiculous but exhibits the core of the following proof of the Leibniz formula.
Now fix some covariant derivative on $\Om$, not necessarily torsion free.
The n-th order covariant derivative on $E$, 
$$
\nabla^n \von \Gamma(E) \nach \Gamma(E\tens\Om^\otimes_n) \qquad(n\in\N_0)
$$
is defined as $\nabla^{0}:=\id_E$ and $\nabla^{n+1} := \nabla^{E\tens\Om^\otimes_n}\nabla^n$ using the total product rule.

Theorem. The Leibniz formula holds for the n-th order covariant derivative on tensor products:
$$\nabla^{n} (e\tens f)  = 
  \sum_{p+q =n} \id_{E\tens F}\tens\shuffle_{p,q}\cdot \nabla^{p}e\boxtens_{E,F}\nabla^{q}f$$
where $\shuffle$ is the shuffle multiplication on $\Om^\otimes_\bullet$.
Proof
The case $n=0$ is trivial. Assume the formula holds for some $n$.
Then
$$
\nabla\nabla^n(e\tens f)   =  
\sum_{p+q =n} 
   \id_{E\tens F}\tens\shuffle_{p,q}\tens\id_\Om \cdot
   \nabla\left(\nabla^{p}e\boxtens_{E,F}\nabla^{q}f\right)
$$
since 
$\homnabla(\id_{E\tens F}\tens\shuffle_{p,q})  = 0$ since this is a sum of permutation operators. Moreover we have
$$
\nabla\left(\nabla^{p}e\boxtens_{E,F}\nabla^{q}f\right)    =  
 \id_{E\tens F}\tens \left(\begin{smallmatrix}\Om^\otimes_p&\Om&\Om^\otimes_{q}\\1&3&2\end{smallmatrix}\right)
    \cdot\nabla^{p+1}e\boxtens_{E,F}\nabla^{q}f
    +  
 \nabla^{p}e\boxtens_{E,F}\nabla^{q+1}f
$$
Reorganizing the summation now gives
$$ \nabla^{n+1} (e\tens f)  =  \\
  \id_{E\tens F}\tens\shuffle_{n,0}\tens\id_\Om \cdot \nabla^{n+1}e\boxtens_{E,F}\nabla^{0}f\\
     +  \id_{E\tens F}\tens\shuffle_{0,n}\tens\id_\Om \cdot \nabla^{0}e\boxtens_{E,F}\nabla^{n+1}f\\
     +   \sum_{\substack{p+q=n+1\\1\le p\le n}}
        \id_{E\tens F}\tens
        \left(
          \shuffle_{p-1,q}\tens\id_\Om\cdot 
      \left(\begin{smallmatrix}\Om^\otimes_{p-1}&\Om&\Om^\otimes_q\\1&3&2\end{smallmatrix}\right)
           + 
          \shuffle_{p,q-1}\tens\id_\Om
        \right)  \cdot
        \nabla^{p}e\boxtens_{E,F}\nabla^{q}f 
$$
By definition of $\shuffle_{p,q}$ this gives the formula for $n+1$.  Q.E.D.

Remarks 1 $\shuffle_{1,1}$ gives symmetrization. Thus the case of $n=2$ and of one bundle trivial proves that curvature is a tensor.
2 A similar formula is proved the same way for homomorphisms. This can be used to prove Gavrilov's "cocycle identity" from the cited paper.
 A: I have not found this exact formula in the physical or mathematical contexts where I would expect to find something like this. (That is, it's not found in any of the several supermanifolds monographs, where such generalizations are most necessary.)
(BTW thanks for the reference to Gavrilov's paper, some interesting things there. Especially in relating this to modern questions.)
In abstract index calculus there is an extended discussion of the general Leibniz law for tensor calculus and related questions in Penrose, R & Rindler, W. 1984. Spinors and spacetime, 1. C:CUP, pp. 191-194. It does not contain this way of defining it, but contains a Leibniz law satisfying generalization.
There is an ultimately equivalent operation (to Gavrilov's formula) that can repeated inside itself thanks to the diagrammatic notation for equivalent generalization: in Penrose, R. The road to reality. London, Cape. P. 300, Fig. 14.6. (This notation is how I think of it, btw, introduced at the end in Penrose and Rindler's monograph.)
Altho the operations expand out to the same objects in the same order, the exact formula with shuffle product does not seem to be present, but an equivalent generalization.
UPDATE:
Further looking in my library. I did not find much discussion of this, despite expecting to find it there, in the more recent comprehensive physical monographs, such as DeWitt, B. 2003. The global approach in quantum field theory, 1. Several good comprehensive texts are missing the general discussion mostly (for example Choquet-Bruhat, Y & DeWitt-Morette, C & Dillard-Bleick, M. 1996. Analysis, manfolds, and physics, 1.)
However the abstract Noll, W. 1987. Finite dimensional spaces, 1. Dordrecht, Nijhoff, chapter 6, section 66, has an approach that is generally relevant here.
A: With sufficient probability the answer is "no".
