The Normal Distribution is the limit, in the sense of distributions, of the scaled sum of $n$ IID variables. This is the Central Limit Theorem.

I posted an outline of a proof of the Central Limit Theorem here:

In what follows, the fourier transform used is
$$
\widehat{f}(\xi)=\int_{-\infty}^\infty f(x)\,e^{-i2\pi x\xi}\;\mathrm{d}x
$$
Suppose we have a probability density function, $\phi$, which has a mean of $0$ and a variance of $1$ (this can be achieved by translation and scaling).
Then, the following are true for the fourier transform of the pdf:

$\widehat{\phi}(0)=1$ (i.e. $\int_{-\infty}^\infty\phi(x)\:\mathrm{d}x=1$; $\phi$ is a probability measure)

$\widehat{\phi}\vphantom{\phi}^\prime(0)=0$ (i.e. $\int_{-\infty}^\infty (-i2\pi x)\phi(x)\:\mathrm{d}x=0$; the mean of $\phi$ is $0$)

$\widehat{\phi}\vphantom{\phi}^{\prime\prime}(0)=-4\pi^2$ (i.e. $\int_{-\infty}^\infty (-i2\pi x)^2\phi(x)\:\mathrm{d}x=-4\pi^2$; the variance of the pdf is $1$)

Thus, to second order, the fourier transform of the pdf looks like
$$
\widehat{\phi}(\xi)=1-2\pi^2\xi^2\tag{1}
$$
The Central Limit Theorem looks at the sum of $k$ random variates contracted by $\sqrt{k}$ and scaled by $\sqrt{k}$. This maintains the unit integral while
compensating for the increased variance due to summation.

The pdf of the sum of random variates is the convolution of the pdf of the variates. The fourier transform of a convolution is the product of the fourier transforms. The fourier transform of a function contracted by $\sqrt{k}$ and scaled by $\sqrt{k}$ is the fourier transform expanded by $\sqrt{k}$.

Thus, the fourier transform of the pdf of the sum of $k$ independent trials contracted and scaled appropriately is
$$
\left(1-\frac{2\pi^2\xi^2}{k}\right)^k\tag{2}
$$
which, for large $k$ approaches
$$
e^{-2\pi^2\xi^2}\tag{3}
$$
the inverse fourier transform of which is
$$
\frac{1}{\sqrt{2\pi}}e^{-x^2/2}\tag{4}
$$
The distribution in $(4)$ is gaussian normal distribution with mean $0$ and variance $1$.

The Normal Distribution is the probability distribution with mean $0$ and variance $\sigma^2$ that maximizes entropy.

I posted a derivation of the distribution which maximizes entropy here:

This answer gives an outline of how to use the Fourier Transform to prove that the $n$-fold convolution of any probability distribution with a finite variance contracted by a factor of $\sqrt{n}$ converges weakly to the normal distribution.

However, in his answer, Qiaochu Yuan mentions that one can use the Principle of Maximum Entropy to get a normal distribution. Below, I have endeavored to do just that using the Calculus of Variations.

**Applying the Principle of Maximum Entropy**

Suppose we want to maximize the entropy
$$
-\int_{\mathbb{R}}\log(f(x))f(x)\,\mathrm{d}x\tag1
$$
over all $f$ whose mean is $0$ and variance is $\sigma^2$, that is
$$
\int_{\mathbb{R}}\left(1,x,x^2\right)f(x)\,\mathrm{d}x=\left(1,0,\sigma^2\right)\tag2
$$
That is, we want the variation of $(1)$ to vanish
$$
\int_{\mathbb{R}}(1+\log(f(x)))\,\delta f(x)\,\mathrm{d}x=0\tag3
$$
for all variations of $f$, $\delta f(x)$, so that the variation of $(2)$ vanishes
$$
\int_{\mathbb{R}}\left(1,x,x^2\right)\delta f(x)\,\mathrm{d}x=(0,0,0)\tag4
$$
$(3)$, $(4)$, and orthogonality requires
$$
\log(f(x))=c_0+c_1x+c_2x^2\tag5
$$
To satisfy $(2)$, we need $c_0=-\frac12\log\left(2\pi\sigma^2\right)$, $c_1=0$, and $c_2=-\frac1{2\sigma^2}$. That is,
$$
\bbox[5px,border:2px solid #C0A000]{f(x)=\frac1{\sigma\sqrt{2\pi}}\,e^{-\frac{x^2}{2\sigma^2}}}\tag6
$$

square-integrablerandom variables tends toward Gaussian. People interested in heavy tails might disagree that this class is "most fundamental". $\endgroup$11more comments