Computational probes of collective excitations in low-dimensional magnetism

The investigation of the dynamics of quantum many-body systems is a concerted effort involving computational studies of mathematical models and experimental studies of material samples. Some commonalities of the two tracks of investigation are discussed in the context of the quantum spin dynamics of low-dimensional magnetic systems, in particular spin chains. The study of quantum fluctuations in such systems at equilibrium amounts to exploring the spectrum of collective excitations and the rate at which they are excited from the ground state by dynamical variables of interest. The exact results obtained via Bethe ansatz or algebraic analysis (quantum groups) for a select class of completely integrable models can be used as benchmarks for numerical studies of nonintegrable models, for which computational access to the spectrum of collective excitations is limited.


Introduction
Notwithstanding the fact that experimental and theoretical studies of condensed matter systems are fundamentally complementary to each other, they share important features, which we wish to illuminate here in the context of quantum many-body dynamics. The object of study is a sample for experimental probes and a model for computational probes. Whether the model is regarded as a mathematical idealization of a real chunk of matter or whether the sample is viewed as a physical realization of a system defined with mathematical rigor may be a matter of perspective, but the influence exerted by theory and experiment is mutual. Any observation of note calls for an explanation, while any prediction of substance brings into motion attempts at verification.
When performing experimental or computational probes in quantum many-body dynamics, the goal of researchers is very much the same, namely trying to make sense of what, in general, is a jumble of fluctuations to which the fundamental degrees of freedom contribute in various combinations and configurations. Experimental probes begin with a measurement and computational probes with a calculation. The results of either probe then require an interpretation in terms and concepts that are common to both approaches.
Qualitatively, the specification of a quantum many-body system involves information on composition, interaction, and environment. Theoretically this information is encoded in the Hamiltonian and in the density operator, whereas experimentally part of this information is manifest or hidden in the sample and other parts are controlled by the instrumental setup.
Any system thus specified is subject to fluctuations, of which we distinguish three kinds. Depending on the circumstances, each kind may govern the dynamical properties of the system. Thermal fluctuations are likely to be dominant in any statistical ensemble at elevated temperatures. Parametric fluctuations are manifestations of quenched random inhomogeneities in composition (e.g. dilution) or interaction (e.g. random bonds or fields), leaving distinct marks in dynamical properties, especially at low temperatures.
In the absence of thermal and parametric fluctuations, quantum fluctuations remain present. They are a direct consequence of the (autonomous) quantum time evolution and the time-delayed projections that are part of any dynamical probe, be it experimental or computational. No zero point motion exists in classical Hamiltonian systems. Here the ground state corresponds to a stable fixed point in phase space, and all dynamical variables are constant.
Consider a (non-random) quantum many-body system in a pure quantum state. The quantum fluctuations then depend on three quantities: (i) the Hamiltonian H, (ii) the state, which we take here to be the ground state |G , and (iii) some dynamical variable, here expressed by operator A. The Hamiltonian governs the dynamics deterministically for as long as the system can be regarded in isolation. It makes no difference whether this time evolution is viewed in the Schrödinger or Heisenberg representation: The statistical element is introduced when we subject the time evolved state |ψ(t) or, equivalently, the time evolved dynamical variable A(t), to an observation. In practice, this step involves the observation or evaluation of a dynamic correlation function, which can then be viewed either as the projection of the time evolved state |ψ(t) or as the product of the dynamical variable A(t) at two different times evaluated in the stationary state |G : zuoz2.1: submitted to World Scientific on March 21, 2022 The sum total of quantum fluctuations in a typical many-body system contains many different dynamical modes no matter in which state the system happens to be. Only once we pick a set of dynamical variables (one variable at a time) do we begin to gain insight into the role of specific modes taken from a multitude of quantum fluctuations. Criteria for choosing dynamical variables include experimental accessibility, computational amenability, or any specific purpose like the study of ordering tendencies or phase transitions.
On the microscopic level, dynamical probes are bound to be somewhat heavy-handed. No observation without perturbation! However, in all situations considered here, we shall assume that the interaction between the probe and the sample takes place in the regime of linear response. 1 This assumption is quite realistic for neutron scattering under most circumstances. Consider the weakly perturbed Hamiltonian H(t) = H 0 − b(t)B, where the external field b(t) couples to the dynamical variable B of the system H 0 . The linear response of any other dynamical variable A to that perturbation, is then determined by Kubo's formula for the response function, Its Fourier transform, the generalized susceptibility χ AB (ω + iǫ), has (for ǫ → 0) a symmetric real part and an antisymmetric imaginary part: χ AB (ω) = χ ′ AB (ω) + iχ ′′ AB (ω). The latter is related to the Fourier transform of the correlation function A(t)A , named the structure function via the fluctuation-dissipation theorem. The spectral representation of this function, in particular the expression resulting in the limit T → 0, where thermal fluctuations are absent, demonstrates that observing quantum fluctuations in the linear response regime means observing collective modes and their transition rates from the ground state. In the following, we look at quantum fluctuations and the associated collective modes from this particular angle for a number of situations as investigated by a variety of methods. a

