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:

$$|psirangle = sum_i^m sum_j^n a_{ij} |irangle |jrangle.$$

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.