I'm looking for as many different arguments or derivations as possible that support the informal claim that Gaussian/Normal distributions are "the most fundamental" among all distributions.
A simple example: the central limit theorem (CLT) shows that the sum of i.i.d. random variables tends towards a Gaussian distribution.
The comments list many reasons why the Gaussian distribution is special, but is it "the most fundamental" among all distributions, as suggested in the OP? I would like to argue that (1) conservation laws are among the most fundamental laws of Nature, and (2) a quantity that obeys a conservation law will naturally follow an exponential -- rather than a Gaussian distribution.
Consider energy: two systems with energies $E_1$ and $E_2$ have total energy $E_1+E_2$, and if they are sufficiently large they will be independent, so the probability distribution must factorize: $P(E_1+E_2)=P(E_1)P(E_2)$, with the exponential distribution $P(E)\propto e^{-\beta E}$ as the unique normalizable solution (assuming $E$ is bounded from below).
This is the Gibbs measure. The Hammersley–Clifford theorem states that any probability measure which satisfies a Markov property is a Gibbs measure for an appropriate choice of (locally defined) energy function. The Gibbs measure is the fundamental measure of statistical physics, but it also applies widely outside of physics.
Economics is one such example. The statistical mechanics of money explains the exponential distribution of money as a direct consequence of the fact that money is conserved in general (there are exceptions). Wealth, in contrast, is not conserved (stocks may rise or fall), hence the non-exponential wealth distribution (Pareto distribution).
If the random vector $(X,Y)$ in the plane has independent coordinates and a rotation-invariant distribution, then it is Gaussian.
There is a whole book which addresses exactly this type of questions: suppose that a distribution has such and such properties, then it must be Gaussian (or sometimes Poisson).
MR0346969 Kagan, A. M.; Linnik, Yu. V.; Rao, C. Radhakrishna Characterization problems in mathematical statistics, John Wiley & Sons, New York-London-Sydney, 1973.
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 $$
Perhaps the reason lies in the characteristic function of the Gaussian distribution. It has several nice properties:
As mentioned above by prof. Yuval Peres it has the rotational invariance property.
Taking the Fourier transform of the normal distribution will again give you a normal distribution, and hence the central limit theorem. (This can also be viewed as a special consequence of the uncertainty principle.)
It has finite moments.
An $\mathbf{R}^n$ valued random variable $X = (X_1, ....., X_n)$ has the multivariate normal distribution if and only if any linear combination $a_1 X_1 + ..... + a_nX_n$ for deterministic real numbers $a_1 ,....... , a_n$ has the univariate normal distribution.
Although the 3rd point seems to be uninteresting in view of the several other distributions having this property, points 2 and 3 above are amongst the reasons why the Gaussian distribution is so common in nature.
When adding two independent random variables $X$ and $Y$, their respective moment generating functions $E\{\exp(t\,X))\}$ and $E\{\exp(t\,Y)\}$ are multiplied to yield the moment generating function of the sum $E\{\exp\bigl(t\,(X+Y)\bigr)\}$. This results in a convolution of the respective probability density functions.
If you take the logarithm of the moment generating functions (of the random variables and their sum), you get the cumulant generating functions that are added for the cumulant generating function of the sum of the independent variables.
Developing this function into a power series yields cumulants as successive coefficients.
The constant term is 0, the linear coefficient is the mean of the respective distributions (means of independent variables add when adding the variables), the quadratic coefficient is the variance, the next two terms are called skew and kurtosis. However, distributions for which all cumulants except mean and variance are zero are normal distributions.
The only cumulant guaranteed to be non-negative in any probability distribution is the variance: that is a significant contributor to the central limit theorem since variances cannot cancel each other.
If a Gaussian random process is stationary in the weak sense, then it is also stationary in the strict sense because the joint distribution of a (centred) Gaussian random vector is completely determined by the covariance matrix. This is not generally true for random processes with other distributions.