Magnons Excited from Valence-Bond-Solid State
Consider the one-dimensional (1D) s = 1 Heisenberg model with bilinear and biquadratic exchange At T = 0 it has a phase with ferromagnetic long-range order (LRO) (−π < θ < −3π/4 or π/2 < θ ≤ π), a phase with dimer LRO (−3π/4 < θ < −π/4), the Haldane phase with hidden topological LRO (−π/4 < θ < π/4), and an obscure phase (π/4 < θ < π/2) that was named trimerized. 2 These ordering tendencies help us choose suitable dynamical variables when we explore the spectrum of collective excitations. All dynamical variables used for this purpose will have the form of fluctuation operators, where the operator A n acts locally at or near lattice site n. The associated structure function (7) is called the dynamic structure factor: where the sum runs over the dynamically relevant excitations |λ with energy ω λ = E λ − E G and spectral weight W A λ = 2π| G|F A Q |λ | 2 . To calculate S AA (Q, ω) for H θ we can employ special methods at the parameter values θ = ±π/4 where it is completely integrable 3 or general methods otherwise. One suitable general method in the present context is the recursion method in combination with a finite-size continued-fraction analysis. 4 Here we introduce four different fluctuation operators F A Q . The local operators A n from which they are constructed are listed in Table 1.
• The spin fluctuations, probed by F S Q , represent Néel order parameter fluctuations for Q = π. They are expected to be strongest at the critical point θ = −π/4, where the Q = π excitations are gapless.
• The dimer fluctuations, probed by F D Q , are also expected to be strongest (for Q = π) at θ = −π/4, which marks the onset of dimer LRO.
• The trimer fluctuations, probed by F T Q , are constructed from projection operators P T n onto local trimer states | [1,2,3] . The state |[1, 2, 3] happens to be the ground state of H θ with N = 3 for arctan 1 3 ≤ θ ≤ π/2, which was interpreted as suggesting that a trimerized phase might exist for N → ∞.
• The center fluctuations, probed by F Z Q , are constructed from a modified spin operator and tune into existing period-three (+, 0, −) or (−, 0, +) patterns of local spin states. Finite-N data indicate that for π/4 ≤ θ ≤ π/2, the fluctuations probed by F Z 2π/3 are the strongest of all the ones listed. Table 1. Local operators An which define the four fluctuation operators F A Q via (9) and the four order parameters P A via (11).
The order parameters associated with the four fluctuation operators defined by Eq. (9) and the entries of Table 1 can be written in the form where Q S = Q D = π, and Q T = Q Z = 2π/3. Each order parameter P A has a set of eigenvectors |Φ A k , k = 1, 2, . . . , K which represent the associated LRO in its purest form. The degree of degeneracy is K = 2 for the Néel and dimer states, K = 3 for the trimer states, and K = 6 for center states. 5 The physical vacuum chosen here is very unlike any of the states |Φ A k . Within the Haldane phase, at the parameter value θ V BS = arctan 1 3 , the ground state of H θ is a realization of the valence-bond solid (VBS) wave function, 6 which is non-degenerate and in which the Néel, dimer, trimer, and center ordering tendencies are all imperceptibly weak. In the VBS state, the spin 1 at each lattice site is expressed as a spin-1/2 pair in a triplet state. The singlet-pair forming valence bond involves one fictitious spin 1/2 from each of two neighboring lattice sites. The VBS state, which has total spin S T = 0, can then be regarded as a chain of valence bonds linking successive spin-1/2 pairs in this manner.
The topological LRO present in the Haldane phase and known to be zuoz2.1: submitted to World Scientific on March 21, 2022 strongest in the VBS state provides an environment, where a specific kind of elementary excitations can propagate freely and where the associated stationary states form a branch with well defined dispersion. We now probe these elementary excitations and composites thereof from different angles by the four fluctuation operators F A Q defined in Table 1.
This produces the selection rules ∆S T = 0 in the dynamic dimer and trimer structure factors. The corresponding selection rules for the dynamic spin and center structure factors are ∆S T = 0, 1 and ∆S T = 0, 1, 2, respectively, with the further restriction that transitions between singlets (S T = 0 states) are prohibited. At the VBS point, S DD (Q, ω) and S T T (Q, ω) thus couple exclusively to the S T = 0 excitation spectrum, and S SS (Q, ω) exclusively to the S T = 1 excitation spectrum, whereas S ZZ (Q, ω) couples to the S T = 1 and S T = 2 spectra.
In Fig. 1 we display ω A λ versus Q of the dynamically relevant spin, center, dimer, and trimer excitation spectra as obtained from the finite-size continuedfraction analysis. 5 The filtered access to the spectrum afforded by the four fluctuation operators gives us valuable clues about the nature of the elementary excitations that thrive in the VBS environment. The low-frequency region at Q π/2 in the spin and center spectra is dominated by a branch of states with S T = 1. They carry more than 95% of the spectral weight in S SS (Q, ω) and S ZZ (Q, ω). These triplets, which have been named magnons, remain invisible in the dimer and trimer spectra. Only composites of the magnon states may be observable in S DD (Q, ω) and S T T (Q, ω).
A dispersion of the general form ω M (Q) = J(a + b cos Q) for the 1-magnon branch can be inferred from the single-mode approximation of S SS (Q, ω). 5,6 Under the assumption that in the VBS environment the magnons are weakly interacting point particles, we can expect the existence of three kinds of 2magnon scattering states formed by pairs of 1-magnon triplets: states with S T = 1, which contribute to S SS (Q, ω) and S ZZ (Q, ω), states with S T = 0, which contribute to S DD (Q, ω) and S T T (Q, ω), and states with S T = 2, which are observable in S ZZ (Q, ω) only. Free 2-magnon states form a two-parameter The predicted coalescence of the 2-magnon continuum into one spectral line at Q = π follows from the symmetry property ω M (Q)+ω M (π−Q) = const of the magnon dispersion. Interestingly, the finite-N dimer spectrum does indeed collapse into a single spectral line at the N -independent excitation energy ω D = √ 10J, which carries all the spectral weight in S DD (Q, ω). We now use the exact 2-magnon excitation energy, ω ± (π) = ω D , and the The dashed lines represent 1-magnon dispersions and the solid lines 2-magnon continuum boundaries as explained in the text. extrapolated value, ω M (π) = 0.66433(2)J, of the 1-magnon excitation gap to fit the parameters a, b. The resulting values, a ≃ 1.581, b ≃ 0.917, used in ω M (Q) and ω ± (Q) yield the dashed and solid lines in Fig. 1. 5 Both the 1-magnon dispersion and the 2-magnon spectral threshold are in very good agreement with all finite-N data shown. Hence the magnon interaction is very weak at the bottom of the 2-magnon region in all three S T subspaces. The finite-N data spilling out on top of the shaded areas in Fig. 1 suggest that at higher energies, the magnon interaction is repulsive, more strongly so in the S T = 0, 2 subspaces than in the S T = 1 subspace.
This raises the possibility that bound 2-magnon states split off the top of the 2-magnon continuum of 2-magnon scattering states in the S T = 0, 2 subspaces. The comparison of panels (a) and (b) at frequencies 3 ω/J 5 indeed suggests that the dynamically relevant finite-N excitations are arranged in contrasting patterns. In panel (a) we have an arrangement of points which is typical of a two-parameter continuum. As N increases, more points are added and spread roughly evenly along the frequency axis. In panel (b), by contrast, the data points are arranged in branches with an almost N -independent separation, which is characteristic for branches of bound states. The spin-1 while not physical realizations of H θ directly, nevertheless realize situations where the spectrum of collective excitations as probed by inelastic neutron scattering 7 shares major features with the magnons excited from the VBS state.

