How does a Hermitian FFT work in Numpy?
You are facing the practical problem of representing a noncausal signal in a programs such as Matlab, Octave, or Python. These programs, in their default, assume that all represented sequences begin at time-index $n=0$.
Therefore, a Hermitian symmetric sequence ( $x[n] = x^*[-n]$ for all $n$ ) ,  which begins at $n=-4$, and ends at $n=4$, will lose this symmetry when interpreted in those programs. The effect is a shift of the signal to the right by $4$ for your particular example.
Hence, what you have compouted is not the FFT of the original symmetric sequence $x[n] = e^{j n}$, but a shifted, non-symmetric, sequence $w[n] = x[n-4]$. Using the DTFT / DFT time-shift property, we see what is actually computed as :
$$ W(e^{j\omega}) = e^{-j4\omega} X(e^{j\omega}) \tag{1} $$ and
$$ W[k] = e^{-j4 \frac{2\pi}{N} k} X[k] ~~~~~ , ~~\text{ for } ~~~k=0,1,...,N-1 \tag{2}$$
Where $N$ is the length of the sequences $x[n]$ and its DFT $X[k]$.
Therefore, eventhough the DTFT / DFT of the desired sequence $x[n]$ has zero imaginary part, as an Hermitian symmetric signal would have, the practically computed $W[k]$ has nonzero imaginary part.
If you wish, you can obtain the non-zero phase samples of the desired DFT $X[k]$ of $x[n]$, from the computed DFT $W[k]$ of $w[n]$, with the following:
$$ X[k] = e^{j 4 \frac{2\pi}{N} k} W[k] ~~~~ ,~~\text{ for } k=0,1,...,N-1 \tag{3}$$
The following MATLAB code demonstrates the issue:
% Demonstrates the Hermitian symmetry and its handling in Matlab/Octave/Python
% environments...

a = -4;
b = +4;
n = [a:b];
x = exp(1j*n);      % Samples of Hermitian symmetrix input: x[k] = x[-k]*

figure,stem( imag( exp(-1j*2*pi*a*[0:(b-a)]/(b-a+1)).*fft(x) ));

The output is: