Bethe ansatz: general comments, Lieb-Liniger and Heisenberg models

Remarkably, there exist exact solutions for the ground state and thermodynamics of a number of one-dimensional problems of quantum-mechanical interacting particles. The idea goes back to Bethe, who constructed an exact solution of the ground state of Heisenberg antiferromagnet in one dimension, but its full mathematical significance was understood only much later. The quantum "integrable models" for which wavefunctions of Bethe type are eigenstates are related to the solvability of certain two-dimensional classical statistical models such as the Ising model, famously solved by Onsager.

Here we consider two examples. The simplest model solvable by Bethe's approach is probably the gas of nonrelativistic bosons with delta-function interactions. This is often called the Lieb-Liniger model, as those authors found its ground state; soon afterward Yang and Yang found its thermodynamics and introduced the "thermodynamical Bethe Ansatz" in a short paper that is highly recommended. The second is the Heisenberg model, where we will present more details.

(Lieb-Liniger was covered in class, following book of Takahashi; will update here later)

We present Bethe's solution of the Heisenberg antiferromagnetic spin-half chain. So as not to be carrying around factors of $J$ everywhere, we take the Hamiltonian in this section to be
$$\begin{equation} \label{heisbethe} H = \frac{1}{2} \sum_{j=1}^N ({\bf \sigma}_j \cdot {\bf \sigma}_{j+1} - 1). \end{equation}$$
This corresponds to $J = 2$ in (\ref{heisenberg}), plus an energy offset. The energy offset means that the ferromagnetic state (all spins up, for example), will have energy 0. A lower bound for the ground state energy can be obtained by imagining that every pair of interacting spins could form a singlet simultaneously, which would then give energy $-2$ per bond. Based on the XX model above, it is useful to divide this into one term that hops a spin-flip, and another term that is an interaction between spin-flips, with a total number $M$ of spin-flips:
$$\begin{equation} H = \sum_{j=1}^N \left[- \sigma^+_j \sigma^-_{j+1} - \sigma^-_j \sigma^+_{j+1} + 2 n_j n_{j+1} \right] - 2 M. \end{equation}$$
(The $-$ sign appears because, following convention, we have performed the canonical transformation $s^x \rightarrow -s^x$, $s^y \rightarrow -s^y$, $s^z \rightarrow s^z$ on every other site. This would change $k$ to $k+\pi$ in our solution of the XX model above but leave all energies and other physical quantities invariant.) Here we have introduced $n_j$ as the "number of spin flips on site $j$", so that the last term in the sum looks like an interaction between particles. The $-2 M$ term is required because of how we have rewritten the $\sigma_z$ term:
$$\begin{equation} \frac{1}{2} \sum_j \sigma^z_j \sigma^z_{j+1} = \frac{1}{2} \sum_j (2 n_j-1)(2 n_{j+1}-1) = -2 M + N/2 + \sum_j 4 n_j n_{j+1} \end{equation}$$
and the constant term eliminates the constant term in (\ref{heisbethe}).

As in the XX model above, we will build up the antiferromagnetic ground state as a system of $M=N/2$ delocalized spin-flips relative to the ferromagnetic state, but now these spin-flips will interact through interactions rather than just through their Fermi statistics. Our trial wavefunction, the "Bethe ansatz", is
$$\begin{equation} \Psi(x_1,x_2,\ldots,x_M) = \sum_P A(P) \exp\left[i \sum_{j=1}^M k_{Pj} x_j\right]. \end{equation}$$
Here the $x_j$ are sites on the one-dimensional lattice. The meaning of this wavefunction in terms of particles is that it gives us the location of the $M$ spin-flips relative to the ferromagnet. It remains to determine the amplitude $A(P)$ for a given permutation $P$ and also the momenta $k_j$.