Magnons Excited from Ferromagnetic State
We now turn to a more quantitative discussion of the difference between scattering states and bound states in the context of completely integrable situations, namely for the 1D s = 1 2 Heisenberg ferromagnet: The ground state is (N + 1)-fold degenerate. We select one of the groundstate eigenvectors, |F ≡ | ↑↑ · · · ↑ , as the physical vacuum for an exploration of collective excitations. As in the VBS case discussed in Sec. 2, the spin fluctuation operator probes transitions between the vacuum state and a branch of 1-magnon states. However, here the magnon dispersion is gapless: Unlike in the VBS case, the 1-magnon states of H F are located in a separate invariant subspace. Only the 1-magnon states can contribute to the spin fluctuations. In the expression for the 1-magnon eigenvectors, the spin fluctuation operator plays the role of a magnon creation operator. That is not to say S − k is a true magnon creation operator. Multi-magnon superpositions, i.e. the states S − k1 · · · S − kr |F , are a redundant set of non-orthogonal and non-stationary states, which are subject to two kinds of interactions: 8 • The kinematical interaction is caused by the restriction on the number of reversed spins at one lattice site.
• The dynamical interaction is caused by the off-diagonal part of H F in the basis of multi-magnon superpositions.
This distinction is quite natural in the framework of the Bethe ansatz 9 as we shall see. The Bethe ansatz is an exact method for the calculation of eigenvectors and eigenvalues of completely integrable quantum many-body systems. The Bethe wave function of any eigenstate in the subspace with r ≡ N/2 − S z T reversed spins relative to the magnon vacuum, has coefficients of the form determined by r magnon momenta k i and one phase angle θ ij = −θ ji for each magnon pair. The sum P ∈ S r is over the permutations of the labels {1, 2, . . . , r}. The consistency requirements for the coefficients a(n 1 , . . . , n r ) inferred from the eigenvalue equation H|ψ = E|ψ and the requirements imposed on the same coefficients by translational invariance can be cast in a set of equations for the momenta k i and phase angles θ ij : Every solution of these equations is specified by a set of r Bethe quantum numbers λ i ∈ {1, 2, . . . , N − 1}. Given a solution, the energy and wave number of the state it describes are In the subspace with r = 1, we thus recover all N 1-magnon states, one for each of the allowed values of λ 1 . There exist N (N +1)/2 distinct 2-magnon superpositions of the 1-magnon states thus identified. However, this set of states must be accommodated in the r = 2 subspace, whose dimensionality is only N (N −1)/2. Inspection shows that of the N (N +1)/2 pairs of Bethe quantum numbers in the allowed range 0 ≤ λ 1 ≤ λ 2 ≤ N − 1, N pairs do indeed not produce a solution of (16). The missing solutions are a consequence of the kinematical interaction between magnons. The consequences of the dynamical 2-magnon interaction are illustrated in Fig. 2, which shows the complete r = 2 spectrum (E − E F )/J versus k for k ≥ 0 and N = 32. 10   The class C 1 contains N states for which one of the two Bethe quantum numbers is zero. This means that one of the magnons has zero wave number. Its effect is a slight rotation of the magnon vacuum, in which the other magnon is as free to propagate as in the original vacuum. There is no dynamical interaction. All states in this class are, effectively, 1-magnon states.
The class C 2 of states has nonzero Bethe quantum numbers λ 2 −λ 1 ≥ 2. There are N (N−5)/2+3 such pairs. All of them yield a solution with real k 1 , k 2 , which makes them 2-magnon scattering states. A measure of the dynamical magnon interaction in Fig. 2 is the vertical displacement of any true scattering state (•) from the nearest free-magnon pair (+). As N increases, the energy correction diminishes for all class C 2 states and vanishes in the limit N → ∞. The 2-magnon scattering states and the free 2-magnon states then form twoparameter continua with identical boundaries E −E F = 2J[1 ± cos(k/2)].
The class C 3 of states has nonzero Bethe quantum numbers which either are equal λ 2 = λ 1 , or differ by unity λ 2 = λ 1 + 1. There exist 2N − 3 such pairs, but only N −3 pairs yield solutions of (16). For the class C 3 states, the effects of the dynamical magnon interaction are much more prominent, and the interaction energy does not disappear when N → ∞. In Fig. 2, these states form a branch of 2-magnon bound states with dispersion E−E F = 1 2 J(1−cos k) below the continuum of 2-magnon scattering states.
The bound state character of the class C 3 states manifests itself in the enhanced probability that the two flipped spins are on neighboring sites of the lattice. This property of the wave function is best captured in the weight distribution |a(n 1 , n 2 )| of basis vectors with flipped spins at sites n 1 and n 2 . In Fig. 3 we have plotted |a(n 1 , n 2 )| versus n 2 − n 1 for a sequence of class C 3 states between k = 0 and k = π. 10 The distribution is peaked at n 2 − n 1 = 1. Its width is controlled by the imaginary parts of k 1 , k 2 in (15). The smallest width is observed in the bound state at k = π. In this case, all coefficients with n 2 = n 1 + 1 are zero, which implies that the two down spins are tightly bound together and have the largest binding energy.
The width of the distribution |a(n 1 , n 2 )| increases as k decreases, and the binding of the two down spins loosens. For finite N , the Bethe ansatz solutions switch from complex to real when the distribution has acquired a certain width. In scattering states the distribution |a(n 1 , n 2 )| is always broad and tends to oscillate wildly. Some scattering states have a smooth distribution with a maximum for n 2 = n 1 + N/2, when the two down spins are farthest apart. The formation of bound states and scattering states of elementary excitations exist in many different contexts. But only in rare cases such as this one can they be investigated on the level of detail presented here.

