intuition on why the density function of a normal law is$\frac{1}{\sqrt{2\pi \sigma^2 }}e^{-\frac{(x-\mu)^2}{2\sigma^2 }}$

What is the intuition of the fact that the density function of a normal law $\mathcal N(\mu,\sigma ^2)$ is given by $$f(x)=\frac{1}{\sqrt{2\pi \sigma^2 }}e^{\large-\frac{(x-\mu)^2}{2\sigma^2 }}?$$ It looks to come from nowhere. For example, the intuition of a poisson law is in fact very natural considering rare event (and binomial law). For the normal law, I really have no idea of the intuition behind, neither why it's so common in the nature ! How did we get to this (incredible) result. How does it work ?

• A start of answer : The coefficient before the exponential is in fact kind of artificial and can be ignored since it is used to "scale" the density function to $1$ : $\int_{-\infty}^{\infty} f(x) dx = 1$ – Zubzub Sep 24 '16 at 16:47
• I corrected it, thank you. – MathBeginner Sep 24 '16 at 17:14
• The limit of binomial coefficients corresponding to the central limit theorem, might be the most natural explanation (and the historically relevant one). – Did Sep 24 '16 at 17:16
• I would say: math.arizona.edu/~jwatkins/505d/Lesson_10.pdf – georg Sep 24 '16 at 17:29

Distributions are often characterised by probabilities or probability densities, but one can also use the characteristic function of a variable $X$, defined as $\varphi_X\left( t\right) := E\left[ e^{itX}\right]$. The pdf of a continuous random variable can be obtained as $\frac{1}{2\pi}\int_{-\infty}^\infty \varphi_X\left( t\right) e^{-itx}dt$; you're basically inverting a Fourier transform. The characteristic function has some conceptually simpler cousins that are not in general well-defined, such as the moment-generating function $E\left[ e^{tX}\right]$, so called because its derivatives' values at $t=0$ are "moments" (means of powers of $X$), and the probability-generating function $E\left[ t^X\right]$ if $X$ has a discrete distribution, in which case the $t^k$ coefficient is the probability that $X=k$ (hence the name).
The point of the technical overview above is that an intuitive reason why a distribution is what it is can be provided not just by deriving probabilities, but also by deriving functions such as those discussed above. I'm sure you're familiar with the former motivation of the Poisson distribution, so let's do it the latter way instead to illustrate how it works. The binomial distribution has probability-generating function $$\sum_{k=0}^n {}^nC_k \left( pt\right)^k q^{n-k}=\left( q+pt\right)^n.$$Fix the mean $np=\lambda$ so as $n\to\infty$ the pgf becomes $$\left( 1+\left( t-1\right)\frac{\lambda}{n}\right)^n=\exp\lambda\left( t-1\right).$$ The $t^k$ coefficient is then $e^{-\lambda}\lambda^k/k!$, but then you already knew that.
If $X_1,\,X_2,\,\cdots X_n$ are independent variables from a distribution of mean $\mu$, standard deviation $\sigma$, then $S_n:=\frac{X-\mu}{\sigma\sqrt{n}}$ has an approximate $N\left( 0,\,1\right)$ distribution for large $n$. That's effectively how one defines Normal distributions; we can restate $Z\sim N\left( \mu,\,\sigma^2\right)$ as $\frac{Z-\mu}{\sigma}\sim N\left( 0,\,1\right)$. While @Did referred to this historical result, there's another way to look at it. Let's use the moment-generating function, $M_X\left( t\right):=E\left[ e^{tX}\right]$. Since $M_{aX+b}\left( t\right)=e^{bt}M_X\left( at\right)$, the moments of each $\frac{X_i-\mu}{\sigma\sqrt{n}}$ vanish as $n\to\infty$, except for the mean and variance, which are $0$ and $1$. So $\frac{X_i-\mu}{\sigma\sqrt{n}}$ has mgf is $1+\frac{t^2}{2n}+o\left( t^2\right)$, and $S_n$ has mgf $\left( 1+\frac{t^2}{2n}+o\left( t^2\right)\right)^n=e^{t^2/2}$.
Specifying this mgf is equivalent to specifying a pdf. Anyone who claims the pdf is $\frac{1}{\sqrt{2\pi}}e^{-x^2/2}$ can prove it simply by checking it gives the right mgf, viz. $$\int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}}e^{tx-x^2/2}dx=\int_{-\infty}^\infty \frac{1}{\sqrt{2\pi}}e^{t^2/2-y^2/2}dy=e^{t^2/2}.$$Alternatively (warning: complex numbers ahead), you could use a slight variant on the above argument to prove the characteristic function is $e^{-t^2/2}$, then compute the pdf as$\frac{1}{2\pi}\int_{-\infty}^\infty e^{-itx-t^2/2}dt$ (although that's a difficult calculation; see here).
From all probability density functions with the same variance, the Gaussian has the maximum entropy. If you seek the probability density function that maximizes entropy, you will arrive at a function of the form $e^{-x^2}$ and the rest is just normalization. There is also the Central Limit Theorem going for it.