I've been trying to find the expansion to the expression given below for some time now.
$$\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^Na[i]a[j]\right)^2$$
So what I am looking for is something like this: $$\left(\sum_{i=1}^Na[i]\right)^2=\sum_{i=1}^Na[i]^2+2\sum_{i=1}^{N-1}\sum_{j=i+1}^Na[i]a[j]$$
Any help will be appreciated.
Thanks
As Eric Synder said in the comments, one can obtain the following formula simply by keeping careful track of the terms: $$ \left(\sum_{i<j} a[i]a[j]\right)^2 = \sum_{i<j} a[i]^2a[j]^2 + 2\sum_{i<j}\sum_{k\neq i, j} a[k]^2a[i]a[j] + 6\sum_{i<j<k<\ell} a[i]a[j]a[k]a[\ell] $$ (For readability, I've omitted the assumption in all sums that the variables are integers between $1$ and $N$. Note that this doesn't contradict, e.g. $i\leq N-1$ on the left-hand side, since $j\leq N$ as well, and $i<j$.)
Note: If $a$ has fewer than four entries, then obviously it is impossible to find four distinct $i,j,k,\ell$ to use for the last summation. This is fine, though; the formula still holds, just with the last summation is being empty.
Theory
However, for more general problems there is a known device for making this bookkeeping easier. This is usually stated in the language of symmetric functions, and so I will give a very brief introduction here. Symmetric functions are polynomial-like expressions in infinitely many variables $a[i]$ (usually $x_i$, but I will keep using your notation), such that any permutation of (finitely many of) the variables does not change the function.
(Aside: It simplifies the theory somewhat to use infinitely many variables; if we only use $N$, simply eliminate anything arising from a tuple with more than $N$ entries. However, based on your notation I assume you are using this for a coding application, in which case your array $a$ is likely to be long enough to not need to worry about this; see e.g. the note above.)
The most "obvious" symmetric functions are the monomial symmetric functions $m_\lambda$, e.g. $$m_{3,1,1} = a[1]^3a[2]a[3] + a[1]^3a[2]a[4]+ \cdots + a[2]^3a[1]a[3] + a[2]^3a[1]a[4] + \cdots +\cdots\cdots.$$ I won't attempt to write down a general definition but it's what you think: start with the mononial term $a[1]^{\lambda_1}a[2]^{\lambda_2}\cdots a[k]^{\lambda_k}$ and then exchange the indices $1,2,\dots, k$ for any other set of natural numbers, throwing out any repeated terms.
If one insists that $\lambda$ are always partitions, that is, a tuple $(\lambda_1\geq \lambda_2\geq \cdots\geq\lambda_k)$, then every function $m_{\lambda}$ is distinct, and every symmetric function is a unique linear combination of the $m_{\lambda}$. (That is, they form a basis for the space of symmetric functions)
Less obvious but quite useful are the elementary symmetric functions $e_\lambda$, which are built recursively. When $\lambda=(n)$ is a 1-tuple, then $$e_n = \sum_{i_1< \cdots < i_k} a[i_1]a[i_2]\cdots a[i_k],$$ and in general $e_\lambda = e_{\lambda_1}e_{\lambda_2}\cdots e_{\lambda_k}$. Again, if we insist that $\lambda$ are always partitions, then every function $e_\lambda$ is distinct and they form a basis for the space of symmetric functions.
Generalization
Your toy example $$\left(\sum_{i=1}^Na[i]\right)^2=\sum_{i=1}^Na[i]^2+2\sum_{i=1}^{N-1}\sum_{j=i+1}^Na[i]a[j]$$ is precisely the statement that $e_{1,1} = m_2 + 2m_{1,1}$, and the expression you were asking about, $\left(\sum_{i=1}^{N-1}\sum_{j=i+1}^Na[i]a[j]\right)^2$, is precisely $e_{2,2}$. Therefore, I interpreted your question as asking for how to write $e_{2,2}$ in the monomial basis, and the formula I gave earlier corresponds to the fact that $e_{2,2} = m_{2,2}+2m_{2,1,1}+6m_{1,1,1,1}$.
One may ask the more general question of how to write an arbitrary $e_\lambda$ in the monomial basis, and the answer is as follows:
$$ e_\lambda = \sum_{A \in \mathcal{P}^\text{min}_{0,1}(\lambda)} m_{(\text{# of 1s in each row of }A)}$$
Here the sum is over the set $\mathcal{P}^\text{min}_{0,1}(\lambda)$, which is the set of matrices such that
This is a somewhat intricate set of conditions, of course, but most people find this easier than keeping track of all the terms! For instance, when $\lambda=(2,2)$ the set $\mathcal{P}^\text{min}_{0,1}(\lambda)$ consists of the following 9 matrices:
$$ \begin{bmatrix} 1&1\\1&1 \end{bmatrix} \begin{bmatrix} 1&1\\1&0\\0&1 \end{bmatrix} \begin{bmatrix} 1&1\\0&1\\1&0 \end{bmatrix} $$ $$ \begin{bmatrix} 1&0\\1&0\\0&1\\0&1 \end{bmatrix} \begin{bmatrix} 1&0\\0&1\\1&0\\0&1 \end{bmatrix} \begin{bmatrix} 1&0\\0&1\\0&1\\1&0 \end{bmatrix} \begin{bmatrix} 0&1\\1&0\\1&0\\0&1 \end{bmatrix} \begin{bmatrix} 0&1\\1&0\\0&1\\1&0 \end{bmatrix} \begin{bmatrix} 0&1\\0&1\\1&0\\1&0 \end{bmatrix} $$ (You can probably imagine how to generate these systematically, although I certainly wouldn't want to write a careful algorithm.)
As you can see the three in the first row correspond to $m_{2,2}$ and the two copies of $m_{2,1,1}$, and the six in the last row all correspond to copies of $m_{1,1,1,1}$, hence the formula at the beginning of the answer.
Thanks to Mathematical Gemstones for this convenient reference. It is the conclusion of a longer exposition on symmetric functions, which beginners may find useful.