I would like to expand the following expression
$$\left(\sum_{i,j=1}^N \,x_i A_{ij} x_j\right)^n$$
where $\mathbf A$ is a symmetric $N\times N$ matrix, $\mathbf {x}$ is an $N$-component vector, and $n$ is a non-negative integer power. The expansion of this expression yields a homogeneous polynomial of order $2n$ in the $x_k$.
What is the coefficient of the term $x_1^{p_1} x_2^{p_2} \cdots x_N^{p_N}$ for $p_1 + p_2 + \cdots + p_N = 2n$ in the expansion of this expression?
Has the formula been worked out before?
The expression can be succintly written as ${\left({\mathbf{x}}^TA\mathbf{x} \right)}^n$, where $\mathbf{x}=(x_1,x_2,\ldots,x_N)$. The inner product $x^TAx$ is some polynomial $p_A(\mathbf{x})$ that is a sum of monomials of degree exactly $2$. Hence, each nonzero monomial in ${\left({\mathbf{x}}^TA\mathbf{x}\right)}^n ={\big[p_A(\mathbf{x})\big]}^n$ is of the form $x_1^{p_1}x_2^{p_2}\dots x_n^{p_N}$ with $p_i\geq 0$ and $\sum_ip_i=2n$.
Let $\rho$ be one such $N$-tuple $\rho=(p_1,p_2,\ldots,p_N)$, and let $\lVert\rho\rVert = \sum_ip_i$. We will use combinatorics to calculate the coefficient of ${\mathbf{x}}^{\rho} = x_1^{p_1}x_2^{p_2}\dots x_n^{p_N}$ in ${\left({\mathbf{x}}^TA\mathbf{x}\right)}^n$. We will consider each different way in which we can choose a term from each copy of $p_A$ in the (ordered) product ${\big[p_A(\mathbf{x})\big]}^n$ so that the end result is a ${\mathbf{x}}^{\rho}$ monomial.
Let $\mho$ be an alphabet $\mho=\{l_1,l_2,\dots,l_N\}$. We associate to each (ordered) choice of monomial $x_iA_{ij}x_j$ in $p_A(\mathbf{x})$ the string $l_il_j$ in $\mho$. This induces a correspondence between choices in the ordered product ${\big[p_A(\mathbf{x})\big]}^n$ and strings of length $\lVert\rho\rVert$ in $\mho$.
Notice that two different strings may give rise to the same coefficient. For instance, $l_2l_2l_4l_1l_1l_3$ is also associated to the coefficient $A_{22}A_{41}A_{13}=A_{13}A_{41}A_{22}$, but this is because the product of the numbers $A_{ij}$ does not depend on the order of the terms. This is not a contradiction with our correspondence, because different strings correspond to different choices (even if they are associated to the same coefficient).
Let $A(s)$ be the coefficient associated to a string $s$. We note that $A(s)$ is in the coefficient of ${\mathbf{x}}^{\rho}$ if and only if $s$ contains $p_i$ $l_i$'s for each $i=1,2\ldots,N$. Hence, with
\begin{equation} S(\rho)=\left\{s\,\middle|\, \begin{array}{l} \text{$s$ is a string of length $\lVert \rho \rVert$ on $\mho$ and}\\ \text{for each $i=1,\dots,N$ the letter $i$ appears $p_i$ times on $s$} \end{array}\right\},\end{equation}
then the coefficient $K(\rho)$ of ${\mathbf{x}}^{\rho}$ is given by
\begin{equation}\tag{1}\label{ksr} K(\rho) = \sum_{s \in S(\rho)}A(s) \end{equation}
We pause briefly to note a rather obvious lemma, that will nonetheless be useful to us.
Proof. In ${\left(\sum_{i,j=1}^N \,x_i A_{ij} x_j\right)}^n$, the only choice that yields $x_k^{2n}$ is taking $i=j=k$ in each term of the $n$ products. $\square$
In the statement of the theorem below, we highlight the following piece of notation. The hat above $c_j$ in $c_1+\ldots+\hat{c_j}+\ldots+c_N$ means $c_j$ does not feature in the sum, and $e_i$ is the vector
$$e_i=\underbrace{(0,0,\ldots,0,1,0,\ldots,0,0)}_{\text{$1$ at position $i$}}$$
With that out of the way:
Notice that the expansion calculates $K(p_1,\ldots,p_N)$ in terms of coefficients $K(\overset{\sim}{p_1},\ldots,\overset{\sim}{p_N})$, where each $\overset{\sim}{p_i}\leq p_i$ and, most importantly, $\overset{\sim}{p_j}=0$. Together with the lemma, this means $K(p_1,\ldots,p_N)$ can be calculated in at most $(N-1)$ applications of the theorem.
Proof. We consider $A$ to be a 'free' symmetric matrix, meaning we assume no relations between the $A_{ij}$ beyond $A_{ij}=A_{ji}$. That said, the theorem can be adapted to general 'free' matrices.
The idea behind the expansion is to rewrite $\eqref{ksr}$ in a sum of the type
\begin{equation}\tag{2}\label{krp} K(\rho)=\sum_{\alpha}\big|\{s \in S(\rho)\,|\,A(s)=\alpha\}\big|\cdot \alpha \end{equation}
where $|X|$ is the cardinality of set $X$. Now, $\big|\{s \in S(\rho)\,|\,A(s)=\alpha\}\big|$ is given by a multinomial coefficient. If $\alpha=\prod_{1\leq i<k\leq N}{\left(A_{ik}\right)}^{n_{ik}}$, and remember $\sum n_{ik}=n$ as per hypotheses, then
$$\big|\{s \in S(\rho)\,|\,A(s)=\alpha\}\big|=\frac{n!}{\prod_{1\leq i<k\leq N}n_{ik}!}$$
Next, we consider how terms $A_{ji}$ (where $j$ is the fixed index choice of the thoerem) show up in the values $\alpha$. Afterewards, the substrings that remain have no $l_j$, so the problem is reduced to an easier one, with a smaller alphabet.
Each $A_{jj}$ corresponds to two $l_j$'s in a string $s$, while each $A_{ji}$ with $j\neq i$ corresponds to a single $l_j$. Hence, if $a$ is the number of $A_{jj}$'s in $A(s)$ and $b$ is the number of $A_{ji}$'s with $j\neq i$, we have that $2a+b=p_j$ and that the total number of terms in $A(s)$ which feature an index $j$ is
$$a+b=p_j-a.$$
The number of terms in $A(s)$ that don't feature an index $j$ is thus $n-(p_j-a)$. Notice that $0\leq a \leq \left\lfloor p_j/2\right\rfloor$.
Now, a choice of $a$ 'well-defines' the terms $A_{jj}$'s, but, even though it defines $b$, there's still room for choice with the $A_{ji}$'s (namely, the indices $i$). To know the terms $A_{ji}$, it suffices to choose a submultiset $M$ from
$$N=\bigcup\limits_{\substack{1\leq k\leq N\\k\neq j}}\left\{l_k^{p_k}\right\}$$
with $|M|=b=p_j-2a$, where above we're using the multiset exponential notation for multiplicities. $M$ represents the indices $i$ (with their multiplicities).
Hence, each such submultiset has the form $M=\bigcup\limits_{\substack{1\leq k\leq N\\k\neq j}}\left\{l_k^{c_k}\right\}$ and corresponds to a solution of
\begin{equation}\tag{3}\label{msys} \left\{ \begin{array}{ll} c_1+\ldots+\hat{c_j}+\ldots+c_N=p_j-2a&\\ 0\leq c_i \leq p_i & (\forall i\in\{1,\ldots,N\}\setminus\{j\}) \end{array}\right. \end{equation}
Besides not having any $l_j$'s, the substring that remains will have $c_i$ fewer $l_i$'s for each $i\in\{1,\ldots,N\}\setminus\{j\}$.
This answer is getting long enough as is, but I'll try to summarize the points made and hope it makes for a satisfactory proof:
If you have any questions about this proof I'll be happy to answer, but I hope this is clear enough. $\square$
Notice the subscript in the summation of $F(\rho,j,a)$. In this form, the number of terms in the summation of $F(\rho,j,a)$ can be found easily via stars and bars. The number is
$$\binom{(p_j-2a)+\left(\left|I^+\right|-1\right)-1}{p_j-2a}= \binom{p_j+\left|I^+\right|-2a-2}{p_j-2a}$$
Proof. For $i \notin I^+$, $p_i=0$. Hence, $0\leq c_i\leq p_i$ implies $c_i=0$, and it's easy to check that there's no loss in removing them from the summation. On the other hand, $\sum_{i\in I^+\setminus{\{j\}}}c_i=p_j-2a$ and $c_i\geq 0$ imply $c_i \leq p_j-2a\leq p_j\leq p_i$, so the sum is as in the theorem. $\square$