Sum of Gaussian pdfs

First of all this has nothing to do with the inflection point of $e^{-\alpha x^2}$. According to Poisson summation formula (see Whittaker, Watson, Modern analysis, chapter 21.51) $$ \sum_{n=-\infty}^\infty e^{-\alpha (x-n)^2}=2{\sqrt{\frac{\pi}{\alpha}}}\left(1+2\sum_{n=1}^\infty e^{-\frac{\pi^2}{\alpha}n^2}\cos2\pi n x\right),\quad \text{Re}~\alpha>0.\tag{1} $$ From this formula one can see that when $\alpha>0$ is small, the Fourier coefficients $a_n=2{\sqrt{\frac{\pi}{\alpha}}}e^{-\frac{\pi^2}{\alpha}n^2}$ of the resulting function decreases rapidly with increasing $n$. When $\alpha=1/2$ one has $$ \frac{a_1}{a_0}=2e^{-2\pi^2}\approx 5.4\times 10^{-9}. $$

I don't know any conceptual explanation for this in mathematics, but there is such an explanation that comes from physics. Consider a quantum particle with mass $m$ on a ring of radius $a$. We assume the ring is pierced with magnetic flux $\phi$. We want to calculate partition function of this system at temperature $T>0$ in two different ways.

On the one hand, it is known that the energy spectrum of the particle is given by $$ E_n=\frac{\hbar^2}{2m^2}\left(n+\frac{\phi}{\phi_0}\right)^2, $$ where $\phi_0$ is the so called magnetic flux quantum. Then partition function is given as the Gibbs sum $$ Z=\sum_{n=-\infty}^\infty e^{-E_n/T}=\sum_{n=-\infty}^\infty e^{-\frac{\hbar^2}{2m^2T}\left(n+{\phi}/{\phi_0}\right)^2}.\tag{2} $$ Here one immediately recognizes the sum analogous to the LHS of $(1)$.

On the other hand, partition function is related to the trace of the density matrix $\rho_{\theta_1,\theta_2}$, i.e. to $\int\rho_{\theta,\theta}d\theta$. It is possible to calculate the density matrix of this system in imaginary time representation by solving a certain differential equation. Details of this calculation can be found for example in this book, chapter 4.3. Consider first the more simple case of unbounded line $-\infty<\theta<+\infty$; then the answer is $$ \rho_{\theta_1,\theta_2}=Ce^{-\frac{mTa^2(\theta_1-\theta_2)^2}{2\hbar^2}}, $$ where $C$ is some normalization constant. Magnetic flux does not enter this expression because there is not any nontrivial loop in the system. On a ring there is a nontrivial loop, and let $n$ be the winding number. It is known that the partition function is a sum over all homotopy classes multiplied by corresponding phase factor. In this case the phase factor comes from the magnetic flux (Aharonov-Bohm phase) and equals $e^{2\pi in\phi/\phi_0}$. As a result we have \begin{align} Z&=\sum_{n=-\infty}^\infty e^{2\pi in\phi/\phi_0}\int\rho_{\theta,\theta+2\pi n}d\theta\\ &=2\pi C \sum_{n=-\infty}^\infty e^{-\frac{2\pi^2mTa^2}{\hbar^2}n^2+2\pi in\phi/\phi_0}.\tag{3} \end{align} Combining $(2)$ and $(3)$ one obtains the transformation in $(1)$.

One can see from this physical interpretation that when the temperature $T$ increases, then the higher Fourier harmonics decrease. This corresponds to the physical intuition that the higher the temperature the more chaotic the system becomes and the effect of the Aharonov-Bohm phase averages out due to thermal fluctuations.

Interesting! Put another way, the fractional part of the standard Gaussian variable closely approximates the standard uniform variable.

You can also closely approximate standard Gaussian from standard uniform: Add 12 independent uniform variables before subtracting 6. This approximation is also shockingly good; see John D. Cook's comment here. (12 isn't "large," but rather uses var(U)=1/12 to get unit variance.)

Notice that the fractional part of this approximate Gaussian is exactly uniform. You might think of this as a dual explanation for the phenomenon you observed.

Too long for a comment...

I'm not sure why this is so surprising. If one sums translates by $\varepsilon \mathbb{Z}$ of the standard Gaussian then for $\varepsilon$ sufficiently small the result ought to be roughly constant. More generally, one can take any sufficiently nice function in lieu of the Gaussian, albeit with the awareness that this will affect $\varepsilon$.

The fact that $\varepsilon = 1$ is "small" for the standard Gaussian is the only thing that should be surprising. So let's set $z := -i\varepsilon^2/4\pi$ and $\tau := i\varepsilon^2/2\pi$ and compare

$\sum_n \exp(-[\varepsilon n]^2/2) = \vartheta(0;\tau) = \vartheta_{01}(z;\tau)$

and

$\sum_n \exp(-[\varepsilon(n-1/2)]^2/2) = \exp(-\varepsilon^2/8) \cdot \vartheta(z;\tau) = \exp(-\pi i[\tau/4+z]) \cdot \vartheta(z;\tau)$.

Now the question becomes when these two functions are approximately equal.