Discrepancy between Gaussian FFT and its Fourier Transform
First consider the theoretical DTFT of $x[n]$
$$ X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j \omega n} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1 - a e^{-j \omega} }$$

Then consider the N-point DFT of the truncated signal:
$$ X[k] = \sum_{n=0}^{N-1} x[n] e^{-j \frac{2\pi}{N} n k} = \sum_{n=0}^{N-1} a^n e^{-j \frac{2\pi}{N} n k} = \frac{1 - a^N}{1 - a e^{-j \omega} }$$

A simple inspection of those DTFT abd DFT yields the following conclusion:
$$ \boxed{ X[k] = (1-a^N) X(e^{j\omega_k}) } $$

Where the DTFT is evaluated at the DFT frequency of $w_k = \frac{2\pi}{N}k$.

And the error between $X[k]$ and $X(e^{j \omega_k})$ is:
$$ E[k] = X(e^{j \omega_k}) - X[k] = a^N X(e^{j \omega_k}) $$

As can be seen, the error due to the finite length DFT truncation is multiplied by the constant  $a^N$ and the frequency dependent variable $X(e^{j \omega_k})$ for each DFT frequency of $\omega_k$. It's obvious from this that the error plot will have a shape of the actual DTFT $X(w)$ scaled by $a^N$.

In order to make this error smaller than the floating point round-off errors (such as $10^{-16}$ for 64 bit IEEE floating point numbers) you have to make sure that for all frequencies $\omega_k$ for $0 \leq k \leq N-1$ you should satisfy $E[k] &lt; 10^{-16}$ that is:$$ E[k]=a^N X(e^{j \omega_k}) &lt; 10^{-16}$$

Note that this error is frequency dependent due to $X(e^{j \omega_k})$ term: the higher the DTFT at that frequency, the higher will be the error and vice versa. So for a randomly fluctuating DTFT $X(e^{j \omega_k})$ the error will also randomly fluctuate.

In the analysis of the Gaussian signal $x[n] = e^{-\alpha n^2} = a^{n^2}$ where $ 0 &lt; a = e^{-\alpha} &lt; 1$ , an analytic expression for the error will be difficult to simplify as you have to evaluate the following sum:

$$ X_g(e^{j\omega}) = \sum_{n=-\infty}^{\infty} a^{n^2} e^{-j \omega n} = \beta a^{\gamma w^2}$$ where $\beta$ and $\gamma$ should be determined.

And the N-point DFT of the truncated (and shifted by $N/2$) Gaussian signal:
$$ X_g[k] = \sum_{n=0}^{N-1} a^{(n-N/2)^2} e^{-j \frac{2\pi}{N} n k} = ? $$

Unfortunately since the series sum may not have a closed form expression (or it will be quite complicated) therefore a simple formulation of the error does not exist as were the case for the simpler signal $x[n] = a^n u[n]$. 

EDIT: The following exposition shows the same results as attained above however from a different perspective which may be less intuitive to follow. I suggest the above perspective for a clear understanding and the below perspective a complementary case.

Consider the following infinite length signal $x[n] = a^n u[n]$, with $|a|&lt;1$, where the discrete-time Fourier tranform (DTFT) $X(e^{j\omega})$ of $x[n]$ is:

$$X(e^{j\omega}) = \text{DTFT} \{x[n] \} = \sum_{n=0}^{\infty} a^n e^{-j \omega n} = \frac{1}{1-a e^{-j\omega} } $$

In order to practically compute the DTFT $X(e^{j\omega})$ using a digital computer, the relative transform, DFT (discrete Fourier transform) is used. For finite length signals it's well known that DFT $X[k]$ will be samples of DTFT $X(e^{j \omega})$. But for infinite length signals (such as $x[n]=a^n u[n]$) the computed N-point DFT (which should practically use a finite length truncated version of the infinite length signal) will be an approximation to the true theoretical DTFT. The question is ; what's the relation between the true DTFT spectrum $X(e^{j\omega})$ and the computed DFT samples  $X[k]$ ?