Spinons Excited from Spin-Fluid State
Turning our attention to the Heisenberg antiferromagnet, we consider the Hamiltonian H A ≡ −H F with H F defined in Eq. (12). All the eigenvectors remain the same, but the energy eigenvalues have the opposite sign. The magnon vacuum |F now is at the top of the excitation spectrum. The ground state |A of H A is located in the invariant subspace with S z T = 0. This subspace also contains the two Néel states |N 1 ≡ | ↑↓↑ · · · ↓ , |N 2 ≡ | ↓↑↓ · · · ↑ , which are not eigenstates of H A . In the framework of the Bethe ansatz, |A can be obtained from |F by exciting r = N/2 magnons with momenta k i and (negative) energies −J(1 − cos k i ). The Bethe quantum numbers for this state are For reasons of computational convenience we rewrite the Bethe ansatz equations (16) in terms of the variables z i ≡ cot(k i /2): The associated Bethe quantum numbers −N/2 < I i ≤ N/2 are integers for odd r and half integers for even r.
with E F = JN/4. For states with real {z i }, Eqs. (18) can be converted into a convergent iterative process: with starting values z (1) i = πI i /N . For the ground state |A we have High-precision solutions {z i } can be obtained with little computational effort. The ground-state energy per site for N = 4096, for example, reproduces the exact result (E A −E F )/JN = − ln 2 of the infinite chain to within 1 part in a million. 11 The distribution of magnon momenta in the ground state |A is broad and peaked at k i = π:

