Introduction to the Bethe ansatz II

Building on the fundamentals introduced in part I, we employ the Bethe ansatz to study some ground-state properties (energy, magnetization, susceptibility) of the one-dimensional s=1/2 Heisenberg antiferromagnet in zero and nonzero magnetic field. The 2-spinon triplet and singlet excitations from the zero-field ground state are discussed in detail, and their energies are calculated for finite and infinite chains. Procedures for the numerical calculation of real and complex solutions of the Bethe ansatz equations are discussed and applied. The paper is designed as a tutorial for beginning graduate students. It includes 10 problems for further study.


I. INTRODUCTION
Quantum spin chains are physically realized in quasione-dimensional (1D) magnetic insulators. In such materials, the magnetic ions (for example, Cu ++ and Co ++ with effective spin s = 1 2 , Ni ++ with s = 1, and Mn ++ with s = 5 2 ) are arranged in coupled linear arrays which are well separated and magnetically screened from each other by large non-magnetic molecules. Over the years, magneto-chemists have refined the art of designing and growing crystals of quasi-1D magnetic materials to the point where physical realizations of many a theorist's pet model can now be custom-made.
Most prominent among such systems are realizations of the 1D spin-1/2 Heisenberg antiferromagnet, the model system which inspired Bethe to formulate the now celebrated method for calculating eigenvalues, eigenvectors, and a host of physical properties. Interest in quasi-1D magnetic materials has stimulated theoretical work on quantum spin chains from the early sixties until the present. Some of the advances achieved via the Bethe ansatz emerged in direct response to experimental data which had remained unexplained by standard approximations used in many-body theory.
In Part I of this series 1 we introduced the Bethe ansatz in the spirit of Bethe's original 1931 paper. 2 The focus there was on the Hamiltonian H F = −H A , where the negative sign makes the exchange coupling ferromagnetic. The ferromagnetic ground state has all the spins aligned and is (N +1)-fold degenerate. The reduced rotational symmetry of every ground state vector relative to that of the Hamiltonian reflects the presence of ferromagnetic long-range order at zero temperature. One vector of the ferromagnetic ground state, |F ≡ |↑ · · · ↑ , in the notation of Part I, serves as the reference state (vacuum) of the Bethe ansatz. All other eigenstates are generated from |F through multiple magnon excitations.
We also investigated some low-lying excitations which involve only a small number of magnons. For example, the Bethe ansatz enabled us to study the properties of a branch of two-magnon bound-states and a continuum of two-magnon scattering states, and gave us a perfect tool for understanding the composite nature of these states in relation to the elementary particles (magnons).
In this column we turn our attention to the Heisenberg antiferromagnet H A . All the eigenvectors remain the same as in H F , but the energy eigenvalues have the opposite sign. Therefore, the physical properties are very different, and the state |F now has the highest energy. Our immediate goals are to find the exact ground state |A of H A , to investigate how the state |A gradually transforms to the state |F in the presence of a magnetic field of increasing strength, and to prepare the ground work for a systematic study of the excitation spectrum relative to the state |A .
As in Part I we will emphasize computational aspects of the Bethe ansatz. The numerical study of finite systems via the Bethe ansatz is akin to a simulation in many respects, and yields important insights into the underlying physics.