The amplitude will reflect the nature of collisions in one dimension: since elastic collisions can not change the momenta of the two particles, the collision's effects are entirely captured by a phase shift. We define this as the phase $\theta(k,k^\prime)$ that appears in the solution of the two-particle problem
$$\begin{equation} \label{twobody} \Psi(x_1,x_2) = e^{i ( k x_1 + k^\prime x_2)} - e^{i (k^\prime x_1 + k x_2) - i \theta(k,k^\prime)}. \end{equation}$$
Based on this, we define $A(P)$ so that, if the only change between permutations $P$ and $P^\prime$ is that the two momenta $k$ and $k_\prime$ are interchanged,
$$\begin{equation} A(P^\prime) / A(P) = - e^{-i \theta(k,k^\prime)}. \end{equation}$$
We will need to compute this phase shift $\theta$ later. But first note that the energy is determined completely by the momentum distribution as long as there is some configuration where no adjacent sites are occupied, as is the case for $M \leq N/2$. Then all the energy is kinetic, and
$$\begin{equation} E = 2 \sum_{j=1}^M (-1 - \cos k_j). \label{partenergy} \end{equation}$$
This is similar in form to the $XX$ result (with the additional offset, and noting that here $M$ is the number of particles, which is half the number of bonds), but the momentum distribution will be different. However, it remains to prove that there is an eigenstate with this energy and to determine the $k$ distribution.

The phase shift $\theta$ describes the scattering of two particles that cannot sit on the same site (since any spin can be flipped at most once) and, if on adjacent sites, have a repulsive interaction energy 2. We will go a step farther and show that our many-particle Bethe ansatz wavefunction is an energy eigenstate, whose energy must then be given by (\ref{partenergy}); this is clear at least for $M \leq N/2$, where it is possible to have an initial state with no particles on adjacent sites. Consider the wavefunction $\Psi(x_1,\ldots,x_M)$ defined above, and allow it to be defined even if any two $x_i$ are identical (call this $\Psi_g$). The next step is somewhat subtle: if we did not require the wavefunction to vanish when any two $x_i$ were identical, then the above kinetic energy (\ref{partenergy}) would apply and
$$\begin{equation} H_{nonint} \Psi_g = E \Psi_g = -2 M \Psi -\sum_{j=1}^M \left[\Psi_g(\ldots,x_j+1,\ldots) + \Psi(\ldots,x_j-1,\ldots)\right], \end{equation}$$
and the state would be an eigenstate of the Hamiltonian without the "interaction energy" term. We will use the interaction energy term to balance the fact that the wavefunction must in reality vanish when two particles are on the same site. In other words, we need that for a wavefunction with one nearest-neighbor pair,
$$\begin{equation} H \Psi = E \Psi \Rightarrow H_{nonint} \Psi + H_{int} \Psi = E \Psi \Rightarrow 2 \Psi + \Psi_g(\ldots,x,x,\ldots) + \Psi_g(\ldots,x+1,x+1,\ldots) = 0. \end{equation}$$
where the second two terms on the right side result because in reality there is no term contributing to the nearest-neighbor $\Psi$ that came from having two particles on the same site. A similar equation, just with more terms, results if there is more than one nearest-neighbor pair.

