Exact two-spin dynamic structure factor of the one-dimensional s=1/2 Heisenberg-Ising antiferromagnet

The exact 2-spinon part of the dynamic spin structure factor $S_{xx}(Q,\omega)$ for the one-dimensional $s$=1/2 $XXZ$ model at $T$=0 in the antiferromagnetically ordered phase is calculated using recent advances by Jimbo and Miwa in the algebraic analysis based on (infinite-dimensional) quantum group symmetries of this model and the related vertex models. The 2-spinon excitations form a 2-parameter continuum consisting of two partly overlapping sheets in $(Q,\omega)$-space. The spectral threshold has a smooth maximum at the Brillouin zone boundary $(Q=\pi/2)$ and a smooth minimum with a gap at the zone center $(Q=0)$. The 2-spinon density of states has square-root divergences at the lower and upper continuum boundaries. For the 2-spinon transition rates, the two regimes $0 \leq Q<Q_\kappa$ (near the zone center) and $Q_\kappa<Q \leq \pi/2$ (near the zone boundary) must be distinguished, where $Q_\kappa \to 0$ in the Heisenberg limit and $Q_\kappa \to \pi/2$ in the Ising limit. The resulting 2-spinon part of $S_{xx}(Q,\omega)$ is then square-root divergent at the spectral threshold and vanishes in a square-root cusp at the upper boundary. In the regime $0<Q_\kappa \leq \pi/2$, by contrast, the 2-spinon transition rates have a smooth maximum inside the continuum and vanish linearly at either boundary. Existing perturbation studies have been unable to capture the physics of the regime $Q_\kappa<Q \leq \pi/2$. However, their line shape predictions for the regime $0 \leq Q<Q_\kappa$ are in good agreement with the new exact results if the anisotropy is very strong.