Now it's obvious that for infinite length signals the values of DFT $X[k]$ will not be the samples of the function $X(e^{j\omega})$. For example given $x[n]=a^nu[n]$ and $X(e^{j \omega}) = \frac{1}{1-a e^{-j\omega}}$ then we see that  $$X[k] \neq X(e^{j\frac{2 \pi}{N}k}) = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$$

First lets consider the following, assume you have artificially obtained $N$ uniform samples from the function $X(e^{j\omega})=\frac{1}{1-a e^{-j\frac{2 \pi}{N}k}}$ as
$$ X_f[k] = \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} $$ where $X_f[k]$ stands for the  DFT samples of a finite length signal $\tilde{y}[n]$ whose N-point DFT would be $X_f[k]$. Therefore we can assert that the inverse DFT of the $X_f[k]$ would yield the N-point discrete-time sequence $\tilde{y}[n]$ as:

$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \frac{1}{N} \sum_{k=0}^{N-1}  \frac{1}{1-a e^{-j\frac{2 \pi}{N}k}} e^{ j \frac{2\pi}{N} k n} $$

It can be shown from the periodicity assumption of the IDFT that the result is also equivalent to:
$$ \tilde{y}[n] = IDFT\{ X_f[k] \} = \sum_{r=-\infty}^{\infty}  x[n-rN] $$

Where the signal $x[n]$ is nothing but the theoretical infinite length inverse of the DTFT $X(e^{j\omega}) =  \frac{1}{1-ae^{-j\omega}}$ which is as we know 
$$ x[n] = IDTFT \{\frac{1}{1-ae^{-j\omega}} \} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \frac{1}{1-a e^{-j\omega}} e^{j\omega n} d\omega = a^n u[n] $$ 

Therefore we have the following:

$$\tilde{y}[n] = \sum_{r=-\infty}^{\infty}  x[n-rN] = \sum_{r=-\infty}^{\infty}  a^{n-rN}u[n-rN] $$

It will be interesting to consider the value of $\tilde{y}[n]$ for the range $0 \leq n \leq N-1$ for which we have:

$$\tilde{y}[n] = a^n \sum_{r=-\infty}^{\infty}  a^{-rN}u[n-rN] = a^n \sum_{r=-\infty}^{0}  a^{-rN} = a^n \sum_{r=0}^{\infty}  (a^N)^r = a^n \frac{1}{1-a^N} $$

Which shows the relation between finite length signal $\tilde{y }[n]$ and the infinite length signal $x[n] = a^n u[n]$ as: 
$$\tilde{y}[n] = \frac{x[n]}{1-a^N}$$ where $x[n] = a^n u[n]$.

One can see that the modification involved on the computed value of finite length approximation to the infinite length signal $x[n]$ is a constant weight of $K = \frac{1}{1-a^N}$ for the signal $x[n]= a^n u[n]$

Now using the linearity properties of DFT one can deduce the following:
Given the pair of N-point sequence $\tilde{y}[n]$ and its N-point DFT $X_f[k]$ which is the uniform samples of the true theoretical DTFT $X(e^{j\omega})$ associated with infinite length signal $x[n]=a^n u[n]$ we have the following 

$$\tilde{y}[n] \leftrightarrow X_f[k] $$
$$\tilde{y}[n] = K x[n] \leftrightarrow X_f[k] \Rightarrow x[n] \leftrightarrow \frac{1}{K} X_f[k] = (1-a^N) X_f[k]$$

This last statment is concluded as follows: given an N-point finite length truncation $x_N[n]$ of the infinite length signal $x[n] = a^n u[n]$ for the range $0 \leq n&lt; N$, then the N-point DFT of the N-point sequence $x_N[n]$ will be $(1-a^N)X[k]$ where $X[k]$ are the samples of the true DTFT spectrum $X(e^{j\omega}) = \frac{1}{1 -a e^{-j \omega}}$ associated with the infinite length signal $x[n]=a^n u[n]$. So the error associated with the DFT computation is the fixed weight $1-a^N$ per sample of the DFT.

You should perform a similar analysis for the Gaussian signal $x[n] = \beta e^{-\alpha n^2}$ to find the error associated per computed DFT sample.