This is a heuristic question that I think was once asked by Serge Lang. The gaussian: $e^{x^2}$ appears as the fixed point to the Fourier transform, in the punchline to the central limit theorem, as the solution to the heat equation, in a very nice proof of the AtiyahSinger index theorem etc. Is this an artifact of the techniques (such as the Fourier Transform) that people like use to deal with certain problems or is this the tip of some deeper platonic iceberg?

Quadratic (or bilinear) forms appear naturally throughout mathematics, for instance via inner product structures, or via dualisation of a linear transformation, or via Taylor expansion around the linearisation of a nonlinear operator. The LaplaceBeltrami operator and similar secondorder operators can be viewed as differential quadratic forms, for instance. A Gaussian is basically the multiplicative or exponentiated version of a quadratic form, so it is quite natural that it comes up in multiplicative contexts, especially on spaces (such as Euclidean space) in which a natural bilinear or quadratic structure is already present. Perhaps the one minor miracle, though, is that the Fourier transform of a Gaussian is again a Gaussian, although once one realises that the Fourier kernel is also an exponentiated bilinear form, this is not so surprising. But it does amplify the previous paragraph: thanks to Fourier duality, Gaussians not only come up in the context of spatial multiplication, but also frequency multiplication (e.g. convolutions, and hence CLT, or heat kernels). One can also take an adelic viewpoint. When studying nonarchimedean fields such as the padics $Q_p$, compact subgroups such as $Z_p$ play a pivotal role. On the reals, it seems the natural analogue of these compact subgroups are the Gaussians (cf. Tate's thesis). One can sort of justify the existence and central role of Gaussians on the grounds that the real number system "needs" something like the compact subgroups that its nonarchimedean siblings enjoy, though this doesn't fully explain why Gaussians would then be exponentiated quadratic in nature. 


I recently came across a strange and beautiful connection between the Gaussian $e^{x^2}$ and the method of least squares. It turns out that the square in $e^{x^2}$ and the square in ``least squares'' is the same square. Let $(x_i,y_i)$ (with $1\leq i \leq n$) be the data set, and assume that for each $x$, the $y$'s are normally distributed with mean $\mu(x)=\alpha x+\beta$ and variance $\sigma^2$. Then, the likelihood of generating our data (assuming that the data points are independent) is $$\prod_{i=1}^n \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(\frac{(y_i\mu(x_i))^2}{2\sigma^2}\right) =\left(\frac{1}{\sqrt{2\pi \sigma^2}}\right)^n \exp\left( \frac{1}{2\sigma^2} \sum_{i=1}^n (y_i  \alpha x+\beta)^2 \right)$$ We would obviously want to choose the parameters $\alpha,\beta$ so that the likelihood is maximized, and this is accomplished by minimizing $$\sum_{i=1}^n (y_i  \alpha x+\beta)^2.$$ In other words, the least squares approximation is the one that makes the data set most likely to happen. 


The sum of two independent normal random variables is again normal, i.e., the shape of the distribution is unchanged under addition except for stretching and scaling. Moreover, the normal distribution is unique among distributions with finite variance in having this property. Many phenomena in nature come from adding together various independent or almost independent terms. Therefore, we would expect the normal distribution to show up a lot in natureinspired mathematics. 


This is just a minor amplification of one of Terry Tao's points. For any prime $p$, the ring $\mathbb{Z}_p$ of $p$adic integers forms an open compact additive subgroup of $\mathbb{Q}_p$, the completion of $\mathbb{Q}$ under the padic metric, and its characteristic function should be viewed as a padic analogue of the Gaussian. It displays many analogues of the nice properties we see from the Gaussian:
The situation with the real line is more complicated because it is connected but not compact, and therefore has no open compact subgroups. There is a maximal compact multiplicative monoid (occasionally called the "ring of integers of $\mathbb{R}$" informally), given by the closed interval $[1,1]$. You can think of the Gaussian as a smoothing of the uniform distribution on $[1,1]$, but it is not clear to me that this particular analogy is very fruitful. 


I'm not an expert, but I believe Stein's method gives a more satisfying connection between the CLT and the heat equation, in particular one that does not involve the Fourier transform. Stein's characterization of the normal distribution and convergence to normality involves an operator closely related to the generator of the OrnsteinUhlenbeck process. On the other hand, the fact that latter has the Gaussian as its invariant measure can be obtained by a trivial transformation of the fundamental solution of the heat equation. 


One of the reasons for the ubiquity of the Gaussian is displayed in what is probably the most electrifying half page of scientific prose ever writtenMaxwell's argument that the distribution of the velocities of molecules in the ideal gas is Gaussian (now known as the MaxwellBoltzmann distribution). The only physical assumptions used are that the density function depends only on the absolute value of the velocity (and not the direction) and that the components in the directions of the coordinate axes are statistically independent. Mathematically, ths means that the only functions in $3$space which depend only on the distance $r$ from the origin and which split as the product of three functions of one variable are those of the form $ae^{br^2}$. Maxwell does this by inspection but it is easy to give a rigorous proof (under very weak smoothness conditions) and the result holds, of course, in any dimension greater then or equal to $2$. Maxwell's reasoning can be found in his collected papers, or, more accesssibly, in Hawking's anthology "On the Shoulders of Giants". 

