I have some random variables $X_1, X_2, X_3, ...$ that are normally distributed with mean 0 and standard deviation $\sigma$. I am trying to model a function of these variables that includes a double-sum of their difference: $$ \sum_{n, m\neq n}^{N,N} (X_n - X_m)^2 $$
I tried to model this sum as a chi-squared distribution with $(N^2-N)$ degrees of freedom (the -N comes from including every term except the terms where $n=m$). Then this sum should have mean $2\sigma^2 (N^2-N)$ and variance $8 \sigma^4 (N^2-N)$. But when I try to confirm that numerically in MATLAB, I find the variance is much larger. I assume this must be because the terms of the sum are not independent, they are using the same variables $X_n$ so they are correlated. But I don't know a way to include the correlations. Is there a neat way to evaluate the distribution of this sum?