zuoz2.1: submitted to World Scientific on March 21, 2022
This state, which has obviously a very complicated structure when described in terms of magnons, will now be configured as the physical vacuum of H A for a different kind of elementary particle called the spinon. 12 A useful way to characterize the new physical vacuum is through the perfectly regular array (21) of Bethe quantum numbers as illustrated in the first row of Fig. 4. The spectrum of H A can then be generated systemati- 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. 4 and rearrange the remaining I i in all configurations over the expanded range |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 N (N+2)/8. A generic configuration consists of three clusters with two gaps between them as shown in the second row of Fig. 4. The position of the gaps between the I i -clusters determine the momentak 1 ,k 2 of the two spinons, which, in turn, add up to the wave number of the two-spinon state relative to the wave number of the vacuum: Q ≡ k−k A =k 1 +k 2 .
A plot of the 2-spinon energies E − E A versus wave number k − k A for N = 16 as inferred via (19) from the solutions {z i } is shown in Fig. 5 (•). The dots, which represent the corresponding data for N = 256, produce a sort of density plot for the 2-spinon continuum which emerges in the limit N → ∞. The exact lower and upper boundaries of the 2-spinon continuum are 13 with s 1 = s 2 = +1/2 (S T = 1, S z T = +1). The 2-spinon triplets with s 1 = s 2 = −1/2 (S T = 1, S z T = −1) are obtained from these states by a spinflip transformation. One set of 2-spinon states with s 1 = −s 2 are the triplets with S z T = 0. They are obtained via Bethe ansatz by a simple modification of the {I i } configurations from the kind shown in the second row of Fig. 4 and yields an equal number of (real) solutions. Symmetry requires that all three triplet components (S z T = 0, ±1) are degenerate. The 2-spinon singlet states (S z T = S T = 0) are characterized by one pair of complex conjugate solutions z 1 = z ⋆ 2 in addition to the real solutions z 3 , . . . , z N/2 . Finding the I i -configurations for a particular set of eigenstates with complex solutions and then solving the associated Bethe ansatz equations is, in general, delicate task. 11,14 The 2-spinon singlets for N = 16 and 0 < Q ≤ π are shown as open circles in Fig. 5. As N grows large, the effect of the complex solutions z 1 = z * 2 on the energy relative to that of the real solutions z 3 , . . . , z N/2 diminishes. In the limit N → ∞, the 2-spinon singlets also form a continuum with boundaries (23). Yet the effect of the complex solutions will remain strong for other quantities including selection rules and transition rates.
The 2-spinon triplets play an important role in the low-temperature spin dynamics of quasi-1D antiferromagnetic compounds such as KCuF 3 , Cu(C 6 D 5 COO) 2 · 3D 2 O, Cs 2 CuCl 4 , and Cu(C 4 H 4 N 2 (N O 3 ) 2 ). They are the elementary excitations of H A which can be directly probed via inelastic neutron scattering. 15 The 2-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. Some quasi-1D antiferromagnetic compounds like CuGeO 3 are susceptible to a spin-Peierls transition, which involves a lattice distortion accompanied by an exchange dimerization. 16 The operator which probes the dimer fluctuations in the ground state of H A couples primarily to the 2-spinon singlets and not at all to the 2-spinon triplets.

