Recursive convolution of a PDF with itself I was playing around with some simulations and I am trying to understand and produce an analytical solution for the behaviour I see. 
I generated $512$ sets of data, each set containing $10000$ data points. All data sets were drawn from a Rayleigh distribution with the same shape parameter $\sigma$. I label each set $X_{i}$, which has elements $X_{i} = \{x_{j},...,x_{M}\}$, where $M = 10000$.
I then take the average of the data sets as 
$$\bar{X}_{N} = \frac{1}{N} \sum_{i}^{N} X_{i} $$
to be explicit, not the average of the $X_{i}$, but the average of the sets. The resulting average still contains $M = 10000$ points. Plotting the averaged data sets gives the histograms, as follows:

What we can clearly see is that as more and more data sets are included in the average, the more normally distributed the resultant data set becomes. 
This is to be expected, I think, because of the central limit theorem. Now to my question:
I would like to show this mathematically. I know that the sum of random variables is described by the convolution from the parent distribution of those variables.
I imagine the result to be a two parameter distribution, $F(\sigma, N)$, $\sigma$ being the same shape parameter of the Rayleigh distribution. And $N$ the number of averages, or convolution iterations as per
$$F(x)_{N} = \int_{-\infty}^{+\infty} F_{N - 1}(X) f(x - X) \ dX$$
where $f(x)$ is the original distribution of interest (in my case the Rayeligh). 
Is it possible to get such a solution?

I have been able to calculate the first and iterations of the convolution. For the $N = 1$ case of the above convolution integral.
$$\int_{-\infty}^{+\infty} R(X) R(x - X) \ dX = \\ \frac{1}{4 \sigma^{3}}\exp\left(\frac{-x^{2}}{2\sigma^{2}}\right) \left( 2 x \sigma + \exp\left(\frac{x^{2}}{4\sigma^{2}}\right) \sqrt{\pi} \left( x^{2} - 2 \sigma^{2} \right) \rm{erf}(x/2 \sigma) \right)$$
Where $R(x)$ is the Rayleigh distribution. Plotting the above gives a reasonable result if I compare it to my data:

 A: Here is a quick outline determining the limiting distribution of the sample mean of $X_1,X_2,\ldots,X_n$ as $n\rightarrow \infty$ using moment generating functions.
All of the $X_i$ are independent and have the same distribution:  $X_i \sim \text{Rayleigh}(\sigma)$.  We have $\mu=\text{E}(X_i)=\sigma \sqrt{\pi/2}$.  Let 
$$S_n=\sum_{i=1}^n X_i$$
We look at the limiting distribution of 
$$Z=\sqrt{n}(S_n/n-\mu)$$
using moment generating functions.  The moment generating function of a Rayleigh distribution is
$$M(t)=\sqrt{\frac{\pi }{2}} \sigma  t \left(\text{erf}\left(\frac{\sigma  t}{\sqrt{2}}\right)+1\right) \exp \left(\frac{\sigma ^2 t^2}{2}\right)+1$$
The moment generating function of $S_n$ is simply $M(t)^n$.  We also know that the moment generating function of $a S_n+b$ (where $a$ and $b$ are known constants) is $e^{b t} M(a t)^n$.  So the moment generating function of $Z=\sqrt{n}(S_n/n-\mu)$ is
$$e^{\sqrt{\frac{\pi }{2}} \left(-\sqrt{n}\right) \sigma  t} \left(\frac{\sqrt{\frac{\pi }{2}} \sigma  t e^{\frac{\sigma ^2 t^2}{2 n}} \left(\text{erf}\left(\frac{\sigma  t}{\sqrt{2} \sqrt{n}}\right)+1\right)}{\sqrt{n}}+1\right)^n$$
Taking the limit as $n->\infty$ we have
$$e^{\frac{1}{4} (4-\pi) \sigma ^2 t^2}$$
for the moment generating function.  As the moment generating function for a normal distrbution with mean $m$ and variance $v$ is
$$e^{m t+\frac{t^2 v}{2}}$$
we see that our limiting moment generating function represents a normal distribution with mean $0$ and variance $\frac{1}{2} (4-\pi ) \sigma ^2$.  What this means is that we can use a normal distribution with mean $\sigma \sqrt{\pi/2}$ and variance $\frac{1}{2} (4-\pi ) \sigma ^2/n$ to approximate the distribution of $\bar{X}=S_n/n$.
As a partial check consider $n=10$.  We can generate a bunch of samples and plot the resulting histogram and approximate density (using Mathematica)
n = 10;
\[Sigma] = 1;
z = RandomVariate[RayleighDistribution[1], {1000000, n}];
z = Mean[#] & /@ z;
Show[Histogram[z, 100, "PDF", Frame -> True, 
  FrameLabel -> (Style[#, Bold, 18] &) /@ {"\!\(\*OverscriptBox[\(X\), \(_\)]\)", "Probability density"}],
 Plot[PDF[NormalDistribution[Sqrt[\[Pi]/2] \[Sigma], Sqrt[(4 - \[Pi]) \[Sigma]^2/(2 n)]], x],
  {x, 0.4, 3}, PlotRange -> All]]


