Why is the ground state energy of the Heisenberg XXZ Model unbounded for some values of $J$? At the moment, I'm looking at numerically studying the Heisenberg XXZ model. The Hamiltonian is given below:
$$
H=\sum_{j=1}^{N-1}\left(J S_j^z S_{j+1}^z+\frac{K}{2}(S_j^+S_{j+1}^-+S_j^-S_{j+1}^+)\right)
$$
where $J$ and $K$ are coupling parameters. I was able to calculate the matrix form of the Hamiltonian fine, but I have a question:
I recently calculated the ground eigenstate $v_0$, and from here calculated the ground state energy $\epsilon_0=\langle v_0|{\bf H}|v_0\rangle$. Note that the Hamiltonian matrix ${\bf H}$ is normalized. When I set $K=1$ and plot $\epsilon_0$ vs. $J$, I find that, for $J<-1$ and $J>1$, the ground state energy is $-1$. However, for $-1<J<1$, I get bizarre behavior - namely, the energy approaches negative infinity as $J\rightarrow 0$. Is this expected behavior, or am I doing something wrong? What's causing this?
 A: This is certainly unexpected, for as Mark Mitchison commented the $J=0$ Heisenberg model is equivalent to free fermions in one dimension. Moreover I suspect that something is amiss even before that, for the $J=-1$ is the ferromagnetic model and $J=+1$ the antifferomagnetic one, and certainly they have different ground state energies.
In fact, for $J\leq -1$ is trivial that the ground state consists of all spins aligned, the state $|\uparrow...\uparrow \rangle $. In this case is easy to verify from the Hamiltonian that the energy of such ground state is $E_F=\frac{JN}{4}$, where I assume you're working with spin one-half.
Now the Heisenberg XXZ model is exactly solvable by means of Bethe Ansatz. Unfortunately the Bethe Ansatz is not particularly nice for analytical calculations, but it is very easy to evaluate the ground state in the thermodynamical limit $N\rightarrow \infty$. I'll just quote the result, which is a straighforward generalization from the one in Giamarchi's book. 
For $-1\leq J\leq 1$ the ground state has total spin zero, and denoting by $E$ the ground state energy and $E_F$ the ferromagnetic ground state energy, we have in the thermodynamical limit that
\begin{equation}
\frac{E-E_F}{N}=-\frac{(J+1)}{2}0.693-\frac{(J-1)}{2}0.307,
\end{equation}
where the constants $0.693$ and $0.307$ are obtained by numerical integration. For long chains, say $N$ larger than $10$, the thermodynamical limit should be a good approximation, so you can use it to compare with your numerical results.
