Sum of Product of Normal Random Variables Consider a collection of $n$ i.i.d.normal random variables $X_i \sim \mathcal{N}(\mu, \sigma^2)$, $i=1,\ldots,n$.  I'm trying to compute the distribution of
$$
\sum_{i=1}^n \sum_{j=i+1}^n X_i X_j.
$$
It seems that each $X_i X_j$ has the normal product distribution, but the sum is elusive for me.

Update based on the answer from @Robert Israel:

Transform $\bf Z$ by an orthogonal matrix...

Let's say ${\bf Z} \in \mathbb{R}^2$ and let the orthogonal matrix be
$$
A = 
\begin{pmatrix}
\frac{1}{\sqrt 2} & -\frac{1}{\sqrt 2} \\
\frac{1}{\sqrt 2} & \frac{1}{\sqrt 2}
\end{pmatrix}.
$$
Applying this transformation gives
$$
A {\bf Z} = 
\begin{pmatrix}
\frac{1}{\sqrt 2}Z_1 -\frac{1}{\sqrt 2}Z_2 \\
\frac{1}{\sqrt 2}Z_1 +\frac{1}{\sqrt 2}Z_2 \\
\end{pmatrix}.
$$
Since we have a ${\bf e}^T {\bf Z}$ term, apply this to the transformed ${\bf e}$:
$$
\begin{pmatrix}
1 & 1
\end{pmatrix}
\begin{pmatrix}
\frac{1}{\sqrt 2}Z_1 -\frac{1}{\sqrt 2}Z_2 \\
\frac{1}{\sqrt 2}Z_1 +\frac{1}{\sqrt 2}Z_2 \\
\end{pmatrix}
=
\sqrt 2 Z_1.
$$
I believe @Robert Israel wrote this as $\sqrt n W = \sqrt n Z_1$, and mentioned

Thus $W$ is the component of ${\bf Z}$ in the direction of the unit vector $n^{−1/2}{\bf e}$, while the sum of the squares of the components for the other $n−1$ basis vectors in an orthonormal basis is $U$.

First, it seems that from my example, $W$ is just the first component of ${\bf Z}$, namely $Z_1$, so I'm not sure how this is also the first component of ${\bf Z}$ in the direction of $n^{−1/2}{\bf e}$.  Secondly, I still don't quite understand the sum of squares of the other basis components becoming $U$.
 A: Your sum is $$S_n = \sum_{i=1}^n \sum_{j=i+1}^n X_i X_j = \dfrac{1}{2}{\bf  X}^T C \bf X$$
where $C$ is the $n \times n$ matrix with all off-diagonal entries $1$ and diagonal entries $0$.  We can write $C = e e^T - I$ where $e = [1,\ldots,1]^T$, and 
 $\bf X = \mu \bf e + \sigma \bf Z$ where $\bf Z$ is multivariate normal with mean $0$ and variance $I$.  Then
$$ S_n = \dfrac{n(n-1)}{2} \mu^2 + (n-1) \mu \sigma {\bf e}^T {\bf Z} + \dfrac{\sigma^2}{2} \left(({\bf e}^T {\bf Z})^2 - {\bf Z}^T {\bf Z}\right) $$  
Changing to an orthonormal basis where the first basis element is $n^{-1/2} \bf e$, this can be written as 
$$ S_n = \dfrac{n(n-1)}{2} \mu^2 + n^{1/2} (n-1) \mu \sigma W + \dfrac{\sigma^2(n-1)}{2} W^2 - \dfrac{\sigma^2}{2} U$$
where $W$ and $U$ are independent, $W$ is standard normal and $U$ is $\chi^2$ with $n-1$ degrees of freedom.
EDIT:
In the case $n=2$, you can take the orthogonal matrix $$M = \pmatrix{1/\sqrt{2} & -1/\sqrt{2}\cr 1/\sqrt{2} & 1/\sqrt{2}}$$
so that $${\bf e} = \sqrt{2} M \pmatrix{1\cr 0\cr}, \ {\bf Z} = M \bf W$$
where $\bf W$ is again multivariate normal with mean $0$ and variance $I$.
Then $${\bf e}^T {\bf Z} = \sqrt{2}\; [1\ 0]\; M^T M {\bf W} = \sqrt{2} W_1$$
where $W_1$, the first entry of $\bf W$, is standard normal, and
$$({\bf e}^T {\bf Z})^2 - {\bf Z}^T {\bf Z} = 2 W_1^2 - {\bf W}^T M^T M {\bf W} = 2 W_1^2 - {\bf W}^T {\bf W} = W_1^2 - W_2^2 $$
$W_1$ is what I called $W$ above, and $W_2^2 = U$.
