The combinatorics of the divergence term in the calculus of variations This question is regarding higher-order, multiple variation problems, but as far as I suspect it is basically a combinatorics problem in disguise that I can't get sorted out by myself.
In fact, the question will be formulated without the calculus-of-variations context. I want to look at operators of the form $$ P[\phi]=\sum_{k=0}^{r}P^{i_1...i_k}\partial_{i_1...i_k}\phi, $$ where $\partial_{i_1...i_k}=\partial_{i_1}\dots\partial_{i_k}$ are repeated partial derivatives with respect to some $m$ independent variables $x^1,\dots ,x^m$ and $\phi$ and the $P^{i_1...i_k}$ are smooth functions of the $x^i$. It will be helpful to use multiindex notation where $I=(i_1\dots i_k)$ is a multiindex of length $k$ ($|I|=k$) and is just an ordered list of $k$ indices. Einstein summation convention is assumed on all repeated indices (but not the multiindices). In terms of the multiindices we have $$ P[\phi]=\sum_{|I|=0}^rP^I \partial_I\phi. $$ My multiindex notation agrees with that of Anderson and Krupka ({1} and {2}), but not with the more combinatorial multiindex notation that is often used in analysis.
As it is well-known (eg. {1}) it is possible to rewrite this operator as $$ P[\phi]=\sum_{|I|=0}^r\partial_I(Q^I\phi)=Q\phi+\partial_i(Q^i\phi)+\partial_{ij}(Q^{ij}\phi)+\dots, $$ where $$ Q^I=\sum_{|J|=0}^{r-|I|}C^{|I|+|J|}_{|I|}(-1)^{|J|}\partial_J P^{IJ}, $$and where $C^k_l=\left(\begin{matrix}k\\l\end{matrix}\right)$ is the binomial coefficient.
I want to write $P[\phi]$ into yet another form, namely $$ P[\phi]=Q\phi+\partial_i\left(\sum_{|J|=0}^{r-1}F^{iJ}\partial_J\phi\right). $$
The problem is that I know that $F^{iJ}$ can be expressed as $$ F^{iJ}=\sum_{|K|=0}^{r-1-|J|}(-1)^{|K|}\partial_K P^{iJK}, \quad(\ast)$$with no combinatorial factors whatsoever because for example this formula is present in eg. {2},{3} and some other texts on the calculus of variations that I can't recall right now, but {3} gives no proof and the proof of {2} is not applicable here ({2} derives the so-called principal Lepage equivalent of a Lagrangian, which involves the expressions $F^{iJ}$, but the proof simply does not work for Lagrangians without introducing Lepage forms).
Moreover, by looking at this decomposition for low orders ($r=2,3$), I can see explicitly that there are no combinatorial factors.
However if I try to use the formula for $Q^I$, I get a very complicated combinatorial expression that does not seem to reduce.
Specifically, we have $$ P[\phi]=\sum_{|I|=0}^{r}\partial_I(Q^I\phi)=Q\phi+\sum_{|I|=0}^{r-1}\partial_i\partial_I(\phi Q^{iI}), $$ and here we can use the higher order Leibniz rule $$ \sum_{|I|=k}\partial_I(\phi Q^{iI})=\sum_{|J|+|K|=|I|}C^{|J|+|K|}_{|K|}\partial_J\phi \partial_K Q^{iJK} $$ and thus $$ \sum_{|I|=0}^{r-1}\partial_I(\phi Q^{iI})=\sum_{|J|+|K|=0}^{r-1}C^{|J|+|K|}_{|K|}\partial_J\phi \partial_K Q^{iJK}. $$ Inserting back the relation for $Q$ gives $$ \sum_{|I|=0}^{r-1}\partial_I(Q^{iI}\phi)=\sum_{|J|+|K|=0}^{r-1}\sum_{|L|=0}^{r-1-(|J|+|K|)}C^{|J|+|K|}_{|K|}C^{1+|J|+|K|+|L|}_{|L|}(-1)^{|L|}\partial_{KL}P^{iJKL}\partial_J\phi. $$
The sum basically goes as $|J|+|K|+|L|=0,1,\dots,r-1$, thus we can write $$ F^{iJ}=\sum_{|K|+|L|=0}^{r-1-|J|}C^{|J|+|K|}_{|K|}C^{1+|J|+|K|+|L|}_{|L|}(-1)^{|L|}\partial_{KL}P^{iJKL}.\quad(\ast\ast)$$
My problem is that comparing this formula with $(\ast)$ I don't see how the combinatorial terms reduce.
It is clear that at this point I should introduce a new multiindex $M=KL$ and rewrite the sums and coefficients in a way that only $M$ and $J$ remain, but for example the binomial terms give $$C^{|J|+|K|}_{|K|}C^{1+|J|+|K|+|L|}_{|L|}=\frac{(|J|+|K|+|L|)!}{|J|!|K|!|L|!}\frac{1+|J|+|K|+|L|}{1+|J|+|K|}, $$and it is really not at all obvious that this "change of variables" is possible. Not to mention the fact that in $(\ast\ast)$ $(-1)^{|L|}$ occurs while a comparision with $(\ast)$ shows that we must have $(-1)^{|K|+|L|}$. But I can't see where the additional factor of $(-1)^{|K|}$ could come from.
Question: How to show that $(\ast)$ and $(\ast\ast)$ are equivalent formulae? If this is not a feasible path, then how to prove explicitly that $$P[\phi]=Q\phi+\partial_i\left(\sum_{|J|=0}^{r-1}F^{iJ}\partial_J\phi\right),$$ where $F^{iJ}$ is given by $(\ast)$?
Remark: Although it should be clear, this is related to the calculus of variations as follows. When one is given a Lagrangian $\mathcal L(x,u,u_{(1)},\dots,u_{(r)})$ of order $r$ (where $u=u(x)$ is the dynamical variable and $u_{(k)}=u_{i_1...i_k}$ is the $k$th derivative of $u$), and one varies $\mathcal L$ with $\delta u=\phi$, then we have $$ P^{i_1...i_k}=\frac{\partial\mathcal L}{\partial u_{i_1...i_k}} $$ and $$ Q=E(L)=\sum_{|I|=0}^r(-1)^{|I|}\partial_{I}\frac{\partial\mathcal L}{\partial u_I}. $$
References:

