Finite Temperature Dynamical Structure Factor of the Heisenberg-Ising Chain
aa r X i v : . [ c ond - m a t . s t r- e l ] F e b Finite Temperature Dynamical Structure Factor of the Heisenberg–Ising Chain
A. J. A. James , , W. D. Goetze and F. H. L. Essler Department of Physics, University of Virginia, Charlottesville, VA 22904-4717, USA Rudolf Peierls Centre for Theoretical Physics, University of Oxford, Oxford OX1 3NP, UK (Dated: November 2, 2018)We consider the spin–1/2 Heisenberg XXZ chain in the regime of large Ising–like anisotropy ∆.By a combination of duality and Jordan–Wigner transformations we derive a mapping to weaklyinteracting spinless fermions, which represent domain walls between the two degenerate groundstates. We develop a perturbative expansion in 1 / ∆ for the transverse dynamical spin structurefactor at finite temperatures and in an applied transverse magnetic field. We present a unifieddescription for both the low–energy temperature-activated response and the temperature evolutionof the T=0 two–spinon continuum. We find that the two–spinon continuum narrows in energywith increasing temperature. At the same time spectral weight is transferred from the two–spinoncontinuum to the low energy intraband scattering continuum, which is strongly peaked around theposition of the (single) spinon dispersion (“Villain mode”). PACS numbers: 75.10.Jm, 75.10.Pq, 75.40.Gb
I. INTRODUCTION
The spin–1/2 Heisenberg XXZ chain is a paradigm inlow–dimensional quantum magnetism. Its Hamiltonianis H = J X n ∆ S zn S zn +1 + S xn S xn +1 + S yn S yn +1 + X n h · S n , (1)where h is an external magnetic field and ∆ controlsthe exchange anisotropy. The one–dimensional case isparticularly significant because for a field h parallel tothe ˆ z direction, the Hamiltonian is integrable and thespectrum of the spin chain may be extracted exactly. Furthermore the Hamiltonian (1) is thought to providea realistic description of several quasi–1D experimen-tal compounds. Examples include Cs CoCl for which∆ = 0 . CsCoBr , CsCoCl and TlCoCl allof which have ∆ ∼ > , h = 0 and T = 0 the XXZ spin chain isin a N´eel phase. The fundamental excitations take theform of gapped fractionalized spin–1/2 quantum solitonsknown as spinons. Strictly in the Ising limit, ∆ = ∞ ,the spinons can be identified simply as domain wallswith gap ∆ /
2. Experiments have established the exis-tence of several manifestations of the XXZ model in theIsing regime.
In some cases these experiments havealso probed the effects of temperature and transversefield on dynamical spin–spin correlations.The measure of dynamical correlations known as thedynamical structure factor is an important quantity inthe study of quantum magnets. This is for two reasons:firstly, the dynamical structure factor is directly mea-surable by inelastic neutron scattering experiments andsecondly the nature of the dynamical response is highlyspecific to the system in question, so that it serves asa characterization tool. A particular feature that onewould like to understand in the case of the XXZ chain isthe finite temperature low energy spin response known as the ‘Villain mode’. This response due to scattering be-tween domain wall pair states has been observed in Refs.[9,13]. As this response occurs only at finite tempera-tures it necessitates a theory that accounts for thermalfluctuations. A recent analysis of the continuum limit oftwo gapped integrable quantum spin chains has shownthat at raised temperatures the effect of thermal fluctu-ations cannot be described in terms of a simple thermaldecoherence or relaxation time picture. Instead markedlyasymmetric thermal broadening of single particle modesis observed. This paradigm has been found to be in agree-ment with theoretical and experimental studies ofthe spin–1/2 Heisenberg chain with strongly alternatingexchange, a model which is gapped but not integrable.In contrast the gapped excitations in the spin–1/2 XXZchain occur only in pairs; it is then interesting to try andunderstand the thermal evolution of the resulting twoparticle response in addition to that of the Villain mode.Despite the advantages afforded by integrability, thetask of calculating correlation functions (and hence dy-namics) for the XXZ chain is still far from simple.First order perturbative treatments around the ∆ = ∞ limit and 1 /S expansions have given some insight.In recent years significant progress has been made for the∆ > and a differentexact technique which works directly with the thermo-dynamic limit. This has lead to an exact expressionfor the transverse dynamical correlations at T = 0. Results at finite temperature are generally still limitedto asymptotically exact thermodynamic quantities .However there have been promising advances that rely ongeneralizing multiple integral representations for time de-pendent correlation functions to finite temperatures.
Currently these methods have not yielded expressions forthe most experimentally relevant quantity, the dynamicalstructure factor.In this paper we present a perturbative calculationfor the transverse spin response in the ∆ ≫ − < ∆ ≤ Theusual perturbative approach to the XXZ chain uses theJordan–Wigner transformation to produce an expan-sion in powers of ∆. This is suitable for investigatingthe XY, | ∆ | ≪ duality transforma-tion to a new Hamiltonian in terms of dual operators.A Jordan–Wigner transformation from these dual oper-ators to spinless fermions then leads to a controlled ex-pansion in 1 / ∆. An equivalent mapping has been usedpreviously, coupled with mean field theory, to find the ap-proximate excitation spectrum in the Ising phase. Wetake an alternative approach resumming certain terms inthe expansion to all orders, in order to take account ofboth quantum and thermal fluctuations.The structure of the paper is as follows. First in Sec.II we discuss symmetries of the Hamiltonian and the dy-namical structure factor. Second, in Sec. III we trans-form the Hamiltonian into a form suitable for the expan-sion. In Sec. IV we describe the perturbative expansionof the transverse spin–spin correlator. In Sec. V we ex-plain how to resum certain terms in this expansion inorder to obtain finite results. In Sec. VI we discuss thebehaviour of the dynamical structure factor for a rangeof parameters. Sec. VIII contains some brief concludingremarks.
II. SYMMETRIES OF THE HAMILTONIANAND THE STRUCTURE FACTOR
We now describe the symmetries of the Hamiltonianand their consequences for spin–spin correlation func-tions. For h = 0 the Hamiltonian, Eq. (1), is invariantunder arbitrary rotations R z ( φ ) around the z-axis as wellas under rotations by π around the x axis, R x ( π ), whichentail the mapping S xj → S xj ,S y,zj → − S y,zj . The two types of symmetry operations do not commute,but we can diagonalize the Hamiltonian simultaneouslywith either S z or with the generator R x ( π ) of the Z symmetry. This in turn implies that all off–diagonalspin correlators vanish for T >
0, by the following argu-ments. When considering the thermal expectation values h S an S zm i with a = x, y we choose a basis of energy eigen-states in which the total S z is diagonal. Then carryingout a rotation by π around the z-axis sends S an → − S an , a = x, y and as a result h S an S zm i = −h S an S zm i = 0 , a = x, y. (2)On the other hand, when considering h S xn S ym i we use abasis of simultaneous eigenstates of H and R x ( π ) to carry out the thermal trace. Under the Z symmetry the ther-mal expectation value is negated, leading to h S xn S ym i = −h S xn S ym i = 0 . (3)This shows that in the absence of the transverse magneticfield all off–diagonal elements of the dynamical structurefactor vanish. In the presence of a finite transverse fieldwe only have the Z to work with. Concomitantly for h > h S xn S ym i = h S xn S zm i = 0 but h S yn S zm i is nolonger required to vanish by symmetry and as a resultacquires a finite value. III. TRANSFORMATIONS OF THEHAMILTONIAN
In order to proceed we aim to re–express the Hamilto-nian in terms of spinless fermions by means of a Jordan–Wigner transformation in such a way that we analyzethe resulting interacting fermion Hamiltonian by stan-dard perturbative methods. In order to achieve this,the ∆ J P n S zn S zn +1 part of the Hamiltonian (1) must bemapped to an expression quadratic in fermions. One ap-proach for doing this is outlined in Appendix A, anotherone is discussed in detail next. In the following we con-sider a transverse magnetic field applied along the ˆ x di-rection. We note that the effects of of a transverse fieldin the critical region of the XXZ–chain − < ∆ ≤ A. Duality Transformation
We work on the infinite chain so that we can ignoreboundary effects. The Hamiltonian, Eq. (1), in terms ofPauli spin matrices is given by H = H ∆ + H h ,H ∆ = J X n σ zn σ zn +1 + J X n (cid:0) σ xn σ xn +1 + σ yn σ yn +1 (cid:1) ,H h = h X n σ xn . (4)The Kramers–Wannier duality transformation is definedby µ xn +1 / = σ zn σ zn +1 , µ zn +1 / = Y j 1, i.e.the Ising model in a transverse field. In this limit themapping gives∆ N − X n =1 σ zn σ zn +1 + h N X n =1 σ xn → h N − X n =1 µ zn +1 / µ zn +3 / + ∆ N − X n =1 µ xn +1 / + hµ z / . (7)Now as h/ ∆ → (cid:10) σ zn (cid:11) = − (cid:10) σ zn +1 (cid:11) . In contrast, for the dual Hamiltonian theground state is given by (cid:10) µ xn +1 / (cid:11) = − n < N .The twofold degeneracy is then maintained by the N thdual spin, µ N +1 / which is free to point in either direc-tion.In the following we will be interested only in bulk cor-relations of operators that are local under the dualitytransformation. Hence the boundary terms do not playa role and will be dropped. B. Fermionic Representation In order to proceed further it is necessary to map thespins to fermions. We perform a rotation of spin axes µ xn +1 / → τ zn , µ zn +1 / → τ yn , (8) with raising and lowering operators τ + n = τ xn + iτ yn , τ − n = τ xn − iτ yn (9)and then use the Jordan–Wigner transformation: τ zn = 2 c † n c n − ,τ + n = c † n e − iπ P j 1. Spin Operators We first consider the transformation of the lattice spinoperators under the mappings. Crucially, the transversespin operator is local under the transformations σ xn = c † n − c n − c † n − c † n + h . c . . (11)On the other hand, both σ z and σ y acquire Jordan-Wigner strings. As a result, our formalism will allow us todetermine the xx –component of the dynamical structurefactor only. It follows from the symmetry considerationsabove that in absence of a transverse magnetic field thissuffices to determine all transverse correlations. 2. Hamiltonian After the Jordan–Wigner transformation the Hamilto-nian takes the form H ∆ = J X n c † n c n + J X n (cid:2) c † n − c n +1 − c † n − c † n +1 − c † n − c † n c n c n +1 + c † n − c † n c n c † n +1 + h.c. (cid:3) + const . (12)and H h = h X n (cid:16) c † n − c n − c † n − c † n + h.c. (cid:17) . (13)We now write the Hamiltonian as a sum of two pieces, H = H + H , containing terms quadratic and quartic in the fermionic operators respectively. The quartic (inter-action) terms are O (∆ ) while the quadratic pieces mixorders O (∆) and O (∆ ). The external field only appearsin the quadratic part of the Hamiltonian. After takingthe Fourier transform of the quadratic part, H , of theHamiltonian we find H = J X k (cid:0) c † k c − k (cid:1) (cid:18) A k iB k − iB k − A k (cid:19) (cid:18) c k c †− k (cid:19) (14)with A k = 2∆+4 cos(2 k )+ hJ cos( k ) and B k = 4 sin(2 k )+ hJ sin( k ). This can be diagonalized by a Bogoliubovtransformation of the form (cid:18) c k c †− k (cid:19) = (cid:18) i cos( θ k ) − sin( θ k )sin( θ k ) − i cos( θ k ) (cid:19) (cid:18) α k α †− k (cid:19) (15)so that H = 12 X k (cid:0) α † k α − k (cid:1) (cid:18) ǫ k − ǫ k (cid:19) (cid:18) α k α †− k (cid:19) (16)with tan(2 θ k ) = B k /A k and ǫ k = J p A k + B k .Now we consider the quartic part of the Hamiltonian: H = − J X n (cid:16) c † n − c † n c n c n +1 + c † n − c † n c † n +1 c n + h.c. (cid:17) . (17)Taking the Fourier transform and manipulating indicesleads to H = − J N X k ,k ,k ,k δ k + k + k + k × (cid:2) f ( k k )( k k ) c † k c † k c − k c − k + 23 (cid:16) ig ( k k k ) c † k c † k c † k c − k + h.c. (cid:17) (cid:3) , (18) where we have defined the new functions f ( k ,k )( k ,k ) = cos( k − k ) − cos( k − k ) − cos( k − k ) + cos( k − k ) , (19) g ( k ,k ,k ) = sin( k − k ) + sin( k − k )+ sin( k − k ) . (20)The function f is antisymmetric under exchange of thetwo momenta within a pair of brackets. The function g is symmetric for cyclic permutations of its momen-tum arguments and antisymmetric otherwise. For exam-ple f ( k ,k )( k ,k ) = − f ( k ,k )( k ,k ) = f ( k ,k )( k ,k ) and g ( k ,k ,k ) = g ( k ,k ,k ) = − g ( k ,k ,k ) . Performing the Bo-goliubov transformation on H is standard but lengthy.By manipulating indices under the sums we arrive at arelatively compact form for the part quartic in Bogoli-ubov operators: H = 1 N X , , , δ k + k + k + k , n V ( k , k , k , k ) h α † k α † k α † k α † k + h . c . i + h V ( k , k , k , k ) α † k α − k α − k α − k + h . c . i + V ( k , k , k , k ) α † k α † k α − k α − k o . (21)The interaction vertices are given by V ( k , k , k , k ) = J X P ∈ S sgn( P ) cos( k P (1) − k P (2) − θ k P (1) + θ k P (2) + θ k P (3) − θ k P (4) ) , (22)with permutation P acting on the set { , , , } , V ( k , k , k , k ) = i J X P ∈ S sgn( P )[sin( k − k P (2) − θ k + θ k P (2) + θ k P (3) − θ k P (4) ) − sin( k P (2) − k P (3) + θ k − θ k P (2) + θ k P (3) − θ k P (4) )] , (23)with permutation P acting on the set { , , } and finally V ( k , k , k , k ) = J (cid:16) X P ∈ S sgn( P ) cos( k − k − θ k + θ k − θ k P (1) + θ k P (2) )+ X P ′ ∈ S sgn( P ′ ) cos( k − k − θ k + θ k − θ k P ′ (3) + θ k P ′ (4) )+ X P ∈ S X P ′ ∈ S sgn( P )sgn( P ′ ) (cid:2) cos( k P (1) − k P ′ (3) − θ k P (1) + θ k P (2) + θ k P ′ (3) − θ k P ′ (4) )+ cos( k P (1) − k P ′ (3) − θ k P (1) − θ k P (2) + θ k P ′ (3) + θ k P ′ (4) ) (cid:3)(cid:17) , (24)where P and P ′ act on { , } and { , } respectively. π/2 π π/2 π k-0.4-0.200.20.4 θ ( k ) ∆=10, h=0 ∆=3, h=0 ∆=10, h=J FIG. 1: (Color Online) The self–consistent Bogoliubov pa-rameter θ k , for J = 1. The parameter scales as ∆ − . New quadratic terms are generated by normal orderingthe quartic piece. We must then include these with theoriginal terms from H and solve for θ k self–consistentlyso that the off–diagonal terms are zero. This requirementmay be recast as a self–consistency condition for every k :tan(2 θ k ) = 2 sin(2 k ) + hJ sin( k ) + N P q Θ ( k, q )∆ + 2 cos(2 k ) + hJ cos( k ) + N P q Θ ( k, q ) , (25)where we have definedΘ ( k, q ) = 2 f ( k,q )( − k, − q ) sin ( θ q ) + g ( k,q, − q ) sin(2 θ q ) , (26)Θ ( k, q ) = f ( k, − k )( q, − q ) sin(2 θ q ) − g ( k, − k,q ) sin ( θ q ) . (27)Clearly the dispersion is also affected by the new π/2 π π/2 π k ε ( k ) ∆=10, h=0 ∆=10, h=1 FIG. 2: (Color Online) The single particle dispersion relation, ǫ k (with J = 1) for ∆ = 10. The dispersion calculated withand without self–consistency (solid and dashed curves respec-tively) and the exact spinon dispersion (dotted curves)are shown. The spinon result is not available in the case ofa finite transverse field. For h = 0 the dispersion is nearlysinusoidal. Note the functions are π periodic for h = 0 and2 π periodic otherwise. quadratic parts, becoming ǫ k = J (cid:16)h ∆2 + cos(2 k ) + hJ cos( k ) + 12 N X q Θ ( k, q ) i + h sin(2 k ) + hJ sin( k ) + 14 N X q Θ ( k, q ) i (cid:17) . (28)Evaluating the self–consistency and dispersion relations(25,28) numerically we can compare the gap (i.e. low-est excitation energy) to the mean field result found byG´omez–Santos. Summing over 100 sites our results forthe physical (two particle) gap are in excellent agreement.The self–consistent Bogoliubov parameter is plotted inFig. 1 for a range of parameters. Figures 2 and 3 showthat ǫ k is an excellent approximation to the spinon dis- π/2 π π/2 π k ε ( k ) ∆=3, h=0 FIG. 3: (Color Online) The single particle dispersion relation, ǫ k (with J = 1) for ∆ = 3. The dispersion calculated withand without self–consistency (solid and dashed curves respec-tively) and the exact spinon dispersion (dotted curves)are shown. persion. It is also apparent that use of a self–consistentBogoliubov parameter is a small effect on the level of thedispersion, except in the presence of a transverse field.It is worth emphasising that the fermions that featurein the diagonalized quadratic Hamiltonian, Eq. (16), arenot the same as the spinons of the exact treatment. Here fermion number is not conserved by the interac-tion vertices. In contrast spinon number is conserved bythe exact solution. Instead the fermions described by α † k should be viewed as propagating domain walls. 3. Properties of the Eigenstates Previously it has been suggested that the fundamentalexcitations of Heisenberg–Ising chains are chiral. Wenow make some remarks on this possibility, in light ofour results. The relevant chiral operator, C x , is definedas C x = ˆ x · X n S n × S n +1 = X n S yn S zn +1 − S zn S yn +1 . (29)For chirality to be a good quantum number for the XXZchain, C x must commute with the Hamiltonian. We find[ H, C x ] = X n i ( S xn − − S xn + S xn +1 − S xn +2 ) × ( S yn S yn +1 + S zn S zn +1 ) , (30)which is an O (1) contribution. A priori this demonstratesthat at least one eigenstate of the Hamiltonian is not aneigenstate of C x . However it is possible to show thatspinon states are generally not chirality eigenstates. First we write the commutator in the fermionic basis[ H, C x ] = − i X n (cid:0) c † n c n − M n − ,n +1 c † n c n M n − ,n +1 (cid:1) × (cid:0) M n − ,n − − M n − ,n + M n,n +1 − M n +1 ,n +2 (cid:1) , (31)where we have defined M n,n +1 = c † n c † n +1 + c n c n +1 − c † n c n +1 − c n c † n +1 . (32)As written, Eq. (31) contains only terms quartic and sex-tic in the creation and annihilation operators. Howeveronce Fourier transformed, rewritten in the Bogoliubovbasis ( α † k ) and normal ordered, quadratic terms will begenerated. A true (one–spinon) excitation of the system, | Ψ i , involves a superposition of domain walls created bythe α † k operators. Schematically | Ψ i = N X i =1 Y k , ··· ,k i l i ( k , · · · , k i )∆ − i α † k i | i . (33)where the l i ( k , · · · , k i ) are c –number functions of themomenta and we have written the small expansion pa-rameter ∆ − i explicitly. It has been shown that at low-est order the excitations are eigenstates of the chiral-ity operator. However from Eqs. (31) and (33) we seethat the expectation h Ψ | [ H, C x ] | Ψ i will generally not bezero, instead having a finite contribution at lower ordersin ∆ − . As ∆ → ∞ these contributions vanish and atthe Ising point the excitations can be chosen as chiralityeigenstates. IV. DYNAMICAL RESPONSE The quantity of interest for inelastic neutron scatteringis the dynamical structure factor S ab ( ω, Q ) given by S ab ( ω, Q ) = 1 N Z ∞−∞ dt π X l,l ′ e iωt e − iQ ( l − l ′ ) h S al ( t ) S bl ′ i . (34)Here S al = σ al is the a component of the spin operatorat site l . The structure factor is related to the retardeddynamical susceptibility χ abR ( ω, Q ) S ab ( ω, Q ) = − π − e − βω Im (cid:2) χ abR ( ω, Q ) (cid:3) , (35)where χ abR ( ω, Q ) = Z β dτ e iω n τ χ ab ( τ, Q ) (cid:12)(cid:12)(cid:12)(cid:12) ω n → η − iω ,χ ab ( τ, Q ) = − N X l,l ′ e − iQ ( l − l ′ ) h T τ S al ( τ ) S bl ′ i . (36)Here we have introduced the Matsubara formalismand the expectation implies a thermal trace h· · · i = Z − P m h m | e − βH · · ·| m i .Following the discussion in Sec. II it is apparent thatfor h = 0, χ ab ( ω, Q ) and hence S ab ( ω, Q ) will be diagonalin the indices a, b but that this is no longer the case for h > 0, in agreement with experiment. A. Dynamical Structure Factor in the FermionicRepresentation As previously discussed, the Jordan–Wigner transfor-mation introduces non–local ‘strings’ which make calcu-lating χ zz complicated. We instead focus our attentionon the transverse susceptibility, χ xx .First the required spin operator must be written interms of the new fermionic operators: σ xQ = 1 √ N X k e i Q (cid:2) k − θ k − θ k + Q + Q/ α † k α k + Q − i sin( k − θ k − θ k + Q + Q/ α † k α †− k − Q − α − k α k + Q ) (cid:3) . (37)We will use Eq. (37) to evaluate the time ordered dy-namical susceptibility, χ xx ( τ, Q ) = − (cid:10) T τ σ xQ ( τ ) σ x − Q (cid:11) . (38)As we aim to calculate (38) in perturbation theory in H we now switch to the interaction picture. In order tosimplify the perturbative calculation of χ xx it is usefulto express (38) in terms of a 3 × βγ ( τ, Q | k, k ′ )(the matrix indices take values β, γ = 1 , , 3) as follows χ xx ( τ, Q ) = 1 N X k,k ′ L β ( k )Π βγ ( τ, Q | k, k ′ ) L † γ ( k ′ ) , (39) where L β ( k ) =( i sin( γ k ) , cos( γ k ) , − i sin( γ k ) ) β , (40) γ k = k + Q/ − θ k − θ k + Q , (41)and the 3 × βγ ( τ, Q | k, k ′ ) is given byΠ βγ ( τ, Q | k, q ) = − (cid:10) T τ X ββ ( τ, Q | k ) X † γγ (0 , Q | q ) U ( β ) (cid:11) , (42) X βν ( τ, Q | k ) = α − k α k + Q α † k α k + Q 00 0 α † k α †− k − Q βν . (43)The imaginary time evolution operator in the interactionpicture is U ( τ ) = T τ exp (cid:18) − Z τ dτ H ( τ ) (cid:19) . (44)The Fourier transform of the matrix Π is given by Π ( iω n , Q | k, q ) = Z β dτ e iω n τ Π( τ, Q | k, q ) . (45) B. Zeroth Order At zeroth order in perturbation theory we replace U ( β )in (42) by 1. All off–diagonal elements then vanish andwe find ( ω nm = ω n − ω m )Π ( iω n , Q | k, q ) = 1 β X iω m G ( iω nm , k + Q ) G ( iω m , − k ) [ δ k, − q − Q − δ k,q ] = n k + n k + Q − iω n − ǫ k − ǫ k + Q [ δ k, − q − Q − δ k,q ] , (46)Π ( iω n , Q | k, q ) = 1 β X iω m G ( iω m , k + Q ) G ( − iω nm , k ) δ k,q = n k − n k + Q iω n + ǫ k − ǫ k + Q δ k,q , (47)Π ( iω n , Q | k, q ) = 1 β X ik n G ( − iω nm , − k − Q ) G ( − iω m , k ) [ δ k, − q − Q − δ k,q ] = n k + n k + Q − iω n + ǫ k + ǫ k + Q [ δ k,q − δ k, − q − Q ] . (48)Here n k = 1 / ( e βǫ k + 1) and the bare Green’s function isgiven by G ( ik n , k ) = 1 ik n − ǫ k . (49)The dynamical susceptibility at zeroth order in perturba-tion theory is then obtained by substituting the matrix Π into (39) and carrying out the momentum sums. Tak-ing the thermodynamic limit and analytically continuingto real frequencies, iω n → ω + iη , we arrive at the follow-ing expression for the zeroth order retarded susceptibility χ xxR, ( ω, Q ) = − Z π − π dk π (cid:20) (1 − cos(2 k + Q − θ k + θ k + Q ])) (cid:18) n k + n k + Q − ω + iη − ǫ k − ǫ k + Q + 1 − n k − n k + Q ω + iη + ǫ k + ǫ k + Q (cid:19) − k + Q − θ k + θ k + Q ))) n k + Q − n k ω + iη − ǫ k + ǫ k + Q (cid:21) . (50)The remaining k –integral cannot be carried out analyt-ically in general as the Bogoliubov parameters θ k needto be determined self–consistently and are therefore onlyknown implicitly. However, in the limit ∆ → ∞ the inte-gral can be taken and simple expressions for χ xxR, ( ω, Q )may be obtained. 1. The ∆ → ∞ , T = 0 Limit In order to evaluate the susceptibility further we take h = 0 and expand ǫ k as a series in 1 / ∆. This gives ǫ k = J ∆2 + J cos(2 k ) + J ∆ sin(2 k ) + . . . , (51) ǫ k + Q + ǫ k = J (∆ + 2 cos(2 k + Q ) cos( Q )) + . . . , (52) ǫ k + Q − ǫ k = − J sin(2 k + Q ) sin( Q ) + . . . . (53)We see that poles occur at ω = 2 J sin(2 k + Q ) sin( Q ) , (54) ω = J (∆ + 2 cos(2 k − + Q ) cos( Q )) , (55) ω = − J (∆ + 2 cos(2 k + + Q ) cos( Q )) . (56)The contribution from θ k only enters at O ( ) so we ne-glect it. We wish to take the imaginary part of the re-tarded susceptibility as this is proportional to the dy-namical structure factor. Defining P ( ω, E ) = Θ( | E | + ω )Θ( | E | − ω ) = (cid:26) | ω | > | E | | ω | ≤ | E | (57) where the Θ’s are Heaviside step functions, we find − Im χ xxR, ( Q, ω ) ≈ (1 + cos(2 k + Q )) ( n k + Q − n k ) P ( ω, J sin( Q ))16 J | cos(2 k + Q ) sin( Q ) |− 12 (1 − cos(2 k + + Q )) ( n k + + Q + n k + − J | sin(2 k + + Q ) cos( Q ) |× P ( ω + ∆ J, J cos( Q ))+ 12 (1 − cos(2 k − + Q )) ( n k − + Q + n k − − J | sin(2 k − + Q ) cos( Q ) |× P ( ω − ∆ J, J cos( Q )) , (58)with 2 k + Q = arcsin (cid:18) ω J sin( Q ) (cid:19) , (59)2 k + + Q = arccos (cid:20) − 12 cos( Q ) (cid:16) ωJ + ∆ (cid:17)(cid:21) , (60)2 k − + Q = arccos (cid:20) 12 cos( Q ) (cid:16) ωJ − ∆ (cid:17)(cid:21) . (61)Using the expansion of ǫ k one can also show that at nextto leading order n k + Q − n k = sinh( βω ) P ( ω, J sin( Q ))cosh( βω ) + cosh( β J ∆ − β cot( Q )((2 J sin( Q )) − ω ) / ) ,n k ± + Q + n k ± − − cosh( βω ) P ( ω ± J ∆ , J cos( Q ))cosh( βω ) + cosh( β tan( Q )((2 J cos( Q )) − ( ω ± J ∆) ) / ) . (62)Inserting these relations into the full expression yields − Im χ xxR, ( Q, ω ) ≈ (cid:20) J sin( Q )) − ω ) / − J sin( Q ) (cid:21) sinh( βω ) P ( ω, J sin( Q ))cosh( βω ) + cosh( β J ∆ − β cot( Q )((2 J sin( Q )) − ω ) / )+ 132 J cos( Q ) (cid:18) J cos( Q ) + ( ω + J ∆)2 J cos( Q ) − ( ω + J ∆) (cid:19) cosh( βω ) P ( ω + J ∆ , J cos( Q ))cosh( βω ) + cosh( β tan( Q )((2 J cos( Q )) − ( ω + J ∆) ) / ) − J cos( Q ) (cid:18) J cos( Q ) − ( ω − J ∆)2 J cos( Q ) + ( ω − J ∆) (cid:19) cosh( βω ) P ( ω − J ∆ , J cos( Q ))cosh( βω ) + cosh( β tan( Q )((2 J cos( Q )) − ( ω − J ∆) ) / ) . (63)This anaysis shows that, at zeroth order, the responsediverges as an inverse square root both for ω = ± (∆ J − J cos( Q )) and ω = ± J sin( Q ). These divergences areassociated with the gapped two–spinon response and thethermally activated Villain mode respectively. C. First Order Perturbation Theory At first order two classes of contributions appear,which may be distinguished by considering their diagram-matic representation. As we have already seen, the threezeroth order diagrams take the form of bubbles. Thefirst class of diagrams in first order consists of two bub-bles connected by an interaction vertex. The second classof diagrams consists of a single bubble with a self–energyinsertion. We detail these contributions in Appendix B.An important feature is that several of the first ordercontributions have stronger singularities as functions ofthe external frequency and momentum than the zerothorder results. This indicates that it is necessary for theexpansion to be resummed before useful physical resultscan be extracted. V. BUBBLE SUMMATION We have shown that at zeroth order the susceptibilitydiverges for certain ω and Q and that this is matchedby stronger divergences in some of the first order terms(see Appendix B). Moreover, it is clear from the firstorder calculation that higher orders in perturbation the-ory will exhibit stronger and stronger divergences. Inorder to get physically meaningful results we thereforeshould sum the most divergent classes of diagrams. Inthe case at hand, the complicated momenta dependenceof the vertices and the self–consistent Bogoliubov trans-formation makes the determination of the most divergentcontributions an impossible task. Instead we resum justthe connected bubble diagrams and justify our choice bycomparing results with the exact T = 0 calculation. Weexplain how to incorporate self–energy corrections in Sec.V D. A. RPA–Like Scheme The RPA scheme consists of carrying out bubble sumsof the type shown in Fig. 4 (cid:1) = (cid:2) + (cid:3) + (cid:4) + · · · + (cid:5) + · · · + (cid:6) + · · · FIG. 4: Bubble summation for one contribution to the dy-namical susceptibility matrix χ xx ( iω n , Q ). The thick lines in-dicate that the single particle propagators may be resummedto include self–energy corrections. To carry out these summations is non–trivial as the linesin the internal bubbles can take any orientation and themomentum dependence of the vertices is very compli-cated. In order to proceed we organize the interactionvertices into a 3 × V ( Q | k, q ) = V ( k + Q, − k, q, − q − Q ) ,V ( Q | k, q ) = − V ( q + Q, − q, k, − k − Q ) ,V ( Q | k, q ) = 6 V ( k + Q, − k, q, − q − Q ) ,V ( Q | k, q ) = 3 V ( k + Q, − k, q, − q − Q ) ,V ( Q | k, q ) = − V ( k + Q, q, − k, − q − Q ) ,V ( Q | k, q ) = − V ( k, q + Q, − q, − k − Q ) ,V ( Q | k, q ) = 6 V ( q + Q, − q, k, − k − Q ) ,V ( Q | k, q ) = 3 V ( q, k + Q, − k, − q − Q ) ,V ( Q | k, q ) = V ( q, − q − Q, k + Q, − k ) . (64)0As is shown in Appendix C, in the thermodynamic limitthe RPA–like bubble summation without taking into ac- count self–energy corrections results in an integral equa-tion of the formΠ RPA αβ ( iω n , Q | k, k ′ ) = Π αβ ( iω n , Q | k, k ′ ) + Z dq π K αγ ( iω n , Q | k, q )Π RPA γβ ( iω n , Q | q, k ′ ) . (65)where the kernel K is defined as the infinite volume limitof K αβ ( iω n , Q | k, q ) = X k ′ Π αγ ( iω n , Q | k, k ′ ) V γβ ( Q | k ′ , q ) . (66) Defining a convolution ∗ by( X ∗ Y ) αβ ( iω n , Q, k, k ′ ) ≡ Z π − π dq π X αγ ( iω n , Q | k, q ) Y γβ ( iω n , Q | q, k ′ ) , (67)this can be rewritten as (cid:0) I − K (cid:1) ∗ Π RPA = Π , (68)The integral equation (68) is then readily solved Π RPA = (cid:0) I − K (cid:1) − ∗ Π . (69)After analytic continuation to real frequencies iω n → ω + iη the quantity of interest is calculated as Z π − π dk dq (2 π ) L α ( k )Π RPA αβ ( ω + iη, Q | k, q ) L † β ( q ) . (70)In practice the solution of the integral equation discussedabove is reduced to a simple matrix inversion problem.We discretize all momentum integrals in terms of sumsover N = 400 points. We find that this value is largeenough to make discretization effects negligible. We set J = 1 and the regulator η = 10 − . The discretized rep-resentation of I − Π ∗ V (which is a matrix both inmomentum space as well as in the 3 × η is finite, thecalculated structure factor will be composed of a number ∼ N of Lorentzian peaks of width η . Finally we convolvethis result with a suitable Gaussian. B. Comparison with Exact Results for T = 0 Given the uncontrolled nature of our bubble summa-tion it is essential to compare it to exact results at zero ω I n t e n s it y ResummedExactIS Q=0, ∆=10, T=0 FIG. 5: (Color Online) The dynamical structure factor asfound by resummation, the exact result and the calcu-lation of Ishimura and Shiba (IS) at T = 0, Q = 0 and∆ = 10. In all cases the curves are convolved with a Gaus-sian in frequency space of full width half maximum 0.12. temperature in order to assess its quality. In Figs.5, 6, 7 and 8 we plot our results against the exact resultsfor the dynamical structure factor for ∆ = 10, T = 0and several momenta. We also include for comparison anearlier result for the DSF due to Ishimura and Shiba. We see that for ∆ = 10, T = 0 and at most wavevectors the resummation is a highly accurate approxima-tion. This suggests that the diagrams we have selectedaccount for most of the spectral weight at low temper-1 ω I n t e n s it y ResummedExactIS Q= π /4, ∆=10, T=0 FIG. 6: (Color Online) The dynamical structure factor asfound by resummation, the exact result and the calcu-lation of Ishimura and Shiba (IS) at T = 0, Q = π/ ω I n t e n s it y ResummedExactIS Q= π, ∆=10, T=0 FIG. 7: (Color Online) The dynamical structure factor asfound by resummation, the exact result and the calcu-lation of Ishimura and Shiba (IS) at T = 0, Q = π and∆ = 10. In all cases the curves are convolved with a Gaus-sian in frequency space of full width half maximum 0.12. ature. It is known from the exact result that thegapped response diverges along its lower energy thresh-old, in a region of momentum centred about π/ → π/ 2, although there isstill qualitative agreement (Fig. 8). For larger values of∆ our approximation becomes even better. ω I n t e n s it y ResummedExactExact Q=0 10 10.2 10.4 ω I n t e n s it y ResummedExact FIG. 8: (Color Online) A comparison of our calculation andthe exact result for Q = π/ T = 0 and ∆ = 10. Leftpanel: π/ Q = 0, in order to demon-strate the scale. Right panel: The same π/ Q = π/ C. Analytical Resummation for T = 0 It is possible and instructive to evaluate the resumma-tion exactly in the limit T = 0, ∆ → ∞ . At zeroth order,positive frequencies and T = 0, the only non–zero dia-gram is the particle–particle propagation bubble. Thisdiagram leads to a response in the region ω ∼ ∆ J . As T → / ∆ we find θ k = O (∆ − ) , (71)so that we can neglect the Bogoliubov phases. The dia-gram then reduces to Z dk π sin ( k + Q/ iω n − ǫ k − ǫ k + Q = A ( iω n , Q ) . (72)Using the approximation ǫ k + ǫ k + Q ≈ ∆ J + 2 J cos( Q ) cos(2 k + Q )we have A ( iω n , Q ) = Z π − π dk π sin ( k + Q/ iω n − ǫ k − ǫ k + Q = 14 J cos( Q ) Z π dk π − cos(2 k + Q )˜ ω − cos(2 k + Q )and ˜ ω = iω n − ∆ J J cos( Q ) . The remaining integral can be takenby standard contour integration methods and analyticcontinuation to real frequencies is then straightforward.With the approximations made here the vertex V takesa particularly simple form and resummation amounts to2the geometric series − Im[ χ xx ( ω, Q )] = − Im (cid:0) A ( ω, Q ) + 4 J cos( Q ) A ( ω, Q ) + (4 J ) cos ( Q ) A ( ω, Q ) + . . . (cid:1) = − Im( A ( ω, Q )1 − J cos( Q ) A ( ω, Q ) ) . (73) Finally the T = 0 result, to lowest order in ∆ is − Im[ χ xx ( ω, Q )] = ( √ (2 J cos( Q )) − ( ω − ∆ J ) J cos ( Q ) | ω − ∆ J | ≤ | J cos( Q ) | T = 0 usinga method based on perturbation theory combined with aLehmann representation. The expression obtained here(74) agrees to lowest order with their result. D. One–Loop Self–Energy Corrections As they stand, the first order tadpole contributions,Eqs. (B10–B16), cannot be incorporated into the RPA.This is because they contain divergences which are notresummed according to the scheme above. Instead wemust first calculate a new single particle Green’s func-tion, using a Dyson equation to resum the tadpole self–energy corrections. This Green’s function is then used tocalculate a new bubble diagram.Including the tadpole corrections generates anomalouspropagators at higher orders. These are most efficientlytaken into account by formulating the propagator as amatrix, g ( iω n , k ) = − Z β dτ e iω n τ g ( τ, k ) , g ( τ, k ) = (cid:28) T τ (cid:18) α k ( τ ) α † k (0) α k ( τ ) α − k (0) α †− k ( τ ) α † k (0) α †− k ( τ ) α − k (0) (cid:19) U ( β, (cid:29) (75)At zeroth order the result is g ( iω n , k ) = (cid:18) G ( iω n , k ) 00 − G ( − iω n , − k ) (cid:19) . (76)The single particle self–energy is defined by the Dysonequation g − ( iω n , k ) = g − ( iω n , k ) − Σ ( iω n , k ) . (77)If only ‘tadpole’ type diagrams are included the self–energy is, to first order in perturbation theory, frequencyindependent Σ = X p n p N (cid:18) V ( k, p, − p, − k ) − V ( p, − p, k, − k )3 V ( p, − p, k, − k ) − V ( k, p, − p, − k ) (cid:19) . (78) The elements of the matrix are O ( n k ∆ ) and hence theireffects will be most pronounced when the gap is small orthe temperature is large.Using the relations g ( iω n , k ) = 1 iω n − ǫ k , (79) g ( iω n , k ) = 1 iω n + ǫ k , (80)Σ = (cid:0) Σ (cid:1) ⋆ = − Σ , (81)Σ = − Σ , (82)one finds that g ( iω n , k ) = 1( iω n ) − ( ǫ k + Σ ) − | Σ | × (cid:18) iω n + ǫ k + Σ − Σ Σ iω n − ǫ k − Σ (cid:19) . We rewrite this as g ( iω n , k ) = h Z − k iω n − E k + Z + k iω n + E k i , (83) g ( iω n , k ) = h Z − k iω n + E k + Z + k iω n − E k i , (84) g ( iω n , k ) = λ k h iω n + E k − iω n − E k i , (85) g ( iω n , k ) = − g ( iω n , k ) , (86)with the definitions E k = p ( ǫ k + Σ ) + | Σ | , (87) Z ± k = 12 (cid:16) ∓ ǫ k + Σ E k (cid:17) , λ k = Σ E k . (88) E. Bubble Summation with Self–EnergyCorrections The one–loop self–energy corrections to the propaga-tors can be taken into acount in the bubble summation3for the dynamical susceptibility as follows. We define a3 × Π S byΠ Sβγ ( τ, Q | k, q ) = − (cid:10) T τ X ββ ( τ, Q | k ) X † γγ (0 , Q | q ) U ( β ) (cid:11)(cid:12)(cid:12)(cid:12) − loop Σ , (89)where only one–loop self–energy corrections are takeninto account. This amounts to calculating the two–point function of X and X † using the (anomalous) propagators(82). The elements of Π S are listed in Appendix D. Wenow follow the same steps as in Appendix C and showthat summing all bubble diagrams of the form shown inFig.4 gives rise to an integral equation obtained from (65)by the replacement Π αβ −→ Π Sαβ , i.e.Π RPA αβ ( iω n , Q | k, k ′ ) = Π Sαβ ( iω n , Q | k, k ′ ) + Z dq π K Sαγ ( iω n , Q | k, q )Π RPA γβ ( iω n , Q | q, k ′ ) . (90) π/2 π k E ( k ) T=0T=2JT=3JT=4JT=10J π/2 π π/2 π k -0.0500.05 Σ Im Σ T=2J 11 12 FIG. 9: (Color Online) The dispersion E k , including tadpoleself–energy corrections for ∆ = 10. The left panel showsthe gradual narrowing of the bandwidth with temperature.At this scale the dispersion for T = J is indistinguishablefrom that for T = 0. The right panel shows the two distinctelements of the self energy matrix, Σ sp at T = 2 J . where the kernel K S is defined as the infinite volumelimit of K Sαβ ( iω n , Q | k, q ) = X k ′ Π Sαγ ( iω n , Q | k, k ′ ) V γβ ( Q | k ′ , q ) . (91)At T = 0 the one–loop corrections to Σ vanish, so thatwe recover our previous result. On the other hand, for T > Σ plays an important role, altering the dispersionand shifting the thresholds of the dynamical response (seeFig. 9).One can also calculate the two–loop corrections to thesingle particle propagator but it is not as simple to in-corporate them into the RPA. For reference they are in-cluded in Appendix E. VI. RESULTS AND DISCUSSION We now turn to a discussion of our results at finite tem-peratures and applied fields. For the perturbative expan-sion to be valid we require ∆ ≫ 1. In selecting diagramsto resum we have favoured those that are most relevantat low temperatures, i.e. those that feature few thermaloccupation factors. The relevant energy scale is the sin-gle particle gap ∼ ∆ J/ 2, so we restrict our discussion totemperatures T < ∆ J/ 2. We calculate the transverse dy-namical structure factor as described in Sec. V A, usingthe matrix Π S found in Sec. V D and setting ∆ = 10, J = 1. For ∆ ≫ ω ∼ ∆ and aresponse for ω ∼ T = 0 should dis-appear. In our approach the thresholds are still presentalthough they are obfuscated by the necessity of convolv-ing the response with a Gaussian. This is a consequenceof the diagrams we have taken into account. We expectthis ‘thermal broadening’ to be a small effect that couldbe taken into account by including certain two–loop di-agrams. These two–loop diagrams connect the responseto decay channels of higher particle number. We discussthe issue further in Appendix E. A. Gapped Response for h = 0 We first consider the gapped (interband) response at fi-nite temperatures and for h = 0. At T = 0 the transverseDSF is dominated by a two–spinon scattering continuumthat occurs at energies around twice the single–particlegap, i.e. ω ∼ ∆ J . In Figs. 10, 11 and 12 we show howthe DSF in this regime of energies changes at finite tem-perature. The most striking feature is a narrowing of theresponse with increasing temperature. This is in agree-ment with inelastic neutron scattering experiments on4 ω I n t e n s it y T=0T=JT=2JT=3J Q=0 FIG. 10: (Color Online) The dynamical structure factor atfinite temperature for ω ∼ ∆ and ∆ = 10 at wave vector Q = 0. ω I n t e n s it y T=0T=JT=2JT=3J Q= π FIG. 11: (Color Online) The dynamical structure factor atfinite temperature for ω ∼ ∆ and ∆ = 10 at wave vector Q = π . TlCoCl . Another notable feature at T = 0 is that theresponse is not symmetric about Q = π/ As can beseen in Figs. 5–7 our calculation captures this behaviour.Accordingly the response develops asymmetrically withtemperature, and in a non–trivial way. In particular themaximum of the response for 0 ≤ Q < π/ π/ < Q ≤ π .For all wave vectors the total spectral weight in thegapped region decreases as temperature increases. Thethresholds of the response shift as temperature increasesdue to the thermal dependence of the single particle dis-persion E k , Eq. (87). This depletion of spectral weight isphysically sensible because the excitations are fermionic;as temperature increases more states are thermally oc- ω I n t e n s it y T=0T=JT=2JT=3J ω I n t e n s it y Q= π/4 Q= π/2 FIG. 12: (Color Online) The dynamical structure factor atfinite temperature for ω ∼ ∆ and ∆ = 10 at wave vectors Q = π/ π/ 2. Note that the scale of the intensity axesdiffers dramatically between the two plots. -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 ω I n t e n s it y Q= π/4 Q= π/2 Q= -2 -1 0 1 2 ω T=1JT=2JT=3J T=1J Q= π/2 FIG. 13: (Color Online) The dynamical structure factor for ω ∼ Q = π with temperature. cupied and concomitantly there are fewer states for thenew pair of fermions to fill. B. Villain Mode for h = 0 At temperatures greater than zero there is a thermalpopulation of spin excitations. Neutron scattering canthen lead to processes that do not change the spinonnumber of a given microstate, thus giving rise to an in-traband respose at low energies, ω ∼ 0. Dynamics ofthis kind were first described by Villain for the case of5 ω I n t e n s it y ω I n t e n s it y h=0h=J/2 Q=0 Q= π/2 FIG. 14: (Color Online) Comparison of the T = 0 dynami-cal structure factor for ∆ = 10 with and without an appliedtransverse field h = J/ 2. Main panel: the divergent responseat π/ Q = 0 is only apparent at the upper threshold. an finite length chain with an odd number of sites. Assuch processes rely on states being thermally occupied,their contribution to the DSF grows with temperature.The principle feature of the response is a well definedresonance or ‘mode’ at ω = 2 J sin( Q ). At lowest or-der in our calculation the intraband response exhibits asquare–root divergence (IV B 1). Taking interactions intoaccount by our bubble summation leads to a smoothingof the divergence, which however occurs in a very smallregion in energy, close to the threshold. This means thatthe zeroth order (especially when convolved with an ex-perimental resolution) is an excellent approximation tothe resummed result.The intraband response is also asymmetric about Q = π/ ≤ Q < π/ Q = π/ 2. For π/ < Q ≤ π it is enhanced. C. Response in a Transverse Field The transverse field h only enters the quadratic part ofthe Hamiltonian (16), hence its influence on the scatter-ing response is through the single particle dispersion andnot the interaction vertices. The first property to takeinto consideration is that h will have an effect on the ex-citation gap. As the magnitude of the field h approaches∆ J the gap collapses and the perturbative expansion isinapplicable. Hence we consider a field small comparedto ∆ J . For fields h ≫ ∆ J the gap opens again with thechain ferromagnetically ordered. A second important fea-ture to note is that a non–zero h causes the period of E k to double (see Fig. 1): the maximum at k = (2 n + 1) π is ω I n t e n s it y Q=0Q= π/4 Q= π/2 ω I n t e n s it y Q=3 π/4 Q= π h=J/2 FIG. 15: (Color Online) Wave vector dependence of the dy-namical structure factor for ∆ = 10 at h = J/ , T = 0. Themain panel shows wave vectors Q = 0 , π/ π/ 2. The insetshows plots for Q = 3 π/ π . -3 -2 -1 0 1 2 3 ω I n t e n s it y Q= π/4 Q= -3 -2 -1 0 1 2 3 ω I n t e n s it y h=0h=J h=J/2, T=J h=J/2, T=J, Q= π/2 FIG. 16: (Color Online) Wave vector dependence of the Vil-lain mode for ∆ = 10 at h = J/ , T = J . The main panelshows wave vectors Q = π/ , π/ π/ 4. The inset showsthe response at Q = π/ h = 0 and h = J/ reduced relative to those at k = 2 nπ . This leads to thedouble peaked structure seen at Q = π/ et al . The asymmetry of the response about Q = π/ T = 0 and ∆ = 10 showsthat the energy thresholds of the gapped response arevery nearly symmetric about π/ 2, this is no longer thecase in a finite field. Instead the upper threshold near Q = 0 is pushed to higher energies but is relatively un-changed near Q = π . We also note that the narrowing in6energy of the response near Q = π/ h = 0. Asymmetry is also seen in the thresholdsof the low energy response (Fig. 16). VII. COMPARISON TO DIAGONALIZATIONOF SHORT CHAINS AT T > As we have seen above, at zero temperature our ap-proach gives good agreement with the exact DSF. In or-der to assess the quality of our approximate DSF at finitetemperatures, we have computed the DSF by means ofnumerical diagonaliztion of the Hamiltonian on finite pe-riodic chains of up to 16 sites. The dynamical suscepti-bility (36) is expressed by means of a Lehmann expansionin terms of Hamiltonian eigenstates | n i of energy E n as χ xx ( ω, Q ) = N X n,m e − βE n − e − βE m ω + iη + E n − E m δ Q + p n ,p m × h n | S x | m i h m | S x | n i , (92)where the sums are taken over the eigenbasis of H . Fora sufficiently small system the eigenstates can be calcu-lated numerically. Due to the exponential increase inthe dimension of the Hilbert space with the number ofsites the method is restricted to short chains. We con-sider systems with N ≤ 16. In order to approximatethe thermodynamic spectrum from a finite number of al-lowed transitions, we take the parameter η in (92) to besufficiently large, so that the finite sum of Lorentziansin (92) becomes a smooth function of (the real part ofthe) frequency. Clearly this procedure can give a mean-ingful approximation of the susceptibility in the thermo-dynamic limit only if η is very small in comparison tothe thermal broadening. Hence the method is restrictedto sufficiently large temperatures. In order to obtain ameasure of the importance of finite–size effects we calcu-late the DSF for system sizes N = 8 , , 16. We find thatthe finite–size effects depend strongly on which region ofenergy and the momentum we consider.In Figs. 17–20 we compare the dynamical structurefactor calculated for a 16–site chain to the results ob-tained by our perturbative approach for several values ofmomentum. For Q = 0, π and T = 2 J (Fig. 17) theagreement of the two methods is good. The differenceat very small freqencies is probably due to finite–size ef-fects in the exact diagonalization results. The finite–sizeeffects for the main peak are found to be quite small.At Q = π/ T = 2 J (Fig. 18) the agreement ofthe two methods is less impressive. The disagreement athigh frequencies is likely due to inaccuracy of the bubblesummation in our perturbative method (we recall thatthe agreement of our resummation with the exact resultat T = 0 was worst for Q = π/ ω I n t e n s it y ResummedED FIG. 17: (Color Online) Comparison of the dynamical struc-ture factor obtained from the perturbative approach to exactdiagonalization on a 16–site chain with η = 0 . J for T = 2 J at momentum Q = 0. In order to facilitate a comparison theperturbative result has been convolved with a Lorentzian ofwidth η . ω FIG. 18: (Color Online) Comparison of the dynamical struc-ture factor obtained from the perturbative approach to exactdiagonalization on a 16–site chain with η = 0 . J for T = 2 J at momentum Q = π . In order to facilitate a comparison theperturbative result has been convolved with a Lorentzian ofwidth η . VIII. CONCLUSIONS We have calculated the transverse dynamical structurefactor of the XXZ spin chain at finite temperature andapplied field. The perturbative method we use is ac-curate at low temperatures and for large ∆ & 10. Inthis case the chain is in the Ising phase and the excita-tions are descended from propagating domain walls. Thescattering response is composed of two distinct parts for7 ω I n t e n s it y ω FIG. 19: (Color Online) Comparison of the dynamical struc-ture factor obtained from the perturbative approach to exactdiagonalization on a 16–site chain for T = 2 J at momentum Q = π/ 2. Left panel: low frequency response with η = 0 . J .Right panel: high frequency response with η = 0 . J . In or-der to facilitate a comparison the perturbative results havebeen convolved with a Lorentzian of width η . ω FIG. 20: (Color Online) Comparison of the dynamical struc-ture factor obtained from the perturbative approach to exactdiagonalization on a 16–site chain with η = 0 . J for T = 10 J at momentum Q = 0. In order to facilitate a comparison theperturbative result has been convolved with a Lorentzian ofwidth η . ∆ & 10, namely a gapped continuum at ω ∼ ∆ J andthe ‘Villain mode’ at ω ∼ 0. Our results pertain to thelow temperature and field dependence of both. Our mainobservations are1. The response associated with the gapped two–spinon continuum at T = 0 narrows in energy astemperature is increased and loses spectral weightto the emerging low frequency Villain mode. 2. The position of the peak (as a function of fre-quency) of the high frequency (gapped) continuumat T > Q = π/ ω ∼ Q = π/ It would be interesting to performquantitative comparisons with polarized neutron scatter-ing data. Acknowledgments We are grateful to Alan Tennant for getting us inter-ested in this problem and to Isaac P´erez–Castillo for hiscollaboration on the early part of this work and for manyhelpful discussions since. In particular we thank him forproviding efficient code to generate the exact solutioncurves at T = 0. We thank I. Affleck, S.E. Nagler andA.M. Tsvelik for useful discussions. This work was sup-ported by the EPSRC under grant EP/D050952/1. APPENDIX A: DIRECT JORDAN–WIGNERTRANSFORMATION OF THE HAMILTONIAN We may write the Hamiltonian (1) as H = H + H ,where H = J X n ∆ σ zn σ zn +1 + σ yn σ yn +1 + h X n σ xn ,H = J X n σ xn σ xn +1 . (A1) H can then be expressed as a quadratic form in spinlessfermions by means of the Jordan–Wigner transformation σ xn = c † n c n − ,σ zn − iσ yn − c † n e − iπ P j In this section we list the first order contributions tothe transverse response and discuss their divergent be-haviour. 1. ‘Connected Bubble’ diagrams The contributions to L ki Π kqij ( ω, Q ) R qj given by connect-ing two bubbles are (the box vertices in the diagrams in-dicate the inclusion of the external factors, L ki and R qj ): (cid:7) = 1 N X k,q sin( γ k ) sin( γ q ) V ( − k, k + Q, q, − q − Q ) × n k + n k + Q − iω n − ǫ k − ǫ k + Q n q + n q + Q − iω n − ǫ q − ǫ q + Q (B1) (cid:8) = 1 N X k,q sin( γ k ) sin( γ q ) V ( q, − q − Q, − k, k + Q ) × n k + n k + Q − iω n + ǫ k + ǫ k + Q n q + n q + Q − iω n + ǫ q + ǫ q + Q (B2) (cid:9) = 4 N X k,q cos( γ k ) cos( γ q ) V ( k + Q, q, − k, − q − Q ) × n k − n k + Q iω n + ǫ k − ǫ k + Q n q − n q + Q iω n + ǫ q − ǫ q + Q (B3) (cid:10) = 6 N X k,q sin( γ k ) sin( γ q ) V ( − k, k + Q, q, − q − Q ) × n k + n k + Q − iω n + ǫ k + ǫ k + Q n q + n q + Q − iω n − ǫ q − ǫ q + Q (B4) (cid:11) = 6 N X k,q sin( γ k ) sin( γ q ) V ( − k, k + Q, q, − q − Q ) × n k + n k + Q − iω n − ǫ k − ǫ k + Q n q + n q + Q − iω n + ǫ q + ǫ q + Q (B5) (cid:12) = 3 N X k,q cos( γ k ) sin( γ q ) iV ( k, q + Q, − q, − k − Q ) × n k − n k + Q iω n + ǫ k − ǫ k + Q n q + n q + Q − iω n + ǫ q + ǫ q + Q (B6) (cid:13) = − N X k,q sin( γ k ) cos( γ q ) iV ( − q, q + Q, − k − Q, k ) × n k + n k + Q − iω n − ǫ k − ǫ k + Q n q − n q + Q iω n − ǫ q + ǫ q + Q (B7) Æ q + Q − q = 3 N X k,q cos( γ k ) sin( γ q ) iV ∗ ( k + Q, − k, q, − q − Q ) × n k − n k + Q iω n + ǫ k − ǫ k + Q n q + n q + Q − iω n − ǫ q − ǫ q + Q (B8)9 (cid:15) = − N X k,q sin( γ k ) cos( γ q ) iV ∗ ( − q − Q, − k, k + Q, q ) × n k + n k + Q − iω n + ǫ k + ǫ k + Q n q − n q + Q iω n − ǫ q + ǫ q + Q (B9) 2. ‘Tadpole’ type diagrams These contributions consist of bubbles in which oneof the propagators features a tadpole type self energyinteraction: " (cid:16) + (cid:17) = 8 N X k,q (cos(2 γ k ) − n q V ( − k, q, − q, k ) iω n − ǫ k − ǫ k + Q × (cid:20) n k + n k + Q − iω n − ǫ k − ǫ k + Q − βn k (1 − n k ) (cid:21) (B10) " (cid:18) + (cid:19) = − N X k,q (cos(2 γ k ) − n q V ( k, q, − q, − k ) iω n + ǫ k + ǫ k + Q × (cid:20) n k + n k + Q − iω n + ǫ k + ǫ k + Q + βn k (1 − n k ) (cid:21) (B11) " (cid:20) + (cid:21) = 2 N X k,q (cos(2 γ k ) + 1) n q V ( k, q, − q, − k ) iω n + ǫ k − ǫ k + Q × (cid:20) n k − n k + Q iω n + ǫ k − ǫ k + Q + βn k (1 − n k ) (cid:21) − N X k,q (cos(2 γ k ) + 1) n q V ( k + Q, q, − q, − k − Q ) iω n + ǫ k − ǫ k + Q × (cid:20) n k − n k + Q iω n + ǫ k − ǫ k + Q + βn k + Q (1 − n k + Q ) (cid:21) (B12) " (cid:22) = − N X k,q sin(2 γ k ) iV ( q, − q, − k, k ) n q iω n + ǫ k − ǫ k + Q × (cid:20) n k + n k + Q − iω n − ǫ k − ǫ k + Q + 2 n ( ǫ k ) − ǫ k (cid:21) (B13) " (cid:23) = − N X k,q sin(2 γ k ) iV ( q, − q, − k, k ) n q iω n + ǫ k + ǫ k + Q × (cid:20) n k − n k + Q iω n − ǫ k + ǫ k + Q + 2 n ( ǫ k ) − ǫ k (cid:21) (B14) " (cid:24) = 3 N X k,q sin(2 γ k ) iV ( q, − q, − k, k ) n q iω n − ǫ k − ǫ k + Q × (cid:20) n k − n k + Q iω n + ǫ k − ǫ k + Q − n ( ǫ k ) − ǫ k (cid:21) (B15) " (cid:25) = − N X k,q sin(2 γ k ) iV ( q, − q, − k, k ) n q iω n − ǫ k + ǫ k + Q × (cid:20) n k + n k + Q − iω n + ǫ k + ǫ k + Q − n ( ǫ k ) − ǫ k (cid:21) (B16)One expects that, just as at zeroth order, these diagramswill have divergences at certain energies. Inspection ofthe pole structure of the connected bubble diagrams (B1–B9) suggests that they will have stronger divergencesthan the zeroth order because they feature products ofpoles. Numerical results support this assumption for Eqs.(B1,B2). On the other hand, it is clear that Eqs. (B4–B9) will be very small (for large ∆) because they containproducts of terms that, individually, are only significantfor different discrete regions in ω . The behaviour of Eq.(B3) is subtler. This diagram gives the most significantfirst order contribution to the response at ω ∼ 0. How-ever it does not contain a stronger divergence than itszeroth order equivalent. The reason is as follows: con-sider the sum X k I ( k, Q ) ω + ǫ k − ǫ k + Q + iη (B17)with I ( k, Q ) an analytic function of k, Q . As shown atzeroth order, because the dispersion ǫ k is bounded, the0imaginary part of the sum (as η → 0) will in turn bebounded. In general there is a divergence at the thresh-old ω = max( ǫ k − ǫ k + Q ) of the corresponding response.Naively one then expects that in the double sum X k,q I ( k, Q ) I ( q, Q ) I ′ ( k, q, Q )( ω + ǫ k − ǫ k + Q + iη )(( ω + ǫ k − ǫ k + Q + iη )) (B18)with I ′ ( k, q, Q ) analytic, the maximum contribution willoccur for k = q and ω = max( ǫ k − ǫ k + Q ). For Eq. (B3)the function V ( k + Q, k, − k, − k − Q ) vanishes at thethreshold ω = max( ǫ k − ǫ k + Q ). This means that the di-vergence is in fact substantially weaker than at zerothorder. This behaviour does not persist when self–energycorrections to the propagator are included. This is be-cause the self–energy corrections shift the thresholds ofthe response. Stronger than leading order divergencesare also found in Eqs. (B10–B16). To take them into ac-count the single particle propagator must be resummedusing a Dyson equation, as in Sec. V D. APPENDIX C: MATRIX STRUCTURE OF THEBUBBLE SUMMATION In this appendix we show that summing all bubble di-agrams results in the integral equation (65) for the ma- trix Π RPA . The proof follows by induction. Our startingpoint is expressions (39) and (42) for the dynamical sus-ceptibility. The n th order contribution to the matrix Πis by definition − (cid:10) T τ X ( τ, Q | k ) X † (0 , Q | q ) U ( n ) (cid:11) , (C1)where U ( n ) = ( − n n ! n Y m =1 Z β dτ m H ′ ( τ m ) . (C2)Out of all possible contractions in (C1) we want to selectonly those that give rise to bubble diagrams. We denotetheir contribution by Π RPA( n ) . We wish to show that Π RPA( n ) ( τ, Q | k, q ) = (cid:16) K ◦ K ◦· · ·◦ K ◦ Π (cid:17) ( τ, Q | k, q ) , (C3)where K is defined as K αβ ( τ, Q | k, q ) = X k ′ Π αγ ( τ, Q | k, k ′ ) V γβ ( Q | k ′ , q ) , (C4)and ◦ denotes a convolution and simultaneous matrixmultiplication( K ◦ K ) αβ ( τ, Q | k, q ) = 1 N X k ′ Z dτ K αγ ( τ − τ , Q | k, k ′ ) K γβ ( τ , Q | k ′ , q ) . (C5)Fourier transforming (C3) then gives Π RPA( n ) ( iω n , Q | k, q ) = (cid:16) K ∗ K ∗ · · · ∗ K ∗ Π (cid:17) ( iω n , Q | k, q ) , (C6)where the convolution ∗ is defined as( K ∗ K ) αβ ( iω n , Q | k, q ) = 1 N X k ′ K αγ ( iω n , Q | k, k ′ ) × K γβ ( iω n , Q | k ′ , q ) . (C7)Summing over n then leads to Π RPA ( iω n , Q | k, q ) = ∞ X n =0 Π RPA( n ) ( iω n , Q | k, q ) . (C8)This can be written as Π RPA = ∞ X n =0 K n ∗ Π = ( I − K ) − ∗ Π , (C9)where K n = n z }| { K ∗ K · · · ∗ K (C10) . The basic identity underlying the inductive proof of(C3) is X ( τ, Q | k ) U ( n ) (cid:12)(cid:12)(cid:12)(cid:12) RPA = ( K ◦ K ◦ · · · ◦ K ◦ X ) ( τ, Q | k ) , (C11)where the contraction notation indicates that only con-tractions compatible with our RPA-like summation havebeen carried out. Eqn (C3) is clearly a direct consequenceof (C11). We now prove (C11) by induction. For n = 1we prove by a lengthy but straightforward calculationthat X ( τ, Q | k ) U (1) (cid:12)(cid:12)(cid:12)(cid:12) RPA = ( K ◦ X ) ( τ, Q | k ) . (C12)The induction step is then straightforward. We have X ( τ, Q | k ) U ( n +1) (cid:12)(cid:12)(cid:12)(cid:12) RPA = − Z dτ n +1 Y H ′ ( τ n +1 ) (cid:12)(cid:12)(cid:12)(cid:12) RPA , (C13)1where Y = X ( τ, Q | k ) U ( n ) (cid:12)(cid:12)(cid:12)(cid:12) RPA . (C14)The combinatorial factor n + 1 cancels exactly againstthe 1 / ( n + 1) in the definition of U ( n +1) . Now using theinduction assumption (C11) in (C14), the final contrac-tion in (C13) reduces to the induction start n = 1, thusestablishing the validity of (C11) for n + 1. APPENDIX D: ELEMENTS OF Π S Introducing the definitions B −− ( k, k + Q ) = − n ( E k ) + n ( E k + Q ) − iω n − E k − E k + Q , (D1) B ++ ( k, k + Q ) = − − n ( E k ) − n ( E k + Q ) iω n + E k + E k + Q , (D2) B − + ( k, k + Q ) = − n ( E k ) − n ( E k + Q ) iω n − E k + E k + Q , (D3) B + − ( k, k + Q ) = − n ( E k + Q ) − n ( E k ) iω n + E k − E k + Q , (D4)the explicit elements of the matrix, Eq. (89), areΠ S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± Z σk Z σ ′ k + Q B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D5)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± σλ k ′ Z σ ′ k ′ + Q B σσ ′ ( k ′ , k ′ + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D6)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± σσ ′ λ k + Q λ k B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D7)Π S ( iω n , Q | k, k ′ ) = − X σ,σ ′ = ± σλ k Z σ ′ k + Q B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D8)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± (cid:0) Z σ ′ k + Q Z − σk δ k,k ′ + σσ ′ λ k λ k + Q δ k, − k ′ − Q (cid:1) B σσ ′ ( k, k + Q ) , (D9)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± σ ′ λ k + Q Z − σk B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D10)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± σσ ′ λ k + Q λ k B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D11)Π S ( iω n , Q | k, k ′ ) = − X σ,σ ′ = ± σ ′ λ k ′ + Q Z − σk ′ B σσ ′ ( k ′ , k ′ + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) , (D12)Π S ( iω n , Q | k, k ′ ) = X σ,σ ′ = ± Z − σk Z − σ ′ k + Q B σσ ′ ( k, k + Q ) (cid:0) δ k,k ′ − δ k, − k ′ − Q (cid:1) . (D13) APPENDIX E: FURTHER CONTRIBUTIONS TOTHE XXZ SPIN CHAIN RESPONSE The calculation described in the main text leads toa response that features sharp thresholds even at finitetemperature. Physically one expects these thresholdsto be absent at T = 0. The diagrams we considerare incapable of capturing this effect because they in-clude poles that only depend on two single particle en- ergies. This limits the response at positive frequenciesto the regions min( E k + E q ) ≤ ω ≤ max( E k + E q ) and0 ≤ ω ≤ max( E k − E q ), for general k, q . Coupling todecay channels involving more than two particles shouldalleviate this problem.One can evaluate the two–loop self energy correction2to the propagator. This includes diagrams of the form (cid:26) , (cid:27) , (cid:28) , . . . (E1)In analogy with the one–loop calculation, we define a self energy matrix Σ . As before the transferred frequencyand momentum are labelled by ik n and k respectively.Defining r = k + p + q , the elements of the matrix Σ ij areΣ = − N X p,q ˜ V ( k, − r, p, q ) ˜ V ( k, − r, p, q ) n ( ǫ q ) + n ( ǫ r ) − ik n + ǫ p + ǫ q + ǫ r ( n ( ǫ p ) − n B ( ǫ q + ǫ r ) − − N X p,q ˜ V ( k, − r, p, q ) ˜ V ( − q, − p, r, − k ) n ( ǫ q ) + n ( ǫ r ) − ik n − ǫ p − ǫ q − ǫ r ( n ( ǫ p ) − n B ( ǫ q + ǫ r ) − N X p,q ˜ V ⋆ ( r, − p, − q, − k ) ˜ V ( r, − k, − p, − q ) n ( ǫ q ) − n ( ǫ r ) ik n + ǫ p + ǫ q − ǫ r ( n ( ǫ p ) + n B ( ǫ r − ǫ q ))+ 12 N X p,q ˜ V ⋆ ( r, − p, − q, − k ) ˜ V ( k, − r, p, q ) n ( ǫ q ) − n ( ǫ r ) ik n − ǫ p − ǫ q + ǫ r ( n ( ǫ p ) + n B ( ǫ r − ǫ q )) . (E2)Σ = 96 N X p,q (cid:0) ˜ V ( k, − r, p, q ) (cid:1) n ( ǫ q ) + n ( ǫ r ) − ik n + ǫ p + ǫ q + ǫ r ( n ( ǫ p ) − n B ( ǫ q + ǫ r ) − N X p,q | ˜ V ( k, − r, p, q ) | n ( ǫ q ) + n ( ǫ r ) − ik n − ǫ p − ǫ q − ǫ r ( n ( ǫ p ) − n B ( ǫ q + ǫ r ) − N X p,q | ˜ V ( r, − q, − p, − k ) | n ( ǫ q ) − n ( ǫ r ) ik n + ǫ p + ǫ q − ǫ r ( n ( ǫ p ) + n B ( ǫ r − ǫ q ))+ 8 N X p,q (cid:0) ˜ V ( k, − r, p, q ) (cid:1) n ( ǫ q ) − n ( ǫ r ) ik n − ǫ p − ǫ q + ǫ r ( n ( ǫ p ) + n B ( ǫ r − ǫ q )) , (E3)Here we have used the boson occupation factor, n B ( ǫ ) =1 / (exp( βǫ ) − = − (Σ ) ⋆ , (E4)Σ = (Σ ) ⋆ . (E5)These contributions have poles that feature three single particle energies. Some of these contributions will lead toa temperature dependent broadening of the response inthe vicinities of ω ∼ ∆ J and ω ∼ 0. Including the two–loop self–energy terms should further increase the qualityof our approximation. However the extra sum over theloop momentum means that such a computation wouldbe a factor of N ∼ 400 slower. R. Orbach, Phys. Rev. , 309 (1958). J. des Cloizeaux and M. Gaudin, J. Math. Phys. , 1384(1966). C. N. Yang and C. P. Yang, Phys. Rev. , 321 (1966). V. E. Korepin, N. M. Bogoliubov and A. G. Izergin, Quan-tum Inverse Scattering Method and Correlation Functions ,Cambridge University Press, Cambridge, England (1993). H. A. Algra, L. J. de Jongh, H. W. J. Bl¨ote, W. J.Huiskamp, and R. L. Carlin, Physica (Utrecht) , 239(1976). P. M. Duxbury, J. Oitmaa, M. N. Barber, A. van der Bilt,K. O. Joung, and R. J. Carlin, Phys. Rev. B , 5149(1981). H. Yoshizawa, G. Shirane, H. Shiba, and K. Hirakawa, Phys. Rev. B , 3904 (1983). S. E. Nagler, W. J. L. Buyers, R. L. Armstrong, andB. Briat, Phys. Rev. Lett. , 590 (1982). S. E. Nagler, W. J. L. Buyers, R. L. Armstrong, andB. Briat, Phys. Rev. B , 3873 (1983). H. Yoshizawa, K. Hirakawa, S. K. Satija, and G. Shirane, Phys. Rev. B , 2298 (1981). A. Oosawa, K. Kakurai, Y. Nishiwaki, and T. Kato, J.Phys. Soc. Jpn. , 074719 (2006). L. D. Faddeev and L. A. Takhtajan, Phys. Lett. , 375(1981). H.-B. Braun, J. Kulda, B. Roessli, D. Visser, K. W.Kr¨amer, H.-U. G¨udel, and P. B¨oni, Nat. Phys. , 159(2005). I. Zaliznyak and S. Lee, Magnetic Neutron Scattering inModern Techniques for Characterizing Magnetic Materials,ed. Y. Zhu, Springer, Heidelberg (2005). J. Villain, Physica , 1 (1975). F. H. L. Essler and R. M. Konik, Phys. Rev. B ,100403(R) (2008). A. J. A. James, F. H. L. Essler and R. M. Konik, Phys.Rev. B , 094411 (2008). H. J. Mikeska and C. Luckmann, Phys. Rev. B , 184426(2006). D. A. Tennant et al., unpublished. N. Ishimura and H. Shiba, Prog. Theor. Phys. , 743(1980). J.-S. Caux, J. Mossel, and I. P´erez-Castillo, J. Stat. Mech. ,P08006 (2008). M. Jimbo and T. Miwa, Algebraic Analysis of Solvable Lat-tice Models , American Mathematical Society, Providence,Rhode Island (1995). A. H. Bougourzi, M. Karbach, and G. M¨uller, Phys. Rev.B , 11429 (1988). M. Takahashi, Prog. Theor. Phys. , 401 (1971); M. Taka-hashi and M. Suzuki, Prog. Theor. Phys. , 2187 (1972);C. Destri and H.J. de Vega, Phys. Rev. Lett. , 2313(1992); A. Kl¨umper, Z. Phys. , 507 (1993); C. Destri andH.J. de Vega, Nucl. Phys. B , 413 (1995); K. Fabricius, A. Kl¨umper and B.M. McCoy, Phys. Rev. Lett. , 5365(1999). F. H. L. Essler, H. Frahm, F. G¨ohmann, A. Kl¨umper,and V. E. Korepin, The One–Dimensional Hubbard Model ,Cambridge University Press, Cambridge (2005). M. Takahashi, Thermodynamics of One–Dimensional Solv-able Models , Cambridge University Press, Cambridge, Eng-land (2005). J. D. Johnson, J. Appl. Phys. , 1991 (1981). H. J. Mikeska, Phys. Rev. B , 2794 (1975). N. Kitanine, J. M. Maillet, N. A. Slavnov and V. Terras, Nucl. Phys. B , 558 (2005) F. G¨ohmann, A. Kl¨umper and A. Seel, J. Phys. A: Math.Gen. , 7625 (2004). F. G¨ohmann, N. P. Hasenclever and A. Seel, J. Stat. Mech. P10015 (2005). K. Sakai, J. Phys. A: Math. Theor. , 7523 (2007). K. Fabricius, U. L¨ow and J. Stolze, Phys. Rev. B , 5833(1997); S. Grossjohann and W. Brenig, arXiv:0811.1956;T. Barthel, U. Schollw¨ock and S. White, arXiv:0901.2342. P. Jordan and E. Wigner, Z. Phys. , 631 (1928). H. A. Kramers and G. H. Wannier, Physical Review ,252 (1941). G. G´omez-Santos, Phys. Rev. B , 6788 (1990). D.V. Dmitriev, V.Y. Krivnov and A.A. Ovchinnikov, Phys.Rev. B , 172409 (2002); D.V. Dmitriev, V.Ya. Krivnov,A.A. Ovchinnikov and A. Langari, JETP , 538 (2002);J.-S. Caux, F.H.L. Essler and U. L¨ow, Phys. Rev. B ,134431 (2003); F. Capraro and C. Gros, Eur. Phys. J. B29 ,35 (2002). R. Hagemans, J.S. Caux and U. L¨ow, Phys. Rev.B , 014437 (2005). J. B. Kogut, Rev. Mod. Phys.51