At what amplitude does spectral leakage start to be noticeable in the phase spectrum?
The issue is that you are not taking a large enough piece of your original x to get an accurate frequency domain representation.

Try setting

x = exp(-linspace(-20,20,N).^2);


If I do that, I get the bottom two plots below (the middle plots are for -10,10).  The top plot is a replication of your values.

In the three pairs below, the left plot is the semilogy plot of x and the right plot is the semilogy plot of its FFT f.





Derivation of the FFT of a windowed Gaussian

To answer your question in the comments, we need to look at how a windowed version of:
$$
x(t) = e^{-t^2/2}
$$
can be analyzed in the frequency domain.

The windowed version is:
$$
x_w(t) = w(t)e^{-t^2/2} 
$$
where $w(t) = \mbox{rect(t/T)}$ if we use a rectangular window of width T (span from $-T/2$ to $+T/2$).

The Fourier transform will be:
$$
X_w(\omega) = \int_{-T/2}^{T/2} e^{-t^2/2}  e^{j\omega t} dt\\
= 2 \int_{0}^{T/2} e^{-t^2/2} \cos(\omega t) dt\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + j \mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} )
 \right] \right |^{T/2}_0\\
= \left . \sqrt{\frac{\pi}{2}} e^{-\omega^2/2} \left[
\mbox{erf}(\frac{t - j\omega}{\sqrt{2}}) + \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}} )
 \right] \right |^{T/2}_0
$$
where this integral is used and the identity
$$
\mbox{erfi}(\frac{ \omega - j t}{\sqrt{2}} ) = \mbox{erf}(j \frac{ \omega - j t}{\sqrt{2}} ) / j = -j \mbox{erf}(\frac{ t + j\omega}{\sqrt{2}})
$$

Unfortunately, I'm not sure where that gets us. I'll leave it for now and see if I can doodle more later.



scilab code below (Sorry, I tend to use scilab rather than Octave).

// Q26236;
N = 101;
x1 = exp(-linspace(-5,5,N).^2);
f1 = fftshift(fft(x1));
x2 = exp(-linspace(-10,10,N).^2);
f2 = fftshift(fft(x2));
x3 = exp(-linspace(-20,20,N).^2);
f3 = fftshift(fft(x3));

clf
subplot(321)
plot2d("nl",abs(x1))
title('x1')
subplot(323)
plot2d("nl",abs(x2))
title('x2')
subplot(325)
plot2d("nl",abs(x3))
title('x3')

subplot(322)
plot2d("nl",abs(f1))
subplot(324)
plot2d("nl",abs(f2))
subplot(326)
plot2d("nl",abs(f3))