# Bethe ansatz eigenfunctions of Hubbard model

I want to calculate the exact eigenfunctions of Hubbard model by using Bethe ansatz. I have solved the periodic boundary condition equations (Lied-Wu equations) and have the values of the momentum "$$k$$"s, for $$4$$ sites and $$2$$ electrons model. Now I use the following formula for the wavefunction (attached as image):

the amplitudes $$\psi(x_1,x_2,..x_{N_e})$$ are complex valued for general $$k_1$$ and $$k_2$$. But when I diagonalize the Hubbard Hamiltonian for $$4$$ sites and $$2$$ electrons, all the amplitudes of the eigenfunctions are real valued. What am I doing wrong?