I. INTRODUCTION
Among all the spin-chain models that are directly relevant for the description of real quasi-one-dimensional (1D) magnetic insulators, the s = 1/2 XXZ model, (σ x n σ x n+1 + σ y n σ y n+1 + ∆σ z n σ z n+1 ), (1.1) is the one whose physical properties have been studied most comprehensively. Today there exist more exact results for this model than for any other model of comparable importance. The early demonstration 1-3 that the XXZ model is amenable to the Bethe ansatz led to a steady stream of advances in our understanding of many of its groundstate properties, 4,5 its thermodynamic properties, [6][7][8] and the structure of its excitation spectrum. [9][10][11][12][13][14][15] The T = 0 phase diagram of the XXZ model, which was rigorously established by these advances, consists of a ferromagnetic phase at ∆ ≥ 1, a critical phase (spin-fluid, Luttinger liquid) at −1 ≤ ∆ < 1, and an antiferromagnetic phase at ∆ < −1.
The mapping between the XXZ model and the exactly solvable 6-vertex and 8-vertex models yielded additional ground-state properties of the former on a rigorous basis, notably the spontaneous staggered magnetization in the antiferromagnetic phase and some critical exponents in the spin-fluid phase. 16,4,17 Until recently, exact results for the T = 0 spin dynamics of the XXZ model were limited to a single nontrivial case, the XX model (∆ = 0). For this case, the spin system is equivalent to a system of free lattice fermions, 18 and the dynamic spin correlation functions can be expressed as fermion density correlations (zz) 19 or as infinite determinants or Pfaffians (xx). 20,21 In the surrounding spin-fluid phase (−1 ≤ ∆ < 1), exact results for the infrared singularities of dynamic structure factors were obtained by field-theoretic approaches. 4,22 A new avenue for the study of the T = 0 dynamics of the XXZ model on a rigorous basis was opened up by important advances in the study of this model and the related vertex models in the framework of the algebraic analysis based on quantum group symmetries. A detailed description of this method with all the results which our calculations build on can be found in a recent book by Jimbo and Miwa. 23 Unlike the Bethe ansatz, this approach considers an infinite chain from the outset and exploits the higher symmetry of the infinite system (compared to the finite system) described by the quan-tum group U q (sl 2 ). 24,25 The algebraic analysis of the XXZ model for the purpose of calculating correlation functions and transition rates (form factors) of local spin operators requires the execution of the following program: (i) span the infinitedimensional physically relevant Hilbert space in the form of a separable Fock space of multiple spinon excitations and to generate the XXZ eigenvectors in this Fock space by products of spinon creation operators (so-called vertex operators) from the XXZ ground state (physical vacuum); (ii) determine the spectral properties (energy, momentum) of the spinon excitations; (iii) express the local spin operators in terms of vertex operators; and (iv) evaluate matrix elements of products of vertex operators in this spinon eigenbasis.
There exist two similar yet distinct programs which operate under different circumstances for essentially the same purpose. One is the fermion representation of the 1D s = 1/2 XY model or the equivalent 2D Ising model, 20,26 and the other is conformal field theory for critical (massless) continuum models. 27 Quantum inverse scattering theory provides yet different ways of calculating some correlation functions and matrix elements for massive relativistic continuum models 28 and for the XXX model. 29 The algebraic analysis 23 operates in the massive phase stabilized by Néel long-range order at ∆ < −1, but the isotropic limit ∆ → −1 − can be performed meaningfully at various stages of the calculation and thus yields equivalent results for the (massless) Heisenberg antiferromagnet. 30,31 In this paper we infer from the diverse ingredients now accessible via the Bethe ansatz and the algebraic analysis, an explicit expression for the exact 2-spinon part of the dynamic spin structure factor at T = 0 and ∆ < −1. The line shapes thus obtained are of direct relevance for the interpretation of existing spectroscopic data obtained via inelastic neutron scattering [32][33][34] and Raman scattering 35 on the quasi-1D magnetic compounds CsCoCl 3 and CsCoBr 3 . In Sec. II we discuss the m-spinon eigenbasis and infer a suitable parametric representation of the energymomentum relation for spinon excitations from it. In Sec. III a closed-form expression for the 2-spinon density of states D(Q, ω) is derived from this spectral information. In Sec. IV we analyze the matrix elements (form factors) between the twofold degenerate ground state and the 2-spinon excitations and derive from them (in Sec. VI), after having solved the 2-spinon energymomentum relations in the appropriate parametrization (Sec. V), a function M (Q, ω) which, when multiplied with D(Q, ω), yields the 2-spinon dynamic structure factor S (2) xx (Q, ω).
A related study was previously undertaken by Weston and Bougourzi. 36 In that study the goal was to calculate the 2-spinon part of the dynamic structure factor S(Q, ω) defined as the Fourier transform of σ n (t) · σ 0 . The result was expressed as an expansion about the Ising limit (∆ → −∞) carried out explicitly to 12 th order. It is much more difficult to calculate this quantity than to calculate S (2) xx (Q, ω), where no expansion is necessary to obtain explicit results.
Finally, it is interesting to note that the exact result for the frequency-dependent spin autocorrelation function Φ xx (ω) ≡ π −π (dQ/2π)S xx (Q, ω) of the case ∆ = 0, which was calculated in the fermion representation, 37 represents all m-spinon contributions for m = 2, 4, . . . simultaneously. There the m-spinon structure of the excitation spectrum is reflected in Φ xx (ω) by an infinite sequence of singularities at the band-edge frequencies ω/J = 0, 1, 2, . . ..