Spinons Excited from Antiferromagnetic State
Although the computational application of the Bethe ansatz yields the exact wave functions of the 2-spinon states for very large systems with little effort, the very structure of the Bethe wave function makes it hard to use this knowledge for the calculation of the transition rates pertaining to any fluctuation operator of interest. Nevertheless, there exist ways to extract useful lineshape information from the Bethe wave function for specific situations. 17 Here we discuss an alternative approach, which exploits the higher symmetry of H A in the limit N → ∞, described by the quantum group U q (sl 2 ). For example, the asymptotic degeneracy of the 2-spinon triplets and singlets noted previously is attributable to this higher symmetry.
The algebraic analysis of a completely integrable spin chain employed for the purpose of calculating dynamic structure factors for specific fluctuation operators requires the execution of the following program: 18 • Span the infinite-dimensional physically relevant Hilbert space in the form of a separable Fock space of multiple spinon excitations.
• Generate the eigenvectors of the Hamiltonian in this Fock space by products of spinon creation operators (so-called vertex operators) from the ground state (physical vacuum).
• Determine the spectral properties (energy, momentum) of the multiplespinon states accessible from the ground state by the selected fluctuation operator.
• Express the fluctuation operator of interest in terms of vertex operators.
• Evaluate matrix elements of products of vertex operators as are needed for the selected fluctuation operator in the spinon eigenbasis.
In the following, we outline how this program was carried out for the calculation of the exact 2-spinon part of one dynamic structure factor for the 1D s = 1 2 XXZ model, At T = 0 H XXZ has a ferromagnetic phase for ∆ ≥ 1, a critical phase (spinfluid) for −1 ≤ ∆ < 1, and an antiferromagnetic phase for ∆ < −1. The algebraic analysis operates in the massive phase stabilized by Néel LRO at ∆ < −1, but the isotropic limit ∆ → −1 − , which is equivalent to H A , can be performed meaningfully.
These ingredients yield the following expression to be evaluated: 21 A compact rendition of the exact result reads 21 where γ = πK ′ /K, t ≡ 2β/K ′ , and ϑ d (x) is a Neville theta function. This function is plotted in Fig. 7 for (Q, ω) ∈ C + . S (2) xx (Q, ω) has a square-root divergence at the portion ω 0 (Q) of the spectral threshold. Along the portion ω + (Q) of the lower boundary and along the entire upper boundary of C + , S (2) xx (Q, ω) has a square-root cusp. The line shapes for (Q, ω) ∈ C − are inferred from the symmetry property S   The isotropic limit ∆ → −1 − is delicate because of its singular nature. As the LRO in the spinon vacuum vanishes, the size of the Brillouin zone changes from (−π/2, +π/2) to (−π, +π). A practical consequence of this phase transition is that we switch our perspective from considering both sheets C ± of 2-spinon excitations over the range (−π/2, +π/2) to considering only the sheet C + , now with boundaries (23), over the extended range (−π, +π). With these subtleties taken into account, Eq. (36) reduces to 22 This function, which is plotted in Fig. 8, has a square-root cusp at ω U (Q) and a square-root divergence with logarithmic corrections at ω L (Q). The expression is nonzero only in the shaded region of the (Q, ω)-plane bounded by ω L (Q) and ω U (Q).