II. GROUND STATE
What is the nature of the ground state state |A ? How do we find it? Is its structure as simple as that of |F ? An obvious candidate for |A is the Néel state. The two vectors reflect antiferromagnetic long-range order in its purest form just as the vector |F does for ferromagnetic longrange order. Henceforth we assume that the number of spins N in (1) is even and that periodic boundary conditions are imposed. Inspection shows that neither |N 1 , |N 2 , nor the translationally invariant linear combinations, |N ± = (|N 1 ± |N 2 )/ √ 2, are eigenvectors of H A . In the energy expectation value H A , the Néel state minimizes S z n S z n+1 but not S x n S x n+1 and S y n S y n+1 (Problem 1). A state with the full rotational symmetry of H A can have a lower energy. Like |N ± , the ground state |A will be found in the subspace with S z T ≡ n S z n = 0. Because of their simplicity, the Néel states are convenient starting vectors for the computation of the finite-N ground state energy and wave function via standard iterative procedures (steepest-descent and conjugate-gradient methods). 3 However, here we take a different route.
In the framework of the Bethe ansatz, all eigenstates of H A with S z T = 0 can be obtained from the reference state |F by exciting r ≡ N/2 magnons with momenta k i and (negative) energies −J(1−cos k i ). The exact prescription was stated in Part I. Each eigenstate is specified by a different set of N/2 Bethe quantum numbers {λ i }. The momenta k i and the phase angles θ ij in the coefficients (I28) of the Bethe wave function (I27) result from the Bethe ansatz equations (I33) and (I35). 4 A finite-N study indicates that the ground state |A has real momenta k i and Bethe quantum numbers (Problem 2) In Part I we worked directly with k i and θ ij , which represent physical properties of the elementary particles (magnons) created from the vacuum |F . At the opposite end of the spectrum, the state |A , generated from |F via multiple magnon excitations, will be configured as a new physical vacuum for H A . The entire spectrum of H A can then be explored through multiple excitations of a different kind of elementary particle called the spinon.
Computational convenience suggests that we replace the two sets of variables {k i } and {θ ij } by a single set of (generally complex) variables {z i }. If we relate every magnon momentum k i to a new variable z i by via the function φ(z) ≡ 2 arctan z, then the relation (I33) between every phase angle θ ij for a magnon pair and the difference z i − z j involves φ(z) again: Here ℜ(x) denotes the real part of x, and sgn(y) = ±1 denotes the sign of y. Relations (4) and (5) substituted into (I35) yield the Bethe ansatz equations for the variables z i : 5 The new Bethe quantum numbers I i assume integer values for odd r and half-integer values for even r. The relation between the sets {λ i } and {I i } is subtle. It depends, via sgn[ℜ(z i − z j )], on the configuration of the solution {z i } in the complex plane. For the ground state |A , we obtain (Problem 3a) Given the solution {z 1 , . . . , z r } of Eqs. (6) for a state specified by {I 1 , . . . , I r }, its energy and wave number are (Problem 3b) with E F = JN/4. The quantity φ(z i ) is called the magnon bare momentum, and ε(z i ) = dk i /dz i = −2/(1 + z 2 i ) is the magnon bare energy. The Bethe wave function (I27) is obtained from the {z i } via (4) and (5).
The ground state |A belongs to a class of eigenstates which are characterized by real solutions {z i } of the Bethe ansatz equations. To find them numerically we convert Eqs. (6) into an iterative process: Starting from z (0) i = 0, the first iteration yields z (1) i = tan(πI i /N ). Convergence toward the roots of (6) is fast. High-precision solutions {z i } can be obtained on a personal computer within seconds for systems with up to N = 256 sites and within minutes for much larger systems (Problem 4). The ground-state energy per site inferred from (8a) is listed in Table I for several lattice sizes. In preparation of the analytical calculation which produces (E A − E F )/JN for N → ∞, we inspect the z iconfiguration for the finite-N ground state |A obtained numerically. All roots are real, and their values are sorted in order of the associated Bethe quantum numbers I i . The 16 circles in Fig. 1(a) show z i plotted versus I i /N for N = 32. The solid line connects the corresponding data for N = 2048. The line-up of the finite-N data along a smooth monotonic curve suggests that the solutions of (6) can be described by a continuous distribution of roots for N → ∞.
We give the inverse of the discrete function shown in Fig. 1(a) the name Z N (z i ) ≡ I i /N and rewrite Eqs. (6) in the form: For N → ∞, Z N (z i ) becomes a continuous function Z(z) whose derivative, σ 0 (z) ≡ dZ/dz, represents the distribution of roots. In Eq. (10) the sum (1/N ) j . . . is replaced by the integral dz ′ σ 0 (z ′ ) . . . Upon differentiation the continuous version of (10) becomes the linear integral equation Applying the inverse Fourier transform to its solution yields the result The (asymptotic) ground-state energy per site as inferred from (8a), is significantly lower than the energy expectation value of the Néel states (Problem 1). The state |A has total spin S T = 0 (singlet). Unlike |F , it retains the full rotational symmetry of (1) and does not exhibit magnetic long-range order. 6 The wave number of |A , obtained from (8b), is k A = 0 for even N/2 and k A = π for odd N/2. The important result (13) was derived by Hulthén 7 in the early years of the Bethe ansatz. In Fig. 1(b) we plot the magnon momenta k i of |A as inferred from (4) versus I i /N for the same data as used in Fig. 1(a). The smoothness of the curve reflects the fact that the state |A for N → ∞ can be described by a continuous k i -distribution (Problem 5)