*

*{1} I. M. Anderson, The Variational Bicomplex

*{2} D. Krupka, Introduction to Global Variational Geometry

*{3} G. W. Bluman, S. Kumei, Symmetries and Differential Equations
 A: I am not certain this gives a direct answer to your question, but from looking over what you have given. This looks to me like you're asking about what my reference calls the higher Euler operators. This is discussed on pg. 365-367 of the reference, which is Applications of Lie Groups to Differential Equations by Peter J. Olver.
Following the reference, he first gives an example to show the origins of the general pattern. Then states the first few operators for the general case before introducing some further multi-index notation to state the full general case.
In particular he gives
$$\mathsf{E^{(0)}}(P)=\frac{\partial P}{\partial u}-D_x\frac{\partial P}{\partial u_x}+D^2_x\frac{\partial P}{\partial u_{xx}}-D^3_x\frac{\partial P}{\partial u_{xxx}}+\cdots$$
Which is the usual Euler operator (for 1 independent variable), and then.
$$\mathsf{E^{(1)}}(P)=\frac{\partial P}{\partial u_x}-2D_x\frac{\partial P}{\partial u_{xx}}+3D^2_x\frac{\partial P}{\partial u_{xxx}}-4D^3_x\frac{\partial P}{\partial u_{xxxx}}+\cdots$$
$$\mathsf{E^{(2)}}(P)=\frac{\partial P}{\partial u_{xx}}-3D_x\frac{\partial P}{\partial u_{xxx}}+6D^2_x\frac{\partial P}{\partial u_{xxxx}}-10D^3_x\frac{\partial P}{\partial u_{xxxxx}}+\cdots$$
Still working w/ 1 independent variable.
Now to proceed to the general case he introduces more notation. Earlier in the text, (Section 2.3, pg. 95) He introduced the unordered multi-index notation, $J=(j_1,j_2,...j_k)$ allowing him to write things like
$$\partial_J f(x)=\frac{\partial f}{\partial x^{j_1}\partial x^{j_2}\cdots \partial x^{j_k}}$$
Now let $I,J$ be indices of the type just described. We will say that $J\subset I$ if every index in $J$ appears in $I$ and write $I\setminus J$ for the remaining indices. For example, if $J=(1,1,2,4)$ and $I=(1,1,1,2,4,4)$ and $J\subset I$ and $I\setminus J=(1,4)$. Given $I=(i_1,\ldots ,i_n)$, let $\tilde{I}=(\tilde{i}_1,\ldots,\tilde{i}_p)$ be the "transposed" ordered multi-index, where $\tilde{i}_j$ equals the number of occurrences of the integer $j$ in $I$. For the above example, $\tilde{I}=(3,1,0,2)$ since there are three $1$'s, one $2$, no $3$'s, and two $4$'s in $I$. Set $I!=\tilde{I}!=\tilde{i}_1!\tilde{i}_2!\cdots \tilde{i}_p!$ (note: $p$ is the number of independent variables), and define the multinomial coefficient
$$\left(\begin{array}{c}
I\\
J
\end{array}\right)=I!/(j!/(I\setminus J)!$$
when $J\subset I$, $0$ otherwise. In the above example, we have $I!=3!\cdot 1! \cdot 0! \cdot 2!=12$ and
$$\left(\begin{array}{c}
I\\
J
\end{array}\right)=12/(2\cdot 1)=6$$
Armed with this notation, he then states the general form of the higer Euler operators in Prop 5.96.
Let $1\leq\alpha\leq q$ (note: $q$ dependent variables), #$J\geq0$ (#$J$ is just the number of indices in $J$). Then
$$\mathsf{E}^J_{\alpha}(P)=\sum_{I\supset J}\left(\begin{array}{c}
I\\
J
\end{array}\right)(-D)_{I\setminus J}\frac{\partial P}{\partial u^{\alpha}_I}$$
for all $P\in \mathscr{A}$ (Note: $\mathscr{A}$ is the algebra of differentiable functions).
I think the conbinatorial thing you're noticing is this business of the multinomial coefficient discuss in the above excerpt from my text. Not sure if you're only interested in the one independent variable case or the general case, but either way, hope that helps you figure out the pattern.
A: In the meantime I managed to solve this problem but I still don't understand how the combinatorial coefficients simplify.
It turns out Krupka's derivation of the coefficients of the principal Lepage extension is actually workable in the "classical" formulation and I was too hasty in dismissing it.
We basicall want to write $$ P[\phi]=\sum_{k=0}^rP^{i_1...i_k}\partial_{i_1...i_k}\phi=F\phi+\partial_i\left(\sum_{k=0}^{r-1}F^{ii_1...i_k}\partial_{i_1...i_k}\phi\right). $$ Calculating the derivative of the parentheses gives $$ P[\phi]=F\phi+\sum_{k=0}^{r-1}\partial_i F^{ii_1...i_k}\partial_{i_1...i_k}\phi+\sum_{k=1}^{r}F^{i_1...i_k}\partial_{i_1...i_k}\phi. $$ The coefficients $P^{i_1...i_k}$ are unique as long as they are symmetric in the indices, so we get the recursion relations $$ F^{i_1...i_r}=P^{i_1...i_r} \\ F^{i_1...i_k}=P^{i_1...i_k}-\partial_i F^{ii_1...i_k},\quad 0\le k\le r-1. $$ We can now proceed by induction. We suspect the formula $$ F^{i_1...i_k}=\sum_{l=0}^{r-1-k}(-1)^l\partial_{j_1...j_l}P^{j_1...j_lii_1...i_k}, $$ which is true for $k=r$, then we can use descending induction on $k$.
I still would like to see however how the combinatorial expressions in $(\ast\ast)$ reduce in my OP.