Solitons Excited from Néel State
Near the Ising limit (∆ → −∞), the exact result (36) can be expanded in powers of the anisotropy parameter κ. To leading order, we obtain with Ω ≡ ω/2I − 1, which is identical (in that order) to the result obtained by Ishimura and Shiba 23 from a first-order perturbation calculation about the Ising limit. In their calculation, the two vacuum states are approximated by the pure Néel states (see Fig. 9). The spinons are kink solitons (domain walls), 24 which produce two adjacent up or down spins in an otherwise unperturbed Néel configuration. The only states which contribute to S xx (Q, ω) in leading order of the perturbation calculation are states which contain two solitons with equal spin orientation. The states with solitons of opposite spin orientation dominate the spectrum of S zz (Q, ω) except for the central peak at ω = 0, Q = π, which involves a transition between the two vacuum states.   Fig. 6, but it does not capture any of the features that characterize the regime Q κ ≤ Q ≤ π − Q κ around the Brillouin zone boundary. The limitations of the perturbation results are much less severe near the center of the Brillouin zone. This illustrated in Fig. 10, which compares the line shapes of Eq. (36) and (42) at Q = π. Here a fairly strong departure from the Ising limit is required before a significant difference between the two results is produced.  Two intensively studied physical realizations of H XXZ at ∆ < −1 are CsCoCl 3 and CsCoBr 3 . Spectroscopic data which probe both the quantum and thermal fluctuations of those materials are available from several neutron scattering experiments. 25