III. MAGNETIC FIELD
In the presence of a magnetic field h, the Hamiltonian (1) must be supplemented by a Zeeman energy: The two parts of H are in competition. Spin alignment in the positive z-direction is energetically favored by the Zeeman term, but any aligned nearest-neighbor pair costs exchange energy. Given that [H A , S z T ] = 0, the eigenvectors are independent of h. The 2S T +1 components (with |S z T | ≤ S T ) of any S T -multiplet fan out symmetrically about the original level position and depend linearly on h.
The largest downward energy shift in each multiplet is experienced by the state with S z T = S T , and that shift is proportional to S T . The state |A , which has S T = 0, does not move at all, whereas the state |F with S T = N/2 descends more rapidly than any other state. Even though |F starts out at the top of the spectrum, it is certain to become the ground state in a sufficiently strong field. The saturation field h S marks the value of h where |F overtakes its closest competitor in the race of levels down the energy axis.
The pattern in which levels with increasing S z T become the ground state of H as h increases depends on their relative starting position along the energy axis. From the zero-field energies of this set of states, we will now determine the magnetization m z ≡ S z T /N of the ground state as a function of h. 8 The Bethe quantum numbers of the lowest state with quantum number S z T = N/2 − r ≥ 0 are 5 as can be confirmed by finite-N studies of all states in the invariant subspaces with r = 1, . . . , N/2. Using the iterative process (9), we can determine the energies of these states with high precision (Problem 6). The red circles in Fig. 2(a) represent the quantity [E(S z T ) − E F ]/JN for N = 32. The solid line connects the corresponding results for N = 2048. For the finite-N analysis it is important to note that both the level energies E(S z T ) and the level spacings At h = 0, all of these levels experience a downward shift of magnitude hS z T . All spacings between adjacent levels shrink by the same amount h. The first level crossing occurs between the state |A with S z T = 0 and the state with S z T = 1, which thereby becomes the new ground state. Next, this state is overtaken and replaced as the ground state by the state with S z T = 2 and so forth. The last of exactly N/2 replacements involves the state with S z T = N/2 − 1 and the state |F with S z T = N/2. Their energy difference in zero field is 2J independent of N (see Part I). Consequently, the saturation field is h S = 2J. The magnetization m z grows in N/2 steps of width 1/N between h = 0 and h = h S . In Fig. 2(b) we plot m z versus h for two system sizes based on data obtained numerically. The blue staircase represents the results for N = 32. For N = 2048 the step size has shrunk to well within the thickness of the smooth curve shown.
An important observation is that the midpoints of the vertical and horizontal steps of the m z (h) staircase (red dots) lie very close to the limiting curve. This behavior made it possible to extract quite accurate magnetization curves for various spin chain models from very limited data. 9 Different scenarios are conceivable. For example, if the levels were arranged in the same sequence as in Fig. 2(a) but with spacings increasing from top to bottom, then the state with S z T = 0 would be replaced directly by the state with S z T = N/2. In the thermodynamic limit, the energy per site of the lowest level with given S z T becomes the internal energy density at zero temperature, From (17) we obtain, via the thermodynamic relations, The function m z (h) shown in Fig. 2(b) is the inverse of the derivative of the function u(m z ) plotted in Fig. 2(a). The slope of m z (h) represents the longitudinal susceptibility χ zz . In a finite system, where m z = S z T /N varies in steps of size 1/N , Eqs. (18) are replaced by .
The data in Fig. 2(b) indicate that χ zz (h) has a nonzero value at h = 0, grows monotonically with h, and finally diverges at the saturation field h = h S . The initial value, 8,5 turns out to be elusive to a straightforward slope analysis because of a logarithmic singularity which produces an infinite curvature in m z (h) at h = 0 (Problem 7). The divergence of χ zz (h) at h S is of the type (Problem 8) .
The characteristic upwardly bent magnetization curve with infinite slope at the saturation field is a quantum effect unreproducible by any simple and meaningful classical model system. The Hamiltonian (1), reinterpreted as the energy function for coupled three-component vectors, predicts a function m z (h) which increases linearly from zero all the way to the saturation field.

