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.