Jun 24, 2008 - can be calculated exactly with the Lanczos algorithm [17]. ..... [17] Cornelius Lanczos, An iteration method for the solution of the ei...

0 downloads 9 Views 304KB Size

No documents

arXiv:0806.3953v1 [cond-mat.str-el] 24 Jun 2008

Herman E. Boos∗ , Jens Damerau‡ , Frank G¨ohmann†, Andreas Kl¨umper§ Fachbereich C – Physik, Bergische Universit¨at Wuppertal, 42097 Wuppertal, Germany Junji Suzuki¶ Department of Physics, Faculty of Science, Shizuoka University, Ohya 836, Suruga, Shizuoka, Japan Alexander Weißek Institut f¨ur Physik, Universit¨at Greifswald, 17487 Greifswald, Germany

Abstract Recent studies have revealed much of the mathematical structure of the static correlation functions of the XXZ chain. Here we use the results of those studies in order to work out explicit examples of short-distance correlation functions in the infinite chain. We compute two-point functions ranging over 2, 3 and 4 lattice sites as functions of the temperature and the magnetic field for various anisotropies in the massless regime −1 < ∆ < 1. It turns out that the new formulae are numerically efficient and allow us to obtain the correlations functions over the full parameter range with arbitrary precision. PACS: 05.30.-d, 75.10.Pq

∗ e-mail: ‡ e-mail: † e-mail: § e-mail: ¶ e-mail: k e-mail:

[email protected] [email protected] [email protected] [email protected] [email protected] [email protected]

2

1 Introduction In recent years considerable progress was achieved in the understanding of the static correlation functions of local operators in the XXZ spin chain. This concerns mostly the ground state correlation functions at zero magnetic field for which a hidden fermionic structure (outside the free Fermion point ∆ = 0) was identified in [4–6]. The exponential of an operator, bilinear in two types of annihilation operators and acting on the space of local operators on the spin chain, generates all local correlation functions, while the space of local operators on the chain is generated by the corresponding creation operators. For the inhomogeneous chain this structure implies the factorization of arbitrary correlation functions into sums over products of neighbour correlators. A generalization of the exponential form to finite temperatures and finite longitudinal magnetic field was suggested in [2]. Without magnetic field the same Fermi operators as in the ground state can still be used, and the modification of the ground state result concerns only certain functions multiplying the annihilation operators in the exponent. For finite magnetic field, on the other hand, it seems that the algebraic structure must be altered as well and a new Fermi-like operator appears linearly in the exponent. Here we shall work out the finite temperature formula of [2] in detail and compute the longitudinal and transversal spin-spin correlation functions for up to four lattice sites. We shall concentrate on the massless regime −1 < ∆ < 1, since the massive regime ∆ > 1 requires the calculation of different functions in the analytical part of our work, which we postpone to a separate publication.

2 The density matrix of a chain segment We consider the XXZ chain with Hamiltonian N

HN = J

∑

j=−N+1

y y σxj−1 σxj + σ j−1 σ j + ∆(σzj−1 σzj − 1)

(1)

in the antiferromagnetic regime (J > 0 and ∆ > −1)† . Here the operators σαj , j = −N + 1, . . . , N, α = x, y, z, act locally as Pauli matrices, J is the exchange coupling, and ∆ = ch(η) = 12 (q + q−1 ) is the anisotropy parameter. The XXZ Hamiltonian preserves the z-component of the total spin SzN

=

1 2

N

∑

σzj .

(2)

j=−N+1

results of this and of the following section are valid for all ∆ > −1. Later on in sections 4 and 5 we restrict ourselves to the massless regime −1 < ∆ < 1. † The

3 The equilibrium properties of the chain at temperature T and external magnetic field h are determined by the statistical operator z

e−(HN −hSN )/T . ρN (T, h) = z tr−N+1,...,N e−(HN −hSN )/T

(3)

We are interested in the system in the thermodynamic limit L = 2N → ∞. To perform this limit in a sensible way we fix a positive integer n and define Dn (T, h) = lim tr−N+1,...,−1,0,n+1,n+2,...,N ρN (T, h) , N→∞

(4)

the density matrix of the segment [1, n] of the infinite chain. By construction it satisfies the reduction property trn Dn (T, h) = Dn−1 (T, h) . (5) Denote by W the space of operators with non-trivial action only on a finite number of positive lattice sites. Because of the translational invariance of the Hamiltonian it is sufficient to consider this space. For O ∈ W we write its restriction to the first n lattice sites as O[1,n] . Then D∗T,h (O) = lim tr1,...,n Dn (T, h)O[1,n] (6) n→∞

is well defined because of (5) and determines the thermal average of the operator O in a magnetic field h, which we shall also denote hOiT,h. Mathematically (6) defines D∗T,h as a functional on W. Along with D∗T,h we often study its inhomogeous generalization (see [2]) which we denote by the same letter. When acting on an operator of length n it depends on n so-called inhomogenieties ξ j , j = 1, . . ., n. In our previous work [2] we conjectured a rather explicit formula for D∗T,h which was inspired in part by [5]. The conjecture says that D∗T,h (O) = tr eΩ1 +Ω2 (O) , (7) where tr is the trace functional defined by tr(O) = . . . 21 tr1

1 1 2 tr2 2 tr3 . . . (O) ,

(8)

and the Ω j : W → W, j = 1, 2 are linear operators of the form‡

Ω2 = − lim

Z

α→0 Γ

α→0

‡ We

dζ21 dζ22 ω(µ1 , µ2 ; α)b(ζ1 , α)c(ζ2 , α − 1) , 2 2 Γ 2πiζ1 2πiζ2

Ω1 = − lim

Z Z

dζ21 ϕ(µ1 ; α)h(ζ1 , α) . Γ 2πi

(9a) (9b)

have slightly changed our notation for operators. Instead of those of [2] we rather employ the conventions of the more recent paper [6].

4 Here µ j = ln ζ j , j = 1, 2. As a consequence of taking the limit α → 0 we shall only need the values of the function ω(µ1 , µ2 ; α) and of its α-derivative at α = 0. Similarly we only need ϕ(µ1 ; 0). All the functions that are necessary for the actual evaluation of (9) will be introduced in section 4 below. They depend on the temperature and on the magnetic field and describe what we call the physical part of the problem. By way of contrast, the operators b(ζ, α), c(ζ, α) and h(ζ, α) are independent of T and h. These operators are rational functions of ζ2 and constitute the algebraic part of the problem. We will explain their construction in the next section. The closed contour Γ encircles the poles of the operators at the inhomogeneities ξ2j (see below) and excludes all other singularities of the integrand.

3 The algebraic part of the construction We should first explain that our present status of understanding is different for Ω1 and Ω2 in (9). The algebraic part of Ω1 , built up from the operators b(ζ, α) and c(ζ, α), is the same as for the ground state at vanishing magnetic field. First versions of these operators were introduced in [5]. They were optimized and further studied in [4, 6]. Their analytic and algebraic structure is by now well understood. In particular, they satisfy a reduction property and anticommutation relations. The operator h(ζ, α) was introduced in [2] in order to account for correlation functions odd under spin reversal. This operator is less well studied. We plan to further elaborate on its mathematical structure elsewhere. Below, when we work out examples of correlation functions, we use some of its properties for α → 0 that we verified by explicit calculations for up to four lattice sites. In this section we follow to a certain extent the recent paper [6] to which we refer the reader for further details. For the definition of the operators b(ζ, α), c(ζ, α) we first of all generalize the space of local operators. Instead of local operators we consider operators of the form q