II. SPECTRUM
The 2 N -dimensional Hilbert space of the XXZ model for a chain of N sites becomes non-separable in the limit N → ∞. However, for the infinite chain, a separable subspace F can be constructed, and all physical properties of the XXZ model can, in principle, be derived exactly from it. The classification of the XXZ spectrum in terms of m-spinon excitations, which is instrumental in the quantum group analysis, had already been established by Faddeev and Takhtajan, 10,15 for ∆ = −1 in the framework of the algebraic Bethe ansatz.
The (infinite-dimensional) space F is spanned by vectors |ξ m , ǫ m ; . . . ; ξ 1 , ǫ 1 j with m = 0, 1, . . . and j = 0, 1, which represent multiple spinon excitations. In the regime of interest here, the twofold degenerate vacuum state is represented by the two vectors |0 0 , |0 1 . These states break the translational symmetry of H. The translation operator T (shift by one lattice site) transforms the two vectors into each other: (2.1) In the Ising limit (∆ → −∞), they become the pure Néel states, | . . . ↑↓↑↓ . . . , | . . . ↓↑↓↑ . . . . Each spinon excitation is characterized by a (complex) spectral parameter ξ l and a spin orientation ǫ l = ±1. The subspaces of F with even and odd numbers of spinon excitations are disconnected in all matters of concern here. They describe the physics of chains with even and odd N asymptotically for N → ∞. 10,15 The completeness relation for the spinon basis in F reads 23 with the respective eigenvalues determined by in terms of the spectral parameter ξ and the anisotropy parameter Here q is the deformation parameter of the quantum group U q (sl 2 ), 24 and For most of the analysis to be carried out later, it is convenient to express ξ in terms of the alternative spectral parameter β: (2.7) The energy and momentum of a spinon are then expressed in terms of Jacobian elliptic functions, The anisotropy parameter (2.5) is related to the nome and thus determines the moduli k, k ′ ≡ √ 1 − k 2 of the elliptic integrals K ≡ K(k), K ′ ≡ K(k ′ ). The spinon energy-momentum relation resulting from (2.8), is equivalent to the corresponding relation obtained via Bethe ansatz. 7,15,38 For the calculation of S (2) xx (Q, ω) from the 2-spinon density of states and the 2-spinon matrix elements we introduce here translationally invariant vacuum states, which have wave numbers (total momenta mod 2π) 0 and π, respectively, in the extended Brillouin zone (−π, +π). 39 In the isotropic limit, the state |0 is a singlet (S T = 0), and the state |π is the vector with S z T = 0 of a triplet (S T = 1). The corresponding linear combinations of 2-spinon states, are then also translationally invariant, Since the 2-spinon momenta and energies are independent of the spin orientations ǫ 1 , ǫ 2 = ±1, all 2-spinon states at fixed P will be at least fourfold degenerate. In the isotropic limit, this degeneracy involves a singlet state (S z T = S T = 0) and the three vectors with S z T = 0, ±1 of a triplet state (S T = 1). The four sets of 2-spinon excitations are readily identified in the framework of the Bethe ansatz. In a finite system (N < ∞), the singlet-triplet degeneracy is removed, and for anisotropic coupling (∆ < −1), the triplet levels are split up as well. The fourfold degeneracy emerges only asymptotically for N → ∞, and thus reflects the higher U q (sl 2 ) symmetry of the infinite system, which is used in the algebraic analysis.
The 2-spinon density of states 41 was evaluated before in closed form: 12 for (Q, ω) ∈ C ± , respectively, and where With the auxiliary quantity the result (3.6b) can be written more compactly: .
Note the reflection symmetry: The 2-spinon density of state has square-root divergences all along the lower and upper boundaries of each sheet. At the zone center, expression (3.9) turns into The function D + (Q, ω) is plotted in Fig. 2 for two values of anisotropy.

IV. MATRIX ELEMENTS
The m-spinon eigenbasis provides a useful framework for the separate analysis of the m-spinon contributions (m = 0, 2, 4, . . .) to any zero-temperature dynamical quantity of interest if a means of calculating the relevant matrix elements can be found. Here we focus on the 2-spinon matrix elements of the dynamic spin structure factor S xx (Q, ω) at T = 0.

V. ENERGY-MOMENTUM RELATIONS
Performing the integrals over the spectral parameters in expression (4.6) brings the 2-spinon dynamic structure factor into the form: where the numerator, is governed by the 2-spinon transition rates and the denominator, by the 2-spinon density of states. In these expressions, the spectral parameters now have fixed values (β c 1 , β c 2 ). These values are the solutions of the 2-spinon energymomentum relations arising from the two products of δ-functions in (4.6): for the future analysis. For fixed σ and at a generic point (Q, ω) within the range of the 2-spinon continuum, there exists exactly one distinct solution per sheet C ± . Every such solution has multiplicity 8, accounted for by the permutation symmetry β 1 ↔ β 2 of (5.5) [factor 2] and the periodicity of the elliptic functions [factor 4]. Now we use addition theorems 42 to convert (5.5) into we obtain These solutions yield and reduce (5.6a), effectively, into a quadratic equation for dnβ − with σ-independent solutions where ω 0 (Q) and T (Q, ω) are given in (3.2a) and (3.7), respectively. Finally, a Landen transformation (k → κ) converts (5.9) and (5.10) into more explicit solutions in terms of incomplete elliptic integrals, where W c (Q, ω) is given in (3.8), and the new label c = ± indicates that (Q, ω) ∈ C ± .
The function β + − (Q, ω) = β − − (π − Q, ω), which alone enters the final result, is plotted in Fig. 3 for two values of anisotropy. It is finite along the lower boundary of C + , decreases monotonically with increasing ω at fixed Q, and vanishes in a square-root cusp at the upper boundary. Note the different behavior along the portions ω + (Q) and ω 0 (Q) of the lower boundary of C + , β + − (Q, ω + ) = K, β + − (Q, ω 0 ) < K, which will give rise to different singularities in S (2) xx (Q, ω) in the two parts of the spectral threshold.