IV. TWO−SPINON EXCITATIONS
Returning to zero magnetic field, let us explore the spectrum of the low-lying excitations. From here on, the ground state |A (with S z T = S T = 0) replaces |F (with S z T = S T = N/2) as the new reference state for all excited states. The Bethe quantum numbers (7), which characterize |A , describe a perfectly regular array on the I-axis as illustrated in the first row of Fig. 3. This array will be interpreted as a physical vacuum. The spectrum of H A can then be generated systematically in terms of the fundamental excitations characterized by elementary modifications of this vacuum array. In the subspace with S z T = 1, a two-parameter set of states is obtained by removing one magnon from the state |A . In doing so we eliminate one of the N/2 Bethe quantum numbers from the set in the first row of Fig. 3 and rearrange the remaining I i in all configurations over the expanded range − 1 4 N ≤ I i ≤ 1 4 N . Changing S z T by one means that the I i switch from half-integers to integers or vice versa. The number of distinct configurations with I i+1 − I i ≥ 1 is 1 8 N (N + 2). A generic configuration consists of three clusters with two gaps between them as shown in the second row of Fig. 3.
The two gaps in the otherwise uniform distribution of Bethe quantum numbers can be interpreted as elementary particles (spinons) excited from the new physical vacuum. The position of the gaps between the I iclusters determine the momentak 1 ,k 2 of the two spinons, which, in turn, add up to the wave number of the twospinon state relative to the wave number of the vacuum: What remains to be done for the finite-N analysis is straightforward. We solve the Bethe ansatz equations via (9) for all the I i -configurations just specified. A plot of the two-spinon energies E − E A versus wave number k − k A for N = 16 as inferred from the solutions {z i } via (8a) and (8b) is shown in Fig. 4 (red circles). Also shown are the Bethe quantum numbers for each excitation. The pattern is readily recognized and extended (via reflection I i ↔ −I i ) to the other half of the Brillouin zone. We have drawn (in blue) the corresponding two-spinon excitations for N = 256. The 1 4 N ( 1 4 N + 1) = 4160 dots in the range 0 < q ≤ π produce a density plot for the two-spinon continuum which emerges in the limit N → ∞. Compare this set of two-spinon scattering states with the set of two-magnon states plotted in Fig. 2 of Part I. There we found that pairs of magnons form a continuum of two-magnon scattering states and a branch of twomagnon bound states. Two-spinon bound states exist as well and will be discussed in a later column.
A derivation of the exact lower and upper boundaries of the two-spinon continuum, 10 starts from the Bethe ansatz equations (10). We set r = N/2 − 1 and use the Bethe quantum numbers from the second row of Fig. 3. When we replace the sum by an integral, we must account for the gaps between the I i -clusters by two 1/N -corrections. The result (after differentiation) is an integral equation which differs from (11) by two extra terms related to the I i -gaps: If we write σ(z) = σ 0 (z) + σ 1 (z) + σ 2 (z), where σ 0 (z) is the solution (12) of Eq. (11), the two-spinon corrections are solutions of Because Eqs. (24) have the same integral kernel K as Eq. (11), the solutions of all three equations can be expressed by the same resolvent R, where g 0 (z) and g 1 (z), g 2 (z) are the inhomogeneities of Eqs. (11) and (24), respectively. The resolvent, which does not depend on the inhomogeneity, satisfies Because g l (z) = (1/N )K(z −z l ) for l = 1, 2 makes (24) equivalent to (26), we can write σ l (z) = (1/N )R(z −z l ). This is all we need to know about the solutions. Using Eq. (8a) for the two-spinon energies, we must again correct for the two I i -gaps when we convert the sum into an integral: Now we must relate the I i -gaps, that is, the valuesz 1 ,z 2 to the spinon momentak 1 ,k 2 . Starting from the configuration in the first row of Fig. 3, we remove one Bethe quantum number at or near the center. The accompanying integer ↔ half-integer switch of the remaining I i opens a double gap at or near the center. Shifting several I i from the left and right toward the center then produces the generic three clusters. This change in configuration applied to (8b) yields the following two-spinon wave number relative to the vacuum: or sink l = sech(πz l /2), l = 1, 2. These relations substituted into Eq. (28) produce the two-spinon spectrum = πJ| sin(q/2) cos(p/2)| (30) for 0 ≤ p ≤ q and 0 ≤ q ≤ 2π − p. The spectrum is a two-parameter continuum with boundaries (22). Like magnons, spinons carry a spin in addition to energy and momentum. 11 Unlike magnons, which are spin-1 particles (see Part I), spinons are spin-1 2 particles. In a chain with even N , where all eigenstates have integer-valued S z T , spinons occur only in pairs. The spins s l = ± 1 2 , l = 1, 2 of the two spinons in a two-spinon eigenstate of H A can be combined in four different ways to form a triplet state (S T = 1, S z T = 0, ±1) or a singlet state (S T = S z T = 0). They are described by distinct configurations of Bethe quantum numbers.
We have already analyzed the spectrum of the twospinon triplet states with s 1 = s 2 = + 1 2 (S T = 1, S z T = +1), as specified by I i -configurations of the kind shown in the second row of Fig. 3. The two-spinon states with s 1 = s 2 = −1 (S T = 1, S z T = −1) are obtained from these states by a simple spin-flip transformation. They have exactly the same energy-momentum relations.
The remaining two sets of two-spinon states have s 1 = −s 2 , that is, S z T = 0. One of these sets contains triplets (S T = 1) and the other singlets (S T = 0). The former is obtained if we shift all Bethe quantum numbers in the second row one position to the right (I i → I i + 1 2 , i = 1, . . . , r − 1) as shown in the third row and add one nonmovable Bethe quantum number I r = 1 4 N + 1 2 to the set (Problem 9). The number of distinct configurations is then the same as in row two and the corresponding states have the same wave number and energy. Symmetry requires that all three components (S z T = 0, ±1) of the triplet states are degenerate.
Whereas all two-spinon triplets are described by real solutions {z i }, the two-spinon singlet states (S z T = S T = 0) are characterized by one pair of complex conjugate solutions z 1 ≡ u + iv, z 2 ≡ u − iv in addition to the real solutions z 3 , . . . , z N/2 . 12 Finding the I i -configurations for a particular set of eigenstates with complex solutions (such as indicated in row four of Fig. 3) and then solving the associated Bethe ansatz equations turns out to be a much more delicate task than was the identification and the calculation of purely real solutions.
A straightforward adaptation of the iterative process (9) to complex z i will, in general, not converge toward a physically meaningful solution. Making the Bethe ansatz equations amenable to root-finding algorithms which are adequate for this task requires that we convert Eqs. (6) into a set of equations with real solutions. For the twospinon singlets we arrive at the following set of N/2 equations for {u, v, z 3 , . . . , z N/2 } (Problem 10): where ϕ(x) ≡ 2 atanh(x). Here we have replaced the first two Bethe quantum numbers I 1 , I 2 of Eq. (6) by a single (integer) Bethe quantum number I (2) associated with z 1 = z * 2 . The Bethe quantum numbers associated with real solutions, now renamed I (1) i , are constrained to the range (− 1 4 N − 1 2 ≤ I (1) i ≤ 1 4 N + 1 2 ) as indicated in the fourth row of Fig. 3. This notation is used in preparation of a general classification of Bethe ansatz solutions (string hypothesis) which is valid for N → ∞ and will be discussed in a later column.
The dependence of the wave number on the Bethe quantum numbers for all states with one pair of complex solutions is then where r = N/2 for the two-spinon singlets considered here. Figure 5 depicts the energy (8a) versus the wave number (32) of the two-spinon singlets (red circles) for N = 16 and 0 < q ≤ π. We use the same representation as in Fig. 4 and show the two-spinon triplets again for comparison. The number of two-spinon singlets identified here is smaller, namely 1 8 N (N −2) over the entire Brillouin zone, than the number 1 8 N (N + 2) of two-spinon triplets identified previously. Although the I i -configurations of the singlets are more complicated than those of the triplets, the pattern is still recognizable to the extent that it may serve as a useful guide for applications to other system sizes. The reader will find it easy to identify in the tables of Part I one state for N = 4 and three states for N = 6 which belong to this class of singlet excitations. In Table  II we present the two-spinon singlet data for N = 8. Note that the pattern of {I 3 , . . . , I N/2 } is akin to the pattern of {I 1 , . . . , I N/2 } for the triplets of a shorter chain. As N grows larger, the effect of the complex solutions z 1 = z * 2 on the energy (8a) relative to that of the real solutions z 3 , . . . , z N/2 diminishes and disappears in the limit N → ∞. The two-spinon singlets will then form a continuum with the same boundaries (22). 12 Yet the effect of the complex solutions in the two-spinon singlets will remain strong for quantities inferred from the Bethe ansatz wave function (for example, selection rules and transition rates).
How do we solve Eqs. (31)? At q = π the real roots z 3 , . . . , z N/2 of the two-spinon singlets are symmetrically distributed about z = 0. This observation paves the way for finding the complex roots exactly. Equation (31c) admits a solution in the limit |u| → 0 + , |v| → 1 + along a path with |u| ∝ (|v| − 1) 1/N . This limit matches the divergence on the left with that in the first term on the right. The z i in the non-singular terms on the right have no effect on this solution.
Significant computational challenges arise in the determination of two-spinon singlet states at q = π, where the real roots are no longer symmetrically distributed and the complex roots have u > 0, v > 1. Now there is no way around solving Eqs. (31) simultaneously. For this task a good root-finding algorithm 13 combined with carefully chosen starting values will be needed.
The two-spinon triplets play an important role in the zero-temperature spin dynamics of quasi-1D antiferromagnetic compounds. They are the elementary excitations of (1) which can be directly probed via inelastic neutron scattering. The two-spinon singlets, in contrast, cannot be excited directly from |A by neutrons because of selection rules. The singlet excitations are important nevertheless, but in a different context.
Several of the quasi-1D antiferromagnetic compounds are susceptible to a spin-Peierls transition, which involves a lattice distortion accompanied by an exchange dimerization. The operator which probes the dimer fluctuations in the ground state of (1) couples primarily to the two-spinon singlets and not at all to the two-spinon triplets. In a forthcoming column of this series, the focus will be on transition rates between the ground state of (1) in zero and nonzero magnetic field and several classes of excited states that are important for one reason or another. 3. (a) Use the derivation of (6) from (I33) and (I35) to establish the relation between the λ i and the I i . Given that all λ i are integers, show that all I i are integers for odd r and half-integers for even r. For the given (real) solution {z 1 < · · · < z N/2 } determine the one-to-one correspondence between the elements of (3) and (7) Table I. Extrapolate your finite-N data for (E A − E F )/JN and compare your best prediction and its uncertainty with the exact result (13). (4) and (12) to calculate the distribution of the magnon momenta in the ground state |A via ρ 0 (k) = ∞ −∞ dz σ 0 (z)δ(k − π − φ(z)). (b) Express the energy per site of |A as an integral of the magnon energy, −J(1 − cos k), weighted by ρ 0 (k), and evaluate it to confirm the result (13). 6. In Problem 5 we interpreted |A as a composite of N/2 magnons with momentum distribution (14). Increasing m z from zero means reducing the number of magnons in the ground state to N ( 1 2 − m z ). From the solution {z i }, calculate the finite-N distribution ρ N (k i ) = [N (k i −k i+1 )] −1 of magnon momenta for m z = 0, 0.125, 0.25, 0.375. Use data for N = 32 and for a much larger system, for example, N = 2048. Show that for h = 0, ρ N (k i ) converges very rapidly toward (14). Show that increasing h does not only deplete magnons, it also constricts the range of the magnon momenta.

(a) Use
7. Generate data for E(S z T ) via (9). Calculate h(m z ) from (19a) and χ zz (m z ) from (19b). Then plot χ zz versus h for 0 < h < 0.25 in relation to the exact value 1/(π 2 J) at h = 0. Show that the slope of χ zz (h) increases without bounds as h → 0. Determine the χ zz at h = 0 via extrapolation of the h = 0 data points, and compare your prediction with the exact result (20).