Why does the probability transition from a Dirac distribution to a Gaussian? I'm trying to follow this paper in Section 3.1, but I'm having trouble with a comment they make in the final paragraph of that section.
First, we start with a pure state:
$$|\psi\rangle = \sum_i^m \sum_j^n a_{ij} |i\rangle |j\rangle.$$
Then, they consider the complex numbers $a_{ij}$ being uniformly distributed on a hypersphere, so they have probability density:
$$P(a) \sim \delta \left( \sum_i^m \sum_j^n \lvert a_{ij} \rvert^2 - 1 \right),$$
where $\delta()$ is the Dirac distribution. Then, in the last paragraph before Section 3.2, they write:

Second, notice that the exact result of Lubkin can be estimated by relaxing the normalization constraint in the distribution, and replacing it with a product of independent Gaussian distributions, $P(a) = \prod_{i,j} (nm / \pi) \exp \left( -nm \lvert a_{ij} \vert^2 \right)$, with $\langle a_{ij} \rangle = 0$ and $\langle \lvert a_{ij} \rvert ^2 \rangle = 1/nm$.

They then go on to say that the central limit theorem says that this distribution tends to a Gaussian in $\sum_i^m \sum_j^n \lvert a_{ij}\rvert^2$ centered at $1$, with variance $1/\sqrt{nm}$.
Basically, I don't see why we get this distribution. From what I can tell, we need a sum of random variables (which we have when you take the product of those exponentials), but I'm not seeing the mean of $1$, nor the (few) steps being taken to get this end result. If someone could explain, that would be great. Also, I would like to understand what the assumption behind replacing our initial normalization with the Gaussians was.
 A: I think that the gaussian distribution just emerges from the definition of the Dirac delta, but I'm only able to understand the approximation in a reversed way (starting from the delta and ending up with the product of gaussians). The Dirac function is defined as:
$$\delta(x)=\lim_{\sigma\rightarrow 0}\frac{1}{\sqrt{2\pi}\sigma}e^{-x^2/2\sigma^2},$$
where $\sigma^2$ is the variance (or width) of the gaussian distribution, centered around zero in this case.
Using this definition, you can have:
$$P(\alpha)\sim \delta (\sum^m_i\sum^n_j|\alpha_{ij}|-1)=\lim_{\sigma\rightarrow 0}\frac{1}{\sqrt{2\pi}\sigma}e^{-(\sum^m_i\sum^n_j|\alpha_{ij}|-1)^2/2\sigma^2},$$
where $\sigma$ must be $\sigma=1/mn$ since it is the parameter that controls the width of the distribution (in the limit of $mn\rightarrow \infty$, your system should have a perfect delta). Then they "relax" the normalization assuming a finite but very large dimension, so:
$$\delta (\sum^m_i\sum^n_j|\alpha_{ij}|-1)\simeq\frac{nm}{\sqrt{2\pi}}e^{-nm(\sum^m_i\sum^n_j|\alpha_{ij}|-1)^2/2}.$$
Finally, they apply the limit central theorem to write the distribution as a product of distributions. For me, this derivation is a bit weird, but it is the only way I find to explain it. I hope it can help you!
