Distribution of a certain functional of iid $N(0,1)$ random variables 
Suppose that $X_1,\ldots,X_n$ are iid $N(0,1)$ random variables. Consider the random variable given by
$$
\xi_n
=\Bigl|\frac1{\sqrt{n}}\sum_{t=1}^nX_t\Bigr|^2-\frac1n\sum_{t=1}^nX_t^2
=\frac1n\sum_{s\ne t}X_sX_t.
$$
What is the distribution of $\xi_n$?

My observations are as follow.

*

*We have that $n^{-1/2}\sum_{t=1}^nX_t\sim N(0,1)$ and $|n^{-1/2}\sum_{t=1}^nX_t|^2\sim\chi^2(1)$, where $\chi^2(k)$ is the chi-square distribution with $k$ degrees of freedom.

*$\sum_{t=1}^nX_t^2\sim\chi^2(n)$ and $n^{-1}\sum_{t=1}^nX_t^2\sim\Gamma(k=n/2,\theta = 2n^{-1})$ (see here).

*$|n^{-1/2}\sum_{t=1}^nX_t|^2$ and $n^{-1}\sum_{t=1}^nX_t^2$ are not independent since the covariance between them is equal to $2/n$ (although the covariance goes to $0$ as $n\to\infty$).

*$n^{-1}\sum_{t=1}^nX_t^2\to 1$ in probability as $n\to\infty$ using the law of large numbers and hence it follows that $\xi_n\to \chi^2(1)-1$ in distribution as $n\to\infty$ using Slutsky's theorem. However, I am interested if the exact distribution has some manageable form when $n$ is some fixed positive integer.

If $|n^{-1/2}\sum_{t=1}^nX_t|^2$ and $n^{-1}\sum_{t=1}^nX_t^2$ were independent, then we could try to use the answers to this question to derive the exact distribution. But they are not independent. Maybe Cochran's theorem might be of use here.
Does the exact distribution of $\xi_n$ have some manageable expression?
Any help is much appreciated!
 A: One has
$$\xi_n=(n-1)\bar X^2-\frac{n-1}n\,S^2=U+V,$$
where
$$\bar X:=\frac1n\,\sum_1^n X_i,\quad S^2:=\frac1{n-1}\,\sum_1^n(X_i-\bar X)^2
=\frac1{n-1}\Big(\sum_1^n X_i^2-n\bar X^2\Big),$$
$$U:=\frac{n-1}n Z^2,\quad Z:=\bar X\sqrt n,\quad V:=-\frac1n\,\chi_{n-1}^2,\quad
\chi_{n-1}^2:=(n-1)S^2.$$
Further, $Z\sim N(0,1)$, $\chi_{n-1}^2=(n-1)S^2$ has the chi-squared distribution with $n-1$ degrees of freedom, and $\bar X$ and $S^2$ are independent. So, $U$ and $V$ are independent random variables with easily found pdf's, and the pdf of $\xi_n$ is the convolution of the pdf's of $U$ and $V$.
A: For $X_1,\dots, X_n$ Gaussian random variables with variance-covariance matrix $M$, and $L$ a $m \times n$ matrix, the variables $Y_1,\dots, Y_m$ given by $Y_i = \sum_{j=1}^n L_{ij} X_j$ have variance-covariance matrix $L M L^T$.
This follows immediately from writing the covariance as an expectation of products of sums, exchanging the sum and product, and then exchanging the sum and expectation.
In particular, if $M$ is the identity matrix, and $L$ is an $n \times n$ orthogonal matrix, then the $Y_i$ have variance-covariance matrix the identity, hence are independent standard Gaussians (because Gaussians are determined by their variance-covariance matrix (and mean, which is zero throughout this discussion and hence can be ignored)).
Given any quadratic function $\sum_{i=1}^n \sum_{j=1}^n C_{ij} X_i X_j$ of the $n$ Gaussian random variables, determined by a matrix $C$, which without loss of generality we can assume is symmetric, we can understand the distribution of $C$ by changing variables from the $X_i$ to the $Y_i$. This has the effect of sending $C$ to $ L^{-1} C L^{-T}$, which if $L$ is orthogonal is $L^T C L$.
Since every symmetric matrix can be diagonalized, we can make $C$ a diagonal matrix by this approach.
In your case $C$ is the matrix with off-diagonal entries $\frac{1}{n}$ and diagonal entries $0$. By viewing this as the all-$\frac{1}{n}$s matrix minus a diaogonal matrix, we see it has eigenvalues $\frac{n-1}{n}, \frac{-1}{n},\dots, \frac{-1}{n}$, and thus diagonalizing it turns it into a diagonal matrix with entries $\frac{n-1}{n}, \frac{-1}{n},\dots, \frac{-1}{n}$.
To be explicit we can find a specific orthonormal basis on which the matrix diagoanlizes, of which one vector must be $\frac{1}{\sqrt{n}},\dots, \frac{1}{\sqrt{n}}$ and the rest can be any orthonormal basis orthogonal to it - one of which I wrote down in the comment.