If we look at our form for the wavefunction above, there are $M!$ terms. Given two particles $x_a, x_b$, these terms can be grouped into pairs where each pair assigns the same two momenta to the same two particles; the terms within a pair differ by an interchange of which particle gets which momentum. So each pair has the two-particle eigenstate $\Psi(x_a,x_b)$ as a common factor; in this sense the Bethe ansatz wavefunction is built up from two-particle wavefunctions. From the above, the many-particle wavefunction will be an eigenstate if the two-particle wavefunction satisfies
$$\begin{equation} \Psi(x,x) + \Psi(x+1,x+1) + 2 \Psi(x,x+1)=0. \end{equation}$$
Inserting the two-body wavefunction from (\ref{twobody}), we get
$$\begin{equation} 1 + e^{i k + i k^\prime} + 2 e^{i k^\prime} = e^{-i \theta(k,k^\prime)} \left[1 + e^{i k + i k^\prime} + 2 e^{ik} \right] \end{equation}$$
so
$$\begin{equation} - i \theta(k,k^\prime) = \log\left[{1 + e^{i k + i k^\prime} + 2 e^{i k^\prime} \over 1 + e^{i k + i k^\prime} + 2 e^{ik} } \right] \end{equation}$$
Since $\arctan z = {i \over 2} \log [(1-iz)/(1+iz)]$, we can rewrite this if desired as
$$\begin{equation} \theta(k,k^\prime) = i \log\left[{\cos((k+k^\prime)/2) + 2 e^{i (k^\prime-k)/2} \over \cos((k+k^\prime)/2) + 2 e^{i(k-k^\prime)/2}}\right] = -2 \arctan\left[\sin((k^\prime-k)/2) \over \cos((k+k^\prime)/2) + \cos((k^\prime-k)/2) \right] \end{equation}$$
The phase shift takes a simpler form if we make the change of variables
$$\begin{equation} e^{ik} = {1 + 2 i \alpha \over 1 - 2 i \alpha} \Leftrightarrow k(\alpha) = 2 \arctan (2 \alpha), \end{equation}$$
so
$$\begin{equation} \theta(k,k^\prime) = i \log \left[{(1-2 i \alpha)(1-2 i \alpha^\prime) + (1 + 2 i \alpha)(1 + 2 i \alpha^\prime) + 2 (1 + 2 i \alpha^\prime)(1 - 2 i \alpha) \over (1-2 i \alpha)(1-2 i \alpha^\prime) + (1 + 2 i \alpha)(1 + 2 i \alpha^\prime) + 2 (1 + 2 i \alpha) (1 - 2 i \alpha^\prime)} \right] = i \log \left[{4-4 i (\alpha-\alpha^\prime) \over 4 +4 i (\alpha-\alpha^\prime)}\right], \end{equation}$$
or simply
$$\begin{equation} \theta(k,k^\prime) = 2 \arctan (\alpha - \alpha^\prime). \end{equation}$$

The "fundamental equation" of the Bethe ansatz is derived from the rule the total phase shift in moving the particle around the system, a closed loop of N sites, must be a multiple of $2 \pi$. In other words,
$$\begin{equation} e^{i k N} = \prod_{k^\prime \not = k} \left[-e^{i \theta(k,k^\prime)} \right], \end{equation}$$
where the negative sign appears because of how we defined the phase shift above. Taking the logarithm of both sides, we obtain
$$\begin{equation} k N = 2 \pi I(k) + \sum_{k^\prime \not = k} \theta(k,k^\prime). \label{fundamental} \end{equation}$$
This is consistent that for free particles as in the XX model, we have momenta $k_j = 2 \pi I_j / N$, $I_j = j - (M+1)/2$. There is a subtle difference between bosons and fermions here; for an odd number $M$ of bosons, or for fermions, we get $I(k)$ integer, while for an even number of bosons, $I(k)$ becomes half-integer. See Sutherland, Beautiful Models, chapter 6, for more details. We are constructing a wavefunction for bosons in this section, because the coordinates $x_j$ label spin flips that should commute, but describing "hard-core bosons" or "fermions" in one dimension is a difference of only an overall factor $(-1)^P$.

Now we come to the key mathematical step: we make the change of variables above to make a difference kernel for the integral equation that results from the continuum limit of (\ref{fundamental}). Being able to make such a change of variables is a large part of what separates solvable 1D models from others. For future reference, we note the inverse transformation
$$\begin{equation} k(\alpha) = 2 \arctan(2 \alpha), \quad k^\prime(\alpha) = {4 \over 1+4 \alpha^2}, \end{equation}$$
and that the energy for a given momentum is
$$\begin{equation} E_k = 2(-1 - \cos k) = - 2 \theta^\prime(2 \alpha) = - k^\prime(\alpha). \end{equation}$$
The density of momenta $\rho(k)$ in the continuum limit then becomes a different function $R(\alpha)$ in terms of $\alpha$:
$$\begin{equation} R(\alpha)\,d\alpha = \rho(k)\,dk \Rightarrow R(\alpha) = k^\prime(\alpha) \rho(\alpha). \end{equation}$$
The fundamental equation in terms of $R$, and dividing by $N$, is
$$\begin{equation} k(\alpha) = 2 \pi f(\alpha) + \int_{-\infty}^{\infty} \theta(\alpha-\alpha^\prime) R(\alpha^\prime)\,d\alpha^\prime, \end{equation}$$
where $f(\alpha)$ is the number of momenta less than $\alpha$, as the integers $I(k)$ are assumed increase by one for each momenta, as in the free case. Here the boundary conditions run to $\pm \infty$ after converting the original boundary condition $k = \pm \pi$ to the $\alpha$ variables.

