Discretized function becomes complex while free propagating a real function when sampled at even number of points using FFT and IFFT in Python

My code involves propagation of a real function using Fourier transform and inverse Fourier transform. Specifically, the function evolves as

∂ψ(z,t)/∂t - v ∂ψ(z,t)/∂z =0

I solve this problem by Fourier transforming the above equation given by

∂ψ(k,t)/∂t + ikv ψ(k,t)=0

which has a solution given by

ψ(k,t)=e^(-ikvt)ψ(k,0)

Now performing inverse Fourier transform of this will give me the propagated function

ψ(z-vt,0)

However, when I do this in Python (see below for code), the final function (after applying FFT and IFFT) seems to have a tiny complex part, particularly when the function is sampled at even number of points. This error builds up over the course of my simulation resulting in less accurate answer. However, this issue seems to go away when I sample the function at odd number of points. Can somebody please help with what is happening here?

Below is a simple code that illustrates my point using the norm of the function. Upon running this code, we can see that when the function is sampled at 10 points, the norm of the real part of the function after doing FFT and IFFT is accurate sometimes only up to 2 or 3 decimal points (while the norm of the whole function (real+complex) is still conserved). On the other hand, the norm is accurate up to 16 decimals in the case of odd sampling. It would be great if someone could explain to me what's happening and how would I get better accuracy using this function while sampling at even number of points

import numpy as np

def free_prop(vectr,Nd,vel=0.92,dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2.0*np.pi*np.fft.fftfreq(Nd, d=dz)
    return np.fft.ifft(np.exp(-1j*dt*k_vals*vel)*psi_k)

vectr_even=np.random.rand(10) ## Even case
vectr_odd=np.random.rand(11) ## Odd case

print('Norm of input array (with even sampling):',np.linalg.norm(vectr_even))
print('Norm of output complex array (with even sampling):',np.linalg.norm(free_prop(vectr_even,10)))
print('Norm of real part of output array (with even sampling):',np.linalg.norm(np.real(free_prop(vectr_even,10))))
print()
print('Norm of input array (with odd sampling):', np.linalg.norm(vectr_odd))
print('Norm of output complex array (with odd sampling):',np.linalg.norm(free_prop(vectr_odd,11)))
print('Norm of real part of complex array (with odd sampling):',np.linalg.norm(np.real(free_prop(vectr_odd,11))))

As can be seen from here, the norm in the even case (when only real part is considered) differs at the second decimal, and I would be very surprised if it a floating point error.

Look at k_vals:

With Nd=6:

[0., 6.28318531, 12.56637061, -18.84955592, -12.56637061, -6.28318531]

With Nd=5:

[0., 6.28318531, 12.56637061, -12.56637061, -6.28318531]

When Nd is even, there's one negative frequency that doesn't have a positive frequency counterpart. As we know, the input to the IFFT must be conjugate symmetric for the output to be real-valued. With an odd-length array, the computed exponential will automatically be conjugate symmetric. With the even-length array, there's the one frequency that doesn't have a positive counterpart. We need to make this frequency component real-valued for the exponential to be conjugate symmetric, and the inverse FFT to be real-valued.

For example, this modification to your function accomplishes this:

def free_prop(vectr, Nd, vel=0.92, dt=1):
    dz=1/Nd
    psi_k = np.fft.fft(vectr)
    k_vals = 2 * np.pi * np.fft.fftfreq(Nd, d=dz)
    k_kernel = np.exp(-1j * dt * k_vals * vel)
    if np.mod(Nd, 2) == 0:
       k_kernel[Nd // 2] = k_kernel[Nd // 2].real
    return np.fft.ifft(k_kernel * psi_k)

Whether this still matches the physics I don't know. One could also argue that that frequency bin should be 0, because it contains aliased data by definition (it's the Nyquist frequency, and to properly sample signal it must has 0 energy at and above the Nyquist frequency).