2αS(0)

O,

S(0) =

1 2

0

∑

σzj ,

(10)

j=−∞

where O is local and α ∈ C. These operators are called quasi-local with tail α. They span a vector space Wα . The parameter α is introduced in order to regularize the various objects introduced below. Some of them are rather singular for α → 0, and the limit is, in general, only well defined for appropriate combinations. In [6] the operators b(ζ, α), c(ζ, α) were defined in two steps. In step one endomorphisms b[kl] (ζ, α) and c[kl] (ζ, α) acting on M[k,l] = End(Vk ⊗ · · · ⊗Vl ) were defined, where every V j ∼ = C2 is the space of states of a spin in the infinite chain. In step two a certain reduction property was used to extend their action inductively to Wα . Here we shall only review step one, since we will be dealing with small finite chain segments anyway, and we refer the reader to [6] for the reduction property. The operators b[kl] (ζ, α) and c[kl] (ζ, α) are constructed from weighted traces of the

5 c2 ). These monodromy maelements of certain monodromy matrices related to Uq (sl trices are products of two types of L-matrices with two-dimensional and with infinite dimensional auxiliary space, respectively. The L-matrices with two-dimensional auxiliary space are directly related to the R-matrix of the six-vertex model, 1 0 0 0 0 β(ζ) γ(ζ) 0 (11) R(ζ) = 0 γ(ζ) β(ζ) 0 , 0 0 0 1 where β(ζ) =

(1 − ζ2 )q , 1 − q2 ζ2

γ(ζ) =

(1 − q2 )ζ . 1 − q2 ζ2

(12)

Fixing an auxiliary space Va isomorphic to C2 we define La, j (ζ) = ρ(ζ)Ra, j (ζ), where ρ(ζ) is a scalar factor to be specified below. This is the standard L-matrix of the sixvertex model. The corresponding monodromy matrix is Ta,[k,l] (ζ) = La,k (ζ/ξk ) . . . La,l (ζ/ξl ) .

(13)

It acts on Va ⊗ Vk ⊗ · · · ⊗ Vl . We are interested in operators acting on M[k,l] or on Ma ⊗ M[k,l] , where Ma = End(Va ). Such type of operators are naturally given by the adjoint action of operators acting on Vk ⊗ · · · ⊗Vl or on Va ⊗Vk ⊗ · · · ⊗Vl . An example is the L-operator La, j (ζ) defined by La, j (ζ) X[k,l] = La, j (ζ)X[k,l] L−1 a, j (ζ)

(14)

for j ∈ {k, . . ., l} and for all X[k,l] ∈ Ma ⊗ M[k,l] . These L-operators generate the first type of monodromy matrices entering the constructions of b(ζ, α), c(ζ, α), Ta (ζ, α)X[k,l] = La,k (ζ/ξk ) . . . La,l (ζ/ξl )qασ X[k,l] z

(15)

for all X[k,l] ∈ Ma ⊗ M[k,l] . A second type of monodromy matrix acts adjointly on an infinite dimensional bosonic auxiliary space, more precisely on a representation space of the q-oscillator algebra Osc defined in terms of generators a, a∗ , q±D and relations qD aq−D = q−1 a ,

qD a∗ q−D = q a∗ ,

aa∗ = 1 − q2D+2 ,

a∗ a = 1 − q2D .

Let OscA a copy of Osc. We define −D q A 0 1 − ζ2 q2DA +2 −ζaA LA, j (ζ) = σ(ζ) 1 −ζa∗A 0 q DA j j

(16)

(17)

6 which acts on OscA ⊗Vk ⊗· · ·⊗Vl . The scalar factor σ(ζ) is a convenient normalization and is required to solve the functional equation σ(ζ)σ(q−1ζ) =

1 . 1 − ζ2

(18)

It also fixes the scalar factor in the definition of La, j (ζ), ρ(ζ) =

q−1/2 σ(q−1 ζ) . σ(ζ)

(19)

For the L-operator (17) we define again its adjoint action LA, j (ζ) X[k,l] = LA, j (ζ)X[k,l]L−1 A, j (ζ) ,

(20)

for j ∈ {k, . . ., l} and for all X[k,l] ∈ End(OscA ) ⊗ M[k,l] , and also the corresponding monodromy matrix TA (ζ, α)X[k,l] = LA,k (ζ/ξk ) . . . LA,l (ζ/ξl )q2αDA X[k,l] .

(21)

The monodromy matrices (15) and (21) are the basic ingredients in the definition of the fermionic operators entering the exponential form of the density matrix. In addition to the monodromy matrices we shall need certain operators acting on the spin. We introduce the adjoint action of the z-component of the total spin on M[k,l] , z S X[k,l] = [S[k,l] , X[k,l]] ,

z S[k,l] =

1 2

l

∑ σzj ,

(22)

j=k

and the spin reversal operator defined by J X[k,l] =

h

i h l i x x X σ σ ∏ j [k,l] ∏ j . l

j=k

(23)

j=k

Somewhat loosely we say that an operator X[k,l] has spin s if S X[k,l] = sX[k,l] .

(24)

We call X[k,l] spin-reversal even if J X[k,l] = X[k,l] and spin-reversal odd if J X[k,l] = −X[k,l] . The annihilation operators b(ζ, α) and c(ζ, α) as well as a couple of other interesting fermionic operators (see [6]) can be derived from a master object k(ζ, α) defined as n z o α−S −2S[k,l] k(ζ, α) X[k,l] = trA,a σ+ T (ζ, α)T (ζ, α)ζ q X . (25) a A [k,l] a Here the bosonic trace is understood in such a way that it agrees with the trace in the bosonic representation W + = ⊕k≥0 C|0i for |q| ≤ 1 but not a root of unity (with

7 qD |ki = qk |ki, a|ki = (1 − q2k )|k − 1i, a∗ |ki = |k + 1i). For more details see appendix A of [6]. We shall give a description of the operators b(ζ, α) and c(ζ, α) in terms of certain residua of k(ζ, α) which is appropriate for implementing them on a computer. Let us fix an operator X[k,l] of spin s. Note that |s| ≤ l − k + 1 = n. The following statements about the analytic structure of k(ζ, α) were worked out in [6]. (i) ζ−α+s−1 k(ζ, α)X[k,l] is rational in ζ2 , (ii) Its only singularities in the finite complex plane are poles at ζ2 = ξ2j , q±2 ξ2j , (iii) ζ−α−s+1 k(ζ, α)X[k,l] is regular at ζ2 = ∞. It follows that kskal (ζ, α)X[k,l] = ζ−α−s−1 k(ζ, α)X[k,l] is regular at ζ2 = 0 if s ≤ 0 and has at most a pole of order s in ζ2 at ζ2 = 0 if s > 0. Furthermore kskal (ζ, α)X[k,l] ∼ ζ−2 for ζ → ∞. Hence, kskal (ζ, α)X[k,l] =

(ε)

ρ j (α)

n

∑ ∑

2 2ε 2 j=1 ε=0,± ζ − q ξ j

s

+∑ζ j=1

−2 j

κ j (α) X[k,l]

(26) (ε)

(where the second sum over j is zero by definition for s < 1). The residua ρ j (α) and the Laurent coefficients κ j (α) determine all those operators that are needed for vanishing magnetic field. They are defined by equation (26). First of all there are the operators c¯ (ζ, α), c(ζ, α) and f(ζ, α) defined in section 2.7 of [6]. In terms of the residua and Laurent coefficients they read α+s+1

c¯ (ζ, α)X[k,l] = ζ

α+s+1

c(ζ, α)X[k,l] = ζ

α+s+1

f(ζ, α)X[k,l] = ζ

(0)

ζ2 + ξ2j ρ j (α)X[k,l] ∑ 2 2 · 2ξ2 , j=1 ζ − ξ j j n

(27a)

ζ2 + ξ2j qε(α+s−1) ρ j (α)X[k,l] , ∑ ∑ 2 2· 4ξ2j j=1 ε=± ζ − ξ j (ε)

n

(27b)

ζ−2 j κ j (α)X[k,l] ∑ α+s+1−2 j − q−α−s−1+2 j j=0 q s

(ε) ζ2 + ξ2j ε qε(α+s−1) ρ j (α)X[k,l] −∑ ∑ 2 · , 2 4ξ2j j=1 ε=± ζ − ξ j n

(27c)

where by definition n

κ0 = − ∑

∑

j=1 ε=0,±

(ε)

ρ j (α) 2q2ε ξ2j

.

(28)

In terms of these operators k(ζ, α) is decomposed as k(ζ, α)X[k,l] = c¯ (ζ, α) + c(qζ, α) + c(q−1ζ, α) + f(qζ, α) − f(q−1ζ, α) X[k,l] . (29)

8 The operator b(ζ, α) that is needed in the construction of the density matrix is defined as b(ζ, α) = q−1 (qα−S −1 − q−α+S +1 ) J c(ζ, −α) J , (30) where J is the spin-flip in the adjoint representation. Taking into account that J S J = − S we find b(ζ, α)X[k,l] = ζ where

(ε)

ρ j (α)X[k,l] ζ2 + ξ2j q−ε(α+s+1) e , ∑ ∑ ζ2 − ξ2 · 4ξ2j j=1 ε=± j n

−α−s+1

(31)

(ε) (ε) e ρ j (α)X[k,l] = q−1 (qα−s − qs−α ) J ρ j (−α) J X[k,l] ,

(32)

for ε = ±. For later convenience we introduce the operators b j (α) =

∑

ε=±

(ε)

q−ε(α+S +2) e ρ j (α) 2ξ2j

,

c j (α) =

∑

qε(α+S −2) ρ j (α)

ε=±

(ε)

2ξ2j

,

(33)

for j = 1, . . ., n. In terms of these operators we have 2 2 1 ζ +ξj · b j (α)X[k,l] , 2 − ξ2 2 ζ j j=1 n

b(ζ, α)X[k,l] = ζ−α−s+1 ∑ α+s+1

c(ζ, α)X[k,l] = ζ

2 2 1 ζ +ξj ∑ 2 2 · c j (α)X[k,l] . j=1 2 ζ − ξ j

(34a)

n

(34b)

Let ω(µ1 , µ2 ; α) denote the function that determines the physical part of the problem. ω(µ1 , µ2 ; α) is a known function of the temperature and the magnetic field (see section 3 of [2]). Here we only need its property ω(µ2 , µ1 ; α) = ω(µ1 , µ2 ; −α). We e (µ1 , µ2 ; α) = ω(µ1 , µ2 ; α)(ζ2 /ζ1 )α (recall that ζ j = eµ j ) which has the same also set ω symmetry. Then, by definition, dζ21 dζ22 Ω1 (α)X[k,l] = − ω(µ1 , µ2 ; α)b(ζ1 , α)c(ζ2 , α − 1)X[k,l] 2 2 Γ 2πiζ1 Γ 2πiζ2 e (λi , λ j ; α)bi (α)c j (α − 1) + ω e (λ j , λi ; α)b j (α)ci (α − 1) X[k,l] (35) =− ∑ ω Z

Z

1≤i< j≤n

(where λ j = ln ξ j ) for all X[k,l] with S(X[k,l]) = 0. The closed contour Γ in (35) encircles all the ξ2j (and excludes the singularities of ω). In the second equation we used that bi (α)ci (α − 1) = 0 ,

for i = 1, . . . , n .

(36)

9 For the physical correlation functions the limit α → 0 remains to be done. We shall write Ω1 = Ω1 (0). We claim that ′ − (37) Ω1 X[k,l] = − ∑ ω(λi , λ j )Ω+ i, j + ω (λi , λ j )Ωi, j X[k,l] , 1≤i< j≤n

where

Ω+ i, j = lim bi (α)c j (α − 1) + b j (α)ci (α − 1) , α→0

Ω− = lim α b (α)c (α − 1) − b (α)c (α − 1) , i j j i i, j α→0

and

e (λi , λ j ; 0) , ω(λi , λ j ) = ω

(38a) (38b)

e (λi , λ j ; α) α=0 . ω′ (λi , λ j ) = ∂α ω

(39)

The existence of the limits in (38) was proved in [2]. As was already mentioned in the introduction to this section, the operator h(ζ, α) introduced in [2] is a less well understood than b(ζ, α) and c(ζ, α). Here we adapt the definition of [2] to the novel conventions of [6] and set n z o α−S −2S[k,l] h(ζ, α) X[k,l] = (1 − qα ) trA,a σ+ T (ζ, α)T (ζ, α)ζ q a X . (40) A [k,l] a a

This operator has some similarity with the operator k(ζ, α) introduced above. The only differences are the prefactor (1 − qα ) and the insertion of an extra bosonic annihilation operator on the right hand side. For this reason its analytic structure as a function of ζ is rather similar to that of k(ζ, α). Following the same reasoning as in section 2.5 of [6] and using that [S, h(ζ, α)] = 0 one sees that for an operator X[k,l] of spin s (i) ζ−α+s h(ζ, α)X[k,l] is rational in ζ2 , (ii) its only singularities in the finite complex plane are poles at ζ2 = ξ2j , q±2 ξ2j , (iii) ζ−α−s+2 h(ζ, α)X[k,l] is regular at ζ2 = ∞. It follows that h(ζ, α)X[k,l] is of the form " # (ε) n s h (α) h (α) j j h(ζ, α)X[k,l] = ζα+s ∑ 2 +∑ 2 + ∑ ζ−2 j θ j (α) X[k,l] . (41) 2 2 2ε ε=± ζ − q ξ j j=1 ζ − ξ j j=1

We denote the function that determines the physical part of Ω2 by ϕ(µ; α) (see [2]) e(µ; α) = ϕ(µ; α)ζα. Then, for every X[k,l] of spin 0, and set ϕ Ω2 X[k,l] = − lim

α→0

dζ21 ϕ(µ1 ; α)h(ζ1 , α)X[k,l] = Γ 2πi

Z

n

n

e(λ j ; α)h j (α)X[k,l] = − ∑ ϕ(λ j )h j X[k,l] , (42) − lim ∑ ϕ α→0 j=1

j=1

10 where h j = lim h j (α) , α→0

e(λ j ; 0) . ϕ(λ j ) = ϕ

(43)

Since only the residua at the simple poles at ξ2j enter the formula for the correlation functions, there is a certain degree of arbitrariness in the definition of the operator h(ζ, α). Equations (37) and (42) can be used in order to calculate short-distance correlators explicitly on a computer. Examples will be shown in section 5 below. The problem clearly divides into two separate tasks. The first task is the calculation of Ω± i, j and h j . For this purpose one has to generate the operators k(ζ, α) and h(ζ, α) for a finite number of lattice sites and then calculate their residua. The second task is the actual calculation of the functions ω(µ1 , µ2 ), ω′ (µ1 , µ2 ) and ϕ(µ). They will be described in the following section. Still, the knowledge of the operators Ω1 and Ω2 in the form (37) and (42) would not be very useful, if these operators would not be nilpotent, which causes the exponential series in (9) to terminate after finitely many terms. As was first shown in [3] we have − − + + − Ω+ i, j Ωk,l = Ωi, j Ωk,l = Ωi, j Ωk,l = 0 ,

if {i, j} ∩ {k, l} = 6 0/ ,

− − − + + [Ω+ i, j , Ωk,l ] = [Ωi, j , Ωk,l ] = [Ωi, j , Ωk,l ] = 0 ,

if {i, j} ∩ {k, l} = 0/

(44a) (44b)

which is a consequence of the fermionic nature of the operators b(ζ, α) and c(ζ, α). It follows for chain segments of length n that [ n ]+1 Ω12 = 0,

(45)

where [x] stands for the integer part of x. Again for Ω2 our knowledge is so far re− stricted to what be learned from explicit examples. We calculated Ω+ i, j , Ωi, j and h j explicitly in the inhomogeneous case for n = 1, 2, 3, 4. The reader might wonder why we do not present the matrices here. The reason is that they are very big. Already for n = 4 they are so big that they can not be handled any longer by an all-purpose computer algebra program like Mathematica. We used the program FORM [19] instead which turned out to be efficient. We found that for n up to 4 the h j mutually anticommute h j hk + hk h j = 0 (46) and that they commute with the Ω± i, j , [Ω± i, j , hk ] = 0 .

(47)

± Ω± i, j hk = hk Ωi, j = 0 if k ∈ {i, j}.

(48)

We also have the stronger property

We conjecture that (46)-(48) are true in general. This would mean, in particular, that [Ω1 , Ω2 ] = 0 (49)

11 and that Ω2 is nilpotent,

Ω22 = 0 .

(50)

As a further consequence the exponential form (9) would factorize, eΩ1 +Ω2 = eΩ1 eΩ2 = eΩ1 (1 + Ω2 ) = (1 + Ω2 )eΩ1 .

(51)

We verified this explicitly for n = 1, 2, 3, 4.

4 The physical part of the construction As we have seen above we need the functions ω(µ1 , µ2 ), ω′ (µ1 , µ2 ) and ϕ(µ) for the description of the physical correlation functions in the limit α → 0. A relatively simple description of these functions in terms of an auxiliary function a and a generalized magnetization density G was given in [2]. We call it the a-formulation. The a-formulation is useful for deriving multiple integral formulae [9] and for studying the high-temperature expansion of the free energy and the correlation functions [18]. For the purpose of the actual numerical calculation of thermodynamic properties [15, 16] or short-distance correlators from the multiple integral [7] one has to resort to a different formulation we refer to as the bb-formulation. It needs pairs of functions, but the defining integral equations in the massless regime |∆| < 1 involve only the real axis as integration contour and become trivial in the zero temperature limit. It turns out that the bb-formulation is numerically extremely stable (see below). Switching from the a- to the bb-formulation is a rather standard procedure which basically requires the application of the Fourier transformation (see e.g. [7]). For this reason we present only the result. We restrict ourselves to the massless regime |∆| < 1 and set γ = −iη ∈ R. Let us define a pair of auxiliary functions as the solution of the non-linear integral equations πh 2πJ sin(γ) ln b(x) = − − + 2(π − γ)T T γ ch(πx/γ)

Z ∞ dy

F(x − y) ln(1 + b(y)) 2π dy − F(x − y + η− ) ln(1 + b(y)) , −∞ 2π

ln b(x) =

πh 2πJ sin(γ) − + 2(π − γ)T T γ ch(πx/γ) −

with kernel

−∞ Z ∞

Z ∞ dy

−∞ 2π Z ∞ dy −∞ 2π

F(x − y) ln(1 + b(y)) F(x − y − η− ) ln(1 + b(y))

sh ( π2 − γ)k eikx F(x) = dk . −∞ 2 sh (π − γ) 2k ch γk 2 Z ∞

(52a)

(52b)

(53)

Note that the physical parameters temperature T , magnetic field h, and coupling J enter only through the driving terms of equations (52) into our formulae. In this formulation

12 (52) is valid for 0 < γ ≤ π/2 which means 0 ≤ ∆ < 1. In order to access negative ∆ we have to reverse the sign of J and conjugate the local operators whose correlation functions we are interested in by σz on every second site of the spin chain (see [21]). Except for the auxiliary functions b and b we need two more pairs of functions (±) ′ gµ and g′ (±) µ in order to define ω and ω . Both pairs satisfy linear integral equations involving b and b, (+) π(x−µ) gµ (x) = iπγ sech γ +

(−)

gµ (x) =

Z ∞ dy F(x − y) −∞

2π 1 + b−1 (y)

π(x−µ) iπ γ sech γ

+

(+)

gµ (y) −

Z ∞ dy F(x − y)

gµ (y) −

−∞

(−)

2π 1 + b−1 (y)

Z ∞ dy F(x − y + η− ) (−) −∞

2π 1 + b−1 (y)

gµ (y) , (54a)

Z ∞ dy F(x − y − η− ) (+) −∞

2π 1 + b−1 (y)

gµ (y)

(54b)

and (+)

g′ µ (x) =

(−) g′ µ (x) =

iπ π γ (x − µ) − 2

+γ

π(x−µ) γ

Z ∞ dy D(x − y)

(+) gµ (y) − γ

Z ∞ dy D(x − y + η− ) (−)

gµ (y) 2π 1 + b−1 (y) Z ∞ Z ∞ dy F(x − y + η− ) ′ (−) dy F(x − y) ′ (+) g (y) − g µ (y) , (55a) + µ −1 −1 −∞ 2π 1 + b (y) −∞ 2π 1 + b (y) −∞

2π 1 + b−1 (y)

iπ π γ (x − µ) + 2

+γ +

where

sech

sech

π(x−µ) γ

Z ∞ dy D(x − y) −∞

2π 1 + b−1 (y)

Z ∞ dy F(x − y) −∞

2π 1 + b−1 (y)

(−) gµ (y) − γ

(−) g′ µ (y) −

−∞

Z ∞ dy D(x − y − η− ) (+) −∞

2π 1 + b−1 (y)

Z ∞ dy F(x − y − η− ) ′ (+) g µ (y) , (55b) −1 −∞

2π 1 + b

(y)

ch ( π2 − γ)k D(x) = dk . −∞ 4 sh2 (π − γ) 2k ch2 γk 2 Z ∞

sin(kx) sh

gµ (y)

πk 2

(56)

The functions ω(µ1 , µ2 ), ω′ (µ1 , µ2 ) and ϕ(µ) that determine the explicit form of the inhomogeneous correlation functions of the XXZ chain can be written as integrals (±) involving b, b, gµ and g′ (±) µ . The function (+) (−) gµ (x) gµ (x) dx − , ϕ(µ) = −1 −1 −∞ 2(π − γ)i 1 + b (x) 1 + b (x) Z ∞

(57)

13 also determines the magnetization m(T, h) = − 12 ϕ(0)

(58)

which is the only independent one-point function of the XXZ chain. The function k sh (π − γ) 2 cos(k(µ1 − µ2 )) dk ω(µ1 , µ2 ) = − 21 K(µ1 − µ2 ) − γk −∞ ch i sh πk 2 2 Z ∞

− with

Z ∞ dx

γ

−∞

sech

π(x−µ2 ) γ

(−) g(+) (x) gµ1 (x) µ1 + (59) 1 + b−1 (x) 1 + b−1 (x)

K(µ) = cth(µ − η) − cth(µ + η)

(60)

also determines the internal energy y

y

e(T, h) = Jhσxj−1 σxj + σ j−1 σ j + ∆(σzj−1 σzj − 1)iT,h = J sh(η)ω(0, 0) − ∆ of the XXZ chain. The function ω′ (µ1 , µ2 ) is defined as ω (µ1 , µ2 ) = ′

η (+) (µ1 − µ2 ) + 2K

−

Z ∞ dx

+

−∞

γ

sech

Z ∞ dx −∞

γ

Z ∞

dk

−∞

π(x−µ2 ) γ

(61)

γ sin(k(µ1 − µ2 )) 2 γk 2i th πk 2 ch 2

(−) f (+) (x) fµ1 (x) µ1 + 1 + b−1 (x) 1 + b−1 (x)

(x − µ2 ) sech

π(x−µ2 ) γ

(−) g(+) (x) gµ1 (x) µ1 + , (62) 1 + b−1 (x) 1 + b−1 (x)

where K (+) (µ) = cth(µ − η) + cth(µ + η) ,

(±)

(±)

(±)

fµ (x) = g′ µ (x) ∓ iγ2 gµ (x) .

(63)

For the function ω′ (0, 0) itself we did not find a simple interpretation in terms of thermodynamic quantities. Its derivative ω′x = ∂µ ω′ (µ, 0)|µ=0, however, appears in the homogeneous density matrix D2 (T, h) and, hence, is related to the neighbour correlators (see [2]). From the symmetry ω(µ1 , µ2 , α) = ω(µ2 , µ1 , −α) (see [2]) we conclude that ω(µ1 , µ2 ) = ω(µ2 , µ1 ) ,

ω′ (µ1 , µ2 ) = −ω′ (µ2 , µ1 ) .

(64)

This property was crucial for the derivation of (37). For the calculation of the physical correlation functions one has to perform the homogeneous limit. This limit is rather singular, because the coefficients coming from the

14 algebraic part in general have poles when two inhomogeneities coincide. These poles are canceled by zeros stemming from certain symmetric combinations of the functions ω(λ1 , λ2 ), ω′ (λ1, λ2 ) and ϕ(λ) in the numerators. In order to perform the limit one has to apply l’Hˆospital’s rule, finally leading to polynomials in the functions ω(0, 0), ω′ (0, 0), ϕ(0) and derivatives of these functions with respect to the inhomogeneity parameters evaluated at zero. We shall denote derivatives with respect to the first argument by subscripts x and derivatives with respect to the second argument by subscripts y and leave out zero arguments for simplicity. Then e.g. ω′xyy = ∂x ∂2y ω′ (x, y)|x,y=0 etc. For the examples in the next section the non-linear integral equations for b and b (±) as well as their linear counterparts for gµ and g′ (±) µ were solved iteratively in Fourier space, utilizing frequently the fast Fourier transformation algorithm. The derivatives (±) of gµ and g′ (±) µ with respect to µ, needed in the computation of the respective deriva′ tives of ω, ω and ϕ, satisfy linear integral equations as well, which were obtained (±) as derivatives of the equations for gµ and g′ (±) µ . Taking into account derivatives is particularly simple in Fourier space.

5 Examples of short-distance correlators Here we mostly concentrate on the longitudinal and transversal two-point functions hσz1 σzn iT,h, hσx1 σxn iT,h for n = 2, 3, 4. The algebraic part for n = 2, 3 was already calculated in the more general inhomogeneous case in our previous paper [2]. In the inhomogeneous case for n = 2 we obtained 1 1 tr12 eΩ1 (1 + Ω2 )σz1 σz2 = tr12 Ω1 σz1 σz2 4 4 1 ′ − z z = − tr12 [ω(λ1 , λ2 )Ω+ + ω (λ , λ )Ω ]σ σ 1 2 1,2 1,2 1 2 4 cth(λ1 − λ2 ) ′ = cth(η)ω(λ1, λ2 ) + ω (λ1, λ2 ) (65) η and in the homogeneous limit λ j → 0, hσz1 σz2 iT,h = cth(η)ω +

ω′x . η

(66)

Note that the limit exists, because ω′ (λ1 , λ2 ) is antisymmetric. Similarly for the transversal correlation functions, 1 1 tr12 eΩ1 (1 + Ω2 )σx1 σx2 = tr12 Ω1 σx1 σx2 4 4 1 = − tr12 [ω(λ1 , λ2 )Ω+ + ω′ (λ1 , λ2 )Ω− ]σx1 σx2 1,2 1,2 4 ch(η) ch(λ1 − λ2 ) ω(λ1 , λ2 ) − ω′ (λ1, λ2 ) (67) =− 2 sh(η) 2η sh(λ1 − λ2 )

15 which in the homogenous limit becomes hσx1 σx2 iT,h = −

ω ch(η)ω′x − . 2 sh(η) 2η

(68)

The two-point functions hσz1 σzn iT,h , hσx1 σxn iT,h are even under spin reversal. Therefore we do not expect any non-vanishing contribution from Ω2 , which would come with functions ϕ, i.e. odd functions of the magnetic field h. This is indeed what we observe from our calculation for n = 1, 2, 3, 4. The above two-point correlation functions depend only on ω, ω′ and their derivatives. As an example of a correlation function which is of no definite symmetry under spin reversal and which does depend on ϕ and its derivatives we shall discuss the emptiness formation probability below. We obtained the formula for n = 3, 4 in the same way as for n = 2 by first calculating the inhomogeneous generalization of the correlation function and performing the homogeneous limit at the end. The formula for hσz1 σz3 iT,h in the inhomogeneous case was presented in [2]. For space limitations we do not repeat it here and also do not show the inhomogeneous formulae in the remaining cases. Instead we merely list the results for the physical correlation functions. For n = 3 we already obtained in [2] hσz1 σz3 iT,h

2 ω′x th(η)(ωxx − 2ωxy ) sh (η)ω′xxy = 2 cth(2η)ω + + − , η 4 4η

hσx1 σx3 iT,h = −

(69a)

ch(2η) ′ 1 ω− ωx sh(2η) 2η ch(2η) th(η)(ωxx − 2ωxy ) sh2 (η)ω′xxy + . − 8 8η

(69b)

The length of the formulae grows rapidly with the number of lattice sites. For n = 4 we found hσz1 σz4 iT,h =

1 768q4 (−1 + q6 ) (1 + q2 ) η2

2 384q4 1 + q2 5 − 4q2 + 5q4 η2 ω 2 4 6 8 12 4 η2 ωxy − 8 1 + q 52 + 64q − 234q + 64q + 52q + q 2 + 192q4 −1 + q2 1 + 4q2 + q4 η2 ωyy i 4 h + −1 + q2 1 + q4 1 + 4q2 + q4 η2 −4ωxyyy + 6ωxxyy − 768q4 −1 − q2 + q6 + q8 ηω′ y 2 4 6 8 10 2 3 1 + 6q + 11q + 11q + 6q + q ηω′ xyy + 16 −1 + q 5 1 + 2q2 + 2q4 + q6 ηω′ xxyyy − 2 −1 + q2

16 i h 3 + 8 −1 + q2 1 + q2 1 + 6q2 + 34q4 + 6q6 + q8 η2 ω2y − ωωxy h 2 4 6 10 12 14 16 + −1 − 4q − 22q − 12q + 12q + 22q + 4q + q η2 −6ω2yy i + 12ωyy ωxy + 4ωy ωyyy − 12ωy ωxyy − 4ωωxyyy + 6ωωxxyy i 4 2 h + 16 −1 + q2 1 + q2 1 + q2 + q4 η ωyy ω′ y − ωy ω′ yy + ωω′ xyy h 2 4 6 8 4 2 1 + 5q + 6q + 5q + q η 4ωxyyy ω′ y − 6ωxxyy ω′ y − 2ωyyy ω′ yy + −1 + q i + 6ωxyy ω′ yy + 2ωyy ω′ yyy − 4ωxy ω′ yyy − 6ωyy ω′ xyy + 4ωy ω′ xyyy − 2ωω′ xxyyy i h ′ ′ ′ ′ ′ ′ 4 3 2 4 (70) + 3 −1 + q 1 + q + q ω yyy ω xyy − ω yy ω xyyy + ω y ω xxyyy and in the transversal case, 1 hσx1 σx4 iT,h = 3072q5 (−1 + q6 ) η2 − 768q6 1 + 10q2 + q4 η2 ω 2 2 2 2 4 6 8 + 16q −1 + q 31 + 56q − 30q + 56q + 31q η2 ωxy 2 3 + 5q2 − 4q4 + 5q6 + 3q8 η2 ωyy − 96q2 −1 + q2 i 4 h + q2 −1 + q2 1 + 4q2 + q4 η2 8ωxyyy − 12ωxxyy + 192q2 −3 − q2 − q4 + q8 + q10 + 3q12 ηω′ y 2 3 2 4 6 8 10 1 − 12q − 25q − 25q − 12q + q ηω′ xyy + 8 −1 + q 5 1 + 2q2 + 2q4 + q6 ηω′ xxyyy + 2 −1 + q2 i h 3 17 + 7q2 + 7q4 + 17q6 η2 ωωxy − ω2y + 16q2 −1 + q2 2h 2 2 4 8 10 12 + q −5 − 4q − 13q + 13q + 4q + 5q η 12ω2yy − 24ωyy ωxy i − 8ωy ωyyy + 24ωy ωxyy + 8ωωxyyy − 12ωωxxyy h i ′ ′ ′ 2 4 6 8 2 4 1 − 9q − 8q − 9q + q η ωyy ω y − ωy ω yy + ωω xyy + 8 −1 + q h 2 + −1 + q4 1 + 5q2 + 6q4 + 5q6 + q8 η −4ωxyyy ω′ y + 6ωxxyy ω′ y

i +2ωyyy ω′ yy −6ωxyy ω′ yy −2ωyy ω′ yyy +4ωxy ω′ yyy +6ωyy ω′ xyy −4ωy ω′ xyyy +2ωω′ xxyyy i h ′ 4 3 2 4 ′ ′ ′ ′ ′ + 3 −1 + q 1 + q + q −ω yyy ω xyy + ω yy ω xyyy − ω y ω xxyyy . (71)

17

0 0.2 −0.1

〈 σz1 σz2 〉

−0.3 −0.2 −0.4 ∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995

−0.4 −0.6 1e−03

1e−01

0.3

1e+01 T/J

1e+03

∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995 1e+05

1e−03

1e−01

∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995

0.25 0.2

1e+01 T/J

1e+03

〈 σx1 σx2 〉

−0.2

0

−0.5 −0.6 1e+05

∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995

0.5 0.4 0.3

0.1 0.05

0.2

〈 σx1 σx3 〉

〈 σz1 σz3 〉

0.15

0 0.1 −0.05 0

−0.1 1e−03 1e−01 1e+01 1e+03 1e+05 T/J

1e+01 T/J

1e+03

1e+05

∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995

0.2

〈 σz1 σz4 〉

1e−01

0 −0.1

0.1

−0.2

0

−0.3 ∆ = −0.995 ∆ = −0.707 ∆ = 0.010 ∆ = 0.707 ∆ = 0.995

−0.1 −0.2 1e−03

1e−01

1e+01 T/J

1e+03

1e+05

1e−03

1e−01

1e+01 T/J

1e+03

〈 σx1 σx4 〉

0.3

1e−03

−0.4 −0.5 1e+05

Figure 1: Two-point correlators for n = 2, 3, 4 as a function of temperature for various values of ∆ in the massless regime at zero magnetic field.

18 ∆

n=2

(L = 18)

n=3

(L = 18)

n=4

(L = 18)

−0.1 −0.2 −0.3 −0.4 −0.5 −0.6 −0.7 −0.8 −0.9

4.96645 2.43157 1.56079 1.10294 0.806967 0.588818 0.412795 0.262355 0.129195

(4.966) (2.432) (1.561) (1.103) (0.807) (0.589) (0.413) (0.265) (0.137)

3.32288 1.64332 1.07081 0.771287 0.577718 0.434179 0.316321 0.211402 0.111127

(3.323) (1.643) (1.071) (0.771) (0.578) (0.434) (0.318) (0.215) (0.118)

2.56077 1.27520 0.839169 0.611558 0.4641a 0.354030 0.262606 0.179803 0.098055

(2.561) (1.275) (0.839) (0.612) (0.464) (0.355) (0.264) (0.184) (0.104)

Table 1: Numerical values of the cross-over temperature T0 (n; ∆) on the infinite chain and for a chain of 18 sites in units of J/2 (values in brackets, taken from [8]). a

The prefactor on the rhs of equation (70) has a pole as a function of q when q is a third root of unity corresponding to ∆ = −0.5. We estimated the value of T0 (n; ∆) for ∆ = −0.5 by interpolating between two values closely above and below −0.5.

Using the representations of the previous section for the functions ω and ω′ we can determine high-precision numerical values for the various two-point correlators. Figures 1-3 show selected examples. For n = 2, 3, 4 the longitudinal correlation functions hσz1 σzn iT,h are exhibited in the left panels, while the right panels show the transversal correlation functions hσx1 σxn iT,h. We observe a rich scenario for their temperature and field dependence, in particular non-monotonous behaviour at intermediate temperatures. The numerical precision in all figures is between 3 and 4 significant digits. In figure 1 we display our data for vanishing magnetic field and various values of ∆ ∈ (−1, 1). The transversal correlations are alternating in sign with n for all values of ∆. But for the longitudinal correlations their qualitative behaviour depends on the sign of ∆. They alternate for positive ∆. For negative ∆ they are always negative at sufficiently low temperatures, change their sign at a ‘crossover temperature’ T0 (n; ∆), reach a positive maximum and tend to zero from above in the high temperature limit. This phenomenon was studied numerically on finite chains of 16 and 18 sites in [8] and was interpreted there as a ‘quantum-to-classical crossover’. In [8] it was estimated that the crossover temperature T0 (n; ∆) remains finite in the thermodynamic limit. From our formulae we can calculate it with high precision (see table 1). Note the peculiar low temperature behaviour close to ∆ = −1, where one would expect simultaneously isotropy (hσz1 σzn i = (−1)n+1 hσx1 σxn i) as well as full polarization (hσz1 σzn i + 2(−1)n+1 hσx1 σxn i = 1). However, this is not seen, not even for ∆ as close to the isotropic point as −0.995, since the full polarization is still forced into the xy-plane by the slightly dominant transversal exchange. At elevated but not too high temperature the system is still close to full polarization, but the small anisotropy is irrelevant and longitudinal as well as transversal correlation attain values close to ±1/3. At much

19

1

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.8

−0.1

−0.3

0.2 0

−0.4 h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

−0.2 −0.4 −0.6 1e−03

1e−01

1e+01 T/J

1

1e+03

1e+05

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.8 〈 σz1 σz3 〉

〈 σx1 σx2 〉

−0.2

0.4

1e−03

1e−01

1e+01 T/J

1e+03

−0.5 −0.6 1e+05

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.6

0.35 0.3 0.25 0.2 0.15

0.4

〈 σx1 σx3 〉

〈 σz1 σz2 〉

0.6

0

0.1 0.2

0.05 0

0 1e−03

1e−01

1e+01 T/J

1

1e+03

1e+05

1e−03 1e−01 1e+01 1e+03 1e+05 T/J

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.8

0 −0.05 −0.1

0.4 −0.15 0.2

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0 −0.2 1e−03

1e−01

1e+01 T/J

1e+03

1e+05

〈 σx1 σx4 〉

〈 σz1 σz4 〉

0.6

−0.2 −0.25

1e−03 1e−01 1e+01 1e+03 1e+05 T/J

√ Figure 2: Two-point correlators for n = 2, 3, 4 and ∆ = 1/ 2 as a function of temperature for various magnetic fields.

20 higher temperatures all correlations tend to zero. As was explained in [8] the negative sign of the longitudinal correlations at low temperatures is a pure √ quantum effect. In figure 2 we plotted the correlators for a fixed value of ∆ = 1/ 2 for various magnetic fields. The curves show non-monotonous behaviour for intermediate fields. In particular, a sign change is possible for positive ∆ if the magnetic field is switched on, as can be seen in the upper and lower left panels. In general, we observe a competition between the effect of the magnetic field, which tends to align the spins in z-direction, the longitudinal exchange interaction in the Hamiltonian, which enforces antiparallel alignment of neighbouring spins, and the transversal exchange interaction, which can be interpreted in terms of quantum fluctuations. The temperature sets the scale for the relative and absolute relevance of these terms. At low temperatures the quantum fluctuations dominate, decreasing e.g. the value of hσz1 σz3 i from 1 to a smaller positive value. When increasing the temperature, at intermediate field the quantum fluctuations are suppressed and the field is less hindered to align the spins, before at very high temperature the thermal fluctuations destroy all order. For large enough fields (light blue curves) the correlation functions reach saturation. This can be also observed in figure 3, where the correlators are shown as functions of the magnetic field for relatively low temperature, T /J = 0.1, and where the saturation roughly sets in at its zero temperature value h = (1 + ∆)4J. An interesting feature is the decrease of hσz1 σz3 i for small fields which we interpret as due to a weakening of the antiferromagnetic neighbour correlations. As an independent test of our results, we numerically calculated the correlation functions for finite systems. Since the Hilbert space dimension grows exponentially with the chain length, this required us to characterise the spectrum and eigenvectors of high-dimensional sparse matrices. The complete diagonalisation of such matrices is too time consuming or simply impossible due to memory limitations. We therefore used Krylov space methods, in particular, Chebyshev expansion [20]. Given the D eigenstates |ki of the system, we split the expectation value of an operator A into a low-energy and a high-energy part, −H/T

Z = Tr(e

C−1

)=

∑e

−Ek /T

k=0

1 1 hAi = Tr(A e−H/T ) = Z Z

"

D−1

+

∑ e−Ek /T ,

(72)

k=C

C−1

D−1

k=0

k=C

∑ hk|A|ki e−Ek /T + ∑ hk|A|ki e−Ek/T

#

,

(73)

where H is the Hamiltonian of the system, Z the partition function, and Ek the energy of the eigenstate |ki. The C low-energy states and the corresponding matrix elements can be calculated exactly with the Lanczos algorithm [17]. The contribution of the remaining states is estimated via Chebyshev expansions of the density of states, D−1

∑e

k=C

−Ek /T

=

Z∞ −∞

ρ(E) e−E/T dE ,

(74)

21

1

0

0.8

−0.1

0.6

−0.2

−0.4

0 −0.2

∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995

−0.4 −0.6 −0.8 0

2

4

6

10

0

2

4

h

6

−0.6 −0.7 8

−0.8 10

h ∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995

1 0.8 0.6 〈 σz1 σz3 〉

−0.5

∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995 8

0.5 0.4 0.3

0.4 0.2

0.2 ∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995

0 −0.2 −0.4 0

2

4

〈 σx1 σx2 〉

−0.3

0.2

6

〈 σx1 σx3 〉

〈 σz1 σz2 〉

0.4

0.1 0 8

10

0

2

4

h

6

8

10

h 0

1

−0.05

0.8

−0.1 −0.15 0.4

−0.2

0.2

∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995

0 −0.2 0

2

4

6 h

−0.25

∆ = 0.010 ∆ = 0.309 ∆ = 0.707 ∆ = 0.951 ∆ = 0.995 8

10

0

2

4

6

〈 σx1 σx4 〉

〈 σz1 σz4 〉

0.6

−0.3 −0.35 8

10

h

Figure 3: Two-point correlators for n = 2, 3, 4 and T /J = 0.1 as a function of magnetic field for various values of ∆.

22

1

L = 20 L = 24 L = 28 L = 32 L→∞

0.8

0 −0.05 −0.1

0.4 −0.15 0.2

h = 0.0 h = 2.0 h = 4.0 h = 6.0 h = 8.0

0 −0.2 1e−02

1e+00 T/J

1e+02

1e−02

1e+00 T/J

〈 σx1 σx4 〉

〈 σz1 σz4 〉

0.6

−0.2 −0.25 1e+02

√ Figure 4: Two-point correlators for n = 4 and ∆ = 1/ 2 as a function of temperature, comparison with data from direct diagonalization of finite chains. ρ(E) =

D−1

1

∑ δ(E − Ek ) = πps2 − (E − E) ¯ 2

k=C

i h N−1 ρ ρ E−E¯ µ0 + 2 ∑ µn Tn ( s ) ,

¯ µρn = gN n Tr(PTn [(H − E)/s]) with P = 1 −

(75)

n=1 C−1

∑ |kihk| ,

(76)

k=0

and of the diagonal matrix elements of A, D−1

∑ hk|A|ki e

−Ek /T

k=C

=

Z∞

a(E) e−E/T dE ,

−∞

i h N−1 1 a a E−E¯ p ) µ + 2 hk|A|kiδ(E − E ) = µ T ( n k ∑ ∑ n s , ¯ 2 0 π s2 − (E − E) n=1 k=C ¯ µan = gN n Tr(PATn [(H − E)/s]) . D−1

(77)

a(E) =

(78) (79)

¯ s denote scaling factors Here Tn are the Chebyshev polynomials of the first kind and E, that shift the spectrum of H into the domain of the polynomials [−1, 1]. The coefficients gN n account for a special kernel [12, 20] which improves the convergence of ρ the truncated Chebyshev series. The traces defining the moments µn and µan can be estimated stochastically, Tr(X ) ≈

D−C R

R−1

∑ hr|X |ri

r=0

with R ≪ D random states |ri .

(80)

Thus the most time consuming part of the moment evaluation is the recursive calcula¯ tion of Tn [(H − E)/s]|ri for a few dozen states |ri. The whole procedure is feasible for Hilbert space dimensions of 107 to 108 . Employing Sz and momentum conservation,

23

1

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.8

1

h = 0.0 h = 2.0 h = 5.0 h = 6.5 h = 8.0

0.8

0.4

0.4

0.2

0.2

P(3)

0.6

P(2)

0.6

0

0 1e−03

1e−01

1e+01 T/J

1e+03

1e+05

1e−03

1e−01

1e+01 T/J

1e+03

1e+05

√ Figure 5: Emptiness formation probability for n = 2, 3 and ∆ = 1/ 2 as a function of temperature for various magnetic fields. with moderate effort we can reach a chain length of L = 32 (about 430 CPU hours on a 2.66 GHz Xeon Cluster). We also performed tests against other exact results. We compared the high temperature expansions of equations (70), (71) with the high temperature expansions for the multiple integrals from [10] up to third order in 1/T . We compared with the known ground state results of [13, 14]. We compared the zz-correlation functions in the XX limit with the explicit result (see e.g. [11]), and we computed the correlation functions for ∆ very close to 1 and compared with the results obtained for the isotropic limit and n = 2, 3 in [1]. As an example of a correlation function which has no definite symmetry with respect to the reversal of all spins and hence must depend on ϕ and its derivatives we consider the emptiness formation probability P(n) = 2−n h∏nj=1 (1 + σzj )iT,h. We know from [2] that 1 ϕ cth(η)ω ω′x P(2) = − + + (81) 4 2 4 η and 1 + q2 + q4 ω −1 + q2 (ωyy − 2ωxy ) 3ω′y 1 3ϕ P(3) = − + + − 8 8 2 (−1 + q4 ) 32 (1 + q2 ) 8η 2 ′ 2 −1 + q ωxyy 3 1 + q2 (−ωϕxx + 2ωy ϕx + ωyy ϕ − 2ωxy ϕ) + + 128q2 η 32 (−1 + q2 ) 1 + 10q2 + q4 (ω′xyy ϕ − ω′yy ϕx + ω′y ϕxx ) + . (82) 128q2 η P(2) and P(3) are displayed in figure 5. They show the characteristic high temperature behaviour P(n) ∼ 2−n and saturate for small temperature and h > (1 + ∆)4J. Like the

24 two-point functions they exhibit non-monotonous behaviour at intermediate magnetic fields.

6 Conclusion The aim of this work was to demonstrate that short-distance correlation functions of the XXZ chain can be computed efficiently from the exponential form of the density matrix proposed in [2]. We found that an accurate computation is possible for arbitrary temperatures and magnetic fields in the massless regime −1 < ∆ < 1, even beyond the phase boundary to the fully polarized state. This has to be contrasted with our experience [7] in computing the same correlation functions from the multiple integral representations obtained in [9,10], which appeared to be much harder. In fact, in [7] we were unable to proceed beyond n = 3. In the present approach the difficulty for larger n comes mainly from the algebraic part of the calculation. We believe that ranges of the order of n = 10 may be reached if one works directly with the homogeneous transfer matrices (which would mean to calculate higher order residua in (9)). Acknowledgment. HB, JD and FG gratefully acknowledge financial support by the Volkswagen Foundation and by the DFG through Graduiertenkolleg 1052. JS was supported by a Grant-in-Aid for Scientific Research #20540370 from the Ministry of Education in Japan.

References [1] H. Boos, F. G¨ohmann, A. Kl¨umper, and J. Suzuki, Factorization of multiple integrals representing the density matrix of a finite segment of the Heisenberg spin chain, J. Stat. Mech. (2006), P04001. [2]

, Factorization of the finite temperature correlation functions of the XXZ chain in a magnetic field, J. Phys. A 40 (2007), 10699.

[3] H. Boos, M. Jimbo, T. Miwa, F. Smirnov, and Y. Takeyama, Reduced qKZ equation and correlation functions of the XXZ model, Comm. Math. Phys. 261 (2006), 245. [4]

, Fermionic basis for space of operators in the XXZ model, preprint, hepth/0702086, 2007.

[5]

, Hidden Grassmann structure in the XXZ model, Comm. Math. Phys. 272 (2007), 263.

[6]

, Hidden Grassmann structure in the XXZ model II: creation operators, preprint, arXiv:0801.1176, 2008.

25 [7] M. Bortz and F. G¨ohmann, Exact thermodynamic limit of short range correlation functions of the antiferromagnetic XXZ chain at finite temperatures, Eur. Phys. J. B 46 (2005), 399. [8] K. Fabricius and B. M. McCoy, Quantum-classical crossover in the spin-1/2 XXZ chain, Phys. Rev. B 59 (1999), 381. [9] F. G¨ohmann, A. Kl¨umper, and A. Seel, Integral representations for correlation functions of the XXZ chain at finite temperature, J. Phys. A 37 (2004), 7625. [10]

, Integral representation of the density matrix of the XXZ chain at finite temperature, J. Phys. A 38 (2005), 1833.

[11] F. G¨ohmann and A. Seel, XX and Ising limits in integral formulae for finite temperature correlation functions of the XXZ chain, Theor. Math. Phys. 146 (2006), 119. [12] Dunham Jackson, On approximation by trigonometric sums and polynomials, Trans. Amer. Math. Soc. 13 (1912), 491–515. [13] G. Kato, M. Shiroishi, M. Takahashi, and K. Sakai, Next-nearest-neighbour correlation functions of the spin-1/2 XXZ chain at the critical region, J. Phys. A 36 (2003), L337. [14]

, Third-neighbour and other four-point correlation functions of spin-1/2 XXZ chain, J. Phys. A 37 (2004), 5097.

[15] A. Kl¨umper, Free energy and correlation length of quantum chains related to restricted solid-on-solid lattice models, Ann. Physik 1 (1992), 540. [16]

, Thermodynamics of the anisotropic spin-1/2 Heisenberg chain and related quantum chains, Z. Phys. B 91 (1993), 507.

[17] Cornelius Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Res. Nat. Bur. Stand. 45 (1950), 255–282. [18] Z. Tsuboi and M. Shiroishi, High temperature expansion of the emptiness formation probability for the isotropic Heisenberg chain, J. Phys. A 38 (2005), L363. [19] J. A. M. Vermaseren, New features of FORM, preprint, arXiv:math-ph/0010025, 2000. [20] Alexander Weiße, Gerhard Wellein, Andreas Alvermann, and Holger Fehske, The kernel polynomial method, Rev. Mod. Phys. 78 (2006), 275–306. [21] C. N. Yang and C. P. Yang, Ground-state energy of a Heisenberg-Ising lattice, Phys. Rep. 147 (1966), 303.