VI. DYNAMIC STRUCTURE FACTOR
A. Exact result for S (2) xx (Q, ω) The 2-spinon dynamic structure factor (5.1) for both sheets C ± of the 2-spinon spectrum will now be evaluated With the solutions (5.11) of the energy-momentum relations, the numerator (5.2) and the denominator (5.3) can be evaluated in the form: l sinh(2lγ) cosh(γl) e ∓γl (6.4) with γ = πK ′ /K, and ϑ d (x) is a Neville theta function. 42 With the exact results (3.9), (6.2), and (6.3), the physically motivated factorization (6.1) of S (2) xx (Q, ω) = c=± S c (Q, ω) can now be established: . (6.6) The exact 2-spinon transition rate function thus obtained from (6.2) substituted into (6.6) is β c − is given in (5.11b), and t ≡ 2β/K ′ . The final result for the exact 2-spinon dynamic structure factor reads: B. Line shapes and singularity structure The function M + (Q, ω), which represents the 2-spinon transition rates for (Q, ω) ∈ C + , is plotted in Fig. 4 for two values of anisotropy. The product of M + (Q, ω) with the 2-spinon density of states D + (Q, ω) (already shown in Fig. 2 for the same two cases) yields the spectral-weight distribution S + (Q, ω) of the 2-spinon dynamic structure factor for (Q, ω) ∈ C + . This function is plotted in Fig. 5 for four values of anisotropy, including the values chosen in Figs. 2-4.
The transition rate function M + (Q, ω) exhibits qualitatively different properties in the two regimes Q κ ≤ Q < π − Q κ and π − Q κ < Q ≤ π, where the spectral threshold is given by ω 0 (Q) and ω + (Q), respectively (see Fig. 1). In the first regime, M + (Q, ω) is nonzero at the lower boundary, decreases monotonically with increasing ω, and approaches zero linearly at the upper boundary. In the second regime, by contrast, M + (Q, ω) approaches zero linearly at both boundaries and has a smooth maximum in between. The transition rates for (Q, ω) ∈ C − are the exact mirror image: M − (π − Q, ω) = M + (Q, ω).
These properties of M ± (Q, ω) imply that the 2-spinon dynamic structure factor S (2) xx (Q, ω) diverges all along the portion ω 0 (Q) of the spectral threshold and that the leading singularity there is the square-root divergence of the 2-spinon density of states. Here the function S + (Q, ω) decreases monotonically from infinity at the lower boundary to zero at the upper boundary. Along the portion ω + (Q) of the lower boundary and along the entire upper boundary of C + , the linear behavior of the transition rates removes the square-root divergence of the density of states in the product and replaces it by a square-root cusp in the dynamic structure factor. Here the spectralweight distribution at fixed Q has a smooth maximum between the band edges. For strong anisotropy, the line shapes are broad and featureless. At moderate to weak anisotropy, the line shapes are distinctly asymmetric with the maximum positioned close to the spectral threshold.
The function S xx (Q, ω), which is symmetric about Q = π/2, is equal to one or the other of the two func- tions S ± (Q, ω) except for Q κ ≤ Q ≤ π − Q κ . Here the two sheets overlap, and the 2-spinon dynamic structure factor has two square-root cusp singularities, one at the upper boundary of each sheet. The two upper boundaries coincide only at Q = π/2.

C. Isotropic limit
When we take the isotropic limit ∆ → −1 − in the results for the 2-spinon density of states, transition rates, and dynamic structure factor for the purpose of linking up with results of calculations that were performed for the Heisenberg model in the (non-degenerate) critical ground-state, we must heed the fact that the size of the Brillouin zone changes from (−π/2, +π/2) to (−π, +π) as the Néel long-range order in the ground state of the infinite system vanishes.
In practice, this means that we switch our perspec-tive from considering both sheets C ± of 2-spinon excitations over the range (−π/2, +π/2) to considering only the sheet C + over the extended range (−π, +π). The term with c = − in all expressions that contain a sum c=± is then be omitted. These contributions are now accounted for in the term with c = + over the extended range of physically distinct wave numbers.
In the isotropic limit, the boundaries of C + turn into the familiar sine curves, 10) and the 2-spinon density of states becomes (6.11) A major simplification occurs in the 2-spinon matrix elements, X 1 (ξ 2 , ξ 1 ) → − X 0 (ξ 2 , ξ 1 ), (6.12) which implies that all terms with σ = − in (5.1) and (6.2) disappear in the isotropic limit. The results presented previously for the isotropic case 30,31 can be recovered from the more general result presented here, if we replace the spectral parameter β in (2.7) by the scaled spectral parameter α = −βπ/2K ′ , and then take the limit q → −1 − . In this limit, the auxiliary expressions (4.5) can be simplified considerably: sinh 2x cosh x e ∓x   (6.14) represents (6.4) with x = γl in the limit γ → 0. The magnitude of the resulting 2-spinon matrix element becomes 44 . were discussed in considerable detail.