Taking the derivative of both sides with respect to $\alpha$, we obtain
$$\begin{equation} k^\prime(\alpha) = 2 \pi R(\alpha) + \int_{-\infty}^{\infty} \theta^\prime(\alpha-\alpha^\prime) R(\alpha^\prime)\,d\alpha^\prime. \end{equation}$$
We can solve this for $R(\alpha)$ by noting the convolution form of the integral equation. Dividing through by $2 \pi$, we have an equation of the form
$$\begin{equation} ({\bf 1} + K) R = {k^\prime \over 2 \pi}, \end{equation}$$
where ${\bf 1}$ is the identity operator and $K$ is the integral operator with kernel $\theta^\prime / (2 \pi)$. The solution of this equation is easily expressed in terms of Fourier components, since that converts the convolution integral to a product:
$$\begin{equation} {\tilde R}_0 = {1 \over 2 \pi} {{\tilde k^\prime} \over 1+{\tilde K}}. \end{equation}$$

Now we have to do a few integrals. Since $K$ is the integral operator with kernel $(1/\pi) (1 + (\alpha-\alpha^\prime)^2)^{-1}$, its Fourier transform is a simple exponential:
$$\begin{equation} {\tilde K}(s) = \int_{-\infty}^\infty {1 \over \pi (1 + \alpha^2)} e^{i s \alpha}\,d\alpha = \exp(-|s|),\quad {\tilde R}_0(s) = {\exp(-|s|/2) \over 1+{\tilde K}(s)} = {1 \over 2 \cosh(s/2)}. \end{equation}$$
where we have used also that the Fourier transform of $k^\prime / (2 \pi) = (2 / \pi) (1 + 4 \alpha^2)^{-1}$ is $\exp(-|s|/2)$.

From this expression for ${\tilde R}_0$ , we can Fourier transform back to obtain the ground state distribution of $\alpha$'s,
$$\begin{equation} R_0(\alpha) = {1 \over 2 \pi} \int_{-\infty}^\infty ds\,e^{-i s \alpha} {1\over 2 \cosh(s/2)} = {1 \over 2 \cosh(\pi \alpha)} \end{equation}$$
and from that the average energy per site:
$$\begin{equation} E_0 = \int E_k \rho(k) dk = - \int_{-\infty}^\infty \theta^\prime(2 \alpha) R_0(\alpha)\,d\alpha = - \int_{-\infty}^\infty\,k^\prime(\alpha)\,R_0(\alpha)\,d\alpha = -2 \int_{-\infty}^\infty\,{1 \over \cosh(\pi \alpha) (1 + 4 \alpha^2)} = -2 \log 2. \end{equation}$$
As expected, this is significantly higher than $-2$, reflecting the fact that a spin cannot simultaneously form singlets with both its neighbors. To compare it to the XX and Ising values, we restore $J$, remove the offset, and change to the convention in the XX model above. We then get
$$\begin{equation} E_{\rm Ising} = - {J \over 4}, \quad E_{XX} = -{J \over \pi} \approx -0.318 J,\quad E_{\rm Heisenberg} = -J {4 \log 2 - 1 \over 4} \approx -0.443 J. \end{equation}$$
With more work, one can confirm that the Heisenberg chain, like the XX one, is gapless and critical. Does this apply for a spin-1 chain as well? For the nearest-neighbor chain, it does not...