D. Ising limit
When we analyze the exact 2-spinon dynamic structure factor for very strong anisotropy (∆ → −∞), a convenient expansion parameter about the limiting Ising model is κ as defined in (3.3). Expressed in terms of this parameter, the exchange anisotropy (2.5) and the amplitude of the 2-spinon dispersion (2.9) become Here we set J = 1. To lowest order in κ, the overlap region of the two sheets C ± that make up the 2-spinon continuum (see Fig. 1) collapses to a single point at Q = π/2. The continuum boundaries, as expressed in terms of the reduced frequency Ω ≡ ω/2I − 1, (6.21) are now described by the curves The regime Q κ ≤ Q ≤ π − Q κ of the 2-spinon continuum is thus squeezed to zero width at Q = π/2, where the 2-spinon bandwidth is zero as well in lowest order. The expansion of the 2-spinon dynamic structure factor can be carried out in the final result (6.9) by using yielding the closed-form expression which is identical (in lowest order) to the result obtained by Ishimura and Shiba 13 from a first-order perturbation calculation about the Ising limit. It is instructive to perform the expansion at an earlier stage of the calculation, namely in (4.6), which then becomes × δ(Q +P ) + δ(Q − π +P ) (6.28) with Ē (β 1 , β 2 ) → 2 κ 1 + κ − κ(sin 2 β 1 + sin 2 β 2 ) , (6.29a) and where we have used The asymptotic 2-spinon energy-momentum relations (5.5) are − ω + 2 + 2/κ = 2(sin 2 β 1 + sin 2 β 2 ), (6.31a) −σ sin Q = sin(β 1 + β 2 ), (6.31b) and the solutions (6.32b) have multiplicity 8 for (Q, ω) within the boundaries (6.22) of the asymptotic 2-spinon continuum. Performing the integrals in (6.28) then yields the asymptotic version of (5.1): from which (6.26) is recovered upon evaluation.

E. Experiments
Two of the most intensively studied physical realizations of the 1D s = 1/2 Heisenberg-Ising antiferromagnet are the quasi-1D magnetic compounds CsCoCl 3 and CsCoBr 3 . A very comprehensive set of spectroscopic data, which probe diverse aspects of the lowtemperature spin dynamics of these materials, is now available from several inelastic neutron scattering experiments performed over the course of 15 years at the Brookhaven 32 , Chalk River 33,45 , Laue-Langevin 46 , and Rutherford 34 Laboratories. This impressive collection of data is complemented by a set of high-resolution spectroscopic data for the long-wavelength spin dynamics obtained via Raman scattering. 35 The basis for the interpretation of all the experimental data that involve the frequency range now known under the name 2-spinon continuum has been the perturbation calculation about the Ising limit of (1.1) carried out to first order by Ishimura and Shiba 13 (see also Ref. 47), which yielded the explicit expression (6.27) for the T = 0 dynamic structure factor S xx (Q, ω). Whereas this calculation does reproduce the 2-spinon continuum boundaries correctly to first order in the expansion parameter, the reliability of its line-shape prediction -a broad peak with the maximum near the center of the band and steep drops near both boundaries -has remained very much in question.
The fact is that the line shapes observed in all experiments turned out to be highly asymmetric with a high concentration of intensity near the spectral threshold and a tail extending to the upper continuum boundary.
Various attempts have been made at reconciling the discrepancy between theory and experiment by considering a second order perturbation calculation 48 for the pure XXZ model, and by considering the impact of biaxial anisotropy, 49 next-nearest-neighbor coupling, 50 interchain coupling, 51 and exchange mixing, 34 all within the framework of a first-order perturbation treatment.
In the range −10 < ∼ ∆ < ∼ −7 of anisotropies, which best describes CsCoCl 3 and CsCoBr 3 according to some indicators, the line-shape predictions obtained via perturbation calculation are in fair agreement with the exact 2-spinon result for wave numbers near the zone center. However, for wave numbers near the zone boundary, the exact 2-spinon spectral-weight distribution is more asymmetric, consistent with the experimental data. For Q κ < Q < π − Q κ , S (2) xx (Q, ω) even diverges at the spectral threshold and decreases monotonically toward the upper continuum boundary. The conspicuous line-shape asymmetry found in the experimental data near the zone center, which is not at all reproduced by the perturbation calculation, does also exist in S (2) xx (Q, ω), but only for weaker exchange anisotropy (see Fig. 5).